Parametric Instability Rates in Periodically-Driven Band Systems

This work analyses the dynamical properties of periodically-driven band models. Focusing on the case of Bose-Einstein condensates, and using a meanfield approach to treat inter-particle collisions, we identify the origin of dynamical instabilities arising due to the interplay between the external drive and interactions. We present a widely-applicable generic numerical method to extract instability rates, and link parametric instabilities to uncontrolled energy absorption at short times. Based on the existence of parametric resonances, we then develop an analytical approach within Bogoliubov theory, which quantitatively captures the instability rates of the system, and provides an intuitive picture of the relevant physical processes, including an understanding of how transverse modes affect the formation of parametric instabilities. Importantly, our calculations demonstrate an agreement between the instability rates determined from numerical simulations, and those predicted by theory. To determine the validity regime of the meanfield analysis, we compare the latter to the weakly-coupled conserving approximation. The tools developed and the results obtained in this work are directly relevant to present-day ultracold-atom experiments based on shaken optical lattices, and are expected to provide an insightful guidance in the quest for Floquet engineering.


I. INTRODUCTION
Periodically-driven systems have been widely studied over the past decades, revealing rich and exotic phenomena in a large class of physical systems. Some of the original applications include examples from classical mechanics, e.g. the periodically-kicked rotor [1] and the Kapitza pendulum [2], which display fascinating effects, such as dynamical stabilization/destabilization [3] or integrability to chaos transitions. More recently, studies of their quantum counterparts have led to remarkable applications in a wide range of physical platforms, such as ion traps [4], photonic crystals [5], irradiated graphene [6][7][8] and ultracold gases in optical lattices [9,10].
A major reason for the renewed interest in periodicallydriven quantum systems is the fact that they constitute the essential ingredient for Floquet engineering [6,[9][10][11][12]. To see this, recall that periodically-driven quantum systems obey Floquet's theorem [13], which implies that the time-evolution operatorÛ (t, 0) of a system, whose dynamics is generated by the HamiltonianĤ(t + T ) =Ĥ(t), can be written as [11,14] U (t, 0) = T exp −i = e −iK kick (t) e −itĤ eff e iK kick (0) , where T denotes time-ordering. Here we introduced the time-independent effective HamiltonianĤ eff , as well as the time-periodic "kick" operatorK kick (t + T ) =K kick (t), which has zero average over one driving period [11]. It then follows that the stroboscopic time-evolution (t = N T , N ∈ Z) of such systems is governed bŷ U (N T, 0) = e −iK kick (0) e −iN TĤ eff e iK kick (0) = e −iN TĤF , H F = e −iK kick (0)Ĥ eff e iK kick (0) .
Hence, up to a change of basis [15], the stroboscopic time-evolution is completely described by the effective HamiltonianĤ eff , which can be designed by suitably tailoring the driving protocol [11,12,[16][17][18]. Moreover, we note that the dynamics taking place within each period of the drive, the so-called micro-motion, is entirely captured by the kick operatorK kick (t). On the theoretical side, perturbative methods to compute both the effective Hamiltonian and the micromotion operators have been developed in the form of an inverse-frequency expansion [11,14,16,19], variants of which include the Floquet-Magnus [12,20], the van Vleck [17], and the Brillouin-Wigner [18] expansions. On the experimental side, Floquet engineering has been used to explore a wide variety of physical phenomena, ranging from dynamical trapping of ions in Paul traps [4] to one-way propagating states in photonic crystals [5]. In the context of ultracold quantum gases, this concept soon resonated with the idea of quantum simulation [21], as it became clear that Floquet engineering could offer powerful schemes to reach novel models and properties, typically inaccessible in conventional condensed matter materials [9,10]. For instance, shaken optical lattices have been used to explore dynamical localization [22][23][24][25][26], photon-assisted tunneling [27,28] and driven-induced superfluid-to-Mott-insulator transi- shaken BEC trapped in a one-dimensional lattice (along x), and allowed to move along a continuous (tube) transverse direction (y), develops unique resonance structures in the presence of a periodic drive. Dominating the dynamics in the early stage of evolution, they signal the onset of instabilities and provide a clear experimental signature of the parametric resonance phenomenon explained in this work. The momentum distribution is shown after 60 driving periods (t = 60T ). Here, the distribution nq is divided by the system's volume Vol; ξ denotes the BEC healing length and ax is the lattice spacing constant. The precise model and the corresponding parameters are the same as in Fig. 17.
Today, the theory of periodically-driven quantum systems (including the resulting artificial gauge fields and topological band structures [9,10]) is well-established for single particles, and there is a strong need for incorporating inter-particle interactions in the description of such systems. Indeed, inter-particle interactions constitute an essential ingredient in simulating condensed-matter systems [21], and the relevance of many theoretical propos-als strongly relies on the extent to which the systems of interest remain stable away from the non-interacting limit. Moreover, in their quest for realizing new topological states of matter, experiments based on driven quantum systems are more than ever eager to consider interacting systems [59]; indeed, building on the successful optical-lattice implementation of flat energy bands with non-trivial Chern numbers [39], mastering interactions in this context would allow for the engineering of intriguing strongly-correlated states, such as fractional Chern insulators [60]. In this framework, we point out that a useful tool to derive the effective low-energy physics of strongly-correlated systems, namely the Schrieffer-Wolff transformation, was generalized to periodically-driven systems [58].
However, until now, the presence of inter-particle interactions has led to complications. Recent experiments on time-modulated ("shaken") optical lattices [36][37][38][39]61] have reported severe heating, particle loss and dissipative processes, whose origin presumably stems from a rich interplay between the external drive and inter-particle collisions. In this sense, understanding the role of interactions in driven systems, and its relation to heating processes and instabilities, is both of fundamental and practical importance. On the theory side, several complementary approaches have been considered to tackle the problem. On the one hand, a perturbative scattering theory, leading to a so-called "Floquet Fermi Golden Rule" [62,63], has been developed to estimate heating and loss rates, and showed good agreement with the band-population dynamics reported in Ref. [39]. Extensions to Bose-Einstein condensates (BEC) have been proposed [64,65], however, to the best of our knowledge, the resulting loss rates have not yet been confirmed by experiments nor through numerical simulations. On the other hand, various studies analyzed the dynamical instabilities that occur in BEC trapped in moving, shaken or timemodulated optical lattices [66][67][68][69][70][71][72]. For instance, Ref. [72] analyzed such dynamical instabilities through a numerical analysis of the Bogoliubov-de Gennes equations; however, no link was made with heating processes and the approach lacked a conclusive comparison with analytics. Importantly, several studies pointed out the parametric nature of those instabilities [69,70,73,74]. Interestingly, we note that parametric instabilities, and their impact on energy absorption and thermalization processes, have also been studied theoretically in a wider context, ranging from Luttinger liquids [75] to cosmology [76][77][78]. In the context of driven optical lattices, the parametric amplification of scattered atom pairs, through phasematching conditions involving the band structure and the drive, has also been identified [79][80][81]. Altogether, there is a strong need for a combined analytical and numerical study of driven optical-lattice models, which would provide analytical estimates of instability and heating rates supported by numerical simulations.

Scope of the paper
The object of the present paper is to explore the physics of periodically-driven bosonic band models, where inter-particle interactions are treated within meanfield (Bogoliubov) theory. Indeed, as in Refs. [64,65,72,74], we will assume that the driven atomic gas forms a Bose-Einstein condensate (BEC), and we will analyze the properties of the corresponding Bogoliubov excitations in response to the external drive. To do so, we develop a generic and widely-applicable method to investigate the short-time dynamics of periodically-driven bosonic systems and to extract the corresponding instability rates, both numerically and analytically.
We will first focus on a shaken one-dimensional (1D) lattice model, and extensively explore the occurrence of instabilities and heating in this system, before considering more elaborate models. In a previous study on a similar model [72], a dynamical instability of the condensate was found to occur above a critical interaction strength, which was numerically calculated as a function of the drive amplitude K and frequency ω. This allowed the author to numerically map out the stability boundaries between the stable and unstable phases. More recently, Ref. [73,74] suggested that such boundaries could be understood in terms of a simple energy-conservation criterion, which could be formulated within a parametricinstability analysis. In this work, we offer a significant advance in the description and understanding of the instabilities that affect bosonic periodically-driven systems, both qualitatively and quantitatively, as we summarize below: 1. We show how to obtain the instability rates as a function of the model parameters, both numerically and analytically, and we establish their relation to the growth of physical quantities in the system, such as the energy, thus providing a link between dynamical instability and heating. Moreover, we identify different instability regimes in the system, which are associated with different timescales and characterized by a different behaviour in the growth of physical quantities. We introduce several observables that reveal clear signatures of these instabilities, among which the non-condensed (depleted) fraction and the momentum distribution of quasiparticle modes [see Fig. 1]; we discuss how they could be observed in current ultracold-atom experiments.
2. We present a full numerical solution of the meanfield problem, determining the stability diagram of the model from first principles [see Fig. 2]. We stress that the quantitative character and versatility of our procedure make it general enough to be applicable to a wide class of periodically-driven systems, including those involving resonant modulations, higher dimensional lattices, various geometries (e.g. square or honeycomb lattices, continu- (2)] and interaction strength g, for a driving frequency ω = 5J. The vertical lines and letters summarize the agreement and disagreement regions, as discussed in Sec. IV C: A: Away from the zeros of J0(K/ω) and J2(K/ω), the analytics successfully capture the instability rates extracted from numerics. B: Agreement zone near the (close) zeros of J0(K/ω) and J2(K/ω), where both analytics and numerics predict a stable behaviour. C/D: Close to a zero of J0(K/ω), the analytical perturbative approach breaks down. E: Potential quantitative disagreement when the contribution of the second harmonic has a weak amplitude: the instability is then partly due to higher harmonics (see Sec. IV B).
ous space), and spin-dependent lattices. Moreover, it can also be enriched to incorporate a full experimental sequence, e.g. including adiabatic state preparation [61]. Our results are expected to determine experimentally-favorable regimes by providing ab initio numerical (and analytical) estimates for the instability rates in a variety of experimental configurations.
3. We perform an analytical treatment of the problem based on the existence of parametric resonances in the Bogoliubov-De Gennes equations. In Ref. [73], following a weak-coupling conserving approximation, it was argued that a driven-lattice model was stable against parametric resonance provided that the drive frequency satisfies ω > 2W eff , with W eff the bandwidth associated withĤ eff in the Bogoliubov approximation. Intuitively, this can be understood by recalling that the elementary excitations of the system (i.e. the Bogoliubov phonons) are always created in pairs, and thus, whenever the drive frequency exceeds twice the maximum possible excitation energy, stability is ensured by energy conservation. However, our rigorous analytical derivation reveals that this simple criterion needs to be revised: As we demonstrate, an accurate qualitative explanation requires taking into account both higher-order photon-absorption processes, as well as the detuning from resonance within the parametric-resonance treatment [82]. This analysis allows us, for the first time, to derive the functional dependence of the instability growth rate on the model parameters, and ultimately, to understand all features of the stability diagram, see Fig. 2.
4. We extend the stability analysis to two-dimensional (2D) models and study the effects due to transverse directions, both considering the case of a transverse lattice (i.e. a discrete transverse degree of freedom), and that of transverse tubes (i.e. a continuous transverse degree of freedom). In the case of continuous degrees of freedom, a theoretical simplification in the equations allows one to obtain very simple analytical formulas for instability rates, which are in perfect agreement with our full numerical simulations. Moreover, this study provides a clear physical picture of how instabilities are enhanced by the presence of transverse modes in the system, as already anticipated in Refs. [63,65]. More generally, we argue that our theory should provide guidance for experiments; it illustrates, for instance, the advantage to work at high frequency and the necessity of using a strong transverse confinement to reduce parametric instabilities. We also include a discussion of finite-size effects, making a link between the physics of double-wells and that of optical lattices.

