Heating Rates in Periodically Driven Strongly Interacting Quantum Many-Body Systems

TL;DR

This paper investigates heating rates in strongly interacting quantum lattice systems under periodic driving, revealing exponential regimes and connections to eigenstate thermalization, with implications for experimental probing.

Contribution

It introduces a numerical linked cluster expansion approach to calculate heating rates and links these rates to the eigenstate thermalization hypothesis in both integrable and nonintegrable systems.

Findings

Heating rates follow an exponential regime.

Heating rates agree with Fermi's golden rule.

Heating rates can be experimentally probed via drive frequency.

Abstract

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.

Equations

$${\dot{E}}_m(\tau )=2\pi g_m^2\sum _{\begin{array}{c}i,f\end{array}}$$

$${\hat{H}}_0=\sum _i[\left(-t{\hat{b}}_i^{†}{\hat{b}}_{i+1}-t^{\prime }{\hat{b}}_i^{†}{\hat{b}}_{i+2}+h{\hat{b}}_i^{†}\right)+\text{H.c.}$$

$$+V\left({\hat{n}}_i-\frac{1}{2}\right)\left({\hat{n}}_{i+1}-\frac{1}{2}\right)+V^{\prime }\left({\hat{n}}_i-\frac{1}{2}\right)\left({\hat{n}}_{i+2}-\frac{1}{2}\right)],$$

$$\hat{K}=-\sum _i\left({\hat{b}}_i^{†}{\hat{b}}_{i+1}+\text{H.c.}\right),$$

$$\displaystyle K^{(s)}_{i,f}=K^{(s)}(E)\delta_{i,f}+[D^{(s)}(E)]^{-1/2}f^{(s)}_{K}(E,\omega)R_{i,f},\ \$$

$${\Gamma }_m=$$

$${\Gamma }_m^{\infty }=\frac{2\pi (m\Omega g_m{)}^2}{\text{Tr}({\hat{H}}_0^2)}\sum _s|f_K^{(s)}(E_{\infty },m\Omega ){|}^{2}Z(s),$$

$${\dot{E}}_{\Delta E,m}(\tau )=\frac{2\pi g_m^2}{\Delta E}\sum _i{P}_i^0(\tau )\times \sum _{|E_f^0-E_i^0\pm m\Omega |\le \Delta E/2}{\left|⟨E_f^0|\hat{K}|E_i^0⟩\right|}^2\left(E_f^0-E_i^0\right),$$

$${\Gamma }_{\Delta E,m}(\tau )=\frac{{\dot{E}}_{\Delta E,m}(\tau )}{E_{\infty }-E(\tau )},$$

$${\dot{E}}_m(\tau )=2\pi g_m^2\sum _s\left(\sum _{i,f\in s}|K_{i,f}^{(s)}{|}^2(E_f^0-E_i^0)P_i^0(\tau )\delta (E_f^0-E_i^0\pm m\Omega )\right)$$

$${\dot{E}}_m(\tau )=2\pi m\Omega (g_m{)}^2\sum _{(s)}$$

$${\dot{E}}_m(\tau )=4\pi m\Omega (g_m{)}^2\sinh \left[\frac{\beta (\tau )m\Omega }{2}\right]\sum _s{\int }_{E_{\text{min}}^{(s)}+m\Omega /2}^{E_{\text{max}}^{(s)}-m\Omega /2}\frac{dE{e}^{-\beta (\tau )E}}{Z(\tau )}{\left|f_K^{(s)}(E,m\Omega )\right|}^2\frac{D^{(s)}(E+m\Omega /2)D^{(s)}(E-m\Omega /2)}{D^{(s)}(E)}.$$

$$\sinh \left[\frac{\beta (\tau )m\Omega }{2}\right]\simeq \frac{\beta (\tau )m\Omega }{2},e^{-\beta (\tau )E}\simeq 1,Z(\tau )\simeq \text{Tr}[1],\text{and}[E_{\infty }-E(\tau )]\simeq \frac{\beta (\tau )\text{Tr}[{\hat{H}}_0^2]}{\text{Tr}[1]}.$$

$$\overline{|f_K^{(s)}(E_{\alpha },{\omega }_{\alpha }){|}^2}=D^{(s)}(E_{\alpha })\left(\frac{{\sum }_{j,k}|⟨E_j^0|\hat{K}|E_k^0⟩{|}^2}{{\sum }_{j,k}1}\right),$$

