Floquet prethermalization in the resonantly driven Hubbard model
Abstract
We demonstrate the existence of longlived prethermalized states in the Mott insulating Hubbard model driven by periodic electric fields. These states, which also exist in the resonantly driven case with a large density of photoinduced doublons and holons, are characterized by a nonzero current and an effective temperature of the doublons and holons which depends sensitively on the driving condition. Focusing on the specific case of resonantly driven models whose effective timeindependent Hamiltonian in the highfrequency driving limit corresponds to noninteracting fermions, we show that the time evolution of the double occupation can be reproduced by the effective Hamiltonian, and that the prethermalization plateaus at finite driving frequency are controlled by the nexttoleading order correction in the highfrequency expansion of the effective Hamiltonian. We propose a numerical procedure to determine an effective Hubbard interaction that mimics the correlation effects induced by these higher order terms.
Introduction. The properties of materials can be tuned by chemical substitution, applied pressure, or static external fields. In the theoretical description, these modifications correspond to changes in the Hamiltonian parameters or the addition of extra terms describing the applied fields. In recent years, “Floquet engineering” has emerged as a versatile tool which enables new levels of control [1]. The idea is to apply periodic perturbations to a system which lead to modified parameters or new terms in the effective static Hamiltonian describing the “stroboscopic” evolution from one period to the next. The theoretically predicted phenomena range from hopping renormalizations [2, 3] and modified exchange couplings [4, 5, 6] to synthetic gauge fields [1, 7] and topological phase transitions [8, 9, 10]. Some of these effects have been demonstrated in cold atom systems [11, 12, 13, 14] and laserirradiated topological insulators [15].
Since the Hamiltonian of a periodically driven system satisfies , where is the period, the time evolution operator can be written as , where is a timeindependent static Hamiltonian and the socalled kick operator [1]. In the regime of large driving frequency one can derive by a highfrequency expansion, which is truncated at a given order in [1, 16, 17]. While the resulting effective model may contain new terms associated with interesting physical effects or novel phases, lowenergy properties of this model can only be realized if heating effects are small, and the driven system is stuck in a longlived quasisteady state exhibiting the desired properties. When the driving frequency is large enough compared to the characteristic energy scales of the system, the effective Hamiltonian evaluated with the high frequency expansion describes the system for exponentially long times [16, 18]. However, for nonintegrable systems and not too large driving frequency, the validity of these assumptions is not a priori clear. Such systems are expected to heat up in the presence of periodic driving, and it is an interesting question whether, and for how long, a quasisteady “Floquet prethermalized state” (FPS) different from the trivial infinite temperature state can be established.
Several recent theoretical works have demonstrated the existence of longlived FPSs in interacting models [16, 19, 20, 21, 22, 23], or questioned the general belief that heating to infinite temperature occurs in such systems. In this work we consider the Mott insulating singleband Hubbard model in timeperiodic electric fields. We will show that FPSs exist under various driving conditions, even when the driving frequency is smaller than the Hubbard interaction and comparable to the bandwidth. Focusing specifically on the case of resonant driving, where the local interaction is a multiple of the driving frequency and a violent heating might be naively expected, we demonstrate that the system can be trapped in longlived states with a suppressed number of double occupations and/or a nonzero current and kinetic energy. In the driving regime where the leading order term in describes free fermions, we can interpret the FPS as a state resulting from a quench to a weakly interacting effective Hubbard model.
Model. The Hamiltonian of the driven Hubbard model is with the creation operator for an electron of spin at site , the number operator, the nearestneighbor hopping, the onsite interaction, the electron charge (which is set to 1), and the electric field with polarization . To solve the model we use the nonequilibrium dynamical mean field theory (DMFT) [24] in combination with a strong coupling perturbative impurity solver (noncrossing approximation, NCA) [25, 26]. We consider an infinitedimensional hybercubic lattice with a Gaussian density of states and apply the electric field along the body diagonal, . In a gauge with pure vector potential , the electric field enters the calculation as a timedependent shift of the noninteracting dispersion, . The implementation and NCA treatment of this problem has been discussed in Refs. [27, 24]. We start at in the equilibrium state at inverse temperature and switch on the electric field as , where is the field amplitude and the driving frequency. Energy is measured in units of and time in units of .
Results. Figure 1 shows the time dependence of the current , the double occupation and the kinetic energy in a model with for indicated values of the field amplitude and driving frequency. (For the DMFT measurement of these quantities, see Ref. [24].) All these results correspond to resonant driving (, with integer). In the nonresonant case, the time evolution is slow, and we cannot reach the timescales needed to observe a saturation in a prethermalized state, or a heating to infinite temperature. In the following, we will thus focus on resonantly driven systems, where the doublonholon production is strong and a (quasi)steady FPS is rapidly reached.
Panels a) and e) show that for small driving frequency and large field amplitude, the system quickly approaches the infinite temperature state characterized by , and . Especially for the stronger field () the main mechanism for doublonholon production in this lowfrequency driving regime is fieldinduced tunneling, which creates an almost flat (infinite temperature) energy distribution of the photo carriers [28, 29]. On the other hand, for and , the system can be trapped in a longlived quasisteady state with a suppressed double occupation and a strongly oscillating current and kinetic energy (panels b), f), c) and g)). Note that the drift of the double occupation and kinetic energy, and hence the heating rate in the quasisteady state, depends in a nontrivial way on the field amplitude and driving frequency. There are also examples of FPSs with , but (panels d) and h)). In general, the time scale on which the steady state (either infinite temperature state or FPS) is reached becomes longer as we decrease the field strength.
The effective temperature of the doublons and holons in the trapped state can be estimated by calculating the timedependent spectral functions from the lesser and retarded local Green’s function as and . Averaging these spectral functions over one period yields and , which can be used to define the “nonequilibrium distribution function” . Panels i)l) of Fig. 1 show and obtained by averaging the Fourier transforms starting from . The almost flat distribution function in panel i) confirms the expectation from the quasistatic picture. Fitting the slopes of the distribution functions in panels j)l) with a Fermi function in the energy range of the Hubbard bands yields , for , , for , and , for . Hence, the effective temperature of the longlived trapped state can be lower than the temperature of the initial equilibrium state (). For driving frequencies somewhat above the resonant value we furthermore find negative effective temperatures, i.e. inverted populations of doublons and holons in the Hubbard bands.
We also note that satisfies . This property can be explained by analyzing the Green’s functions in terms of the kick operator and the effective Hamiltonian assuming that the system has reached a thermal state of the effective Hamiltonian at and that the characteristic energy scale of the effective Hamiltonian is smaller than . These assumptions imply that around behaves like the Fermi distribution function at , see Supplemental Material for details.
While the response of the Mott insulating Hubbard model to the electric field driving is complex and depends sensitively on the parameters of the driving field, our DMFT results clearly demonstrate the existence of longlived FPSs, even for resonant driving, moderate frequencies () and for large field amplitudes. An interesting question is how well the FPSs and the relaxation into these states are described by the timeindependent effective Hamiltonian.
To shed some light on this issue, we consider the case of resonant driving with even , where the effective Hamiltonian becomes simple. As shown by Bukov, Kolodrubetz and Polkovnikov [7], the leadingorder effective Hamiltonian has two terms, describing doublon/holon hopping, and doublonholon production/recombination. For even and suitably chosen driving amplitude, corresponds to free fermions. In the specific case of a hypercubic lattice with hopping , the leading term in the highfrequency expansion becomes , with , and (). Here, denotes the th order Bessel function, the operator describing the hopping of doublons and holons, and the operator for doublonholon production.^{1}^{1}1To be precise, this effective Hamiltonian is for . In our case of , the leading order effective Hamiltonian from the van Vleck highfrequency expansion is identical to the free Hamiltonian after a further unitary transformation, which can be absorbed into a redefinition of the kick operator. If is chosen such that the two amplitudes are equal, i.e. , and is even so that , then the driven system (in the highfrequency limit) is expected to behave like a noninteracting model with hopping amplitude .
We demonstrate this behavior in Fig. 2, where we compare the time evolution of the double occupation of the driven system with interaction to the double occupation in a Hubbard model after a quench from to and a simultaneous reduction in the hopping amplitude. We choose and , so that , and we quench the hopping amplitude from to . For (panel a)) we are in the high frequency driving regime, where the leading order effective Hamiltonian of Bukov et al. should be valid. Indeed, the time dependence of the double occupation shows the behavior expected for a quench from the undriven to and the double occupation increases to a value close to . As the interaction (and hence the driving frequency) is reduced (panels c) and e)), larger deviations between the driven system and the effective static description appear. While the quench to the noninteracting inevitably leads to a saturation of the double occupation at , the dynamics of the driven system shows a trapping of at a value below 0.25.
The observed deviations from the effective model description must be due to the higher order corrections in . We have calculated (see Supplemental Material), and it contains a large number of terms involving up to three different sites. One can see that induces correlations and modifications of the bandwidth, which represent the O() corrections to the leading order noninteracting Hamiltonian . For example, some of these terms describe the hopping of a pair of electrons from (to) the same site. This acts like a local interaction whose strength is determined by some average kinetic energy squared times a prefactor . Another term describes a three body interaction , which may also act as a local interaction whose strength is determined by the average of the occupancy on neighbouring sites. In addition there are correlated hopping terms which cannot be reduced to an effective Hubbard interaction, but which may change the effective bandwidth. It is thus an interesting question to what extent the effect of the additional terms in can be captured by a simple Hubbard Hamiltonian with modified interaction and bandwidth. In the following, we will demonstrate that to a large extent, acts as a local Hubbard interaction, and we will use this insight to interpret the FPSs with shown in Fig. 2.
First of all, we note that a significant change of the effective bandwidth would manifest itself in a change of the timescale on which the double occupation grows after the electric field quench. However, Fig. 2 shows that the quench correctly reproduces this growth rate not only in the highfrequency regime (), but also for and . This implies that the band widening effect of is not significant.
Instead, our numerical analysis shows that the effects of are, to a large extent, mimicked by a local Hubbard interaction . To determine , we drive the system with , choose corresponding to the noninteracting condition for , and vary the interaction of the driven system as . In the highfrequency driving regime, this model should behave like a static Hubbard model with interaction and a rescaled hopping parameter. At finite driving frequency, produces additional interaction effects, so that the effective static model is . The measured timeaveraged double occupation in the trapped state is plotted as a function of in the top panel of Fig. 3. The double occupation reaches a maximum very close to the noninteracting value at some negative value of . There is apparently a cancellation between the interaction coming from and the local interaction at this particular point. Conversely, the resonantly driven system with has an effective Hubbard interaction of strength , which originates from . Indeed scales with (see panel c)), which is consistent with the interpretation of an effective interaction coming from the next leading order.
In panel b) of Fig. 3, we plot the double occupation in the FPS as a function of . If the effective Hubbard interaction would perfectly capture the effect of higherorder terms in , and the absorbed energy were independent of and the kick operator, we would expect a collapse of the curves for different . The shifted data show a rather good agreement for , and , but small deviations remain. These deviations indicate that some of the correlations induced by cannot be captured by an effective Hubbard interaction, and they provide a rough estimate of these beyondHubbard effects on .
An interesting question is why the double occupation shows a parabolic maximum near as the interaction is varied near the resonance condition (Fig. 3a)). In fact, the expression for suggests that below the effective noninteracting point , the driven system should behave like an attractive Hubbard model, which usually yields . The observed suppression for occurs because in this regime, the driven system exhibits an inverted population. In panel a) of Fig. 4 we plot the effective inverse temperature of the model with as a function of . This figure shows that the effectively noninteracting driven system has an infinite temperature distribution, while the effectively attractive system has a negative effective temperature. The nonequilibrium distribution functions for , , are illustrated in panel b). As discussed in Ref. [3], the Hubbard model with interaction and negative effective temperature can be mapped onto a Hubbard model with interaction and positive temperature. This explains why the doublon occupation does not exceed 0.25 even when the effective model shows an attractive interaction. On the other hand, the resonantly driven model () with effective interaction coming from has a positive temperature distribution and an effective inverse temperature comparable to the initial equilibrium state (). For driving frequencies below the resonance (), the driven system can be effectively colder than the initial state.
We also note that the parabolic dependence of on and the scaling imply that the deviation of the double occupation in the FPS from 0.25 is proportional to . This at first sight surprising scaling is the result of an effective interaction proportional to and the fact that the system approaches an infinite temperature state in the highfrequency limit as , see Fig. 4e).
It is interesting to check how well the effective Hubbard model explains the observed values of the double occupation in the FPSs of the resonantly driven system (Fig. 2). An interaction quench to does not reproduce the plateau value very well (black dashed line). This is because the absorbed energy, and hence the effective temperature of the trapped or quenched state, depends sensitively on the details of the transient evolution, i.e. the kick operator. It is thus more meaningful to extract the effective inverse temperature of the driven system from a Fermi function fit of and to compare the double occupation of a Hubbard model with interaction and inverse temperature to the double occupation in the FPS. The corresponding results are indicated by the red dashed lines in the right hand panels of Fig. 2, and they are in rather good agreement with the timeaveraged double occupation. The remaining deviation to the FPS is comparable to the deviations evident in Fig. 3b), and may be attributed to “beyondHubbard” interaction effects.
We finally comment on the question whether the FPS observed here is a thermal or prethermal state in terms of . Due to the vicinity to the integrable noninteracting limit, the FPS might be expected to be a longlived prethermalized state [30], where the properties of nonlocal observables are different from those of the thermalized system described by . Though the direct simulation of and the analysis of nonlocal observables are beyond the scope of this study, we confirmed that the quench to the effective Hubbard model with shows a fast thermalization of local observables such as the double occupation, kinetic energy, and distribution function .
Summary. In this study, we have analyzed the properties of the resonantly driven Mott insulating Hubbard model. Contrary to naive expectations, and despite an efficient doublonholon production in the resonant regime, this nonintegrable system can be trapped in longlived Floquet prethermal states characterized by a suppressed double occupation and a nonzero current. While such trapping phenomena are found under various driving conditions, we have focused on the case , with even, where the leading order effective Hamiltonian reduces to a noninteracting fermion model. In this driving regime, the longlived trapped states can be understood as states resulting from a quench to a weakly interacting effective Hamiltonian. While these interactions originate from higher order terms in the highfrequency expansion, i.e., multisite correlated hopping terms, their effect on local observables can to a large extent be captured by an effective Hubbard repulsion . We have demonstrated a numerical procedure for evaluating and showed that it scales with , as expected for an interaction resulting from the nextleading order. At driving frequencies comparable to the bandwidth, there are however nonnegligible beyondHubbard interaction effects.
We have also demonstrated that driving below the resonance can lead to a cooling of the doublons and holons, and that for driving above the resonance, the driven Mott insulators can exhibit inverted doublon and holon populations. In the latter case, a sign change in the interaction terms of the effective Hamiltonian is required to correctly describe the properties of the Floquet prethermalized states by means of the effective static model with a positive temperature.
Acknowledgments The calculations have been performed on the Beo04 cluster at the University of Fribourg. This work was supported by the European Research Council through ERC Starting Grant 278023 and Consolidator Grant 724103, and by the Swiss National Science Foundation through NCCR Marvel.
References
 [1] Bukov M., D’Alessio L., and Polkovnikov A., Advances in Physics, 64, (2015) 139.
 [2] Eckardt A., Weiss C., and Holthaus M., Phys. Rev. Lett., 95 (2005) 260404.
 [3] Tsuji N., Oka T., Werner P. and Aoki H., Phys. Rev. Lett., 106 (2011) 236401.
 [4] Mentink J., Balzer K., and Eckstein M., Nat. Commun., 6, (2015) 6708.
 [5] Kitamura S. and Aoki H., Phys. Rev. B, 94, (2016) 174503.
 [6] Eckstein M., Mentink J., and Werner P., arXiv:1703.03269.
 [7] Bukov M., Kolodrubetz M., and Polkovnikov A., Phys. Rev. Lett., 116 (2016) 125301.
 [8] Oka T. and Aoki H., Phys. Rev. B, 79, (2009) 081406.
 [9] Lindner N., Refael G., and Galitski V., Nature Physics, 7, (2011) 490.
 [10] Kitagawa T., Oka T., Brataas A., Fu L., and Demler E., Phys. Rev. B, 84, (2011) 235108.
 [11] Struck J., et al., Phys. Rev. Lett., 108, (2012) 225304.
 [12] Aidelsburger M., Atala M., Lohse M., Barreiro J. T., Paredes B., and Bloch I., Phys. Rev. Lett., 111, (2013) 185301.
 [13] Jotzu G., Messer M., Desbuquois R., Lebrat M., Uehlinger T., Greif D., and Esslinger T., Nature, 515 (2014) 237.
 [14] Görg F., Messer M., Sandholzer K., Jotzu G., Desbuquois R., and Esslinger T., arXiv:1708.06751.
 [15] Wang Y. H., Steinberg H., JarilloHerrero P., and Gedik N., Science, 342 (2013) 453.
 [16] Mori T., Kuwahara T., and Saito K., Phys. Rev. Lett., 116, (2016) 120401.
 [17] Mikami T., Kitamura S., Yasuda K., Tsuji N., Oka T., and Aoki H., Phys. Rev. B, 93, (2016) 144307.
 [18] Kuwahara T., Mori T., and Saito K., Annals of Physics, 367, (2016) 96.
 [19] D’Alessio L. and Polkovnikov A., Annals of Physics, 333, (2013) 19.
 [20] Citro R., Dalla Torre E. G. , D’Alessio L. , Polkovnikov A., Babadi M., Oka T., and Demler E., Annals of Physics, 360, (2015) 694.
 [21] Bukov M., Gopalakrishnan S., Knap M., and Demler E., Phys. Rev. Lett., 115, (2015) 205301.
 [22] Canovi E., Kollar M., and Eckstein M., Phys. Rev. E, 93, (2016) 012130.
 [23] Weidinger S. and Knap M., Scientific Reports, 7, (2017) 45382.
 [24] Aoki H. et al., Rev. Mod. Phys., 86, (2014) 779.
 [25] Keiter H. and Kimball J. C., Intern. J. Magnetism, 1, (1971) 233.
 [26] Eckstein M. and Werner P., Phys. Rev. B, 82 (2010) 115115.
 [27] Eckstein M. and Werner P., Phys. Rev. B, 84, (2011) 035122.
 [28] Oka T., Phys. Rev. B, 86, (2012) 075148.
 [29] Eckstein M. and Werner P., Journal of Physics: Conf. Series, 427, (2013) 012005.
 [30] Moeckel M. and Kehrein S., Phys. Rev. Lett., 100 (2008) 175702.
