Parametric instabilities of interacting bosons in periodically-driven 1D optical lattices

Periodically-driven quantum systems are currently explored in view of realizing novel many-body phases of matter. This approach is particularly promising in gases of ultracold atoms, where sophisticated shaking protocols can be realized and inter-particle interactions are well controlled. The combination of interactions and time-periodic driving, however, often leads to uncontrollable heating and instabilities, potentially preventing practical applications of Floquet-engineering in large many-body quantum systems. In this work, we experimentally identify the existence of parametric instabilities in weakly-interacting Bose-Einstein condensates in strongly-driven optical lattices through momentum-resolved measurements. Parametric instabilities can trigger the destruction of weakly-interacting Bose-Einstein condensates through the rapid growth of collective excitations, in particular in systems with weak harmonic confinement transverse to the lattice axis.

Floquet engineering has proven to be a powerful technique for the design of novel quantum systems with tailored properties, unattainable in conventional static systems [1][2][3]. It is based on the design of time-periodic systems, whose stroboscopic evolution is governed by an effective time-independent Hamiltonian featuring the desired properties. Floquet engineering is captivating due to its conceptual simplicity and its potentially farreaching applications for engineering novel states of matter. For instance, it has been used to manipulate the electronic properties of solid-state systems [4][5][6][7], to realize time crystals [8,9], to engineer artificial magnetic fields and topological Bloch bands in cold atoms [10][11][12][13][14][15], photonics [16][17][18] and superconducting circuits [19], to generate density-dependent gauge fields [20,21] and to explore the rich physics of lattice gauge theories [22].
The complex interplay between periodic driving and interactions poses theoretical and experimental challenges. In time-periodic systems, energy conservation is relaxed due to the possibility to absorb and emit energy quanta from the drive, and any driven ergodic system is expected to eventually heat up to infinite temperatures [23,24]. Recent experiments [25][26][27][28][29] have addressed this problem for interacting atoms in shaken optical lattices. In particular, it has been shown [27,[30][31][32][33] that heating rates are well captured by a Floquet Fermi's Golden Rule (FFGR) approach, if they are evaluated at sufficiently long times. This approach suggests that the long-time dynamics is dominated by incoherent two-body scattering processes [32]. In contrast, the onset of heat-ing in bosonic systems is expected to be triggered by coherent processes [34][35][36][37].
First evidence for parametric instabilities in amplitude-modulated optical lattices has been found indirectly via spectroscopic measurements [34,38]. In this work, we directly reveal the existence and nature of parametric instabilities by measuring the momentum distribution of weakly-interacting bosons in a periodically-driven one-dimensional (1D) optical lattice (Fig. 1a). These instabilities exist whenever the energy quantum ω associated with the drive matches the energy of a collective excitation, as dictated by the effective Bogoliubov dispersion E B eff (q, J eff ), where q denotes the momentum of the excitation and J eff is the effective tunnel coupling renormalized by the drive [39] (Fig. 1b). In a strictly 1D lattice, the strongest instability occurs at the two-photon resonance, 2 ω = 2E B eff (q x res , J eff ). In the single-band approximation, there exists a stable parameter regime without parametric instabilities for ω > ω sat ≈ W eff / , where W eff = 4|J eff |(4|J eff | + 2g) is the bandwidth of the effective Bogoliubov dispersion and g is the interaction energy (Ref. [37] and App. A). In a 3D system with weak harmonic confinement transverse to the lattice axis (Fig. 1a), however, no stable parameter regime exists [27,[31][32][33][34]37] and parametric instabilities occur via the closely-spaced transverse modes [34,37,38] [regime (II) in Fig. 1c]. Moreover, our numerical simulations indicate that a harmonic confinement along the lattice axis further prevents the existence of a stable region, even in a true (a) Schematic of the driven 1D optical lattice with tunneling J, generated by two laser beams with frequencies ω1,2, and weak harmonic transverse confinement. Modulating ω2(t) periodically with frequency ω generates a force F (t) with period T = 2π/ω. (b) Illustration of the effective 1D Bogoliubov dispersion E B eff (q x , J eff ), with width W eff reduced according to the Floquet-renormalized tunnel coupling J eff (blue) compared to the static lattice case (gray). Parametric resonances at E B eff = ω induced by the modulation lead to instabilities centered around the most unstable mode q x mum ; = h/(2π) is the reduced Planck's constant. (c) Momentum of the most unstable mode qmum = (q x mum , q ⊥ mum ) as predicted by Bogoliubov theory, which shows a clear separation between lattice (I) and transverse (II) degrees of freedom, that occurs at the saturation frequency ωsat ≈ W eff ; d is the lattice constant.
1D geometry (App. B). The appearance of instabilities is common in static weakly-interacting bosonic lattice systems. For instance, Landau instabilities can occur when the condensate is prepared at a finite quasimomentum, where the effective mass in the band structure is negative; such configurations can also display dynamical instabilities, where Bogoliubov excitations grow exponentially [40][41][42][43][44][45][46]. We emphasize that the origin of such instabilities is different compared to those revealed in this work. These instabilities exist even for a condensate initially at rest, in contrast parametric instabilities originate from the timedependent nature of the drive. While most experiments have studied the appearance of instabilities via decoherence or loss-rate measurements. We note that instabilities are not necessarily detrimental but can result in interesting phenomena, such as parametric amplification, four-wave mixing [44][45][46][47] and pattern formation [48][49][50][51].
The experiment starts by loading an almost pure Bose-Einstein condensate (BEC) of about N = 3.7(4) × 10 5 39 K atoms within 100 ms into a 1D optical lattice aligned along the x axis, with lattice constant d = 425 nm and depth V lat = 11.0(3) E R , where E R = h 2 /(8md 2 ) = h × 7.1 kHz is the recoil energy and m the mass of an atom. Additional confinement is provided by an optical dipole trap. The harmonic trapping frequencies of the combined potential are ω r /(2π) = 26 (2) Hz in the xyplane and ω z /(2π) = 204(3) Hz in the vertical direction. The lattice is created by interfering two laser beams with λ = 736.8 nm under an angle of 120 • . Its position is modulated by varying the frequency of one lattice laser beam (Fig. 1a), ω 2 (t) = ω 1 +2πν sin(ωt+ϕ), where ϕ is the phase of the drive. The lattice modulation is turned on suddenly in order to be able to observe the presence of collective excitations after few modulation periods (we verified that ramping up the modulation amplitude within five cycles does not modify our main results). We hold the atoms in the modulated lattice for integer multiples of the driving period and determine the momentum distribution by performing bandmapping (the modulation is turned off abruptly followed by a 100µs-long linear rampdown of the lattice) and subsequent time-of-flight (TOF) imaging (Fig. 2a).
In the reference frame of the lattice, the modulation leads to a time-varying force F (t) = F 0 cos(ωt + ϕ), with F 0 = mdνω [39] and the time-dependent tight-binding Hamiltonian describing the dynamics takes the form: whereâ † i,r ⊥ andâ i,r ⊥ are the bosonic creation and annihilation operators on lattice site i and transverse position r ⊥ ,n i,r ⊥ is the corresponding number operator, J/h = 108(7) Hz is the tunnel coupling, U is the on-site interaction, K = F 0 d,Ĥ ⊥ denotes the kinetic energy along the transverse direction andĤ harm the 3D harmonic confinement. The Feshbach resonance of 39 K at 403.4(7) G enables us to work in the weakly-interacting regime at the scattering length a s = 20a 0 (with a 0 the Bohr radius), where a mean-field approach is expected to be valid. The Feshbach resonance is crucial in order to find experimental parameter regimes, where parametric instabilities could be clearly identified in momentumspace images. Additional nonlinear effects rapidly counteract the exponential growth of the instabilities ( [37] and App. C), leaving only a small window of suitable parameters to observe it.
In the non-interacting limit, Floquet theory predicts that the dynamics is well described by a timeindependent Hamiltonian, with renormalized tunnel coupling J eff = JJ 0 (α) [1][2][3]39]. Here J ν is the νth-order Bessel function of the first kind and α = K/( ω). Note, that α is independent of the modulation frequency ω. The measurements are performed in a strong-driving regime, 1 < α < 2, where the effects of the drive are nonperturbative, while the minimum of the effective disper-  Difference image at t = 14 T , obtained by subtracting the condensate profile at t = 0. The lower panels show the 1D profiles along q x , resulting from summation of the difference image perpendicular to the lattice axis for the left and right excitation peak together with asymmetric Lorentzian fits (black solid lines), used to determine the amplitude A x and position q x max (dashed lines) of the excitation peaks [53]. The left and right plot show the profiles resulting from summation along the lattice axis with Lorentzian fits to determine the transverse full-width at half maximum ∆q y (dashed lines) of the peaks. Here, we integrate the density profile in a region of interest containing the left or right peak only. sion remains at q x = 0 [39,46]. To avoid single-particle inter-band resonances [52], the modulation frequency is chosen well below the first single-particle bandgap of the static lattice ∆ 21 /h = 41.6(5)kHz. Indeed, we do not observe excitations to higher bands during the measurements, which would be visible in the TOF images (Fig. 2a).
In order to identify instabilities in the system, we monitor the appearance of collective excitations by measuring the momentum distribution of the atoms as a function of modulation time for various modulation frequencies and amplitudes. We find that after few modulation periods a small fraction of atoms is excited into additional momentum components that are distinct from the initial condensate at q x = 0 (Fig. 2a). The amplitude of these modes grows as a function of the modulation time and excitation peaks eventually start to broaden after long modulation times (see t = 17 T in Fig. 2a) due to satu-ration and nonlinear effects (App. C).
We attribute these excitations to parametric instabilities that occur whenever a resonance exists between the energy associated with one drive quantum ω and the Bogoliubov spectrum E B eff (q) [see Eq. (A2) in App. A] of collective excitations [34,35,37]. In the presence of transverse modes there is no stable parameter regime, since the transverse kinetic energy has no upper bound. There is always a set of resonant excitations determined by the resonance condition which causes the system to be necessarily unstable. Each collective mode grows at a different rate Γ q and we denote the one associated with the dominating growth rate Γ = max q Γ q the most unstable mode ("mum") with the corresponding momentum q mum . For simplicity, we henceforth neglect the 3D harmonic confinementĤ harm in our theoretical analysis and set the transverse kinetic term equal to a free-particle kinetic en- 2mn j,q ⊥ . For weak harmonic confinement, where the transverse modes are closely spaced, we expect this to be a good approximation. We note, however, that the harmonic confinement along the lattice axis largely modifies the stability of a strictly 1D system, where it prevents the existence of a true stable parameter regime, as the energy spectrum becomes unbounded (App. B).
Following the approach developed in [35,37] based on the time-dependent Bogoliubov-de-Gennes (BdG) equations of motion [53], we find approximate analytic solutions for the momentum q mum and the growth rate Γ of the most unstable mode (Ref. [37] and App. A). To lowest order in this perturbative treatment, there is a clear separation between the lattice and transverse degrees of freedom, resulting in two distinct regimes ( Fig. 1c): (I) ω < ω sat : q x mum < π/d and |q ⊥ mum | = 0, (II) ω > ω sat : q x mum = π/d and |q ⊥ mum | > 0. In the first regime (I), the modulation mainly couples to excitations along the lattice direction and the instability is not affected by the transverse modes (|q ⊥ mum | = 0); here g = nU denotes the interaction energy and n is the mean density (see App. B for a discussion of the inhomogeneous density profile of a harmonicallytrapped gas). When reaching the second regime (II), in contrast, the transverse degrees of freedom dominate; q ⊥ mum becomes finite and grows according to To quantify the position of the most unstable mode experimentally from the TOF images (Fig. 2a), we subtract the mean initial condensate profile at t = 0 from  Position of the peak maximum q x max along q x as a function of modulation time. The most unstable mode q x mum is defined as the weighted mean of all q x max for t ≤ ts. Its value is shown as the solid black line; the shaded blue bars denote its standard error of the weighted mean. The shaded gray area defines the range of hold times over which the transverse width ∆q y is averaged. Each data point is an average over ∼ 10 individual experimental realizations, error bars indicate the standard deviation. To minimize systematic deviations in the bandmapped images, we average < ∼ 5 realizations for a modulation phase ϕ = 0 and < ∼ 5 with ϕ = π (Fig. S1).
each individual TOF image, measured at t > 0. A typical result is shown in Fig. 2b (a more detailed description of the data analysis can be found in the Supplementary Material [53]). To characterize the excitations along the lattice direction q x mum , we integrate the 2D difference profile perpendicular to the lattice axis (lower panels in figure 2b). The resulting 1D profile is divided into two parts excluding the negative part at small |q x | (purple), which arises due to the small depletion of the condensate. We determine the position q x max and amplitude A x of the excitation peaks by fitting asymmetric Lorentzian functions to each peak [Eq. (S.1) in the Supplementary Material].
The transverse momentum component q ⊥ mum of the most unstable mode is masked by the initial momentum spread of the condensate and the width of the parametric resonance. Nevertheless, it manifests itself in a broadening, which we monitor by extracting the transverse width ∆q y of the excitation peaks: We integrate the 2D difference profile along the lattice direction using a region of interest that contains only the left or the right peak. We fit the resulting 1D profiles (left and right panel in figure 2b) with a symmetric Lorentzian and extract the full width at half maximum. The fit parameters obtained for the left and right peak are then averaged over all images for each modulation time t > 0.
In Fig. 3 we show a typical time series for the position q x max and amplitude A x of the dominant mode (full data set is shown in the Supplementary Material [53]). While there are remnants of a moderate exponential growth at short times t < t s we predominantly find an explicitly time-dependent growth rate. This is due to incoherent processes caused by interactions between collective excitations and between the excited modes and the condensate [37]. This behavior is confirmed by our numerical simulations beyond BdG theory discussed in App. C. In order to quantitatively compare our experimental data with analytical expressions obtained using Bogoliubov theory we restrict the analysis to a small time window t < t s before saturation and nonlinear effects start to dominate. This is motivated by the good agreement between numerical simulations and analytical formulas for short times (Fig. 9 in App. C). The parameter t s is evaluated by fitting a piecewise function consisting of a linear and a constant part to ln(A x ) individually for the left and right excitation peak (Fig. 3a) and averaging the two results [53]. The corresponding growth rates (for t < t s ) are in good agreement with Bogoliubov theory (Fig. 10 in App. C), validating this approach. The position of the most unstable mode q x mum is then defined as the average over the peak positions q x max for all modulation times t < t s . We observe a decrease of q x max at modulation times t > t s indicating the onset of additional scattering events not captured by Bogoliubov theory. In order to study the behavior along the transverse direction, we analyze the peak width ∆q y at the end of the short time regime, where the amount of transverse excitations is expected to be maximal (full data set is shown in the Supplementary Material [53]). In order to reduce statistical errors, we average ∆q y over modulation times between t s − 1 and t s + 1 in units of the modulation period ( Fig. 3b) to obtain the experimental value for the width ∆q y lin [53]. The position of the most unstable mode along the lattice axis q x mum and the transverse width ∆q y lin are measured for various modulation parameters around the saturation frequency, where we expect a crossover between parametric instabilities dominated by excitations along the lattice axis [regime (I)] and those that are facilitated by the presence of transverse modes [regime (II)]. The results are shown in Fig. 4.
We find that q x mum indeed increases with ω until it saturates at q x mum ≈ π/d, the edge of the Brillouin zone, at a frequency ω sat that depends on the driving amplitude α (Fig. 4 upper panel). The small deviation of the measured positions from π/d is mainly due to the short TOF used in the experiment [53]. The saturation frequency ω sat calculated from Eq. (3) matches the experimental data well for an interaction parameter g = 11.5 J. From our measured in-situ density profiles we obtain an interaction parameter g max ≈ 8 J in the center of the trap. The deviation is most likely due to a systematic uncertainty in the atom number calibration [53]. At the same time, we observe that the width ∆q y lin (ω) transverse to the lattice axis starts to increase with frequency for ω > ω sat (α), simultaneously with the saturation of q x mum at the BZ edge ( Fig. 4 lower panel), as expected from lowest-order perturbation theory Eqs. (3)-(4). Moreover, the transition point ω sat (α) decreases for larger driving parameters α, in line with the reduction of J eff .  We further observe that the shape of the q x mum -curve differs from the theoretical prediction for low modulation frequencies in particular for the large modulation amplitude α = 1.78. We attribute this deviation mainly to the following effects: First, the short-time window, during which we can clearly identify the most unstable mode strongly depends on the modulation parameters. The time t s , which marks the end of the short-time window, decreases with decreasing modulation frequency and increasing modulation amplitude (Fig. S2). This complicates the identification of the peak maxima. If t s becomes too small, there is not enough time for the most unstable mode to dominate over the other excited modes. Second, the initial width of the momentum distribution of the condensate at t = 0 poses a fundamental limitation that prevents us from measuring positions smaller than ≈ 0.4π/d. Moreover, for large modulation amplitudes higher-order corrections to the analytical formulas Eqs (3)-(4) become relevant (App. A).
In summary, we have demonstrated the first direct evidence for parametric instabilities in shaken optical lat-tices via momentum-resolved measurements. Our experiments were performed with a large 3D system of weaklyinteracting bosons, where exact numerical calculations including the harmonic trap are not feasible. Tuning the scattering length with a Feshbach resonance we were able to identify and address experimental parameter regimes, where the short-time dynamics is well-described by BdG theory. While the instability rates are time-dependent due to competing processes that dominate on different timescales, the momentum of the most unstable mode turned out to be a reliable observable. The obtained results are in agreement with the analytic approach derived in Ref. [37], which enables the development of an intuitive understanding and allows us to identify stable parameter regimes in driven lattice models from simple energetic arguments. We were able to verify the existence of key bottlenecks in current experimental settings with weak transverse confinement [10,11,14,26,[54][55][56]] that need to be overcome by freezing the transverse degrees of freedom and generating a box-type longitudinal confinement [57,58]. Parametric instabilities can indeed lead to a depletion of the condensate, where subsequent scattering events result in large heating rates. Our results are of strong interest for future experiments based on Floquet engineering [3], as they indicate the necessity to engineer full 3D lattice systems, where stable regimes can be found [35,37]. Parametric resonances are expected to be present whenever the BdG equations of motion include time-periodic features, and hence may turn out to play an important role in a wide family of Floquet-engineered systems, such as periodically-driven superfluids [3] and superconductors [59], photonic devices [60,61], but also in the context of cosmology [62].
We note that during the completion of this work, shorttime heating rates, which are expected to be dominated by parametric instabilities, have been investigated in modulated 2D lattices [63].
We thank T. Boulier Quantum Simulation MURI, AFOSR-MURI Photonic Quantum Matter (award FA95501610323). We used QuSpin [64,65] to perform the numerical simulations. The authors are pleased to acknowledge that the compu-tational work reported on in this paper was performed on the Shared Computing Cluster which is administered by Boston University's Research Computing Services. Here, we would like to recall the analytical method developed in Ref. [35] to extract the instability properties of the system within the Bogoliubov approximation. As mentioned in the main text, for simplicity, we neglect the harmonic confinement, so that the transverse kinetic energy is determined by the free-particle dispersion relation. It has been shown that the Bogoliubov equations of motion [Eq. (S.5) in Sect. S3 A of the Supplementary Material [53]], can be mapped to a parametric oscillator model [36,74], a seminal model of periodically-driven harmonic oscillator known to display parametric instabilities as soon as the drive frequency approaches twice the natural frequency. To see that, one should perform a series of suitable changes of basis and reference frames [37]. Applying the Rotating Wave Approximation (RWA), we keep the leading-order harmonic and recast the Bogoliubov equations into the form: where1 is the identity matrix, denotes the effective (time-averaged) Bogoliubov dispersion, and we introduced the amplitude For the sake of simplicity, we have here restricted ourselves to the dominant harmonic of the drive, and dropped all terms in Eq. (A1) that are irrelevant regarding the occurrence of instabilities -a simplification that was rigorously established in Ref. [37]. Perturbation theory-As can be seen from a RWA treatment of Eq. (A1), each momentum mode q will display a dynamical instability (characterized by an exponential growth of its population), whenever the drive frequency ω approximately matches its effective Bogoliubov energy E B eff (q), i.e., ω ≈ E B eff (q). The analytical method to extract the associated instability rate is detailed in Refs. [37,74] and relies on a perturbation theory formulated in both A q (the small parameter of the expansion, which is assumed to be smaller than 1 for the expansion to converge), and the detuning from the resonance δ q ≡ ω − E B eff (q). In brief, the findings are the following: 1. Zeroth order: the instability occurs only on resonance, i.e., ω = E B eff (q) and the instability rate is given by for all q fulfilling the resonance condition ω = E B eff (q) and is zero for all other modes. Interestingly, not all resonant modes necessarily have the same instability rate. We focus on the most unstable mode, which has the largest rate Γ = max q Γ * q , which results in the following rates for the two regimes introduced in the main text: (I) ω < ω sat : (II) ω > ω sat : The corresponding momentum q mum of the most unstable mode to lowest order in this pertubative treatment is defined in Eq. (3) and (4) in the main text.
2. First order: we find that instability does not only occur on resonance but still arises in a finite window around the resonance point. The associated instability rate is given by if the argument of the square root is positive, and Γ q ≡ 0 otherwise. It is maximal at resonance and decreases with the distance from resonance, until it vanishes at the edges of the instability domain [37,74].
3. Second order: We no longer have explicit analytical expressions, but we find that Γ q is the solution of an implicit equation which can be numerically solved [37,74]. While the width of the resonance is unaffected, we obtain that the instability rate is no longer maximal on resonance, but that the maximum instability point is slightly shifted from the resonance.
In Fig. 5 we show a comparison between the analytic formulas Eq. (3) and (4), the perturbation theory up to second order and a numerical simulation of the exact time-dependent BdG equations Eq. (S.5). We find that, for large modulation amplitudes, higher-order corrections become important because A q [Eq. (A3)] is not a small parameter any more; more specifically J 2 (α) is no longer small compared to the effective Bogoliubov dispersion E B eff (q). We believe that this explains why the measured q x mum -curve in Fig. 4 deviates more from the zero-order BdG theory for large modulations amplitudes, i.e. α = 1.78, as compared to the rather good agreement observed for α = 1.44. We stress that the BdG analysis we present here does not rely on the inversefrequency expansion, which is routinely used to determine the effective Hamiltonian in the context of a timeperiodic Schrödinger equation. In this section, we analyze modifications to the translationally-invariant Bogoliubov-de-Gennes (BdG) theory [Sect. S3 A in the Supplementary Material [53]] due to the harmonic trap along the lattice axis. We constrain the discussion to 1D systems for simplicity. Energy (J) 35