$${\Gamma }_m^L=\frac{2\pi (m\Omega g_m{)}^2}{\text{Tr}({\hat{H}}_0^2)}\sum _s\sum _{\alpha }\Delta E\overline{|f_K^{(s)}(E_{\alpha },m\Omega ){|}^2}\frac{D^{(s)}(E_{\alpha }+m\Omega /2)D^{(s)}(E_{\alpha }-m\Omega /2)}{D^{(s)}(E_{\alpha })},$$

$${\Gamma }_m^{\infty }=\frac{2\pi (m\Omega g_m{)}^2}{\text{Tr}({\hat{H}}_0^2)}\sum _s\overline{|f_K^{(s)}(0,m\Omega ){|}^2}Z(s).$$

$${\delta }_{\text{NLCE}}^l(\tau )=\frac{|e_l(\tau )-e_{l=17}(\tau )|}{|e_{l=17}(\tau )|}\text{and}{\delta }_{\text{ED}}^L(\tau )=\frac{|e_L(\tau )-e_{l=17}(\tau )|}{|e_{l=17}(\tau )|}.$$

$${\Delta }_l(\tau )=\frac{|e_l(\tau )-\text{fit}(\tau )|}{|\text{fit}(\tau )|}\text{and}{\Delta }_L(\tau )=\frac{|e_L(\tau )-\text{fit}(\tau )|}{|\text{fit}(\tau )|},$$

Full text

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 Dalibard et al. (2011); Kitagawa et al. (2011); Aidelsburger et al. (2013); Goldman et al. (2014), as well as to engineer topologically nontrivial band structures Oka and Aoki (2009); Lindner et al. (2011); Rechtsman et al. (2013); Cooper et al. (2019) and novel time-crystalline phases Else et al. (2016, 2017); Khemani et al. (2016); Yao et al. (2017a); *Yao_2017_Erratum; Zhang et al. (2017); Choi et al. (2017). 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 D’Alessio and Rigol (2014); Lazarides et al. (2014).

Driving at high frequencies, because of prethermalization Moeckel and Kehrein (2008); *moeckel_kehrein_2009; Eckstein et al. (2009); Kollar et al. (2011); Tavora and Mitra (2013); Tavora et al. (2014); Nessi et al. (2014); Essler et al. (2014); Bertini et al. (2015); *bertini2016prethermalization; Canovi et al. (2016); Fagotti and Collura ; Lange et al. (2018); Lenarčič et al. (2018); Reimann and Dabelow (2019); Mallayya et al. (2019), has been proposed to slow down heating Abanin et al. (2015); Else et al. (2017); Machado et al. ; Weidinger and Knap (2017); Kuhlenkamp and Knap . It results in initial fast prethermal dynamics towards time-periodic steady states (prethermal states) of effective local Hamiltonians Abanin et al. (2017a, b); Mori et al. (2016); Kuwahara et al. (2016), before thermalization dynamics eventually results in featureless “infinite-temperature” states D’Alessio and Rigol (2014); Lazarides et al. (2014); Prosen (1998); D’Alessio and Polkovnikov (2013). Prethermalization is a universal phenomenon that occurs during dynamics in isolated Mallayya et al. (2019) and open Lange et al. (2018); Lenarčič et al. (2018) systems whenever conservation laws are weakly broken. Numerical studies of prethermalization and thermalization, or, in general, of energy absorption in driven strongly interacting systems with many particles (or spins) are challenging. Progress has been achieved using massively parallel Krylov subspace methods Machado et al. , density matrix truncation Ye et al. , and t-DMRG Kollath et al. (2006), but there is a dearth of computational techniques to study generic models in arbitrary dimensions.

Here, we report on the implementation of a numerical linked cluster expansion (NLCE) for driven systems. NLCEs can be used to study arbitrary interaction strengths in arbitrary dimensions. They were originally introduced to study thermal equilibrium ensembles Rigol et al. (2006); *rigol2007numerical; *rigol2007numerical2, where they outperform full exact diagonalization calculations Iyer et al. (2015). NLCEs were recently implemented to study thermalization Rigol (2014); *rigol_16 and quantum dynamics under time-independent Hamiltonians in one Mallayya and Rigol (2018); Mallayya et al. (2019) and two White et al. ; Guardado-Sanchez et al. (2018) dimensions, and combined with dynamical quantum typicality Richter and Steinigeweg (2019). We use them to determine heating rates in strongly interacting 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.

