{'fraction_non_alphanumeric': 0.05711400359066427, 'fraction_numerical': 0.02053411131059246, 'mean_word_length': 4.122126436781609, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heating Rates in Periodically Driven Strongly Interacting Quantum Many-Body Systems

We study heating rates in strongly interacting quantum lattice systems in the thermodynamic limit. Using a numerical linked cluster expansion, we calculate the energy as a function of the driving time and find a robust exponential regime. The heating rates are shown to be in excellent agreement with Fermi's golden rule. We discuss the relationship between heating rates and, within the eigenstate thermalization hypothesis, the smooth function that characterizes the off-diagonal matrix elements of the drive operator in the eigenbasis of the static Hamiltonian. We show that such a function, in nonintegrable and (remarkably) integrable Hamiltonians, can be probed experimentally by studying heating rates as functions of the drive frequency.

Periodic perturbations are a ubiquitous tool to excite and probe quantum systems and study their response functions. Recent developments in theory and experiments have expanded the scope of periodic driving to generate effective magnetic fields [1][2][3][4], as well as to engineer topologically nontrivial band structures [5][6][7][8] and novel time-crystalline phases [9][10][11][12][13][14]. However, under periodic driving, generic many-body interacting systems are expected to heat up and (for a bounded spectrum, typical of lattice fermions and spins) equilibrate at long times to states that are effectively at infinite temperature [15,16].
Here, we report on the implementation of a numerical linked cluster expansion (NLCE) for driven systems. NL-CEs can be used to study arbitrary interaction strengths in arbitrary dimensions. They were originally introduced to study thermal equilibrium ensembles [43], where they outperform full exact diagonalization calculations [44]. NLCEs were recently implemented to study thermalization [45] and quantum dynamics under time-independent Hamiltonians in one [30,46] and two [47,48] dimensions, and combined with dynamical quantum typicality [49]. We use them to determine heating rates in strongly in-teracting one-dimensional (1D) lattices in the thermodynamic limit. The numerically obtained rates are shown to agree with Fermi's golden rule predictions. We argue that, in addition to helping quantify the stability of prethermal states, heating rates can be used to probe the structure of the off-diagonal matrix elements of the drive operator in the eigenstates of the static Hamiltonian.
The obvious conservation law broken by g(τ )K is energy conservation. For sufficiently small g in the thermodynamic limit, we expect prethermalization to occur (independently of the value of Ω), wherein the system quickly relaxes to the equilibrium state ofĤ 0 described by a (generalized) Gibbs ensemble [up to O(g) corrections]. The relaxation towards infinite temperature can be described by a slowly evolving (generalized) Gibbs ensemble ofĤ 0 , characterized by the instantaneous expectation values of the conserved quantities ofĤ 0 [30]. The dynamics of those quantities is described by autonomous equations, with drifts given by Fermi's golden rule [30].
where standard notation was used [52]. We drive the system with a square wave g(τ ) = g sgn[sin(Ωτ )], and set t = V = 1 (our unit of energy and frequency).Ĥ 0 is integrable for t = V = h = 0 (and mappable to the spin-1/2 XXZ Hamiltonian [52]), and nonintegrable for nonvanishing t , V , and h. We study integrable and nonintegrable (with t = V = 0.8 and h = 1.0) cases, and selectĤ I to have the same terms asĤ 0 [Eq. (2)] but with different nearest neighbor coupling parameters (t I = 0.5 and V I = 2.0). We implement a NLCE to calculate the energy per site e(τ ) = E(τ )/L at stroboscopic times in the thermodynamic limit (L → ∞). Within NLCEs, e(τ ) is expressed as a sum over the contributions of all connected clusters (c) that can be embedded on the lattice, e(τ ) = c M (c) × W e c (τ ), where M (c) is the number of "embeddings" (per site) of cluster c, and W e c (τ ) is the weight of e(τ ) in cluster c. W e c (τ ) is obtained recursively using the inclusion-exclusion principle: W e c (τ ) = E c (τ ) − c ⊂c W e c (τ ), where c denotes the connected subclusters of c and E c (τ ) = Tr[Ĥ c 0ρc (τ )] is the energy in cluster c [Ĥ c 0 is the static Hamiltonian, andρ c (τ ) is the density matrix at time τ , both in cluster c]. The series starts with the smallest cluster (a site) for which W c (τ ) = E c (τ ). For each cluster, E c (τ ) is calculated numerically using full exact diagonalization. We use maximally connected clusters (clusters with contiguous sites and all possible bonds) as they are optimal to study dynamics in chains in the presence of nearest and nextnearest neighbor interactions [45,46,54]. The order of the NLCE is set by the number of sites of the largest cluster considered. For nonintegrableĤ 0 , we compute 17 orders of the NLCE (after exploiting all symmetries, the dimension of largest sector of the Hamiltonian is 32 896). WhenĤ 0 is integrable, due to particle number conservation, we are able to compute 18 orders of the NLCE (the dimension of the largest sector in this case is 21 942).
In the main panels of Fig. 1, we show NLCE results for |e(τ )| vs τ for (a) the nonintegrable and (b) the integrable static Hamiltonians, for three strengths g = 0.05, 0.2, and 0.8 of the drive, for an initial thermal equilibrium state ofĤ I at an inverse temperature β I = (30) −1 . The exponential fits, which exclude the short-time transient dynamics and long times at which the NLCE does not converge, make apparent that the approach of e(τ ) to the infinite-temperature energy (E ∞ /L = 0) is exponential. The rates obtained from such fits are plotted in the insets of Fig. 1 vs g, for the two highest orders of the NLCE. They agree with each other, indicating that the fits are robust. The rates are ∝ g 2 and are in excellent agreement with Fermi's golden rule [Eq. (1)], evaluated numerically using full exact diagonalization in chains with periodic boundary conditions [53].
To illustrate this, in the main panel of Fig. 2 we plot |e(τ )| (normalized by its initial value |e(0)|) for various initial inverse temperatures β I ∈ [0.01, 0.5]. The normalized energies e(τ )/e(0) for β I = 0.033 and 0.01 exhibit a nearly identical exponential decay (within the times at which the NLCE has converged) implying that Γ is independent of β I [hence, of β(τ )] when β I 0.03. For β I 0.2, one can still use exponentials to fit e(τ ), but the rates obtained depend on β I . In the inset in Fig. 2, we report the rates obtained from such fits vs β I using two orders of the NLCE and for two values of g. The rates from the two orders of the NLCE agree with each other and agree well with Fermi's golden rule predictions. (A worse agreement is seen for g = 0.5 than for g = 0.2 due to the effect of higher order corrections.) The increase in the rate seen in the inset in Fig. 2 with decreasing β I is the one expected to occur as a function of driving time for initial states that are not in the regime β I 1. Next, we focus on the dependence of the heating rates on Ω. In nonintegrable systems, the eigenstate thermalization hypothesis [50,51,55,56] allows one to compute Γ m (τ ). After resolving all symmetries of the static Hamiltonian, the eigenstate thermalization hypothesis ansatz for the matrix elements K is the density of states ofĤ 0 in sector s at energy E, and R i,f is a random variable with zero mean and unit variance. K (s) (E) and f (s) K (E, ω) are smooth functions of their arguments.
Using Eqs. (1) and (4), changing sums over eigenstates by integrals over energy, replacingρ DE (τ ) byρ GE (τ ) and assuming high temperature [β(τ ) 1], one obtains the following expression for the heating rate [53] where E In Fig. 3(a), we compare heating rates (for the nonintegrable case and normalized by g 2 ) obtained from dynamics evaluated with NLCE (see inset) and the ones predicted by Eq. (5) [53]. NLCE results are not reported for small and large values of Ω because the time interval in which the NLCE converges is not sufficiently long to produce robust exponential fits. The normalized rates for g = 0.2 and g = 0.3 are nearly identical to one another, and are well described by Eq. (5). For high Ω, we find that the evaluation of Eq. (5) results in heating rates that can be well described by an exponential in Ω. This is consistent with rigorous bounds [10,31,36].
WhenĤ 0 is integrable (the spin-1/2 XXZ limit), the prethermal states are described by a generalized Gibbs ensemble (GGE)ρ GGE (τ ) [62][63][64]. Whenρ I is a thermal state with β I 1 (or in general after long driving times),ρ GGE (τ ) ρ GE (τ ) with β(τ ) 1 [65]. In this regime, Eq. (5) gives the heating rates for the integrable static Hamiltonian provided that there is a well defined |f Fig. 3(b), we show the equivalent of Fig. 3(a) but for the integrable case. Despite the differences between the dependence of the heating rates on Ω in the nonintegrable and integrable cases, the heating rates in the latter are described by Eq. (5) and, for high Ω, they are well described by an exponential in Ω.
The previous results show that heating rates can be used to probe the function f where Z(s) is the Hilbert space dimension of sector s. Thus, the rate for Fourier mode m = 1, which Fig. 4 shows to be in excellent agreement with the heating rates obtained from the NLCE dynamics for a wide range of values of Ω, gives the average |f K (E ∞ , Ω)| 2 over all sectors of the Hamiltonian in the thermodynamic limit [53].
In summary, we studied heating in strongly interacting driven lattice systems and showed that, at sufficiently high effective temperatures ([β(τ )] −1 2), it can be well characterized by rates no matter whether the system is nonintegrable or integrable. We also showed that the rates agree with Fermi's golden rule predictions for both nonintegrable or integrable cases. We then argued that heating rates can be used to probe the structure of off-diagonal matrix elements of the operator used to drive the system, in the eigenstates of the static Hamiltonian. Our results suggest that there is a well defined |f (s) K (E, Ω)| 2 in integrable interacting systems. This has been confirmed in a recent full exact diagonalization study of the spin-1/2 XXZ chain [66], and needs to be further explored to place it on equal footing with what is known for quantum chaotic systems [51,[67][68][69][70][71][72].
This work was supported by the National Science Foundation under Grant No. PHY-1707482. We are grateful to W. De Roeck and S. Gopalakrishnan for motivating discussions. The computations were carried out at the Institute for CyberScience at Penn State.
Supplemental Material: Heating rates in periodically driven strongly interacting quantum many-body systems  (1) in the main text is evaluated using full exact diagonalization of chains with L sites and periodic boundary conditions. Defining a small energy window ∆E, Eq. (1) is modified to the following expression (which is amenable to numerical evaluation) With this coarse graining procedure, Γ ∆E,m (τ ) for Fourier mode m is calculated as where E(τ ) is also evaluated using full exact diagonalization, and E ∞ = 0 for our model. Γ ∆E (τ ) = m>0 Γ ∆E,m (τ ) is the relaxation rate of E(τ ). In Fig. S1, we show Γ ∆E (τ ) vs τ for three values of ∆E when g = 0.05, 0.2, and 0.8, for the nonintegrable [L = 18, Fig. S1(a)] and the integrable [L = 20, Fig. S1(b)] static Hamiltonians (the period of the drive is T = 1.0). The initial thermal state has β I = (30) −1 . It is apparent in Fig. S1 that Γ ∆E (τ ) is nearly constant, with a slight drift at long times (apparent for g = 0.8), and that it is independent of the value of ∆E. We identify a range of ∆E and τ where Γ ∆E (τ ) is (nearly) constant [and where the dynamics of e(τ ) is exponential and robust against finite-size effects, see Sec. S3] and compute the average of Γ ∆E (τ ) in this range.
In the main text, all the rates reported in Figs From the eigenstate thermalization hypothesis (ETH) it follows that the results from the diagonal ensemble and the Gibbs ensemble agree [50,51,55,56], so one can replace P 0 is the partition function, and the inverse temperature β(τ ) is set by the energy E(τ ).
In Fig. S2, we also show results of the numerical evaluation of Eq. (6) in the main text using full exact diagonalization of chains with periodic boundary conditions and L sites. The coarse grained Eq. (6), using Eq. (S7), has the form It is apparent in Fig. S2, both for the nonintegrable and the integrable static Hamiltonians, that the results for Γ ∞ m=1 calculated using Eq. (S9) do not agree with the ones for Γ L m=1 using Eq. (S8). This is because of strong finite-size effects in Γ ∞ m=1 . We note that the disagreement increases as Ω increases. The fact that finite-size effects in Γ ∞ m=1 increase with increasing Ω is also apparent in the increasing discrepancy with increasing Ω between the results for the two chain sizes shown in Fig. S2. This is in contrast to the results for Γ L (similar to Γ L m=1 at large Ω) evaluated from Eq. (S8) and reported in Fig. 3 in the main text for two systems sizes. The strong finite-size effects in Γ ∞ m=1 are not surprising as the assumptions made to derive Eq. (6) are not valid for the small system sizes studied in this work.

S3. CONVERGENCE OF NLCE AND EXACT DIAGONALIZATION
In Fig. S3(a), we plot the absolute value of the energy per site |e(τ )| of the nonintegrable Hamiltonian for g = 0.5, period T = 1, and β I = 30 −1 . The results shown for |e(τ )| were obtained within the last four orders l of the numerical linked cluster expansion (NLCE), and for the four largest chain sizes L (with periodic boundary conditions) studied using exact diagonalization (ED). At short times (τ 10), both NLCE and ED give nearly identical results (all lines are indistinguishable in the scale of the figure). The small finite-size effects of the ED calculations for τ 10 are the reason we evaluate Eq. (1) of the main text with ED in the range 3 ≤ τ ≤ 10 (see Fig. S1). All the curves in this range are well described by an exponential, so we fit an exponential to the results from the largest chain calculated with ED (L = 18) in this range of τ [shown in Fig. S3(a) as a black line]. For τ 10, ED calculations deviate from the exponential, and with increasing L the curves monotonically approach the exponential fit. On the other hand, NLCE results (l = 16 and 17) remain exponential up to τ ∼ 20 [apparent in the inset in Fig. S3(a)]. Increasing the order l of the NLCE significantly improves the convergence towards the exponential at longer times. We remark here that, purely from the ED calculations, it is difficult to identify the exponential regime (in order to accurately predict the rates) as the results from ED smoothly drift away from an exponential and the discrepancies between L = 17 and 18 are small at most times. We use the NLCE results as reference in order to identify the time interval in which the ED results exhibit the "correct exponential". That time interval is then used to compute robust Fermi's golden rule predictions.
Next, we quantify the convergence errors of NLCE and the finite-size errors of ED calculations. At the times at which the highest order (l = 17) NLCE results are very close to the exponential fits to the ED (and NLCE) data, the NLCE results serve best as reference to quantify the errors at lower orders of the NLCE and finite-size errors of ED. For e(τ ) evaluated with the l th NLCE order [e l (τ )] and with the L-site periodic chain ED [e L (τ )], we define the convergence errors for l < 17 and L, respectively, as the relative differences from e l=17 (τ ) given by δ l NLCE (τ ) = |e l (τ ) − e l=17 (τ )| |e l=17 (τ )| and δ L ED (τ ) = |e L (τ ) − e l=17 (τ )| |e l=17 (τ )| .
(S10) Fig. S3(b) reports δ l NLCE (τ ) and δ L ED (τ ) vs l and L, respectively, at τ = 10 and τ = 12, for the same e(τ ) as in Fig. S3(a). The plots make apparent that the errors decrease with increasing l and L, that the errors at τ = 12 are greater than the corresponding ones at τ = 10, and suggest that the NLCE convergence errors decrease faster with increasing l than the ED finite-size errors with increasing L, a known fact in equilibrium calculations [44].
At τ = 10, the estimated of error for the ED calculation [δ L ED (10)] for L = 18 in Fig. S3(b) is less than 0.2%. For all calculations with ED in this paper, τ = 10 is the largest time considered for the nonintegrableĤ 0 .
As argued in the main text, e(τ ) in the thermodynamic limit is essentially a single exponential in τ for β I 30 −1 . Hence, for τ 20, we can estimate errors via the deviation of the NLCE and ED results from an exponential fit to the shorter-time results. To avoid any bias towards the NLCE results, here we compute the exponential fit from the ED results of e L=18 (τ ) for 3 ≤ τ ≤ 10. The normalized deviation from the exponential fit for the l th NLCE order [∆ l (τ )] and the L-site periodic chain solved with ED [∆ L (τ )], are given by ∆ l (τ ) = |e l (τ ) − fit(τ )| |fit(τ )| and ∆ L (τ ) = |e L (τ ) − fit(τ )| |fit(τ )| , where fit(τ ) is the value of the exponential fit at τ . In Fig. S3(c), ∆ l (τ ) and ∆ L (τ ) are plotted vs l and L, respectively, for τ = 22 and τ = 25. It is apparent that the NLCE results converge faster towards the exponential [at τ = 22, ∆ l=17 (τ ) is less than 0.5%]. The errors in the ED calculations are an order of magnitude higher, and vanish more slowly with increasing L [as in Fig. S3(b)]. The excellent convergence of the NLCE results (for τ ≤ 20 the nonintegrablê H 0 ) was essential for the accuracy of the relaxation rates computed via fits to the NLCE data that were reported in the main text.