Statistical Floquet prethermalization of the Bose-Hubbard model

The manipulation of many-body systems often involves time-dependent forces that cause unwanted heating. One strategy to suppress heating is to use time-periodic (Floquet) forces at large driving frequencies. For quantum spin systems with bounded spectra, it was shown rigorously that the heating rate is exponentially small in the driving frequency. Recently, the exponential suppression of heating has also been observed in an experiment with ultracold atoms, realizing a periodically driven Bose-Hubbard model. This model has an unbounded spectrum and, hence, is beyond the reach of previous theoretical approaches. Here, we study this model with two semiclassical approaches valid, respectively, at large and weak interaction strengths. In both limits, we compute the heating rates by studying the statistical probability to encounter a many-body resonance, and obtain a quantitative agreement with the exact diagonalization of the quantum model. Our approach demonstrates the relevance of statistical arguments to Floquet perthermalization of interacting many-body quantum systems.

The study of periodically driven systems has a long history, tracing back to the work of Floquet on classical systems governed by linear equations of motion [1]. Floquet showed that these equations can be solved using a time-independent unitary matrix, U F , which captures the evolution over one period of the drive, τ . Remarkably, because the time evolution of quantum systems is determined by a linear equation (namely, the Schrödinger equation), Floquet theory can be used to study any quantum system, even in the presence of interactions. The practical applicability of Floquet theory is hindered by the fact that finding U F , and diagonalizing it, is generically very difficult. This difficulty is especially acute for many-body quantum systems, where the size of U F grows exponentially with the number of degrees of freedom. Nevertheless, at large driving frequencies, U F can be derived using a controlled analytical approximation, the Magnus expansion [2]. The first term of this expansion is U F ≈ e −iHavτ , where H av = τ −1 τ 0 H(t)dt is the time-averaged Hamiltonian. The other terms are integrals of commutation relations of the Hamiltonian at different times [3].
Using the Magnus expansion, Refs. [4][5][6][7][8] were able to obtain rigorous constraints on the time evolution of periodically driven quantum many-body systems. These rigorous theorems apply to quantum spin systems that satisfy a local norm bound: their Hamiltonians consist of sums of local operators whose matrix elements are smaller than a given energy scale J. For these systems, the heating rate Φ was shown to be exponential suppressed at large driving frequencies Ω = 2π/τ , according to where is the Plank's constant, A and B are unitless constant. This exponential suppression was observed in several numerical studies [9][10][11][12] and in an experiment with dipolar spin chains [13]. The rigorous bound of Eq. (1) can be understood using a perturbative argument [4]: Due to the local norm bound, a single application of the driving field can change the energy of the system by J, at most. On the other hand, the absorption of a quantum of energy from the pump injects energy Ω. Hence, the absorption of energy from the pump requires the product of n = Ω/J operators and is governed by the nth order perturbation theory. Refs. [5][6][7][8] used the Magnus expansion to extend this argument and demonstrate that Eq. (1) is a rigorous bound, valid to all orders. Interestingly, in the limit of → 0, this bound applies to classical systems with a bounded spectrum [14,15].
Many physical systems escape the regime of validity of the aforementioned rigorous bounds. For example, massive particles with momentum p have a kinetic energy p 2 /2m that is unbounded from above. Ref. [16] demonstrated that systems of interacting particles can, nevertheless, show an exponential suppression of heating. They considered a canonical model of coupled kicked rotors [17][18][19][20] and showed that, for appropriate initial conditions, the system shows an exponentially long-lived prethermal plateau with vanishing energy absorption. This effect was explained in Ref. [21] using the following statistical argument: At large driving frequencies, the heating rate is small and the time-averaged energy of the system is (quasi) conserved. If the system is ergodic, the state of the system can be approximated by the Boltzmann distribution function, where Z is the partition function, k B is the Boltzmann constant, and the temperature T is determined by the initial energy of the system, measured with respect to the time-averaged Hamiltonian H av . If other quantities, such as the total momentum or the total number of particles are conserved, the appropriate Lagrange multipliers need to be taken into account. The resulting distribution can then be used to estimate the heating rate by computing the probability to incur into a many-body resonance [22]. Under physical assumptions, this probability is exponentially small, leading to a statistical Floquet prethermalization [21].
Having introduced the concepts of rigorous and statistical Floquet prethermalization, we now move to the focus of this article, namely the periodically driven Bose-Hubbard model, described by with J(t) = J 0 +δJ cos(Ωt). Here, b i and b † i are canonical bosonic operators, n i = b † i b i is the number of particles on site i and i, j are nearest neighbors. The U term describes onsite repulsion and the J term hopping. Importantly, the U term is unbounded from above, making the rigorous bounds of Ref. [4][5][6][7][8] unapplicable. The Hamiltonian of Eq. (3) conserves the total number of particles in the system, N = i n i andn denotes the average number of particles per site.
Floquet prethermalization in the Bose-Hubbard model was studied theoretically in Ref.
[23] using a selfconsistent quadratic approximation. This work employed the concept of many-body parametric resonance [24] to predict the existence of a frequency threshold above which the system does not absorb energy. However, in practice, terms that are neglected in the quadratic approximation lead to finite heating rates at all frequencies.
Ref. [4] predicted that at large driving frequency, the heating rate should be rigorously bounded by a stretched exponential [25]. In the limit of a large number of particles per site (n 1), the model can be mapped to a system of classical rotors, where the heating rate is exponential suppressed [21].
Recently, the heating rate of the Bose-Hubbard model with one particle per site (n = 1) was studied by Ref.
[26], using three methods: (i) the numerical calculation of the linear response of the model; (ii) the experimental measurement of single-site excitations (doublons or holes); (iii) the experimental measurement of the system's temperature. The experiments were performed using ultracold atoms in one and two-dimensional optical lattices. The time-periodic drive was obtained by modulating the intensity of the laser fields that generate the lattice [27]. The findings of Ref. [26] demonstrate that the heating rate is exponentially suppressed as a function of Ω in all dimensions. As explained, this observation cannot be accounted by the available theoretical methods.
In this article, we present two semiclassical approximations that capture the exponential suppression of the heating in two opposite limits. The first limit is strong interactions (U J), where we link the heating suppression to the low probability of finding many particles on a single site. The second limit is weak interactions (U J), where we can perform a controlled expansion of the heating rate in orders of U . For both cases, we use a statistical approach to compute the heating rate to lowest order in the strength of the periodic drive (∼ δJ 2 ) and compare it with the exact numerical diagonalization of the model.
Strong interactions (U J) -In the regime of large interactions, U J, we can describe the system in terms of semiclassical particles hopping on a lattice. The periodic drive moves one particle from one site to a neighboring one. This process changes the value of the on-site interaction by where the upper (or lower) sign refers to a particle hopping from site j to site i (or vice versa). Following Ref.
[21], we need to identify the many-body resonances of the model. Here, a resonance occurs when Eq. (4) equals to an integer multiple of the frequency of the drive (in units of Schrödinger's equation constant ), or ∆E = m Ω, where m is an integer. For high-frequency drives, the heating rate is dominated by the lowestorder available resonance, which corresponds to m = ±1. Without loss of generality, we assume that n i > n j , such that when a particles moves from j to i (or vice versa) the interaction energy increase (decreases). The resonance condition ∆E = ± Ω becomes ±(n i − n j ) + 1 = ±n Ω , or where we defined n Ω = Ω/U . Here, the upper (or lower) sign refers to the absorption (or emission) of energy. Note that this condition can be matched only if n Ω is integer. If the maximal occupation of each site is limited to n i ≤ 2, such as in the case of spin-1/2 fermions, the resonant condition can be satisfied only for n Ω = 1 [28]. In contrast, for bosons n i is unbounded and energy can be resonantly absorbed at arbitrarily high frequencies.
Because the probability to find sites with large n i is exponentially small, so is the probability to satisfy the resonance condition, leading to suppressed heating rates. The goal of this article is to put this intuitive argument on solid mathematical ground. The probability to satisfy Eq. (5) is determined by P i,j (n i , n j ), the joint distribution function to find n i and n j particles in sites i and j, according to P ± (Ω) = n P i,j (n, n − n Ω ± 1) .
This expression needs to be multiplied by a factor of 2 to take into account the case of n i < n j . In a d dimensional square lattice, we need to further multiply the result by the coordination number d [29]. Evaluating the distribution function P i,j (n i , n j ) in a (pre)thermal state described by Eq. (2) is a formidable task in many-body quantum physics. In what follows, we focus on the regime of large temperatures T J, where we can neglect quantum fluctuations and describe the prethermal state by with Here, in addition to the quasi-conservation of the energy in the prethermal state, we took into consideration the conservation of the total number of particles, through the chemical potential µ. The values of Z 0 and µ are determined by the constraints n Z i (n) = 1 and n nZ i (n) =n. These constraints, along with the numerical solution of Eqs. (6)-(8) enable us to compute the semiclassical heating rate of the Bose-Hubbard model, Φ. The total heating rate is given by the probability to incur into a resonance (P + −P − ), times the heating rate of an individual resonance. According to the linear response theory, one obtains where the delta function δ( Ω − ∆E) imposes the relevant resonance condition. To regularize this function, one needs to take into account the effects of small, but finite, J/U : the hopping term in Eq. (3) transforms the single particle states into "conduction bands" of width Λ = 4dJ. To model this effect, we substitute the delta function in Eq. (9) by a square function of width 2Λ, namely δ( ω) = [Θ( ω > −Λ)−Θ( ω > Λ)]/(2Λ), where Θ is the Heaviside function. In Fig. 1, we plot the resulting heating rates in d = 1, obtained from the numerical solution of our semiclassical approach, Eq. (6)-(9), for different values of the temperature [30]. We find that the heating rate is exponentially suppressed for all temperatures and, at large temperatures, inversely proportional to the temperature. To gain physical insight into this result, we now develop an analytical high-temperature expansion. In the limit of T → ∞, the distribution function is solely determined by the conservation laws and with Z 0 = 1 − z and z =n/(1 +n) [31]. By combining Eqs. (6) and (10), we obtain Note that the two sums have different lower limits because P + can occur only if n j ≥ 1, while P − requires only n j ≥ 0. Because P + = P − the net energy absorption is zero, Φ = 0. This result is not surprising: infinite temperature ensembles do not absorb energy! We can use this result as the starting point of a perturbative analysis. By approximating Eq.
leading to (see symbolic script in appendix B) In particular, atn = 1 (z = 1/2), we obtain Eq. (15) shows that, in the regime of U J and at very high temperatures, the heating rate of the Bose-Hubbard model is an exponential function of the ratio between the driving frequency and the onsite interaction. At intermediate temperatures, the heating rate is additionally suppressed by the fact U n 2 /(2k B T ) in Eq. (8), leading to a faster-than-exponential decay of Φ(Ω), see Fig. (1). Hence, Eq. (15) can be considered as an upper bound of the heating rate at all temperatures.
We now compare the results of our semiclassical approximation with the exact diagonalization of the Bose-Hubbard model. At finite temperatures, linear response gives [26] Here |ψ n and E n are, respectively, the eigenstates and eigenvalues of the average Hamiltonian H av atn = 1 and is the time-dependent perturbation. We evaluate this quantity numerically for N = 9 particles on a one dimensional lattice with L = 9 sites (n = N/L = 1) and open boundary conditions [33]. To mitigate the effects due to the finite dimension of the lattice, we have regularized the delta function of Eq.  using the above-mentioned square function with Λ = 2J. Because the maximal number of particles per site is always smaller or equal to the total number of particles N , we need to restrict ourselves to frequencies Ω, such that n Ω < N , or Ω < N U [34]. As shown in the Fig. 1, for all temperatures T > U the results of our numerical calculations are well approximated by the semiclassical description.
Weak interactions (U J) -We now turn to the other extreme limit, where the interactions are small in comparison to the kinetic energy and can be treated perturbatively. The periodically driven Bose-Hubbard model of Eq. (3) can be written as the sum of a timeindependent part H = H 0 + H int and a periodic drive, δJ cos(Ωt)V , with x e ikx b x . Using this notation, the heating rate of Eq. (9) takes the form where the square brackets denote a commutator and ... T is the expectation value with respect to a thermal state of the time-independent Hamiltonian H 0 + H int at temperature T . This expression can be computed numerically using path integrals techniques, either as the analytic continuation of an imaginary-time correlator, or as a real-time (Keldysh) response function.
In what follows, we present a semiclassical approach, aimed at computing Φ for U J. As we will show below, our approach captures the correct scaling laws of Φ and highlights its exponential suppression at large frequencies. We treat the eigenstates of H 0 as classical particles (quasi-particles), generated by the interaction term H int . The probability to observe a process involving the nth order of H int is given by Here, the factor n! derives from the nth order Taylor expansion of the exponent used in the perturbation theory. At zero temperature, this process creates up to n qp = n/2 + 1 quasiparticles. This relation is justified by the diagrams shown in Fig. 2, which demonstrate that the leading order contribution to the creation of n qp quasiparticles involves n = 2n qp −2 vertexes. From a semiclassical perspective, this relation indicates that the second order perturbation creates two quasiparticles (n qp = 2 for n = 2), and that the number of quasiparticles increases by one for every two additional orders of perturbation. We now use the statistical approach of Eq. (9) to compute the heating rate. A many-body resonance condition is satisfied when the total energy is conserved, namely if Ω = nqp j=1 ε kj < 4Jn qp . Hence, the lowest order resonance is obtained for n * qp = Ω/4J , where ... is the ceil function. The heating rate is, then, given by In Fig. 3(a) we compute Φ T =0 (Ω) as a function of U/J FIG. 2. Representative diagrams for the leading contributions in the second (n = 2) and fourth (n = 4) orders of the perturbative expansion. The annihilation (creation) operators are denoted by outgoing (incoming) arrows, from left to right. The maximal number of simultaneous quasiparticles is nqp = 2 and nqp = 3 for the second and fourth orders, respectively. All other diagrams of the same order (i.e. with the same number of vertexes) create less quasiparticles. The generalization to higher order diagrams is straightforward and shows that the nth-order perturbation can create up to nqp = n/2 + 1 quasiparticles.
using the exact diagonalization of a finite-size system (L = N = 9) and show that Eq. (20) captures the correct scaling behavior. As one increases the driving frequency Ω, the heating rate is dominated by higher orders of perturbation theory in U/J. Hence, at a fixed U/J < 1, the heating rate decreases exponentially with Ω. To see this effect, we consider a smooth version of Eq. (20) by approximating x ≈ x + 1/2 and substituting n! → Γ(n), leading to This expression is found to be in quantitative agreement with the exact diagonalization calculations, see Fig. 3(b). We now study the temperature dependence of the heating rate by considering the statistical properties of the aforementioned semiclassical quasiparticles. For simplicity, we approximate the band structure ε(k) as two plateaus, one at ε k=0 = −2J and one at ε k=π = 2J. In this simplified model, the creation of a quasiparticle involves an energy jump of ∆ε = 4J. This event is possible only if the k = 0 state is full and the k = π state is empty. The probability to excite n qp quasiparticles simultaneously is, then, Bose-Einstein distribution function and the chemical potential µ is determined by the condition f q=0 + f q=π = 1. The resulting heating rate is where Φ T =0 is given in Eq. The numerical results are compared with Eq. 21, which is the analytical continuation of the heating rate for all frequencies (dashed lines). Each range of energies is corresponded to a specific resonance condition, where increasing the external drive requires an additional order in the perturbative expansion of the heating rate. From the linear slope it is evident that for each such range, Φ scales like U 2 . In particular, the blue circles corresponded to the regime of 0 < Ω/J < 8, and further increasing the external drive in an additional 4 Ω/J, yields the power-law dependent for the other regions. (b) The scaling of the heating rate, Eq. (21) as a function of the external drive. This expression (dashed lines) gives a good estimation for the exponential suppression as obtained numerically from ED.
temperatures, when k B T /J approaches J/U , the subleading orders of our perturbative approach become nonnegligible and the analytical expression deviates from the exact numerical results. Note that as the temperature increases, the thermal weight in the square brackets of Eq. (22)  bound for the exponential suppression, which persists at all temperatures.
Conclusion -To summarize, we discussed the differences between rigorous [4-8] and statistical [21] Floquet prethermalization. The former approach relies on the boundedness of quantum operators and applies to spin models only. See also Ref. [35], where it was shown that the rigorous approach applied to systems of interacting particles with an unbounded spectrum does not lead to exponential bounds on diffusion rates. The latter approach relies on the statistical description of the prethermal state and applies to a wider range of models, including interacting particles in a lattice and in the continuum [36]. A key difference between these two approaches is that, while the rigorous approach is independent on the initial state, the statistical approach depends on the initial state, through its (quasi)conserved quantities, such as energy and particles' number.
In this article, we applied the statistical argument to the periodically driven Bose-Hubbard model, which was recently realized experimentally [26]. We developed two semiclassical descriptions of Floquet prethermal states, valid in two extreme regimes. The first limit corresponds to strong interactions and large temperatures (U > k B T J), where the suppressed heating rate is the outcome of the low probability to find many particles on a single site. The second limit corresponds to low temperatures and weak interactions (k B T < U J) and is relevant to the experiment of Ref. [26]. Here, the exponential suppression results from the low probability to create simultaneously many quasiparticles in momentum space. In both limits, we described the system semiclassically and applied statistical arguments to derive an analytical expressions for the heating rate Φ as a function of the driving frequency Ω and of the temperature T . These expressions are found to match the results of the exact diagonalization of the model, without any fitting parameter. Importantly, we demonstrated that in both regimes, the exponential suppression of the heating persists at all temperatures.
In this aspect, the Bose-Hubbard model differs from the coupled rotors model of Refs. [16][17][18][19][20][21], where the exponential suppression of heating disappears at large temperatures, eventually leading to a runaway from the prethermal regime. This fundamental difference stems from the nature of the conserved quantities of the two models: In the rotor model, the conserved quantity, namely the momentum of the rotors p i , is a continuous variable and can acquire both positive and negative values. At large temperatures, the fluctuations of p i diverge making the exponential suppression of heating ineffective. In contrast, in the Bose-Hubbard model, the conserved quantity, namely the particles' number n i , is non-negative. If the expectation value of n i is kept fixed, the fluctuations of this quantity remain finite and the heating rate is suppressed at all temperatures. The prediction of the two models coincide when the average number of particles per site is taken to infinity (n → ∞).
Our semiclassical approach disregards effects associated with quantum coherence. In the case of a single kicked rotor, quantum coherence strongly suppresses heating through the dynamical localization in energy space [37,38]. Accordingly, it was recently shown that dynamical localization can lead to ergodicity breaking in many-body kicked models, such as coupled rotors [39] and the Bose-Hubbard model [40]. However, as conjectured in Ref.
[41], dynamical localization is probably restricted to kicked models and, hence, is not relevant to the present study, where we considered a sinusoidal time dependence.
We  nn=np . a r r a y ( r a n g e ( 1 , 4 5 ) ) 123 p l t . s e m i l o g y ( nn , 1 / 3 * np . power ( 1 / 2 , nn/U) /2/gamma , ' k . : ' , l a b e l="Eq . ( 1 3 ) " ) 124 nn=np . a r r a y ( r a n g e ( 1 ,Nmax+1) ) 125 f o r i i n r a n g e ( sh  Fig. 1 for N particles on L sites. No fitting parameter is used. The semiclassical approximation matches the exact results for frequencies Ω < N/U and temperatures kBT > U .