We consider a time-periodic Hamiltonian of the form $\hat{H}(\tau)=\hat{H}_{0}+g(\tau)\hat{K}$, where $\hat{H}_{0}$ is the static Hamiltonian and $g(\tau)\hat{K}$ is a weak time-periodic perturbation of strength $g$, period $T=2\pi/\Omega$, and zero time average. The system is initialized (at $\tau=0$) in a state $\hat{\rho}_{I}=\exp[-\beta_{I}\hat{H}_{I}]/\text{Tr}\{\exp[-\beta_{I}\hat{H}_{I}]\}$ that is a thermal equilibrium state of an initial static Hamiltonian $\hat{H}_{I}$ at an inverse temperature $\beta_{I}$. At stroboscopic times $\tau=nT$ $(n=0,1,2,\dots)$, the density matrix $\hat{\rho}(\tau)$ can be written as $\hat{\rho}(\tau)=(\hat{U}_{F})^{n}\hat{\rho}_{I}(\hat{U}_{F}^{\dagger})^{n}$, where $\hat{U}_{F}=\mathcal{T}\exp[-i\int_{0}^{T}\hat{H}(t)dt]$ is the (time ordered $\mathcal{T}$) Floquet evolution operator (we set $\hbar=1$). We assume that $\hat{H}_{I}$, $\hat{H}_{0}$, and $\hat{K}$ are translationally invariant sums of local operators, and that they are mutually noncommuting (nontrivial dynamics occurs even if $g=0$).