Unbounded energy spectrum
In most optical lattice experiments there is a harmonic confinement along the lattice axis, and therefore the energy spectrum of the system is no longer bounded from above. Figure 6 (orange dots) shows the energy spectrum of the effective time-averaged non-interacting Hamiltonian in the presence of the harmonic trap, and the corresponding drive-renormalized Bogoliubov dispersion (blue dots). Details about the numerical simulations of the inhomogeneous system can be found in Sec. S4 B of the Supplementary Material [53]. Notice the occurrence of states above the Bogoliubov bandwidth (blue dashed line), which can be excited resonantly by the drive. In momentum space, these states occupy modes in the vicinity of q x = π/d, and typically have weight over a finite range of momenta. Therefore, the presence of a weak harmonic confinement allows the system to absorb energy even for ω > ω sat , where unconfined systems are shown to be stable [37]. As a result, we do not expect a truly stable parameter regime to exist in harmonically confined systems, even in the absence of transverse modes.

Modified saturation frequency ωsat
The expressions for the saturation frequency ω sat and the effective bandwidth W eff , derived using translationally-invariant Bogoliubov theory, depend on the density-renormalized interaction parameter g. To zeroth-order (App. A) the effective bandwidth for the 1D lattice is obtained from Eq. (A2) by setting q ⊥ = 0. In the presence of a harmonic confinement, the condensate profile obeys Thomas-Fermi theory and is no longer uniform, which induces a position dependence in g → g(x). Since the energy spectrum is unbounded, there is no natural energy bandwidth to use. Hence, we need to determine the appropriate value for g used to compute the saturation frequency ω sat . We observe a correlation between the energy at which the states in the Bogoliubov spectrum rapidly become quasidegenerate [cf. Fig. 6 inset], and the energy scale ω sat (g max ) with g max = max The same scale also coincides with the energy, starting from which the Bogoliubov states attain a significant occupation of the q x = π/d mode. This suggests that, in the presence of a harmonic trap, the maximum value g max can be used to estimate the saturation frequency ω sat .
To test this conjecture, we performed numerical BdG simulations of 1D lattices comparing the position of the most unstable mode with and without harmonic confinement. Indeed, we find a qualitative behavior similar to the homogeneous system, i.e., the momentum q x mum (ω) increases with ω until it saturates at q x mum ≈ π/d. We can also quantitatively compare the α-dependence of the saturation frequency ω sat (α) for the homogeneous and trapped systems as follows: (1) we fix a trap frequency and determine the value of g max = U max x n TF (x) from the Thomas-Fermi profile of the condensate wavefunction n TF (x). (2) we simulate the dynamics of a homogeneous system using the same value of g. If Eq. (B1) for g = g max provides a correct description, this procedure will result in the same value of ω sat for the trapped and homogeneous simulations by construction. The values of ω sat (α) are extracted numerically from the q x mum vs. ω curves for every fixed value of α, as follows: (i) we compute numerically the momentum distribution profile of the Bogoliubov modes as a function of momentum q x and time t, evolved under the periodic drive. We do this for a grid of various α and ω points. (ii) we evolve the system for 20 driving cycles which is enough to single out the most unstable mode that grows exponentially according to the BdG equations [53]. (iii) we extract the fastest growing mode q x mum from the latest time slice for every point on the (α, ω) grid. (iv) for every fixed α, we determine the saturation frequency, by finding the frequency for which q x mum reaches π/d for the first time upon increasing ω. Figure 7 shows that the behavior of the saturation frequency ω sat (α) agrees with analytic predictions of Eq. (B1) for g = g max . We observe that the agreement gradually becomes limited in the regime of large α where the effective kinetic energy of the system is parametrically reduced for strong interactions g. This regime of large α and large g/J is precisely where higher-order corrections to Eq. A3 become pronounced (see the discussion in App. A). This explains the observed mismatch in the g/J = 10.2 curves at large α in Fig. 7.

