Theory of AC quantum transport with fully electrodynamic coupling
NNoname manuscript No. (will be inserted by the editor)
Theory of AC quantum transport with fully electrodynamic coupling
Timothy M. Philip · Matthew J. Gilbert
Received: November 14, 2018/ Accepted: date
Abstract
With the continued scaling of microelectronic de-vices along with the growing demand of high-speed wire-less telecommunications technologies, there is increasingneed for high-frequency device modeling techniques that ac-curately capture the quantum mechanical nature of chargetransport in nanoscale devices along with the dynamic fieldsthat are generated. In an effort to fill this gap, we developa simulation methodology that self-consistently couples ACnon-equilibrium Green functions (NEGF) with the full so-lution of Maxwell’s equations in the frequency domain. Weapply this technique to simulate radiation from a quantum-confined, quarter-wave, monopole antenna where the length L is equal to one quarter of the wavelength, λ . Classi-cally, such an antenna would have a narrower, more di-rected radiation pattern compared to one with L (cid:28) λ , butwe find that a quantum quarter-wave antenna has no direc-tivity gain compared to the classical solution. We observethat the quantized wave function within the antenna sig-nificantly alter the charge and current density distributionalong the length of the wire, which in turn modifies the far-field radiation pattern from the antenna. These results showthat high-frequency radiation from quantum systems can bemarkedly different from classical expectations. Our method,therefore, will enable accurate modeling of the next genera-tion of high-speed nanoscale electronic devices. Keywords
Quantum Transport · High-frequency · An-tenna · Radiation
This work was supported by the National Science Foundation underCAREER ECCS-1351871.T. M. PhilipDepartment of Electrical and Computer Engineering, University of Illi-nois at Urbana-Champaign, Urbana, IL 61801, USAE-mail: [email protected]. J. GilbertE-mail: [email protected]
The non-equilibrium Green function (NEGF) methodologyhas found great utility in modeling DC electron transportin a wide variety of quantum systems such as quantumdots [1,2,3], resonant tunneling diodes [4,5], and metal-oxide-semiconductor (MOS) transistors [6,7,8]. The suc-cess of the method derives not only from its ability to ac-curately model quantum coherent transport but also from itssystematic framework within which to introduce incoher-ence and interactions [9,10,11,12]. The standard, DC for-mulation of the NEGF method, however, fails to capturetransient features of transport and electrodynamic couplingthat are important in phenomena such as cross-talk in high-frequency circuits [13]. Although, such full-wave couplinghas been treated semi-classically in the context of the time-domain Boltzmann transport and the drift-diffusion equa-tions [14,13,15,16,17], fully dynamic electromagnetic cou-pling to quantum transport has not yet been achieved.The time-domain NEGF formulation can be used to un-derstand dynamic phenomena and has successfully been ap-plied to elucidate phenomena ranging from the AC conduc-tance in graphene [18] to the transient properties of quantuminterferometers [19,20,21]. Transient NEGF calculations,however, suffer from the memory effect, whereby the entirehistory of the Green function must be stored in order to accu-rately evolve the Green function [22,23]. Although some ap-proximations, such as the wide bandwidth limit, can greatlyreduce the storage requirements, time-dependent calcula-tions are still demanding for time scales longer than a fewpicoseconds since a sub-femtosecond time step is requiredfor numerical stability and the computation time scales atleast linearly with the total number of time steps. There-fore, to accurately model device responses at sub-terahertzoperating frequencies, the DC NEGF formulation has beenextended to simulate AC steady state transport, thereby pro- a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n Timothy M. Philip, Matthew J. Gilbert viding a platform for studying quantum transport at any fre-quency [24,25,26,27].Although AC NEGF provides a simple extension of DCNEGF to understand high-frequency transport, care mustbe taken to accurately account for displacement currentand thereby satisfy gauge invariance and charge conserva-tion [24,27]. The AC NEGF methodology on its own doesnot satisfy either of these conditions, and so they must be en-forced through additional measures. Fortunately, this prob-lem can be alleviated by self-consistently solving the trans-port equations with the Coulomb potential at the Hartreelevel of approximation [24,27]. To date, however, electro-magnetic interactions with AC NEGF transport have beenassumed to be quasi-static, where a solution of Poisson’sequation is sufficient to capture the electrostatics of theproblem. The quasi-static approximation fails, however, insituations that involve dynamic electromagnetic fields, suchas electromagnetic induction or field radiation, which arefeatures that are vital to characterize the properties of high-frequency devices [28].In this work, we present a simulation technique that self-consistently solves the AC NEGF equations with the full so-lution of Maxwell’s equations in the frequency domain toconsider the influence of both static and dynamic fields inthe presence of fully quantum mechanical electron transportin three dimensions. In Sect. 2, we introduce details of thesimulation methodology, including the tight-binding form ofHamiltonian and the governing equations of the method forquantum transport, which are used to calculate charge andcurrent density within the system. We then present a finite-difference frequency-domain (FDFD) method that, insteadof solving for the electric and magnetic fields, solves di-rectly for the scalar and vector electromagnetic potentialson a Yee cell. By solving for the potential formulation ofMaxwell’s equations rather than the electric and magneticfields that standard treatments obtain, we can directly cou-ple the output potentials into the Hamiltonian. The outputcharge and current densities from the AC NEGF are coupledwith the FDFD potentials to obtain a self-consistent solu-tion. In Sec. 3, we apply this technique to model the electro-dynamic radiation emitted by a quantum monopole antenna,that is an antenna that possesses quantum-confined states.We begin with an overview of the classical operation of amonopole antenna, with emphasis on the expected radiationpatterns at different operating frequencies. We then comparethese classical expectations to the the AC NEGF/FDFD sim-ulations of a quantum monopole. The quantization withinthe quantum wire significantly alters the charge and currentdensity distribution, which ultimately changes the macro-scopic radiation pattern. Despite driving the antenna at thequarter-wave frequency, we find that the quantum monopoleantenna fails to achieve the directivity gain associated with a classical quarter-wave antenna. In Sec. 4, we summarizeand discuss the implication of our results. H ( r ) = ∑ r (cid:34) ψ † r H ψ r + ∑ δ (cid:0) ψ † r H δ ψ r + δ + H.c. (cid:1)(cid:35) , (1)where ψ r is the electron annihilation operator at position r , ψ † r is the electron annihilation operator at position r , δ arethe distances between nearest neighbor atoms in the lattice, H is the on-site term, and H δ is the nearest-neighbor hop-ping term. For the case of the cubic lattice we consider inthis work, δ has the form δ = ( ± a ˆx , ± a ˆy , ± a ˆz ) , where a isthe lattice constant.2.2 AC NEGFThe governing equations for the AC NEGF derive from thetime-domain Keldysh equation [29,12]. The formal expres-sion for the full, two-time Green function for the total sys-tem including the leads is given by the Dyson equation [29,30] G γ ( t , t ) = g γ ( t , t )+ (cid:90) dt (cid:90) dt (cid:2) g γ ( t , t ) Σ γ ( t , t ) δ ( t − t ) G γ ( t , t ) , + g γ ( t , t ) U ( t ) δ ( t − t ) G γ ( t , t ) (cid:3) , (2)where G γ ( t , t ) is the fully-dressed retarded or advancedGreen function ( γ = { r , a } ), g γ ( t , t (cid:48) ) is the Green function forthe isolated (i.e. uncontacted) device Hamiltonian includingany time-independent potentials, Σ γ ( t , t (cid:48) ) is the self-energythat describes the time-independent and time-dependent in-fluence of the leads, and U ( t ) is the time-dependent poten-tial energy. Without loss of generality, Eq. (2) can be Fouriertransformed to the energy domain by using the two-timeFourier transform defined as [25] H ( E , E ) = (cid:90) dt π (cid:90) dt π e iE t / ¯ h e iE t / ¯ h h ( t , t ) , (3)where h ( t , t ) is any generic time-domain function of times t and t , and H ( E , E ) is the Fourier transform of h at en-ergies E and E . Since the AC bias implies that the sys-tem is no longer time-independent, we cannot utilize a sin-gle time/energy Fourier transform based on the difference t − t as is traditionally considered in DC NEGF. heory of AC quantum transport with fully electrodynamic coupling 3 The resultant energy-domain form of Eq. 2 is written G γ ( E , E ) = G γ ( E ) δ ( E − E )+ (cid:90) dE G γ ( E ) Σ γ ( E − E ) G γ ( E , E ) (cid:90) dE G γ ( E ) U ( E − E ) G γ ( E , E ) , (4)where G γ ( E ) is the DC contribution to the Green functiongiven by G r ( E ) = [( E + i η ) I − H − U − Σ r ( E )] − . (5)Here, I is the identity matrix, η is an infinitesimal positivenumber that pushes the poles of the Green function into thecomplex plane, allowing for integration along the real en-ergy axis [11], U is the DC potential energy profile, and Σ r ( E ) is the DC contact self-energy that integrates out theinfluence of the semi-infinite leads. In general, the contact-self-energy is non-analytic, but decimation techniques canbe used to efficiently obtain self-energies for semi-infiniteleads [31,32].To obtain the system response at frequency ω , we eval-uate the Green function in Eq. (4) at E = E + = E + ¯ h ω and E = E . For a AC bias of the form V ( t ) = V AC cos ω t ap-plied to the leads, we simplify the expression for the Greenfunction to G r ( E ) = G r ( E ) + g r ω ( E ) , (6)where g r ω ( E ) is first-order response of the system to the ACbias. Similarly, the total self-energy can be expressed as afunction of DC component and a small-signal AC compo-nent: Σ γ ( E ) = Σ γ ( E ) + σ γω ( E ) ( γ = r , < ) , (7)where σ γω ( E ) is the AC self-energy due to a the small-signal AC bias. Previous work has shown that this ACself-energy can be perturbatively expanded in powers of eV AC / ¯ h ω , where e is the electron charge, which provides asystematic framework to evaluate the system’s response [26]Using this result, the AC self-energy to first-order in the biasamplitude is given simply as σ γω ( E ) = eV AC ¯ h ω [ Σ γ ( E ) − Σ γ ( E + )] ( γ = r , < ) . (8)Taking advantage of this simplified form, the AC retardedGreen function, g r ω ( E ) , can be determined by substitutingEqs. (6)-(8) into Eq. (4) and only retaining terms that arelinear in the bias amplitude, V AC . The resulting expressionis a product of DC Green functions at energies E and E + = E + ¯ h ω [24]: g r ω ( E ) = G r ( E + ) [ U ω + σ r ω ( E )] G r ( E ) , (9)where U ω is the AC potential energy profile. The AC re-tarded Green function is sufficient to calculate density of H q,q H q,q-1 H q-1,q H q-1,q-1 H q+1,q+1 H H H H N,N H N-1,N-1 ... ... H q-1,q H q,q+1 H N,N-1 H N-1,N H H H H . . .. . . . . .. . . . . . . . .. . . . . .. . . . . . H H H H H H H H q − − H q − H q,q − H q,q H q,q H q H q,q H N − − H N − H N ,N − H N ,N (a)(b) +1,q +1 { { { N b N b N b N Fig. 1 (a) The device region can be considered diagonal matrix slices H q , q connected to the nearest neighbor lattice H q + , q + through off di-agonal matrix blocks H q , q + and H q + , q (b) The tight-binding Hamil-tonian matrix of this layered hopping structure is a block tridiagonalmatrix. The AC recursive Green algorithm exploits this matrix struc-ture to efficiently calculate the AC Green functions. states and, in the wide-band limit when the contact self-energy is assumed to be independent of energy, terminalcurrents, but the AC lesser Green function is necessary tocalculate spatially-resolved particle and current density.The AC lesser Green function gives one direct accessto charge flow within the system by convolving the densityof states information contained in the AC retarded Greenfunction with the occupancy of the leads that is stored inthe lesser self-energy, Σ < ( E ) . In the energy-domain, thisconvolution can be simply written by the Keldysh equation G < ( E ) = G r ( E ) Σ < ( E ) G r ( E ) † [24]. After Eqs. (6) and (7)are substituted into this relation, the expression for the AClesser Green function to first-order in the AC bias amplitudeis given as g < ω ( E ) = G r ( E + ) Σ < ( E + ) g r ω ( E ) † + G r ( E + ) σ < ω ( E ) G r ( E ) † + g r ω ( E ) Σ < ( E ) G r ( E ) † . (10)Again, we see that the AC linear response of the system is aproduct of DC Green functions at energies E and E + .Once the AC retarded Green function and AC lesserGreen function are computed, observables can be calculatedin a fashion similar to DC NEGF [11,12]. The spatially-resolved local density of states (LDOS) is only dependenton the AC retarded Green function and is calculated byLDOS ( r ) = i π (cid:0) g r ω ; r , r ( E ) − g r ω ; r , r ( E ) † (cid:1) . (11) Timothy M. Philip, Matthew J. Gilbert
The frequency-dependent electron density n ω ( r ) is given as n ω ( r ) = − i (cid:90) dE π g < ω ; r , r ( E ) . (12)While the electron density is important to understand chargedynamics, the AC particle current density must be deter-mined to self-consistently calculate the dynamic magneticfield within the system and is given by J α , ω ( r ) = − e ¯ h (cid:90) dE π (cid:104) H r + a ˆ α , r g < ω ; r , r + a ˆ α ( E ) − g < ω ; r + a ˆ α , r ( E ) H r , r + a ˆ α (cid:105) ( α = x , y , z ) . (13)Additionally, the AC terminal currents can be computed as[24] I β , ω = − e ¯ h (cid:90) dE π Tr (cid:2) g < ω ( E ) Σ r β ( E ) † − Σ r β ( E + ) g < ω ( E ) g r ω ( E ) Σ < β ( E ) − Σ < β ( E + ) g r ω ( E ) † G r ( E + ) σ < ωβ ( E ) − σ < ωβ ( E ) G r ( E ) † G < ( E + ) σ r ωβ ( E ) † − σ r ωβ ( E ) G < ( E ) (cid:3) . (14)where β = L , R for the left and right contact currents.2.3 AC Recursive Green Function AlgorithmAlthough the AC NEGF equations described in the previoussection can be directly solved to simulate electron transport,it is a computationally intensive task since two matrix in-versions are required to obtain G r ( E ) and G r ( E + ) at eachstep of the energy integration necessary to calculate the par-ticle density and current density in Eqs. (12) and (13). For amatrix of side length N tot , the number of operations neededto invert it scales as O ( N ) [11]. Since the Hamiltonianside length is a product of the number of lattice sites andthe number of orbitals, the AC NEGF algorithm becomesprohibitively expensive for large systems due to the com-putational burden of inverting large matrices. Calculatingparticle and current density using Eqs. (12) and (13), how-ever, does not require all of the matrix elements of the entireGreen function. In fact, only the main diagonal blocks andthe first off-diagonal blocks of the lesser Green function arenecessary to calculate these observables. With this in mind,we develop a AC recursive Green function (RGF) algorithmfor the AC NEGF equations that efficiently solves for onlythe main diagonal and first off-diagonal of the AC lesserGreen function while avoiding full matrix inversions by ex-ploiting the sparse nature of the Hamiltonian matrix. OurAC RGF algorithm extends the well-known DC RGF algo-rithm that is commonly used in modern DC NEGF transportsimulations to the AC NEGF Eqs. (9) and (10) [9,11]. Since our algorithm requires some parts of the DC RGFalgorithm, we begin by reviewing the RGF algorithm forthe DC NEGF method before extending it to the AC case.Equation (5) for the DC retarded Green function, G r , can berewritten as A G r = I , (15)where A = (( E + i η ) I − H − U − Σ r ) is the inverse ofthe Green function, H is the matrix representation of theHamiltonian, U is the DC potential energy profile, and Σ r is the self-energy that captures the influence of external leadsand and scattering mechanisms. In order to avoid invertingthe entire A matrix to obtain the Green function, we mustconsider the structure of the Hamiltonian. The RGF algo-rithm is most efficiently applied to layered systems, wherethe Hamiltonian matrix parallel to the transport directionconsists of N b orbitals/lattice sites that can be considered asingle layer or slice. Each N b × N b slice of the Hamiltonianis connected only connected to the next slice in the transportdirection though an off-diagonal N b × N b block. Aside fromthe these main diagonal block and first off-diagonal blocks,the rest of the Hamiltonian matrix H is zero. Figure 1schematically depicts the layered system and the resultingHamiltonian matrix. We adopt a matrix notation where H q , s refers to a N b × N b block of the Hamiltonian, H , for blockrow index q and block column index s : H q , s = H ( q − ) N b + qN b , ( s − ) N b + sN b (16)For example, when N b = H , refers to the matrix blockcomprising of rows 21 to 30 and columns 31 to 40. InFig. 1a, we depict the each layer of N b orbitals and latticesites as blue slices labeled by H q , q , where q runs from 1 to N total number of slices. Coupling between the layers inthe transport direction is indicated by the arrows betweenthe slices labeled H q + , q and H q + , q . Figure 1b depicts thematrix representation of this layered structure as a sparse,block tridiagonal matrix with side length N tot = N b N , whereall nonzero blocks are labeled by their block indices. The A matrix in Eq. (15) for these layered structures is, therefore,also sparse and has the same block tridiagonal structure asthe Hamiltonian. The DC RGF algorithm exploits the spar-sity and block tridiagonal form of this matrix to calculate therelevant blocks of the Green function using N inversions of N b × N b and a number of computationally trivial matrix mul-tiplications instead of inverting a single matrix of side length N b N . This reduces the computational complexity of calcu-lating the Green function from O ( N b N ) to O ( N b N ) andprovides significant speedup in obtaining the Green functionfor systems that have many layers [9,11]. The RGF algorithm begins by inverting the first N b × N b block of the A matrix, which is the first block of the left- heory of AC quantum transport with fully electrodynamic coupling 5 connected Green function, g r , L . The name of this Greenfunction derives from the fact that it only contains informa-tion about the left-lead through the contact self-energy thatis part of the first block of A . To obtain the elements of thefull Green function, we must build the rest of g r , L until it isconnected to the right contact self-energy. The blocks of agiven Green function will be subscripted with the same no-tation as the A matrix. The first block of the left-connectedGreen function is given simply as g r , L , = A − , . (17)Subsequent diagonal blocks of g r , L are connected back to thefirst block recursively through the following relationship: g r , Lq , q = (cid:16) A q , q − A q , q − g r , Lq − , q − A q − , q (cid:17) − . (18)Thus, g r , L is constructed by started on the top left of the A matrix and iterating over q until q = N , where the right con-tact is located. Since this last block, g r , LN , N , includes the influ-ence of the contact self energy at both q = q = N , it isexactly the last diagonal block of the fully-connected Greenfunction G rN , N .We also build an analogous right-connected Green func-tion that begins at q = N and is connected only to the rightcontact though the contact self-energy that is present in A N , N . The q = N block is given as g r , RN , N = A − N , N . (19)The blocks of the right-connected Green function are foundby the recursive relationship g r , Rq , q = (cid:16) A q , q − A q , q + g r , Rq + , q + A q + , q (cid:17) − . (20)Therefore, the blocks of g r , R are computed by iterating back-wards from q = N to q =
1. Again, we note that similar tothe last block of g r , L being the last block of the full Greenfunction, the q = g r , R , , is exactly G r , . In the DC NEGF formulation, one typ-ically only needs to compute one of either g r , L or g r , R to ob-tain the desired elements of the full Green function, but forthe AC extension of the recursive algorithm described in thefollowing sections, we will need both to obtain the leftmostand rightmost columns of the full Green function. Althoughthese columns can be calculated using just g r , L and matrixmultiplication, the large number of off-diagonal componentsof the Green that must be computed makes the approachcomputationally inefficient. It is, therefore, more efficient tocalculate both g r , L and g r , R to obtain the desired columns.Once the left-connected and right-connected Greenfunctions have been constructed, any off-diagonal block of the full G r matrix can be computed using multiplicative re-lationships. Elements above the main diagonal of the Greenfunction are recursively related by G rq , s (cid:12)(cid:12) q < s = − g r , Lq , q A q , q + G rq + , s = − G rs , q − A q − , q g r , Rq , q , (21)and elements below the main diagonal follow the relation-ship G rq , s (cid:12)(cid:12) q > s = − G rs , q + A q + , q g r , Lq , q = − g r , Rq , q A q , q − G rq − , s . (22)The AC RGF implementation requires the left and rightcolumns of G r , so we use Eqs. (21)-(22) to obtain them bymultiplying away from the main diagonal blocks: G r (rgt) q , N = − g r , Lq , q A q , q + G r (rgt) q + , N , (23) G r (lft) q , = − g r , Rq , q A q , q − G r (lft) q − , . (24)For clarity, we add the superscript (rgt) to the elements ofthe rightmost column of the Green function and (lft) to theelements of the leftmost column. Once the left and rightcolumns of the DC Green function are obtained, the ACGreen function can be computed. For the AC RGF algorithm described below, we assume co-herent transport such that AC self-energies, σ r ω and σ < ω , areonly nonzero in the first ( q =
1) and last ( q = N ) diago-nal sub-blocks due to the presence of contacts. Under suchan assumption, we show in the proceeding section that thelesser AC Green function, g < ω ; q , q , blocks can be calculatedusing simple matrix multiplication using just the left andright block columns of G r , G r + , and g r ω .The AC NEGF method involves calculating the prod-uct of Green functions at two different energies, previouslydenoted as G r ( E ) and G r ( E + ¯ h ω ) . As such, there will bea number left-connected and right-connected Green func-tions associated with each full matrix. To simplify notation,we will denote functions associated with G r ( E ) with a sub-script 0 and those with G r ( E + ¯ h ω ) with subscript +. Forexample, g r , L refers to the left-connected Green function of G r ( E ) = G r , while g r , L + refers to that of G r ( E + ¯ h ω ) = G r + .The AC Green function will continue to be subscripted with ω , so, g r , L ω refers to left-connected Green function of the ACretarded Green function, g r ω .We can rewrite Eq. (9) in this abbreviated notation as g r ω = G r + ( σ r ω + U ω ) G r , (25) A + g r ω = ( σ r ω + U ω ) G r , (26)where A + = (( E + ¯ h ω + i η ) I − H − Σ ( E + ¯ h ω )) . Equa-tion (26) resembles the equation of the DC lesser Green Timothy M. Philip, Matthew J. Gilbert function G < [10], and so, we develop an RGF algorithm for g r ω , similar to the recursive algorithm used to calculate theDC lesser Green function [11]. To obtain g r ω , we first usethe recursive algorithm from the previous section to obtainthe left- and right-connected Green functions for both G r and G r + in addition to the left and right columns of their fullmatrices.To calculate the relevant blocks of the AC retardedGreen function, we first must obtain the diagonal blocks ofthe left-connected AC retarded Green function, g r , L ω ; q , q . Thefirst, q =
1, block of g r , L ω ; q , q can be trivially calculated fromthe left-connected Green functions g r , L and g r , L + : g r , L ω ;1 , = g r , L + ;1 , ( U ω ;1 , + σ r ω ;1 , ) g r , L , . (27)The subsequent diagonal blocks of of the g r , L ω can be com-puted through the multiplicative, recursive relationship g r , L ω ; q , q = g r , L + ; q , q (cid:2) U ω ; q , q + σ r ω ; q , q + A + ; q , q − g r , L ω ; q − , q − A + ; q − , q (cid:105) g r , L q , q . (28)In general, for any pair of block indices q and s where q (cid:54) = s , A + ; q , s = A q , s , so no additional matrix elements must becomputed to evaluate Eq. (28). We simply retain the sub-scripted + to keep consistent notation with Eq. (26). We alsocompute the right-connected AC retarded Green function, g r , R ω , via analogous expressions: g r , R ω ; N , N = g r , R + ; N , N [ U ω ; N , N + σ r ω ; N , N ] g r , R N , N (29) g r , R ω ; q , q = g r , R + ; q , q [ U ω ; q , q + σ r ω ; q , q + A + ; q , q + g r , R ω ; q + , q + A + ; q + , q ] g r , R q , q . (30)To calculate the lesser AC Green function, we need theleftmost and rightmost columns of g r ω ( E ) . We can obtainthem using the minimum number of matrix multiplicationthrough the following relations: g r (rgt) ω ; q , N = − g r , L ω ; q , q A + q , q + G r (rgt)0; q + , N − g r , L + ; q , q A + q , q + g r (rgt) ω ; q + , N , (31) g r (lft) ω ; q , = − g r , R ω ; q , q A + q , q − G r (lft)0; q − , − g r , R + ; q , q A + q , q − g r (lft) ω ; q − , . (32)Again, we add the superscript (rgt) to the elements of theright column of the Green function and (lft) to the elementsof the left column. Using this methodology, we simplifythe AC Green function computation from two inversions of N b N × N b N matrices to 2 N inversions of N b × N b matricesalong with a handful of computationally trivial matrix multi-plications. This reduces the computational complexity from O ( N b N ) to O ( N b N ) , which can greatly increase computa-tional speed. Once the leftmost and rightmost columns of the AC retardedGreen function are obtained using Eq. (31)-(32), we can cal-culate the AC lesser Green function, which is then used cal-culate relevant observables. In the abbreviated notation wehave established, Eq. (10) for the AC lesser Green functionis written as g < ω = G r + Σ < + g r † ω + G r + σ < ω G r †0 + g r ω Σ < G r †0 , (33) = g < ω + g < ω + g < ω . (34)In the second line, we label each term in the first line witha superscripted number. It is more convenient to calculateeach of these terms individually and then add them togetherto obtain the total AC lesser Green function.If the two-terminal transport is coherent, the self-energies Σ < / + / σ < ω are zero except for the q = q = N block diagonal due to the presence of con-tact self-energies. With such simplification, the three termsof g < ω can be calculated using block matrix multiplicationinvolving the leftmost and rightmost columns of G r / + / g r ω .In the block notation of the RGF algorithm, the three termsof g < ω are calculated as g < ω ; q , q = G r (lft) + ; q , Σ < + ;1 , (cid:16) g r (lft) ω ; q , (cid:17) † + G r (rgt) + ; q , N Σ < + ; N , N (cid:16) g r (rgt) ω ; q , N (cid:17) † , (35) g < ω ; q , q = G r (lft) + ; q , σ < ω ;1 , (cid:16) G r (lft)0; q , (cid:17) † + G r (rgt) + ; q , N σ < ω ; N , N (cid:16) G r (rgt)0; q , N (cid:17) † , (36) g < ω ; q , q = g r (lft) ω ; q , Σ < , (cid:16) G r (lft)0; q , (cid:17) † + g r (rgt) ω ; q , N Σ < N , N (cid:16) G r (rgt)0; q , N (cid:17) † , (37)Which when added together to obtain g < ω ; q , q , allows us tocompute the particle density using Eq. (12). To compute cur-rent density using Eq. (13), we need the off-diagonal ele-ments g < ω ; q , q + and g < ω ; q + , q . The off-diagonal componentsof g < ω are computed in a similar fashion to the diagonal ele-ments: g < ω ; q + , q = G r (lft) + ; q + , Σ < + ;1 , (cid:16) g r (lft) ω ; q , (cid:17) † + G r (rgt) + ; q + , N Σ < + ; N , N (cid:16) g r (rgt) ω ; q , N (cid:17) † , (38) g < ω ; q + , q = G r (lft) + ; q + , σ < ω ;1 , (cid:16) G r (lft)0; q , (cid:17) † + G r (rgt) + ; q + , N σ < ω ; N , N (cid:16) G r (rgt)0; q , N (cid:17) † , (39) g < ω ; q + , q = g r (lft) ω ; q + , Σ < , (cid:16) G r (lft)0; q , (cid:17) † + g r (rgt) ω ; q + , N Σ < N , N (cid:16) G r (rgt)0; q , N (cid:17) † , (40) heory of AC quantum transport with fully electrodynamic coupling 7 g < ω ; q , q + = G r (lft) + ; q , Σ < + ;1 , (cid:16) g r (lft) ω ; q + , (cid:17) † + G r (rgt) + ; q , N Σ < + ; N , N (cid:16) g r (rgt) ω ; q + , N (cid:17) † , (41) g < ω ; q , q + = G r (lft) + ; q , σ < ω ;1 , (cid:16) G r (lft)0; q + , (cid:17) † + G r (rgt) + ; q , N σ < ω ; N , N (cid:16) G r (rgt)0; q + , N (cid:17) † , (42) g < ω ; q , q + = g r (lft) ω ; q , Σ < , (cid:16) G r (lft)0; q + , (cid:17) † + g r (rgt) ω ; q , N Σ < N , N (cid:16) G r (rgt)0; q + , N (cid:17) † , (43)If transport on incoherent or other blocks of the self-energiesare nonzero, Eqs. (35)-(43) are no longer valid. Instead, g < ω must be constructed from left-connected and right-connected Green functions for each term g < / / ω follow-ing the recursive framework used to calculate the elementsof g r ω , as presented in the previous section. Evaluating theabove equations to obtain the elements of g < ω only requiresmatrix multiplications, which are significantly less compu-tationally demanding than matrix inversions. As such thecomputational complexity for this AC NEGF algorithm toobtain g < ω remains O ( N b N ) . In summary, the RGF procedure to obtain the relevantblocks of the AC retarded and lesser Green functions withinthe limit of coherent, two-terminal transport can be calcu-lated with the following steps:1. Evaluate the blocks of left-connected DC Green functionusing Eq. (17) at E , the energy of interest, to obtain g r , L , and at E + = E + ¯ h ω to obtain g r , L + ;1 , .2. Using the recursive relationship for the left-connectedDC Green function defined in Eq. (18), construct theremaining blocks of g r , L and g r , L + by iterating q = , , . . . N − , N .3. Evaluate the blocks of the right-connected DC Greenfunction using Eq. (19) at E to obtain g r , R N , N and at E + to obtain g r , R + ; N , N .4. To obtain the the remaining blocks of the right-connected DC Green functions g r , R and g r , R + , iterate over q = N − , N − , . . . , , G r and G r + using Eq. (23).5. Now that all the blocks of the left- and right-connectedDC Green functions g r , L , g r , L + , g r , R and g r , R + , have been ob-tained, the AC left- and right-connected Green functions, g r , L ω and g r , R ω , can be computed. Start by using Eq. (27) tocalculate the first block of the left-connected AC Greenfunction, g r , L ω ;1 , . xz y ρ /V { J x /A x ;E x B y a J z A z ;E z J y A y ;E y B x B z Fig. 2
The FDFD equations are solved on a cubic Yee cell with sidelength a , which matches the lattice constant of the tight-binding Hamil-tonian. The tight-binding Hamiltonian lattice where the charge density, ρ , is located by the solid squares at the corners of the grid; the poten-tial, V , is naturally defined at the same positions. The current density, J α ( α = x , y , z ) , is located between lattice sites in the ˆ α direction, sothe same component of the vector potential A α is also at these po-sitions. The definition of electric field, E = − ∇ φ + i ω A , naturally co-locates E α at the same sites. The definition of the magnetic flux density, B = ∇ × A , places its components at offset positions from the electricfield.
6. Build the remainder of the left-connected AC Greenfunction, g r , L ω , by iterating Eq. (28) over q = , , . . . , N − , N . It is also convenient to calculate the left columns ofthe fully-connected DC Green functions G r and G r + us-ing Eq. (24) on this iterative loop.7. Next, begin calculating the blocks of the right-connectedAC Green function g r , R ω by using Eq. (29) to obtain g r , R ω ;1 , .8. Construct the remainder of the right-connected ACGreen function g r , R ω by iterating Eq. (30) over q = N − , N − , . . . , ,
1. The right column of the fully-connected AC Green function g r ω can also be calculatedduring the loop using Eq. (31).9. Finally, one last iterative loop over q = , , . . . , N − , N is necessary to calculate the the leftmost column of thefully-connected AC Green function g r , R ω using Eq. (31).10. Now that the leftmost and rightmost columns of DCGreen function at E , G r , the DC Green function at E + , G r + , and the AC Green function, g r ω , have been calcu-lated, Eqs. (35)-(43) are used to obtain the main diago-nal and first off-diagonal blocks of the AC lesser Greenfunction, g < ω .This procedure must be repeated over a discretized en-ergy range to compute particle density using Eq. (12) withthe main diagonal of g < ω and the current density usingEq. (13) with the first off-diagonal of g < ω . Timothy M. Philip, Matthew J. Gilbert
Charge transport through a system is strongly influenced byelectric and magnetic fields, and thus must be accounted forwhen solving the AC NEGF equations. In DC NEGF, theself-consistent electrostatic potential is sufficient to capturethe the effect of the electromagnetic environment. At fre-quencies when the ratio of the frequency to the speed oflight, ω / c , is not negligible, however, this electrostatic as-sumption is not adequate to account for the dynamic chargemotion, so the full solution of Maxwell’s equation must besolved to full understand the influence of electrodynamiccoupling.Standard treatments for solving Maxwell’s equationscalculate the electric field, E , and magnetic field, B [33],but the quantum mechanical wave function is dependent,however, on the scalar potential V and vector potential A [34]. We, therefore, reformulate the finite-differencefrequency-domain (FDFD) method [35,36] to solve directlyfor the electromagnetic potentials. In the frequency-domain,Maxwell’s equations in the Lorenz gauge, where ∇ · A = − i ω c V , take the form (cid:18) ∇ + ω c (cid:19) V ω = − ρ ω ε , (44) (cid:18) ∇ + ω c (cid:19) A ω = − µ J ω , (45)where ω is the frequency of interest, c is the speed of light, ε is the electric permittivity, and µ is the magnetic perme-ability. The charge density ρ ω and current density J ω are ex-tracted from AC NEGF using Eqs. (12) and (13). The elec-tric field and magnetic flux density components, E α and B α ( α = x , y , z ) , respectively, are numerically calculated fromthe potentials using the frequency domain relations: E ω = − ∇ V ω + i ω A ω , (46) B ω = ∇ × A ω . (47) Solving for the scalar and vector potential via Eqs. (44) and(45) on a discrete lattice requires care because the chargedensity and current density components from AC NEGF arenot co-located on the same grid. The charge density is lo-cated at the lattice sites of the Hamiltonian, but the currentdensity components, calculated from the particle flux trans-fered between lattice sites, is located between lattice sites.Therefore, Eqs. (44) and (45) must be solved on the stag-gered Yee cell illustrated in Fig. 2, where the scalar potentialis defined at the same sites as the tight-binding lattice, whilethe components of the vector potential A α ( α = x , y , z ) are offset in position from the lattice to be co-located with thecurrent density components J α . As a result, when Eqs. (46)and (47) are evaluated using these staggered potentials, theelectric and magnetic field naturally are defined locationsas would be expected in the standard finite-difference time-domain (FDTD) Yee cell discretization [35].Equations (44) and (45) both describe Helmholtz equa-tions, which in rectangular coordinates can be genericallywritten (cid:18) ∂ ∂ x + ∂ ∂ y + ∂ ∂ z + ω c (cid:19) u ( r ) = − f ( r ) . (48)To solve for the scalar potential, f ( r ) is replaced with ρ ω ( r ) / ε and u ( r ) with V ω ( r ) . Similar replacements can bemade for the different components of the vector potential.We discretize the Laplacian in Eq. (48) in rectangular coor-dinates using a second-order, central finite difference: ∂ u ∂ x ≈ u i + , j , k − u i , j , k + u i − , j , k ∆ x , (49)where we adopt the notation u i , j , k = u ( i ∆ x , j ∆ y , k ∆ z ) . Foruniform, cubic grid-spacing where ∆ x = ∆ y = ∆ z , the dis-crete form of the Helmholtz equation is written u i + , j , k + u i − , j , k + u i , j + , k + u i , j − , k + u i , j , k + + u i , j , k − + (cid:0) r ω − (cid:1) u i , j , k = − ∆ x f i , j , k , (50)where r ω = ∆ x ω / c and the ∆ x is equal to the lattice con-stant, a , of the lattice Hamiltonian. To minimize discretiza-tion error, the condition r ω < n × A ω ( r ) = V ω ( r ) = n is the unit normal to the surface of the bound-ary. The PEC boundary conditions perfectly reflect any in-cident electromagnetic waves, which accurately capturesthe behavior of highly conductive surfaces. When PECboundary conditions are applied to all faces of the simula-tion domain, however, an electromagnetic cavity is createdwhereby waves with wavelengths commensurate with thedomain side length are enhanced and others are attenuated. heory of AC quantum transport with fully electrodynamic coupling 9 Such behavior must be avoided when modeling field radia-tion, as these cavity modes will significantly alter the radia-tion profile. As such, radiative boundary conditions that donot reflect incident waves must be applied for any bound-ary that is not explicitly metallic. To this end, stretched-coordinate, convolutional perfectly-matched-layers are at-tached to all non-metallic boundaries to absorb any radi-ating fields and inhibit the development of artificial cavitymodes [38,39].
Once Eqs. (44) and (45) are solved using the charge andcurrent density obtained from the AC NEGF calculations,the output scalar and vector potentials must be reinsertedinto the transport simulation. The scalar potential enters asan on-site potential which modifies the on-site term in theHamiltonian described in Eq. 1: H → H − eV ( r ) (53)The vector potential is coupled to the tight-binding Hamil-tonian through a Peierl’s phase on the off-diagonal hoppingterms [40]: H δ → H δ exp (cid:18) − ie ¯ h (cid:90) r + δ r A ( r ) · d (cid:96) (cid:19) . (54)Based on the Yee cell discretization of the FDFD equationsdescribed in the previous section, the components of the vec-tor potential are piece-wise constant between lattice sites.Therefore, the modified hopping amplitude is simplified to H α → H α exp (cid:18) − ie ∆ α ¯ h A α ( r + ∆ α / ˆ α ) (cid:19) , (55)where α = { x , y , z } . After these substitutions are made, theAC NEGF and FDFD equations must be solved again re-peatedly until their solutions stabilize and self-consistency isattained. We mark self-consistency in our simulations whenthe change in the scalar potential between successive it-erations is less than 1 µ V. Reaching self-consistency withthis iterative method can fail to converge and is often time-consuming, so we utilize the Anderson mixing scheme tostabilize and accelerate self-consistency [41]. L is mounted Ground Plane (a) (b) (c) (d)
Feed Point
Fig. 3
Schematic of the device design and the operation of a classi-cal monopole antenna. (a) A monopole antenna comprises of a currentcarrying wire of length L that is extended above a conductive groundplane. (b) The image charge and currents in the ground plane resultsin radiation from the monopole antenna that mimics that from dipoleantenna of length 2 L . (c) A short antenna is one in which the antennalength is much less than wavelength of operation ( L (cid:28) λ ), which gen-erates a current profile that linearly drops from the feed point to theend of the antenna. When the antenna length L = λ /
4, the currentdistribution forms a sinusoidal, standing wave current profile. (d) Thesinusoidal current profile of the quarter-wave monopole creates in amore directed classical radiation pattern than that of a short monopole. on a conductive ground plane, which reflects the fields radi-ated from the antenna. The reflected fields from the groundplane can be modeled as image charges and currents belowthe monopole that result in radiation that is equivalent to adipole antenna of length 2 L , as depicted in Fig. 3(b).Figure 3(c) illustrates the classically expected currentdensity in an electrically short monopole antenna, that iswhere L (cid:28) λ and λ is the operating wavelength, and aquarter-wave monopole antenna, where L = λ /
4. In a shortmonopole, the current drops linearly from the feed point tothe end of the antenna with I short ( x ) = iI ( L − x ) . When theantenna is extended to a length L = λ /
4, however, the an-tenna operates at its resonant frequency and the current dis-tribution forms a sinusoidal standing wave along the antennawith the form I λ / ( x ) = iI cos π x λ . Figure 3(d) shows theexpected radiation pattern for both a short and quarter-wavemonopole. We see that these altered current profiles alongthe length of the antenna result in different classical radi-ation patterns. The resultant far-field electric field for theshort antenna has the form [42] | E short | = I ε cr L λ sin θ , (56)where θ is the angle measured from the x axis. The angu-lar dependence of radiation follows a simple sin θ relation- Absorbing EM Boundaries Perfect Electric Conductor
Closed
NEGF Boundary Open
NEGF
Boundary
Fig. 4
Schematic of simulated quantum monopole antenna. The ACNEGF technique is used to calculate transport for an end-fed 1D metalwire with L = µ m. Transport is coupled self-consistently with theelectrodynamic scalar and vector potentials in an FDFD domain withside length (cid:96) = .
73 mm. Perfect electrical conductor electromagnetic(EM) boundary conditions, indicated with the solid gray face, are ap-plied to the bottom y - z plane while absorbing EM boundary conditions,indicated by the dashed edges, are applied on the other faces to inhibitthe development of cavity modes. ship that creates the circular lobes in the radiation pattern de-picted in Fig. 3(d). When the antenna is driven at the quarter-wave frequency, the far-field electric field profile becomes | E λ / | = I πε cr cos ( π / θ ) sin θ . (57)This more complicated angular dependence results in lesspower delivered at small angles and more at angles closeto θ = π , which is illustrated by the oval shaped lobes ofthe λ / θ = π .3.2 Quantum CaseHaving reviewed the expected radiation patterns from a clas-sical short and quarter-wave monopole antenna, we nowinvestigate the radiation behavior of a monopole antennathat possesses quantized energy states, which we will call aquantum monopole antenna. We use the self-consistent ACNEGF/FDFD technique described here to model the radi-ation of the quantum monopole antenna shown in Fig. 4.The antenna is modeled by a 20 lattice site one-dimensionalmetal Hamiltonian H ( r ) = ∑ i t (cid:2) ψ † x ψ x − (cid:0) ψ † x ψ x + a + H.c. (cid:1)(cid:3) , (58)with hopping amplitude t = . L = µ m [10]. Although the length of this antenna is relatively Fig. 5
The calculated local density of states of the quantum wire an-tenna reveals quantized energy states at µ =
34 meV and µ = large and seemingly classical, the hopping parameter valuewe use creates quantized states within the wire. Since cur-rent is only injected at one end of the monopole antenna,only one end of the antenna has an open NEGF boundarycondition, while the other end has a closed boundary con-dition. The quantum antenna is placed in a cubic electro-magnetics domain with side length (cid:96) = .
73 mm, with per-fect electric conductor boundary conditions applied to the y - z plane to provide the ground plane needed for operation ofthe monopole antenna. Absorbing electromagnetic bound-ary conditions are applied to all other faces of the electro-dynamics domain to allow for field radiation away from theantenna. With such a small electrodynamics simulation do-main, the field surrounding the antenna is largely dominatedby the non-radiative near field. To understand the far-fieldradiation pattern, we perform a near-to-far-field transforma-tion on the self-consistent electrodynamic fields around theantenna [43]. This technique allows us to understand themacroscopic radiative characteristics of the antenna withoutsimulating a large electromagnetic domain outside the nearfield. We drive the antenna at the quarter-wave frequency of100 GHz to understand the differences between the classicalmodel of a monopole antenna using our quantum coherentAC NEGF/FDFD simulation methodology.Figure 5 shows the computed LDOS of the quantum an-tenna. The energies of the states within the one-dimensionalwire match those of an infinite square well within 5 meV,thereby demonstrating the quantum nature of transportwithin the antenna. We self-consistently calculate transportthrough the first and second quantized state by setting thechemical potential µ =
34 meV and µ =
140 meV, re-spectively. In Fig. 6(a), we see that the AC quantum chargedensity differs significantly from the classically expectedcharge distribution. Rather than being maximized at the endof the antenna as is anticipated in the classical charge dis- heory of AC quantum transport with fully electrodynamic coupling 11 -90-60-300306090
Quantum = Quantum = Classical x ( m) R e [ n ] ( c m - ) Quantum = Quantum = Classical x ( m) I m [ J x ] ( m A / c m ) Quantum = Quantum = Classical (a)(b)(c)
Fig. 6
Simulation results of a quantum quarter-wave antenna throughthe first and second quantized energy level at the chemical potentials µ =
34 meV and µ =
140 meV, respectively. (a) The charge densitydistribution along the length of the antenna differs from the classicalexpectation due to the wave function of the quantized states. (b) Thequantum confinement also alters the current density and does not cre-ate the expected sinusoidal profile. (c) The quantization effects on thecharge and current density result in a modification to the macroscopicradiation pattern. tribution, the wave function of the quantized electronic statemaximizes the the charge distribution where the anti-nodesof the wave functions are seen in the LDOS. Because ofthe current-conserving nature of the AC NEGF/FDFD tech-nique, the AC current density is intimately related to chargedensity via the continuity equation, and it too is altered fromthe expected classical arrangement. Despite being driven onthe quarter-wave, resonant frequency, the current distribu-tion no longer forms the standing wave along the lengthof the antenna. These non-ideal charge and current densi-ties due to the spatial form of the quantum-confined wavefunction result in the distorted macroscopic radiation pat-tern in Fig. 6(c). Therefore, instead of having the directivityof gain associated with a classical quarter-wave monopole,the quantum quarter-wave monopole has little to no gain at either chemical potential and radiates identically as a shortantenna.
In this work, we have detailed a novel quantum transportsimulation methodology that goes beyond the quasi-staticapproximation by coupling AC NEGF with the full, dynamicsolution of Maxwell’s equations. This unique technique al-lows us to explore electrodynamic phenomena that cannotbe understood via traditional DC or transient methods. Wedemonstrate the utility of this formulation by modeling theradiation from a quantum, quarter-wave, monopole antenna.We find that the quantum-confined wave function withinthe wire constrains the charge and current density profiles,which disallows the formation of the standing wave alongthe length of the wire despite being driven on resonance.Because of the altered current profile within the antenna, weobserve that the macroscopic radiation pattern is altered bythe quantization, resulting in no directivity gain compared toa short antenna. Our results illustrate that this new simula-tion methodology can uncover unexpected behavior in high-frequency quantum mechanical systems and will enable thecharacterization of the of next-generation, nanoscale, high-frequency devices, where quantum effects dominate trans-port phenomena.
The authors acknowledge support from the NSF under CA-REER Award ECCS-1351871. T.M.P. thanks M.J. Park andY. Kim for helpful discussions.
References
1. G. Klimeck, S.S. Ahmed, N. Kharche, M. Korkusinski, M. Usman,M. Prada, T.B. Boykin, IEEE Trans. Electron Devices (9), 2090(2007). DOI 10.1109/TED.2007.9048772. L. Huang, Y.C. Lai, D.K. Ferry, R. Akis, S.M. Good-nick, J. Phys. Condens. Matter (34), 344203 (2009).DOI 10.1088/0953-8984/21/34/344203. URL http://stacks.iop.org/0953-8984/21/i=34/a=344203?key=crossref.06e69231502c0dec272d7605e1646d29
3. K. Balzer, M. Bonitz, R. van Leeuwen, A. Stan, N.E. Dahlen,Phys. Rev. B (24), 245306 (2009). DOI 10.1103/PhysRevB.79.245306. URL http://link.aps.org/doi/10.1103/PhysRevB.79.245306
4. R. Lake, S. Datta, Phys. Rev. B (12), 6670 (1992). DOI 10.1103/PhysRevB.45.66705. V.N. Do, P. Dollfus, V. Lien Nguyen, J. Appl. Phys. (9) (2006).DOI 10.1063/1.23640356. M. Luisier, A. Schenk, W. Fichtner, J. Appl. Phys. (4) (2006).DOI 10.1063/1.22445227. G. Fiori, G. Iannaccone, IEEE Electron Device Lett. (8), 760(2007). DOI 10.1109/LED.2007.9016802 Timothy M. Philip, Matthew J. Gilbert8. S.O. Koswatta, S. Hasan, M.S. Lundstrom, M.P. Anantram, D.E.Nikonov, IEEE Trans. Electron Devices (9), 2339 (2007). DOI10.1109/TED.2007.9029009. R. Lake, G. Klimeck, R.C. Bowen, D. Jovanovic, J. Appl.Phys. (12), 7845 (1997). DOI 10.1063/1.365394. URL http://scitation.aip.org/content/aip/journal/jap/81/12/10.1063/1.365394
10. S. Datta, Superlattices Microstruct. (4), 253 (2000). DOI10.1006/spmi.2000.0920. URL http://linkinghub.elsevier.com/retrieve/pii/S0749603600909200
11. M. Anantram, M. Lundstrom, D. Nikonov, Proc. IEEE (9),1511 (2008). DOI 10.1109/JPROC.2008.927355. URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=4618725
12. M. Pourfath,
The Non-Equilibrium Green’s Function Methodfor Nanoscale Device Simulation . Computational Micro-electronics (Springer Vienna, Vienna, 2014). DOI 10.1007/978-3-7091-1800-9. URL http://link.springer.com/10.1007/978-3-7091-1800-9
13. R.O. Grondin, S.M. El-Ghazaly, S. Goodnick, IEEE Trans. Mi-crow. Theory Tech. (6 PART 2), 817 (1999). DOI 10.1109/22.76931514. A. Witzig, C. Schuster, S. Member, P. Regli, W. Fichtner, (6),919 (1999)15. M. Sirbu, S.B. Lepaul, F. Aniel, IEEE Trans. Microw. TheoryTech. (9), 2991 (2005). DOI 10.1109/TMTT.2005.85422816. K.J. Willis, S.C. Hagness, I. Knezevic, Appl. Phys. Lett. (6)(2010). DOI 10.1063/1.330849117. K.J. Willis, S.C. Hagness, I. Knezevic, J. Appl. Phys. (6)(2011). DOI 10.1063/1.362714518. E. Perfetto, G. Stefanucci, M. Cini, Phys. Rev. B (3), 035446(2010). DOI 10.1103/PhysRevB.82.035446. URL http://link.aps.org/doi/10.1103/PhysRevB.82.035446
19. B. Gaury, X. Waintal, Nat. Commun. (May), 1 (2014).DOI 10.1038/ncomms4844. URL
20. M.W.Y. Tu, A. Aharony, W.M. Zhang, O. Entin-Wohlman, Phys.Rev. B - Condens. Matter Mater. Phys. (16), 1 (2014). DOI10.1103/PhysRevB.90.16542221. M.W.Y. Tu, A. Aharony, O. Entin-Wohlman, A. Schiller, W.M.Zhang, Phys. Rev. B - Condens. Matter Mater. Phys. (12)(2016). DOI 10.1103/PhysRevB.93.12543722. P. My¨oh¨anen, A. Stan, G. Stefanucci, R. van Leeuwen, EPL (Eu-rophysics Lett. (6), 67001 (2008). DOI 10.1209/0295-5075/84/67001. URL http://iopscience.iop.org/0295-5075/84/6/67001http://stacks.iop.org/0295-5075/84/i=6/a=67001?key=crossref.7f76c45bb53d32a8d0f63a9c53ac212b
23. M. Puig von Friesen, C. Verdozzi, C.O. Almbladh, Phys. Rev. B (15), 155108 (2010). DOI 10.1103/PhysRevB.82.155108. URL http://link.aps.org/doi/10.1103/PhysRevB.82.155108
24. Y. Wei, J. Wang, Phys. Rev. B (19), 195315 (2009). DOI 10.1103/PhysRevB.79.195315. URL http://link.aps.org/doi/10.1103/PhysRevB.79.195315
25. D. Kienle, M. Vaidyanathan, F. L´eonard, Phys. Rev. B (11),115455 (2010). DOI 10.1103/PhysRevB.81.115455. URL http://link.aps.org/doi/10.1103/PhysRevB.81.115455
26. O. Shevtsov, X. Waintal, Phys. Rev. B (8), 085304 (2013).DOI 10.1103/PhysRevB.87.085304. URL http://link.aps.org/doi/10.1103/PhysRevB.87.085304
27. J.Q. Zhang, Z.Y. Yin, X. Zheng, C.Y. Yam, G.H. Chen, Eur. Phys.J. B (10) (2013). DOI 10.1140/epjb/e2013-40325-728. J. Larsson, Am. J. Phys. (3), 230 (2007). DOI 10.1119/1.2397095. URL http://scitation.aip.org/content/aapt/journal/ajp/75/3/10.1119/1.2397095
29. L.P. Kadanoff, G. Baym,
Quantum Statistical Mechanics (W. A.Benjamin, New York, 1962)30. G. Stefanucci, R. van Leeuwen,
Nonequilibrium Many-Body The-ory of Quantum Systems: A Modern Introduction (Cambridge Uni-versity Press, 2013)31. M.P.L. Sancho, J.M.L. Sancho, J.M.L. Sancho, J. Rubio, J. Phys. FMet. Phys. (4), 851 (2000). DOI 10.1088/0305-4608/15/4/00932. M.P.L. Sancho, J.M.L. Sancho, J. Rubio, J. Phys. F Met. Phys. (5), 1205 (2000). DOI 10.1088/0305-4608/14/5/01633. A. Taflove, S. Hagness, Computational Electrodynmics (ArtechHouse, 2005)34. W.C. Chew, Prog. Electromagn. Res. (August), 69 (2014).DOI 10.2528/PIER14060904. URL
35. Kane Yee, IEEE Trans. Antennas Propag. (3), 302 (1966). DOI10.1109/TAP.1966.1138693. URL http://ieeexplore.ieee.org/document/1138693/
36. R. Luebbers, F.R. Hunsberger, K.S. Kunz, R.B. Standler,M. Schneider, IEEE Trans. Electromagn. Compat. (3), 222(1990). DOI 10.1109/15.5711637. Y. Wong, G. Li, Int. J. Numer. Anal. Model. Ser. B (1), 91 (2011)38. Y. Rickard, N. Georgieva, Wei-Ping Huang, IEEE Microw.Wirel. Components Lett. (5), 181 (2002). DOI 10.1109/7260.1000196. URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1000196
39. W. Shin, S. Fan, J. Comput. Phys. (8), 3406 (2012). DOI10.1016/j.jcp.2012.01.013. URL http://dx.doi.org/10.1016/j.jcp.2012.01.013http://linkinghub.elsevier.com/retrieve/pii/S0021999112000344
40. M. Graf, P. Vogl, Phys. Rev. B (8), 4940 (1995). DOI10.1103/PhysRevB.51.4940. URL http://journals.aps.org/prb/abstract/10.1103/PhysRevB.51.4940{%}5Cnhttp://link.aps.org/doi/10.1103/PhysRevB.51.4940http://link.aps.org/doi/10.1103/PhysRevB.51.4940
41. D.R. Bowler, M.J. Gillan, (July), 473 (2000). DOI10.1016/S0009-2614(00)00750-8. URL http://arxiv.org/abs/cond-mat/0005521{%}0Ahttp://dx.doi.org/10.1016/S0009-2614(00)00750-8
42. W.H. Hayt, J.A. Buck,
Engineering Electromagnetics (McGraw-Hill New York, 2001)43. P. Petre, T.K. Sarkar, IEEE Trans. Antennas Propag.40