Phonon excitations of Floquet-driven superﬂuids in a tilted optical lattice

Tilted lattice potentials with periodic driving play a crucial role in the study of artiﬁcial gauge ﬁelds and topological phases with ultracold quantum gases. However, driving-induced heating and the growth of phonon modes restrict their use for probing interacting many-body states. Here, we experimentally investigate phonon modes and interaction-driven instabilities of superﬂuids in the lowest band of a shaken optical lattice. We identify stable and unstable parameter regions and provide a general resonance condition. In contrast to the high-frequency approximation of a Floquet description, we use the micromotion of the superﬂuids to analyze the growth of phonon modes from slow-to fast-driving frequencies. The model describes phonon excitations in both resonantly and nonresonantly driven systems, with or without a tilted potential. Our observations enable the prediction of stable parameter regimes for quantum-simulation experiments aimed at studying driven systems with strong interactions over extended time scales.

Especially, the combination of periodic driving with a constant force that tilts the lattice potential has been instrumental for the implementation of artificial gauge fields.The constant force suppresses tunneling between lattice sites, while other mechanisms, such as laser-assisted tunneling [17,18] or near-resonant driving [19,20], reintroduce the coupling between lattice sites with the desired properties.These so-called Floquet-engineered lattice potentials have been successfully applied to investigate single-particle effects or weakly interacting quantum gases.Simulating many-body states of interacting particles, however, remains challenging as interactions create instabilities and heating on short time scales comparable with the modulation period and quickly destroy the coherence of the system [21,22].Developing an understanding of these instabilities and finding an optimal window for the driving frequency [23] are crucial to include interactions in quantum-simulation experiments with periodic driving.
In this article, we experimentally and theoretically investigate instabilities due to the spontaneous formation of phonon modes in a tilted lattice potential.We examine the role of interactions for the time evolution of a BEC of cesium atoms confined in a one-dimensional (1D) lattice potential that is periodically shaken and subject to a constant force F B [Fig. 1(a)].In such a system, phonon modes can grow exponentially in time and eventually destroy the BEC.To understand the origin of this phenomenon, often referred to as Floquet heating [23], we experimentally probe the time evolution of the many-body ground state of the driven system and measure its stability with respect to driving strength and frequency.Our findings are used to develop a model that predicts the stability of the system based on resonant phonon excitations.
In addition to phonon instabilities, the system exhibits a multitude of other driving resonances [Fig.1(b)].Resonant excitations to higher lattice bands occur when the driving frequency matches the energy gap between the bands or fractions thereof [24,25].In analogy to electrons in solids that are driven by electric fields, these higher-order excitations are often referred to as multiphotonlike [26].Furthermore, tunneling resonances between neighboring lattice sites occur when the driving frequency matches the energy shift between sites.Multiples of this driving frequency couple sites at further distances, while fractions of it allow for multiphoton-assisted tunneling [19,20,27].For tilted potentials, phonon instabilities coincide with these interband and tunneling resonances and add another manifold of higher-order driving resonances to the system.Finding a frequency window free of resonances is experimentally challenging and requires a detailed understanding of the underlying mechanisms [23].Within the lowest lattice band, periodically driven quantum gases show an exponential growth of phonon modes due to modulational and parametric instabilities.Parametric instabilities are caused by oscillating system parameters [28], while modulational instabilities result from properties of the medium that exist even without periodic driving, such as attractive interactions or a negative effective mass for certain momenta of the superfluid in the lattice [29,30].They also occur in driven systems when the driving force accelerates the medium on a micromotion in momentum space.In a previous study [31], we demonstrated the existence of parametric and modulational instabilities in driven superfluids independently and studied their properties.However, in the presence of a tilted potential, the micromotion of the superfluid always covers the complete first Brillouin zone, and it is experimentally challenging to separate both mechanisms.In this article, we aim to explain the origin of Floquet heating due to resonant phonon excitations that occur when driving frequencies are comparable with phonon energies.
Floquet theories provide a very successful approach to analyze quantum gases in the fast-driving limit when driving frequencies exceed all other frequencies of the system.In a high-frequency approximation, the system is well-described stroboscopically using a time-independent Floquet Hamiltonian [2].Recently, this description was extended to lower frequencies to analyze parametric resonances [28,[32][33][34][35].We were not able to map our experimental data to those models.Instead, we took the opposite approach and extended the description of phonons and modulational instabilities in nondriven systems to the case of resonant driving frequencies.Our approach offers an intuitive explanation of our measurement results, and the description remains valid in the limits of both slow and fast driving.For intermediate driving frequencies, resonant excitations are caused by the periodic growth of phonon modes whenever the micromotion passes through a modulationally unstable region of the Brillouin zone.Band excitations, which occur for very large driving frequencies [25], and trap excitations, which require long observation times [36], have little impact on our results, and we omit them in our discussion.
We experimentally study phonon modes in three settings: nondriven with a constant force F B (Sec.III), resonantly driven with a frequency that matches the energy shift between adjacent lattice sites hν D = F B d L , at distance d L (Sec.IV), and off-resonantly driven with detuning h ν = F B d L − hν D (Sec.V).For all settings, the superfluid cycles through unstable regions of the first Brillouin zone which causes the growth of excitation modes.As a result, resonant growth of excitations does not depend directly on the driving frequency ν D but on the frequency ν c of those crossings into critical regions.Time-averaging, e.g., to calculate the time-averaged Bogoliubov energy of a phonon mode within a Floquet description, results in the loss of information about the shape of the micromotion.Instead, we use the number of crossings and the time intervals between them to determine ν c and the resonance condition.

II. EXPERIMENTAL SETUP
Our starting point for the experiment is a BEC of ∼50 000 cesium (Cs) atoms in a vertical optical lattice potential [Fig.1(c)] that was created using two counterpropagating laser beams with wavelength λ = 1064 nm, lattice constant d L = λ/2, and lattice momentum k L = π/d L [37].Typical lattice depths V were between 2 and 14 E r , where E r is the recoil energy of Cs atoms.The BEC was confined by a cross-beam optical dipole trap with trapping frequencies ω x,y,z = 2π × (11,18,14) Hz and levitated by a magnetic field gradient [38,39].We controlled the atomic interaction strength by tuning the s-wave scattering length a s with a magnetic Feshbach resonance before loading the wave packet into the lattice potential in 150 ms.The values of a s and V were varied for different measurements to provide experimentally convenient time scales for the growth of excitation modes.
We created the constant force F B and the driving force F (t ) with different experimental techniques.The constant force was applied by reducing the magnetic levitation in 1.5 ms to a fraction of the gravitational force, which results in Bloch oscillations with period T B = h/(F B d L ) and frequency ν B = 1/T B that are directly observable in the superfluid momentum distribution in the laboratory frame.The driving force F (t ) = F 0 cos(2πν D t ) was applied by periodically shifting the sites of the lattice using two acoustic-optical modulators that create a frequency difference between the laser beams [40,41].This method provides precise control over driving frequency and strength.However, F (t ) is an inertial force in the lattice frame, and the resulting micromotion is only indirectly detectable in the laboratory frame by measuring the total momentum of the superfluid or the weight of the reciprocal lattice peaks in momentum profiles [42].By combining the two techniques, we were able to directly measure the Bloch period in the laboratory frame while facilitating fast-driving frequencies in the kilohertz regime.
After loading the superfluid into the lattice potential, we applied F B and F (t ) for a time t that was adjusted to the closest multiple of the driving period to reduce effects of the micromotion on the final momentum distribution.To minimize the impact of the trapping potential, which determines the state of the superfluid on long time scales [36], we kept t smaller than the trap period and studied phonon excitations along the lattice direction.The lattice potential was switched off instantly or ramped off in 1.2 ms to determine the real momentum or the quasimomentum distribution with absorption imaging after an expansion time of typically 75 ms.For measurements with a clearly identifiable carrier wave packet, we applied Gaussian fits to determine the atom number in the carrier wave N C and in excitation modes N E = N tot − N C , where N tot is the total atom number.This approach was challenging for measurements with a significant fraction of atoms in excitation modes, and we instead counted the number of atoms in fixed momentum intervals where we expected the phonon modes or the carrier wave (see Appendix C).

III. EXCITATIONS WITHOUT PERIODIC DRIVING
To demonstrate the main concepts, we first studied excitation modes for only a constant force F B .Without driving, the superfluid Bloch oscillates in the first Brillouin zone with its micromotion k(t ) following straight lines in quasimomentum space [41].Absorption images, showing the momentum distribution in quasimomentum space [Fig.2(a)], were used to determine the fraction of atoms in excitation modes N E /N tot [blue circles in Fig. 2(b)].Whenever the superfluid crosses into critical regions of the Brillouin zone with negative effective mass |k(t )| > k L /2, small perturbations of the stationary state grow exponentially in time [30,43], while excitations are steady or decay in stable regions [31] (see Appendix A).
We quantified the stability of the system for a varying constant force and lattice depth by measuring N E /N tot after approximately t = 20 ms.Yellow and blue colors in the (V, T B )-stability diagram [Fig.2(c)] indicate stable and unstable regions with weak and strong growth of phonon modes, respectively.We expect the superfluid to become stable in deep lattices and for large values of F B due to the suppression of tunneling when the energy gap between neighboring sites F B d L exceeds the width of the lattice band 4J, with tunneling matrix element J.However, we find that the decoupling of lattice sites is far less abrupt for a 1D superfluid with large atom numbers per site than for single atoms [44], and the superfluid remains phase coherent beyond F B d L = 4J [dashed red line in Fig. 2(c)].This allows us to observe the reduction of modulational instabilities when the system becomes stable in the fast-cycling regime.This onset of stability shifts for increasing interaction strength toward smaller values of T B .
For our analysis, we assume that phonon modes with momentum hq are resonantly excited when the frequency ν c of the carrier wave crossing into critical regions of the Brillouin zone matches twice the time-averaged Bogoliubov energy of the mode E q .This resonance condition for modulational instabilities is similar to the condition for parametric instabilities hν D = 2 E q in Refs.[28,33,35,45], but it relies on the frequency ν c .For strong driving forces, several crossings of critical regions can occur per Bloch period, and we introduce the parameter α to relate the resulting frequencies ν c to ν B .The parameter α provides the fraction of time between two consecutive crossings per Bloch period with ν c = ν B /α.As a result, our resonance condition predicts a series of values T B or ν B with strong growth of excitation modes for where m p is an integer that indicates higher-order phonon resonances (see Sec. IV).
The time-averaged energy of a phonon mode E q [k(t )] with micromotion k(t ) can be approximated by two methods, both of which are based on the Bogoliubov-de Gennes (BdG) equations [28,33,46,47].Please see Eq. (A5) for more information.Time-averaging the single-particle energy in the BdG equations provides an approximation E f q (K, k 0 ) that depends on the driving strength K and on the initial momentum k 0 = k(0) [33,43,46,48], see Appendix A. For fast-driving frequencies, this analytic expression maps directly to the energy of a phonon E q (k) in a nondriven system when we identify k with k 0 and the tunneling matrix element J with J eff (K ) (see Appendix A).This mapping indicates that modulational instabilities also exist in periodically driven systems for |k 0 | > 0.5k L .However, the approximation assumes fast-driving frequencies and neglects the short time intervals of phonon growth when the micromotion cycles through the Brillouin zone.
For a low-frequency approximation, we numerically timeaverage the energy E q [k(t )] of an existing phonon over k(t ) [31]: Due to modulational instabilities, the phonon energy E q (k) is a complex number with an imaginary component that indicates the growth rate of the mode q and a real component that provides the frequency of the phase oscillations of the mode [49,50].The same complex description applies to E s q but not for the approximation E f q .Only the real component of E s q is used in Eq. ( 1) to predict resonant cycling frequencies.Without driving, the micromotion crosses once per Bloch period into critical regions with α = 1 [Fig.2(a)], resulting in a periodic increase of the phonon mode occupation [Fig.2(b)].Resonances occur when the frequency ν c of these crossings matches the natural oscillation frequency of the phonon density 2 E q (with m p = 1).For faster frequencies, the system becomes stable when hν c exceeds the energy of all phonon resonances, i.e., when the cycling due to Bloch oscillations becomes faster than possible response times of the phonons.We find that the transition to this fast-cycling regime (or fast-driving regime in later sections) is well approximated by Eq. ( 1) [solid red line in Fig. 2(c) with α = m p = 1 in Eq. (1)].An even better prediction for the transition from an unstable to a stable system is provided by T B = h/(2.96√ 2JU ) in Ref. [43] [black line in Fig. 2(c)], which approximates a numerical solution of the BdG equations.The equation describes the transition line between stable and unstable parameter regions with real and complex phonon energies.
We also observed a second stable region at shallow lattice depths V < 2 E r [Fig.2(c)].For repulsive interactions, modulational instabilities are induced by the lattice potential, and we expect the superfluid to become stable in the limit of vanishing lattice depths.We omit the discussion of instabilities in shallow lattice potentials, as identifying unstable intervals of the micromotion is more challenging beyond the tight-binding regime [49].
Gray patches in Fig. 2(c) indicate strong atom loss due to excitations to higher bands.Band excitations can be driven by Landau-Zener tunneling at the edge of the Brillouin zone, which typically occurs for small lattice depth and large accelerations, or by resonant excitations when the driving frequency matches the energy of the band gap.Except for Fig. 2(c), we did not observe excitations to higher bands in our experimental data, and we focus our discussion in the following sections on the lowest band.

A. Experimental measurements
In a second series of experiments, we added a resonant driving force F (t ) with a driving period that matches T B .Despite a strongly tilted lattice potential, resonant driving reintroduces coupling between the lattice sites with an ef-  1) when approximating the phonon energy by E s q and E f q .
fective tunneling matrix element J eff (K ) = JJ 1 (K ) [51].Here, J 1 (K ) is the first-order Bessel function, and is the dimensionless driving strength.For increasing driving strength, the straight lines of the micromotion [K = 0, dotted black line in Fig. 3(a)] start bending upward or downward [red and blue lines in Fig. 3(a)], depending on the initial directions of F B and F (0). Thick lines in Fig. 3(a) indicate time intervals during which the micromotion of the superfluid crosses into critical regions of the Brillouin zone.We used opposite initial directions of the forces to reduce these time intervals [red and orange lines in Fig. 3(a)], as identical initial directions extend them (blue dashed line).
To develop a better understanding of the stability of the system, we calculated the growth rate q (T B , K ) of phonons with momentum q by integrating and diagonalizing the BdG equations [33,46,48], see Appendix A. The resulting (q, T B )-stability diagram for K = 2 [Fig.3(b)] shows a sequence of unstable regions with strong growth (blue color) which belong to fundamental and higher-order phonon resonances m p = 1-7.The resonance condition in Eq. ( 1) provides a good approximation of these resonances (red lines with α = 1).Successive higher-order resonances are separated by a stable region (yellow color).Position, width, and growth rate for these stable and unstable regions can be related to the time that the micromotion spends in critical regions of the Brillouin zone (see Appendix B).However, this q dependence of is usually not directly visible in experiments due to the stochastic nature of modulational instabilities.Phonon modes require small initial modulations of the medium to start growing [43].Those seeds are of random nature, caused by thermal, quantum, or technical fluctuations, which result in the occurrence of modes with varying values of q [52].
We experimentally determined the growth of phonon modes by measuring N E /N tot after a driving period of ∼30 ms.Two different initial quasimomenta, hk 0 = 0 and hk L , were used to prepare the superfluid in the ground state of the timeaveraged lattice band [panels in Fig. 3(c)].The micromotion of the ground state starts with k 0 = 0 for K < 3.83 when J eff (K ) is positive and with k 0 = k L for 3.83 < K < 7.02 for negative values of J eff (K ).To push the superfluid to the edge of the Brillouin zone k 0 = k L , we applied a weak magnetic field gradient before starting the drive (see Appendix C).
We compare the measured (T B , K)-stability diagram in Fig. 3(c) to the calculated growth rate (T B , K ) in Fig. 3(d).The statistical fluctuations of q are included by averaging (T B , K ) over 11 modes with equally spaced values of q = 0-k L .Our measurements and numerical calculations show good qualitative agreement for characteristic parameter regions that are indicated by Roman numerals in the diagrams, except for a residual structure of curved horizontal lines which are caused by higher-order phonon resonances.Our measurement data do not have the resolution along the T B axis and sufficient repetitions to resolve this structure.
Particularly important for future experiments is the stable region at high driving frequencies [region (i) in Fig. 3(d)] which, for weak driving strengths, starts approximately at the first resonance with m p = 1.Both E s q and E f k L provide a good approximation for E q in Eq. ( 1), except for K values close to the zeros of J 1 (K ) [red and white lines in Fig. 3(d)].Stronger driving strengths, e.g., 1.5 < K < 3.0, shorten the time that the superfluid spends in critical regions and reduce growth rates and widths of the resonances.As a result, even large driving periods show an increased stability, both in our data and in calculations [region (ii) in Figs.3(c) and 3(d)].Even stronger driving strengths, e.g., with K ≈ 3, cause three crossings of the micromotion into critical regions of the Brillouin zone [see also Fig. 4(b)].This reduces the α parameter and shifts the stable fast-driving regime toward larger values of T B [region (ii) in Figs.3(d)].Even stronger driving strengths, 3.9 < K < 5.9, for which the ground state has an initial momentum k 0 = k L , cause multiple crossings into critical regions with varying values of α.Despite the resulting complex structure of stable and unstable regions in our calculation [region (iv)], we find good agreement between calculations and experiment within our measurement resolution.