Modified BdG instability rates
In Fig. 8 we show numerically evaluated instability rates of 1D lattice models using BdG simulations. We compare the drive frequency dependence for the inhomogeneous vs. homogeneous lattices. Ideally, a homogeneous 1D system becomes stable for ω > ω sat , since there are no states above the lattice bandwidth [75]. The simulation reveals that the presence of the harmonic trap does not lead to any further shift in the position of ω sat , provided one compares a homogeneous system of interaction strength g to a harmonically confined system with g max = g. In the latter case, ω sat is approximately given by the Bogoliubov bandwidth evaluated at g max [ Fig. 7]. For larger modulation frequencies ω, the trapped system displays a distinctly different behavior compared to the homogeneous system. As already anticipated, due to the unbounded nature of the energy spectrum, there is no true stable (i.e. Γ ≡ 0) parameter regime and the confined system can always absorb energy via resonant processes. Nonetheless, we find a decrease of the instability rates with increasing drive frequency.
Due to lack of translation invariance, it is not feasible to carry out an analytical calculation to determine the precise decay law Γ(ω) for ω > ω sat in the presence of the harmonic trap. However, assuming the Local Density Approximation for a weak-enough harmonic confinement and keeping in mind the unbounded structure of the spectrum, we can apply Eqs (A5)-(A6) in all spatial regions of approximately constant density, according to which Γ ∝ ω −1 for ω > ω sat ; while this does not predict the exact functional form of Γ(ω) observed in Fig. 8, it presents an approximate argument for the observed rate decay.
As a result, for trapped systems one can recover an approximately stable regime, where the instability rates are small compared to the duration of the experiment. An alternative way to mitigate the problem of non-vanishing rates at large drive frequencies in experiments could be provided by the use of uniform box traps [57,58].