Outline
This paper is organized as follows. In Sec. II, we derive the general mean-field equations that will constitute the basis of our analysis. Section III is devoted to the numerical solution of those equations, including details on the general procedure used and a presentation of the results, among which the stability diagram of the model under consideration, the identification of several timescales and instability regimes in the problem and the dynamics of various physical observables. In Sec. IV, we perform an analytical treatment of the problem: mapping the Bogoliubov equations on a parametric oscillator (Sec. IV A), we build an effective model from which analytical instability rates can be inferred (Sec. IV B). The analytical results are presented in Sec. IV C, including a discussion of the validity regimes of the approach. The case of finite-size systems is presented in Sec. V. We discuss in Sec. VI the case where a transverse direction is present (a lattice or a continuous one); this includes simple analytical formulas for instability rates as well as a physical understanding of the enhancement of instabilities by transverse modes. Finally, we discuss in Sec. VII the application of our results to the weakly-interacting Bose-Hubbard model in the meanfield regime; in order to understand the role of non-linear processes and study the regime of validity of our linearized analysis, we employ a weak-coupling conserving approximation to study the leading-order (in the interaction strength) features of particle-conserving dynamics; we identify clear signatures of the instabilities (among which the non-condensed fraction and the momentum distribution of quasiparticles) that could be directly probed by current ultracold-atom experiments in modulated optical lattices.

II. MEAN-FIELD EQUATIONS AND STABILITY ANALYSIS
Consider a Bose gas in a shaken 1D lattice, with meanfield interactions, governed by the time-dependent Gross-Pitaevskii equation (GPE) [10,72] i∂ t a n = −J (a n+1 + a n−1 ) + K cos(ωt)na n + U |a n | 2 a n , (2) where n labels the lattice sites, J > 0 denotes the tunneling amplitude for nearest-neighbour processes, U > 0 is the on-site interactions strength, and where the timeperiodic modulation has an amplitude K and a frequency ω = 2π/T . Note that we set = 1 throughout this work.
Equation (2) describes a wide variety of physical systems. On the one hand, it is expected to provide a good description for the shaken 1D Bose-Hubbard model [10], when treated in the weakly-interacting regime: in this case, the time-modulated system in Eq. (2) can be realized by mechanically modulating an optical lattice filled with weakly-interacting bosonic atoms [10]; see also Sec. VII for further discussion. On the other hand, some physical systems are "exactly" described by Eq. (2): for instance, non-linear optical systems [83], including helical photonic crystal [5]. In all cases, Eq. (2) defines a close self-consistent problem, which constitutes the core of the present study.
Similar to the analysis of Ref. [72], we study the dynamical instabilities of a BEC described by Eq. (2). To do so, we proceed along the following steps: (i) First, we determine the time-evolution of the condensate wavefunction a (0) n (t), by solving the full time-dependent GPE and assuming that the initial state a (0) n (t = 0) forms a BEC. More precisely, our choice for the initial state corresponds to the solution of the static (effective) GPE: where we introduced the effective tunneling amplitude J eff renormalized by the drive [10], and where J 0 denotes the zeroth-order Bessel function. For the model under consideration, we find that a (0) n (t = 0) is the Bloch state e ip0n of momentum p 0 = 0 if J 0 (K/ω) > 0 (homogeneous condensate), and p 0 = π for J 0 (K/ω) < 0. Note that such a choice takes into account the initial kick due to the launching of the modulation [11] [84]. Altogether, this choice for the initial state features both the effects of tunneling renormalization and initial launching of the drive, which are both present in our model (2). We emphasize that this prescription for the initial state is the only step of our calculations that relies on the existence of a well-defined high-frequency limit, as provided by the inversefrequency expansion [11,14,16]; indeed, all subsequent results are based on the full time-dependent equations. Note that it is not the purpose of this paper to discuss how to prepare the system in this initial state; one possibility is to use Floquet adiabatic perturbation theory, see Ref. [61,85], which is also the experimentally-preferred strategy.
(ii) Given the time-dependent solution for the condensate wave function a (0) n (t), we analyse its stability by considering a small perturbation a n (t) = a (0) and linearizing the Gross-Piteavskii equation (2) in δa n ; we refer to Appendix A for a discussion on the specific parametrization chosen in Eq. (4). This yields the time-dependent Bogoliubov-de Gennes equations, which take the general form where we introduced the operator F n (t), whose action on δa n is defined by and where the discrete operatorL is defined bŷ L(·) n ≡ (·) n+1 + (·) n−1 .
At this stage, let us emphasize that the Bogoliubov equations (5) contain the complete time dependence of the problem, which includes the effects related to the micro-motion (note that the BEC wavefunction a (0) n (t) is computed exactly, not stroboscopically). Importantly, and as will become apparent below in Section IV, it is the micro-motion (and not the time-averaged dynamics) that determines the stability of the system.
The Bogoliubov equations of motion (5) can be timeevolved over one driving period T , which allows one to determine the associated "time-evolution" (propagator) matrix Φ(T ). From this, we extract the "Lyapunov" exponents q , which are related to the eigenvalues λ q of Φ(T ) by λ q = e −i q T . The appearance of Lyapunov exponents with positive imaginary parts thus indicates a dynamical instability [70,72], i.e. an exponential growth of the associated modes at the rate s q = Im q , where s q denotes the growth rate of the momentum mode q. As we shall see later, a quantitative indicator of the instability is the maximum growth rate of the spectrum, which, in the following, will be referred to as the instability rate Γ. This choice will be justified in Sec. III B.
For the model under investigation, one can considerably simplify the problem by working in a rotating frame, in which translational invariance is manifest. This is achieved by applying the following gauge transformation [86] a n (t) −→ e −in(K/ω)sin(ωt) a n (t). (8) In this frame, and going to momentum space, the system of coupled Bogoliubov-de Gennes equations reduces to a 2×2 matrix [72], and the dynamics of the mode with momentum q ∈ BZ (first Brillouin zone) is described by [87] i∂ t where ε(q, t) = 4J sin with p 0 denoting the momentum of the initial BEC wavefunction. Throughout this paper, we set the lattice constant to unity (a x = 1). Here, we introduced the parameter g ≡ ρU , where ρ = |a (0) (0)| 2 denotes the condensate density; the latter enters the normalization of the solution a Equation (9) is numerically easier to integrate, and will also constitute the starting point of the analytical analysis (see Sec. IV).

III. NUMERICAL SIMULATION
The general procedure outlined above can be implemented and solved numerically, regardless of the precise form of the physical band model. In more involved models, the initial state a  n (t) is then determined by solving the time-dependent GPE with this initial condition, through real-time propagation (e.g. using a Crank-Nicolson integration scheme). Finally, the Bogoliubov equations are also solved over one driving period by real-time propagation, yielding the operator Φ(T ), which is then exactly diagonalized (e.g. using a Lanczos algorithm). In the present case, this procedure can be shortcut by directly solving Eqs. (9), and we have checked that this produces the same results for all physical quantities. We now present the numerical results obtained for our model.