The obvious conservation law broken by $g(\tau)\hat{K}$ is energy conservation. For sufficiently small $g$ in the thermodynamic limit, we expect prethermalization to occur (independently of the value of $\Omega$), wherein the system quickly relaxes to the equilibrium state of $\hat{H}_{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 $\hat{H}_{0}$, characterized by the instantaneous expectation values of the conserved quantities of $\hat{H}_{0}$ Mallayya et al. (2019). The dynamics of those quantities is described by autonomous equations, with drifts given by Fermi’s golden rule Mallayya et al. (2019).

We study the evolution of the energy defined by the static Hamiltonian, which is also the time-averaged Hamiltonian $\overline{\hat{H}(\tau)}=\hat{H}_{0}$, $E(\tau)=\text{Tr}[\hat{H}_{0}\hat{\rho}(\tau)]$. We consider general time-periodic perturbations, which can be Fourier decomposed as $g(\tau)\hat{K}=\sum_{m>0}2g_{m}\sin(m\Omega\tau)\hat{K}$. After a short initial transient dynamics, in the linear response regime, the system absorbs energy independently from each Fourier mode $m$. The average rate of energy absorption over a cycle is $\dot{E}(\tau)=\sum_{m>0}\dot{E}_{m}(\tau)$ with, as expected from Fermi’s golden rule,

[TABLE]

where $|E^{0}_{i}\rangle$ ($|E^{0}_{f}\rangle$) are the eigenkets of $\hat{H}_{0}$ with eigenenergies $E^{0}_{i}$ ($E^{0}_{f}$), and $P_{i}^{0}(\tau)=\langle E^{0}_{i}|\hat{\rho}(\tau)|E^{0}_{i}\rangle$ is the projection of $\hat{\rho}(\tau)$ into the basis of $\hat{H}_{0}$. The latter defines the so-called diagonal ensemble (DE) at time $\tau$ Rigol et al. (2008), $\hat{\rho}_{\text{DE}}(\tau)=P_{i}^{0}(\tau)|E_{i}^{0}\rangle\langle E_{i}^{0}|$. $\hat{\rho}_{\text{DE}}(\tau)$ is expected to characterize the equilibrated state under $\hat{H}_{0}$ at time $\tau$ D’Alessio et al. (2016). We define the rate $\Gamma(\tau)=\sum_{m>0}\Gamma_{m}(\tau)$, where $\Gamma_{m}(\tau)=\dot{E}_{m}(\tau)/[E_{\infty}-E(\tau)]$ is the rate for Fourier mode $m$, and $E_{\infty}$ is the energy at infinite temperature. Only when it is sufficiently small does one expect $|E_{\infty}-E(\tau)|$ to be an exponential function, and $\Gamma(\tau)$ to be meaningful.

We focus on 1D lattice system of hard-core bosons, with $\hat{H}_{0}$ and $\hat{K}$ given by

[TABLE]

where standard notation was used Cazalilla et al. (2011). We drive the system with a square wave $g(\tau)=g\,\text{sgn}[\sin(\Omega\tau)]$, and set $t=V=1$ (our unit of energy and frequency). $\hat{H}_{0}$ is integrable for $t^{\prime}=V^{\prime}=h=0$ (and mappable to the spin-1/2 $XXZ$ Hamiltonian Cazalilla et al. (2011)), and nonintegrable for nonvanishing $t^{\prime}$, $V^{\prime}$, and $h$. We study integrable and nonintegrable (with $t^{\prime}=V^{\prime}=0.8$ and $h=1.0$) cases, and select $\hat{H}_{I}$ to have the same terms as $\hat{H}_{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(\tau)=E(\tau)/L$ at stroboscopic times in the thermodynamic limit $(L\rightarrow\infty)$. Within NLCEs, $e(\tau)$ is expressed as a sum over the contributions of all connected clusters ($c$) that can be embedded on the lattice, $e(\tau)=\sum_{c}M(c)\times W^{e}_{c}(\tau)$, where $M(c)$ is the number of “embeddings” (per site) of cluster $c$, and $W^{e}_{c}(\tau)$ is the weight of $e(\tau)$ in cluster $c$. $W^{e}_{c}(\tau)$ is obtained recursively using the inclusion-exclusion principle: $W^{e}_{c}(\tau)=E_{c}(\tau)-\sum_{c^{\prime}\subset c}W^{e}_{c^{\prime}}(\tau)$, where $c^{\prime}$ denotes the connected subclusters of $c$ and $E_{c}(\tau)=\text{Tr}[\hat{H}_{0}^{c}\hat{\rho}_{c}(\tau)]$ is the energy in cluster $c$ [$\hat{H}^{c}_{0}$ is the static Hamiltonian, and $\hat{\rho}_{c}(\tau)$ is the density matrix at time $\tau$, both in cluster $c$]. The series starts with the smallest cluster (a site) for which $W_{c}(\tau)=E_{c}(\tau)$. For each cluster, $E_{c}(\tau)$ 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 next-nearest neighbor interactions Rigol (2014); *rigol_16; Mallayya and Rigol (2018, 2017). The order of the NLCE is set by the number of sites of the largest cluster considered. For nonintegrable $\hat{H}_{0}$, we compute 17 orders of the NLCE (after exploiting all symmetries, the dimension of largest sector of the Hamiltonian is 32 896). When $\hat{H}_{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(\tau)|$ vs $\tau$ 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 $\hat{H}_{I}$ at an inverse temperature $\beta_{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(\tau)$ to the infinite-temperature energy ($E_{\infty}/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 $\propto 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 sup .

It follows from eigenstate thermalization for nonintegrable Hamiltonians Deutsch (1991); Srednicki (1994); Rigol et al. (2008); D’Alessio et al. (2016) (generalized eigenstate thermalization for integrable Hamiltonians Cassidy et al. (2011); Vidmar and Rigol (2016)) $\hat{H}_{0}$ that the predictions of $\hat{\rho}_{\text{DE}}(\tau)$ for few-body operators agree with those of the thermal (generalized Gibbs) ensemble D’Alessio et al. (2016); Vidmar and Rigol (2016); Essler and Fagotti (2016); Caux (2016). We first focus on the case in which $\hat{H}_{0}$ is nonintegrable with no local conservation law. In this case, the inverse temperature $\beta(\tau)$ alone characterizes the thermal (grand canonical) ensemble at $\tau$, $\hat{\rho}_{\text{GE}}(\tau)=\exp[-\beta(\tau)\hat{H}_{0}]/\text{Tr}\{\exp[-\beta(\tau)\hat{H}_{0}]\}$, where $\beta(\tau)$ is determined by the condition $\text{Tr}[\hat{H}_{0}\hat{\rho}_{\text{GE}}(\tau)]=\text{Tr}[\hat{H}_{0}\hat{\rho}(\tau)]$. Only when $\beta(\tau)\ll 1$ is that one expects $\Gamma(\tau)$ to become independent of $\beta(\tau)$, and $E(\tau)$ to approach $E_{\infty}$ as a single exponential.

To illustrate this, in the main panel of Fig. 2 we plot $|e(\tau)|$ (normalized by its initial value $|e(0)|$) for various initial inverse temperatures $\beta_{I}\in[0.01,0.5]$. The normalized energies $e(\tau)/e(0)$ for $\beta_{I}=0.033$ and $0.01$ exhibit a nearly identical exponential decay (within the times at which the NLCE has converged) implying that $\Gamma$ is independent of $\beta_{I}$ [hence, of $\beta(\tau)$] when $\beta_{I}\lesssim 0.03$. For $\beta_{I}\gtrsim 0.2$, one can still use exponentials to fit $e(\tau)$, but the rates obtained depend on $\beta_{I}$. In the inset in Fig. 2, we report the rates obtained from such fits vs $\beta_{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 $\beta_{I}$ is the one expected to occur as a function of driving time for initial states that are not in the regime $\beta_{I}\ll 1$.

Next, we focus on the dependence of the heating rates on $\Omega$. In nonintegrable systems, the eigenstate thermalization hypothesis Deutsch (1991); Srednicki (1994); Rigol et al. (2008); D’Alessio et al. (2016) allows one to compute $\Gamma_{m}(\tau)$. After resolving all symmetries of the static Hamiltonian, the eigenstate thermalization hypothesis ansatz for the matrix elements $K^{(s)}_{i,f}=\langle E^{0}_{i}|\hat{K}|E^{0}_{f}\rangle$ of the operator $\hat{K}$ (used as drive) in each block diagonal sector $s$ of $\hat{H}_{0}$ has the form D’Alessio et al. (2016); Srednicki (1999)

[TABLE]

where $E=(E_{i}+E_{f})/2$, $\omega=E_{f}-E_{i}$, $D^{(s)}(E)$ is the density of states of $\hat{H}_{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,\omega)$ are smooth functions of their arguments.

Using Eqs. (1) and (4), changing sums over eigenstates by integrals over energy, replacing $\hat{\rho}_{\text{DE}}(\tau)$ by $\hat{\rho}_{\text{GE}}(\tau)$ and assuming high temperature [$\beta(\tau)\ll 1$], one obtains the following expression for the heating rate sup

[TABLE]

where $E^{(s)}_{\text{min}}$ ($E^{(s)}_{\text{max}}$) is the minimum (maximum) energy in sector $s$, and $m\Omega$ is smaller than $E^{(s)}_{\text{max}}-E^{(s)}_{\text{min}}$ (otherwise there is no linear response heating for that mode).

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) sup . NLCE results are not reported for small and large values of $\Omega$ 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 $\Omega$, we find that the evaluation of Eq. (5) results in heating rates that can be well described by an exponential in $\Omega$. This is consistent with rigorous bounds Abanin et al. (2015, 2017b); Else et al. (2017).

When $\hat{H}_{0}$ is integrable (the spin-1/2 $XXZ$ limit), the prethermal states are described by a generalized Gibbs ensemble (GGE) $\hat{\rho}_{\text{GGE}}(\tau)$ Wouters et al. (2014); Pozsgay et al. (2014); Ilievski et al. (2015). When $\hat{\rho}_{I}$ is a thermal state with $\beta_{I}\ll 1$ (or in general after long driving times), $\hat{\rho}_{\text{GGE}}(\tau)\simeq\hat{\rho}_{\text{GE}}(\tau)$ with $\beta(\tau)\ll 1$ He and Rigol (2012). In this regime, Eq. (5) gives the heating rates for the integrable static Hamiltonian provided that there is a well defined $|f_{K}^{(s)}(E,\omega)|^{2}$. In 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 $\Omega$ in the nonintegrable and integrable cases, the heating rates in the latter are described by Eq. (5) and, for high $\Omega$, they are well described by an exponential in $\Omega$.

The previous results show that heating rates can be used to probe the function $f_{K}^{(s)}(E,m\Omega)$ in nonintegrable and integrable systems. Still, Eq. (4) involves the density of states. For large system sizes, since $E$ is extensive but $\Omega$ is not, $D^{(s)}(E+m\Omega/2)D^{(s)}(E-m\Omega/2)\simeq[D^{(s)}(E)]^{2}$ and $E^{(s)}_{\text{min,max}}\pm m\Omega/2\simeq E^{(s)}_{\text{min,max}}$. Using the saddle point approximation to compute the integral in Eq. (5), and using that $D^{(s)}(E_{\infty})$ is maximal, the heating rate for Fourier mode $m$ in the thermodynamic limit $({\Gamma}^{\infty}_{m})$ can be written as

[TABLE]

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 $\Omega$, gives the average $|f_{K}^{(s)}(E_{\infty},\Omega)|^{2}$ over all sectors of the Hamiltonian in the thermodynamic limit sup .

In summary, we studied heating in strongly interacting driven lattice systems and showed that, at sufficiently high effective temperatures ($[\beta(\tau)]^{-1}\gtrsim 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_{K}^{(s)}(E,\Omega)|^{2}$ in integrable interacting systems. This has been confirmed in a recent full exact diagonalization study of the spin-1/2 $XXZ$ chain LeBlond et al. , and needs to be further explored to place it on equal footing with what is known for quantum chaotic systems Khatami et al. (2013); Steinigeweg et al. (2013); Beugeling et al. (2015); D’Alessio et al. (2016); Luitz and Bar Lev (2016); Mondaini and Rigol (2017); Jansen et al. (2019).

S1 S1. Numerical evaluation of Eq. (1) in the main text

Equation (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 $\Delta E$, Eq. (1) is modified to the following expression (which is amenable to numerical evaluation)

[TABLE]

where $|E^{0}_{i}\rangle$ $(|E^{0}_{f}\rangle)$ are eigenkets of $\hat{H}_{0}$ with eigenenergies $E^{0}_{i}$ ($E^{0}_{f}$), and $P_{i}(\tau)=\langle E^{0}_{i}|\hat{\rho}(\tau)|E^{0}_{i}\rangle$. With this coarse graining procedure, $\Gamma_{\Delta E,m}(\tau)$ for Fourier mode $m$ is calculated as

[TABLE]

where $E(\tau)$ is also evaluated using full exact diagonalization, and $E_{\infty}=0$ for our model. $\Gamma_{\Delta E}(\tau)=\sum_{m>0}\Gamma_{\Delta E,m}(\tau)$ is the relaxation rate of $E(\tau)$.

In Fig. S1, we show $\Gamma_{\Delta E}(\tau)$ vs $\tau$ for three values of $\Delta 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 $\beta_{I}=(30)^{-1}$. It is apparent in Fig. S1 that $\Gamma_{\Delta E}(\tau)$ is nearly constant, with a slight drift at long times (apparent for $g=0.8$), and that it is independent of the value of $\Delta E$. We identify a range of $\Delta E$ and $\tau$ where $\Gamma_{\Delta E}(\tau)$ is (nearly) constant [and where the dynamics of $e(\tau)$ is exponential and robust against finite-size effects, see Sec. S3] and compute the average of $\Gamma_{\Delta E}(\tau)$ in this range.

In the main text, all the rates reported in Figs. 1 and 2 for which Eq. (1) was used were obtained averaging over $\Delta E/L=\{0.002,0.004,\dots,0.04\}$ and $\tau=\{3,4,\dots,10\}$ (a total of 160 values) for the nonintegrable $\hat{H}_{0}$, and $\Delta E/L=\{0.002,0.004,\dots,0.02\}$ and $\tau=\{2,3,\dots,6\}$ (a total of 50 values) for the integrable $\hat{H}_{0}$. The standard deviation of the averages were reported as error bars.

S2 S2. Derivation of Eq. (5) in the main text

Equation (1), accounting for the block diagonalization of $\hat{H}_{0}$ in symmetry sectors $\{s\}$, has the form

[TABLE]

where $K^{(s)}_{i,f}=\langle E^{0}_{i}|\hat{K}|E^{0}_{f}\rangle$, for $i,f\in s$, and $P_{i}^{0}(\tau)=\langle E^{0}_{i}|\hat{\rho}(\tau)|E^{0}_{i}\rangle$. From the eigenstate thermalization hypothesis (ETH) it follows that the results from the diagonal ensemble and the Gibbs ensemble agree Deutsch (1991); Srednicki (1994); Rigol et al. (2008); D’Alessio et al. (2016), so one can replace $P_{i}^{0}(\tau)$ by $\exp[-\beta(\tau)E^{0}_{i}]/Z(\tau)$, where $Z(\tau)$ is the partition function, and the inverse temperature $\beta(\tau)$ is set by the energy $E(\tau)$.

Using the Gibbs ensemble, the ETH ansatz for $K^{(s)}_{i,f}$ (see main text), and replacing sums by integrals, Eq. (S3) can be written as

[TABLE]

where $E^{(s)}_{\text{min}}$ ($E^{(s)}_{\text{max}}$) is the minimum (maximum) energy in sector $s$, and we used that $|f_{K}^{(s)}(E,-\omega)|=|f_{K}^{(s)}(E,\omega)|$. A change of variable $E\rightarrow E+m\Omega/2$ in the first integral, and $E\rightarrow E-m\Omega/2$ in the second integral, allows one to rewrite the expression above as

[TABLE]

At high temperatures, when $\beta(\tau)\ll 1$, one has to lowest order in $\beta(\tau)$

[TABLE]

Using Eqs. (S5) and (S6), the heating rate $\Gamma_{m}=\dot{E}_{m}(\tau)/[E_{\infty}-E(\tau)]$ reduces to Eq. (5) in the main text.

S2.1 Numerical evaluation of Eq. (5) in the main text

Like Eq. (1), Eq. (5) in the main text is evaluated using full exact diagonalization of chains with $L$ sites and periodic boundary conditions. We define a small energy window $\Delta E$, which we use to bin the spectrum of $\hat{H}_{0}$ in each symmetry sector $s$. Each bin $\alpha$, with energy $E_{\alpha}$, includes all eigenstates with eigenenergies $E_{i}^{0}\in(E_{\alpha}-\Delta E/2,E_{\alpha}+\Delta E/2)$. The density of states at energy $E_{\alpha}$ is then $D^{(s)}(E_{\alpha})=n_{\alpha}/\Delta E$, where $n_{\alpha}$ is the number of energy eigenstates in bin $\alpha$. The function $|f_{K}^{(s)}(E_{\alpha},\omega_{\alpha})|^{2}$, with $\omega_{\alpha}>0$, after coarse graining is given by

[TABLE]

where $E^{0}_{j}$ and $E^{0}_{k}$ are such that bin $\alpha_{j}$ containing $E_{j}^{0}$ and $\alpha_{k}$ containing $E_{k}^{0}$ satisfy $(E_{\alpha_{j}}+E_{\alpha_{k}})/2\in(E_{\alpha}-\Delta E/2,E_{\alpha}+\Delta E/2)$ and $|E_{\alpha_{j}}-E_{\alpha_{k}}|=\omega_{\alpha}$. This coarse graining procedure modifies Eq. (5) in the main text to

[TABLE]

where the inner sum is over all the bins $\alpha$ whose energy $E_{\alpha}\in(E^{(s)}_{\text{min}}+m\Omega/2,E^{(s)}_{\text{max}}-m\Omega/2)$.

In contrast to Eq. (S1), Eq. (S8) does not involve calculating the time evolution of the system. As a result, we are able to evaluate Eq. (S8) in chains with $L=19$ ($L=21$) for the nonintegrable (integrable) static Hamiltonian. The dimension of the largest symmetry resolved sector is 13,797 (16,796) for the nonintegrable (integrable) $\hat{H}_{0}$.

In Fig. S2, we show heating rates for the $m=1$ mode, $\Gamma^{L}_{m=1}$, evaluated at $\Omega=\Delta E,2\Delta E,\dots$ for two values of $\Delta E$ for the nonintegrable and the integrable static Hamiltonians. Our values of $\Delta E/L$ are such that the spectrum of $\hat{H}_{0}$ is divided into $10L$ bins [(a) $\Delta E/L\sim 0.014$ and (b) $\Delta E/L\sim 0.006$] and $40L$ bins [(a) $\Delta E/L\simeq 0.004$ and (b) $\Delta E/L\simeq 0.002$]. The results obtained can be seen to be robust against the choice of $\Delta E$. For the results reported in Fig. 3 of the main text, we use $\Delta E/L\simeq 0.004$ ($40L$ bins) to evaluate $\Gamma^{18}$, $\Gamma^{19}$, and $\Gamma_{m=1}^{19}$ in Fig. 3(a) for the nonintegrable static Hamiltonian, and $\Delta E/L\simeq 0.002$ ($40L$ bins) to evaluate $\Gamma^{20}$, $\Gamma^{21}$, and $\Gamma_{m=1}^{21}$ in Fig. 3(b) for the integrable static Hamiltonian.

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

[TABLE]

It is apparent in Fig. S2, both for the nonintegrable and the integrable static Hamiltonians, that the results for $\Gamma^{\infty}_{m=1}$ calculated using Eq. (S9) do not agree with the ones for $\Gamma^{L}_{m=1}$ using Eq. (S8). This is because of strong finite-size effects in $\Gamma^{\infty}_{m=1}$. We note that the disagreement increases as $\Omega$ increases. The fact that finite-size effects in $\Gamma^{\infty}_{m=1}$ increase with increasing $\Omega$ is also apparent in the increasing discrepancy with increasing $\Omega$ between the results for the two chain sizes shown in Fig. S2. This is in contrast to the results for $\Gamma^{L}$ (similar to $\Gamma^{L}_{m=1}$ at large $\Omega$) evaluated from Eq. (S8) and reported in Fig. 3 in the main text for two systems sizes. The strong finite-size effects in $\Gamma^{\infty}_{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 S3. Convergence of NLCE and exact diagonalization

In Fig. S3(a), we plot the absolute value of the energy per site $|e(\tau)|$ of the nonintegrable Hamiltonian for $g=0.5$, period $T=1$, and $\beta_{I}=30^{-1}$. The results shown for $|e(\tau)|$ 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 ($\tau\lesssim 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 $\tau\lesssim 10$ are the reason we evaluate Eq. (1) of the main text with ED in the range $3\leq\tau\leq 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 $\tau$ [shown in Fig. S3(a) as a black line]. For $\tau\gtrsim 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 $\tau\sim 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(\tau)$ evaluated with the $l^{\text{th}}$ NLCE order [$e_{l}(\tau)$] and with the $L$-site periodic chain ED [$e_{L}(\tau)$], we define the convergence errors for $l<17$ and $L$, respectively, as the relative differences from $e_{l=17}(\tau)$ given by

[TABLE]

Fig. S3(b) reports $\delta^{l}_{\text{NLCE}}(\tau)$ and $\delta^{L}_{\text{ED}}(\tau)$ vs $l$ and $L$, respectively, at $\tau=10$ and $\tau=12$, for the same $e(\tau)$ as in Fig. S3(a). The plots make apparent that the errors decrease with increasing $l$ and $L$, that the errors at $\tau=12$ are greater than the corresponding ones at $\tau=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 Iyer et al. (2015).

At $\tau=10$, the estimated of error for the ED calculation [$\delta^{L}_{\text{ED}}(10)$] for $L=18$ in Fig. S3(b) is less than 0.2%. For all calculations with ED in this paper, $\tau=10$ is the largest time considered for the nonintegrable $\hat{H}_{0}$.

As argued in the main text, $e(\tau)$ in the thermodynamic limit is essentially a single exponential in $\tau$ for $\beta_{I}\lesssim 30^{-1}$. Hence, for $\tau\gtrsim 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}(\tau)$ for $3\leq\tau\leq 10$. The normalized deviation from the exponential fit for the $l^{\text{th}}$ NLCE order [$\Delta_{l}(\tau)$] and the $L$-site periodic chain solved with ED [$\Delta_{L}(\tau)$], are given by

[TABLE]

where fit($\tau$) is the value of the exponential fit at $\tau$. In Fig. S3(c), $\Delta_{l}(\tau)$ and $\Delta_{L}(\tau)$ are plotted vs $l$ and $L$, respectively, for $\tau=22$ and $\tau=25$. It is apparent that the NLCE results converge faster towards the exponential [at $\tau=22$, $\Delta_{l=17}(\tau)$ 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 $\tau\leq 20$ the nonintegrable $\hat{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.