Appendix C: Instability rates
Besides the position of the most unstable mode discussed in the main text, a characteristic property of parametric instabilities are the associated instability rates Γ, which after sufficiently long times are dominated by the growth of the occupation of the most unstable mode n qmum . Here, we present an experimental and numerical study of these rates.

Influence of mode coupling
In order to provide more insight on the timedependence of the instability rates found in the experiment ( Fig. 3a and Fig. S2 in the Supplementary Mate- rial [53]), we performed numerical simulations (see [53] for details) on a homogeneous hybrid (translationallyinvariant) 2D system [76], composed of one lattice and one continuum direction, based on three different approximation methods [53]: (i) The linearized BdG equations, which capture the parametric instability at short times, i.e., before saturation effects, such as particle-number conservation, and non-linear effects associated with the Gross-Pitaevskii equation, become significant. (ii) The weak-coupling conserving approximation (WCCA) [36], where particle-number conservation is restored, and which couples the condensate mode to the excitations to leading order in the interaction strength U . This method keeps track of the number of atoms scattered into finite-momentum modes, as well as the back action of the Bogoliubov quasiparticles onto the condensate. The WCCA, however, does not capture collisions between quasiparticles. Hence, it does not offer any insight on the thermalization dynamics at longer times, during which the system heats up steadily to an infinitetemperature state. The BdG and WCCA approaches have already been compared in Ref. [37] and their validity at short times is further confirmed in this work by comparing their predictions to a third approach: (iii) the Truncated Wigner approximation (TWA), which produces thermalizing dynamics and, even though quantum effects are only partially captured in this semi-classical approximation, it has recently been demonstrated that classical Floquet systems thermalize in a very similar way to quantum models [68][69][70]. Figure 9 shows the time evolution of the occupation of the most unstable mode. In both, the WCCA and TWA simulations, the obtained curves directly reflect the condensate depletion dynamics. At short times, we find an agreement between the three theoretical approaches, which further improves as one decreases the on-site interaction strength U , while keeping g = nU fixed. At later times, the three approximations exhibit different behaviors: due to the lack of particle conservation, the BdG curves grow exponentially in an unphysical and indefinite manner. In contrast, the WCCA and TWA curves show clear manifestations of saturation effects, indicating that the instability rate is a truly time-dependent physical quantity. After an intermediate transient, all momentum modes are equally populated, which is expected for an infinite-temperature state that is reached at long times. These results offer a qualitative explanation for the saturation of the peak growth observed in the experiment, indicating the importance of saturation effects at intermediate modulation times, and further highlight the advantage of momentum-resolved measurements for revealing parametric instabilities (Fig. 4).