B. Interpretation of the stability diagram
The structure of stable and unstable regions in our calculation of can be explained using Eq. ( 1).For a direct comparison with Eq. (1), we omitted the averaging and calculated q (K, T B ) for a single mode q = k L and increasing values of K [Fig.4(a)].The local maxima of k L were determined using a peak finding algorithm (dashed black lines), and we compared them with the predicted resonances with parameters (m p , α) (solid red and blue lines in Fig. 4).
We find that unstable regions on the left side of each panel [regions (L) in Fig. 4(a)] show a mostly regular pattern of higher-order resonances.The vertical distance between two resonances increases for each panel with 1/α (blue lines) as predicted by Eq. ( 1).Here, α was determined using the micromotion for each value of K [Figs.4(b)-4(d) for each panel in Fig. 4(a)].The largest values of α [arrows in Figs.4(b)-4(d)] provide a good prediction for the resonance positions in regions (L).However, we find a more complex pattern of intersecting resonances on the right side of each panel in regions (R), and our prediction of the resonance position works less well [red lines in Fig. 4(a)].The difference between regions (L) and (R) is caused by the different number of crossings of the micromotion (colored patches), which have similar but slightly different time intervals between them in regions (R).We speculate that the structure of resonances is caused by this multitude of slightly different α values and their corresponding cycling frequencies ν c .Equation (1) provides a general condition for resonances in lattice systems with modulational instabilities and a cycling micromotion, irrespective of whether the cycling is caused by driving, a tilted lattice potential, or both.For systems without tilt, ν B must be replaced by the driving frequency ν D , and α increases from 0 for weak driving strength, when the micromotion never crosses into critical regions of the Brillouin zone, to 1  2 for moderate driving strength.The resonance condition for this special case with α = 1 2 and m p = 1 is hν D = E q , which matches the condition suggested in Ref. [28].

