Density matrix approach to photon-assisted tunneling in the transfer Hamiltonian formalism
DDensity matrix approach to photon-assisted tunnelingin the transfer Hamiltonian formalism
Paul S. Davids ∗ and Joshua Shank Sandia National Laboratories, Albuquerque, NM 87185, USA ∗ (Dated: July 16, 2018)The transfer Hamiltonian tunneling current is derived in a time-dependent density matrix for-mulation and is used to examine photon-assisted tunneling. Bardeen’s tunneling expression arisesas the result of first order perturbation theory in a mean-field expansion of the density matrix.Photon-assisted tunneling from confined electromagnetic fields in the forbidden tunnel barrier re-gion occurs due to time-varying polarization and wavefunction overlap in the gap which leads to anon-zero tunneling current in asymmetric device structures, even in an unbiased state. The photonenergy is seen to act as an effective temperature dependent bias in a uniform barrier asymmetrictunneling example problem. Higher order terms in the density matrix expansion give rise to multi-photon enhanced tunneling currents that can be considered an extension of non-linear optics wherethe non-linear conductance plays a similar role as the non-linear susceptibilities in the continuityequations. I. INTRODUCTION
Tunneling is a quantum phenomena that has been ex-tensively studied mainly using semiclassical approachesbased on the Wentzel Kramers Brillioun (WKB)approximation[1, 2]. Bardeen introduced a many par-ticle transfer Hamiltonian [3] approach to examine tun-neling currents between superconductors and normalmetals through an insulating gap. The model treatedthe tunneling heuristically to match the observation ofgaps in quasiparticle tunneling seen in the superconduc-tor to superconductor[4] and superconductor to normalmetal[5] tunneling current voltage data of Giaever. Amore detailed transfer Hamiltonian approach based onthe second quantized operator method and canonicaltransformations of the operators was subsequently devel-oped [6, 7]. The extension of Bardeen’s transfer Hamil-tonian approach to include tunneling between semicon-ducting half-spaces forms the basis for the independentparticle or mean-field type approach[8, 9]. This approachhas been applied to a wide range of tunneling phenomenain physics from superconducting tunneling[3, 6, 7, 10–12]to semiconductor devices[13–18] to scanning tunnelingmicroscopy[19, 20].The observation of photon enhanced tunneling in su-perconducting tunnel junctions under microwave illumi-nation led to theoretical investigations into dynamic tun-neling phenomena[10]. Tien-Gordon proposed a theorybased on time-dependent interaction of the microwavephotons in the tunnel barrier of a superconducting tun-nel junction[10] where they applied the transfer Hamilto-nian to the multi-photon tunneling processes in the insu-lator. This multi-photon enhanced tunneling process insuperconducting tunnel junctions results in plateaus inthe current voltage characteristics in these devices andhas been used as a sensitive direct detector for millime- ∗ [email protected] ter waves[12]. Direct detection based on rectification ofdisplacement currents in tunneling devices have been pro-posed and examined theoretically[21].Photon-assisted tunneling has been experimentally ob-served in antenna-coupled semiconductor superlatticedevices[22], and quantum dots [23, 24]. These meso-scopic systems show evidence of multi-photon assistedtunneling in the current voltage characteristics underhigh frequency microwave illumination at cryogenic tem-peratures. Modified Tien-Gordon models have been de-veloped and well describe the data [11]. Furthermore,Boson-assisted tunneling [25] in single charged defectpoint-contact barriers and quantum dots have shownthe Kondo effect and resonant tunneling in these struc-tures [26, 27]. These results have been modeled us-ing the Anderson model for a single charged or mag-netic impurity, and the tunneling currents computedusing Bardeen’s transfer Hamiltonian. These new as-sisted tunneling observations have led to the develop-ment of advanced modeling methods, including pertur-bative expansions and non-equilibrium Green’s functiontechniques[28] which are closely related techniques to ourderived density matrix expansion at finite temperature.Recently, photon-assisted tunneling in resonant quantumdot systems has been examine as a means of thermoelec-tric energy conversion[29, 30] similar to enhanced tunnel-ing in nanoantenna-coupled MOS tunnel diodes[31, 32].Another key application of tunneling is in examina-tion of leakage currents and gate control in metal-oxide-semiconductor (MOS) device structures owing to its tech-nological importance. The starting point for tunnel-ing calculations is the Bardeen transfer tunnel currentexpression in the static approximation and the barriertransmission expression computed within the WKB ap-proximation. In the WKB approach, the semi-classicalbarrier transmission expression depends on the energy ofthe carriers on the side under bias and is computed forelectron or hole distributions driven out of equilibrium bya bias voltage applied to the device. Parabolic bands inthe effective mass approximation are typically assumed a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n for the metal and semiconductors and the transverse mo-mentum integration is performed to give the usual supplyfunction expression for the tunneling current[13, 17]. Theasymmetirc effective masses in the metal and semicon-ductor give rise to large current asymmetries under bias.Tunneling calculations for non-analytic barriers can becalculated using a transfer matrix approach[33].Recently, photon enhanced tunneling in nanoantenna-coupled MOS tunnel diodes have been studied as a newinfrared direct conversion devices[31, 32] These largearea distributed tunnel diodes concentrate infrared time-varying transverse electromagnetic fields in tunnelingbarriers and have been shown to lead to measurable tun-neling currents and generation of electrical power froma thermal source. In this paper, we examine the time-dependent transfer Hamiltonian in a mean-field densitymatrix perturbative expansion. The transfer Hamilto-nian problem is expanded into complete and orthogo-nal bases on the right and left half-spaces, which are ex-tended over the entire domain but not complete[34, 35].The Bardeen tunneling current formula is found as thefirst order expression in the mean-field current expan-sion. Higher order terms in the density matrix expan-sion are seen to give rise to resonant and non-resonantmulti-photon enhanced tunneling currents. A uniformbarrier tunneling problem is examined under incoherentinfrared illumination between dissimilar materials andsingle-photon enhanced tunneling currents are calculatedfor the asymmetric tunneling problem. II. TRANSFER HAMILTONIAN APPROACH
The transfer Hamiltonian approach to tunneling wasoriginally formulated by Bardeen as a semi-empirical ap-proach to calculate tunneling currents between a super-conductor and a normal metal state through a very smallgap. Since then it has been applied to a host of generaltunneling problems. The transfer Hamiltonian approachis formulated as two half-spaces separated by a thin bar-rier interaction that connects the two half-spaces suchthat we can describe the system phenomenologically as H = H l + H r + H T ( t ) , (1)where r and l refer to right and left half-spaces and H T ( t ) = λ V cos(Ω t ) is the transfer Hamiltonian. Thesplitting of the problem into left right half-spaces im-poses conditions on the completeness of the states. Theeigenstates on the left and right should be complete ontheir respective support and exist in the entire range ofthe total Hamiltonian. This requirement is not compati-ble with the requirement that the left and right Hamilto-nian commute. We therefore give up on the notion thatthe left and right energy is observable and the states ofthe unperturbed Hamiltonian are not left-right productstates [34].With these caveats in mind, the eigenvalue problem inthe two half spaces can be readily stated, for the left half space x <
0, we have | Ψ l (cid:105) = (cid:80) k a k ( t ) | φ k (cid:105) such that H l | φ k (cid:105) = E lk | φ k (cid:105) (2)and we have the half-space orthogonality (cid:104) φ k (cid:48) | φ k (cid:105) l = δ k,k (cid:48) . Similiarly, the right half space x ≥ t ox , we have | Ψ r (cid:105) = (cid:80) q b q ( t ) | ψ q (cid:105) with H r | ψ q (cid:105) = E rq | ψ q (cid:105) , (3)and (cid:104) ψ q (cid:48) | ψ q (cid:105) r = δ q,q (cid:48) . The left and right eigenvectors areorthogonal on their respective half-spaces but extend intothe gap region, 0 < x ≤ t ox and beyond. Many particleGreen’s function approaches have been used to examinedtunneling in the transfer Hamiltonian approach[28, 35].Density matrix approaches are closely related to theseGreen’s function approaches but give simplified time-dependent perturbation expansions and are amenable tomean-field approximations at finite temperature.In the following sections, we will show that the Bardeenexpression for the tunnel current corresponds to a first or-der term in a systematic density matrix expansion in thestrength of the interaction. This expansion in the den-sity matrix will lead to a generalized non-linear currentanalogous to the non-linear polarization from non-linearoptics. These higher order terms correspond to resonantand non-resonant enhanced multi-photon tunneling cur-rents. The transfer interaction will be considered to arisefrom the electromagnetic coupling term confined with inthe tunnel barrier and we will derive the photon-assistedtunneling current induced from a broadband incoherentsource in a simple uniform barrier model. III. MULTI-PHOTON DENSITY MATRIXEXPANSION
The transfer interaction for a confined electromagneticfield in the gap can be examined by considering a semi-classical electromagnetic field interaction that takes onthe familiar form, H T ( t ) = i (cid:126) emc A ( t ) · ∇ , (4)where A ( t ) is the spatially-varying real vector poten-tial in the gap region. The transfer interaction from abroadband source can be expressed in terms of the time-dependent electric field, H T ( t ) = (cid:90) ∞−∞ dω π e − iωt e (cid:126) mω E ( ω ) · ∇ , (5)where we have used the Fourier transform of the vectorpotential to derive the electric field along with the con-straints of real valued fields. This transfer interaction isshown to transfer a quasi-particle from right to left andand vice versa, and for the above form of the interactionis driven be the induced polarization from the electro-magnetic field in the gap. The states of the system arenot product states, but single particle states that are de-fined by projection specified with the general index n .The two non-interacting half-spaces satisfy H | n (cid:105) = E n | n (cid:105) , (6)where H = H l + H r , where E n = { E n r , E n l } dependingon whether projected into right or left half-space. The density matrix operator can be defined in termsof the single particle states and is given by ρ = (cid:88) n,m ρ nm | n (cid:105) (cid:104) m | . (7)The Liouville equation gives the equation of motion forthe density matrix is i (cid:126) ∂ρ nm ∂t = ( E n − E m ) ρ nm + (cid:88) n (cid:48) m (cid:48) (cid:0) (cid:104) n | H T n (cid:48) (cid:105) δ m,m (cid:48) − (cid:104) m (cid:48) H ∗ T | m (cid:105) δ n,n (cid:48) (cid:1) ρ n (cid:48) m (cid:48) . (8)where (cid:104) m | H T n (cid:105) means that H T operates on the nth state.(see Appendix A for general derivation and Fourier anal-ysis of the electromagnetic transfer interaction). The Li-ouville equation for the density matrix can be solved us-ing a formal perturbation expansion in the strength ofthe interaction, H T ( t ) → λH T ( t )), where λ is treatedas a formal parameter in perturbation. Here the densitymatrix expansion is expanded in the interaction strengthto give ρ nm = ρ (0) nm + λρ (1) nm ( t ) + λ ρ (2) nm ( t ) + λ ρ (3) nm ( t ) . . . , (9)where the unperturbed density matrix is ρ (0) nm = f n δ nm and is assumed to be time-independent and diagonal and f n is the static non-equilibrium distribution function.We can perform the same Fourier analysis for the time-dependent density matrix and is derived in Appendix B,to give ρ ( i ) nm ( t ) = ∞ (cid:90) −∞ dω π e − iωt ˜ ρ ( i ) nm ( ω ) . (10)By grouping the terms in the strength of the interaction,we obtain the series expansion terms in the freqency do-main. We have the density matrix expansion terms˜ ρ (0) nm ( ω ) = f n δ n,m δ ( ω ) (11)˜ ρ (1) nm ( ω ) = G nm ( ω ) (cid:88) n (cid:48) m (cid:48) ∞ (cid:90) −∞ dω (cid:48) π ˜ T nm ; n (cid:48) m (cid:48) ( ω − ω (cid:48) )˜ ρ (0) n (cid:48) m (cid:48) ( ω (cid:48) )˜ ρ (2) nm ( ω ) = G nm ( ω ) (cid:88) n (cid:48) m (cid:48) ∞ (cid:90) −∞ dω (cid:48) π ˜ T nm ; n (cid:48) m (cid:48) ( ω − ω (cid:48) )˜ ρ (1) n (cid:48) m (cid:48) ( ω (cid:48) )˜ ρ (3) nm ( ω ) = G nm ( ω ) (cid:88) n (cid:48) m (cid:48) ∞ (cid:90) −∞ dω (cid:48) π ˜ T nm ; n (cid:48) m (cid:48) ( ω − ω (cid:48) )˜ ρ (2) n (cid:48) m (cid:48) ( ω (cid:48) )... ...where G nm ( ω ) = ( (cid:126) ω − E n + E m ) − is the free carrierpropagator, and we define˜ T nm ; n (cid:48) m (cid:48) ( ω ) = (cid:0) V n,n (cid:48) δ m,m (cid:48) − V m (cid:48) ,m δ n,n (cid:48) (cid:1) , (12)and V n,n (cid:48) ( ω ) = e (cid:126) mω (cid:90) d x E ( ω ) · (cid:0) φ ∗ n ( ∇ φ n (cid:48) (cid:1) . (13)The iterated density matrix can be expressed as an in-tegral equation in the frequency domain and is derived in Appendix B. The first order correction to the densitymatrix is ˜ ρ (1) nm ( ω ) = G nm ( ω ) V n,m (cid:0) f m − f n (cid:1) , (14)where we have assumed a uniform transverse electricfield.The hierarchy of density matrices in the dynamic tun-neling case allows for multi-photon processes to enhancethe tunneling current within the device structure. Thedynamic tunneling current and the photon-assisted staticcurrent are derived in Appendix C. The dynamic current -‐0.6 -‐0.4 -‐0.2 0 0.2 0.4 0.6 0.8 1 1.2 -‐50 -‐30 -‐10 10 30 50 70 90 -‐1.5 -‐1 -‐0.5 0 0.5 1 1.5 -‐50 -‐30 -‐10 10 30 50 70 90 FIG. 1.
Illustration transfer interaction for a simple 1D potential barrier. (a) The schematic shows the right andleft wave functions, φ k and ψ q , respectively. These become evanescent in the potential barrier region. The overlap in thisregion leads to the tunneling current. (b) Energy level diagram showing the one photon tunneling process at zero bias. Thequasiparticle occupancy of the left and right half-spaces are shown schematically. The schematic impact of applying a bias(dashed black lines) leads to shift in the bands on the right side of the barrier with a spatially varying potential in the barrier. is given by I ( ω (cid:48) ) = i e (cid:126) (cid:88) nm ∞ (cid:90) −∞ dω π V m,n ( ω ) ˜ ρ nm ( ω (cid:48) − ω ) , (15)and the static tunneling current is I = − i e (cid:126) (cid:88) nm ∞ (cid:90) −∞ dω π ˜ V m,n ( ω ) ˜ ρ nm ( ω ) . (16)Here we have˜ V m,n (cid:48) ( ω ) = e (cid:126) mω (cid:90) d x E ∗ ( ω ) · (cid:0) φ ∗ m ( ∇ φ n (cid:48) (cid:1) . (17)The expansion of the density matrices in the strengthof the interaction leads to a similar expansion of the staticand dynamic currents. We have that I = I (0) + λI (1) + λ I (2) + λ I (3) . . . (18)from the expansion of ˜ ρ . This expression gives rise to anexpansion for the non-linear conductances in the same fashion as the non-linear susceptibilities arise in non-linear optics. Note, the static current at zero order isnot included in the expansion, the static part of the den-sity matrix is included in the dynamic zeroth order termand has the frequency dependence of the electric fieldonly.The first order term in the static current expressioncan be written explicitly, I (1) = − i e (cid:126) (cid:88) nm ∞ (cid:90) −∞ dω π ˜ V m,n V n,m G nm ( ω ) (cid:0) f m − f n (cid:1) , (19)where we have used Eq. 14 in obtaining the above cur-rent. The first order current can be simplified if we con-sider the special case of a uniform transverse electric fieldin the gap and use the principal value expression for theGreen’s function energy denominator G nm ( ω ) = P (cid:18) (cid:126) ω − E n + E m (cid:19) + iπδ ( (cid:126) ω − E n + E m ) , (20)to obtain, I (1) = 2 πe (cid:126) ∞ (cid:90) −∞ dω π (cid:18) e (cid:126) mω (cid:19) | E z ( ω ) | (cid:88) nm u mn u nm (cid:0) f m − f n (cid:1) δ ( (cid:126) ω − E n + E m ) , (21)where u nm = (cid:104) n |∇ z m (cid:105) ,and which represents the first or-der real-valued static current induced by a frequency-dependent transverse electric field. Eq. 21 is a general-ized Bardeen transfer expression for photon-induced tun- neling from a broadband radiation source. The frequencyintegral can be restricted to positive frequencies leadingto two contributions to the photon-assisted current cor-responding to absorption and emission of a photon.The higher order currents can be derived by iteratingthe density matrix expansion in Eq. 12, but lead to much more complex expressions. The second order static cur-rent can be obtained in the uniform field case, I (2) = − i e (cid:126) ∞ (cid:90) −∞ dω π (cid:18) e (cid:126) mω (cid:19) ∞ (cid:90) −∞ dω (cid:48) π (cid:18) e (cid:126) mω (cid:48) (cid:19) (cid:18) e (cid:126) m ( ω − ω (cid:48) ) (cid:19) E ∗ z ( ω ) E z ( ω (cid:48) ) E z ( ω − ω (cid:48) ) (cid:88) ln (cid:48) m (cid:48) u m (cid:48) l u ln (cid:48) u n (cid:48) m (cid:48) [ f m (cid:48) − f n (cid:48) ] G n (cid:48) m (cid:48) ( ω (cid:48) ) ( G m (cid:48) l ( ω ) − G ln (cid:48) ( ω )) , (22)where u nm = (cid:104) n |∇ z m (cid:105) . The 2nd order current is ob-tained from the real value of the above expression andis cubic in the electric field strength and has a complexfrequency dependence. In the following, we will considera special example of a uniform barrier with different con-tact materials on either side of the barrier. This sim-ple problem will be evaluated using our 1st order staticcurrent result (Eq. 21) and will show photon-assistedtunneling in an simple example problem with technicalrelevance. IV. UNIFORM BARRIER EXAMPLE
It is instructive to examine a simple example. We con-sider tunneling through a uniform barrier in the pres-ence of a transverse uniform time varying electromagneticfield. The electric field is from an incoherent black-bodysource with a small and finite bandwidth. We consider asimple transfer Hamiltonian of the form, H = (cid:18) p m ox + V (cid:19) + H T ( t ) , (23)where V is a real-valued constant barrier height in theregion 0 < x ≤ t ox , and last term is the imaginary poten-tial that arises from time-varying vector potential. Eq.21 can be applied to the uniform barrier problem to eval-uate the first order tunneling current. The wavefunctionsin the barrier from the left and right regions are given by φ k = √ p l exp( i k · r − kx ) (24) ψ q = √ p l exp( i q · r − q ( t ox − x )) (25)where k and q are the transverse momentum, and k = (cid:112) m l ( V − E l ) / (cid:126) and q = (cid:112) m r ( V − E r ) / (cid:126) , where V is the barrier height. Here p l,r = 1 /V l,r are the real-valued volume normalizations in the left and right half-spaces. Figure 1 (a) shows schematically the wavefunc-tions in the two half-spaces and the overlap in the barrier. The wavefunction in the barrier region exponentially de-cays away from the left and right interfaces respectively.The tunneling occurs due to the dipole polarization in-duced by the field confined in the gap and the wavefunc-tion overlap in this region.The matrix element of the transfer interaction arereadily evaluated u qk = − k T qk (2 π ) δ ( q − k ) , (26)and similiarly for the transpose u kq = q T kq (2 π ) δ ( k − q ) , (27)where transverse momentum conservation is explicit, k = q , and T kq = T qk is symmetric with T qk = t ox √ V l V r e − ( k + q ) t ox / sinh(( k − q ) t ox / k − q ) t ox / . (28)The imaginary part of the transfer interaction matrix el-ements occurs due to the momentum operator matrix el-ements with the evanescent decaying wavefunctions. Theimaginary transfer interaction is needed to give rise to anon-zero current through the barrier and implies that weare not in an equilibrium situation but a dynamic situ-ation in which quasiparticles are injected and removedfrom a reservior at the left and right contacts.The first order photon transfer tunneling current isgiven by Eq. 21, and the sums are over both transverseand longitudinal wavevectors. These sums are convertedto integrals (cid:88) k → V l (cid:90) dk π (cid:90) d k (2 π ) (29)over longitudinal and transverse momenta and similarlyfor right half-space momenta q . The current density canbe obtained by integrating over the transverse momen-tum in the left and right half-spaces. This can be doneanalytically if the parallel effective masses in the left andright half-spaces are assumed to be equal, m e , and theevanescent x wavevectors do not depend on the trans-verse momenta. The single photon tunneling current is V=1.0 eV V=1.25 eV V=1.5 eV V=1.75 eV -10 -9 -8 -7 -6 -5 -4 -3 -2 tox=1nm tox=1.3nm tox=1.5nm tox=2.0nm tox=2.5nm -10 -9 -8 -7 -6 -5 -4 -3 -2 Temperature ( o C) J ( A / c m ) J ( A / c m ) Temperature ( o C) FIG. 2.
Tunneling current density for various barrier parameters. (a) Computed tunneling current for varying barrierenergy V for fixed barrier thickness t ox = 1 . t ox for fixed barrier energy, V = 1 . m l = m e , m r = 0 . m e , and m ox = 0 . m e where m e is the bare electron mass. The photon energy (cid:126) ω = 0 .
17 eV corresponding toa wavelength of approximately 7.3 microns. The field in the barrier, E , is obtained from a blackbody distribution for thespecified temperature at the above central wavelength with dν = 1 THz bandwidth. J (1) = 4 πem e h β (cid:18) m r m l m ox (cid:19) ∞ (cid:90) dω π (cid:18) et ox (cid:126) ω (cid:19) | E ( ω ) | T ( ω ) , (30)where the transmittance is defined as T ( ω ) = β (cid:90) dE l t ( E l , E l + (cid:126) ω ) log (cid:20) − βE l )1 + exp( − β ( E l + (cid:126) ω ) (cid:21) + ( ω → − ω ) , (31)and where β = 1 /k B T , and the photon energy (cid:126) ω actsas an effective bias in the supply function. There is asmaller back-flow current contribution that correspondsto photon emission in the barrier, E l − (cid:126) ω . This con-tribution to the transmittance has been ignored. In Eq.30, the first term has the units for the current density, AT where A = 120 Amps / cm K is Richardson’s con-stant. The second expression in parenthesis is the ratioof the x-directed perpendicular effective masses for theright and left half-spaces and oxide barrier. The thirdterm is dimensionless interaction strength and dependson the amplitude of the field in the barrier, | E | and thetransmittance, T ( ω ) . This integral is over all positivefrequencies, but can be limited to a finite bandwidth dν near the photon energy. The transmittance is the scaledintegral over the barrier transmittance and supply func-tion give the overall transmittance. The barrier trans-mittance is t ( E l , E r ) = e − ( k + q ) t ox sinhc (cid:18) ( k − q ) t ox (cid:19) , (32)where the first term has the form of the standard tun-neling exponential and the second term is the squarehyperbolic sinc function. The evanescent wavevectors in the barrier region are k = (2 m l ( V − E l ) / (cid:126) ) / and q = (2 m r ( V − E r ) / (cid:126) ) / as previously defined.Figure 2(a) shows the computed first order current forfixed barrier thickness and (b) for fixed barrier heightas the device temperature is varied from 25 o C to 425 o C.A key component of the calculation is the photon fieldamplitude in the barrier, | E | . The photon field amplitudeis derived from a black-body spectral exitance, M ν , givenby d Φ dν = M ν ( T ) = 2 πhν c βhν ) − ν = ω/ π is thephoton frequency, and dν is the bandwidth. The electricfield amplitude in the barrier is | E ( ω ) | = 2 Z M ν ( T ) (34)where Z is the impedance of free-space. In fig. 2 (a) wesee that there is strong variation of the current densitywith barrier height with increasing current with decreas-ing barrier height as expected. There is also a strong tem-perature dependence with the current density sharply in-creasing above a critical temperature. This critical tem-perature decreases with lowered barrier and is due to thethermal occupancy of the left and right energy bandsand the increase of the thermal photon electric field am-plitude in the barrier with increasing temperature. Fig-ure 2(b) shows the dependence of the the tunneling cur-rent density on the barrier thickness. The tunneling cur-rent increases with decreasing barrier thickness as ex-pected. Interestingly, the current density above a criticaltemperature asymptotically approaches an common ex-ponential growth for increasing temperatures. The bar-rier transmittance and the dimensionless field amplitudeterms in Eq.30 depend exponentially and quadraticallyon the barrier thickness, respectively, and are the sourceof this asymptotic behavior. Short circuit tunneling cur-rent in large area nanoantenna coupled tunnel diodes il-luminated a thermal source have been reported experi-mentally [31, 32]. In these devices, high field concentra-tion in the tunnel barrier from a thermal infrared sourceslead to high short currents with strong temperature de-pendence and may lead to new radiative thermoelectricdirect conversion of thermal radiation into direct electricpower. V. CONCLUSIONS
We have derived the transfer Hamiltonian tunnelingcurrent using a time dependent mean-field density matrixperturbation expansion. The density matrix is computedto all orders in perturbation theory and the generalizedtunneling current is obtained. The Bardeen tunnelingcurrent expression is obtained to first order in the densitymatrix expansion. Photon-assisted tunneling arises dueto confined electromagnetic fields in the tunnel barrier re-gion. These fields and wavefunction overlap in the barrierproduce a time-varying polarization which leads to a netcurrent flow. This complex interaction in the barrier is adirect consequence of the imaginary momenta in the bar-rier and results in an imaginary transfer interaction thatleads to the tunneling current. The dependence on thephoton energy in the density matrix hierarchy gives rise to photon-assisted tunneling through the barrier, wherethe photon energy acts as an effective bias.To examine the effective of photon-assisted tunneling,we examine a uniform barrier model with asymmetric ef-fective masses in the two half-spaces. We consider theconfined field in the tunnel gap arising from black-bodyirradiation consistent with the experimental observationof short-circuit tunneling current observed in infrarednanoantenna-coupled tunnel diode rectfiiers[31, 32]. Inthese devices, electric field confinement in the tunnel gapresults from strong photon-phonon interaction near a lon-gitudinal optical phonon mode in the oxide tunnel bar-rier. We find that the tunneling current has a tempera-ture dependent turn on as the source temperature is in-creased, and is seen to saturate at different temperaturesfor varying barrier thicknesses. Multi-photon assistedtunneling arises from higher order terms of eE t ox / (cid:126) ω ,which is small and represents a minor contribution to thefirst order Bardeen tunneling current computed for theuniform barrier model considered. In this paper, we havepresented a generalized density matrix approach to thecomputation of multi-photon assisted tunneling currentsand examined a simple 1D asymmetric tunneling model.The effect of applying a bias and the role of surface statesare challenging problems under consideration in photon-assisted tunneling devices and need to be solved withina consistent many-body framework. ACKNOWLEDGMENTS
Funding for this work was provided by Sandia’s Labo-ratory Directed Research and Development (LDRD) pro-gram. Sandia National Laboratories is a multi-missionlaboratory managed and operated by National Technol-ogy and Engineering Solutions of Sandia, a wholly ownedsubsidiary of Honeywell International Inc., for the UnitedStates Department of Energy’s National Nuclear SecurityAdministration under contract DE-NA0003525. [1] E. L. Murphy and R. H. Good, “Thermionic emission,field emission, and the transition region,” Phys. Rev. , 1464–1473 (1956).[2] J. L. Dunham, “The wentzel-brillouin-kramers methodof solving the wave equation,” Phys. Rev. , 713–720(1932).[3] J. Bardeen, “Tunnelling from a many-particle point ofview,” Phys. Rev. Lett. , 57–59 (1961).[4] Ivar Giaever, “Electron tunneling between two supercon-ductors,” Physical Review Letters , 464 (1960).[5] Ivar Giaever, “Energy gap in superconductors measuredby electron tunneling,” Phys. Rev. Lett. , 147–148(1960).[6] John Bardeen, “Tunneling into superconductors,” Phys.Rev. Lett. , 147–149 (1962). [7] M. H. Cohen, L. M. Falicov, and J. C. Phillips, “Su-perconductive tunneling,” Phys. Rev. Lett. , 316–318(1962).[8] Walter A. Harrison, Solid state theory (Dover, 1979).[9] Walter A. Harrison, “Tunneling from an independent-particle point of view,” Phys. Rev. , 85–89 (1961).[10] A. H. Dayem and R. J. Martin, “Quantum interaction ofmicrowave radiation with tunneling between supercon-ductors,” Phys. Rev. Lett. , 246–248 (1962).[11] P. K. Tien and J. P. Gordon, “Multiphoton process ob-served in the interaction of microwave fields with the tun-neling between superconductor films,” Phys. Rev. ,647–651 (1963).[12] John R. Tucker and Marc J. Feldman, “Quantum de-tection at millimeter wavelengths,” Rev. Mod. Phys. , , 2272–2285 (2001),http://dx.doi.org/10.1063/1.1337596.[14] Leo Esaki, “Long journey into tunneling,” Rev. Mod.Phys. , 237–244 (1974).[15] L.B. Freeman and W.E. Dahlke, “Theory of tunnelinginto interface states,” Solid-State Electronics , 1483 –1503 (1970).[16] Evan O. Kane, “Theory of tunneling,” Jour-nal of Applied Physics , 83–91 (1961),http://dx.doi.org/10.1063/1.1735965.[17] Simon M Sze and Kwok K Ng, Physics of semiconductordevices (John wiley & sons, 2006).[18] Claude Weisbuch and Borge Vinter,
Quantum semicon-ductor structures: fundamentals and applications (Aca-demic press, 2014).[19] C Caroli, R Combescot, P Nozieres, and D Saint-James,“Direct calculation of the tunneling current,” Journal ofPhysics C: Solid State Physics , 916 (1971).[20] C. Julian Chen, “Tunneling matrix elements in three-dimensional space: The derivative rule and the sumrule,” Phys. Rev. B , 8841–8857 (1990).[21] A. Sanchez, C. F. Davis Jr., K. C. Liu, and A. Ja-van, “The mom tunneling diode: Theoretical estimateof its performance at microwave and infrared frequen-cies,” Journal of Applied Physics , 5270–5277 (1978),http://dx.doi.org/10.1063/1.324426.[22] B. J. Keay, S. J. Allen, J. Gal´an, J. P. Kaminski,K. L. Campman, A. C. Gossard, U. Bhattacharya, andM. J. W. Rodwell, “Photon-assisted electric field domainsand multiphoton-assisted tunneling in semiconductor su-perlattices,” Phys. Rev. Lett. , 4098–4101 (1995).[23] L. P. Kouwenhoven, S. Jauhar, J. Orenstein, P. L.McEuen, Y. Nagamune, J. Motohisa, and H. Sakaki,“Observation of photon-assisted tunneling through aquantum dot,” Phys. Rev. Lett. , 3443–3446 (1994).[24] L. P. Kouwenhoven, S. Jauhar, K. McCormick, D. Dixon,P. L. McEuen, Yu. V. Nazarov, N. C. van der Vaart,and C. T. Foxon, “Photon-assisted tunneling through aquantum dot,” Phys. Rev. B , 2019–2022 (1994).[25] D. C. Ralph and R. A. Buhrman, “Kondo-assisted and resonant tunneling via a single charge trap: A realizationof the anderson model out of equilibrium,” Phys. Rev.Lett. , 3401–3404 (1994).[26] J¨urgen K¨onig, Herbert Schoeller, and Gerd Sch¨on,“Zero-bias anomalies and boson-assisted tunnelingthrough quantum dots,” Phys. Rev. Lett. , 1715–1718(1996).[27] J¨urgen K¨onig, J¨org Schmid, Herbert Schoeller, and GerdSch¨on, “Resonant tunneling through ultrasmall quantumdots: Zero-bias anomalies, magnetic-field dependence,and boson-assisted transport,” Phys. Rev. B , 16820–16837 (1996).[28] Gloria Platero and Ram´on Aguado, “Photon-assistedtransport in semiconductor nanostructures,” Physics Re-ports , 1 – 157 (2004).[29] Andrew N. Jordan, Bj¨orn Sothmann, Rafael S´anchez,and Markus B¨uttiker, “Powerful and efficient energyharvester with resonant-tunneling quantum dots,” Phys.Rev. B , 075312 (2013).[30] Bj¨orn Sothmann, Rafael S´anchez, and Andrew N Jordan,“Thermoelectric energy harvesting with quantum dots,”Nanotechnology , 032001 (2015).[31] Paul S Davids, Robert L Jarecki, Andrew Starbuck,D Bruce Burckel, Emil A Kadlec, Troy Ribaudo, Eric AShaner, and David W Peters, “Infrared rectification in ananoantenna-coupled metal-oxide-semiconductor tunneldiode,” Nature nanotechnology , 1033 (2015).[32] Emil A. Kadlec, Robert L. Jarecki, Andrew Starbuck,David W. Peters, and Paul S. Davids, “Photon-phonon-enhanced infrared rectification in a two-dimensionalnanoantenna-coupled tunnel diode,” Phys. Rev. Applied , 064019 (2016).[33] Yuji Ando and Tomohiro Itoh, “Calculation of transmis-sion tunneling current across arbitrary potential barri-ers,” Journal of Applied Physics , 1497–1502 (1987),http://dx.doi.org/10.1063/1.338082.[34] Richard E. Prange, “Tunneling from a many-particlepoint of view,” Phys. Rev. , 1083–1086 (1963).[35] Joel A. Appelbaum and W. F. Brinkman, “Theory ofmany-body effects in tunneling,” Phys. Rev. , 464–470 (1969). Appendix A: Generalized density matrix
The transfer Hamiltonian under consideration explicitly in this paper is given by the electromagnetic interaction inthe tunnel barrier, we have H T ( t ) = i e (cid:126) mc A ( t ) · ∇ . (A1)The real-valued vector potential can be written in terms of a Fourier transform A ( t ) = (cid:90) ∞−∞ dω π e − iωt A ( ω ) , (A2)and we require A ∗ ( − ω ) = A ( ω ) for real-valued fields. The vector potential in the Coulomb gauge can be expressed interms of the electric field, we have E ( ω ) = iω/c A ( ω ) and we have E ∗ ( − ω ) = E ( ω ). The resulting transfer interactionis given by H T ( t ) = (cid:90) ∞−∞ dω π e − iωt e (cid:126) mω E ( ω ) · ∇ , (A3)and we find that H ∗ T ( t ) = − H T ( t ) but care must be taken regarding order of operations. In the following, we willexamine the equations of motion for the density matrix in great detail for clarity.The density matrix can derived for the transfer interaction. The Hamiltonian for the system is H = H + H T ( t ),with H | n (cid:105) = E n | n (cid:105) where H is the unperturbed Hamiltonian with orthonormal states, (cid:104) m | n (cid:105) = δ nm . Schrodinger’sequation of motion is i (cid:126) ∂ t | Ψ (cid:105) = ( H + H T ( t )) | Ψ (cid:105) , (A4)and the total wavefunction can be expanded in the unperturbed basis | Ψ (cid:105) = (cid:80) n c n ( t ) | n (cid:105) , to give i (cid:126) ˙ c n = E n c n + (cid:88) n (cid:48) (cid:104) n | H T n (cid:48) (cid:105) c n (cid:48) , (A5)and a similar expression for the conjugate gives i (cid:126) ˙ c † m = − E m c † m − (cid:88) m (cid:48) (cid:104) m (cid:48) H ∗ T | m (cid:105) c † m (cid:48) . (A6)The generalized density operator is ˆ ρ = | Ψ (cid:105) (cid:104) Ψ | = (cid:88) n,m c n ( t ) c † m ( t ) | n (cid:105) (cid:104) m | , (A7)where ρ nm = c n ( t ) c † m ( t ). The density matrix equation of motion is i (cid:126) ∂ρ nm ∂t = ( E n − E m ) ρ nm + (cid:88) n (cid:48) m (cid:48) (cid:0) (cid:104) n | H T n (cid:48) (cid:105) δ m,m (cid:48) − (cid:104) m (cid:48) H ∗ T | m (cid:105) δ n,n (cid:48) (cid:1) ρ n (cid:48) m (cid:48) . (A8)It is convenient to define the transfer matrix T nm ; n (cid:48) m (cid:48) = (cid:0) (cid:104) n | H T n (cid:48) (cid:105) δ m,m (cid:48) − (cid:104) m (cid:48) H ∗ T | m (cid:105) δ n,n (cid:48) (cid:1) . (A9)The matrix elements can be evaluated in the transfer matrix, we have T nm ; n (cid:48) m (cid:48) = (cid:90) d x (cid:90) ∞−∞ dω π e − iωt e (cid:126) mω E ( ω ) · (cid:0) φ ∗ n ( ∇ φ n (cid:48) ) δ m,m (cid:48) − φ ∗ m (cid:48) ( ∇ φ m ) δ n,n (cid:48) + ∇ ( φ ∗ m (cid:48) φ m ) (cid:1) , (A10)where the last term can be converted to a divergence which vanishes on the bounding surface and a term proportionalto ∇ · E = 0. Therefore, we obtain T nm ; n (cid:48) m (cid:48) = (cid:90) ∞−∞ dω π e − iωt e (cid:126) mω (cid:90) d x E ( ω ) · (cid:0) φ ∗ n ( ∇ φ n (cid:48) ) δ m,m (cid:48) − φ ∗ m (cid:48) ( ∇ φ m ) δ n,n (cid:48) (cid:1) , (A11)where the equation of motion becomes i (cid:126) ∂ρ nm ∂t = ( E n − E m ) ρ nm + (cid:88) n (cid:48) m (cid:48) T nm ; n (cid:48) m (cid:48) ( t ) ρ n (cid:48) m (cid:48) . (A12) Appendix B: Time-dependent perturbation expansion
The expansion of the density matrix in the strength of the interaction is given by Eq. 9 in the text, we obtain0 i (cid:126) ˙ ρ (1) nm − ( E n − E m ) ρ (1) nm = (cid:88) n (cid:48) m (cid:48) ∞ (cid:90) −∞ dω π e − iωt T nm ; n (cid:48) m (cid:48) ρ (0) n (cid:48) m (cid:48) (B1) i (cid:126) ˙ ρ (2) nm − ( E n − E m ) ρ (2) nm = (cid:88) n (cid:48) m (cid:48) ∞ (cid:90) −∞ dω π e − iωt T nm ; n (cid:48) m (cid:48) ρ (1) n (cid:48) m (cid:48) (B2) i (cid:126) ˙ ρ (3) nm − ( E n − E m ) ρ (3) nm = (cid:88) n (cid:48) m (cid:48) ∞ (cid:90) −∞ dω π e − iωt T nm ; n (cid:48) m (cid:48) ρ (2) n (cid:48) m (cid:48) (B3)... ...and we recall that the zeroth order term is time independent and diagonal. We can perform the same Fourier analysisfor the time-dependent density matrix, ρ ( i ) nm ( t ) = ∞ (cid:90) −∞ dω π e − iωt ˜ ρ ( i ) nm ( ω ) . (B4)The first order correction to the density matrix is˜ ρ (1) nm ( ω ) = 1 (cid:126) ω − E n + E m (cid:18) e (cid:126) mω (cid:19) E ( ω ) · (cid:104) n |∇ m (cid:105) (cid:0) f m − f n (cid:1) , (B5)where we have assumed a uniform transverse electric field. The second order density matrix and higher order termsare ˜ ρ (2) nm ( ω ) = 1 (cid:126) ω − E n + E m (cid:88) n (cid:48) m (cid:48) ∞ (cid:90) −∞ dω (cid:48) π ˜ T nm ; n (cid:48) m (cid:48) ( ω − ω (cid:48) )˜ ρ (1) n (cid:48) m (cid:48) ( ω (cid:48) ) (B6)˜ ρ (3) nm ( ω ) = 1 (cid:126) ω − E n + E m (cid:88) n (cid:48) m (cid:48) ∞ (cid:90) −∞ dω (cid:48) π ˜ T nm ; n (cid:48) m (cid:48) ( ω − ω (cid:48) )˜ ρ (2) n (cid:48) m (cid:48) ( ω (cid:48) ) (B7)... ...where we note that ˜ T nm ; n (cid:48) m (cid:48) ( ω ) = e (cid:126) mω (cid:90) d x E ( ω ) · (cid:0) φ ∗ n ( ∇ φ n (cid:48) ) δ m,m (cid:48) − φ ∗ m (cid:48) ( ∇ φ m ) δ n,n (cid:48) (cid:1) . (B8)The density expansion can be considered a consequence of an iterated integral equation for the total density matrix.We have ˜ ρ nm ( ω ) = ˜ ρ (0) nm ( ω ) + 1 (cid:126) ω − E n + E m (cid:88) n (cid:48) m (cid:48) ∞ (cid:90) −∞ dω (cid:48) π ˜ T nm ; n (cid:48) m (cid:48) ( ω − ω (cid:48) )˜ ρ n (cid:48) m (cid:48) ( ω (cid:48) ) , (B9)which has the form of a scattering equation. The matrix element ˜ T acts as a self-energy term in the non-equilibriumdensity matrix and the above equation can be formally inverted for the non-equilibrium density matrix. Appendix C: Current continuity
The transfer Hamiltonian for the system is given by H = H + H T ( t ) , (C1)where in general H = − (cid:126) m ∇ + U, (C2)1with U the spatially dependent real-valued potential. The expansion of the wavefunction is also as given above, wehave Ψ = (cid:88) n c n ( t ) φ n ( x ) . (C3)From the equation of motion and the definition of the density matrix given in pervious Appendices, we obtain (cid:88) n,m (cid:18) ∂ρ nm ∂t φ ∗ m φ n + ∇ · j nm ρ nm (cid:19) = 2 e (cid:126) (cid:88) nm ∞ (cid:90) −∞ dω π e − iωt ω E ( ω ) · j nm ρ nm , (C4)where the current is defined as j nm = (cid:126) mi (cid:0) φ ∗ m ( ∇ φ n ) − ( ∇ φ ∗ m ) φ n (cid:1) . (C5)If we now integrate over an arbitrary volume, we can apply Green’s theorem to convert to a surface integral over thebounding volume to obtain ∂N∂t + (cid:90) (cid:90) dA ˆ n · j = 2 e (cid:126) ∞ (cid:90) −∞ dω π e − iωt ω (cid:90) V d x E ( ω ) · j , (C6)where the current is j = (cid:80) nm j nm ρ nm , and number of quasi-particles is N = (cid:88) nm ρ nm (cid:90) V d x φ ∗ m φ n . (C7)The charge current is obtained by multiplying by the − e electronic charge and the electrical current density is J = − e j .The resulting expression for the electrical current is I = 2 e (cid:126) ∞ (cid:90) −∞ dω π e − iωt ω (cid:90) V d x E ( ω ) · J + e ∂N∂t , (C8)where I = − e (cid:90) (cid:90) dA ˆ n · J . (C9)Eq. C8 is the general equation for current continuity in the presence of dissipation. If we ignore the time-varyingcharge density term, we obtain I = 2 e (cid:126) ∞ (cid:90) −∞ dω π e − iωt ω (cid:90) V d x E ( ω ) · J . (C10)The current from Eq. C10 can be expressed in terms of the transfer Hamiltonian. Explicitly, we have I = i (cid:88) nm ∞ (cid:90) −∞ dω π e − iωt e mω (cid:90) V d x E ( ω ) · (cid:0) φ ∗ m ( ∇ φ n ) − ( ∇ φ ∗ m ) φ n (cid:1) ρ nm . (C11)Using the previous integration by parts, we get I = i e (cid:126) (cid:88) nm ∞ (cid:90) −∞ dω π e − iωt e (cid:126) mω (cid:90) V d x E ( ω ) · (cid:0) φ ∗ m ( ∇ φ n ) (cid:1) ρ nm , (C12)or I ( ω (cid:48) ) = i e (cid:126) (cid:88) nm ∞ (cid:90) −∞ dω π e (cid:126) mω (cid:90) V d x E ( ω ) · (cid:0) φ ∗ m ( ∇ φ n ) (cid:1) ˜ ρ nm ( ω (cid:48) − ω ) , (C13)2and since we are mainly interested in the DC current that results from the dynamic electromagnetic field, we evaluateat ω (cid:48) = 0 I = − i e (cid:126) (cid:88) nm ∞ (cid:90) −∞ dω π e (cid:126) mω (cid:90) V d x E ∗ ( ω ) · (cid:0) φ ∗ m ( ∇ φ n ) (cid:1) ˜ ρ nm ( ω ) , ., .