Correlated electron systems periodically driven out of equilibrium: Floquet + DMFT formalism Naoto Tsuji, Takashi Oka, and Hideo Aoki Department of Physics, University of Tokyo, Hongo, Tokyo 113-0033, Japan (Dated: January 23, 2009) Wepropose tocombinetheFloquet formalism for systems in ac fieldswith thedynamical mean- field theory to study correlated electron systems periodically driven out of equilibrium by external fieldssuch asintense laser light. This approach hasa virtue that wecan nonperturbatively include both the correlation effects and nonlinear effects due to the driving field, which is imperative in 9 analyzingrecentexperimentsforphotoinducedphasetransitions. Insolvingtheproblem,weexploit 0 a general theorem that the Hamiltonian in a Floquet matrix form can be exactly diagonalized 0 for single-band noninteracting systems. As a demonstration, we have applied the method to the 2 Falicov-Kimball model in intense ac fieldsto calculate thespectral function. The result shows that photoinducedmidgap states emerge from strong ac fields, triggering an insulator-metal transition. n a J PACSnumbers: 71.27.+a,78.20.Bh,71.15.−m 3 2 I. INTRODUCTION method (FGFM)8,9,10,11 is plugged into the dynamical mean-field theory (DMFT).12 This formulation provides ] l us with a promising way to treat both of the nonlinear e Controlling phases of matter is a central issue in the - physics of strongly correlated electron systems, where a effectoftheelectricfieldandthecorrelationeffectsimul- r taneously in a nonperturbative manner. We then apply t rich variety of phases are realized for various physical s the method to study the responseofthe Falicov-Kimball . degrees of freedom such as spin and charge. Since a t (FK) model,13,14 one of the simplest lattice models for a phase transitiondramatically alters macroscopicproper- m ties of the system, it is ofgreatimportance to know how correlated electrons, to ac electric fields. A particular emphasis is put on the technical basis of FGFM. the phase changes as external parameters (temperature, - d pressure, band filling of the system, etc.) are varied, in Duringthe preparationofthe presentstudy, wenotice n equilibrium. thatJouraet al.15 adoptedatechniquesimilartoFGFM o Now, the last decade witnessed that a controlling to solve the Dyson equation for a system in dc fields. c [ factor, intense laser fields, can trigger “phase transi- Here we present a general framework of DMFT out of tions” in correlated electron materials.1,2 One repre- equilibrium for arbitrary time-periodic fields in a more 2 sentative class of materials is perovskite manganites, transparent viewpoint exploiting the Floquet formalism. v in which an insulator-to-metal transition is induced by In this context we should emphasize that FGFM is not 9 photoexcitation.3,4Recentexperimentsalsoindicatethat just a numerical technique but offers a fruitful physical 7 3 a ferromagnetic spin alignment can emerge in the in- picture for nonequilibrium systems as revealed here. 0 duced metallic phase in manganites.5 Such phenom- Thepaperisorganizedasfollows: inSec.II,wereview . ena are called photoinduced phase transitions (PIPTs), theFloquettheoremandthe Floquetmatrixmethod,on 8 0 where irradiation of photons allows one to change elec- whichourtheoreticaldescriptionisbased. Therestofthe 8 tronic, magnetic, optical, or structural properties of the paper is devoted to our original formulation. In Sec. III 0 system. we define the Floquet matrix form of Green’s function, : The photoinduced phase transition, however, is dis- which is our starting point of FGFM. Then in Sec. IV, v i tinct from conventional phase transitions in equilibrium we derive a general expression for Green’s function and X in that the photon field drives the system out of equilib- its inverse for noninteracting electrons in a Floquet ma- r rium. In other words, the change in phases in nonequi- trixform. We calculatenoninteractingGreen’sfunctions a libriumchallengesourunderstandingofphasetransitions for several cases to discuss their physical implications. whichisnormallyconceivedinequilibrium. InthePIPT, InSec.V,wederiveageneraltheoremthatidentifiesthe we haveto consider,ontopofthe nonlinearelectric-field eigenvaluesandtheeigenvectorsofaFloquetmatrixform effect,theelectroncorrelationeffect. Thenonlineareffect of the Hamiltonian for noninteracting electrons. While appearsasathresholdbehavior,i.e.,amacroscopictran- Secs.IVandVaddressnoninteractingcases,wemoveon sition only occurs when the intensity of the driving field to correlated electron systems by incorporating FGFM exceedsacertainstrength. Ontheotherhand,important inDMFT inSec.VI. WeformulatetheobtainedGreen’s correlation effects such as Mott’s metal-insulator transi- function in a gauge-invariant manner in Sec. VII. We tionarenonperturbativeeffects. Thus,wecannotemploy then apply our method to the FK model in Sec. VIII, thelinear-responsetheory6 noramean-fieldtreatmentof and calculate the spectral functions for dc and ac fields, the electron-electron interaction. where we discuss how the Mott-insulating state is trans- Here we propose a theoretical approach7 for photoin- formedinto the metallic states inthe externalfields. We duced phenomena, where the Floquet-Green function summarizethepaperandgivefutureproblemsinSec.IX. 2 II. FLOQUET THEOREM AND FLOQUET dependent, the energy is not conserved in general. How- MATRIX ever,Eq.(4)showsthattheenergyisconserveduptoan integermultiple ofΩ,correspondingto theabsorptionor Anexternalfielddrivesanelectronhavinganenergyε emissionof the photon with the energy Ω. Eachelement intoanotherstatewithadifferentenergy,wherethereare in the Floquet matrix Hmn corresponds to the proba- many scatteringchannels innonequilibrium. However,if bility amplitude of the transition from the mth Floquet the driven system is periodic in time with a frequency mode to the nth one, so that off-diagonal components Ω,the allowedchannelsarelimited tothe processessuch represent excitations driven by the external field while that ε ε+nΩ, where n is an integer. This greatly the diagonal ones the probability to remain in the same reduces→thechannels’degreesoffreedomtobedealtwith. mode. One can take advantage of such a simplification through The consequence of the Floquet theorem is remark- the Floquet matrix method,16,17,18 of which we give an able: Eq.(4)resemblesthestaticSchr¨odingerequationin overview in this section. equilibrium except for the presence of the Floquet mode The method has been used as an effective approach index n, which means that we have no longer to solve toward photoexcited systems. The concept of the Flo- the time-dependent Schr¨odinger Eq. (1), in favor of the quet matrix originates from the Floquet theorem19 for a time-independent Eq.(4). This is the greatadvantage of periodically driven system, an analog of the Bloch theo- the Floquet matrix method, which also plays a crucial rem applied to a spatially periodic system. Floquet the- role in Green’s-function approach. orem is a general theorem for differential equations of a form dx(t)/dt = C(t)x(t) with C periodic in t, which III. FLOQUET REPRESENTATION OF include equationsofmotionforsystemssubjecttoexter- GREEN’S FUNCTION naldrivingforcesthatperiodicallyoscillateintime. One representative example is the Mathieu equation which describes a parametric resonance phenomenon. Here we Besides the Floquet analysis of the Hamiltonian, restrictourselvestoaquantumsystemwhosedynamicsis we can alternatively describe nonequilibrium states determinedbythetime-dependentSchr¨odingerequation, in Green’s-function approach based on the Keldysh formalism.20,21 The approach of the Floquet matrix d proves its own worth when it is used within Green’s- i Ψ(t)=H(t)Ψ(t), (1) dt function formalism. The idea of FGFM was first intro- ducedbyFaisal,8followedbyseveralgroups.9,10,11 Inthis where Ψ(t) is a state vector of the system, and H(t) is section we give another way to define a Floquet matrix thetime-dependentHamiltonian,whichisassumedtobe formofGreen’sfunction,whichweshalluseinthispaper. periodic in t, A Green’s function has two independent arguments of H(t+τ)=H(t), (2) time,tandt′,asdenotedbyG(t,t′). Wedefinevariables t t t′andt (t+t′)/2. Inequilibriumthesystem rel av ≡ − ≡ with a period τ. The Floquet theorem states that there is invariant against continuous time translation, so that exists a solution of Eq. (1) which is an eigenstate of the Green’s functions depend on t and t′ only through t , rel time translation operation t t+τ, implying which enables us to Fourier transform them into func- → tions of ω. However, when the system is driven out of Ψ (t)=e−iεαtu (t) (3) α α equilibrium, they generally depend on both t and t . rel av Since the periodic system that we consider in this pa- with e−iεατ as an eigenvalue of the time translation, α perhasthe discretetimetranslationinvariance[Eq.(2)], as aset ofquantumnumbers,andu (t)=u (t+τ)as a α α it is guaranteed that Green’s function is also invariant periodicfunctionoft. HencewecanFourierexpandu (t) α against t t +τ. For an arbitrary function G(t,t′) as uα(t) = ne−inΩtunα with the frequency Ω = 2π/τ, (not limiatved→toaGvreen’s function) that satisfies the peri- whereun iscalledthenthFloquetmodeofFloquetstate α P odicitycondition,G(t+τ,t′+τ)=G(t,t′),wecandefine (3). We can then Fourier transform Eq. (1) to have the Wigner transformation of G as Hmnunα =(εα+mΩ)umα, (4) ∞ 1 τ/2 Xn Gn(ω)= dtrel τ dtav eiωtrel+inΩtavG(t,t′). Z−∞ Z−τ/2 where (6) We call G (ω) the Wigner representation of the func- 1 τ/2 n H dtei(m−n)ΩtH(t) (5) tion G. Using the Wigner representation, we define the mn ≡ τ Z−τ/2 Floquet matrix form of G as is the Floquet matrix form of the Hamiltonian. The fac- m+n G (ω) G ω+ Ω , (7) tor ε +mΩ appearing on the right-hand side (rhs) of mn m−n α ≡ 2 Eq.(4) is calledquasienergy, whichformsa ladder ofen- (cid:18) (cid:19) ergies with a spacing Ω. Since the Hamiltonian is time and call G (ω) the Floquet representation. Hereafter mn 3 a function with one index n should be understood as a tobe takeninnumericalcalculationsis usuallysmallbe- Wigner representation, while two indices m and n mean cause sufficiently high-order processes should tend to be a Floquet representation. In the Floquet representation, irrelevant when the driving field is not so large, which we use the reduced zone scheme, i.e., the range of ω is also supports the usefulness of FGFM. restricted to the “Brillouin zone” on the frequency axis: Ω/2 < ω Ω/2. We can readily check that definition − ≤ (7)is equivalentto the one givenby Refs.8,9,10,11. The IV. NONINTERACTING ELECTRONS Floquet representationis used during calculations,while the Wigner representationis used when we interpret the Having defined the Floquet representation of Green’s result, since the Wigner representationhas a clear phys- function in Eq. (7), we then compute the Floquet- icalinterpretationthat Gn is the nthoscillatingmode in representedGreen’sfunctionfornoninteractingelectrons. tav of G(t,t′). Although FGFM has been used by several authors to Actually, the Floquet representation Gmn(ω) has a study noninteractingelectronsdrivenoutofequilibrium, one-to-one correspondence with the Wigner representa- there are still further developments yet to be explored. tion Gℓ(ω′). Gℓ(ω′) Gmn(ω): the integers m and n This has motivated us to present an exact and unified → should obey the conditions descriptionofFGFMthatcanbeappliedtogeneralnon- interacting single-band systems in this section. m n=ℓ, (8) In Sec. IVA we provide a general expression of the − Ω m+n Ω FloquetrepresentationofGreen’sfunctionforanysingle- <ω′ Ω , (9) −2 − 2 ≤ 2 band model. Then in Sec. IVB we derive the inverse of Green’s function, which will be used to build DMFT in for ℓ and ω′. There are two consecutive integers k and the Floquet matrix formin Sec. VI. After that, we show k +1 which can be equal to m+n satisfying Eq. (9). several examples of Green’s function for the hypercubic Either k or k+1 is congruent to ℓ modulo 2. Thus m+ (Sec. IVC) and other lattices (Sec. IVD). Through- n is uniquely determined via Eq. (9) and the condition out the paper we restrict our discussionto a single-band m+n m n ℓ (mod 2). Together with Eq. (8), model, and omit spin degrees of freedom for simplicity. ≡ − ≡ m and n are uniquely determined. For such m and n, ω is given by ω′ (m+n)Ω/2. G (ω) G (ω′): ℓ mn ℓ − → and ω′ are uniquely determined by ℓ = m n and ω′ = A. General lattices and fields − ω+(m+n)Ω/2. We can immediately realize the advantage of the Flo- Let ǫk be a band dispersion of the system. We make quet representation in multiplications of two Floquet- the system subject to a homogeneous time-dependent represented functions. As shown in Appendix A, the electric field periodic in t. Here we choose the temporal mapping from G(t,t′) to Gmn(ω) preserves a multipli- gauge or the Hamiltonian gauge in which the scalar po- cation structure, tentialφ=0. Replacingthemomentumkwithk eA(t) [A(t): a vector potential] in ǫk gives the noninte−racting dt′′A(t,t′′)B(t′′,t′)=C(t,t′), Hamiltonian, Z Amℓ(ω)Bℓn(ω)=Cmn(ω). H(t)= ǫk−eA(t)c†kck, (11) ⇔ k Xℓ X As an example, the Floquet representation of the Dyson where c†k and ck are the creation and the annihilation equation is simplified into operators of the electrons, respectively, and we treat the electric field as a classical one. The retarded Green’s (Gk)mn(ω)=(G0k)mn(ω)+ (G0k)mm′(ω)(Σk)m′n′(ω) function for noninteracting electrons reads m′n′ (Gk)n′n(ω),X (10) GRk0(t,t′)=−iθ(t−t′)h[ck(t), c†k(t′)]+i0 × t where G and G0 are, respectively, the full and the non- =−iθ(t−t′)exp −i dt′′ [ǫk−eA(t′′)−µ] , interacting Green’s functions and Σ is the self-energy. (cid:18) Zt′ (1(cid:19)2) Note that each function has the additional 2 2 matrix GR GK × where θ(t) represents the step function, [ , ] is the an- structure, G = in the Keldysh space (with + 0 GA ticommutation relation, is the statistical average thethreelinearly(cid:18)independ(cid:19)entcomponents: theretarded, withrespecttotheinitialhd·e·n·is0itymatrixρ =e−βH(A=0) 0 theadvanced,andtheKeldyshone). Thankstotheusual (where the system is assumed to be in equilibrium with multiplicationruleofthelinearalgebra,onecansolvethe the temperature β−1 at t= ), and µ is the chemical DysonEq.(10)asGk(ω)=[G0k−1(ω) Σk(ω)]−1. Inad- potential of the system. We−ca∞n transform Eq. (12) into − dition, a typical size of a Floquet matrix that is needed the Wigner representationvia Eq.(6), and then into the 4 Floquet representation through Eq. (7). The details of analytically invert Green’s function as GRk0−1(ω)= Λk · the calculation are described in Appendix B. The final Q−1(ω) Λ†,whichcanbeevaluatedexactlyaspresented k k result is inAppe·ndixD. Thederivedexpressionfortheinverseof 1 Green’s function is (GR0) (ω)= k mn Xℓ ω+ℓΩ+µ−(ǫk)0+iη (GRk0−1)mn(ω)=(ω+nΩ+µ+iη)δmn−(ǫk)m−n. (18) π dx π dy ei(m−ℓ)x+i(ℓ−n)y One can also prove Eq. (18) from the Dyson equation × 2π 2π for GR0(t,t′) ina straightforwardmanner. Relation (18) Z−π Z−π k i x means that Green’s function is the kernel of Eq. (4), ×exp −Ω dz [ǫk−eA(z/Ω)−(ǫk)0] , or that the Floquet representation of Green’s function is (cid:18) Zy (cid:19) equivalent to the inverse of the Floquet matrix form of (13) the quasienergy minus the Hamiltonian. One can use Eq. (18) for any single band Hamiltonian with a homo- where η is a positive infinitesimal, and (ǫk)0 is the time- geneous electric field periodic in time. In Secs. IVC and averageddispersion,whichcoincideswiththezerothFlo- IVD,wepresentexamplesofthecalculationsthatutilize quet mode of ǫk, relation (18). π dz (ǫk)m−n = 2π ei(m−n)zǫk−eA(z/Ω). (14) Z−π C. Hypercubic lattice Equation (13) is the general Floquet representation of Green’sfunctionforthe noninteractingsystemdrivenby Asafirstexample,letusconsiderasimplecubiclattice a periodic field. What is notable about expression (13) in d dimensions, whose energy dispersion is given by is thatit canbe decomposedinto well-behavedmatrices. d Let us define two Floquet matrices, ǫsc = 2t cosk , (19) k i − π dx Xi=1 (Λk)mn = ei(m−n)x where t is the hopping and we set the lattice constant 2π Z−π a = 1. For simplicity we assume that the vector poten- ×exp(cid:18)−Ωi Z0xdz [ǫk−eA(z/Ω)−(ǫk)0](cid:19), (15) Atiai(ltA)=(t)Ai(st)p.aSraulblesltittout(i1n,g1k,.i.w.,it1h)kwiitheAea(cth)icnomEqp.o(n1e9n)t, − we have and ǫsc =ǫsccoseA(t)+ǫ¯scsineA(t), (20) 1 k−eA(t) k k (Qk)mn(ω)= δmn, (16) ω+nΩ+µ (ǫk)0+iη where we have defined − d whereΛk isunitaryasshowninAppendix C,andQk(ω) ǫ¯sc = 2t sink , (21) is a diagonalmatrix. The physicalmeaning of these ma- k − i trices will be given later in Sec. V. Here let us just note Xi=1 a strikingly simple decomposition, after Turkowski and Freericks.22 Note that ǫ¯sc has the k odd time-reversal symmetry. Every equation including GRk0(ω)=Λk·Qk(ω)·Λ†k. (17) ǫskc and ǫ¯skc must be consistent against the time-reversal operation. Forinstance,oneusuallyfindsthefactorǫsc+ k where we denote a multiplicationof a Floquet matrix by iǫ¯sc,whichistime-reversaleven,sincetheimaginaryunit k “ ,”andomitFloquetindicesinEq.(17). Thedecompo- · i is time-reversalodd. sition [Eq. (17)] is essentially used to derive the inverse An integral over k is performed through of Green’s function in Sec. IVB. The Floquet representation of the advanced Green’s ρ(ǫ,¯ǫ)= δ(ǫ ǫsc)δ(ǫ¯ ǫ¯sc), (22) k k functionisequaltotheHermitianadjointoftheretarded k − − one: (GA0) (ω)= (GR0†) (ω) =(GR0)∗ (ω). Using X k mn k mn k nm whichiscalledthe jointdensityofstates(JDOS).22 This Eq. (17), we have GAk0(ω)=Λk Q†k(ω) Λ†k. contrasts with the equilibrium cases in which an inte- · · grand depends on k only through ǫsc, so that we can k replace the integral variable k with ǫsc accompanied by k B. Inverse of Green’s function the usual density of states ρ(ǫ). The analytic expression forJDOSinarbitraryddimensionsissummarizedinAp- When one solves a Dyson equation such as Eq. (10) pendix E. Inparticular,ininfinite dimensionsthe JDOS to include effects of interaction, the noninteracting part becomes a Gaussianfunction [Eq. (E2)].22 In the follow- appears as an inverse, GR0−1(ω), rather than GR0(ω) ing, we consider two kinds of electric fields, a dc field k k itself. Usingrelation(17)andtheunitarityofΛk,wecan (Sec. IVC1), and an ac field (Sec. IVC2). 5 1. Hypercubic lattice in a dc field in infinite dimensions. First, we transform Eq. (25) into the Wigner representation, and then integrate it over k Ahomogeneousdcfieldisgivenbythevectorpotential with the JDOS [Eq. (E2)]. After taking the imaginary proportional to time, part, we arrive at eA(t)a=Ωt, (23) A (ω)=δ a δ(ω+ℓΩ+µ). (28) n n,0 ℓ where Xℓ Ω= eEa (24) Here the coefficients in front of the delta functions are − a = e−1/2Ω2I (1/2Ω2), where I (z) is the modified ℓ ℓ ℓ is the Bloch frequency. The system with the field Bessel function of the first kind. Note that the coeffi- [Eq. (23)] fulfills the periodicity condition [Eq. (2)] with cients satisfy the normalization condition, a =1. ℓ ℓ the period 2π/Ω because of the periodic potential of the Equation (28) shows that the local spectral function lattice. WithEq.(13),Green’sfunctionhasthefollowing P has only the n=0 component, indicating that the spec- Floquet representation, tral function evolves into a time-independent function 1 [Eq. (28)] for a sufficiently long time elapsed after the (GRk0)mn(ω)=ei(m−n)θk dc fieldbeganto drivethe system.22 LaterinSec.VIIA, ω+ℓΩ+µ+iη ℓ we shall show that the components with n = 0 vanish X 6 ζk ζk due to a symmetry in the system. J J , (25) m−ℓ ℓ−n × Ω Ω The spectral function on the hypercubic lattice, dis- (cid:18) (cid:19) (cid:18) (cid:19) played in Fig. 1, consists of a set of delta functions with where J (z) is the nth-order Bessel function, and n a spacing Ω, namely, the Wannier-Stark ladder. The width of each peak (approximately effective hopping) is ζk = (ǫskc)2+(ǫ¯skc)2, (26) infinitesimal due to the Bloch oscillation, where an elec- tanθkq=ǫ¯skc/ǫskc. (27) tronisnotfreetoruninalattice,butonlyoscillateswith the frequency Ω. To see how Green’s function (25) behaves, we calculate The inverse of the Floquet-represented Green’s func- thelocalspectralfunctionA (ω)= 1Im (GR0) (ω) tion in this case is given via relation (18) by n −π k k n P (GRk0−1)mn(ω)=(ω+nΩ+µ+iη)δmn 1[(ǫskc+iǫ¯skc)δm−n,1+(ǫskc iǫ¯skc)δm−n,−1] − 2 − 1 =(ω+nΩ+µ+iη)δmn ei(m−n)θkζk(δm−n,1+δm−n,−1). (29) − 2 The Hamiltonian part of Eq. (29) is explicitly written down in a tridiagonal matrix form as ... ... ... 0 ζkeiθk 0 ζke−iθk 0 ζkeiθk 1 ζke−iθk 0 ζkeiθk . (30) 2 ζke−iθk 0 ζkeiθk 0 ζke−iθk 0 ... ... ... As remarked in Sec. II, each component of the Floquet cannotstaystationarybutarealwaysexcitedbythefield. matrix [Eq. (30)] represents a probability amplitude of We also note that the off diagonal components do not a transition from one Floquet mode to another. Note dependonΩ whichis proportionaltothe strengthofthe that for the case of the dc field, Hamiltonian (30) has field. The dependence of Ω is only taken into account no diagonalcomponents,whichmeansthatthe electrons through the quasienergy part of Eq. (29). 6 0.2 Bloch frequency) introduced in Eq. (24) for the dc field. FollowingEq. (13), we derive the Floquet representation of Green’s function as æ æ æ æ æ 1 a{0.1 (GRk0)mn(ω)= ω+ℓΩ+µ ǫscJ (A)+iη æ æ k 0 ℓ − X æ æ π dx π dy ei(m−ℓ)x+i(ℓ−n)y æ æ × 2π 2π 0-æ10æ æ æ æ -5 0 5 æ æ æ æ 1æ0 Z−π Zi−π x { exp dz ǫskc[cos(Asinz) J0(A)] × − Ω { − (cid:18) Zy FIG. 1: The coefficients aℓ of the delta functions in the local +ǫ¯skcsin(Asinz) . (33) } spectral function of the noninteracting electrons on the hy- (cid:19) percubiclattice with thedcfield Ω=0.25 are plotted by the circles. Thedeltafunctionswith thespacing Ωareschemati- Taking the imaginary part of Eq. (33) and integrating cally shown by the solid lines. The broken line is a guide for over k with the JDOS [Eq. (E2)], we obtain the local theeyes. spectral function. We depict it for several A and Ω = 1 inFig.2. Onecanseethatthespectralfunctionhasnar- row peaks at ω = nΩ (n = 0, 1, 2,...) known as the ± ± dynamical Wannier-Stark ladder. The peaks at ω = Ω 2. Hypercubic lattice in an ac field ± correspond to one-photon absorption or emission, and the peaks at ω = 2Ω correspond to two-photon one, ± Let us move on to the case of the ac field on the hy- etc. The width of each peak, or the effective hopping, is percubic lattice. The vector potential is defined by renormalized by the zeroth-order Bessel function J (A) 0 as seen in Eq. (33), and even vanishes making the elec- eA(t)a=AsinΩt, (31) trons completely localized when J (A) = 0. In Fig. 2 0 one finds that the widths of the peaks shrink as A ap- where Ω is the frequency of the ac field, and proaches the first zero (z = 2.40483 ) of J (A) until 0 ··· finally the peaks become the delta functions. This scal- eEa A= (32) ing has been known as dynamical localization since the − Ω proposal by Dunlap and Kenkre.23 is its amplitude divided by the frequency. Although we The inverse of Green’s function in the ac field can be usethesymbolΩ,thisshouldnotbeconfusedwithΩ(the calculated via Eq. (18). The result is (GRk0−1)mn(ω)=(ω+nΩ+µ+iη)δmn−(ǫisk¯ǫcskJcJmm−−nn(A(A)) mm−nn::oevdedn. (34) − The explicit Floquet matrix form of the Hamiltonian part of Eq. (34) reads ... ... ǫkJ0(A) iǫ¯kJ1(A) ǫkJ2(A) iǫ¯kJ3(A) ǫkJ4(A) iǫ¯kJ1(A) ǫkJ0(A) iǫ¯kJ1(A) ǫkJ2(A) iǫ¯kJ3(A) − ǫkJ2(A) iǫ¯kJ1(A) ǫkJ0(A) iǫ¯kJ1(A) ǫkJ2(A) , (35) ··· − ··· −iǫ¯kJ3(A) ǫkJ2(A) −iǫ¯kJ1(A) ǫkJ0(A) iǫ¯kJ1(A) ǫkJ4(A) −iǫ¯kJ3(A) ǫkJ2(A) −iǫ¯kJ1(A) ǫkJ0(A) ... ... where we omit the label “sc” attached to ǫk and ǫ¯k for care because of definition (32) which is singular at Ω = simplicity. Taking the dc limit Ω 0 requires a special 0. Therefore a quantity calculated for a system in the → 7 2 A=0 d α−1 1.5 4 2 ǫfkcc = coskαcoskβ. (36) 2 d(d 1) LΩ1 2.2 − αX=2βX=1 H 2.3 A p In the infinite-dimensional limit (d ) the dispersion of the fcc lattice ǫfcc is related to t→hat∞for the sc lattice k [Eq. (19)] through,24 0 -3 -2 -1 0 1 2 3 1 Ω ǫfkcc =(ǫskc)2 . (37) − 2 FIG. 2: (Color online) The local spectral functions of the UsingEqs.(18)and(37),wederivetheinverseofGreen’s noninteracting electrons on thehypercubiclattice coupled to function for the infinite-dimensional fcc lattice in the dc theac field with Ω=1 and A=0,1.5,2,2.2, and 2.3. field, (GRk0−1)mn(ω)=(ω+nΩ+µ+iη)δmn presence of the ac field with a finite frequency Ω = 0 1 doesnotnecessarilyreproducethe resultcalculatedf6ora ei(m−n)θk[ζk2(δm−n,2+δm−n,−2) − 4 system with the dc field discussed in Sec. IVC1. +2(ζ2 1)δ ]. (38) Let us examine the physical meaning of the Floquet k− mn matrix [Eq. (35)]. The (m,n) component of Hamilto- We see that the Hamiltonian part of the inverse of nian (35) is proportional to J (A). Since J (A) A|m−n| if A is sufficiently smma−lln, the transitimon−nproba∝- Green’sfunction[the secondtermonthe rhsofEq.(38)] is written in a pentadiagonal matrix form. In the same bility m n is proportional to E2|m−n|. For m < n, → way we can obtain Green’s function on the ac field via the process corresponds to stimulated absorption, while Eq. (18). The result reads it corresponds to stimulated emission for m > n. The process of spontaneous emission is not included since we assume that the electric field is classical so that there (GRk0−1)mn(ω)=(ω+nΩ+µ+iη)δmn issumnoptqiounanitsumappflruocptruiaattieonasoflopnhgoatosnthneuminbteernss.ityThoifstahse- 1 ζk2cos(2θk)Jm−n(2A)+(ζk2−1)δmn m−n:even. electric field considered here be so strong as a pulsed − 2(iζk2sin(2θk)Jm−n(2A) m−n:odd laser. (39) WenoticethatthefactorsJ (2A)andJ (0)=δ m−n m−n mn D. Application to other lattices appearinEq.(39)[whileJ (A)appearsonthesclat- m−n tice; seeEq.(34)]. Equations(38)and(39)indicate that So far we have assumed that the vector potential A Green’sfunctionsdependonkonlyviathetwofunctions points to the specific direction (1,1,...,1) in the hyper- ǫskc [Eq. (19)] and ǫ¯skc [Eq. (21)], so that we can integrate cubiclattice. Wecanmoregenerallycalculatetheinverse over k using the JDOS [Eq. (22)]. of the noninteracting retarded Green’s function GR0−1 bymakinguse offormula(18)forarbitrarylattice struc- turesandthevectorpotentials. Thedisadvantageofsuch 2. bcc lattice in infinite dimensions a general case is that, since the k dependence of GR0−1 is not so simple as to be only through ǫsc and ǫ¯sc, the The bcc lattice in d dimensions (d 3) is defined by k k integral over k becomes computationally heavier. the dispersion relation, ≥ Here we note that there is a group of lattice models that give a simple k-dependence of GR0−1 in infinite 8 d α−1β−1 dimensions. Among them are the face-centered cubic ǫbkcc = coskαcoskβcoskγ. −2 d(d 1)(d 2) (fcc) (Sec. IVD1) and body-centered cubic (bcc) − − αX=3βX=2γX=1 (Sec. IVD2) lattices with A parallel to (1,1,...,1) in p (40) We can take the limit d in the same way as in the infinite dimensions. It is instructive to give examples →∞ case of the fcc lattice, and the dispersion converges to other than the hypercubic lattice. 2 1. fcc lattice in infinite dimensions ǫbcc = (ǫsc)3 ǫsc, (41) k k k 3 − The energy dispersion of the fcc lattice generalized to fromwhichwecanderiveGreen’sfunctioninthedcfield, arbitrary d dimensions (d 2) is given by ≥ 8 (GRk0−1)mn(ω)=(ω+nΩ+µ+iη)δmn 1 ei(m−n)θkζk[ζk2(δm−n,3+δm−n,−3)+3(ζk2 2)(δm−n,1+δm−n,−1)]. − 12 − (42) In this case the Hamiltonian part of the inverse of Green’s function [the second term on the rhs of Eq. (42)] becomes a heptadiagonal matrix. Similarly, Green’s function in the ac field is written as (GRk0−1)mn(ω)=(ω+nΩ+µ+iη)δmn− 61(ζikζk[ζ[k2ζk2cosisn(3(3θθkk)J)Jmm−−nn(3(3AA))++33(ζ(ζk2k2−22))csoisnθθkkJJmm−−nn((AA))]] mm−nn::oevdedn, − − (43) where the factor J (3A) newly appears besides One should note that the theorem cannot be applied to m−n J (A). Again, Green’s functions depend on k only a multiband system, where contributions of a geometri- m−n throughǫsc andǫ¯sc,whichmakesanintegraloverk com- cal phase may survive and interband transition could be k k putationally quite efficient with the JDOS [Eq. (22)]. caused by the electric field. The theorem provides a unified description of period- ically driven systems: the original band structure ǫk is V. DIAGONALIZATION OF FLOQUET renormalizedintothetime-averagedone(ǫk)0 duetothe MATRICES field, and the renormalized band splits into its replicas with the spacing Ω. On the (hyper)cubic lattice, the Here we examine how to diagonalize a Floquet matrix dc field changes the band into the Wannier-Stark lad- (ω+nΩ)δ H appearingintheoriginalSchr¨odinger der withthe infinitesimalbandwidth [Fig.3(a)], while in mn mn − Eq. (4) for the noninteracting electrons. We have shown the ac field, the bandwidth scales with the factor J (A) 0 in Sec. IVA that the Floquet matrix form of Green’s [Fig. 3(b)]. functionisdiagonalizedintoQk(ω)bytheunitarytrans- To see how the theorem actually works, let us apply formation Λk [see Eq. (17)]. In Sec. IVB, we have men- it to the d-dimensional fcc lattice model in the dc field tioned that the inverse of Green’s function is equivalent as an example. For convenience, we restrict our discus- totheFloquetmatrix(ω+nΩ+µ+iη)δmn Hmn. Com- sion in the limit d . The Floquet matrix form of biningthesetwofacts,weidentifytheeigenv−aluesandthe the Hamiltonian is g→ive∞n by Eq. (38). According to the eigenvectors of the Floquet matrix. We summarize the theorem, its eigenvalues are equal to the diagonal com- statement below. ponents: ω +nΩ 1(ζ2 1). This means that the dc The eigenvalues of a Floquet matrix Hmn−nΩδmn for fieldmodifiesthee−ne2rgykd−ispersionfromǫfkcc [Eq.(37)]to a single-band noninteracting system subject to a homoge- ˜ǫfcc 1(ζ2 1). Thelatterdispersionisequivalenttothe neous electric field periodic in time are k ≡ 2 k− one for a (d 1)-dimensional hyperplane perpendicular − to (1,1,...,1), the direction of the electric field. We can −(Q−k1)nn(0)=(ǫk)0−nΩ interpretthis fact asfollows: the dc field makesthe elec- =Hnn nΩ (n=0, 1, 2,...), (44) trons localize along the direction of the field due to the − ± ± Blochoscillation,buttheelectronsarefreetomovealong and for each n, the corresponding eigenvector is given by the directions perpendicular to the field because along those directions the force of the dc field does not act on umk−n =(Λk)mn (m=0, 1, 2,...). (45) the electrons. As a result, the motion of the electrons ± ± is confined on the hyperplane, and the energy dispersion This completely solves Floquet matrix problems for a becomes ǫ˜fcc. k single-band system of noninteracting electrons. We note Anotherexampleisprovidedbyanacfield. Ifonesees that the fact that Eq. (44) gives the eigenvalues is es- the Floquet matrix form of the Hamiltonian for the sys- sentially a consequence of Shirley’s relation [Eq. (5) of tem in an ac field on the fcc [Eq. (39)] or bcc [Eq. (43)] Ref. 16]. Since (ǫk)0τ is a dynamical phase, the above lattices,onefindsthatvariouskindsofbandrenormaliza- statementindicatestheabsenceofanadditionalgeomet- tionbesidesthe J (A)scalingonthe simple cubic lattice rical phase (see Appendix D). 0 are implied by the theorem. From Eqs. (17) and (45), the noninteracting Green’s Finally, let us consider electrons in a parabolic band function can be written with the Floquet states as (e.g., a conduction-band bottom in a semiconductor) um−ℓ(un−ℓ)∗ with ǫk = k2/2m∗. If we apply an ac field on it, the (GRk0)mn(ω)= ℓ ω+ℓΩk+µ−k(ǫk)0+iη. (46) d[Fisigp.er3si(ocn)].isreTnhoirsmraelsizueldtsinftroom(ǫkt)h0e=dǫykn+ame2iEca2l/4Fmra∗nΩz2- X 9 grees of freedom except for a representative site i = o, we have the local partition function, Z [G ] = loc 0 [dco][dc†o]eiSloc[co,c†o]. Herethe localactionreadsSloc = dt dt′c†(t)G−1(t,t′)c (t′) + S [c ,c†] [G (t,t′): the R o 0 o int o o 0 Weiss function]. If one ignores the nonlocal fluctua- R R tions,Z andZ giveacommonsite-diagonalself-energy, loc Σ (t,t′)=δ Σ(t,t′). This fact enables us to build a set (a) ij ij of self-consistent closed equations, which can be solved ² with an iterative numerical calculation. Although ne- k glecting the spatial fluctuation is generally an approx- imation, the nonlocal corrections become rigorously ir- ²k (²k) 0 rheolpevpainngt pinartahmeetliemriits socfaliendfinaistet=dimt∗e/n2s√iodns(,t∗w: hfiexreed)t.h28e DMFT isalsoapplicableto nonequilibriumsystemsas recentlystudied29,30,31 basedonnonequilibriumGreen’s- (b) function formalism. Similar to the equilibrium case, one hasself-consistentequationsforGreen’sfunctionandthe self-energy. Now, the present proposal is that if the driven system is periodic in time, one can rewrite the equations in the Floquet matrix form, δZ [G ] loc 0 (G ) (ω)= , (47) loc mn δ(G−1) (ω) 0 nm (G−1) (ω)=(G−1) (ω) Σ (ω), (48) (c) loc mn 0 mn − mn k (G−k1)mn(ω)=(G0k−1)mn(ω)−Σmn(ω), (49) (Gloc)mn(ω)= (Gk)mn(ω). (50) k FIG. 3: The band renormalization ǫk (ǫk)0 on the (hy- X → per)cubic lattice in (a) the dc field and (b) the ac field and To solve Eqs. (47)-(50) self-consistently, we first in- (c) on a parabolic band in the ac field. The solid line repre- put the inverse of the noninteracting Green’s function sents ǫk, while the dashed line respresents (ǫk)0. The circles given by Eq. (18) into Eq. (49). After the initial self- are occupied states movingalong thearrows. energy is properly chosen, the calculation is iterated un- til Green’s function converges. For the lattices discussed in Secs. IVC and IVD, the integralover k in Eq. (50) is Keldysh effect,25,26 which is, in the present formalism, performed via the JDOS [Eq. (22)]. naturallyunderstoodthroughthediagonalizationpicture As remarkedinSec.III,the sizeofthe Floquetmatrix of the Floquet matrix. that needs to be taken in a calculation is usually small ( 5-30, depending on Ω), which, with the analytic ex- ∼ pression of the inverse of Green’s function (18), makes VI. DYNAMICAL MEAN-FIELD THEORY our computational costs dramatically small. WITH THE FLOQUET-GREEN FUNCTION METHOD VII. GAUGE-INVARIANT GREEN’S NowweareinpositiontocombineFGFMwithDMFT FUNCTION for treating interacting systems in external fields. At the basis of DMFT lies the fact that a lattice prob- Before applying our method to a model, we exam- lem of correlated many-body systems can be approx- ine the gauge invariance of Green’s function. Let us imately mapped to a problem of an impurity embed- write the coordinates xν = (t,r) and the vector po- ded into the environment of an effective medium when tential Aν = (φ,A) in the four-vector form. The one ignores spatial fluctuations but takes fully into ac- gauge transformation, A (x) A (x) + ∂ χ(x), puts ν ν ν count on-site dynamical correlation.12,27 The mapping a phase factor to the creationa→nd the annihilation oper- is given as follows: let Z = [dci][dc†i] eiS[ci,c†i] be ators as c†(x) → e−ieχ(x)c†(x) and c(x) → eieχ(x)c(x). a partition function in terms of the original action, Accordingly Green’s function changes as G(x,x′) S = dt dt′ c†(t)G0 −1(t,t′R)c (t′)+ S [c ,c†], eie[χ(x)−χ(x′)]G(x,x′), i.e., the usual Green’s functio→n ij i ij j i int i i where the interaction term is assumed to be a sum is not gauge invariant. It is known32,33 that one can R R P P of the local terms. Integrating out each site’s de- makeGreen’sfunctiongaugeinvariantwithanadditional 10 phase factor as Using Eq. (53), one can evaluate the gauge-invariant Green’s function (52) as x G˜(x,x′)=exp(cid:18)−iZx′ dyν eAν(y)(cid:19)G(x,x′). (51) G˜Rm(ζ,θ,ω)=δm,0 GRn(ζ,θ,ω), (54) n X G˜ depends on the path of the line integral in the expo- where we can see that every mode of Green’s function nent. Here we adopt the conventional straight line con- equally contributes to G˜R. necting x with x′ as the path of the integral. Suppose thatGreen’sfunctiondepends onk onlythroughǫsc and k ǫ¯skc. Theninthetemporalgauge(φ=0)the Wignerrep- B. ac field resentation of the modified Green’s function G˜ becomes Fortheacfield,whichisoneofthekeyquestionsinthe G˜ (ζ,θ,ω)= dω′ 1 τ/2 dt dt present paper, Green’s function has a more complicated m av rel 2π τ θ dependence. To evaluate the gauge-invariant Green’s n Z Z−τ/2 Z X function (52), here we expand the original Green’s func- ×ei(m−n)Ωtav+i(ω−ω′)trel tion with respect to eA(t) in a Taylor series: GRn(ζ,θ+ 1/2 dλ eA,ω) = 1( dλ eA)ℓ ∂ℓGR(ζ,θ,ω). Then Gn ζ,θ+ dλeA(tav+λtrel),ω′ , theWignerrepresenℓtℓa!tionofthegauθge-ninvariantGreen’s × Z−1/2 ! fRunction is exprePssed asR (52) 2Aℓ G˜R(ζ,θ,ω)= dω′ ∂ℓGR(ζ,θ,ω′) where ζ and θ are defined in Eqs. (26) and (27). Note m ℓ!Ω θ n thatG˜ iscalculatedbyshiftingthevariableθ intheorig- Xℓn Z 2 inal Green’s function. This suggests that Green’s func- X(ℓ) Y(ℓ) (ω ω′) , (55) tion integrated in terms of θ is definitely gauge invari- × m−n Ω − (cid:18) (cid:19) ant, so that the local Green’s function kGk(ω) is also where gauge invariant.33 In the following, we derive the gauge- invariantGreen’sfunctionG˜ forthedcfiPeldinSec.VIIA π dx X(ℓ) ei(m−n)xsinℓx and the ac field in Sec. VIIB. m−n ≡ 2π Z−π ℓ 1 ℓ = ( 1)r δ , (56) (2i)ℓ r − m−n,2r−ℓ A. dc field r=0(cid:18) (cid:19) X and To obtain G˜ for the dc field, we first note that the Hamiltonianinthe dc field[Eq.(23)]hasthetime trans- ∞ dx sinx ℓ lation symmetry. If one makes a time translation t Y(ℓ)(k) eikx t+δt,thevectorpotentialchangesasA(t) A(t)+Ωδ→t. ≡Z−∞ 2π (cid:18) x (cid:19) → ℓ Since the change can be absorbed by the gauge trans- 1 ℓ = ( 1)r(k+ℓ 2r)ℓ−1θ(k+ℓ 2r) formation with χ = Ωδt d xi, the Hamiltonian is 2ℓ(ℓ 1)! r − − − − i=1 − r=0(cid:18) (cid:19) invariantagainsttime translation. Then weassumethat X (57) P inthelong-timelimitafterthedcfieldisswitchedonthe retarded Green’s function becomes independent of the for ℓ 1. When ℓ = 0, we have Y(ℓ)(k) = δ(k). The initial correlations. This assumption seems to be valid34 ≥ Floquet representation of Eq. (55) reads as numerically checked in Sec. VIIIA. As a result, a gauge-invariant quantity that is calculated from the re- ∞ 2 A ℓ tardedGreen’sfunctionshouldbeindependentoftheav- G˜R (ζ,θ,ω)=GR (ζ,θ,ω)+ mn mn ℓ!Ω 2i erage time tav. For instance, the local Green’s function Xℓ=1 (cid:18) (cid:19) in the Wigner representation (GRloc)n(ω), which is gauge ℓ ℓ Ω/2 ℓ−1 ω ℓ invariant as shown above, vanishes for n = 0, so that ( 1)r dω′ + dω′ GRloc(t,t′) does not depend on tav. In the s6ame way the ×Xr=0 − (cid:18)r(cid:19) Zω Xk=0 Z−Ω/2 kX=1! self-energy ΣR(ω) also vanishes for n=0. 2 n 6 ∂ℓGR (ζ,θ,ω′)Y(ℓ) (ω ω′) 2k+ℓ . SincetheFloquet-representedself-energyΣRmn(ω)isdi- × θ m−r+k,n+r+k−ℓ Ω − − (cid:18) (cid:19) agonal due to the symmetry, we can identify the θ de- (58) pendence of the Floquet representation of the retarded Green’s function as The derivative with respect to θ in Eqs. (55) and (58) is simplified when the system is on the hypercubic lat- GR (ζ,θ,ω)=ei(m−n)θGR (ζ,θ =0,ω). (53) tice. The Floquet representation of the noninteracting mn mn