A. Stability diagrams
The stability diagram of the model in Eq. (2), which displays the behavior of the instability rate Γ as a function of the interaction strength g = U ρ and modulation amplitude K/ω, is shown in Fig. 2, for a reasonably large driving frequency ω = 5J. The stability boundary is similar to the one previously reported in [72] [88]; in particular, the system is found to be stable in regions where J 0 (K/ω) = 0, which can be attributed to the fact that the dynamics is frozen by the cancellation of the effective tunneling [72]. At the transition to instability, the instability rate builds up continuously from zero, and then increases when going further in the unstable regime. Figure 3 shows similar stability diagrams for two other values of the driving frequency, ω = 10J and ω = J. In the first case, we observe that the diagram is mostly unaffected by a change of ω provided g is rescaled as g ∝ ω 2 . More generally, except in what we shall refer to as the "low-frequency regime" (defined by the condition ω < 4|J eff | for the present model; see the next paragraph), transitions to instability always occur at some finite g, and the corresponding rate is found to mostly depend on the quantities K/ω and g/ω 2 [89]. As we shall see in Sec. IV, this can be understood from the fact that the instability rate depends on a competition between ω and the Bogoliubov dispersion associated with the linearised effective GPE [i.e. the Bogoliubov dispersion stemming from the linear analysis of Eq. (3)], E av (q) = 4|J eff | sin 2 (q/2)(4|J eff | sin 2 (q/2) + 2g), (12) where the two cases p 0 = 0, π are taken into account through the absolute value |J eff |. As soon as the transition occurs at sufficiently large g, the term 4|J eff | sin 2 q/2 in this dispersion becomes negligible compared to g, At large ω (i.e. outside the "low-frequency" regime, see text), the instability rates depend mostly only on K/ω and g/ω 2 , while in the low frequency regime, instabilities can occur at infinitesimal interaction strength.
which results in E av ∝ √ g, explaining the scaling indicated above. This is no longer true in the "low-frequency regime" (see Fig. 3 for ω = J), which we generically define through the criterion according to which ω is smaller than the effective free-particle bandwidth (i.e. ω < 4|J eff | for the model under consideration). In that case, we observe that instabilities may occur at any finite interaction strength g, which is related to the fact that ω is smaller than the bandwidth of the effective Bogoliubov dispersion (12) at g = 0; see also the analytical analysis in Sec. IV for more details.
Finally, in the singular case where J eff = 0, the effective Bogoliubov dispersion E av (q) in Eq. (12) vanishes, and thus, the system is necessarily stable. In the stable regime (top, g = 4), the energy displays variations due to micromotion, but no long-term growth, while in the unstable regime (bottom, g = 7), the energy curve features exponential growth, defining a heating rate compatible with the maximal Lyapunov exponent Γ [see Fig. 5]. The parameters are Aq = 1 and Bq = 0 (for the initial condition, see text), ω = 5J, and K/ω = 7.5.

B. Dynamics of physical observables
Instead of solving the Bogoliubov equations [Eq. (5)] over one driving period only, we can also use these to compute the full time-evolution of physical quantities, hence revealing both their long-time and micromotion dynamics. To do so, one has to choose an initial condition for the fluctuation term, δa n (t = 0). In this section, we will consider a generic (small) perturbation δa n (t = 0), which has Fourier components over all Bogoliubov modes, i.e. of the generic form δa n (0) = α q A q e iqn + B q e −iqn with A q , B q being complex numbers and α a small amplitude. In Sec. VII, we discuss a more physical initial condition, in the framework of the periodically-driven Bose-Hubbard model. Altogether, given such an initial condition, one can numerically evolve Eqs. (5) using real-time propagation, which yields δa n (t) (and thus also a n (t) = a in the linearized approximation), hence revealing the behavior of physical quantities.

Energy growth and heating
As an illustration, we first evaluate the energy in the rotating frame, which is given at time t by In the linear approximation, one can explicitly recast this expression as a function of δa n and only keep terms of lowest order in δa n . This yields where the energy of the fluctuations E fluct (t) is given to lowest order by The behavior of the energy of the fluctuations in Eq. (14) over several driving periods is shown in Fig. 4, in the stable and the unstable regimes [red curves], for an initial condition corresponding to A q = 1 and B q = 0. In the stable regime, it displays modulations due to micromotion, but no long-term growth. Conversely, as anticipated for a parametric instability, the energy displays an exponential growth in the unstable regime. By fitting this growth, we find that its rate is given by 2Γ [the factor 2 stems from the square moduli (|δa n | 2 , ...) in Eq. (14); see also below, Eq. 17], which expresses the fact that the maximally unstable mode dominates the growth of the energy in the system. This unstable behaviour, and the validity range of the related regime, will be further addressed in the next paragraphs. Figure 5 shows the diagram obtained by extracting the heating rates from energy curves, and is found to be in excellent agreement with the stability diagram of Fig. 2. Therefore, we conclude that the instability rate Γ can be used as a quantitative estimator of heating in such systems.

Validity regimes of the calculation
The computation of the full dynamics also highlights in which regimes our previous calculation of Γ is expected to hold. On the one hand, the calculation of the rate Γ is based on a Floquet treatment of the Bogoliubov equations Eq. (5), and thus describes the dynamics at times t T.
In turn, the dynamics within one driving period, appearent in Fig. 4, cannot be captured by our analysis.
On the other hand, the analysis presented so far is based on the linear approximation of the GPE. To investigate non-linear effects, we compared the time-evolution obtained from the linearized Eq. (5) (red curves in Fig. 4, as previously discussed) with the time-evolution of the original non-linear GPE in Eq. (2) (grey dashed lines in Fig. 4), for the same initial condition. While the two graphs agree reasonably well in the stable regime, we find a significant deviation at longer times, within the unstable regime, due to the growth of the fluctuation term. This illustrates the intuitive fact that our linear-approximation-based treatment only holds at short times, imposing a second validity condition on our theory: At longer times, non-linear corrections to the evolution damp the exponential growth of the energy, which leads to a slowdown (saturation) of the heating dynamics. Altogether, we find two natural time scales in the problem, and our approach thus requires the condition T t lin to be relevant. Such a condition is fulfilled in a wide window of realistic system parameters.

Instability regimes
A third natural time scale arises when studying the growth of physical observables. In the linear treatment of , for ω = 5J, K/ω = 7.5 and g = 5.7. In this regime where Γ is small, three regimes clearly show up: at very short times t T , the initial growth is not described by our Floquet approach, as previously explained; at intermediate times T t Γ −1 , the growth is linear and described by Eq. 18 (red dotted line); at large times t Γ −1 , the growth is exponential and described by Eq. 17 (red short-dotted line).

the GPE, a generic physical observable can be expressed as
where O[u q (t), v q (t)] is a functional, which is quadratic in the Bogoliubov modes u q (t) and v q (t); for instance, this is straightforward for the energy when substituting Eq. (11) into Eq. (14). This also applies to the non-condensed fraction when considering a weakly-interacting Bose-Hubbard model in the meanfield interacting regime [see Eq. (45) in Sec. VII]. Using the fact that in our Floquet treatment of the Bogoliubov equation, each mode q stroboscopically evolves according to the rate s q , one can rewrite Eq. (15) as where the time t N ≡ N T denotes an integer multiple of the driving period T . By introducing Γ = max q s q , as considered above, two cases may arise: a) if Γt N 1, the growth of the maximally unstable mode dominates in Eq. (16) and we find that the observables stroboscopically grow up exponentially with the rate 2Γ, as already observed for the energy in Sec. III B 1.
Remarkably, this rate only involves the most unstable mode, and it is the same for all physical quantities. As announced above, this justifies the use of Γ as a global instability rate in our theoretical analysis.
. (18) In this regime, the energy increases linearly in time, with a rate given by the slope The latter now involves a summation over all modes, since no single mode contribution can be singled out at these early times (t N 1/Γ). Note also that, in this regime, the rate depends on the physical observable considered (through the functional O).
For most values of the system parameters [such as those used in Fig. 4], ΓT 1, so that the linear regime [Eq. (18)] is hidden in the first oscillation (which, as stated above, is not accurately captured by our Floquet approach). Yet, very close to the stability boundary, where Γ is small, this marginal regime can become apparent in a very narrow range of parameters, as we illustrate in Fig. 6.
Here, three regimes clearly show up: at very short times t T , the initial growth is not described by our Floquet approach, as previously explained; at intermediate times T t Γ −1 , the growth is linear and accurately described by Eq. (18); at longer times t Γ −1 , the growth is exponential and well described by Eq. (17). Note that the linear regime is expected to be hard to probe in experiments, since it appears at very short times, typically a few milliseconds (see also discussion in Sec. VIII).