Measured instability rates
Motivated by the numerical analysis discussed above, we quantitatively study the exponential growth rate of the most unstable mode at short times, before saturation effects dominate the dynamics. To determine the saturation point of the amplitudes and the instability rates from the experimental data, we fit a piecewise function consisting of a linear and a constant part to to the logarithm of the peak amplitude A x [53]. In BdG theory, the peak amplitudes along the x-lattice direction grow in time according to A x ∝ e 2Γt . Hence, from the fitted slopes m (in the regime t < t s ), we extract the instability rates as Γ = m/2 and average over the left and right peak. In Fig. 10, we show a comparison between the experimentally observed rates and the theoretical expectation based on the analytical BdG predictions [Eqs (A5)-(A6)]. We further display the long-time heating rates, which have shown to be captured by Floquet Fermi's Golden Rule (FFGR) in Ref. [27]. For completeness the relevant equations are summarized in the subsequent section.
We find that the experimental rates [dots in Fig. 10] are of similar magnitude as compared to the BdG predictions and about several orders of magnitude larger than the FFGR predictions [32]. The inset of Fig. 10 shows all rates on a linear scale, to illustrate the dependence on the modulation parameters. From BdG theory, it is expected that the rates have a maximum at the saturation frequency ω sat [where ω equals the effective bandwidth]. This maximum value is expected to increase with larger modulation amplitude.
Although the measured rates lie in the same order of magnitude as the BdG rates, the BdG-predicted parameter dependence is not directly confirmed by the measurements. Given the narrow window determined by t s during which the most unstable mode dominates, this is not surprising. More importantly, the discrepancy between the measured rates and the FFGR predictions shows that we are indeed probing the system on sufficiently short timescales, where instabilities are mostly driven by coherent processes.