V. EXCITATIONS WITH NONRESONANT DRIVING
Finally, we examine the stability for off-resonant driving frequencies with detuning ν = ν B − ν D .Off-resonant driving in tilted lattice potentials can induce super-Bloch oscillations with large amplitudes in position space [20,53,54] that make the study of phonon modes challenging.To reduce the amplitude of those oscillations, we chose large detunings with ν > 20 Hz and a strong force with ν B = 1 kHz.This Bloch frequency also ensures that the system is well in the fast-driving regime, and any growth of phonon modes can be attributed to ν.In momentum space, super-Bloch oscillations show the characteristic straight lines of Bloch oscillations in the Brillouin zone (as in Fig. 2) when evaluating k(t ) stroboscopically at integer multiples of the driving frequency.The period of those oscillations is T = 1/ ν.
To understand the growth of phonon modes for nonresonant driving, it is helpful to study the complete micromotion k(t ).The shape of the micromotion gradually shifts, resulting in predominantly stable and predominantly unstable time intervals [red patches in Figs.5(a) and 5(b)], and it repeats itself after a period T tot = s/ν B = r/ν D for integer values s, r and rational ratios ν D /ν B .The period T tot depends on the choice of s, r and is less relevant for the phonon growth than the period between predominantly stable and unstable intervals T .We found that the growth of phonon modes shows a similar time dependence for super-Bloch oscillations as for Bloch oscillations in Sec.III.For instance, we observe for T = 10 ms that excitation modes grow in the predominately unstable time intervals at 5 and 15 ms [blue circles in Fig. 5(a)], and the system once again attains stability for larger detunings, e.g., for ν = 800 Hz [blue circles in Fig. 5(b)].
As in Sec.IV, we quantified the growth of phonon modes by measuring N E /N tot after ∼30 ms of driving.The resulting (K, ν)-stability diagram shows a stable region for zero detuning [yellow color in Fig. 5(c)], as expected in the fastdriving regime.Small values of ν cause instabilities as the superfluid spends longer time intervals in predominantly unstable regions (blue colors), while we find the system to be again stable for large detunings with | ν| > 360 Hz.The overall shape of the unstable region resembles an ellipse with a stable region in the form of a horizontal line at the center.
We compared our measurement results to a numerical calculation of the phonon growth rate q that was again based on the BdG equations.Instead of integrating the BdG equations over T B , as in Sec.IV, we used the complete driving cycle T tot .To match our experimental parameters and to provide rational values of ν B /ν D , we chose ν B = 1 kHz, T tot = sT B , and T D = (s/r)T B , with s = 500 and integer val- ues of r.We again used the phonon momentum q = k L that provides the strongest growth [Fig.6(a)].The shape of the unstable region and the stable line at resonance match our experimental results well.As for normal Bloch oscillations in Fig. 2(c), the transition lines between stable and unstable parameters regions are again well predicted by the approximation h| ν| = 2.96 √ 2J eff (K )U [43], where we replaced J by the renormalized J eff (K ) [black line in Fig. 6(a)].
To explain the shape of the unstable regions in Fig. 6(a), we use a resonance condition like the one in Eq. ( 1).We replace the Bloch frequency with the super-Bloch frequency ν and use α = 1 for Bloch oscillations: Due to the detuning, the averaged phonon energy continues to oscillate when we time-average over the Bloch period, and we instead approximate E q with a long time average E tot q calculated over T tot .For a direct comparison, we use the momentum q = k L instead of averaging over different values of q.Equation (3) provides a good approximation for the resonances [blue regions and dashed red lines in Fig. 6(a)]; however, E tot q is independent of K, and our approach only predicts resonance positions close to the maximum of J eff (K ), i.e., at K ≈ 1.85.There, we recovered the full K dependence of the resonance position by scaling the black transition line to E tot k L , with Gray lines in Fig. 6(a) indicate the predicted positions of the phonon resonances with m p = 1, 2, 3.As in Fig. 3(d), we observe stable regions that form thin stable lines between neighboring resonances.For completeness, we also compared the interaction dependence of the growth rate for K = 1.85 with Eq. ( 3) [Fig.6(b)] and found excellent agreement.The lines in Fig. 6(b) refer to the same resonances as in Fig. 6(a).We extended our calculation of k L to larger driving strengths and detuning to demonstrate the connection between phonon resonances and tunneling resonances [Fig.6(d)].As in Figs.3(c) and 3(d), we used the momentum of the ground state of the superfluid in the fast-driving limit as the initial momentum k 0 .The driving frequency ν D is provided on the vertical axis instead of ν to allow for easy recognition of tunneling resonances.For reference, the unstable region m p = 1 in Fig. 6(a) matches the first region at ν D = 1 kHz in the top left corner of Fig. 6(d).The unstable region is repeated for increasing K values with a shape that is given by the Bessel function in E tot q (K ).Solid black lines indicate the transition lines between stable and unstable regions.
We found that phonon resonances also occur at other driving frequencies due to additional tunneling resonances [Fig.6(d)].In strongly tilted lattice potentials, phonon resonances can only exist on top of a tunneling resonance, as phonons typically spread over several lattice sites which requires tunneling and coherence between the sites.Periodic driving restores the tunneling in tilted potentials when the driving frequency is resonant to the energy shift hν B .For driving frequencies ν D = (n t /m t )ν B , the tunneling resonances are of order m t and couple lattice sites at a distance n t d L [19,20].Combining all contributions, the complete resonance condition for phonon resonances of order m p and for tunneling resonance (n t /m t ) is Higher-order tunneling resonances require the regularization of J with the m t -order Bessel function J m t (K ) [54] which must be included when approximating the average phonon energy with E tot q (K ).We indicate the predicted driving frequencies for transition lines in Fig. 6(d) with black lines and find good agreement with the unstable regions in our calculation.For the tunneling resonances with n t = 2, we used products of Bessel functions to provide a guide to the eye.
To experimentally demonstrate the existence of phonon resonances at the next tunneling resonance, we measured N E /N tot after 30 ms of driving with a constant value K = 1.9 [dashed gray line in Fig. 6(d)].We observe two regions with strong growth of phonon modes and stable center points at ν D = 1.0 and 0.5 kHz [data points with red and blue colors in Fig. 6(d)].The first region matches our previous measurement in Fig. 5(c) for the tunneling resonance m t = n t = 1, while the second region occurs at tunneling resonance m t = 2, n t = 1.For our measurement parameters in Figs.5(c) and 6(c), the unstable regions of phonon resonances start merging, and we used a smaller interaction strength of U = 4J and a larger lattice depth V = 10 E r for the calculations to clearly show disjunct regions with phonon growth.