Supplemental material: Floquet prethermalization in the resonantly driven Hubbard model
A.1. General idea of the SchriefferWolff transformation for periodically driven systems
In our study, we follow Ref. [7] to evaluate the effective Hamiltonian for the AC driven Hubbard model in the resonant or nearly resonant regime. In a rotating frame, the Hamiltonian becomes time dependent with the characteristic frequencies of the interaction () and the excitation frequency (). We can then define a common energy scale of these two and perform a highfrequency expansion if is much higher than the kinetic energy. Without the driving this procedure leads to the SchriefferWolff transformation for the  model or the Heisenberg model.
To be more precise, our Hamiltonian is the Hubbard model under a periodic excitation of the form
(0) 
The rotating frame with respect to , where , is chosen to obtain
(0) 
where
(0) 
Here , represents the hopping of doublons and holons and the creation of doublonholon pairs.
The second step is to regard this Hamiltonian as a timeperiodic Hamiltonian and to perform a highfrequency expansion. When the Hamiltonian is timeperiodic, from the Floquet theorem, the full time evolution can be decomposed as . Here is the timeindependent effective Hamiltonian and is the kick operator, which is time periodic with period . Usually, these operators can be approximately evaluated with highfrequency expansions. We note that the expansion is not unique and several methods exist [17]. One possible procedure is the van Vleck highfrequency expansion [7].
To apply the highfrequency expansion, we introduce a common frequency so that and and write
(0) 
From the van Vleck highfrequency expansion, the leading terms of become
(0)  
(0) 
where .
We also note that
in these expressions,
we can reduce some terms
under the assumption
that the number of doublons and holons is small.
A.2. Application to the resonantly driven Mott insulator
Now we apply this theory to our case of the infinitedimension hypercubic lattice under AC driving near the resonance. This is the case when . Then one can resonantly excite the system using multiphoton processes. Therefore one cannot neglect the terms involving doublons or holons.
The driving term explicitly reads
(0) 
where . For , . From this one can see that
(0) 
Here and and denotes the th order Bessel function. For , the effective Hamiltonian is identical to the case of after the unitary transformation .
The leading order becomes
(0) 
Hence, the leading term describes correlated hoppings, and one can see that at some special points the system becomes effectively free.
We now calculate the next order correction . One can identify three different groups of terms.

Group 1
(0) These terms vanish in our case of AC field driving because .

Group 2
(0) 
Group 3
(0) Here we can directly show that .
Let us write down the explicit form of the terms appearing in these expressions.
The coefficient of the terms can be evaluated by substituting Eq. ( ‣ Floquet prethermalization in the resonantly driven Hubbard model), but here we just want to show
what kind of processes are generated by .
Group 2 :
Group 2 :
Group 3 :
B. Spectral function and distribution function in the Floquet prethermal state
In the numerical simulation, we observed a periodic behavior of the timeaveraged distribution function in the prethermal states, . Here we give an argument why this periodicity appears and how it is connected to the effective temperature assuming that the prethermal state is an equilibrium state of the effective Hamiltonian.
In the following we assume the resonant condition, . We first note that our effective Hamiltonian is for the rotating frame, . Therefore, we also need to transform the expression of the Green’s function into the rotating frame.
First we consider the greater component,
(0) 
Here , and
(0) 
By introducing the timeperiodic operators and and ,