Parametric Heating in a 2D Periodically Driven Bosonic System: Beyond the Weakly Interacting Regime

We experimentally investigate the effects of parametric instabilities on the short-time heating process of periodically driven bosons in 2D optical lattices with a continuous transverse (tube) degree of freedom. We analyze three types of periodic drives: (i) linear along the x -lattice direction only, (ii) linear along the lattice diagonal, and (iii) circular in the lattice plane. In all cases, we demonstrate that the Bose-Einstein condensate (BEC) decay is dominated by the emergence of unstable Bogoliubov modes, rather than scattering in higher Floquet bands, in agreement with recent theoretical predictions. The observed BEC depletion rates are much higher when shaking along both the x and y directions, as opposed to only x or only y . We also report an explosion of the decay rates at large drive amplitudes and suggest a phenomenological description beyond the Bogoliubov theory. In this strongly coupled regime, circular drives heat faster than diagonal drives, which illustrates the nontrivial dependence of the heating on the choice of drive.


I. INTRODUCTION
An area of increasing interest in ultracold atoms concerns the engineering of novel states of matter using highly controllable optical lattices [1]. In this context, a promising approach relies on applying time-periodic modulation to the system, in view of designing an effective time-independent Hamiltonian featuring the desired properties [2][3][4]. This Floquet engineering has emerged as a promising and conceptually straightforward way to expand the quantum simulation toolbox, enabling appealing features such as suppressed [5,6] or laser-assisted [7] tunneling in optical lattices, enhanced magnetic correlations [8], state-dependent lattices [9], and subwavelength optical lattices [10], as well as synthetic dimensions [11,12], synthetic gauge fields [13,14], and topological band structures [15].
Despite these promising applications, progress in Floquet engineering for interacting systems has been hindered by heating due to uncontrolled energy absorption. This heating, which stems from a rich interplay between the periodic drive and interparticle interactions, is a particularly challenging problem in interacting systems, where it is known to occur due to the proliferation of resonances between many-body Floquet states, not captured by the inverse-frequency expansion [4,16]. This problem constrains the applicability of Floquet engineering to regimes where heating is slower than the engineered dynamics [9,[17][18][19][20][21]. A deeper understanding of the underlying processes is essential to determine stable regions of the (large) parameter space, where the system is amenable to Floquet engineering. Additionally, interaction-mediated heating is itself an interesting nontrivial quantum many-body process. Energy absorption and entanglement production in periodically driven systems have recently been the focus of theoretical studies [16,[22][23][24][25][26][27][28][29] and experimental investigations [5,9,20,21]. It was predicted that, whenever the drive frequency is larger than all single-particle energy scales of the problem, heating succumbs to a stable long-lived prethermal steady state, before it can occur at exponentially long times [25,[30][31][32][33][34][35]. However, it is unclear whether this physics is accessible in interacting bosonic experiments.
A perturbative approach to understanding drive-induced heating is to analyze the underlying two-body scattering processes using Fermi's golden rule (FGR) [21,23,[36][37][38]. In the weakly interacting limit, interactions provide a small coupling between noninteracting Floquet states. However, Floquet states cannot be treated as noninteracting when the Floquet-modified excitation spectrum is itself unstable [22,[24][25][26]. These instabilities indicate that heating can occur on a shorter timescale than expected from the scattering theory alone.
For Bose-Einstein condensates (BECs) in optical lattices, increased decay rates arise due to the emergence of unstable collective modes. The resulting parametric heating can be described using a Floquet-Bogoliubov-de Gennes (FBdG) approach [24], and the short-time dynamics is dominated by an exponential growth of the unstable excited modes in the BEC. The depletion time of the condensate fraction provides an experimental window to observe this and related effects. Qualitatively different behavior is expected between scattering and parametric instability rates, most notably, different power laws as a function of the interaction strength, tunneling rate, and drive amplitude.
We experimentally explore these predictions in a 2D lattice subject to 1D and 2D periodic drives, by measuring the decay of the BEC condensed fraction. We provide strong experimental evidence that parametric instabilities dominate the short-time dynamics over FGR-type scattering processes, which are responsible for long-time thermalization [21]. Our experiment reveals effects beyond Floquet-Bogoliubov predictions and points out limitations in the applicability of the FBdG theory.
The experiments are performed on a BEC of 87 Rb atoms loaded into a square 2D optical lattice [39,40] with principal axes along x and y, formed by two pairs of counterpropagating laser beams with a wavelength of λ ¼ 814 nm. The total atom number is N ≃ 10 5 (AE20% systematic uncertainty). Two piezoactuated mirrors [41] sinusoidally translate the lattice along x and y with any desired amplitude, relative phase, and angular frequency: rðtÞ ¼ fΔx sin ðωtÞ; Δy sin ðωt þ ϕÞg. We consider the effect of three drive trajectories on the decay rate: translation along x only (Δy ¼ 0), diagonal translation along x and y (Δy ¼ Δx and ϕ ¼ 0), and circular translation (Δy ¼ Δx and ϕ ¼ π=2). Therefore, the driving is 1D (x only) or 2D (diagonal or circular), in a 2D system (2D array of tubes), as shown in Fig. 1. We express the amplitude Δx in terms of the drive-induced maximum effective energy offset between neighboring lattice sites in the comoving frame, K 0 ¼ ΔE=ℏω, where ΔE ¼ mω 2 aΔx, a is the lattice spacing, and m is the 87 Rb mass. The physical displacement is Δx ¼ ℏK 0 =aωm. The lattice depth V 0 is held constant during shaking and is measured in units of lattice recoil energy E R ¼ h 2 =mλ 2 . The lattice tunneling energy ℏJ and the interaction strength g are controlled via V 0 . The value of JðV 0 Þ and gðV 0 Þ are calculated from the band structure, peak atom density, and scattering length (see Appendix B). This calculation results in the following periodically driven Bose-Hubbard Hamiltonian: FIG. 1. Schematic of the lattice driving. (Top) Lattice translation is performed along the x direction (green), along the diagonal (blue), and in circles (red). The normalized drive amplitude projected along the x axis, K 0 , is used to characterize the drive strength for all trajectories. (Bottom) Example of a periodic drive for circular driving.
whereâ ð †Þ ij ðzÞ is the annihilation (creation) operator at lattice site ði; jÞ and transverse position z and κ ¼ 0 for x-only and κ ¼ 1 for 2D drives. The interaction U is defined such as U=V P i;j R z hâ † i;j ðzÞâ i;j ðzÞi ¼ g with V the volume of the system.
In order to avoid micromotion effects during a drive period [10], the experiments are performed at integer multiples of the period T ¼ 2π=ω (see Appendix A). The drive amplitude is ramped up smoothly [10,42,43] in a fixed time (minimum 2 ms) corresponding to an integer number of periods (Fig. 1). The shaking is then held at a constant amplitude for a time τ. Finally, the amplitude is ramped down to zero in a few periods. We check the effect of adiabaticity of this procedure (or lack thereof) for the three drive trajectories and describe it in the Appendix A. Once the lattice is at rest, we turn it off in 300 μs to determine the atomic distribution. We use absorption imaging after the time of flight to measure the condensate fraction as a function of τ.
For most conditions, the condensate decay agrees with an exponential decay [NðtÞ ¼ Nð0Þe −Γ cf t ], whose rate Γ cf we extract from a least-square fit (see Appendix C). We measure Γ cf for the three drives at different values of ω, K 0 , and V 0 . FBdG predicts an undamped parametric instability, characterized by exponential growth of unstable modes, i.e., accelerated condensate loss. This behavior is inconsistent with the measured exponential decay of the BEC. Hence, the undamped FBdG regime does not last long compared to the typical BEC lifetime for our parameters, and interactions between the excited unstable modes and the BEC play a significant role in the observed heating process. Nonetheless, as we discuss below, the magnitude and scaling of Γ cf are well captured by a FBdG description.
Since the Floquet-renormalized hopping is J eff ¼ JJ 0 ðK 0 Þ [J ν ðK 0 Þ is the νth order Bessel function], J eff <0 for K 0 > 2.4 and the lowest Floquet band is inverted [5]; the BEC then becomes dynamically unstable at q ¼ ð0; 0Þ [44], but a stable equilibrium occurs at the band edge, and a sudden change in the stability point can trigger a dynamical transition [45] between the two equilibrium configurations. In the band-inverted regime, the stable quasimomenta are q ¼ ðAEπ; 0Þ for a 1D drive along x and q ¼ ðAEπ; AEπÞ for a 2D drive along x and y, where the components of the crystal momentum q are measured in units of the inverse lattice spacing a −1 . FBdG assumes an initial macroscopic occupation of these modes [24]. Unless stated otherwise, for data taken at K 0 > 2.4, we first accelerate the BEC to the appropriate stable point while simultaneously turning on the Floquet drive (see Appendix B).

II. RESULTS
A. Lattice depth scans: Γ cf ðV 0 Þ A major difference between FGR and the FBdG theory is the scaling of the decay rate with the hopping J and interaction strength g. Whereas FGR predicts Γ cf ∝ ðgJÞ 2 , the parametric instability rate is expected to be linear (Γ cf ∝ gJ) [24]. Figure 2 shows the condensed fraction decay rate Γ cf measured at different V 0 (between 7.3E R and 17.0E R ) and ω, plotted as a function of gJ=ω, for the 2Ddiagonal drive and 1D x-only drive. The solid lines show the FBdG theory, and the dashed lines are linear fits to the data.
For the x-only drive, the magnitude and slope extracted from the data are well described by the FBdG theory, and while some deviations appear for the 2D drive, both results are clearly inconsistent with a quadratic dependence.
B. Amplitude scans: Γ cf ðK 0 Þ Figure 3 shows the decay rate as a function of the drive amplitude K 0 , at lattice depth V 0 ¼ 11E R (which gives J ¼ 2π × 50 Hz, g ¼ 2π × 700 Hz, and a gap ΔE v =h ¼ 21 kHz to the next vibrational level) and a drive frequency FIG. 2. Γ cf ðJg=ωÞ: Measured decay rates for K 0 ¼ 2.1, ω ¼ 2π × 4 kHz (filled circles), and ω ¼ 2π × 2.5 kHz (empty circles) and various lattice depths, compared to the FBdG theory. The horizontal axis is Jg=ð2 × ωÞ, since the FBdG theory predicts the mode growth rate to be linear in J, g, and 1=ω, as derived in the Appendix E. The dashed lines are linear fits to the data. ω ¼ 2π × 2.5 kHz. The same data are shown on a full range (top) and enlarged (bottom).
Consider first the x-only drive. The instability growth rate predicted from the FBdG theory [16] Γ mum ¼ 8JJ 2 ðK 0 Þg=ω, which we detail in the Appendix E, agrees with the measured decay rates (Fig. 3). The appearance of the Bessel function J 2 ðK 0 Þ in the expression for the rate can be traced back to the parametric resonance condition, requiring the excitation spectrum to match the drive frequency; cf. the Appendix. In contrast, the FGR scattering approach prediction is too low by a factor of about 30, and the predicted scaling, ∝ jJ 2 ðK 0 Þj 2 to leading order, does not describe the data as well. The agreement with the FBdG theory, despite evidence for effects beyond simple undamped parametric instability, is consistent over a range of parameter space in K 0 , ω, and V 0 . As expected, the decay dramatically increases when J eff < 0 (K 0 > 2.4) for q ¼ ð0; 0Þ, while it is partially stabilized when accelerating the BEC to q ¼ ðπ; 0Þ. We note that significant heating occurs during the drive turn-on and acceleration phase for the q ¼ ðπ; 0Þ data, resulting in partial BEC losses.
For the two 2D drives (circular and diagonal), we observe decay rates that follow roughly the same functional form as the 1D drive but about 3× larger. Additionally, the sudden increase in the decay rate Γ cf for the 2D drives consistently occurs at a critical amplitude K c 0 below 2.4. At 11E R , K c 0 ≃ 2.15 (Fig. 3). Above K c 0 , both 2D rates increase dramatically beyond any prediction, and the circular rate increases faster than the diagonal rate. This drive dependence and drastic rate increase for K c 0 < K 0 < 2.4 suggest effects beyond the FBdG theory, distinct from the simple parametric instability, and are discussed at the end of this paper. For K 0 ≳ 2.4, the rates are essentially unmeasurable since Γ cf > ω. As with the 1D drive, accelerating the BEC to q ¼ ðπ; πÞ partially stabilizes the decay, just enough to be measurable.

C. Frequency scans: Γ cf ðωÞ
The FBdG theory predicts distinct behavior at low and high drive frequencies. In the low-frequency regime, the momenta of the maximally unstable mode q mum evolve as one increases the drive frequency, until it saturates at the Bogoliubov band edge at q mum ¼ fðπ; 0Þg, fðπ; 0Þ; ð0; πÞg, and fðπ; πÞg for the x-only, circular, and diagonal drive, respectively. The periodic drive results in interference which depends on the relative phase of the two components of the drive, leading to different most unstable modes for the three shaking geometries. The saturation frequency ω c ¼ E Bog eff ðq mum Þ marks the onset of the high-frequency regime [24]. Note that, unlike in 1D lattices [24], the energy of the maximally unstable mode can be lower than the full effective bandwidth (see Appendices B and E). The rate is predicted to increase quasilinearly for ω ≤ ω c , while Γ mum ∝ ω −1 for ω ≥ ω c [21,24], resulting in a cusp in the rate at ω c . A detailed derivation of the expressions for the most unstable modes and the instability rates are given in the Appendix. Figure 4 shows the experimental values for Γ cf ðωÞ compared with the FBdG theory. Since there is no observed rate explosion for the x-only drive, we use the cusp in Γ x cf ðωÞ to calibrate our experimental value for g, which agrees to within 20% with an estimate calculated from the lattice parameters (see Appendix B). Using this value of g, the prediction for the diagonal drive ω diag c matches the experiment. While the FBdG theory predicts the same ω c FIG. 3. Γ cf ðK 0 Þ: Comparison between circular (red), diagonal (blue), and x-only (green) drives for V 0 ¼ 11E R and ω ¼ 2π × 2.5 kHz. The rates are in units of J ¼ 2π × 50 Hz, shown full scale (top) and enlarged (bottom). The Floquet-Bogoliubov-de Gennes theory [24] (see Appendix E) Γ mum is shown for each drive as solid lines, and the 1D FGR-based scattering theory is shown as the green dashed line (enlarged only). Filled circles indicate data taken at q ¼ ð0; 0Þ, while open circles indicate data taken at q ¼ ½π; πð0Þ (see the text) to keep the BEC in the stable region of the band (illustrated with the bottom plot cartoon). A dramatic increase in the rate occurs at K c 0 ≳ 2.4 for the circular and diagonal drives, highlighted by the light red region.
for circular and x-only drives, the measured cusp for the circular drive lies between the cusps of the two linear drives.
The observed behavior qualitatively fits the FBdG theory, with rates generally higher than predicted for ω ∼ ω c . For the x-only drive, the measured rates are slightly above the prediction below 2π × 1.5 kHz, and the agreement is excellent above 2π × 1.5 kHz, as observed with Γ cf ðK 0 Þ at 2π × 2.5 kHz (Fig. 3). The 2D rates show a larger discrepancy at low frequencies and a decent quantitative agreement for ω ≳ 2π × 2 kHz. This result is related to the rate explosion appearing for K 0 > K c 0 . As we discuss below, the observed value of K c 0 increases with the frequency. This increase implies a similar rate explosion should happen when decreasing ω at fixed K 0 . This explosion is especially visible with the diagonal drive ( Fig. 4): For ω < 2π × 1 kHz, the data abruptly depart from the prediction. This increased rate at low frequencies for 2D drives is likely responsible for the discrepancy between the experiment and theory. The presence of the cusps in the rate explosion region is still expected, since, for ω low enough, some modes are energetically inaccessible, and the limit Γ cf ðω → 0Þ → 0 must be fulfilled.

D. Rate explosion
Beyond a critical amplitude K c 0 , we observe a sudden increase of the 2D-driven decay rate (Fig. 3). The dependence of K c 0 on the frequency for circular driving Fig. 5. We observe that K c 0 → 2.4 as ω → ∞, suggesting that the giant heating arises from a finite-frequency effect. Because the effect appears for both diagonal and circular drives, whereas a perturbative argument using the inverse-frequency expansion shows that there are finite 1=ω corrections only for the circular drive, we surmise that the effect is likely nonperturbative. Assuming the anomalously strong heating results from an interplay of correlated physics beyond the Bogoliubov regime and recalling that resonance effects beyond the infinite-frequency approximation lead to a 1=ω dependence in the instability rates [2][3][4], we make the following scaling argument. Viewing the system as an effective Bose-Hubbard model, the strongly correlated regime is reached for g=J eff ≳ 1. On the other hand, we note that corrections to the infinite-frequency Floquet Hamiltonian scale as J=ω. We make the phenomenological observation that the dimensionless ratio ðg=J eff ÞðJ=ωÞ should be relevant to a combination of beyond-mean-field and finite-frequency effects. When J eff is low in all lattice directions, ðg=J eff ÞðJ=ωÞ is large, and, therefore, these effects should be large. The simple scaling relation g=ωJ 0 ðK 0 Þ ¼ 1 gives K c 0 ðωÞ ¼ J −1 0 ðg=ωÞ and is shown as the red line in Fig. 5 with the experimental data. The agreement is surprisingly good for such a simple argument. The quantum many-body nature of the rate explosion calls for more extensive study, that promises new insights into periodically driven strongly correlated quantum lattice systems.
We presented a detailed investigation of heating for interacting bosons in a periodically driven 2D lattice. The observed decay rates are substantially larger than expected from a scattering theory based on Fermi's golden rule [21] and scale as expected for interaction-driven parametric instabilities [24]. The observed exponential decay of the condensed fraction suggests that interactions between these excited modes and the BEC, not captured by the FBdG theory, play an important role in the dynamics. Nonetheless, the linear scaling of the condensate loss rate with gJ=ω is indicative of direct, interaction-induced instabilities. Importantly, these instabilities arise from collective modes and involve coherent processes, unlike scattering in a purely FGR approach. In addition, for 2D driving, there exist regions where the heating is even larger than predicted by FBdG, which is not explained by current theories. Altogether, our observations provide important insight into a major heating mechanism in bosonic systems subject to a position-modulation drive, a valuable knowledge for future many-body Floquet engineering schemes.
We note that complementary signatures of parametric instabilities have been recently investigated with bosonic atoms in periodically driven 1D optical lattices [46].  [47,48] to perform the numerical simulations. The authors are pleased to acknowledge that the computational work reported on in this paper was performed on the Shared Computing Cluster which is administered by Boston University's Research Computing Services.

APPENDIX A: THERMALIZATION OF EXCITED ATOMS
The condensate fraction can decay either by direct Floquet-driven loss or by heating due to relaxation of energetic excitations. The latter mechanism occurs over a thermalization timescale and can be probed by observing relaxation of out-of-equilibrium states in an undriven, static lattice. Special attention to the drive turn-off is required to avoid unwanted excitation due to micromotion during a drive period. Abruptly turning off the drive induces a kick large enough to create a significant out-of-equilibrium population by interband excitation.
To determine the relaxation timescale and test for this additional condensate loss mechanism, we measure the evolution of the condensed fraction when holding the atomic cloud in a static lattice, immediately after an abrupt stop of the drive. While the condensed fraction is initially unchanged, upon letting the static system evolve for a time t hold , we observe a subsequent decrease of the condensed fraction as excited atoms thermalize with the rest of the sample. Figure 6 shows an example of thermalization for a diagonal drive at ω ¼ 2π × 2 kHz. When there is sufficient initial excitation (blue data in Fig. 6), we observe a characteristic thermalization time of the order of 2 ms that does not depend on the initial nonequilibrium population.
We observe that the condensed fraction increases slightly with the hold time in the cases where minimal atom excitation occurs (red and green data in Fig. 6). We attribute this increase to rethermalization along the tube axis: Entropy added in the degree of freedom (d.o.f.) associated with the lattice direction can be redistributed along the z (tube) axis, which is visible as a slightly decreased temperature in the x-y plane.
The amplitude of the condensed fraction decrease is indicative of the energy of the initial nonequilibrium population. We can measure this amplitude by comparing the condensed fraction at t hold ¼ 0 and t hold ¼ 6 ms, a time sufficient for the thermalization to have occurred. Figure 7 FIG. 6. Thermalization of the nonequilibrium population. The lattice is shaken for a few periods and then held static at t hold ¼ 0. There exists at t hold ¼ 0 a nonequilibrium population that thermalizes with the BEC in a few milliseconds, resulting in a clear decrease of the condensed fraction. The amplitude of this effect depends on the initial size of the nonequilibrium population. It is maximum if atoms are kicked to higher bands by an abrupt stop of the drive (blue). If no such population is created by the stopping of the drive, for example, by choosing an end phase that minimizes the kick (red) or by rapidly ramping down the drive amplitude (green), no condensed fraction decay is observed.
shows the remaining condensate fraction after t hold ¼ 6 ms as a function of the phase at which the drive is stopped, for an abrupt stop (blue) and for a smooth but rapid (1 period) turn-off (green). When the stop is abrupt, the fraction of atoms excited depends on the end phase: Abruptly immobilizing the lattice can induce an effective force which depends on the lattice velocity immediately prior to immobilization. Therefore, the excited fraction is maximized when _ qðtÞ is maximally discontinuous (stop phases of 0½π rad) and minimized when smooth (stop phases of ðπ=2Þ½π rad). In all our decay rate data where we smoothly turn off the drive (in at least one period), no turn-offinduced heating is observed.
Abruptly stopping the drive is not the only potential source of nonequilibrium population. The unstable Bogoliubov modes studied in the main text could themselves produce a population that thermalizes and that could potentially modify the measured decay rates. This possibility can be tested with the same stop-and-hold measurement, if turn-off-induced populations are avoided. As visible in Fig. 7, this possibility is realized for stop phases of ðπ=2Þ½π rad or when smoothly turning off the drive. For our data, avoiding a turn-off-induced population when stopping the drive results in no visible decay, which is shown as the green data in Fig. 6. The result is identical for the whole region of parameter space studied here. We deduce that the energy carried by the nonequilibrium population created by the unstable modes is not enough to significantly impact the measured rates.

APPENDIX B: EXPERIMENTAL CONSIDERATIONS
Accelerating the BEC to the band edge.-In order to accelerate our BEC from q ¼ ð0; 0Þ to a desired q ¼ ðq x ; q y Þ, we apply a constant force F ¼ _ q for a fixed time in the lattice plane. The force is generated with a constant magnetic gradient, acting on the BEC in the jF ¼ 1; m F ¼ −1i ground state. Bias coils in the three spatial directions control the gradient direction.
The BEC becomes unstable for quasimomenta about halfway to the band edge, which is the well-known static instability [44]. On the other hand, as mentioned in the main text, ramping up the drive amplitude beyond K 0 ¼ 2.4 reverses the band smoothly as J eff becomes negative: q ∼ ð0; 0Þ becomes unstable, and the band edges (or corners, depending on the drive trajectory) become stable. In order keep the BEC in a stable region (in an effectively static dynamical stability sense) at any given time, we synchronize the BEC acceleration with the ramping on of the drive, such that the BEC crosses the static instability point q ¼ 0.6π=a when K 0 ¼ 2.4.
We accelerate the BEC to q ¼ ðπ; πÞ for the two 2D drives (diagonal and circle) and to q ¼ ðπ; 0Þ for the x-only 1D drive, since these become stable whenever K 0 > 2.4.
Calibration of translation amplitude.-The piezoactuated mirrors are each roughly calibrated offline using an interferometric technique [41]. The final calibration is performed on the atoms by measuring the tunnelingdependent magnetization decay rate of 2D staggered spin magnetization [49] as a function of drive strength K 0 . The drive strength at which the tunneling rate J eff vanishes is determined by the condition K 0 ¼ 2.4.
Calibration of tight-binding parameters.-The lattice depth is calibrated via Raman-Nath diffraction. In the tightbinding limit, the tunneling rate is derived through the modeled 1D dispersion as The on-site interaction g is calibrated from the x-only drive cusp (see Fig. 4). The point at which the rates go from increasing to decreasing, ω diag c , is given by Knowing J, the experimental value of ω x c offers a calibration for g. For V 0 ¼ 11E R , J ¼ 50 Hz and the measured value of ω x c ¼ 444 Hz gives g ¼ 700 Hz. As an additional check, this value of g is then used to predict the diagonal drive cusp, The observed value of approximately 650 Hz is in good agreement with this prediction.
The interaction strength g depends upon the atom number, the dipole trap, and the lattice parameters. To confirm that the experimental calibration matches these FIG. 7. Amplitude of the thermalization process. The difference between the condensed fraction at t ¼ t end (t hold ¼ 0) and at t ¼ t end þ 6 ms (t hold ¼ 6 ms), as a function of the end drive phase, for a 2π × 2 kHz diagonal (2D) drive. For an abrupt stop (blue), a large population can be transferred out of the BEC for some end phases, which results in a large thermalization event and a clear drop of the condensed fraction. The effect is minimal when stopping the drive such that qðtÞ is a smooth function. If, however, the drive is stopped by ramping down the amplitude (green), no condensed fraction decay can be observed at any phase. known experimental parameters, we also estimate g through tight binding and Thomas-Fermi assumptions: In the lattice plane, the wave function ψðxÞ is taken to be well approximated by a Mathieu function, while we use a Thomas-Fermi profile in the tube direction. The interaction energy g ∝ ∬ a=2 −a=2 jψðrÞj 4 d 2 r is then calculated from the known experimental parameters, including the density profile due to the dipole trap (frequencies fω x ; ω y ; ω z g ¼ f11; 45; 120g Hz). For V 0 ¼ 11E R , we find g ¼ 850 Hz, similar to the calibrated value of 700 Hz. Note that the systematic 20% uncertainty in the atom number can easily explain the small offset between the estimation and calibration.
Bandwidths.-It is important to note that ω c is, in general, different from the effective bandwidth B. For a 2D (diagonal or circular) drive in our 2D lattice, and, with a 1D drive in the 2D lattice, which is, in general, different from ω c , as observed in the main text: Only in the case of the diagonal drive do we find that the maximally unstable mode had the maximum ground band energy, and therefore ω c ¼ B.
Background rates.-All theoretical plots take into account the background decay rate, predominantly due to lattice photon scattering. We experimentally determine this rate by setting K 0 ¼ 0 and measuring the resulting rate with the same procedure as in the main text. This constant rate y 0 ∼ 1 s −1 for V 0 ¼ 11E R is then added to the FBdG formula for comparison with the experimental data. Since multiphoton resonances to higher bands [50] could complicate the interpretation of the measured decay rates, we perform heating measurements up to drive frequencies of ω ¼ 2π × 21 kHz to identify excited band resonances. Population transfer to higher bands is directly visible on the band-mapping images. The lowest frequency at which resonant higher band excitation is observed is ω ¼ 2π × 6.25 kHz for V 0 ¼ 11E R , and we therefore limit our heating measurements to below 2π × 4 kHz to avoid these effects.
Heating in the absence of a Floquet drive.-To rule out heating mechanisms that do not depend on the drive but are instead related to the change in J eff , it would be useful to compare the heating observed at a given K 0 to the heating observed in a static lattice with a depth chosen to have an equivalent J ¼ J eff (for J eff > 0). Here, such a direct comparison is experimentally challenging under the same conditions as the Floquet experiment, due to the change in the trap confinement when increasing the lattice depth. In addition to a potential change in gravitational sag (which we minimize by aligning the optical lattice beams directly on the dipole-trapped BEC), the change in trap frequency when changing the lattice depth can excite breathing motion along the nonlattice tube direction, as well as cause a redistribution of atoms within the trap transverse to the lattice direction. These trap-changing effects do not occur for the Floquet modification of the tunneling. To avoid this motion or redistribution, one needs to increase the lattice adiabatically, i.e., on a slower timescale than the Floquet drive is turned on. Nonetheless, we do a type of comparison to the data in Fig. 2 by turning on the lattice slowly (over 200 ms) to a final lattice depth to obtain a given static J and observe the condensate fraction. As an example, we increase the lattice depth to 18.2E R , which changes the tunneling from J ¼ 2π × 50 Hz at 11E R to J ¼ 2π × 12 Hz at 18.2E R , corresponding in the Floquet case to K 0 ¼ 2. The fact that there is still condensate remaining after the 200 ms turn-on is already an indication that the condensate decay rate in the unshaken lattice is slow. Measuring the condensate decay rate at this depth gives Γ=2π ¼ 2 s −1 , which is slower than the equivalent Floquet rate by a factor of 9 and slower than all the data at K 0 ¼ 2.1 shown in Fig. 2 by at least a factor of 7.

APPENDIX C: EXTRACTING THE DECAY RATES
Rate extraction.-Our data consist of a series of measured condensed fractions after various driving times. A time series typically presents an exponential-looking decay. An example for such a decay is given in Fig. 8. Since we FIG. 8. Typical decay plot. Each data point (blue dots) is an experimental realization, where the condensed fraction is measured after a variable shake time (all else kept constant). The resulting decay curve is then fitted with an exponential function (red thick line). This particular example uses K 0 ¼ 2, ω=2π ¼ 4 kHz, and V 0 ¼ 11E R . Each experimental rate in the manuscript is extracted from such a fit, and its error bar is given by AE1 standard deviation. focus on the early time decay rate, greater weight is given to earlier data points. We fit to an exponential with no offset (two fit parameters): fðtÞ ¼ Ae −Γ cf t , where A is the t ¼ 0 condensed fraction (typically, A ≥ 0.5). The rates presented in this manuscript are the extracted fit parameter Γ cf , and the uncertainties correspond to AE1 standard deviation.
A small fraction of the experimental data does not show an exponential decay. In such cases, the condensed fraction versus time looks linear, which is indicative of an acceleration relative to an exponential decay. Whenever fitting to an exponential is impossible, we take the slope, focusing on early times where the condensed fraction is more than half that of t ¼ 0. In any case, the measured decay rate under all conditions corresponds to the decay rate at initial times.
For these rates to be directly compared to the BdG prediction, additional considerations must be taken. First, let us consider a maximally unstable Bogoliubov mode with an amplitude predicted to grow as e Γ mum t . An experiment will actually detect a rate 2Γ mum , as it measures an amplitude squared (typically, the number of atoms in the unstable mode). Second, since the experiment measures how many atoms leave the BEC (to populate the modes) per unit of time, it is sensitive to the number of simultaneous maximally unstable modes, as each is a decay channel. If two modes are equally and maximally unstable, as is possible in 2D, then an additional factor of 2 is needed in the theory to compare to the experiment. This multiplemodes factor is 1 for the x-only drive and 2 for the circle and diagonal drives. All these extra factors are added to the theory plots throughout this paper: In total, the x-only drive theory is 2× and the circle and diagonal drives are 4× larger than the bare BdG rates predicted in Ref. [24]. These points are expanded in the derivation of the FBdG rate Γ mum later in this Appendix.
Rate explosion: Circular versus diagonal.-We observe that the circular decay rate increases faster than the diagonal rate for K 0 ≳ K c 0 . For example, Fig. 9 shows an enlarged version of Fig. 3, where V 0 ¼ 11E R and ω ¼ 2π × 2.5 kHz, to make this observation clearer. The rates go from being approximately equal below K c 0 to Γ circ cf > 2Γ diag cf when tunneling is suppressed. Since breaking time-reversal symmetry is necessary for many proposed schemes, this observation may be of interest to the Floquet engineering community.
Difficulties associated with dynamical rates.-In the main text, we compare the decay rates measured in the experiment to the instability rates predicted by the FBdG theory. Here, we elaborate on some intrinsic difficulties in the procedure which may affect the extracted values.
As explained in Ref. [24], for drive frequencies below the effective drive-renormalized Floquet-Bogoliubov bandwidth, there exists an entire manifold of resonant modes. While all of them contribute to expectation values of observables at very short times, only the maximally unstable mode q mum dominates the long-time BdG dynamics, and the rate associated with q mum sets the parametric instability rate. Thus, at any finite time, the FBdG dynamics is in a crossover between these two regimes, which shrinks exponentially with time. Yet the time width of this crossover also depends on the drive frequency: The higher the frequency, the smaller the instability rate, and the longer it takes for the exponential behavior to become visible.
When extracting the rates from data, effects due to this crossover become relevant. To test this, we perform exact numerical simulations of the BdG equations of motion and compute the dynamics of the excited fraction of atoms n q ðtÞ over a finite number of driving cycles, which increases suitably with the drive frequency. We then extract the instability rates using least-square fitting as the slope of log n q ðtÞ over the last eight driving cycles, to maximally eliminate transient effects. A comparison between the numerically extracted rates and the analytic theory prediction is shown in Fig. 10 for the three types of drives. Note that the agreement becomes worse at larger ω, since this decreases the rate and pushes the exponential regime to later times. This disagreement is a source of error, which is certainly relevant for the experimental determination of the rates.
Additionally, in the experiment there are strong beyond Bogoliubov effects, not captured by the FBdG theory. Because of the nonlinearity of the Gross-Pitaevskii equation which leads to saturation of the condensate depletion, the exponential BdG regime mentioned above crosses over into a third, scattering-dominated regime. In this regime, the population transferred coherently to the maximally unstable modes by the parametric resonance starts decaying into the surrounding finite-momentum modes at high FIG. 9. Decay rates in the critical region. Enlarged plot of the data shown in Fig. 3, focusing on the region between K c 0 and K 0 ¼ 2.5. The difference between the circular and the diagonal drive is clearly visible, with the circular drive rate as much as twice that of the diagonal drive rate. energy, leading to uncontrolled irreversible heating [25]. This result suggests that the instability rates can change in time and additionally obfuscates the comparison of the experiment with the FBdG theory.

APPENDIX D: HEATING DYNAMICS OF THE TRUNCATED WIGNER APPROXIMATION
While the FBdG theory is valid in the short-time regime of the dynamics, it has some serious deficiencies. Perhaps the most notable of these, when it comes to out-ofequilibrium dynamics, is the lack of particle-number conservation: The condensate is assumed to be an infinite reservoir which supplies particles to indefinitely increase the occupation of pairs of modes with finite and opposite momenta. In equilibrium, this description works well and captures the physics in the superfluid phase. Away from equilibrium, however, condensate depletion processes such as the parametric instabilities studied in this work lead to the significant depletion of the BEC, and the mean-field Bogoliubov description ultimately breaks down under typical observation times.
Particle conservation is obeyed in the truncated Wigner approximation (TWA), which also includes nonlinear interactions modeling collisions between Bogoliubov quasiparticles [51,52], and is capable of describing thermalization at later stages, due to the continuous pumping of energy into the system.
The starting point for the TWA is the Gross-Pitaevskii equation which, in the comoving real-space frame, reads (ℏ ¼ 1) i∂ t a r ðtÞ ¼ −J½a rþae x ðtÞ − a r−ae x ðtÞ − J½a rþae y ðtÞ − a r−ae y ðtÞ − ∂ 2 z 2m a r ðtÞ þ ωK 0 r · ½sinðωtÞe x þ κ d sinðωt þ ϕÞe y a r ðtÞ þ Uja r ðtÞj 2 a r ðtÞ; ðD1Þ where a r ðtÞ models the bosonic system at time t and position r ¼ ðx; y; zÞ. The kinetic energy reflects the lattice d.o.f. in the ðx; yÞ plane and the continuous transverse mode along the z axis. The periodic drive is in the ðx; yÞ plane with frequency ω and amplitude K 0 . κ d ¼ 0 for the x-only drive, and κ d ¼ 1 for both 2D drives. ϕ is the relative drive phase between x and y (ϕ ¼ 0 for a diagonal drive and ϕ ¼ −ðπ=2Þ for a circular drive). Finally, the on-site interaction strength is denoted by U (such that U=V P i;j R dzja r j 2 ¼ g). Since at time t ¼ 0 the system forms a BEC of N 0 atoms in a volume V, assuming a macroscopic occupation n 0 ¼ ffiffiffiffiffiffiffiffiffiffiffi ffi N 0 =V p in the uniform (q ¼ 0) condensate, the field a r can be decomposed as where u q and v q are the Bogoliubov modes which solve the time-independent BdG equations at t ¼ 0 [24]. Here, γ q is a complex-valued Gaussian random variable (associated with the quantum annihilatorγ q of Bogoliubov modes) with the mean and variance set by the corresponding quantum expectation values in the Bogoliubov ground state [51]. Hence, in the TWA, one draws multiple random realizations of γ q , each of which corresponds to a different initial state. One then evolves every member of this ensemble according to Eq. (D1), computes the observable of interest, and takes the ensemble averageð ·Þ in the end. For instance, one can compute the total number of excited atoms as which, due to particle number conservation, also reflects the dynamics of the condensate depletion.
In general, we expect that the condensate depletion curve shows two types of behavior: At short times, the FBdG theory applies and n ex ðtÞ ∼ expð2Γ mum tÞ grows exponentially in time. Hence, the condensate depletion curve ja q¼0 ðtÞj 2 ¼ V½n 0 − n ex ðtÞ is concave. At long times, nonlinear interaction effects in the Gross-Pitaevskii equation become important, leading to saturation, and the curvature of condensate depletion changes sign. Therefore, in the longtime regime, the curve is concave. The opposite behavior is true for the evolution of the excitations n ex ðtÞ. The curvature of the experimental data (cf. Fig. 8) suggests that the system enters well into the long-time regime. Yet, the measured decay rates appear consistent with the short-time Bogoliubov theory (main text).
To shed light on this intriguing observation, we perform TWA simulations on a periodically driven homogeneous system in (2 þ 1) dimensions and extract the short-time and long-time rates from the numerical data. We use a comparison with the BdG simulations, to separate the shorttime regime (where agreement between BdG equations of motion and TWA is expected) from the longer-time regime (Fig. 11, left). For the sake of comparison with experiments, we fit the long-time TWA growth to an exponential, even though we find that it follows a more complicated functional form. Figure 11 (right) shows a scan of the TWA rates over the effective interaction parameter g. We find that both the short-time and long-time rates are of similar strength. More importantly, they do not show a quadratic scaling in g, as predicted by Fermi's golden rule. This behavior is consistent with the experimental observations. Note the mismatch between the FBdG theory (black) and the short-time BdG simulations (blue), which arises since the most unstable mode does not yet dominate the dynamics at such short times (see Fig. 10 and the corresponding discussion). Indeed, we find an excellent agreement between BdG numerics and the FBdG theory if we extract the rates from the long-time regime. The short-time BdG rates agree qualitatively with the short-time TWA rates, as expected from the agreement seen in Fig. 11 (left). The rates are extracted from a least-square fit over the last five consecutive driving cycles of the short-time region of agreement between BdG and TWA. Since the rates are dynamical, i.e., change depending on the time window used to extract them, the curves in Fig. 11 (right) are not smooth.
We also do frequency and amplitude scans of the longtime TWA rates to look for signatures of the Bessel function J 2 ðK 0 Þ and the cusp at the critical frequency ω c , as expected from the FBdG theory and found experimentally. Unfortunately, we do not see clear signatures of such behaviors in our TWA simulations. Thus, we cannot conclude that the TWA captures the long-time thermalization dynamics of driven bosonic cold atom systems accurately. More interestingly, the rate explosion (see the main text) is also beyond the TWA dynamics, suggesting that quantum effects, such as loss of coherence, are FIG. 11. Rate comparison: TWA versus BdG. Left: Numerical simulation of the excitations growth using TWA (solid lines) and FBdG (dashed lines). The TWA curves change the curvature beyond the regime of validity of BdG. Right: Excitation (i.e., condensate depletion) growth rate against effective interaction strength g for the FBdG theory (black), BdG short-time evolution (blue), TWA shorttime evolution (red), and TWA long-time evolution (green). The parameters are K 0 ¼ 2.1, ω=J ¼ 20.0, and n 0 ¼ 50.0. We use a system of 80 × 80 × 101 momentum modes in the ðx; y; zÞ direction, respectively. The TWA data are averaged over 50 independent realizations, and the error bars (shaded area) are computed using a bootstrapping approach. important for describing this phenomenon. Another possible reason for the disagreement is the single-band approximation, as its validity for the Floquet system has not been fully understood so far.

APPENDIX E: MOST UNSTABLE MODES FOR LINEAR, DIAGONAL, AND CIRCULAR LATTICE DRIVES WITHIN FBDG THEORY
In this Appendix, we briefly revisit the derivation for the most unstable mode within the FBdG theory [26] and extend the results to diagonal and circular drives. The takehome message is that the critical saturation frequency ω c , which defines the position of the cusp in the instability rate curves, coincides for the linear and circular drives; for the diagonal drive, the value is twice as large on the square lattice. It is straightforward to extend the analysis below to other lattice geometries.
The starting point is the BdG equations of motion (EOM) for the mode functions u q ðtÞ and v q ðtÞ in the rotating frame: where εðq; tÞ ≥ 0 denotes the time-dependent free lattice dispersion relation (shifted by the chemical potential so that it is non-negative).
To analyze the effects of parametric instabilities, we first isolate the time-average term and write εðq; tÞ ¼ ε eff ðqÞ þ gðq; tÞ, which will separate the right-hand side of the BdG EOM into an effective time-averaged term and a time-periodic perturbation. Following Ref. [24], we now apply a phase rotation, followed by a static Bogoliubov transformation with respect to ε eff ðqÞ: where coshð2θ q Þ ≡ ½ε eff ðqÞ þ g=E Bog eff ðqÞ and sinhð2θ q Þ ≡ g=E Bog eff ðqÞ with E Bog eff ðqÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ε eff ðqÞ½ε eff ðqÞ þ 2g p the time-averaged Bogoliubov dispersion. The Bogoliubov transformation diagonalizes the effective static timeaverage term, while the phase rotation makes it easier to identify the parametric resonant terms (see Appendix C in Ref. [24] and Supplemental Material in Ref. [25] for an application to the paradigmatic parametric oscillator). The BdG EOM now read where W q ðtÞ ¼ g q ðtÞcosh 2 ðθ q Þ þ g −q ðtÞsinh 2 ðθ q Þ 0 0 −g −q ðtÞcosh 2 ðθ q Þ − g q ðtÞsinh 2 ðθ q Þ ! ; h q ðtÞ ¼ 1 2 ½g q ðtÞ þ g −q ðtÞ: ðE3Þ By the definition of g q ðtÞ, the diagonal matrix W q ðtÞ has a vanishing period average and, hence, does not contribute to the parametric resonance to leading order. At the same time, the off-diagonal term in Eq. (E2) will be dominant, whenever h q ðtÞ interferes with the phase term e 2iE Bog eff ðqÞt constructively. This latter condition gives the resonant frequencies. In a sense, h q ðtÞ can be thought of as an effective periodic drive, the amplitude of which, multiplied by the prefactor sinh 2θ q ¼ g=E Bog eff ðqÞ, determines properties of the resonance, such as the maximally unstable mode.
Let us now analyze the function h q ðtÞ for the three types of drives we study in the main text.
Because q circ mum is again located at the side edge of the 2D Brillouin zone, the saturation frequency for circular drive ω circ c equals half the effective 2D Bogoliubov bandwidth (similar to a linear drive).
The above analysis is performed for a pure 2D lattice system. However, adding the transverse tube dimension is straightforward and does not change these conclusions [24]; cf. Fig. 10 for a numerical simulation in the presence of a tubelike transverse direction.
Finally, note that the instability rate γ derived above corresponds to the growth rate of the bosonic operators b q ðtÞ and b † q ðtÞ, which are not observable. Observables, consisting of bilinears of the bosonic operators, will thus display a parametric instability rate for each maximally unstable mode. Furthermore, an atom number variation per unit time (an experimental rate) depends on how many such modes are present simultaneously. For the linear drive, there exists a single maximally unstable mode ðπ; 0Þ, while there exist two equally unstable modes for the diagonal drive [fðπ; πÞ; ð−π; πÞg] and for the circular drive [fðπ; 0Þ; ð0; πÞg]. Therefore, we expect an additional factor of 2 in Γ mum for the two 2D drives. In total, As a final remark, we note that, while a choice of small K 0 explores the 2D quasiparabolic dispersion found at small BEC momenta, the above derivation shows that the problem still cannot be treated as rotationally symmetric: The resonant modes are determined by the drive frequency through the resonance condition. Therefore, even for small drive amplitude K 0 , resonant Bogoliubov modes can be in the middle of the dispersion, where rotational symmetry is lost. Hence, the distinctions between 1D and 2D drives are expected, and observed, to hold even in the small K 0 limit.