VI. CONCLUSIONS
In conclusion, we have studied the growth of phonons in a superfluid with a tilted 1D lattice in three scenarios: nondriven, resonantly driven, and nonresonantly driven.To determine the phonon growth, we measured the momentum distribution of the system after a fixed hold time and obtained the fraction of atoms that were not in the ground state.For all settings, we found stable and unstable parameter regimes which can be explained by analyzing the micromotion of the superfluid through the first Brillouin zone.Due to the tilted potential, the micromotion always crosses into critical regions of the Brillouin zone, and modulational instabilities make the superfluid unstable for short time intervals.We found that the duration and multitude of those crossings per cycle determine the growth rate of phonon modes.Time-averaging over the micromotion, e.g., within a Floquet description, loses this information and makes it challenging to predict the stability of the system.
To determine the resonance condition, we matched the phonon energy with the frequency at which the micromotion crosses critical regions of the Brillouin zone.This cycling frequency allowed us to predict the fundamental and higherorder phonon modes of modulational instabilities in all three scenarios.For off-resonant driving, we replaced the driving frequency with the detuning ν between Bloch frequency and driving frequency.In all cases, a stable, fast-driving limit is reached when the cycling frequency exceeds the energy of the fundamental phonon mode.In addition to phonon resonances, band excitations, and tunneling excitations, a complete stability analysis must include the role of the trapping potential [36].Understanding the joint effects of these excitation mechanisms will be instrumental for quantum simulation experiments that study interacting many-body states with periodic driving over long time scales [21].
The data used in this publication are openly available at the University of Strathclyde Knowledge Base [55].