A. The linearised Gross-Pitaevskii equation as a parametric oscillator
In momentum space, the Bogoliubov-de Gennes equations (9) constitute a set of uncoupled equations, which independently govern the time evolution of excitations with a given momentum q: the mean-field analysis effectively reduces to a single-particle problem, and one can study how instabilities occur independently in each of those modes. The general mechanism underlying the appearance of parametric instabilities was suggested in Ref. [73], where a weakly-interacting Bose-Hubbard model was investigated through a weak-coupling particlenumber-conserving approximation; in that study, parametric resonances were shown to appear in the Bogoliubov equations, whenever the drive frequency was reduced below twice the single-particle bandwidth. A similar approach was considered in Ref. [90], which also provides an intuitive picture of the underlying mechanisms, in terms of the so-called Krein signature associated with the Bogoliubov modes.
The same phenomenon occurs in the present framework, as can be seen by performing a change of basis that recasts Eq. (9) into the following form [see Appendix C for details] Here E av (q) is the Bogoliubov dispersion associated with the effective (time-averaged) GPE, within the Bogoliubov approximation [see Eq. (12)]; W q (t) is a diagonal matrix of zero average over one driving period, which will play no role in the following [see Appendix C for its exact expression]; and h q (t) is a (real-valued) function which can be Fourier expanded as with J l (z) the l-th Bessel function of the first kind. Importantly, the specific form of Eq. (19) allows one to clearly distinguish between the contribution due to the time-averaged dynamics, as captured by the effective dispersion E av (q), and the contribution due to the micromotion [W q (t), h q (t)]. As will be shown below, it is the micro-motion contribution that governs the existence of instabilities in the system, through the properties of the real-valued function h q (t).
Casting the equations of motion in the form (19) makes it possible to directly identify any parametric resonance effects. To see this, we recall the reasoning of Ref. [73], which was based on applying a rotating-wave approximation (RWA) to Eq. (19). If one disregards for now the expression in Eq. (20), and if one naively assumes that the lowest frequency appearing in h q (t) is ω, then a dominant contribution to the dynamics is expected when the resonance condition ω = 2E av (q) is fulfilled, resulting in the time-independent non-diagonal term in the matrix displayed in Eq. (19). Keeping this resonant term only yields a time-independent 2 × 2 matrix which can be diagonalized, and whose eigenvalues (i.e. the Lyapunov exponents) are found to exhibit an imaginary part, yielding a non-zero instability rate. This argument was invoked in Ref. [73] to justify the intuitive stability criterion ω > 2W eff (with W eff = 4|J eff |(4|J eff | + 2g) the bandwidth of the effective Bogoliubov dispersion), which states that instability arises from the absorption of the energy ω to create a pair of Bogoliubov excitations on top of the condensate. However, a more careful examination shows that the instability rates inferred from this idea are not consistent with our numerical simulations.
Moreover, to get a flavour of the additional dilemma one is faced with, we note that the expression for h q (t) in Eq. (20) actually only contains even harmonics of the modulation, and therefore, following the reasoning above, the resonance condition should then read 2ω = 2E av (q). Yet, such a criterion provides a wrong estimate of the stability-instability boundaries. Hence, it appears that this simple explanation needs to be thoughtfully revised.
In the following, we present a more rigorous analytical evaluation of the instability rates for our system. Similar to Ref. [73], our derivation relies on that the Bogoliubov equations (19) precisely take the form of a so-called parametric oscillator. In Appendix B we briefly recall some useful results about this paradigmatic model [3], which describes a harmonic oscillator of eigenfrequency ω 0 driven by a weak sinusoidal perturbation of frequency ω and amplitude α. In brief, such a model displays a parametric instability for ω ≈ 2ω 0 . Importantly, the instability is maximal when this resonance condition is fulfilled, but occurs in a whole range of parameters around this point. As detailed in Landau and Lifshits [3], the width of the resonance domain and the instability rates in the vicinity of the resonance can be calculated perturbatively by introducing the detuning δ = ω − 2ω 0 and solving the equations perturbatively in δ and amplitude α (see Appendix B).
More specifically, the Bogoliubov equations (19) exactly take the form of the parametric oscillator [Eq. (B2)], with the frequency ω 0 of the unperturbed oscillator being identified with the dispersion E av (q). Thus, it immediately follows that the model is equivalent to a set of independent parametric oscillators, one for each mode q. Two points should be emphasized here: first, all those parametric oscillators depend on q in a different way, and will thus exhibit different resonance conditions; second, the function h q (t) in Eq. (19) is not a pure sinusoid as in Eq. (B2), but rather contains all (even) harmonics of the driving frequency ω. Altogether, given those two remarks, resonances are expected as soon as one of the harmonics of the modulation, of energy lω, is close to twice the energy of any of the (effective, time-averaged) Bogoliubov modes, 2E av (q).
Before digging into more technical details, let us acquire some intuition about the general mechanisms behind this phenomenon. To do so, consider the stability diagram on Fig. 2 (i.e. outside of the "low-frequency regime" [91]). In this case, in the lowest part of the stability diagram (small g), ω is generically large compared to the bandwidth of the effective dispersion E av (q), so that no resonance can occur. When increasing g at fixed K/ω, the first mode to exhibit a resonance will be that of maximal E av (q), i.e. q = π. Naively, the first resonance would be due to the first harmonic l = 1 and would occur for 2E av (q = π) = ω, which was also the argument in [73]. However, as already discussed, the function h q (t) only contains even harmonics, and the first resonance to occur is, therefore, the one corresponding to l = 2. At first sight, this might seem to be in contradiction with the intuitive stability criterion 2E av (π) = ω, since the l = 2 resonance is centered around E av (π) = ω; however, as we shall discuss in more detail below, the key point is that the resonance domain has a finite width: although centered around E av (π) = ω, its boundaries are in fact close to the point where 2E av (π) = ω. Resonances due to higher harmonics (fourth, sixth, etc.) become important only for higher values of g, and we can, to first approximation, restrict the analysis to the second harmonic only.
To confirm this intuition, we have verified that restricting our analysis to the second harmonic only is generally sufficient to estimate the stability boundary, and to recover the instability rate in its vicinity [92]. Deeper in the unstable region, the fourth, and eventually the sixth harmonics are progressively required to recover the agreement with the numerical results [see Fig. 7]

B. Effective model
To simplify the calculations, we restrict our analysis to the second harmonic (see discussion above); we stress that higher-order resonances can be treated along the same lines. The equations of motion [Eq. (19)] then read with Therefore, for a fixed q, these equations are now strictly equivalent to the equations of motion for the parametric oscillator of Eq. (B2), with the following iden-tifications:  (20)], the odd harmonics do not contribute, while the first even ones successfully describe the full numerical calculation. Bottom : Comparing the analytical and numerical results for the instability rate Γ(g), using the same parameters as above. The red crosses correspond to the full numerical solution, while the blue dashed lines are numerically obtained by only taking the second harmonic of the drive into account. As generically observed in the A zones of the stability diagram [ Fig. 2], the numerical instability rates are well captured by the perturbative analytical approach, implemented here up to order 1 (red dotted-dashed line) and 2 (blue dotted line).
As a result, each mode q exhibits a resonance domain centred around E av (q) = ω. The instability rate is generically maximal around the resonance point (defining the maximal rate s max q ) and it decreases until it cancels at the edges of the resonance domain. In order to obtain analytical estimates of the instability rate and the width of the resonance domain, one can apply the procedure detailed in Appendix B, which amounts to solving the problem perturbatively in the detuning δ = 2ω − 2E av (q) and amplitude α q . To zeroth order, the instability only arises when the resonance condition is precisely fulfilled, which defines the rate "on resonance" s * q , which is given by To first order, the instability rate reads [see Eq. (B3) with the substitutions (22)] whenever the argument in the square root is positive, and s q = 0 otherwise. The associated instability rate s q is thus maximal on resonance (hence, s max q = s * q ) and it decreases when going towards the boundaries of the resonance domain. At second order, s q is given by the solution of the implicit equation (B4) with the substitutions (22), which can be found numerically; this can be performed to improve the accuracy of the analytical rate's value. In particular, at this order, we note that the actual maximal instability point is then slightly shifted from the resonance point (i.e. s max q = s * q ). Altogether, the total instability rate is given by [Eq. (7)] As a function of g [which enters E av (q), see Eq. (12)], the instability rate of a given mode q forms a "bumped curve" [see Fig. 8], which to first approximation simply follows from Eq. (24). Since the resonance domain for each mode is centered around a different energy [and therefore a different g], the curves corresponding to various modes are slightly shifted, and the total instability rate is given by the envelope of all those curves [see Fig. 8].
In particular, the curve centered around the smallest values of g (i.e. the curve most located to the left in Fig. 8, near the transition point) is the one associated with the mode of highest E av , namely q = π. Let us denote by g max the value of g where this curve is maximal: at first order, this corresponds to the solution of the resonance condition ω = E av (π) = 4|J eff |(4|J eff | + 2g). Then, we identify two cases: (a) For g > g max , there always exists, in the thermodynamic limit, a single mode q which is maximally unstable, and the total instability rate, Eq. (25) is thus given by the rate of this particular mode. At first order [see discussion above], this mode is the one fulfilling the resonance condition ω = E av (q), and the instability rate is given by the rate on resonance of this particular mode q res , so that : where we have explicitly re-expressed q res using the expression ω = E av (q res ). (b) Conversely, for g < g max , the instability rate is only governed by the mode q = π [see Fig. 8], and it cannot be estimated from the knowledge of the associated rate on resonance. Therefore, to capture the behaviour of the instability rate near the transition point, as well as the stability boundary itself, one can indeed restrict our analysis to the mode q = π, but one has in turn to resort to the perturbative expansion around the resonance point discussed in Appendix B. In other words, in that case, one has where s π is given at lowest order by Eq. (24) with q = π. The stability boundary can also be computed at any required order following the procedure indicated in Appendix B, which yields, up to second order, (28) Interestingly, in the high-frequency regime, since the transition occurs at sufficiently large g, so that E av ∝ √ g, one recovers that the boundary is a function of the combination g/ω 2 , as previously observed numerically in Sec. III.
When decreasing the frequency, the transition occurs at lower and lower values of g: since all corrections of any order tend to vanish at small g in Eq. (28), the transition point at vanishing g is simply obtained for ω = 4J eff , recovering the criterion for entering the "low-frequency" regime, previously discussed. In this regime, which in some sense corresponds to a negative g max , one is always in case (a) and the rate Γ is given by Eq. (26).
Altogether, to capture the behaviour of the instability rate in the whole range of parameters, one should use Eqs. (24)- (25), which include all modes as well as the first correction due to the finite detuning. Note that the perturbative expansion yielding those expressions is expected to hold provided α q is not too large, i.e E av (q) is not too small [see below for a discussion of the breaking points of the approach].

C. Results based on the analytical approach
We now show the results obtained from this analytical approach for the instability rate, Eq. (25), with s q being generically computed to second order in the detuning (unless explicitly specified). The stability diagram is very similar to the one obtained numerically [see Fig. 2]. Let us investigate in more detail this agreement and comment on the small differences between the numerical and analytical diagrams.

Agreement Regions
For values of K/ω which are not too close to the zeros of the Bessel functions J 0 and J 2 [93] (i.e. not too close to the edges of the lobes in the stability diagram, zones A on Fig. 2), the analytical calculation gives a very good estimate of both the transition point and the instability rate, up to moderate interactions g [ Fig. 7]. At higher g, more harmonics become important (presumably, in a complex and "coupled" way: we verified that building independent effective models for separate harmonics does not reproduce the numerics accurately). Zones where J 4 vanishes are "favorable" since the first correction to the effective model [Eq. (21)] is absent, and the agreement between analytics and numerics survives for larger values of g.
In the vicinity of the "common" zeroes of J 2 (K/ω) and J 0 (K/ω) (i.e. "between" the lobes of the stability diagram, zones B on Fig. 2), both the analytics and numerics agree and predict a stable regime.