Floquet Fermi's Golden Rule (FFGR)
Heating on long timescales dominated by incoherent processes is well described by FFGR as studied in detail in Ref. [27]. The exponential decay of the BEC atoms caused by the modulated lattice can be described aṡ where Γ FFGR = κN 2/5 . The loss rate κ is independent of the atom number and given by: u l = 0.75 for l even 0.15 for l odd , , w 0 (x) denoting the wannier function of the lowest band and J ν being the ν-th Bessel function of the first kind describing ν-photon scattering processes. For the model parameters used in our study, the rates are dominated by processes up to second order, keeping contributions up to l = 2. The FFGR rates are displayed in Fig. 10 [bright solid lines]. It is evident that the FFGR rates are several orders of magnitude smaller than the BdG rates [dark solid lines]. This failure of FFGR to capture the short-time heating dynamics can be traced back to the exponentially dominating short-time parametric instability effect which, in contrast to FFGR, is sustained by coherent dynamics. We note in passing that such parametric instabilities can only occur in bosonic systems [or, more precisely, in arbitrary systems with bosonic elementary excitations].
As described in the main text, the excitation peaks are evaluated from difference images (Fig. 2b), which are obtained by subtracting the mean condensate profile at t = 0 from each single image at time t > 0 for a certain set of modulation parameters. The mean image at t = 0 is obtained as described below. Due to the high density of the Bose-Einstein condensate (BEC) in the center, the imaging light is completely absorbed there, leading to a saturation of the pixel values, especially at t = 0, when the number of atoms in the BEC is maximal.
In order to minimize statistical fluctuations in the mean position of the atom cloud, we shift the pixels such that the center of the BEC peak is at (x, y) = (0, 0). The center is extracted by fitting a two-dimensional (2D) Gaussian function to the BEC peak, which is truncated at some value I 0 to account for the saturation in the center: For each parameter set, we determine a saturation value I 0 from all images at t = 0 and apply it to all images in this data set. We start by fitting a truncated Gaussian to each image at t = 0 separately and take the mean over all I j 0 as a starting value for the following fit. Then, we fit a truncated Gaussian to all images simultaneously with the constraint that the saturation value I 0 has to be the same in all images, whereas all other parameters can vary. This procedure results in the final value I * 0 . To determine the center positions, we fit a truncated Gaussian g * (x, y) with saturation value I * 0 to all images and use these to shift the pixels in each image such that all center positions are at coordinates (0, 0) and truncate the pixel counts on all images at I * 0 . The centered and truncated images are then averaged to give a mean t = 0 image for each parameter set.
For the subtraction, we take each single image at every t > 0 and first fit the previously defined truncated Gaussian g * (x, y) with the fixed I * 0 to determine the center of the BEC peak. Then, the profiles are again shifted to move their center to (0, 0) and the pixel values are truncated at I * 0 . From the centered, truncated image we then subtract the mean t = 0 image in order to obtain a difference image as shown in Fig. 2b in the main text. Due to the centering of the images the statistical fluctuations due to drifts of the BEC position are suppressed. Hence, fluctuations of the position of the excitation peaks that are detected for t > 0 correspond to the uncertainty of this peak position relative to the BEC position and is reflected in the standard error of the corresponding mean value at this time.
The conversion from pixel values to quasimomentum is achieved by taking an absorption image after switching of the lattice abruptly and a time-of-flight (TOF) of 6 ms. This results in interference peaks, which are separated by 2π/d. We determine the positions of these by Gaussian fits and average over several experimental realizations, which gives the width of the first Brillouin zone (BZ) in pixels. From the positions of the Bragg peaks, also the exact direction of the lattice on the camera can be determined. In Fig. 2, the images were rotated by 6 • to match the lattice axis with the horizontal direction for better visualisation. In the data analysis, we sum up the pixels along an axis (or perpendicular to) that is turned by 6 • instead of rotating the images, which would lead to false empty pixels.
As discussed in the main text, we expect that for modulation frequencies above the saturation frequency ω sat , the excitation peaks in momentum space should be located at π/d with a tail towards smaller quasimomenta, resulting from excitations to other modes and from the width of the condensate at t = 0, which has a finite (Gaussian) width of σ 0 ≈ 0.2π/d. As shown in Fig. 4, the measured saturation value of the most unstable mode along the lattice axis q x mum is slightly smaller than π/d. This is most likely due to the short TOF used in the experiments, i.e., we measure a convolution of the momentum profile with the insitu density distribution of the BEC. This leads to a slight shift of the peak maxima in the TOF images towards the center of the BZ.

B. Details of the Peak analysis
To analyze the peaks along the lattice direction, we integrate all difference images (Fig. 2b) along the direction perpendicular to the lattice to obtain a 1D profile, as described in the main text. To extract the peak amplitude A x and peak position we fit an asymmetric Loretzian function where we define the peak amplitude as A x =Ã − |A 0 |. The asymmetry of the profiles arises from the thermal background inside the first BZ which effectively increases the profile height at the inner side of the peaks. To determine the transverse width ∆q y of the peaks, we integrate the difference profile along the lattice direction inside a region of interest that contains only the left or right peak (excluding the negative pixels in the center) and fit a symmetric Lorentzian to each of the peaks: L(q y ) =B (∆q y ) 2 (∆q y ) 2 + 4(q y max − q y ) 2 + B 0 .

(S.2)
The end of the short-time regime t s (Fig. 3) is obtained by fitting the following piecewise function to the logarithm of the peak amplitudes ln(A x ) for the left and right peak independently: The final value for t s is defined as the average over the fit results for the left and right peak (rounded to integer multiples of the modulation period).
We find that for parameters, where q x mum saturates at the band edge, q x ≈ π/d, the observed momentum profiles are slightly asymmetric, i.e. either the left or right excitation peak has a larger amplitude (Fig. S1). We further observed that the asymmetry depends on the initial phase of the drive ϕ. It is known that for momentum components at the band edge any residual gradients or finite ramp times during the switching off of the lattices can result in an asymmetry in the bandmapped images. In order to reduce this systematic deviations we average over individual experimental realizations with ϕ = 0 and ϕ = π. The final results show no residual systematic asymmetry (Figs S2-S4).

C. Atom number calibration
To determine the atom number, we measure the trapping frequencies in the absence of the lattice [ω r /(2π) = 23.2(3) Hz,ω z /(2π) = 189(3) Hz] and the Thomas-Fermi radii as a function of the scattering length a s and fit the number of atoms according to the Thomas-Fermi prediction .

(S.4)
We obtain N = 3.7(4) × 10 5 , where the error bar mainly stems from atom number fluctuations. Note, that for the atom number calibration, the harmonic confinement is reduced compared to the measurements with the lattice that are presented in the main text. The trapping frequencies are determined by monitoring the center of mass motion of the BEC in-situ after a displacement in the trap. This gives the harmonic trapping frequencies along the axes of the dipole trap beams. In the horizontal plane the trap is rather symmetric withω x /(2π) = 24.0(5) Hz andω y /(2π) = 22.3(3) Hz, so we consider their mean valueω r . The same is true for the situation with the lattice being on.

S2. PEAK PARAMETERS VS. MODULATION TIME
In Fig. S2 the logarithm of the fitted peak amplitudes versus modulation time is shown for all modulation frequencies and modulation amplitudes α = 1.44 and α = 1.78 together with the fitted piecewise functions [defined in Eq. (S. 3)]. The dashed lines mark again the transition time t s where the short time regime ends, which is the mean value of the fitted kink positions for the left and right peak rounded to integer multiples of the modulation period. In general, the peak amplitudes decrease with the modulation frequency, also the errorbars become smaller. We attribute this to the enhanced influence of incoherent processes at small modulation frequencies, which are closer to the bandwidth of the static 1D lattice (4J/h = 0.43 kHz). The fitted peak positions q x max are shown in Fig. S3 where the black dashed lines again mark the end of the short-time regime, t s . The black lines denote q x mum which is the mean over all positions up to t s . The data points in the upper panel of Fig. 4 are the average over q x mum for the left and right peak for each set of modulation parameters. In Fig. S4 we plot the transverse width ∆q y versus modulation time, the grey shaded area denotes the data range, which is averaged to obtain the data points in the lower panel of Fig. 4. In general, the transverse width is nearly constant for frequencies below the saturation frequency, as discussed in the main text. For frequencies larger than 750 Hz for α = 1.44 and 625 Hz for α = 1.78, the peak widths increase but they also fluctuate more as indicated by the larger error bars.    Figure S4. Transverse full-width at half maximum (FWHM) versus modulation time. Each point is an average over ∼ 10 individual experimental realizations, < ∼ 5 realizations for a modulation phase ϕ = 0 and < ∼ 5 with ϕ = π. The error bars indicate the standard error. The black, dashed line marks the end of the short-time regime ts and the grey shaded area extends from ts − 1 to ts + 1. The points in this range were averaged to obtain the data points in Fig. 4 (lower panel).