ACKNOWLEDGMENTS
We acknowledge support by the EPSRC through a New Investigator Grant (No. EP/T027789/1), the Programme Grant DesOEQ (No. EP/P009565/1), and the Quantum Technology Hub in Quantum Computing and Simulation (No. EP/T001062/1).CEC was supported by the Spanish MICINN through Grant No. PID2022-139288NB-I00 and by the Universidad Complutense de Madrid through Grant No. FEI-EU-19-12.

APPENDIX A: DESCRIPTION OF EXCITATION MODES
In this Appendix, we summarize the description of phonon modes in superfluids for lattice potentials with a constant force F B and a driving force F (t ) = F 0 cos(2πν D t ).These forces accelerate the superfluid on a periodic motion within the first Brillouin zone in momentum space.This micromotion k(t ) provides the basis for the subsequent analysis.
Without forces F B = F (t ) = 0, a superfluid with phonon excitations can be described by a carrier wave with quasimomentum hk and a weak perturbation δφ k (z, t ) [50,56]: (A1) where μ k is the chemical potential, and φ k is the solution of the stationary Gross-Pitaevskii equation.The perturbation is expressed as a superposition of Bogoliubov modes, each with quasimomentum hq, amplitudes u kq , v kq , and energy E q = hω q (k).Due to momentum conservation, those excitation modes occur in pairs of opposite quasimomenta (+ hq, − hq).
The energy of a single phonon mode with momentum hq and carrier momentum hk is given by [49,56] E q (k) = 2J sin(kd L ) sin(qd L ) ± 2 4J 2 cos 2 (kd L ) sin 4 where J and U are the tunneling matrix element and the interaction energy.The real part of E q (k) describes the energy and the phase oscillations of the mode in Eq. (A2), while the imaginary part provides its growth rate = Im[E q ]/h.The largest value of Re[E q (k)] for any combination of q and k in the Brillouin zone is W = 2 √ 4J 2 + 2JU .Adding a constant force F B suppresses tunneling when the energy shift F B d L between adjacent lattice sites is much larger than the single-particle band width 4J.However, resonant driving with frequency hν D = F b d L restores the coupling between lattice sites, and the wave packet shows a periodic micromotion through the first Brillouin zone: Here, K = F 0 d L /( hω D ) is the dimensionless driving strength [57].The phase ϕ is set by the switch-on procedure of the forces and their initial directions.We used opposite directions for F B and F (0) in calculations and experiments to minimize the growth of excitations for k 0 = 0.An effective dispersion for the driven lattice system −2J eff (K ) cos(d L k 0 ) can be derived by time-averaging the single-particle energy over one period of the micromotion k(t ) [51].Here, J eff (K ) = JJ 1 (K ) is an effective tunneling element with the first-order Bessel function J 1 (K ).The sign of J eff (K ) is negative in the interval 3.8 < K < 7.0, resulting in an initial momentum k 0 = k L for the ground state.We indicate this change of ground states with separate panels in Figs.3(c), 3(d) and 4(a).
The BdG equations provide the time evolution of a phonon mode [43,48].The components u and v of the mode in Eq. (A2) evolve according to where In the limit of fast-driving frequencies hν D W , the BdG equations can be simplified by time-averaging ± (q, t ) over one driving period and diagonalizing Eq. (A5) [48].The resulting eigenvalues provide the energies of phonon and anti-phonon modes: We included k 0 in the description to extend the time-averaged Bogoliubov equation in Ref. [33].As a result, Eq. (A7) directly matches the Bogoliubov dispersion relation for nondriven systems in Eq. (A3) when replacing J with J eff and k with the initial momentum k 0 .This condition is almost identical for a nontilted system [31], except that the first-order Bessel function in J eff is replaced by the zeroth order.
For slow-driving frequencies hν D W , we time-average E q [k(t )] over one driving period [Eq.( 2)].Here, E s q provides the averaged energy of an existing phonon with momentum q that follows a particular micromotion, while E f q describes the energy of a phonon in a system with time-averaged parameters.Note that E f q does not include the periodic growth of phonon modes during a driving cycle due to modulational instabilities, which is the main subject of this article.The values of both energies Re[E s q ] and E f q diverge close to the zeros of J 1 (K ), but they agree well in between [red and white lines in Fig. 3(d)].