Disagreement Regions
A detailed description of the small disagreement regions is provided in Appendix D. In brief, the analytics breaks down in two main regions: (i) around the first zero of J 0 (i.e. the only zero of J 0 where J 2 is significant; zones C on Fig. 2), the perturbative approach breaks down due to the vanishing of the effective dispersion E av , and the analytical approach fails. (ii) Near the common zeros of J 0 and J 2 (i.e. "between" the lobes of the stability diagram), we find that the numerics and the analytics disagree in the way the stable zones close at large g. We identified two main causes for this effect: the breakdown of the perturbative approach due to the vanishing of E av , and the vanishing of the second harmonics near the zeros of J 2 . The consequences of these factors taken together are very involved, since they can both compensate for each other and compete. In brief, at the left border of each lobe (zones D on Fig. 2), the numerics displays a transition to instability (although showing up only at large g), which is not captured by the analytics. At the right border of the lobes (zones E), the instability predicted by the full numerics is most likely due to a joint effect of the second and higher-order harmonics, and is therefore not accurately captured by the analytical approach [see Fig. 21 Nevertheless, the disagreement regions constitute very narrow zones in the stability diagram. Figure 2 summarizes the zones where the analytical and numerical calculations agree (in green) and the ones where they disagree qualitatively (in red) or just quantitatively (in orange).

Scaling in limiting cases
The analytical solutions provide a general scaling behaviour for the instability rates as a function of the model parameters, which we analyze for three limiting cases.
a. Weak interactions: If ω is larger than the effective free-particle bandwidth of the model ("highfrequency" regime), the system is stable at g ≈ 0, and features a transition to an unstable phase at some finite g c , with a scaling [see Eq. (24)] Conversely, if ω is smaller than the effective free-particle bandwidth ("low-frequency" regime), the system is typically unstable at very small g, and [see Eq. (26)] Γ ∝ g.
b. Weak driving amplitudes: The analysis of the transition to the unstable phase from K = 0 (where the system is stable) to finite K is subtle, since this limit does not commute with the thermodynamic limit: indeed, one has to note that the resonance domain of each mode has a vanishing width when K → 0; therefore, for any finite-size system, instabilities at vanishing K occur only for discrete values of the parameters (fulfilling one of the resonance conditions associated with a specific mode); it is only when increasing K that those instability regions grow in parameter space, eventually merging and forming the first lobe of the stability diagram. In the thermodynamic limit, one nevertheless finds the general scaling Γ ∝ K 2 , which stems from the second order Bessel function which dominates at low K in Eq. (26).
c. Low driving frequency: In this limit, one finds [see Eq. (26)] where f is a generic function of K/ω only.

V. FINITE-SIZE SYSTEMS: FROM A DOUBLE WELL TO THE ENTIRE LATTICE
The analytical approach, which treats different modes independently, allows one to readily predict the behaviour of finite-size systems under periodic modulation, close to a parametric resonance. In this case, the "continuous" dispersion E av (q) is replaced by discrete energy modes. For a small number of lattice sites (and thus of modes), the total instability rate Γ = max q s q , which is given by the envelope of the curves associated with individual modes (see Fig. 8), may not be a monotonic curve as a function of g. Instead, it is rather composed of disjointed bumps if the resonance domains of two consecutive modes are disconnected. Figure 9 shows how individual modes contribute to the total growth rate when increasing the number of sites, both in and outside the lowfrequency regime. Interestingly, the situation is slightly different in these two cases. Outside the low-frequency regime, i.e. for ω J eff , the onset of instability is governed by the mode q = π [see Sec. IV], and therefore, the critical frequency, below which the system becomes unstable, happens to be the same for both an infinite and a two-site system. Since the mode q = π remains the "most unstable" one, all the way up to large values of g [see Fig. 8], the stability diagram of a two-site system is in fact very similar to the one associated with an infinite system (whereas a one-site system is trivially stable) [94]. Conversely, in the low-frequency regime (ω J eff ), instability occurs already at vanishingly small g, induced by a certain resonant mode q res fulfilling the resonance condition [see Sec. IV]. Increasing the number of sites therefore lowers the instability boundary as a function of g, since increasing the number of points in momentum space allows for the states to get closer to the maximally unstable mode. Moreover, if the number of sites N is too small, the discretisation in momentum space is so rough that the resonance domains of consecutive modes may be disconnected, resulting in "islands" of instability in the phase diagram. These islands begin to merge with increasing N , as the resonance domains of adjacent modes overlap, only gradually approaching the phase diagram of an infinite system. Both our analytical and numerical methods capture this behaviour, and agree for any  FIG. 9. Total instability rate and contributions of a few modes as a function of the interaction strength g, for ω = J ("lowfrequency" regime) and ω = 5J (outside the "low-frequency" regime), and for an increasing number of lattice sites N = 2, 4, 6, 64. The total instability rate is the envelope of the curves associated with "available" modes in the system (for N = 2, only q = 0, π; for N = 4, only q = 0, π/2, π, 3π/2...). The curve gets smoother when increasing the number of sites. Outside the low-frequency regime, the transition to instability is governed by the mode π and already captured by a two-site system. In the low-frequency regime, the instability rate at low g, which is governed by a certain mode qres, is probed with more and more accuracy when gradually increasing N : here, the 2-site system is always stable since qres > π, the 4-site system displays a transition since qres > π/2, and then this transition lowers when probing qres.
number of lattice sites. Interestingly, since the instability rate is determined by its maximal value over all modes (i.e. the envelope of all curves in Fig. 9), its estimate obtained for finite-size systems is already very decent, and further increasing the number of sites just smoothens the curve of the instability rate as a function of g.
We point out that for very small systems, there might be minor corrections due to edge states: indeed, the previous treatment, formulated in momentum space, tacitly implies periodic boundary conditions. However, in the lab frame, the drive actually breaks translational invariance and the boundary conditions are intrinsically open, possibly yielding different selection rules for states near the edges. Although such a difference is expected to play a negligible role in large systems, it may lead to minor but noticeable deviations in small systems.

VI. THE EFFECT OF TRANSVERSE DIRECTIONS: 2D LATTICES VS TUBES
So far, our study focused on the origin of parametric instabilities that occur in a driven system, which satisfies the 1D GPE on a lattice (2). As further discussed in Sec. VII, this model can be used to describe the physics of weakly-interacting bosonic atoms, trapped in a 1D optical lattice. However, from an experimental point of view, it is intriguing to determine the effects of transverse directions on the stability diagram and instability rates. Indeed, optical-lattice experiments involving weakly-interacting bosons [39,41] typically feature continuous transverse degrees of freedom, commonly referred to as "tubes" or "pancakes". Besides, experiments realizing artificial magnetic fields [95,96] involve twodimensional optical lattices, and hence, it is relevant to study the fate of instabilities as one transforms a 1D optical lattice into a full 2D lattice, by adding sites (and allowing for hopping processes) along a transverse direction.
In this section, we extend our previous analysis to study the effects of (i) a secondary tight-binding-lattice direction (resulting in a ladder or a full 2D lattice), and (ii) an additional continuous ("tube'") degree of freedom. In the following, we assume that the periodic modulation remains exclusively aligned along the original tightbinding lattice dimension.

A. Two-dimensional lattice geometry
In this Section, we consider the addition of a transverse lattice direction, aligned along the y-axis; by doing so, we keep the same time-dependent modulation as before, which is hence exclusively aligned along the x-axis. Our aim is to study how an increase in the number of sites along the transverse direction affects the instability rates, previously evaluated for the 1D configuration. Theoretical studies [63,65] have reported that heating and collisional processes are enhanced by the presence of a transverse direction, which provides crucial information for current experimental studies. For this extended model, the time-dependent GPE reads i∂ t a m,n = − J(a m,n+1 + a m,n−1 + a m+1,n + a m−1,n ) + K cos(ωt)ma m,n + U |a m,n | 2 a m,n , where each site of the underlying (2D) square lattice is now labelled by two integers (m, n). As for the 1D case [Sec. II], the condensate wavefunction at time t = 0 is again given by a Bloch state e i(pxm+pyn) , with momentum p x = 0 if J 0 (K/ω) > 0 and p x = π if J 0 (K/ω) < 0; for this choice of the time-modulation, the momentum along the y direction is necessarily p y = 0. The timedependent Bogoliubov-de Gennes equations, which govern the time evolution of the mode q = (q x , q y ), take the form where which is a straightforward generalization of Eq. (9). Therefore, both the numerical and the analytical methods presented in the previous sections can be readily applied to Eq. (30).  Figures 10 and 11 show the stability diagrams, obtained numerically by increasing the number N y of transverse lattice sites (from 4 to 16); they correspond to the drive frequency ω = 10J (i.e. "high-frequency" regime, where ω is greater than the effective bandwidth) and FIG. 11. Numerical instability rate as a function of interaction strength g = U ρ and modulation amplitude K/ω, for ω = 5J (i.e. in the "low-frequency" regime) and for two different number of lattice sites in the transverse direction: Ny = 4 (top) and Ny = 16 (bottom). ω = 5J (i.e. "low-frequency" regime), respectively. Since the motional degrees of freedom along the x and y directions are independent of each other, the problem being separable, the situation is analogous to that previously discussed in the case of finite-size systems.
Let us summarize the results: away from the "lowfrequency" regime [ Fig. 10], namely for ω > 4(J + J eff ) [the drive frequency is again compared to the modified (effective) free-particle bandwidth], the system is stable for small g and the onset of instability is governed by the mode of maximal energy (i.e. the first mode to exhibit a resonance), which now corresponds to q = (π, π). Therefore, as soon as two sites are present in the transverse direction, the transition point to the unstable regime is well-captured. Moreover, both the instability boundary and the instability rates in its vicinity are expected to remain unaffected when increasing N y , since they are dominated by the mode (π, π) [see Eq. (27)]. Conversely, in the low-frequency regime, there is an instability at vanishingly small g, which is driven by a certain mode (q x , q y ): in this case, adding the number of transverse sites lowers the instability boundary, as increasing the resolution in momentum space allows one to get closer to this maximally unstable mode.
Using Eq. (30), one can readily extend the analytical approach of Sec. IV. Indeed, all the results of Sec. IV remain applicable in the present case, provided that the dispersion E av (q) is now replaced by The resulting stability diagrams are shown in Fig. 12 for the same parameters as in Fig. 10, and they reveal an excellent agreement with the latter. Interestingly, the agreement between analytics and numerics is even better than for the 1D configuration, since the presence of transverse modes in the effective dispersion Eq. (31) prevents E av (q) from vanishing when J 0 (K/ω) = 0 [97]. Analytical instability rate as a function of interaction strength g = U ρ and modulation amplitude K/ω, for ω = 10J (i.e. outside the "low-frequency" regime) and for two different number of lattice sites in the transverse direction: Ny = 4 (top) and Ny = 16 (bottom). To be compared with Fig. 10.

B. Continuous transverse degrees of freedom: the case of "tubes"
We now consider the case where the transverse degree of freedom corresponds to free-motion along a continuum, in one or more transverse directions. In this case, the Bogoliubov-de Gennes equations still take the form of Eq. (30), now with so that the previous analysis still holds, provided that the time-averaged Bogoliubov dispersion E av (q) is now replaced by Interestingly, since this dispersion is unbounded, ω will always be smaller than the total bandwidth, and hence, there will always be (at least) one mode that is precisely set on resonance: in other words, the system necessarily falls into the "low-frequency regime" previously introduced, where it is unstable at any finite g. As explained in Sec. IV B, the instability rate can be evaluated, at lowest order, by calculating the rate on resonance associated with the mode(s) satisfying the condition ω = E av (q res ) [see Eq. (26)]: Importantly, although there are generically several resonant modes q res (corresponding to different q x and q ⊥ ), they do not have the same instability rate s * q res : the total instability rate Γ is then attributed to the most unstable mode, namely where we introduced the notation q mum to denote the momentum of the most unstable mode. We then identify two cases: (i) If ω > 4|J eff |(4|J eff | + 2g), which is often the case in realistic configurations, there exists a resonant mode at q mum x = π, which is thus the most unstable one; q mum ⊥ is then adjusted so as to respect the resonance condition. The maximally unstable mode thus reads In this case, the total instability rate is given by the simple analytical formula (ii) If ω < 4|J eff |(4|J eff | + 2g), the most unstable mode, which obeys the resonance condition ω = E av (q), is necessarily reached at q mum In this case, the total instability rate is given by Interestingly, these results do not depend on the number of transverse dimensions, since all that matters is the existence of an unbounded continuum. Figure 13 compares the numerical stability diagram to the analytical formula Eq. (36), for realistic (experimental) parameters. As visible on Fig. 14, which shows cuts through the diagrams of Fig. 13, the agreement is excellent at low interaction/modulation amplitude, and the scaling of the instability rate is remarkably simple: Γ increases linearly with g and quadratically with K. Such simple results could promisingly be compared with experiments. Note that in experiments, one typically has J/ ≈ 100Hz − 1kHz, so that the inverse instability rates are typically of the order of 1 − 10ms.
At this stage, it is worth comparing those results with the more traditional Fermi-Golden-Rule (FGR) approach: first, in the present case, and contrary to what is expected for a FGR regime [62], the final  FIG. 14. Top: Instability rate as a function of the interaction strength g for ω = 32J and for two values of the modulation amplitude. As captured by the analytical expression Eq. (36), the instability rate increases linearly with g. Bottom: Instability rate as a function of the modulation amplitude K/ω for ω = 32J and for two values of the interaction parameter g. As captured by the analytical expression Eq. (36), the instability rate increases quadratically with K/ω at low amplitudes.
rate can be inferred from the contribution of a single (well-identified) mode, and not from a sum over all modes (which would then involve the density of states in the description). Second, the form of the final formula in Eq. (36) is significantly different from that resulting from a FGR argument, where the rate would typically scale as g 2 (as dictated by the squared matrix element of the perturbation [62]). When considering a driven weakly-interacting degenerate Bose gas, one expects the short-time dynamics to be dominated by the parametric instability identified in this work; at longer times, as the condensate significantly depletes, heating rates are expected to be dominated by a FGR behavior [62,63], not captured in the present analysis. We believe that those two distinct mechanisms and instability regimes could be probed by current ultracold-atom experiments.
As a practical remark, we emphasize that the results presented in this section suggest that working with a lattice seems to be more favorable than with tubes/pancakes, in view of reducing parametric instabilities in periodically-driven Bose systems.

VII. THE PERIODICALLY-DRIVEN BOSE-HUBBARD MODEL
In this section, we discuss the relation between the Gross-Pitaevskii equation introduced in Eq. (2) and the (full quantum) driven Bose-Hubbard model [10]. Our goal is to clarify the validity of our analysis in view of describing parametric instabilities and heating in the context of this quantum model.

A. Effective Hamiltonian and micromotion operator
For the sake of concreteness, here we focus on the 1D model of Eq. (2), but the discussion applies to all models investigated in the paper (in particular the twodimensional models of Section VI). Consider a system of weakly-interacting bosons, trapped in a shaken 1D lattice, as described by the periodically-driven Bose-Hubbard Hamiltonian [10] where J > 0 denotes the tunnelling amplitude of nearestneighbour hopping, and U > 0 is the on-site interaction strength. The on-site potential term describes a timeperiodic modulation of amplitude K and frequency ω = 2π/T . As we stated already in the introduction, the dynamics is stroboscopically governed by the effective Floquet Hamiltonian, see Eq. (1). In the absence of interactions, it is exactly [98] given by [10][11][12] Besides, the kick operator [Eq. (1)] reads [16] K kick (t) = i log R (t)e −iK rot kick (t) , whereR(t) = e −iK/ω sin(ωt)X is the unitary transformation to the rotating frame (withX ≡ Σ n nâ † nân the position operator on the lattice). HereK rot kick denotes the kick operator in the rotating frame [16] , which is explicitly given bŷ Whenever interactions are present, Eq. (40) is no longer exact, and the effective Hamiltonian is a much more complicated (possibly nonlocal) object [59,99]. In the high-frequency regime, it can be approximated using the inverse-frequency expansion [11], which to lowest order yieldŝ Hence, in the infinite-frequency limit, the periodic drive merely renormalizes the hopping matrix element J → J eff . Due to the Bessel function taking both positive and negative values, it is possible to tune the amplitudeto-frequency ratio such that J 0 (K/ω) = 0, in which case tunneling is completely suppressed [22,23,29]. Therefore, even a weakly-interacting lattice system can effectively behave as a strongly interacting one under periodic driving, in the sense that the interaction strength U potentially dominates over the tunneling when J eff ≈ 0. We recall that the lattice dispersion flips sign when J 0 (K/ω) < 0, in which case the stable minimum appears at quasi-momentum q = π: namely, at equilibrium, bosons are expected to condense in a finite-momentum state in this regime. We point out that corrections to the Hamiltonian (43) become non-negligible away from the high-frequency limit, and that these are expected to manifest themselves over longer time scales.

B. Relation to the mean-field approach and observables
The mean-field approach presented in this paper [based on Eq. (2)] is expected to provide a good description for the shaken 1D Bose-Hubbard model [Eq. (39)] at short times, as far as the weakly-interacting regime is concerned. More specifically, if one assumes that the initial state is that of fully condensed bosons, the system in Eq. (39) can be treated within mean-field theory [100]. In this approximation, the annihilation and creation operatorsâ n ,â † n are merely replaced by classical fields a n , a * n , and the Heisenberg equations of motion forâ n ,â † n lead to the time-dependent GPE in Eq. (2).
Throughout this work, our choice for the initial state corresponds to the Bogoliubov ground-state of the effective Hamiltonian (43) [see Section II]. In order to compute heating rates for such a model, following the procedure detailed in Sec. III B, we also had to fix the initial condition for the fluctuation term δa, which in our GPE-based approach, could be taken to be an arbitrary small perturbation on top of the condensate.
However, our choice for the initial condensate wave function (i.e. ground state ofĤ eff ) suggests a natural choice for the initial fluctuation term δa(t = 0): indeed, one could assume that at t = 0, this fluctuation term is precisely given by the Bogoliubov wave function of a condensate in the ground-state ofĤ eff , as obtained from a numerical diagonalization of the corresponding (effective) Bogoliubov equations, namely, where u q , v q correspond to the solution of A relevant observable to quantify heating and losses, which can be probed in experiments, is the noncondensed density, which, in the Bogoliubov approximation, reads [101] ρ nc (t) = 1 with Vol the volume of the system. Given the above initial condition, it is straightforward to apply the numerical tools developed in Sec. III B [102], so as to describe the time evolution of the non-condensed density Eq.(45) in the system. For instance, Fig. 15 shows the exponential growth rate of the non-condensed fraction (referred to as the loss rate) obtained for the 2D lattice model of Sec. VI A, which is therefore expected to describe particle losses in the associated 2D driven Bose-Hubbard model. Note that, unsurprisingly, this is very similar to the stability diagram Fig. 11(b), which corresponds to the same parameters [see Sec. III B]. The aforementioned aspects are intrinsic to the meanfield Bogoliubov approximation in the frame of the model under consideration. More generally, there are a few important points one has to keep in mind, when adopting a mean-field approach to study out-of-equilibrium ergodic (non-integrable) bosonic systems subject to periodic driving: • As already alluded to above, the Floquet Hamiltonian becomes an increasingly nonlocal operator, the further the drive frequency deviates from the infinite-frequency limit. An intriguing and interesting part of this intrinsic nonlocality is due to Floquet many-body resonances [103] appearing because of hybridisation of many-body states induced by the drive. These resonances appear as a result of energy absorption and provide shortcuts between states in energy space separated by multiples of the drive frequency. Hence, in the unstable phase, they are expected to affect the heating rates at times t ∼ U −1 , when interaction effects between quasiparticles become important.
• In fact, one can anticipate (using the inversefrequency expansion in the rotating frame) that the Floquet Hamiltonian of the driven Bose-Hubbard model contains complicated three-and higher-body interaction terms, starting from order ω −3 . It is currently an open problem how these terms modify and limit the application of Bogoliubov's meanfield theory associated with the effective (Floquet) Hamiltonian.
Despite the open character of these potential issues related to interacting bosonic systems, we expect that our GPE-based analysis captures the dominant contribution to the short-time evolution of the periodically-driven Bose-Hubbard model in Eq. (39), in particular the onset of instability.
C. Time evolution beyond the linearised regime: the conserving approximation Since the mean-field Bogoliubov approximation assumes a macroscopic occupation of the condensate mode, it remains valid provided the number of non-condensed atoms remains small compared to the total number of atoms throughout the entire subsequent evolution. As expected, this assumption fails in the unstable regions, where the depletion of the condensate grows exponentially and the mean-field approach typically holds at short times. Thus, the time scale for reaching a sizeable dynamical depletion sets a natural upper bound on the validity of mean-field approaches.
One way to avoid the problems associated with particle conservation, is to apply the weak-coupling conserving approximation (WCCA) [73]. Based on a Keldysh field theory formalism, the WCCA is the minimal extension of Bogoliubov theory, which includes the proper effective interactions between the Bogoliubov quasiparticles and the condensate, order by order in the original interaction strength U , while ensuring particle conservation at all times. Truncated to linear order in U , the WCCA is equivalent to the Bogoliubov-Hartree-Fock approximation [104], and in the following we shall restrict to this case. The unknown variables in the WCCA formalism are [105] where φ(t) = a 0 n (t)e −i t 0 µ(t )dt √ Vol is the q = 0 mode of the spatially-homogeneous (i.e. nindependent) condensate wave function, and where µ(t) = −zJ cos(K/ω sin ωt) + g is the chemical potential, which is irrelevant for U(1)-conserving dynamics. The subscript c stands for connected correlators: . The equal-time correlator F 11 (t; q) is, apart from an additive constant, the phonon density. In momentum space, it is related to the momentum distribution function of the quasiparticles n q (t) by n q (t) = F 11 (t; q) − 1/2. Note that n q (t) does not include the condensate delta function peak at q = 0. The WCCA EOM for the equal-time correlators represent a simple system of non-linear, non-local in space equations where * denotes complex conjugation.
The conserved quantity itself is the total number of particles (condensed and excited): One recognises the GPE equation for the spatially homogeneous condensate density φ(t), and identifies the additional phonon feedback terms. Quite generally, if one eliminates all terms containing integrals over q, the WCCA EOM decouple into the familiar GPE for the condensate, see Eq.
Opening up a channel between the quasiparticles and the condensate, we find that the phonon build-up rate decreases (and consequently, the condensate depletion slows down) compared to the linearised BdG equations, in agreement with intuitive expectations. In other words, the WCCA deviates from the exponentially growing linearised solution, which sets a natural scale on the validity of the mean-field approximation.
To justify the applicability of the mean-field approach at short times, we compare the dynamical build-up of quasi-particle excitations in the parametrically unstable regime, predicted by the linearised EOM [see Eq. (9)], and the WCCA (to order U ). Here, we initialise the system in the Bogoliubov ground state corresponding to the non-driven Hamiltonian H(K = 0). Figure 16 shows the time evolution of the quasiparticle excitations (top) and the condensate density (bottom) for the 2D model of Sec. VI B. One can clearly identify two stages of evolution: in stage I, which continues for about t 220T driving periods, the dynamics is governed by the singleparticle mean-field physics. The dotted vertical line at t * ≈ 60T marks the time up to which the dynamics predicted by the linearised mean-field equations agrees with the WCCA. Thus, for t ≤ t * , this initial regime features an exponential build-up of "phonons" at the momenta satisfying the resonant condition ω = E av (q). For 60T t 220T , the time evolution is characterised by a different growth rate than the one predicted by the mean-  Fig. 16. Here, we set the lattice constant ax = 1 to unity, while ξ = 1.47ax denotes the BEC healing length.
field equations for this system configuration. By the time t ≈ t * , the back-action effect of the quasiparticles onto the condensate becomes sizeable and the dynamics no longer follows the initial linearised exponential growth. Based on the available data, it is not possible to determine whether the growth remains exponential or follows another law. In stage II, t 220T , the population of the resonant q-modes slows down significantly, although it never really saturates, as can be seen from the condensate depletion at longer times, see Fig. 16 (bottom panel). Due to the unbounded character of the dispersion along the q y -direction, the presence of resonant modes set by the condition ω = E av (q) does not allow for the formation of a Floquet steady state. Nevertheless, the growth rate diminishes significantly and the dynamics slows down, as expected for a pre-thermal regime. Interestingly, the con-densate evolution curve [ Fig. 16 (bottom panel)] changes curvature in stage II, compared to stage I. Curiously, we note that a very similar change of behaviour was reported in a cosmological context [76].
Important features of the quasiparticles dynamics are conveniently reflected in their momentum distribution function n q (t) [107]. Figure 17 shows a comparison between the momentum distributions, as obtained by timeevolving the linearised and WCCA EOM over a long duration (t = 1000T ), for the same 2D model. The visible lines, which materialize largely-populated modes, correspond to on-resonance modes [E av (q) = ω], which have indeed been shown to be unstable. In the linearised dynamics, the largest peaks are visible around the maximally unstable modes, as predicted by our theory [see Eqs. (35) and (37)]; for the considered parameters, the most unstable mode indeed corresponds to q mum x = ±π and q mum y is adjusted so that the resonance condition ω = E av (q) is fulfilled [see Sec. VI B]. We stress that the large width of the peaks around q x = ±π, in the top panel of Fig. 17, is an artifact due to the significant population of momentum modes at long times, which yields saturation on this plot; we refer to Fig. 1 for a more striking illustration of these peaks, as calculated at shorter times (for the same model parameters). In contrast, the WCCA couples all momenta through the q -integrals in Eq. (47), allowing the population to further spread along the xaxis. All modes fulfilling the resonance condition now exhibit peaks of similar magnitudes, precisely materializing the curve of equation ω = E av (q); note though that the populations of the modes along this curve is not completely homogeneous, presumably due to the non-linear character of the WCCA EOM and/or the differences in the instability rates of resonant modes. Moreover, as we discussed in Sec. IV C, modes which are not precisely on resonance, but sufficiently close to it, are also affected by the instability.
A summary of the various behavior types of the momentum distribution, which are expected for the different parameters regimes, are presented in the concluding Section VIII, see Figs. 19-19 and Table II.
We conclude this Section by noticing that while the WCCA to leading order in U features a pre-thermal Floquet regime at high driving frequencies [73], the absence of quasiparticle collisions prevents the onset of thermalising dynamics at later times. Including the next-toleading order in U has been done for the periodicallydriven O(N ) model [109,110], where the crossover from a pre-thermal steady state to a thermal regime becomes clearly visible.

VIII. SUMMARY OF THE MAIN RESULTS AND CONCLUDING REMARKS
The analysis presented in this work reveals the crucial role played by parametric instabilities in the short-time dynamics of periodically-driven band models. It provides a quantitative description and a qualitative understanding of the stability diagrams of these periodically-driven systems, capturing both the stability boundaries and the instability rates in its vicinity. This work also identifies clear signatures of these early-stage instabilities, which could be detected in current cold-atom experiments.
On the one hand, the quantitative description of the instability rates will prove helpful to guide experimentalists in finding stable regimes of operation. A first indicator is the instability rates themselves. Table I summarizes the scaling behaviour found for instability rates in the thermodynamic limit as a function of the model parameters, as derived in Sec. IV C. Remarkably, these scaling laws hold in all cases investigated in this work, i.e 1D/2D lattices and high/low-frequency regime.
Limiting scalings of instability rates in the thermodynamic limit, in terms of the interaction strength g, the modulation amplitude K and frequency ω. The "high" and "low" frequency regimes are defined through a comparison between ω and the free-particle effective bandwidth, see Sec. III A.
A second important information concerns the identification of the most unstable mode q mum . Table II recalls the characteristics of this most unstable mode for the various relevant regimes discussed in this work. In the "highfrequency" regime, where ω is larger than the free-particle effective bandwidth, the onset of instability is necessarily driven by the mode q = π (or q x,y,z = π if several lattice degrees of freedom are present). In the "low-frequency" regime, the most unstable mode q mum obeys the resonance condition ω = E av (q). Yet, in dimensions higher than one, we found that this condition is not sufficient to unambiguously determine q mum , since a whole set of resonant modes (associated with different rates) typically satisfy it. Our analysis showed that two situations can occur, and that these could be distinguished by comparing the frequency ω to the effective Bogoliubov bandwidth associated with the shaken-lattice direction alone: if ω > 4|J eff |(4|J eff | + 2g), then q mum x = π and q mum ⊥ = 0 [see Eq. (35)]; conversely, if ω < 4|J eff |(4|J eff | + 2g), we find q mum x < π and q mum ⊥ = 0 [see Eq. (37)]. Let us emphasize that, in principle, the so-called "low-frequency" regime can be reached at arbitrarily large frequencies; for instance, this typically occurs whenever the free-particle effective dispersion is unbounded, as in the case where a continuous direction is present. Figures 18 and 19 show the momentum distribution as calculated for the hybrid-2D band model of Sec. VI B (i.e. a shaken lattice along the x-direction and a continuum along the y-direction); by exploring the different parameter regimes, these results illustrate all the behaviors described above [Table II].
low-frequency regime high-frequency regime q mum ⊥ > 0 (at the onset of instability ) see Eq. (37) see Eq. (35) TABLE II. Characterization of the most unstable mode (mum) for various relevant regimes. The "low-frequency" regime is generically defined through the criterion according to which the drive frequency ω is smaller than the effective free-particle bandwidth (e.g. ω < 4|J eff | for the 1D shaken-lattice model of Section III A). We stress that the effective free-particle dispersion is unbounded for models involving a continuous spatial dimension, in which case the system necessarily corresponds to the "low-frequency" regime (for all ω).
Those predictions on the instability rates can readily be tested in present-day experiments. For typical experimental parameters, the ratio ω/J is generally set in the range 10−100; therefore, in the case of pure-lattice systems, experiments would typically lie in the "highfrequency" regime defined in Section III A, and the most unstable mode is thus expected to be q mum x,y,z = π. In turn, experiments with a continuous degree of freedom fall into the second column of Table II, and we thus expect the maximal instability to occur through q mum x = π and finite q mum ⊥ . We point out that experiments could also investigate the "low-frequency" regime of driven-lattice systems (first column of Table II), by working with lower drive frequencies ω ∼ J (see, e.g., Ref. [28]).
The different observables that we introduced are accessible in current experiments and probing them should permit to have clear experimental signatures of the predicted instability behaviour. In particular, measuring the momentum distribution of the gas, which is routinely performed through time-of-flight techniques, and identifying the positions of the peaks, would give access to the momentum of the most unstable modes, which is one of the most recognizable characteristics of parametric instabilities.
On the other hand, several important conceptual points arise from the detailed analysis presented in this paper: in particular, the intuitive stability criterion ω > 2W eff is not sufficient to properly estimate the instability boundary and rates in the system, as becomes clear from the rigorous derivation we provided. More precisely, the parametric resonance condition for the instability is fulfilled for ω = E av (q), but in addition, one has to take the (possibly) large width of the resonance window into account in order to properly estimate the position of the boundary; we note that similar effects were found in the Quasiparticle momentum distribution function nq after a duration of t = 10T , for the 2D band model from Sec. VI B (with a lattice along the x-direction and continuum along the y-direction), computed using the linearised dynamics, for different values of the relevant parameter 4|J eff |(4|J eff | + 2g)/ω; see Table II. The suggested ellipses correspond to on-resonance modes [Eav(q) = π], which are unstable, while the most unstable modes, which are revealed by the main peaks in the momentum distribution, are accurately described by our theoretical predictions [see Table II].
As a result of our study, the following general guidelines deserve to be highlighted: in lattice systems, working at high frequency (much larger than the effective bandwidth of Bogoliubov excitations) is favourable, since 19. Quasiparticle momentum distribution function nq for the 2D band model from Sec. VI B (with a lattice along the x-direction and continuum along the y-direction), as a function of the parameter 4|J eff |(4|J eff | + 2g)/ω and of qx (upper line)[resp. q ⊥ , (lower line)], after integration over q ⊥ (resp. qx), and computed : (left) after a time t = 10T using the linearised dynamics; (middle) same time but using the WCCA; (right) after a longer time t = 101T using the WCCA. The dotted blue line indicates our theoretical prediction for the most unstable mode [see Table II], which is verified in all regimes: for 4|J eff |(4|J eff | + 2g)/ω > 1, we find q mum ⊥ = 0 and q mum x < π, as predicted by Eq. (37); for 4|J eff |(4|J eff | + 2g)/ω < 1, we find q mum x = π and q mum ⊥ > 0, as given by Eq. (35); let us note that, in the latter regime, the peaks are only apparent after a long time (here 101T , right panel) due to the small values of the corresponding instability rates; see also Fig. 18(a), where these peaks are made visible. As already observed in Fig. 16, both linearized and WCCA agree very well on the position of the peaks, although the WCCA predicts a broadening of the instability zone due to non-linear couplings between the modes. Note that at longer times (right panel), higher order parametric resonances are also visible. The apparent "quantised" character of the resonant peaks for t = 101T is a numerical artifact due to the discretization of the variable 4|J eff |(4|J eff |+2g)/ω used in the numerical simulation.
no parametric resonance can then occur, ensuring a stable evolution. We also analyzed the effect of transverse directions, be it a deep lattice or tubes; although the instability is generically increased in both cases by the fact that more and more modes are available in the system (offering more possibilities to have resonances), we showed that working with a lattice should be much more favorable to find stable regimes and have less heating. We stress that our analysis results from a tight-binding approximation, namely, the restriction of the study to a single Bloch band. This simplification is reasonable for situations where the drive frequency is set away from any inter-band resonances, which is indeed possible in experiment. Interestingly, we note that such inter-band excitations have been generated and studied experimentally in Ref. [111]. The interplay between inter-band transitions, as generated by setting the drive frequency close to resonance with inter-band spacing, and the existence of parametric instabilities, should show rich behaviors, which could be explored through the theoretical approach presented in this work.
Finally, we emphasize that the instability rates derived from parametric resonances (see, e.g., the scalings in Table I) differ quantitatively from the ones computed within a Fermi Golden Rule (FGR) approach [63], and which generically scale as g 2 for low interactions. In particular, we note that contrary to the case of parametric resonances, the rates emanating from the FGR approach strongly depend on the dimensionality of the system, through its density of states, which indicates that the two phenomena are indeed very distinct in essence. Since the system is bosonic, we expect that the parametric instability should dominate the short-time dynamics, due to its exponential growth rate. It is only at later times that Fermi's Golden Rule comes to rule the behaviour of the system, when the dynamics will be dominated by collision processes leading to thermalisation at a higher temperature, a feature which is generically known for quantum Bose fluids [76]. Hence, both phenomena take place in different regimes of the time evolution. More precisely, for typical experimental parameters, the instability rates derived from parametric resonances are found to be of the order of 1 − 10ms, while the Fermi Golden Rule decay usually spreads over times of the order of 100ms [63]. It is within the scope of present-day experiments to observe those two distinct regimes.
where the discrete operatorL is defined byL(·) n ≡ (·) n+1 + (·) n−1 . The equations of motion (A4) are still explicitly time-dependent, and hence, they do not allow for a direct identification of the quasiparticle spectrum.
In the static case, a common way of solving this issue is to use a slightly different parametrization, namely a n (t) = e −iµt [a (0) n (t)+δa n (t)], where we have factored out the dynamical phase µ associated with the free evolution. Using this parametrization, we obtain the Bogoliubov equations in the form where∆δa n ≡ δa n+1 − 2δa n + δa n−1 denotes the discrete Laplacian. The resulting equations of motion correspond to the standard form of the BdGEs, which are indeed time-independent (the spurious complex off-diagonal terms have disappeared), allowing for a direct extraction of the Bogoliubov spectrum.
In the driven-system situation, the time-dependent phase of the BEC wavefunction a n (t) does not reduce to a trivial phase µ, as it also contains an additional dynamical phase induced by the periodic driving. Therefore, factorizing e −iµt alone in the Bogoliubov parametrization would not be sufficient to get rid of all the spurious terms. In order to solve that issue, we generalize the method discussed above by factorizing the entire time-dependent phase of a where we have introduced the modulus and phase of the BEC wavefunction, a n (t) = ρ n (t)e iφn(t) . For the model considered in this work, translational invariance is recovered in the rotating frame, which allows one to further simplify the parametrization: indeed, taking into account the fact that ρ n (t) is uniform (it does not depend on n), we can equivalently write Eq. (A6) in the form [Eq. (4)] a n (t) = a (0) n (t)[1+δa n (t)], where we have now factorized the full condensate wavefunction. Linearizing the Gross-Piteavskii equation in δa n then yields the time-dependent Bogoliubov-de Gennes equations, which take the general form considered in the main text i δ a ṅ δa * n = F n (t) U |a (0) n | 2 −U |a where we introduced the operator F n (t), whose action on δa n is defined by F n (t)δa n ≡ − JL (a (0) δa n ) n a (0) n + 2U |a (0) n | 2 δa n + K cos(ωt)nδa n − iȧ Compared to Eq. (A2), the differences are twofold. First of all, the off-diagonal terms are real in this approach, which is numerically more practical. Then, we note that the last term in Eq. (A9), which involves the logarithmic derivative of the BEC wavefunction a n (t), features a time-dependent chemical potential, which takes into account the entire dynamical phase induced by the periodic driving; in particular, we recover the chemical potential term −µ = −g + 2J in the static limit K → 0.
For the reasons detailed above, we have chosen to consider the convenient parametrization defined in Eq. (A7), both in our numerical and analytical extraction of the instability rates; see also Eqs. (47) for the extension to WCCA. We stress that the parametrization in Eq. (A7) should be revised for more elaborate models (multiband, multicell...), where ρ n (t) is not uniform, using the more general form (A6). Finally, we emphasize that all parametrizations lead to the same physical results, and that choosing one or the other is just a matter of technical convenience.
where ε(q, t) = 4J sin q 2 sin q 2 + p 0 − K ω sin(ωt) , one first introduces the change of basis which diagonalizes the time-averaged part of the equation, where cosh(2θ q ) ≡ (4|J eff | sin 2 (q/2) + g)/E av (q) and sinh(2θ q ) ≡ g/E av (q) with E av (q) = 4|J eff | sin 2 (q/2)(4|J eff | sin 2 (q/2) + 2g) the time-averaged Bogoliubov dispersion. Defining then a second transformation the Bogoliubov equations take the final form where sinh(2θ q ) = g/E av (q), E av (q) is the Bogoliubov dispersion associated with the effective GPE, within the Bogoliubov approximation [see Eq. (12)], W q (t) is a diagonal matrix of zero average over one driving period, with J l (z) the l-th Bessel function of the first kind.

Appendix D: Disagreement regions
The analytical approach is based on a perturbative solution in the detuning and the amplitude of the perturbation, which is performed on an effective model restricted to the second harmonic [see Sect. IV B]. We observe that it breaks down in two main regions: 1. The main discrepancy with the numerical results occurs around the first zero of J 0 (the only one where J 2 is not small as well, zones C on Fig. 2), the analytics fails to capture the numerically observed instability, and we also notice a completely different behavior between the first order and second order analytics [see Fig. 20]. This is coming from the fact that in this case the effective dispersion E av is flat, and the perturbative approach breaks down. More precisely, the first order analytical solution gives a resonance domain centered on ω = E av (q) ≈ 0 but of infinite width (∝ E av (q) −1 ), displaying thus a strong instability. The second order correction displaces the center of the resonance domain as ∝ E av (q) −3 , sweeping away all instability in the considered parameter range. More generally, the perturbative approach does not converge and cannot capture the observed behavior (which seems though to be due to the second harmonic alone). 2. Slight discrepancies also occur "between" the lobes of the stability diagram, as the numerics and the analytics disagree on the way those "black zones" close (or not) at large g. Firstly, near the zeros of J 0 (i.e. at the left border of each lobe, zones D on Fig. 2), the numerics displays a transition to instability which is not captured by the analytics (similarly to what was observed in case C). However, this failure of the analytics is less marked than in case C (in the sense that it appears for much larger g) for two main reasons: on the one hand, the fact that J 2 (K/ω) is also small can partly kill the divergences due to the vanishing of J 0 (K/ω); on the other hand, the stable regime extends here to much larger g coherently with the fact that the second harmonic has a very small amplitude (indeed, as visible on Fig. 21(a), the numerical instability seems rather due to higher harmonics). More precisely, at points where J 2 (K/ω) = 0, the analytics, which is essentially built on a secondharmonic-based model, is always stable and cannot reproduce the numerical instabilities due to higher harmonics (which nevertheless show up only at large g). Near such points (zones E), the numerical instability seems due to a joined effect of the second and next harmonics, and is therefore not accurately captured by the analytical approach [see Fig. 21(b)].