S3. THEORETICAL METHODS AND ANALYTICAL TREATMENT OF INSTABILITIES IN THE PERIODICALLY-DRIVEN BOSE-HUBBARD MODEL
In this section, we give details about the theoretical methods and results used in the main text to discuss the dynamics and stability properties of the periodicallydriven weakly-interacting Bose-Hubbard model (BHM). For simplicity we set = 1, the lattice constant d = 1 and d ⊥ = 1 in the entire section, cf. Sec. S5.
The BHM is defined by the Hamiltonian whereâ i,r ⊥ annihilates a particle at lattice site i and transverse position r ⊥ = l ⊥ d ⊥ê⊥ . J > 0 denotes the tunneling amplitude of nearest-neighbor hopping along 2mn j,q ⊥ describes a free-particle kinetic part along transverse directions, and U > 0 is the repulsive on-site interaction strength. The time-periodic modulation has amplitude K, phase ϕ, and frequency ω = 2π/T with T the driving period.
Transforming the system to the rotating frame [1], we eliminate the external shaking at the expense of introducing a periodically-modulated hopping term with the dimensionless coupling strength α = K/ω. In momentum space, the free-system [U = 0] is governed by the time-dependent dispersion oscillating at the driving frequency ω. Our study relies on three different methods: In the weakly-interacting regime we apply Bogoliubov theory. Weakly-interacting bosons at ultracold temperatures form a Bose-Einstein condensate (BEC) which is macroscopically occupied, whereas the population of all other modes remains small compared to that of the condensate. This scale separation motivates the Bogoliubov approximation: where N BEC (t) denotes the BEC population, q BEC is the mode in which condensation occurs, andb † q creates a particle of momentum q on top of the BEC background.
Plugging this ansatz into the Hamiltonian (S.2), and keeping terms to quadratic order in the operatorsb q , we arrive at the time-dependent Bogoliubov Hamiltonian (S.4) Here g = U n, n = N BEC (t = 0)/V with V the volume of the system. The time-dependent shift ε qBEC (t) arises from the chemical potential µ(t) = ε qBEC (t) + g which is fixed by the condition that the linear term (inb q ,b † q ) in the Bogoliubov expansion vanishes. The same result can be obtained from a linearization of the time-dependent Gross-Pitaevskii equation around the BEC wavefunction [1,2].
In the Bogoliubov approximation, the dynamics of the system is governed by the time-dependent Bogoliubov mode functions u q (t) and v q (t), defined througĥ , we arrive at the Bogoliubov-de Gennes (BdG) equations: The BdG equations (S.5) describe the quasi-particle dynamics, and include micromotion effects.
Due to the time-periodicity of Eq. (S.5), Floquet theory applies. We can thus focus on the stroboscopic funda-mental matrix Φ(T ), which is obtained by time-evolving Eq. (S.5) over a single period T . We denote the eigenvalues of Φ(T ) by q . The appearance of eigenvalues with positive imaginary parts indicates a dynamical instability, i.e. an exponential growth of the corresponding mode(s), characterized by the rate Im q .
We then define the instability rate as the maximum growth rate: It is independent of the reference frame, and governs the stroboscopic dynamics of the mode functions u q and v q . Notice that the number of exited atoms, n(t) ∼ q |v q (t)| 2 , has an instability rate 2Γ. Let us also introduce the most unstable mode which dominates the long-time BdG dynamics because all other modes grow more slowly. In the main text, q mum is shown to provide clear experimental signatures of parametric instabilities. Details on the explicit derivation of the instability rate and the maximally-unstable mode for the BHM can be found, e.g., in Ref [1].

B. Weak Coupling Conserving Approximation (WCCA)
The Bogoliubov approximation holds under the condition that the condensate occupation remains large throughout the duration of interest. A prerequisite for this is the presence of weak enough interaction strength U/J 1, which renders the occurrence of quasiparticle collisions rare. Indeed, the Bogoliubov approximation works well when the latter do not play an important role, and thus can be neglected. This is the case, for instance, for weakly interacting bosons in equilibrium.
Away from equilibrium, excitations may form induced by the drive, leading to fast depletion of the macroscopically occupied condensate, even though the interactions may still remain small (compared to the bare hopping) at all times. The condensate depletion is enhanced when resonant transitions between the condensate mode and higher-energy states occur; such transitions can be stimulated by a periodic drive.
Conceptually, one of the major drawbacks of the timedependent Bogoliubov approximation is the lack of particle number conservation [3]. This means that, even though it may capture the onset of condensate depletion, the approximation is doomed to fail in correctly predicting the dynamics of the quasiparticle distribution during the later stages of the evolution. One reason for this lies in the fact that the BdG EOM assume a constant in time, and thus endless, reservoir of particles provided by the macroscopically occupied condensate mode.
An elegant way to restore particle number conservation is provided by the Weak Coupling Conserving Approximation (WCCA). It represents a minimal extension of the linearized Bogoliubov equations which conserves the corresponding global U (1) symmetry at all times. The WCCA EOM is a set of coupled differential equations for the time-dependence of the condensate wavefunction φ(t) and the quasiparticle correlators F 11 and F 12 , defined as The WCCA equations can be derived by Legendretransforming the action of the BHM on the Keldysh contour with respect to both the order parameter (condensate field) φ, and the quasiparticle propagators F 11 and F 12 [4,5]. The resulting effective action, which contains all two-particle irreducible diagrams, is then expanded to leading order in the interaction strength U . Minimizing this effective action w.r.t. the condensate field and the propagators, one arrives at the following system of coupled integro-differential equations [5]: where we denoted by * the complex conjugation. The initial condition is given by the Bogoliubov ground state, together with |φ(t = 0)| 2 = N BEC (0). Notice that the presence of the chemical potential in the WCCA EOM is irrelevant to any observables, due to the U (1)symmetry being conserved at all times, in contrast to the BdG EOM, for which µ(t) is crucial to recover the correct spectrum of excitations. Hence, for the sake of a proper comparison with the BdG approximation, we may set µ(t) = ε qBEC (t) + g, which also fixes the reference frame.
One can convince oneself that, if one neglects all terms involving q -summations in Eqs. (S.9), the equation for the condensate wavefunction φ(t) decouples from the other two. Moreover, one may further recognize it as the Gross-Pitaevskii equation in the presence of the periodic drive. Using the definition (S.8), it follows that F 11 (t; q) = 1/2(|u q (t)| 2 + |v q (t)| 2 ) and F 12 (t; q) = u q (t)v q (t), and hence the remaining equations for F 11 and F 12 are equivalent to the BdG EOM from Eq. (S.5). Going back to the complete equations, we see that the qsummations represent the essential coupling between the condensate and the quasiparticle excitations, necessary to restore particle number conservation to leading order in U . The conserved quantity itself is the total number of condensed and excited atoms Despite the advantages of the WCCA EOM over the BdG equations discussed above, let us make some remarks about the applicability and usefulness of the WCCA. First, the WCCA EOM are not amenable to simple analytical treatment, due to their non-locality in momentum space. Second, unless condensate depletion is suppressed, e.g. by the existence of a pre-thermal phase [5], the WCCA is also expected to be valid only in the short-time limit, since it misses the collisions between quasiparticles which appear first at order O(U 2 ). Thus, within the WCCA, a periodically driven system cannot thermalize which points out that the WCCA does not capture the later stages of the heating process. Going to order U 2 has been done recently in the large-N limit of the periodically-driven O(N ) model [6,7], and the existence of long-lived pre-thermal plateaus has been shown in the limit where the drive frequency is the largest energy scale in the problem.

