Dissipative Floquet Topological Systems
aa r X i v : . [ c ond - m a t . m e s - h a ll ] N ov Dissipative Floquet Topological Systems
Hossein Dehghani , Takashi Oka , and Aditi Mitra Department of Physics, New York University, 4 Washington Place, New York, NY 10003, USA Department of Applied Physics, University of Tokyo, Hongo 7-3-1, Bunkyo, Tokyo 113-8656, Japan (Dated: November 25, 2014)Motivated by recent pump-probe spectroscopies, we study the effect of phonon dissipation andpotential cooling on the nonequilibrium distribution function in a Floquet topological state. Tothis end, we apply a Floquet kinetic equation approach to study two dimensional Dirac fermionsirradiated by a circularly polarized laser, a system which is predicted to be in a laser induced quan-tum Hall state. We find that the initial electron distribution shows an anisotropy with momentumdependent spin textures whose properties are controlled by the switching-on protocol of the laser.The phonons then smoothen this out leading to a non-trivial isotropic nonequilibrium distributionwhich has no memory of the initial state and initial switch-on protocol, and yet is distinct from athermal state. An analytical expression for the distribution at the Dirac point is obtained that isrelevant for observing quantized transport.
PACS numbers: 73.43.-f, 05.70.Ln, 03.65.Vf, 72.80.Vp
I. INTRODUCTION
Recent years have seen the emergence of topologicalstates of matter which is a new way of characterizingmaterials by the geometric properties of the underly-ing band-structure.
These include time reversal (TR)breaking integer quantum Hall systems, TR preserv-ing spin quantum Hall systems or two-dimensional (2D)topological insulators (TIs), 3D TIs, and their stronglyinteracting counterparts. Another intriguing class of sys-tems are those that can show topological behavior onlyunder out of equilibrium conditions, the main candidatebeing the Floquet TIs which arise under periodic driv-ing.
Consider a time periodic Hamiltonian H ( t ) = H ( t + T )where the periodicity may be due to an external ir-radiation by a laser. Then the time-evolution overone period can be written as U ( t + T, t ) = e − iH F T where H F is the Floquet Hamiltonian: an effective time-independent Hamiltonian that captures the stroboscopictime-evolution over one period. Floquet TIs havebeen mainly described by borrowing concepts from equi-librium where the topological properties are extractedby analyzing the spectrum of H F , with the topologi-cal phase showing non-zero Chern numbers and edge-states, though the precise correspondence between theusual equilibrium definition of the Chern number and thenumber of edge-states does not always work. Exper-imentally Floquet TIs have been realized in a photonicsystem which is effectively in equilibrium because the pe-riodicity in time is replaced by a periodicity in position. A Floquet TI has also been realized in a periodically mod-ulated honeycomb optical lattice of fermionic atoms, where in the limit of high frequency of modulations, theFloquet Hamiltonian maps onto the Haldane model. However Floquet TIs are manifestly out of equilibriumand so raise a unique set of questions that do not arisein systems in equilibrium, one of them being the issue ofthe electron distribution function, critical for determin- ing measurable quantities. Obviously the distributionfunction, at least in ideal closed quantum systems, willdepend on how the periodic driving is switched on where any switching-on protocol breaks time-periodicity.In addition the occupation probability will be very sensi-tive to any coupling to an external reservoir.
Oftenreservoir engineering can even produce topological prop-erties absent in the closed system, which in turn re-quires new measures for topological order in open anddissipative systems.
The main aim of this paper is to understand the elec-tron distribution function of Floquet topological systemsby accounting for the initial switching-on protocol of theperiodic drive and accounting for coupling to a reservoirof phonons. We will derive and solve a kinetic equa-tion for the electron distribution function, and show thatthe combined effect of drive and dissipation can stabilizenon-trivial steady-states. We will discuss the signatureof these states on spin and angle resolved photoemission(ARPES).The outline of the paper is as follows, in Section IIthe model is introduced, in Section III physical quanti-ties are calculated for the closed system and for a quenchswitching-on protocol of the laser. In Section IV we gen-eralize to the open system where the electrons are coupledto a phonon reservoir. The rate or kinetic equation ac-counting for inelastic electron-phonon scattering in thepresence of a periodic drive is derived, the results forphysical quantities at steady-state are obtained and com-pared with results for the closed system. Finally in sec-tion V we present our conclusions. Derivation of generalexpressions for the Green’s functions needed for ARPESis given in Appendix A. Analytic results can be obtainedin the vicinity of the Dirac point for both the closed andthe open system, and these are derived in Appendices Band C respectively.
FIG. 1. (color online) Contour plots for the time-averagedspin density P z ( k x , k y ). Left panel: Without phonons andfor a quench. Right panel: At steady-state with phonons. A / Ω = 0 . λ D ph = 0 . , T = 0 . , Ω = 1.
II. MODEL
We study 2D Dirac fermions coupled to an externalcircularly polarized laser, and also coupled to a bath ofphonons. The Hamiltonian is, H = H el + H ph + H c (1)where (setting ~ = 1) H el = X ~k =[ k x ,k y ] ,σ,σ ′ = ↑ , ↓ c † ~kσ h ~k + ~A ( t ) i · ~σ σσ ′ c ~kσ ′ (2) c † ~kσ , c ~kσ are creation, annihilation operators for the Diracfermions whose velocity v = 1, ~σ = [ σ x , σ y ] are the Paulimatrices, ~A = θ ( t ) A [cos(Ω t ) , − sin(Ω t )] is the circu-larly polarized laser which has been suddenly switchedon at time t = 0, we will refer to this switch-on pro-tocol as a quench. This model plays a central role inthe study of Floquet topological states where the circu-larly polarized laser generates a mass term mσ z in theFloquet Hamiltonian H F , with m = A Ω in the high-frequency limit of A / Ω ≪
1. This implies a Hall con-ductivity σ xy = sign( m ) e / h provided a zero temper-ature equilibrium distribution at half-filling is realized. H F is also the continuum limit of the Haldane model which is an example of a TI (Chern insulator), and wasrecently realized experimentally using optical lattices. Gapless surface states of a 3D TI are also modeled by theDirac Hamiltonian, and were recently studied by pump-probe spectroscopy, while laser induced Hall effectand chiral edge states are being experimentally studiedin graphene.
Dissipation affects the electron distribution and thusthe topological signatures. Here we consider dissipationdue to coupling to 2D phonons H ph = X q,i = x,y h ω qi b † qi b qi i (3)For now we do not specify whether we have acoustic oroptical phonons, and hence the particular form of the dispersion ω qi . We will specify this when presenting ourresults. The electron-phonon coupling is H c = X ~k,q,σ,σ ′ c † ~kσ ~A ph ( q ) · ~σ σσ ′ c ~kσ ′ (4) ~A ph ( q ) = (cid:2) λ x,q (cid:0) b † x,q + b x, − q (cid:1) , λ y,q (cid:0) b † y,q + b y, − q (cid:1)(cid:3) (5)Above we neglect phonon induced scattering betweenelectrons with different momenta. This simplifies the ki-netic equation for the electron distribution function con-siderably, without changing the physics, and is a micro-scopic way of accounting for a Caldeira-Leggett typedissipation. In Section III, we will first discuss the physicsin the absence of the phonons H c = 0, but account-ing for the sudden switch-on protocol of the laser, pre-senting results for the steady-state distribution functionand Green’s functions, quantities that are measured inARPES. In Section IV we will address how the resultsget modified due to coupling to phonons. III. RESULTS FOR THE QUENCH AND INTHE ABSENCE OF PHONONS
Suppose that at t ≤
0, there is no external irradiation,and the electrons are in the ground-state, i.e. , all statesbelow the Dirac point are occupied. Thus the wavefunc-tion right before the switching on of the laser is | Ψ in ( t = 0 − ) i = Y ~k | ψ in ,k i| ψ in ,k i = 1 √ (cid:18) − e − iθ k (cid:19) (6)where θ k = arctan ( k y /k x ). The time-evolution afterswitching on the laser is | Ψ( t > i = U el ( t, | Ψ in i (7)where U el ( t, t ′ ) is the time-evolution operator, i dU el ( t, t ′ ) dt = H el ( t ) U el ( t, t ′ ) ; U el ( t, t ) = 1 (8)Since we neglect any spatial dependence of the laserfield, the system stays translationally invariant. Thusthe dynamics is factorizable between different momenta, U el ( t, t ′ ) = Q k U k ( t, t ′ ) so that, | Ψ( t ) i = Q k | ψ k ( t ) i = Q k U k ( t, | ψ in ,k i , where U k ( t, t ′ ) = X α = ± | ψ kα ( t ) ih ψ kα ( t ′ ) | (9) | ψ kα ( t ) i being the exact solution of the Schr¨odinger equa-tion which may be written in terms of the time-periodicFloquet quasi-modes ( | φ kα ( t + T ) i = | φ kα ( t ) i ) and quasi-energies ( ǫ kα ) as follows, | ψ kα ( t ) i = e − iǫ kα t | φ kα ( t ) i [ H el − i∂ t ] | φ kα i = ǫ kα | φ kα i (10)The quasi-energies ǫ kα represent an infinite ladder ofstates where ǫ kα and ǫ kα + m Ω, for m any integer,represent the same physical state corresponding to theFloquet quasi-modes | φ kα ( t ) i and e im Ω t | φ kα ( t ) i . Con-fusion due to this over-counting can be easily avoidedby noting that in all physical quantities, including thematrix elements for electron-phonon scattering that en-ter the kinetic equation, it is always the combination e − iǫ kα t | φ kα ( t ) i = | ψ kα i that appears, where there areonly two distinct states corresponding to α = ± . How-ever while in a typical two-level system, there is only oneenergy-scale corresponding to the level splitting, in thisproblem, a hierarchy of energy scales ǫ k + − ǫ k − + m Ωare possible, although one needs to take care that not allmatrix elements for inelastic processes at these energy-scales may exist. This will be discussed in more detailbelow when we present our results.We can determine the retarded Green’s function g Rσσ ′ ( k, t, t ′ ) = − iθ ( t − t ′ ) h Ψ in | (cid:26) c kσ ( t ) , c † kσ ′ ( t ′ ) (cid:27) | Ψ in i = − iθ ( t − t ′ ) U k,σσ ′ ( t, t ′ ) (11)and the lesser Green’s function, g <σσ ′ ( k, t, t ′ ) = − i h Ψ in | c † kσ ( t ) c kσ ′ ( t ′ ) | Ψ in i = − i X σ ,σ h Ψ in | c † kσ c kσ | Ψ in i U k,σ σ (0 , t ) U k,σ ′ σ ( t ′ , g R does not depend on the occupation probability(by not depending on the initial state), g < depends on it.We perform a Fourier transformation of the Green’s func-tions g ( k, t, t ′ ) with respect to the time-difference t − t ′ thus moving to the frequency ω space, and all through-out we present results after time-averaging over the meantime T m = ( t + t ′ ) /
2. Thus in what follows, whenever wedenote quantities by the arguments k, ω alone, it shouldbe understood that an averaging over mean time T m has already been performed. General expressions for theGreen’s functions are presented in Appendix A where theaveraging procedure over the mean time is also explained.We refer to ig <σσ ( k, ω ) as the spin resolved ARPESspectrum, a key quantity in this work that can be di-rectly probed in experiments. Note that results for thespectral density A = Im (cid:2) g R (cid:3) have been discussed else-where , however our results for g < even in the absenceof phonons are new. We note that number conservation,absence of momentum mixing, and the fact that we are athalf-filling imply the sum rule R ( dω/ π ) i P σ g <σσ ( k, ω ) =1. The results for the momentum dependent spin-density P z ( k, T m ) = i P σ σg <σσ ( k, T m , T m ) after averaging over T m is shown as a contour plot in the left panel of Fig. 1as well as along the line k y = 0 in Fig. 2. The circu-larly polarized laser induces a strongly momentum de-pendent spin density which shows oscillations each timethe condition for a photon induced resonance between -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 k x -0.8-0.6-0.4-0.200.20.4 P z quench, no phononsphonons, T=0.01phonons, T=1.0 FIG. 2. Spin density P z ( k x , k y = 0) for A / Ω = 0 . , Ω =1 . , λ D ph =0 .
1Ω for three different cases: for the quench withno phonons, steady-state with phonons at temperature T =0 . T = Ω. the Dirac bands | k | ≃ n Ω /
2, where n is an integer, isobeyed. Further, the density is also anisotropic in mo-mentum space. The spin averaged ARPES spectrum ig < tot ( k, ω ) = i P σ g <σσ ( k, ω ) is plotted as an intensity plotin Fig. 3 and its momentum slices in Fig. 4. Note that thedelta-functions are given an artificial broadening which issuch that the heights of the peaks in Fig. 4 equal the pref-actor of the delta function. In other words, Bπδ ( ω − ǫ k )has been approximated by Bγ γ γ +( ω − ǫ k ) with the broad-ening γ arbitrary and chosen so that the plots are visible.Moreover we plot γg (or γG with phonons), so that theheight of the peaks in Fig. 4 equals the prefactor of the δ -function B .The ARPES spectrum clearly shows the appearanceof Floquet bands. Without phonons, the system is free,and the electron distribution is given by the overlap |h φ kα = ± (0) | ψ in , k i| . This is a highly non-thermal statethat retains memory of the initial state | ψ in , k i , and isnot expected to thermalize. Just like the spin resolveddensity, the total density in Fig. 3 shows a clear asymme-try under k x → − k x , where this particular anisotropy isdetermined by the phase of the AC field at t = 0 + . Notethat in our case, initially the gauge field ~A ( t = 0 + ) =[ A ,
0] is entirely along the ˆ x direction.The anisotropy can be understood analytically at k =0(see Appendix B for details), P z ( k = 0 , θ k )= − A Ω∆ cos θ k ;∆ = q A + Ω (13)where θ k is the angle along which k = 0 is approached.This has the same anisotropy as the left panel in Fig.1.For the lesser Green’s function at k =0 we obtain, ig <σσ ( k = 0 , θ k , ω ) = 2 π X α = ± ρ quench k =0 ,αα × (cid:20)(cid:18) ∆ + ασ Ω2∆ (cid:19) δ (cid:18) ω + σ ( α ∆ + Ω)2 (cid:19)(cid:21) ,ρ quench k =0 ,αα = |h φ k =0 ,α (0) | ψ in ,k =0 i| = 12 (cid:18) − αA ∆ cos θ k (cid:19) (14)above σ = + / − for ↑ / ↓ . The analytic expression for g < shows that for k = 0, there are exactly four resonancesfor inelastic scattering, where the two resonances for spin σ occur at ω = − σ (Ω ± ∆) /
2. Naively one would haveexpected infinite number of resonances ǫ k + − ǫ k − + m Ωfor integer m . The fact that at k = 0 there are so few isdue to vanishing matrix elements alluded to above. As k increases, more and more resonances appear, howeverthey are very rapidly suppressed for large | m | .The above location of the resonances also shows thatthe circularly polarized field acts as an effective magneticfield along ˆ z , splitting the energies of the up and downspin electrons. In particular in the high frequency ( A ≪ Ω) limit, the lowest energy excitation is ∆ − Ω ≃ A / Ωand involves flipping a spin from ↓ to ↑ . However, thisis still not a typical two level system, as for k = 0, thereare two energy scales for energy absorption ( ω >
0) (andmore for k = 0), rather than just one energy-scale forenergy absorption encountered in a conventional two levelsystem.The analytic results also show that the weights are farfrom thermal, where by thermal we imply resonances ofthe form δ ( ω − ǫ k ) n F ( ω ), n F being the Fermi distribu-tion function at some temperature T . Rather the heightof the resonances are proportional to amplitude squareof the overlap between the initial wavefunction corre-sponding to the ground state of the Dirac model, and thewavefunctions | ψ kα i . Note that the appearance of onlya couple of Floquet bands, and momentum anisotropy isconsistent with experimental observations. IV. RESULTS IN THE PRESENCE OFPHONONS
The above results for the time-averaged distributionfunctions after a quench are exact and will not evolvein time. However if we turn on the electron-phonon cou-pling, inelastic scattering will cause the distribution func-tions to relax, we now study how this happens, and whatis the resulting steady-state. We first briefly outline thederivation of the kinetic or rate equation in the presenceof phonons within the Floquet formalism (see for gen-eral discussions). Let W ( t ) be the density matrix obeying dW ( t ) dt = − i [ H, W ( t )] (15) FIG. 3. (color online) Intensity plot of the spin-averagedARPES spectrum iG < tot ( k, ω ) / k y = 10 − for the quenchwith no phonons (upper panels) and at steady-state withphonons at temperature T = 0 . ,
1Ω (middle and lowerpanels). A / Ω = 0 . , λ D ph = 0 . , Ω = 1 . It is convenient to be in the interaction representation, W I ( t ) = e iH ph t U † el ( t, W ( t ) U el ( t, e − iH ph t . To O ( H c ),the density matrix obeys the following equation of motion dW I dt = − i [ H c,I ( t ) , W I ( t )] − Z tt dt ′ [ H c,I ( t ) , [ H c,I ( t ′ ) , W I ( t ′ )]] (16)where H c,I is in the interaction representation. Weassume that at the initial time t , the electrons andphonons are uncoupled so that W ( t ) = W el0 ( t ) ⊗ W ph ( t ), and that initially the electrons are in the state | Ψ( t ) i described in Section III, while the phonons arein thermal equilibrium at temperature T . This is justi-fied because phonon dynamics is much slower than elec-tron dynamics. Thus the quench state of Section III canbe achieved within femto-second time-scales, while, thephonons do not affect the system until pico-second time-scales.Thus, W el0 ( t ) = | Ψ( t ) ih Ψ( t ) | = Y k W el k, (17)where W el k, ( t ) = X α,β = ± e − i ( ǫ kα − ǫ kβ ) t | φ kα ( t ) ih φ kβ ( t ) | ρ quench k,αβ (18)with ρ quench k,αβ = h φ kα (0) | ψ in ,k ih ψ in ,k | φ kβ (0) i (19)Defining the electron reduced density matrix as the oneobtained from tracing over the phonons, W el = Tr ph W ,and noting that H c being linear in the phonon operators,the trace vanishes, we need to solve, dW el I dt = − Tr ph Z tt dt ′ [ H c,I ( t ) , [ H c,I ( t ′ ) , W I ( t ′ )]] (20)We assume that the phonons are an ideal reservoir andstay in equilibrium. In that case W I ( t ) = W el I ( t ) ⊗ e − H ph /T / Tr (cid:2) e − H ph /T (cid:3) (we set k B = 1).The most general form of the reduced density matrixfor the electrons is W el I ( t ) = Y k X αβ ρ k,αβ ( t ) | φ k,α ( t ) ih φ k,β ( t ) | (21)where in the absence of phonons, ρ k,αβ = ρ quench k,αβ and aretime-independent in the interaction representation. Thelast remaining assumption is to identify the slow and fastvariables, which allows one to make the Markov approx-imation. We write ρ k,αβ ( t ) = P m =int e im Ω t ρ ( m ) k,αβ ( t )where in what follows we assume that ρ ( m ) k,αβ ( t ) are slowlyvarying on time scales of the period of the AC field andthe relevant phonon frequencies. In addition we onlystudy the diagonal components of ρ ( m ) k,αα , which after theMarkov approximation, obey the rate equation h ˙ ρ ( m ) k,αα ( t ) + im Ω ρ ( m ) k,αα ( t ) i = − X m ′ ,β = ± L m,m ′ k,αβ ρ ( m − m ′ ) k,ββ ( t )(22)The initial condition we will consider corresponds to aquench switch on protocol for the laser ρ ( m ) k,αα ( t = 0) = δ m =0 ρ quench k,αα , with the in-scattering and out-scatteringrates L m,m ′ k,αβ given in Appendix C.Since the rate equation is a weak-coupling quasi-classical approximation in the electron-phonon coupling,the position of the resonances in the spectral density arenot modified, and thus even with phonons, g R is un-changed. The phonons strongly modify the steady-statelesser Green’s function because the distribution functionof the electrons is changed due to inelastic scattering withphonons. In the numerical solutions for the rate equa-tion we assume optical phonons with a uniform phonondensity of states D ph , with a broad band-width so thatinelastic scattering is always possible. We also assumean isotropic electron-phonon coupling λ x = λ y = λ . Theresults can easily be generalized to optical phonons withnarrow band-widths, as for frequencies below or abovethe optical phonon frequencies, the distribution functionwill remain unchanged, and will be given by that for thequench.The time-evolution of the density matrix from aquench-type initial state is shown in Fig. 5, where therate for reaching steady-state is set by the strength ofthe electron-phonon coupling λ D ph . In what follows, wepresent results for G < at long times when a steady-statehas been reached. The solution of the rate equations in Fig. 5 shows that the steady-state is characterizedby some oscillations with time (controlled by λ D ph ),and our results are presented after a time-averaging of ρ k,αα ( t ) = ρ ss k,αα over several cycles. After this time-averaging, the steady-state lesser Green’s function in thepresence of phonons is given by, G <σσ ′ ( k, t, t ′ ) = − i X α = ± (cid:20) ρ ss k,αα h φ k,α (0) | c † kσ ( t ) c kσ ′ ( t ′ ) | φ k,α (0) i (cid:21) (23)where c kσ ( t ) = P σ ′ U kσσ ′ ( t, c kσ ′ (0). Note that due tothe laser field, G < is not time-translationally invariant,and so we average over the mean time ( t + t ′ ) / k = 0, L m,m ′ k =0 ,αβ = δ m,m ′ L mk =0 ,αβ , sothat again analytic results are possible. Here we find forthe spin-density at k = 0, P z ( k = 0; H c = 0) = − (cid:0) ∆ + Ω (cid:1) / ∆ P α = ± (∆ − α Ω) (cid:26) N (∆ + α Ω) (cid:27) (24)where N ( x ) is the Bose distribution function, while iG <σσ ( k = 0 , ω ; H c = 0) = 2 π X α = ± ρ ss k =0 ,αα × (cid:20)(cid:18) ∆ + ασ Ω2∆ (cid:19) δ (cid:18) ω + σ ( α ∆ + Ω)2 (cid:19)(cid:21) (25)where ρ ss k =0 , ++ = P β = ± (∆ − β Ω) N (∆ + β Ω) P α = ± (∆ − α Ω) (1 + 2 N (∆ + α Ω))(26)with P α = ± ρ ss k =0 ,αα = 1. Note that the above results at k = 0 are isotropic in being independent of the angle θ k .Thus the coupling to phonons makes the electrons losememory of the initial state as well as the initial switch-onprotocol for the laser. This results in a symmetric distri-bution of the density in momentum space. This is alsoclearly seen in the contour plot of Fig. 1. Fig. 2 showsthat the spin-density still retains oscillations at momenta k for which the photon frequencies are resonant with theDirac bands, however the magnitude of the oscillationsdecay with increasing temperature of the phonon bath,with the spin-density P z ( k ) approaching zero as the tem-perature increases.The spin averaged ARPES spectrum iG < tot in thesteady-state with phonons is shown as an intensity plotin the middle and lower panels in Fig. 3 and along somemomentum slices in Fig. 4. One finds that as the tem-perature of the phonon bath decreases, the magnitude ofthe resonances at positive frequencies decrease and theones at negative frequencies increase, maintaining thesum rule. While this is also the expected result froma simple thermal Green’s function where the weights ofthe resonances are δ ( ω − ǫ k ) n F ( ǫ k ), yet note that the i g < t o t / i G < t o t / -2 -1 0 1 2 ω/Ω i G < t o t / -1 0 1 ω/Ω -2 -1 0 1 2 ω/Ω k x =0quench k x =0.25quench k x =0.5quenchk x =0phononsT=0.01k x =0phononsT=1 k x =0.25phononsT=1k x =0.25phononsT=0.01 k x =0.5phononsT=0.01k x =0.5phononsT=1 FIG. 4. iG < tot ( k, ω ) for the quench with no phonons (up-per panels) and at steady-state with phonons at temperature T = 0 . ,
1Ω (middle and lower panels) for k y = 10 − and k x = 0 . , . , . A / Ω = 0 . , λ D ph = 0 . , Ω = 1 .
0. Nor-malization such that the peak heights equal the prefactor ofthe δ -functions. precise weights in steady-state are not thermal. This canalso be clearly seen in the analytic solution for k = 0.In particular Eq. (24) implies that in the high frequencylimit P z ( k = 0 , A ≪ Ω) → tanh (cid:2) A / Ω T (cid:3) [1 + ( A / Ω ) tanh ( A / Ω T ) coth (Ω /T )] (27)In this high-frequency limit, the Floquet Hamiltonian is H F ≃ σ x k x + σ y k y + σ z A / Ω, so a naive guess would bethat the thermalized state should have a magnetizationof tanh ( h z / T ) where h z = 2 A / Ω. The result for P z shows deviations from this guess at O ( A / Ω ). Thus, thepresence of the AC drive causes the electrons to reach anonequilibrium steady-state even when the phonon reser-voir to which the electrons are coupled are themselvesalways in thermal equilibrium. Fig. 4 also shows thatas k approaches the photon induced resonance condition | k | ∼ n Ω /
2, the effective temperature is higher, as morefrequencies are excited. This result is clearly reflectedin Fig. 3 (central panel) where even when the phonontemperature is very low, the avoided crossings are char-acterized by a high population density.
V. CONCLUSIONS
In summary we have studied the electron distributionin a Floquet topological system under two circumstances,one is for the closed system, where the results are verysensitive to how the AC field has been switched on, show-ing highly anisotropic distribution functions, the secondis for the open system where the electrons are coupledto a reservoir of phonons. While coupling to phonons Ω t ρ ( t ) ρ -- (k x = 10 -4 ) ρ ++ (k x = 10 -4 ) ρ -- (k x = 0.25) ρ ++ (k x = 0.25) ρ -- (k x = 0.5) ρ ++ (k x = 0.5) FIG. 5. Time-evolution of ρ k,αα = ± from an initial state cor-responding to a quench for k y =0 and k x = 10 − , . , . A / Ω = 0 . , λ D ph = 0 . , T = 0 . , Ω = 1 . causes the system to lose memory of its initial state, yetthe presence of the drive gives rise to non-trivial nonequi-librium steady-states observable in ARPES. Since elec-tron dynamics is much faster than phonon dynamics, theresults for the quench should be observable on short ∼ femto-second time-scales, while the phonons will start re-laxing the system on much longer time-scales. The rateequations show that the system will eventually reach anonequilibrium steady-state with the phonons on time-scales that are inverse of the effective electron-phononcoupling O ( λ D ph ) which is a highly material dependentparameter, and in realistic materials suggests time-scalesof the order of pico-seconds.An important open question is to understand transportphenomena such as the Hall conductance. Since the Hallresponse is dominated by the behavior near k = 0 wherethe Berry curvature is peaked, our results imply that theanisotropic distribution of the closed system will causesignificant deviation from the quantum limit. On theother hand coupling to low temperature phonons inducescooling of Floquet states. The cooling works efficientlynear the Dirac point, which could help the system to ap-proach the quantum limit. However, near resonant points(energy difference ∼ n Ω), our results also show that theeffective temperature stays high due to photo-carriers.Thus how close the system is to the quantum limit willbe a competition between the contribution to the Hall-conductance at the Dirac points, and the role of these ex-cited photo-carriers. A quantitative treatment requiresbeing on the lattice, as in the continuum, the Berry-curvature shows very sharp peaks at the resonances, andbecomes mathematically ill-defined. However it is clearthat to reach the quantum limit, in addition to having alow temperature bath, it would also be helpful to be inthe high frequency regime where Ω is greater than theband-width so as to suppress excited photo-carriers, aregime which is non-existent in the continuum due to theunbounded energy dispersion.
Acknowledgments:
This work was supported by USDepartment of Energy (DOE-BES) under Award No.DE-SC0010821 (HD and AM), and partially by the Si-mons Foundation (academic year support for AM).
Appendix A: Green’s functions
In this section we highlight how the Green’s functionsdefined in Eqns. (11) and (12) can be obtained. Since thequasi-modes are periodic in time, we may write them as, | φ kα ( t ) i = X m ∈ int e im Ω t (cid:18) α kmα β kmα (cid:19) (A1)Thus the time-evolution operator becomes, U k ( t = T m + τ / , t ′ = T m − τ /
2) = X α = ± ,m,m ′ e − iǫ kα τ + m + m ′ Ω τ + i ( m − m ′ )Ω T m × (cid:18) α kmα β kmα (cid:19) (cid:0) α ∗ km ′ α β ∗ km ′ α (cid:1) (A2)Averaging over T m , U k ( τ ) = X α = ± ,m e − iǫ kα τ + m Ω τ (cid:18) α kmα β kmα (cid:19) (cid:0) α ∗ kmα β ∗ kmα (cid:1) (A3)so that on Fourier transforming with respect to the timedifference τ , the retarded Green’s function becomes g R ( k, ω ) = X α,m ω − ( ǫ kα − m Ω) + iδ × (cid:18) α kmα β kmα (cid:19) (cid:0) α ∗ kmα β ∗ kmα (cid:1) (A4)For k = 0, analytic expressions for g R may be obtainedand these are presented in Eq. (B24),For the lesser Green’s functions, using Eq. (12), wehave, g < ↑↑ ( k, t, t ′ ) = − i X α,β = ± ,m,m ′ ,n,n ′ e − iǫ kβ t ′ + iǫ kα t + in Ω t ′ − im Ω t (cid:20) α knβ α ∗ kn ′ β α km ′ α α ∗ kmα h c † k ↑ c k ↑ i + α knβ β ∗ kn ′ β β km ′ α α ∗ kmα h c † k ↓ c k ↓ i + α knβ α ∗ kn ′ β β km ′ α α ∗ kmα h c † k ↓ c k ↑ i + α knβ β ∗ kn ′ β α km ′ α α ∗ kmα h c † k ↑ c k ↓ i (cid:21) (A5) and, g < ↓↓ ( k, t, t ′ ) = − i X α,β = ± ,m,m ′ ,n,n ′ e − iǫ kβ t ′ + iǫ kα t + in Ω t ′ − im Ω t (cid:20) β knβ α ∗ kn ′ β α km ′ α β ∗ kmα h c † k ↑ c k ↑ i + β knβ β ∗ kn ′ β β km ′ α β ∗ kmα h c † k ↓ c k ↓ i + β knβ α ∗ kn ′ β β km ′ α β ∗ kmα h c † k ↓ c k ↑ i + β knβ β ∗ kn ′ β α km ′ α β ∗ kmα h c † k ↑ c k ↓ i (cid:21) (A6)where h c † kσ c kσ ′ i = h ψ in ,k | c † kσ c kσ ′ | ψ in ,k i for the closedsystem with a quench switch-on protocol. For the opensystem in steady-state, h c † kσ c kσ ′ i is the average with re-spect to the steady-state reduced density matrix of theelectrons, which is in turn obtained from solving a kineticequation.Time-averaging over the mean time T m = ( t + t ′ ) / α = β, m = n , so that g < ↑↑ ( k, τ = t ′ − t ) = − i X α = ± ,n e − i [ ǫ kα − n Ω] τ | α knα | ρ k,αα g < ↓↓ ( k, τ = t ′ − t ) = − i X α = ± ,n e − i [ ǫ kα − n Ω] τ | β knα | ρ k,αα (A7)Above ρ k,αα = |h φ k,α (0) | ψ in ,k i| = ρ quench k,αα for the quenchin the closed system, while it is obtained from a kineticequation for the open system. For the latter, inelas-tic scattering causes ρ k,αα to evolve in time, and theMarkov approximation that we will employ requires thatthis time-dependence is slow as compared to all othertime-scales. For the open system, we will then presentresults for the Green’s functions only at long times, wherea steady-state has been reached, where the density ma-trix is replaced by its steady-state value ρ k,αα = ρ ssk,αα .Sometimes, some slow residual oscillations such as thoseshown in Fig 5 persist even at long times, in this casesuch slow oscillations will also be averaged over.Fourier transforming, ig < ↑↑ ( k, ω ) = 2 π X nα δ ( ω − [ ǫ kα − n Ω]) | α knα | ρ k,αα ig < ↓↓ ( k, ω ) = 2 π X nα δ ( ω − [ ǫ kα − n Ω]) | β knα | ρ k,αα (A8)Analytic expressions for the lesser Green’s function for k = 0 are given in Eq. (14) for the quench and in Eq. (25)for the steady-state with phonons. Appendix B: Analytic solution at k = 0 for thequench (no phonons) Let us consider the solution of H el when k = 0. In thiscase, the quasi-modes | φ α i (we suppress the k = 0 label)obey the equation, H el , F ( k = 0) | φ α i = ǫ α | φ α i (B1) H el , F ( k = 0) = ~A · ~σ − i∂ t (B2) | φ α i = (cid:18) φ ↑ α φ ↓ α (cid:19) (B3)where ~A = A (cos Ω t, − sin Ω t ), so that ~A · ~σ = A (cid:18) e i Ω t e − i Ω t (cid:19) . Thus, the φ ↑ , ↓ α obey the coupledequation − i∂ t φ ↑ α + A e i Ω t φ ↓ α = ǫ α φ ↑ α (B4) − i∂ t φ ↓ α + A e − i Ω t φ ↑ α = ǫ α φ ↓ α (B5)Substituting for φ ↓ α = e − i Ω t A [ ǫ α φ ↑ α + i∂ t φ ↑ α ] (B6)into the second equation above gives, ∂ t φ ↑ α − i [2 ǫ α + Ω] ∂ t φ ↑ α + (cid:0) A − Ω ǫ α − ǫ α (cid:1) φ ↑ α = 0(B7)Writing φ ↑± = d ↑± e iλ ± t , one obtains λ ∓ = Ω2 + ǫ ∓ ± ∆2 where ∆ = q A + Ω (B8)Since φ ↑ , ↓ α ( t + T ) = φ ↑ , ↓ α ( t ), λ = m Ω, where m is aninteger. Thus Eq. (B6) gives, φ ↓ α = d ↓ α e i ( m − t ; φ ↑ α = d ↑ α e im Ω t (B9)with ǫ ± = (cid:18) m − (cid:19) Ω ± ∆2 (B10) d ↓± d ↑± = − Ω ± ∆2 A (B11)Thus, d ↑± = √ A p ∆ (∆ ∓ Ω) ; d ↓± = ± √ r ∓ Ω∆ (B12) | φ ± ( t ) i = e im Ω t (cid:18) d ↑± e − i Ω t d ↓± (cid:19) (B13)Note that while there are infinite possible ways to choosethe quasi-modes and the corresponding quasi-energies,where the quasi-energies are related by shifts by integermultiples of the frequency Ω, this degeneracy is absentin the wavefunctions corresponding to the exact solutionsof the Schr¨odinger equation, | Ψ α ( t ) i = e − iǫ α t | φ α ( t ) i . Inparticular the wavefunctions are | Ψ + ( t ) i = e i Ω t/ − i ∆ t/ √ A √ ∆(∆ − Ω) e − i Ω t √ q − Ω∆ (B14) | Ψ − ( t ) i = e i Ω t/ i ∆ t/ √ A √ ∆(∆+Ω) − e − i Ω t √ q Ω∆ (B15) One may also construct the time-evolution operator, U k =0 ( t, t ′ ) = X α = ± e − iǫ α ( t − t ′ ) | φ α ( t ) ih φ α ( t ′ ) | (B16)= X α = ± e − i ( − Ω+ α ∆2 ) ( t − t ′ ) (cid:18) d ↑ α e − i Ω t d ↓ α (cid:19) (cid:0) d ↑ α e i Ω t ′ d ↓ α (cid:1) (B17)If the state just before switching on the AC field is theground state of the Dirac fermions, | ψ in i = 1 √ (cid:18) − e − iθ k (cid:19) (B18)then the wavefunction after the sudden switch-on of theAC field is given by | Ψ( t ) i = U k =0 ( t, | ψ in i = X α = ± C − α p ∆(∆ − α Ω) √ A | Ψ α ( t ) i (B19)where C + = − [ ω − ψ ↑ (0) + A ψ ↓ (0)]( ω + − ω − ) (B20) C − = [ ω + ψ ↑ (0) + A ψ ↓ (0)]( ω + − ω − ) (B21)with ω ± = Ω ± ∆2 (B22) ψ ↑ (0) = − e − iθ k / √ ψ ↓ (0) = 1 / √ | Ψ( t ) i and the time evolution op-erator U ( t, t ′ ) are known, one may compute all the single-time and two-time averages discussed in the main text.Using the above, and Eq. (A4), the expression for theretarded Green’s function is, g R ( k = 0 , ω ) = X α = ± ω − (cid:0) − Ω2 + α ∆2 (cid:1) + iδ (cid:18) d ↑ α
00 0 (cid:19) + X α = ± ω − (cid:0) Ω2 + α ∆2 (cid:1) + iδ (cid:18) d ↓ α (cid:19) (B24)Using Eq. (A8), the lesser Green’s function is given inEq. (14). Appendix C: Rate equations for general k and exactsolution at k = 0 The rate equations after the Markov approximation arefound to be (below N q = N ( ω q ) is the Bose distributionfunction) h ˙ ρ ( m ) k,αα ( t ) + im Ω ρ ( m ) k,αα i = − X q,i = x,y,β = ± ,n ,n (cid:20) πλ iq (cid:18) ǫ i ¯ i (cid:26) C n kαβ C n kβα + C n kαβ C n kβα (cid:27) + C n kαβ C n kβα + C n kαβ C n kβα (cid:19)(cid:21) × (cid:20)(cid:26) (1 + N q ) δ ( ǫ kβ − ǫ kα + ( m − n )Ω + ω qi ) + N q δ ( ǫ kβ − ǫ kα + ( m − n )Ω − ω qi ) (cid:27) ρ ( m − n − n ) k,αα ( t ) − (cid:26) (1 + N q ) δ ( ǫ kβ − ǫ kα + ( m − n )Ω − ω qi ) + N q δ ( ǫ kβ − ǫ kα + ( m − n )Ω + ω qi ) (cid:27) ρ ( m − n − n ) k,ββ ( t )+ (cid:26) (1 + N q ) δ ( ǫ kβ − ǫ kα − ( m − n )Ω + ω qi ) + N q δ ( ǫ kβ − ǫ kα − ( m − n )Ω − ω qi ) (cid:27) ρ ( m − n − n ) k,αα ( t ) − (cid:26) (1 + N q ) δ ( ǫ kβ − ǫ kα − ( m − n )Ω − ω qi ) + N q δ ( ǫ kβ − ǫ kα − ( m − n )Ω + ω qi ) (cid:27) ρ ( m − n − n ) k,ββ ( t ) (cid:21) (C1)where ǫ x ¯ x = 1 , ǫ y ¯ y = −
1, and h φ kα ( t ) | c † k ↑ c k ↓ | φ kβ ( t ) i = X n e in Ω t C n kαβ (C2) h φ kα ( t ) | c † k ↓ c k ↑ | φ kβ ( t ) i = X n e in Ω t C n kαβ (C3)
1. Analytic results for the rate equation at k = 0 At k = 0, the exact expressions for the quasi-modescan be used to show that h φ α ( t ) | c † k =0 , ↑ c k =0 , ↓ | φ β ( t ) i = d ↑ α d ↓ β e − i Ω t (C4) h φ α ( t ) | c † k =0 , ↓ c k =0 , ↑ | φ β ( t ) i = d ↓ α d ↑ β e i Ω t (C5)Thus, the matrix elements entering the rate equation be-come, C ( n )1++ = A ∆ δ n = − ; C ( n )1 −− = − A ∆ δ n = − C ( n )1+ − = − (cid:18) (cid:19) δ n = − ; C ( n )1 − + = 12 (cid:18) − Ω∆ (cid:19) δ n = − C ( n )2++ = A ∆ δ n =1 ; C ( n )2 −− = − A ∆ δ n =1 C ( n )2+ − = 12 (cid:18) − Ω∆ (cid:19) δ n =1 ; C ( n )2 − + = − (cid:18) (cid:19) δ n =1 Let us assume λ xq = λ yq . In this case for k = 0, n + n =0 in the rate equations. So for k = 0, the rate equationssimplify to ∂ t ρ ( m ) k =0 , ++ ρ ( m ) k =0 , −− ! + im Ω ρ ( m ) k =0 , ++ ρ ( m ) k =0 , −− ! = L ( m ) k =0 , ++ L ( m ) k =0 , + − L ( m ) k =0 , − + L ( m ) k =0 , −− ! ρ ( m ) k =0 , ++ ρ ( m ) k =0 , −− ! (C6)where L ( m ) k =0 , ++ = − L ( m ) k =0 , − + ; L ( m ) k =0 , + − = − L ( m ) k =0 , −− . Wenow make the assumption of a uniform phonon density of states D ph so that the rates are, L ( m ) k =0 , ++ = − πλ D ph (cid:18) (cid:19) (cid:20) { N (∆ − Ω − m Ω) } θ (∆ − Ω − m Ω)+ { N (∆ − Ω + m Ω) } θ (∆ − Ω + m Ω)+ N ( − ∆ + Ω − m Ω) θ ( − ∆ + Ω − m Ω)+ N ( − ∆ + Ω + m Ω) θ ( − ∆ + Ω + m Ω) (cid:21) − πλ D ph (cid:18) − Ω∆ (cid:19) (cid:20) { N (∆ + Ω − m Ω) } θ (∆ + Ω − m Ω)+ { N (∆ + Ω + m Ω) } θ (∆ + Ω + m Ω)+ N ( − ∆ − Ω − m Ω) θ ( − ∆ − Ω − m Ω)+ N ( − ∆ − Ω + m Ω) θ ( − ∆ − Ω + m Ω) (cid:21) (C7)and, L ( m ) k =0 , −− = − πλ D ph (cid:18) − Ω∆ (cid:19) (cid:20) { N ( − ∆ − Ω − m Ω) } θ ( − ∆ − Ω − m Ω)+ { N ( − ∆ − Ω + m Ω) } θ ( − ∆ − Ω + m Ω)+ N (∆ + Ω − m Ω) θ (∆ + Ω − m Ω)+ N (∆ + Ω + m Ω) θ (∆ + Ω + m Ω) (cid:21) − πλ D ph (cid:18) (cid:19) (cid:20) { N ( − ∆ + Ω − m Ω) } θ ( − ∆ + Ω − m Ω)+ { N ( − ∆ + Ω + m Ω) } θ ( − ∆ + Ω + m Ω)+ N (∆ − Ω − m Ω) θ (∆ − Ω − m Ω)+ N (∆ − Ω + m Ω) θ (∆ − Ω + m Ω) (cid:21) (C8)Above θ is the Heaviside step function.0 F. D. M. Haldane, Phys. Rev. Lett. , 2015 (1988). M. Z. Hasan and C. L. Kane,Rev. Mod. Phys. , 3045 (2010). X.-L. Qi and S.-C. Zhang,Rev. Mod. Phys. , 1057 (2011). C. L. Kane and E. J. Mele,Phys. Rev. Lett. , 146802 (2005). B. A. Bernevig, T. L. Hughes, and S.-C. Zhang,Science , 1757 (2006). T. Senthil, arXiv:1405.4015 (unpublished). T. Oka and H. Aoki, Phys. Rev. B , 081406 (2009). J. Inoue and A. Tanaka,Phys. Rev. Lett. , 017401 (2010). T. Kitagawa, E. Berg, M. Rudner, and E. Demler,Phys. Rev. B , 235114 (2010). N. H. Lindner, G. Refael, and V. Galitski, Nature Physics , 490 (2011). T. Kitagawa, T. Oka, A. Brataas, L. Fu, and E. Demler,Phys. Rev. B , 235108 (2011). N. H. Lindner, D. L. Bergman, G. Refael, and V. Galitski,Phys. Rev. B , 235131 (2013). A. G´omez-Le´on and G. Platero,Phys. Rev. Lett. , 200403 (2013). Y. T. Katan and D. Podolsky,Phys. Rev. Lett. , 016802 (2013). P. M. Perez-Piskunow, G. Usaj, C. A. Balseiro, and L. E.F. F. Torres, Phys. Rev. B , 121401 (2014). J. H. Shirley, Phys. Rev. , B979 (1965). H. Sambe, Phys. Rev. A , 2203 (1973). M. S. Rudner, N. H. Lindner, E. Berg, and M. Levin,Phys. Rev. X , 031005 (2013). M. Lababidi, I. I. Satija, and E. Zhao,Phys. Rev. Lett. , 026805 (2014). A. Kundu, H. Fertig, and B. Seradjeh, arXiv:1406.1490(unpublished). M. Rechstman, J. Zeuner, Y. Plotnik, Y. Lumer, D. Podol-sky, F. Dreisow, S. Nolte, M. Segev, and A. Szameit, Na-ture (London) , 196 (2013). G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat,T. Uehlinger, D. Greif, and T. Esslinger, arXiv:1406.7874(unpublished). A. Lazarides, A. Das, and R. Moessner,Phys. Rev. Lett. , 150401 (2014). M. A. Sentef, M. Claassen, A. F. Kemper, B. Moritz,T. Oka, J. K. Freericks, and T. P. Devereaux,arXiv:1401.5103 (unpublished). N. Goldman and J. Dalibard, Phys. Rev. X , 031027(2014). Z. Gu, H. A. Fertig, D. P. Arovas, and A. Auerbach,Phys. Rev. Lett. , 216601 (2011). A. Kundu and B. Seradjeh,Phys. Rev. Lett. , 136402 (2013). L. E. F. F. Torres, P. M. Perez-Piskunow, C. A. Balseiro,and G. Usaj, arXiv:1409.2482 (unpublished). S. Diehl, E. Rico, M. A. Baranov, and P. Zoller, NaturePhysics , 971 (2011). J. C. Budich, P. Zoller, and S. Diehl, arXiv:1409.6341(unpublished). A. Uhlmann, Rep. Math. Phys. , 229 (1986). A. Rivas, O. Viyuela, and M. A. Martin-Delgado,Phys. Rev. B , 155141 (2013). O. Viyuela, A. Rivas, and M. A. Martin-Delgado,Phys. Rev. Lett. , 130401 (2014). Y. H. Wang, H. Steinberg, P. Jarillo-Herrero, andN. Gedik, Science , 453 (2013). Y. Onishi, Z. Ren, M. Novak, K. Segawa, Y. Ando, andK. Tanaka, arXiv:1403.2492 (unpublished). J. Karch, P. Olbrich, M. Schmalzbauer, C. Zoth, C. Brin-steiner, M. Fehrenbacher, U. Wurstbauer, M. M. Glazov,S. A. Tarasenko, E. L. Ivchenko, D. Weiss, J. Eroms,R. Yakimova, S. Lara-Avila, S. Kubatkin, and S. D.Ganichev, Phys. Rev. Lett. , 227402 (2010). J. Karch, C. Drexler, P. Olbrich, M. Fehrenbacher,M. Hirmer, M. M. Glazov, S. A. Tarasenko, E. L. Ivchenko,B. Birkner, J. Eroms, D. Weiss, R. Yakimova, S. Lara-Avila, S. Kubatkin, M. Ostler, T. Seyller, and S. D.Ganichev, Phys. Rev. Lett. , 276601 (2011). A. O. Caldeira and A. J. Leggett,Phys. Rev. Lett. , 211 (1981). B. M. Fregoso, Y. H. Wang, N. Gedik, and V. Galitski,Phys. Rev. B , 155129 (2013). S. Kohler, J. Lehmann, and P. H¨anggi, Phys. Rep. ,379 (2005). M. Langemeyer and M. Holthaus,Phys. Rev. E89