APPENDIX B: CALCULATION OF THE GROWTH RATE OF THE PHONON
We calculated the growth rate q of phonon modes with momentum q using Eq.(A5) and the diagonalization procedure in Refs.[46,48].Figures 7(a strength K = 0, 2, 3.5, with yellow and blue colors indicating stable and unstable parameter regions, respectively.Without driving (K = 0) but with a constant force F B , the micromotion passes only once per cycle through a critical region with nonzero growth [α = 1, Fig. 7(d)].The corresponding time-averaged growth rate shows unstable regions that are boomerang shaped and align with the predicted resonances in Eq. ( 1) [Fig.7(a)].
Adding a driving force K > 0 distorts the micromotion of the superfluid [black line in Fig. 7(e)] and reduces the time it spends in critical regions of the Brillouin zone (red patch).This distortion of the micromotion changes the real and imaginary components of the time-averaged energy E s q .Spending less time in the critical region increases Re[E s q ], as the real component of E q [k(t )] is zero in the critical regions, and it decreases Im[E s q ], as the imaginary component is zero in the noncritical regions of the Brillouin zone.As a result, unstable regions shown in Fig. 7(b) move toward smaller values of T B and shrink in size compared with Fig. 7(a).This is the origin of the stable region (ii) in the (K, T B )-stability diagram in Fig. 3(d).
For larger driving strength, e.g., K = 3.5 in Figs.7(c) and 7(f), the micromotion passes three times per cycle through a critical region.Due to these crossings, the relative time interval between two crossings α decreases to 0.4, and the cycling frequency ν c increases.As a result, the fast-driving limit shifts toward larger values of T B , which is the cause for the stable region (iii) in Fig. 3(d).images either after a rapid switch off or after a linear ramp of the lattice potential over 1.2 ms.For the data analysis in Fig. 2, we determined the atom number in the carrier wave N C by fitting the integrated 1D profiles close to the carrier momentum with a Gaussian fit of the form A exp[(k − B) 2 /C 2 ], where A, B, C are fit parameters with N C = AC √ π .
However, this approach was challenging for measurements with a significant fraction of atoms in excitation modes.Instead, we determined N C for the data in Figs. 3 and 5 directly by counting the number of atoms in a momentum interval that enclosed the expected momentum of the carrier wave [red lines in Fig. 8(a)].The fraction of atoms in excitation modes was calculated as N E /N tot = 1 − N C /N tot .Both methods do not account for phonon modes with small q values [e.g., in Fig. 8(b)] because the phonons and the carrier wave overlap in momentum space.However, excitation modes with large q values dominate for strong interactions U > 4J, as used for the measurements in this article, and we expect that the omission of modes with small q values does not change the overall shape of our stability diagrams.
Changing the lattice depth V during a measurement did sometimes introduce small additional forces that shifted the final momentum of the wave packet [e.g., Fig. 8(c)].Instead of balancing those small forces for every value of V , we included the final momentum shift of the wave packet in our data analysis when setting the momentum interval to determine N C .For example, the scan of V in Fig. 8(c) created a total variation of ∼1% of F B , which we expect to have little influence on the stability measurement.
Figure 9 shows the averaged absorption images for the measurement of the (K, T B )-stability diagram [Fig.3(c)].The system was always prepared in its ground state for the fast-driving limit, i.e., with k 0 = 0 for K = 0-3.8[bottom panel in Fig. 9(a)] and with k 0 = k L for K = 3.9-5.9[top panel in Fig. 9(a)].For an easy comparison of both parameter regimes, we show only the peak of the carrier wave in the bottom panel and the sum of the two wave packets at the edge of the Brillouin zone in the top panel [Fig.9(a)].Figure 9(b) provides averaged absorption images for nonresonant driving, [Fig.5(c)].Again, we show only the momentum interval for the main peak of the carrier wave for reference.

FIG. 1 .
FIG. 1. Experimental setup and characteristic energies.(a) Sketch of a superfluid with density modulations (blue wave packets) in a tilted lattice potential (gray line), with tunneling energy J, interaction energy U , lattice spacing d L , constant force F B , and driving force F (t ).(b) Resonant excitations can occur when the driving energy hν D matches multiples or fractions of the energy of the band gap for lattice depth V , the tilt of lattice potential F B d L , phonon modes E q , or trap frequencies hω z .Dark arrows show resonances of the driving frequencies, and light arrows indicate their fractions for higher-order resonances.(c) Experimental setup with optical lattice, trapping laser beams, magnetic field coils, and levitating force.

FIG. 2 .
FIG. 2. Modulational instabilities without periodic driving.(a) Absorption images of Bloch-oscillating atoms and excitation modes in quasimomentum space for T B = 10 ms, V = 6 E r , U/J = 3.(b) Atoms in excited modes N E /N tot (blue circles) and calculated micromotion k(t ) (black line) for parameters in (a).Red patches indicate time intervals when micromotion crosses into unstable regions of the Brillouin zone.(c) Measured ratio N E /N tot after a hold time of t ≈ 20 ms.Lines indicate Bloch periods that match the singleparticle bandwidth 4J/h (dashed red line), the transition line between unstable and stable parameters for a phonon energy 2E s k L (solid red line), and the transition line predicted in Ref. [43] (black line).Gray patches indicate strong atom loss for small V and large values of F B .

FIG. 3 .
FIG. 3. Stability of a superfluid in a tilted potential with a resonant driving force.(a) Calculated micromotion for K = 0, 1 (black and blue lines) and for K = 1, 3 with opposite initial direction of F (0) and F B (red and orange lines).Thick lines with corresponding colors indicate time intervals in critical regions (gray patches) of the Brillouin zone (|k| > 0.5k L ).(b) Numerical calculation of the growth rate for K = 2, V = 10 E r , U/J = 35 together with predicted resonance (red line) with energy 2Re[E s k L ] in Eq. (1).(c) Experimental measurement of atoms in excitation modes after driving for ∼30 ms with strength K and parameters V = 10 E r , N tot ≈ 60 000, a s = 107 a 0 .Left and right panels use initial momentum k 0 = 0 and k L .Average over typically six repetitions.(d) Numerical calculation of with same parameters as in (c).Red and white lines indicate the predicted position of the fundamental resonance in Eq. (1) when approximating the phonon energy by E s q and E f q .

FIG. 4 .
FIG. 4. Growth rate and driving resonances.(a) Calculated growth rate k L for V = 10 E r and U = 30 J. Dashed black lines indicate local maxima of k L that were determined with a peak-finding algorithm.Solid blue and red lines show the predicted resonance positions in Eq. (1).Panels indicate different signs of J eff (K ).(b)-(d) Examples of the micromotion for the three panels in (a), top and bottom panels refer to regions (L) and (R) with blue and red line colors, respectively.(b) K = 0.5 (top), K = 3.5 (bottom); (c) K = 4.0 (top), 6.5 (bottom); and (d) K = 7.5 (top), 10.0 (bottom).Colored patches indicate time intervals of the micromotion in regions of the Brillouin zone with modulational instabilities.Arrows show α for the longest time intervals between consecutive crossings.

FIG. 5 .
FIG. 5. Stability of the superfluid in a tilted potential with offresonant driving.(a) Micromotion (gray line) and measured number of atoms in phonon modes N E (blue circles) for small detuning ν = 100 Hz with parameters K = 1.5, V = 8 E r , k 0 = 0, T B = 0.998 ms, U/J = 15, and (b) for large detuning ν = 800 Hz.Colored patches indicate time intervals in the critical region of the Brillouin zone.(c) Measured stability diagram showing N E /N tot after ∼30 ms of driving and parameters V = 8 E r , N tot ≈ 60 000, a s = 107 a 0 .

FIG. 6 .
FIG. 6. Phonon growth rate for off-resonant driving.(a) Calculated growth rate k L for V = 10 E r , U = 4J, k 0 = 0, s = 500, ν B = 1 kHz.Dashed red lines show the resonances m p = 1, 2, 3 for K = 1.85 calculated using Eq.(3).Other lines are based on the approximation E tot k L (K ) in Eq. (4) predicting the resonances m p = 1, 2, 3 (gray lines), the stable region between m p = 1 and 2 (white line), and the transition between unstable and stable (black line).(b) Interaction dependence of k L for K = 1.85.Other parameters and line colors match (a).(c) Measurement of N E /N tot for the parameters in Fig. 5(c) and K = 1.90 [dashed gray line in (d)].Data colors indicate the tunneling resonances 1 1 (red color) and 1 2 (blue color).(d) Calculated growth rate for a large range of K and ν D values.Phonon resonances occur at tunneling resonances with frequency ν D = (n t /m t )ν B and parameters s = 42, V = 10 E r , U = 4J, q = k L .
FIG. 7. Calculated growth rate for phonon modes with momentum q. (a)-(c) Growth q for resonant driving hν D = F B d L = h × 1 kHz, V = 10 E r , U = 35J, and driving strength with (a) K = 0, (b) K = 2, and (c) K = 3.5.Red lines indicate Bloch periods for phonon resonances with energy (2E s q )/m p and integer values of m p .(d)-(f) Micromotion k(t ) (black lines), phonon energy Re[E s k L ] (red lines), and growth rate Im[E s k L ]/h (blue lines) for (d) K = 0, (e) K = 2, and (f) K = 3.5.
FIG. 8. Measurement of excitation modes.(a) and (b) Integrated absorption profiles in real momentum space for (a) K = 0.8, T B = 4 ms and (b) K = 2.3, T B = 3 ms.Red lines indicate the regions used to determine the number of atoms in the carrier wave for Fig. 3(c).(c) Absorption images for T B = 2 ms in Fig. 2(c).

FIG. 9 .
FIG. 9. (a) Absorption images used to generate the (K, T B )stability diagram with resonant driving in Fig. 3.Only the area close to the peak of the carrier wave is shown.Each image is the average of typically six absorption images with identical parameters.The top panel shows the sum of the two wave packets at k 0 = k L for an easier comparison with the bottom panel with k 0 = 0. (b) Absorption images for (K, ν)-stability diagram with off-resonant driving for Fig. 5.The images show the area close to the central peak of the carrier wave.The panels are rotated by 90 • compared to Figs. 3 and 5.