Floquet prethermalization in a Bose-Hubbard system

Periodic driving has emerged as a powerful tool in the quest to engineer new and exotic quantum phases. While driven many-body systems are generically expected to absorb energy indefinitely and reach an infinite-temperature state, the rate of heating can be exponentially suppressed when the drive frequency is large compared to the local energy scales of the system --- leading to long-lived 'prethermal' regimes. In this work, we experimentally study a bosonic cloud of ultracold atoms in a driven optical lattice and identify such a prethermal regime in the Bose-Hubbard model. By measuring the energy absorption of the cloud as the driving frequency is increased, we observe an exponential-in-frequency reduction of the heating rate persisting over more than 2 orders of magnitude. The tunability of the lattice potentials allows us to explore one- and two-dimensional systems in a range of different interacting regimes. Alongside the exponential decrease, the dependence of the heating rate on the frequency displays features characteristic of the phase diagram of the Bose-Hubbard model, whose understanding is additionally supported by numerical simulations in one dimension. Our results show experimental evidence of the phenomenon of Floquet prethermalization, and provide insight into the characterization of heating for driven bosonic systems.


I. INTRODUCTION
The study of out-of-equilibrium dynamics in Floquet systems is an exciting new frontier in quantum physics [1][2][3][4][5]. By driving a quantum system, it is possible to enhance or stabilize interesting equilibrium phases, or even to create new, inherently non-equilibrium phases without a static analogue, such as the discrete time crystal (DTC) [6][7][8][9][10][11] or the anomalous Floquet insulator [12][13][14]. A seemingly ubiquitous obstruction towards realizing such phases in the many-body setting is thermalization: by absorbing energy from the drive, a quantum system is generically expected to heat up and eventually approach a featureless "infinite temperature" state, which is the maximum entropy state in the absence of any conservation laws [15][16][17]. The only robust exception to this fate is provided by many-body localization (MBL) [18][19][20][21][22][23][24][25][26][27][28], whereby sufficiently strong disorder can prevent this 'heat death' [15,17,[29][30][31][32]. This exception however comes with a number of constraints (e.g. on the presence of disorder, the dimensionality of the system or the range of interactions [33,34]) that may preclude potentially interesting theoretical scenarios or experimental platforms.
An alternative route towards the realization of nonequilibrium phases is to accept the ultimate thermalizing fate of the system, and focus instead on delaying the inevitable by engineering a long-lived 'prethermal' regime. In particular, it has been shown that the timescale for * These two authors contributed equally to this work heating can be bounded from below as t th O(e ω/J eff ) for sufficiently large drive frequency ω, where J eff represents a typical local energy scale of the system [35][36][37][38][39]. Intuitively, such exponentially large heating timescales arise when absorbing one quantum of energy from the drive requires the rearrangement of many local degrees of freedom, which is a high-order process. At times t t th , the system can in principle exhibit rich dynamics, featuring symmetries, quasi-conserved quantities (including an effective Hamiltonian), etc. [37,40,41] The existence of prethermal regimes at large driving frequencies was established in a number of analytical and numerical works [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49]. On the experimental front, the heating mechanisms of interacting driven systems have lately been the focus of several works based on ultracold atoms in optical lattices, tracking condensate decay, atom losses or doublon production [50][51][52][53][54][55]. Recently, signatures of prethermalization were observed in the extreme limit of fast, ultra-strong driving [56] by measuring the lowest-band depletion. However, a clear experimental demonstration of frequency-tuned prethermalization, manifested by exponential-in-frequency heating times, has yet to be provided. Such a demonstration requires overcoming several challenges: the quantum system, while driven, must preserve coherence for long enough times that the exponential scaling of t th with frequency becomes manifest; at the same time, one must be able to tune ω across a wide enough range without exciting high-energy degrees of freedom that lie outside the scope of the original isolated system -e.g. omnipresent higher bands in lattice models; and, finally, the overall exponential trend must be resolved from other system-specific spectral features that may obscure it in the available frequency range.
In this work we show evidence of an exponential-infrequency thermalization time, the main signature of Floquet prethermalization, in a driven Bose-Hubbard system of ultracold bosonic atoms in a square optical lattice. This observation is made possible by two crucial advantages of our experimental setup: (i) the high degree of isolation of our system, which enables us to explore long evolution times [57,58]; and (ii) a sensitive thermometry technique, enabled by quantum-gas microscopy [59], which allows us to measure the heating even for weak drives, thereby suppressing transfers to higher bands, and avoiding parametric instabilities.
By tuning the lattice parameters, we explore the response of the atoms in a broad range of couplings spanning the superfluid and Mott-insulating phases. The exponentially slow heating is observed most clearly on the superfluid side, where it coexists with weak spectral features possibly associated to Bogoliubov quasiparticle excitations. On the Mott-insulating side, the heating rate is non-monotonic in frequency, dominated by resonances with higher-occupation excitations (doublons and triplons). Nonetheless, in both phases, the heating rate drops substantially (by 1 to 2.5 orders of magnitude) even with a modest increase in frequency in the experimentally accessible range. Our understanding of the observed phenomena is aided by numerical simulations in one dimension which, while limited in system size, can explore a broad range of couplings, drive frequencies and timescales.

II. EXPERIMENTAL SETUP
Our experiment begins with the preparation of a twodimensional cloud of ultracold 87 Rb atoms trapped in a single antinode of a vertical optical lattice. Next, the cloud is adiabatically loaded into an in-plane square optical lattice at depth V 0 . We fix the number of atoms such that the density in the central part of the trap is close to one atom per lattice site, typically leading to a total atom number of N at 200. In the prepared state, the atoms only populate the lowest energy band of the lattice potentials, and our system can be well described by a two-dimensional Bose-Hubbard model, with Hamiltonian whereâ i ,â † i andn i respectively denote the annihilation, creation and number operators at a site i of the square lattice [i = (i x , i y )], J is the tunneling amplitude between nearest-neighbour sites i, j , U the on-site interaction energy, and i the harmonic trapping potential (see Appendix A).
At this stage, the atoms are at very low temperature, close to the ground state of Eq. (1). We then start driving the system by periodically modulating the depth of the in-plane lattices as V (t) = V 0 (1 + A cos(ωt)) (see Fig. 1(a)), where A is the relative modulation amplitude, leading to a time-dependent modulation of all the Hamiltonian parameters. However, the tunneling strength, due to its exponential dependence on the lattice depth, dominates the modulation so that with g = δJ/J andÔ drv = J i,j â † iâj [60] . To ensure that during the driving no atoms are excited into higher bands, we keep the driving frequencies well below the bandgap [61] and use a low modulation amplitude A 1 to avoid multi-photon interband transitions [62] (see Appendix A). After driving the system for an integer number of periods, N cyc = ω t drv /2π, we slowly ramp up the lattice depth until the system becomes an atomic-limit Mott insulator. At this stage, all tunneling dynamics is frozen and, if no heating took place during the drive, this results in a Mott insulator near unit filling. Finally, we measure the parity-projected atomic occupation of the lattice sites through fluorescence imaging [59]. Because of parity projection, the growth of the hole density is directly linked to excitation processes that increase the variance of the single-site occupation. Thus, the density of holes is a proxy for the energy density of the cloud and thus for the heating dynamics (see App. D).

III. THERMALIZATION DYNAMICS
In Fig. 1(c), we plot the evolution of the density of holes ρ h under driving for four different frequencies. The depth of the in-plane lattices is V 0 = 8 E r , where E r is the recoil energy, and the modulation amplitude is A = 0.05. While each drive frequency in Fig. 1(c) gives rise to a qualitatively similar time dependence of ρ h , the thermalization times themselves are vastly different -spanning 1.5 orders of magnitude within less than a factor of 2 in ω. Such a striking dependence on the drive frequency is an indication of the exponential slow-down of heating characteristic of Floquet prethermalization.
Theoretically, the phenomenon of prethermalization refers to the quasi-conservation of energy stemming from the existence of a quasilocal time-independent Hamiltonian that captures the dynamics of the system out to an exponentially long time. While one can always formally define a (non-unique) 'Floquet Hamiltonian' from the unitary time evolution operator viaÛ (T ) = T e −i T 0 dtĤ(t)/ ≡ e −iĤ F T / , the operatorĤ F is generally highly non-local in a many-body system (and hence unphysical). Nevertheless, when the frequency ω is large compared to the local energy scales of the arXiv:1912.09443 problem (here denoted collectively by J eff ), one can perform a high-frequency asymptotic expansion forĤ F in powers of 1/ω,Ĥ F = n (1/ω) nĤ (n) F ; the leading-order termĤ (0) F equals the time-averaged Hamiltonian, while higher-order terms are progressively longer-ranged and contribute significantly to the dynamics only at correspondingly later times. While ultimately divergent, this expansion looks convergent out to some optimal order n opt = O(ω/J eff ). Truncating the expansion at this order yields an exponentially accurate approximation to the Floquet time evolutionÛ (T ), which sets the rate of heating to be exponentially small [35,[37][38][39], giving where E(N cyc ) = Ĥ (0) F Ncyc is the energy absorbed by the system.
Thus, in our experiment one expects the energy density to at first increase linearly in time (as is also expected, e.g., from linear response theory), before eventually saturating to its infinite-temperature value; other local observables such as ρ h are expected to follow the same behavior. While more intricate behavior could in principle arise at intermediate times, the simplest ansatz for the stroboscopic time dependence of ρ h takes the form where ρ 0 is the baseline value measured in the absence of the drive, ρ ∞ is the infinite-temperature value, and N th cyc is the heating timescale, predicted to obey N th cyc > e O(ω) for sufficiently high ω.
Fits to the data of the form (4) (solid lines in Fig. 1(c)) in fact reveal a good agreement, and yield thermalization timescales N th cyc between 4 × 10 2 at the lowest frequency and 10 4 at the highest. Note that the long times mea-sured in our data, over 3000 /J, are crucial to the detection of this effect and reflect the high level of isolation of our system. We remark that what may look like a plateau in ρ h at early times is in fact a very slow linear growth as follows from Eq. 3 (and as is true more generally of prethermal "plateaux"). A nonequilibrium initial state, e.g. one with a spatial density imbalance, would first thermalize relative to the leading-order prethermal Hamiltonian (in our caseĤ (2)) on a timescale independent of ω, and then heat to infinite temperature exponentially slowly. However, in our experiment we do not observe any fast transient dynamics before the onset of this slow heating because the initial state we prepare is already in thermal equilibrium relative toĤ While the dynamics of ρ h illustrate qualitatively the phenomenon of Floquet prethermalization, a more precise characterization of the heating-rate dependence can be obtained from a fit of the atomic density profile. This established thermometry method, based on the fit of a grand-canonical model [59,63], provides us with a measure for the temperature of the cloud, from which we can characterize the heating induced by the drive. Due to the high sensitivity of this technique, we can better explore the linear response regime, where the drive amplitude A is small such that interband processes are strongly suppressed (see Appendix B). This weak-drive probing is in contrast to recent measurements of the response of Bose-Einstein condensates in one and two-dimensional optical lattices [54,55], which focused instead on the emergence of parametric instabilities under strong drives.

IV. HEATING IN THE BOSE-HUBBARD MODEL: NUMERICS
Before moving on to the results of the experiment outlined above, it is useful to gain some intuition on the nature of Floquet heating in the Bose-Hubbard phase diagram with the help of numerical simulations. A variety of methods have been applied to the Bose-Hubbard model out of equilibrium [1,45,[64][65][66][67], with particular interest on parametric instabilities of the superfluid condensate in recent years [42,68]. To study the approach to infinite temperature under weak driving, we use numerical exact diagonalization and the Krylov subspace method for time evolution [69] as detailed in the following.
While the general theory of prethermalization applies to arbitrarily strong drives, our experiment considers a weak modulation g 1. In this regime, the energy absorbed per Floquet cycle by the system (simply called 'heating rate' hereafter) is well captured within linear response theory as the dissipative part of a response function, where {E n , |n } label the eigenvalues and eigenvectors of the time-averaged HamiltonianĤ 0 in Eq. (1) (|0 being the ground state, with E 0 = 0), andÔ drv has been introduced in Eq. (2). We note that Φ(ω) is the quantity rigorously bounded by an exponential in Ref. [35]. It has units of energy and, for weak drives (g 1), is proportional to the energy absorbed per Floquet cycle, We compute Φ(ω) in Eq. (5) via numerical exact diagonalization of a one-dimensional chain of L = 9 sites at unit filling. Given the small system size, we replace the harmonic trap potential inĤ 0 with hard-wall (open) boundary conditions (for additional details on the numerics see Appendix D). This approach, while strongly limited in system size, provides complete flexibility in the choice of couplings J/U and frequency ω, while also allowing us to probe extremely long timescales (within linear response). The results, shown in Fig. 2, outline a clear picture of the nature of heating in the two phases. To best highlight each phase's spectral features, we show the heating and the frequency in units of J in the superfluid phase (see Fig. 2(a)), and U in the Mott insulating phase (see Fig. 2(b)).
Deep in the superfluid phase, the system manages to heat efficiently up to frequency ω = 8J/ , then the rates sharply drop in an exponential fashion, with additional kink-like features visible at ω = 16J/ and higher multiple frequencies. This suggests that heating takes place via the excitation of quasiparticles from the condensate. As the drive carries no net momentum, the quasiparticles must come in pairs with opposite momenta ±q.
For ω > Ω 2qp (twice the quasiparticle bandwidth), to absorb a single quantum of energy ω from the drive, the system must concurrently scatter multiple pairs from the condensate into excited states, with each additional scattering event suppressing the amplitude by factors of U/J 1. This gives the expected exponential scaling and explains the threshold features visible in Fig. 2(a) at frequencies commensurate with 8J/ (twice the singleparticle bandwidth). Increasing the interaction strength U gradually washes out the above features, while pushing the value of Ω 2qp upwards approximately in agreement with the Bogoliubov prediction Ω 2qp 8J/ 1 + U/2J -though notice the latter is a mean-field prediction and as such it is expected to work better in higher dimension.
As U increases further, eventually the superfluid spectral features give way to sharp peaks associated to higher on-site occupancy (doublons, triplons, etc.), a characteristic of the Mott-insulating phase. In fact, Fig. 2(b) displays a hierarchy of peaks at integer values of ω/U . The location and height of each peak can be understood by perturbing away from an atomic limit (J = 0) Mott state (see Appendix C). While the envelope of the peaks does obey an exponential bound (in one dimension -see Appendix C for a discussion of higher dimension), the strongly non-monotonic structure of Φ(ω) means that measurements with a limited dynamic range in ω and N th cyc may not be able to identify the overall trend. In Appendix D we show that the above picture remains qualitatively valid even away from the linear response regime by studying the approach of the system to infinite temperature via exact time evolution.

V. EXPERIMENTAL RESULTS: HEATING IN TWO DIMENSIONS
We now turn to experimentally characterizing the slow Floquet thermalization, by extracting the temperatures from the density profile of the cloud, in our twodimensional system. We measure the heating rates at different lattice depths within V 0 = 5 − 11E r over a range of driving frequencies with a fixed relative modulation A = 0.05. We express it as the energy absorbed per Floquet cycle φ(ω), which is related to Φ(ω) through φ(ω) = (2πg) 2 Φ(ω). The results (see Fig. 3(a)) reveal a clear suppression of the heating rate as the frequency is increased, extending over more than two decades in the measured range. This abrupt arrest of Floquet thermalization manifests the presence of a Floquet prethermal regime.
For values of J/U well above the phase transition at J/U 0.06 [70], i.e. in the superfluid phase, all datasets show qualitatively the same behavior -a monotonic and approximately exponential decrease of the heating rate as the frequency is increased. This trend is further indicated, for the two datasets with highest J/U , by the fit of an exponential function φ(ω) = C e − ω/J eff (see dotted lines in Fig. 3(a)). The fit allows us to extract the effective local energy scales, J eff,1 = 5.76(16) J and J eff,2 = 5.9(2) J, which are consistently on the same order of magnitude as the Hamiltonian parameters J and U . As we move away from weak interactions, a visible kinklike feature appears on top of the general exponential trend. Following the line of argument from the previous section, we expect the dominant heating process in the superfluid to be the generation of quasiparticle pairs with opposite momenta. Because of this, the heating is expected to be further reduced for drive frequencies above twice the Bogoliubov quasiparticle bandwidth, which in the 2D case is Ω 2qp,2D = 2 × 8J/ 1 + U/4J [42]. In Fig. 3(a) we indicate with five small arrows the frequency Ω 2qp,2D for the four traces with weakest interactions, which roughly agree with the positions of the kinks observed in the data.
Aside from this feature, we note that below twice the non-interacting bandwidth (here 16J/ ) the heating rate is not flat, in contrast to the results for the 1D numerics. This is to be expected based on the different density of states of a tight-binding model in a square lattice in 1D and 2D: in the latter case, the density of states has a maximum in the middle of the band, which means the excitation of quasiparticles is most efficient, and the heating fastest, near twice the middle of the bandwidth (here ω = 8J/ ) [71]. Finally, we also note that the dynamic range of the driving frequencies is smaller for higher J/U , since higher tunneling strengths J require higher frequencies (on an absolute scale) to see prethermalization, making the limitation posed by interband heating more severe.
Moving to the strongly interaction regime, different features emerge. The dataset with the smallest J/U , in fact the only one in the Mott-insulating phase, shows a distinct non-monotonic behavior in the observed frequency range. To identify the relevant spectral features associated with it, in Fig. 3(b) we show the same data as in (a) but expressed in units of the interaction strength U . Indeed the trace at J/U = 0.06 (dataset in dark orange) shows a peaked structure at ω = U/ and 3 U/ , as expected for the doublon and triplon resonances respectively. As the interaction strength is reduced, these resonant features fade into a continuum associated with the superfluid bandwidth, similar to what one observes for the numerics in Fig. 2(b).
Finally, we also note that in the regime of very high frequencies, which features the smallest heating rates, we reach the sensitivity limit of this experiment. This is caused by the very long measurements times and the then significant contribution of the background heating present in our system. The noise floor is expected to be dependent on J/U due to the change in the excitation spectrum of the system.

VI. EXPERIMENTAL RESULTS: HEATING IN ONE DIMENSION
Our experimental setup allows us also to produce 1D systems. We achieve this by ramping the lattice along the y axis to a depth of V 0,y = 20E r before the start of the drive. The typical atom number in each one of these 1D systems is of N at 15. The measured heating rates in this 1D geometry are shown in Fig. 4, displayed in a similar fashion as in the 2D case, for lattice depths at V 0,x = 3 − 9E r and with a relative lattice modulation of A x = 0.1 (while A y = 0). Here too we observe an exponential suppression of the heating rate as a function of the drive frequency ω. However, for ω < 8J/ the heating rate appears almost constant, in agreement with the numerics in Fig. 2(a), and only beyond this flat part we see a sharp decrease. While this behavior can be solely explained in terms of twice the non-interacting bandwidth of the system, 2 × 4J/ , we also observe a second kink at slightly higher frequency, reminiscent of the 2D case, which shifts to higher frequencies as interactions increase. In Fig. 4 we also use four arrows to indicate the position of twice the Bogoliubov bandwidth, Ω 2qp,1D = 2×4J/ 1 + U/2J, for the first four datasets, showing a reasonable agreement that eventually becomes discrepant for higher interactions. We also include an exponential fit of the first dataset, with J/U = 0.62, taken only above ω < 8J/ and from which we extract an effective local energy scale of J eff = 3.0(3) J. Note that in this case the presence of interactions leads already to a deviation from the simple exponential trend, even for the weakest interaction explored here.
For stronger interactions, the heating rate becomes nonmonotonic, analogously to the 2D results and the 1D numerics, though in this case the associated features are less sharp. This could be explained by the inhomogeneity present in the 1D system, due to much stronger confinement from the transverse lattice. We also note that, while the numerical simulations capture most features of the 1D experimental data, quantitative discrepancies between the timescales are to be expected for various reasons, chiefly the different boundary conditions (hard wall vs harmonic confinement), but also the drive amplitudes (infinitesimal vs finite) and protocols (the experiment naturally includes a weak modulation of U which is not considered in the simulations).

VII. CONCLUSION
We have studied the nature of heating in a system of driven ultracold bosonic atoms in an optical lattice, and found strong evidence that the thermalization time diverges exponentially in the drive frequency -a central prediction in the theory of Floquet prethermalization. When increasing the frequency only by a factor of 2 or 3, we were able to suppress the heating rate by as much as two orders of magnitude. This possibility of driving isolated quantum systems while avoiding heating for exponentially long times is a key step towards the engineering of new forms of matter existing only out of equilibrium.
Furthermore, our results add physical insight to the intuitive picture of Floquet prethermalization -that a quantum system driven at high enough frequency heats slowly because the absorption of a quantum of energy ω from the drive requires many coordinated rearrangements of its local degrees of freedom. Our experimental results not only confirm this picture, but also shed light on the understanding of what those rearrangements look like in a real system. Interestingly, they point to different scenarios for the superfluid and Mott-insulating sides of the Bose-Hubbard phase diagram.
In the future, it will be interesting to explore Floquet heating in our system in the limit of hard-core bosons, where a different prethermalization mechanism, based on weak integrability breaking, may be realized [45]. Another direction is the addition of disorder, where it would be interesting to microscopically characterize the failure of the MBL phase to heat to infinite temperature. Other possible directions include implementing more exotic drives, such as quasiperiodic ones [48,49]; using strong drives to probe the heating rates with our technique beyond the linear response regime; and exploring the dependence on the initial temperature. Most interestingly, our work paves the way for future realizations of novel prethermal Floquet phases of matter with no equilibrium analogs.
Note added.
During the preparation of this manuscript, we became aware of related work providing evidence for the observation of exponential-in-frequency thermalization times in dipolar spin chains [72].
Correspondence address: antonio.rubio@mpq.mpg.de. After the preparation of a two-dimensional degenerate gas of Rubidium-87 atoms, we load the atoms into a square optical lattice, generated by two retroreflected laser beams in the atomic plane, with lattice spacing a lat = 532 nm. We ramp up the in-plane optical lattices with an s-shaped ramp of 75 ms up to a lattice depth V 0 expressed in units of the energy recoil E r = h 2 /8ma 2 lat , where m is the atomic mass. Ramps of the same duration and form are used to ramp to the transition point and to the atomic limit after the driving dynamics, where the energy content of the gas is measured.
In addition to the lattice potential, the atoms experience an overall harmonic trapping potential, given by where ω x and ω y are the harmonic trap frequencies in the plane, which are slightly different for each lattice depth. In the 2D case, they are in the range 2π × 45 Hz < ω x ,ω y < 2π × 55 Hz, while in the 1D case ω x 2π × 70 Hz is roughly constant.

Bose-Hubbard parameters
The values of J and U used in the main text are obtained from a numeric calculation of the band structure, and are based on the calibrated lattice depths V 0 , which are estimated to have an uncertainty of roughly 2%. We show here all the calculated parameters for the relevant depths corresponding to Fig. 3 and Fig. 4. In the 2D case, both in-plane lattices had the same lattice depth and the lattice modulation was A = 0.05 (see Tab. I).
In the 1D case, the lattice along the y axis was fixed to 20E r and the lattice along the x axis was tuned to V 0,x and driven with amplitude A x = 0.1 (see Tab. II). We also plot the modulation of the tunneling strength δJ, obtained as δJ = (J V0−A − J V0+A )/2.

Higher bands
As we mention in the main text, we need to keep the drive at small enough frequencies and low enough amplitudes in order to not populate higher bands. Naively, this should only require staying below the gap to the second excited band, E g,2 = E 2 (q = 0) − E 0 (q = 0), since due to symmetry reasons there is no coupling to the first excited band with gap E g,1 = E 1 (q = π/a) − E 0 (q = 0). However, multi-photon resonances can lead to interband transfer even for frequencies well bellow the gap energies, such that in practice we need to identify the regimes at which interband heating starts to take place and stay below those. In Table III, we plot both E g,1 and E g,2 , also obtained from band-structure numerics, for five different lattice depths within the explored regimes. All the frequencies used to drive the lattice depth in the experiment are well below both E g,1 and E g,2 /3. The experimental heating rates shown in Fig. 3 and Fig. 4 are extracted from the temperature dynamics V0(Er) Eg,1/h(kHz) Eg,1(J) Eg,2/h(kHz) Eg,2(J) 3 3.9 17 9.  within the linear heating regime. The heating rate per Floquet cycle can be expressed as φ(ω) = k B dT /dt × 2π/ω, where k B is the Boltzmann constant, T the temperature, and 2π/ω is the drive cycle period. We verify that we probe the system within linear response by plotting in Fig. 5 sample heating dynamics in 2D at V 0 = 6E r for four different driving frequencies with a relative driving amplitude A = 0.05. The initial temperature of the cloud is typically around 0.1 U 0 /k B , where U 0 /h = 660 Hz is the interaction strength at the atomic limit. In units of the tunneling strength J at 6E r the temperature is roughly k B T = 1.4J. We can see that the increase in temperature is consistent with being in the linear response regime.

Scaling from Fermi's Golden Rule
We have also explored the heating rates for different amplitudes of the drive in 2D, since in the linear regime we expect a scaling of the heating rate predicted by Fermi's Golden Rule, i.e. proportional to the drive amplitude squared. From fitting a power law with the expression φ(A) = cA α we obtain c = 4.0(4) and α = 2.11(4). Fermi-golden-rule scaling. Heating rate as a function of the drive amplitude on a log-log scale. The data was taken at V0 = 8 Er with driving frequency ω = 14.5 J/ . The solid line is a power-law fit. The errorbars denote the s.e.m.

Appendix C: Prethermalization bound for bosons
The bound on the spectral function (linear-response heating rate) proven in Ref. [35], for appropriate constants C and J eff , is strictly speaking only valid for for systems with a bounded local Hilbert space such as spins, fermions or hardcore bosons. The Bose-Hubbard system we study in this work instead allows for unbounded occupation of each site, which was argued in Ref. [35] to generically relax the bound (C1) from exponential to stretched-exponential, Here we discuss the relevance of this relaxed bound to our experimental observations. As the physical reason for relaxing the bound to Eq. (C2) has to do with unbounded energy density, the most natural place to look for violations of the original exponential bound is the Mott-insulating phase. There, Φ(ω) is strongly non-monotonic and exhibits peaks near ω = mU , m ∈ N, as shown in Fig. 2(b). These can be understood via perturbation theory from the Mottinsulator state in the atomic limit J = 0: each ω = mU peak appears at some order p m in perturbation theory, giving a leading contribution ∼ (J/U ) 2pm , or in a form reminiscent of the prethermal bounds.
Generally, the optimal process to absorb the most energy in the fewest moves consists of depleting a whole contiguous region on the lattice and gathering all its particles on a central site. Taking a sphere of radius R in the hypercubic lattice in d dimensions, this process gives an energy absorption of ω ∼ R 2d in p ∼ R d+1 "moves", yielding p m ∼ (m 1/2d ) d+1 . This scaling gives in general a stretched exponential bound like Eq. (C2) with power In one dimension, we recover α = 1, i.e. the exponential bound. On the other hand, in higher dimension this construction gives a series of peaks that violate the exponential bound. However, such peaks occur at very high frequency and are sparsely distributed in the spectrum -in 2D, the lowest such peak is the quintuplon at ω = 10 U (gathering all four nearest neighbors), followed by the "13-uplon" at ω = 78 U (gathering the eight next-nearest neighbors), etc. Such high frequencies are unlikely to be achievable in experiment without exciting other degrees of freedom, and are thus irrelevant in practice.
In the superfluid phase, single-site high-occupancy excitations play a less prominent role than they do in the Mott phase, which means the above physics is even less likely to be relevant to observations, justifying the use of the exponential bound of Eq. (C1) somewhat beyond its domain of mathematical rigor.
We conclude by noting that Ref. [73] found that the quantity J eff in the exponential bound Eq. (C1) is fundamentally connected to many important aspects of quantum dynamics, including operator growth and chaos. This raises interesting questions about such issues in bosonic systems, particularly whether similar bounds on chaos and complexity also hold for a generalized bound like Eq. (C2). subspace method to simulate the dynamics of a state |ψ(N cyc ) at stroboscopic times t = N cyc T (T = 2π/ω is the drive period) and track the value of the "energy" ψ(N cyc )|Ĥ 0 |ψ(N cyc ) ≡ E Ncyc during the evolution to define a thermalization time N th cyc . This method can probe slightly larger sizes than full diagonalization (up to L = 13 sites, Hilbert space dimension in the millions) and, most importantly, is not restricted to drive amplitudes g 1, at the expense of a limited dynamic range in N th cyc . We define a normalized energy density where E ∞ ∝ Tr(Ĥ 0 ) is the infinite-temperature value of the energy. During the dynamics, this obeys 0 ≤ ε(N cyc ) ≤ 1. We define the heating time N th cyc as the lowest N cyc such that ε(N cyc ) is greater than some predefined threshold ε (we use 0.1, though other choices give similar results). We also keep track of the "density One-dimensional chain of L = 12 sites at unit filling in the superfluid phase, driven with amplitude g = 0.5 and subsequently ramped to the atomic limit. Density of holes ρ h plotted as a function of the normalized energy density ε for different values of the frequency of the drive. The frequencies range from ω = 2 J/ (gray) to ω = 12 J/ (dark green).
of holes" , whereρ h,i is a projector onto even occupation of site i , which mimics the fluorescence imaging technique used in the experiment. We choose the initial state |ψ(0) as the ground state ofĤ 0 (obtained via the Lanczos method). We then approximate each Floquet cycle by a sequence of s constant Hamiltonians, {Ĥ(t = T k/s) : k = 0, . . . s − 1}, and time-evolve the state vector for time T /s with each of these Hamiltonians using the Krylov subspace method. In practice, we use s = 32 steps; further increasing s changes the results negligibly.
In addition, to better imitate the experimental procedure, before measuring the observables ε and ρ h we ramp the system into an atomic-limit Mott insulating state, i.e. we arrest the drive and slowly take J → 0 as J(N cyc T + t) = (1 − t/τ ) J(N cyc T ) for 0 ≤ t ≤ τ , where τ is a long time scale (we use τ = 100 /U ). In practice, this is again accomplished by time-evolving with piecewise constant Hamiltonians, keeping the same time step used during the drive. A copy of the wavefunction at t = N cyc T (before the ramp) is stored so that the time evolution can resume after the measurement is taken.
Results for a fairly strong drive amplitude g = 0.5 (see Fig. 7) generally agree with the linear response picture -N th cyc (ω) is approximately flat for ω Ω 2qp , then starts increasing exponentially. The time traces of the density of holes ρ h , shown in Fig. 7(a), are very similar to the experimental data in Fig. 1(c). Comparison with the energy density traces in Fig. 7(b) also confirms that ρ h is indeed a good proxy for the energy density ε. We further confirm the relation between these two quantities by plotting the density of holes vs the normalized energy density (see Fig. 8).
As the method simulates time evolution directly, its cost scales as O(N th cyc ), hence exponentially in ω; this limits us to ω 12J, and in particular prevents us from testing the presence of an additional threshold feature near 2 Ω 2qp as seen in linear response ( Fig. 2(a)). It also makes the method generally less suited to the Mottinsulating phase, where N th cyc has non-monotonic oscillations by many orders of magnitude.