C. Truncated Wigner Approximation (TWA)
Since the exact quantum equations of motion for the BHM are often too hard to analyze, it may be advantageous to perform a semi-classical analysis. This proves particularly useful when the initial condition is deep into the superfluid phase, as the latter is well captured by the GPE, a classical nonlinear wave equation. Recently, it has been demonstrated that thermalization in classical Floquet systems behaves very similarly to their quantum counterparts [8][9][10].
An systematic way of deriving the leading-order quantum corrections in the semiclassical limit, is to apply the Truncated Wigner Approximation (TWA). This method has the advantage that it is both particle-number conserving [since the GPE conserves U (1) symmetry] and, at the same time, it is capable of capturing the thermalizing dynamics induced by the periodic modulation. Since there is extensive literature on the TWA already, below we briefly summarize the main ideas useful to our analysis, and refer the interested reader to Refs. [11,12].
From translational invariance, it follows that the lowest-energy state of the BHM must necessarily be uniform in real space and, therefore, in momentum space the entire weight is carried by the q BEC = 0 mode. Thus, in the classical limit, the uniform condensate field obeys the following time-dependent GPE: Similar to the WCCA EOM, the presence of the chemical potential in the GPE is irrelevant to any observables, due to the U (1)-symmetry of the GPE. For instance |φ(t)| 2 = N BEC ≈ N at all times. In quantum mechanics, however, a small but finite fraction of atoms always occupies the excited modes of the system, even for weak interactions, known as quantum depletion. Let us again denote the Bogoliubov mode functions by u q (t) and v q (t), and recall that the complete atomic field expansion readŝ In the TWA [11], the occupation of momentum mode q is described by a set of independent identically distributed complex-valued Gaussian random variables γ q , analogous to the quantum mechanical operatorsγ q , such that the semiclassical condensate field is decomposed as (S.12) Since any Gaussian distribution is uniquely determined by its only two non-vanishing moments, we can choose the mean and variance of γ q to correctly recover the true quantum mechanical fluctuations within the quadratic Bogoliubov theory [11].
We can now use the TWA ansatz (S.12) as an initial condition for the GPE. Notice that this ansatz wavefunction is no longer uniform in space, due to the random contribution from the excited modes, and thus the time evolution needs to be computed in real space. Thus, within the TWA, the uniform condensate wavefunction is seeded with the correct quantum fluctuations, which ensures that the dynamics is captured correctly at the semiclassical level.
To compute expectation values of observables, we produce an ensemble of initial states a j , time-evolve each state separately up to the time of interest, compute the observable for each realisation a j (t) separately, and perform an ensemble average in the very end. For example, the momentum distribution of the excitations within the TWA is given by (S. 13) where (·) stands for averaging over the random variable γ q .
While the TWA captures the effects of quasiparticle collisions to some extent, it should be noted that this formulation of quantum dynamics is not exact. One of the reasons for this is that it uses Gaussian distributed random variables, whereas in general the true quantum distribution from the BHM has nonzero higher-order moments beyond the Bogoliubov approximation.

A. Homogeneous 2D System
While it is possible to simulate a 3D homogeneous system, the strong harmonic trap along the z-axis present in the experiment effectively confines the system to a quasi-2D geometry. We therefore model the experiment by one lattice and one transverse degree of freedom. While the squeezed continuous z direction is irrelevant for the parametric instability effects within BdG theory, effects of the additional z-confinement are expected to be felt in the TWA simulations only at later times, beyond the regime of agreement between BdG and TWA (see Fig. 6 in the main text).
For the BdG and WCCA simulations, we use a condensate density of n = 20 to initialize the system in the Bogoliubov ground state. To determine the condensate wavefunction for the TWA, we use imaginary time evolution to find the lowest-energy state of the GPE in the presence of weak interactions, and consider an ensemble of 10 3 TWA realizations, performing an additional averaging over the phase of the drive ϕ to reduce the effects of the initial kick.
For all simulations, the momentum in the x direction (lattice) is discretized according to the quantization condition q x = 2πl/(N x sites d) with l ∈ [1, 2, . . . , N x sites ]. To simulate numerically the y (continuous) direction we apply the following auxiliary discretization: q ⊥ = 2πl ⊥ /(N ⊥ sites d ⊥ ) with l ⊥ ∈ [1, 2, . . . , N ⊥ sites ], imposing an upper cutoff on the values of perpendicular momentum at 2π/d ⊥ . We checked that increasing this cutoff does not change the results, since high-momentum states have high energy and remain unaffected by the dynamics.
The volume of the finite 2D system system is V = N x sites N ⊥ sites , where N x sites and N ⊥ sites are the number of sites along the respective axis. Setting N x sites = 81 and N ⊥ sites = 80 modes, we find that the perpendicular momentum grid is dense enough for the system to be in the thermodynamic limit already at d ⊥ = d = 1. We further make sure the results remain unchanged with increasing the number of modes.
Effects due to the weak confining harmonic trap are not considered in these simulations, and are addressed separately [cf. Appendix B]. We use QuSpin [13,14] to simulate the dynamics.

B. Inhomogeneous 1D system
The computational complexity of simulations in the presence of harmonic confinement limits the reachable system sizes. Unlike homogeneous systems, the inhomogeneity induced by the trap leads to a break-down of momentum conservation and couples all momentum modes. Hence, to determine the Bogoliubov mode functions, needed for the initial condition for the EOM, one has to exactly diagonalize a matrix of size N sites × N sites , with N sites the total number of space points. It is the size of this matrix, together with the additional bottlenecks for solving the corresponding equations of motion, which set the current limit on the reachable system sizes. For these reasons, simulating the dynamics of a 3D system in the thermodynamic limit (N sites ∼ 10 6 ) in the presence of a trap and subject to periodic modulation, is not feasible without introducing further approximations. We therefore limit the theoretical discussion of effects due to the harmonic trap to 1D lattices.
To obtain the numerical data in the presence of a harmonic trap potential, we simulate the BdG equations for a 1D lattice of L x /d = 201 sites [15]. We first perform imaginary time evolution in the presence of the trap to determine the exact initial condensate profile. This inhomogeneous profile results in an inhomogeneous effective interaction strength g(x), which is used to compute the Bogoliubov modes in the absence of the periodic drive used as an intial condition [16]. Finally, the timeperiodic BdG dynamics is simulated by solving the realspace variant of Eq. (S.5) with the additional replacement g → g(x).