Parametric Instabilities in Resonantly-Driven Bose-Einstein Condensates

Shaking optical lattices in a resonant manner offers an efficient and versatile method to devise artificial gauge fields and topological band structures for ultracold atomic gases. This was recently demonstrated through the experimental realization of the Harper-Hofstadter model, which combined optical superlattices and resonant time-modulations. Adding inter-particle interactions to these engineered band systems is expected to lead to strongly-correlated states with topological features, such as fractional Chern insulators. However, the interplay between interactions and external time-periodic drives typically triggers violent instabilities and uncontrollable heating, hence potentially ruling out the possibility of accessing such intriguing states of matter in experiments. In this work, we study the early-stage parametric instabilities that occur in systems of resonantly-driven Bose-Einstein condensates in optical lattices. We apply and extend an approach based on Bogoliubov theory [PRX 7, 021015 (2017)] to a variety of resonantly-driven band models, from a simple shaken Wannier-Stark ladder to the more intriguing driven-induced Harper-Hofstadter model. In particular, we provide ab initio numerical and analytical predictions for the stability properties of these topical models. This work sheds light on general features that could guide current experiments to stable regimes of operation.


I. INTRODUCTION
Driving quantum systems periodically in time has been proposed as a versatile tool to generate unusual quantum phases of matter [1][2][3][4][5][6]. In the context of ultracold quantum gases, it was shown that subjecting neutral atoms to an external time-periodic drive could be used to design artificial gauge fields in these systems [6][7][8][9][10], hence opening promising perspectives in the quantum simulation of topological states of matter [11][12][13] and quantum magnetism [10,14,15]; see also the recent work [16]. The underlying concept of Floquet engineering [4,[8][9][10]13] builds on the fact that the dynamics of periodicallydriven systems can be well described by a static effective Hamiltonian, whose properties can be suitably designed by tailoring the driving protocol. On the experimental side, this idea has been extensively used to design artificial gauge fields for neutral atoms [11,12,[17][18][19][20][21][22], some of which lead to non-trivial geometrical and topological properties [13].
A particularly interesting class of periodically-driven setups is that featuring resonant time-modulations [23][24][25][26][27][28], in which the driving frequency ω resonates with an energy separation ∆ ≈ ω that is inherent to the underlying static system. In particular, such schemes can be exploited to finely control the tunneling matrix elements connecting neighboring sites of a lattice, and can be simply implemented by resonantly modulating a superlattice or a Wannier-Stark ladder; see Refs. [29,30] for experimental realizations and Ref. [10] for a review. This so-called photon-assisted tunneling effect constitutes a natural ingredient for the generation of artificial fluxes within cold-atom systems [23,[25][26][27], as exemplified by the recent experimental realizations of the Harper-Hofstadter model [17][18][19][20][21]31], a 2D lattice penetrated by a uniform magnetic flux [32]; the latter band model is of particular interest, due to its rich topological band structure [33,34].
However, a major challenge remains in this context, namely, the addition of inter-particle interactions in view of creating novel strongly-correlated states of matter, e.g. fractional Chern insulators [34][35][36][37][38]. Indeed, the interplay between inter-particle interactions and an external drive, such as lattice shaking, has been shown to lead to significant heating and losses in experiments involving Bose-Einstein condensates (BEC) [31,39]. The complexity of this issue lies on the fact that several underlying mechanisms are believed to be responsible for those undesired effects, and these possibly interplay in a complex manner, as we now explain.
On the one hand, at very short times, time-modulated BECs are believed to be mostly affected by strong parametric instabilities, which are characterized by an exponential growth of collective (Bogoliubov) excitations and are accompanied with a fast decay of the BEC; such processes were thoroughly characterized in our previous work [40], where instability rates, stability diagrams and robust physical signatures of such processes were obtained for simple non-resonant shaken systems. Other approaches have been adopted to characterize dynamical instabilities in shaken BECs [41][42][43][44][45], and a recent study even pointed out the possibility of dynamically stabilizing a modulated BEC, inspired by the Kapitza pendulum [46]; the experimental evidence of staggered-states in time-modulated BECs, whose formation also stems from an instability involving the external drive and collective excitations, was recently reported in Ref. [47]; we note that time-modulating the trapping potential can also be exploited to create correlated excitations in BECs [48].
On the other hand, at longer times, dissipative processes are expected to be dominated by scattering events, which are typically associated with two-body processes and are captured by the so-called Floquet-Fermi golden rule [11,49,50]. Finally, the role of (resonant) interband transitions [39,[51][52][53], and the formation of collective emissions of matter-wave jets upon driving [54], have been analyzed in very recent ultracold-atom experiments.
The aim of this paper is to analyze the onset of parametric instabilities in resonantly-modulated BEC's. Building on the tools developed in Ref. [40], we identify the general features of those instabilities that occur in the resonant-modulations context, and provide useful results that could be readily applied to topical systems, such as the driven-induced Harper-Hofstadter model [17][18][19][20][21]31].
The rest of the paper is organized as follows: In Sec. II, we recall the general method of Ref. [40] to treat parametric instabilities. In Sec. III, we apply this method to a simple resonantly-shaken Wannier-Stark ladder, highlighting the main features of these driven-induced instabilities. We then address the topical case of an optical lattice that is modulated by a secondary moving lattice, and which includes a space-dependent phase [18][19][20]; we first study the 1D case in Sec. IV, and then discuss the 2D configuration (leading to the Harper-Hofstadter Hamiltonian [18][19][20]) in Sec. V. Final remarks are provided in the concluding section VI.

II. GENERAL METHOD
We first briefly summarize the general method of Ref. [40] to study parametric instabilities in periodicallydriven BEC lattice systems; we point out that similar or complementary approaches have been proposed to characterize dynamical instabilities in Refs. [41][42][43][44][45]55].

Linear stability analysis
Consider a weakly-interacting Bose gas, in a generic periodically-driven system of period T = 2π/ω; this discussion disregards the system dimension, and/or the presence of a lattice. To study the stability properties of a potential BEC, the general idea is to assume that the system is initially fully-condensed in some well-defined state and to perform a linear stability analysis around this specific state. Let us denote by a n (t = 0) the condensate wavefunction at initial time t = 0, where n is a generic (possibly multidimensional, discrete or continuous) index for spatial position. This state could be set by hand (based on analytical arguments), or it could be more precisely estimated for a given experimental protocol, by numerically implementing the full preparation sequence or using Floquet adiabatic perturbation theory [56,57].
Given this state a (0) n (t = 0) as an input, we are interested in the stability of the full time-dependent solution a (0) n (t). This property can be evaluated by following the guideline below: (i) In the weakly-interacting regime, this solution is governed by the time-dependent Gross-Pitaevskii Equation (tGPE), whose exact form depends on the precise model under consideration. Thus, we first determine the time-evolution of the condensate wavefunction a  n (t = 0). This can be performed numerically using direct real-time propagation in real space. We stress that this calculation is based on the full time-dependent equations, so that the dynamics of the BEC is exactly computed (including all micromotion effects [8]).
(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 (tGPE) in δa n . This yields the timedependent Bogoliubov-de Gennes equations, which take the general form where L(t) is a T -periodic operator (we set = 1 here and in the following). Based on this time-periodicity, it is convenient to exploit the Floquet theorem and to focus the stability analysis on the "time-evolution" (propagator) matrix Φ(T ), which is obtained by time-evolving Eq. (2) over a single period T . From the knowledge of Φ(T ), we extract the "Lyapunov" exponents q , which are related to the eigenvalues λ q of Φ(T ) through the relation λ q = e −i q T ; here we explicitly introduced the momentum q, which will be used to index the corresponding excitation modes. The appearance of Lyapunov exponents with positive imaginary parts indicates a dynamical instability [43,44], i.e. an exponential growth of the corresponding modes, given by the rate s q = Im q .

Observables
As previously discussed in Ref. [40], a first 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 instability rate is independent of the reference frame (or gauge). It quantifies the parametric instabilities occurring in the system in the sense that it governs the stroboscopic dynamics (t = T × integer) of physical observables in the system. For instance, Γ > 0 indicates that physical observables (e.g. the total energy, the depleted fraction, ...) will exponentially grow up with the rate 2Γ [58].
Another relevant indicator is the most unstable mode The excitations momentum distribution indeed shows a pronounced contribution around q mum and structures of momentum q mum develop in real/momentum space, producing clear signatures of those parametric instabilities. We stress that, due to the factorization of the BEC wavefunction in Eq. (1), all momentum modes q (and hence, the most unstable mode) are always defined relatively to the ground-state (which may not always be associated with a vanishing momentum).
As pointed out in Ref. [40], saturation effects generically alter the instability rates (Γ) predicted by Bogoliubov theory at longer times; these effects are mainly attributed to couplings between the Bogoliubov modes. In this sense, the parametric instabilities investigated in this framework are truly associated with short-time dynamics. While the instability rate must therefore be treated as a dynamical quantity (affected by saturation effects and other, possibly incoherent, mechanisms [49,50]), the most unstable mode q mum and the associated structures developing in real and momentum space (e.g. in the momentum distribution) are found to be very robust; this could, therefore, provide clear signatures of parametric instabilities in realistic experimental configurations.

III. A FIRST EXAMPLE: THE RESONANTLY-SHAKEN WANNIER-STARK LADDER
We first apply this method to a simple one-dimensional (1D) resonantly-driven Wannier-Stark ladder; we note that this toy-model is a direct extension of the shaken 1D lattice studied in Ref. [40]. Beyond its physical interest, this Section aims at demonstrating how the (numerical and analytical) tools that were previously developed can indeed be successfully applied in the context of resonantly-driven models.

A. The model
We consider a system of weakly-interacting bosons, trapped in a shaken 1D Wannier-Stark ladder, as described by the periodically-driven Bose-Hubbard Hamiltonian [10] whereâ n annihilates a particle at lattice site n, J > 0 denotes the tunneling amplitude of nearest-neighbor hopping, and U > 0 is the repulsive on-site interaction strength. The second line in Eq. (5) captures the on-site potential term, which contains two effects: the Wannier-Stark-ladder potential, which introduces an energy shift ∆ > 0 between consecutive sites, and a time-periodic modulation of amplitude K and frequency ω = 2π/T ; this time-modulation simply corresponds to an external shaking of the 1D optical lattice, as viewed from the moving frame [10,27]. The frequency modulation ω is chosen to be resonant with the offset ∆, i.e. ∆ = lω with l denoting some integer; see Refs. [29,30]. It is well known that the main effect of the timemodulation is to restore the tunneling, which is suppressed by the strong offset ∆; this is reflected in the (non-interacting) Floquet effective Hamiltonian associated with Eq. (5), and which reads [10,26,29,59] where the effective tunneling is J eff = JJ −l (K/ω), with J l the l-th order Bessel function, and where we set U = 0.

B. Numerical results
To compute the instability rates of the system, we first numerically implement the procedure detailed in Sec. II. The initial state a (0) n (t = 0) is chosen to be the ground state of the naive "interacting effective Hamiltonian" which is a reasonable approximate for an experimentallyprepared ground-state; in general, we note that this ground-state can be numerically determined using imaginary time propagation. For the specific model under consideration, we find that a n (t = 0) is the Bloch state e ip0n of momentum p 0 = 0 if J eff > 0 (homogeneous condensate), and p 0 = π for J eff < 0: see analytical details in Appendix A. This fact is reminiscent of the non-resonant case; see Refs. [40,44]. Figure 1 displays the behavior of the instability rate Γ as a function of the interaction strength g = U ρ (with ρ the condensate density, which enters the normalization of the initial state) and modulation amplitude K/ω, in the purely resonant case (l = 1) and for a driving frequency ω = 5J. We note that the general features are very similar to what was observed in Refs. [40,44] in the non-resonant case : the stability diagram displays lobes, which are separated by stable regions corresponding to cancellation points of the effective tunneling J eff (here the zeros of the function J −1 (K/ω)). For large enough driving frequencies (such as in Fig. 1), the system is stable at low values of g, and a transition to instability appears at finite g; close to the transition, these instabilities are found to be dominated by the Bogoliubov mode q mum = π, which is the most unstable one; for smaller values of ω (smaller than the effective free-particle bandwidth, i.e. ω < 4|J eff | here), the system would be unstable for any non-zero interaction strength, with the most unstable mode corresponding to q mum < π [40].
We find that this situation is very general: For other values of l, we find similar stability diagrams, except that the positions of the lobes are now governed by the Bessel function J −l (K/ω) (instead of J −1 ); besides, similar conclusions hold for the nature of the most unstable mode.

C. Analytical approach
The numerical results described in the previous Section can be understood using the analytical method developed in Ref. [40], which can indeed be readily transposed to the present model; see Appendix A for the full calculations.
The main idea is that the Bogoliubov equations of motion (2) can be mapped onto a parametric oscillator model [45,60], a seminal model of periodically-driven harmonic oscillator known to display dynamical instabilities as soon as the drive frequency approaches twice its own (intrinsic) frequency. Similarly here, we find (see Appendix A) that each Bogoliubov mode q will display a dynamical instability (characterized by an exponential growth of its population) whenever the drive frequency ω approaches its time-averaged Bogoliubov energy, namely, Note that the latter represents the Bogoliubov dispersion associated with the linearized GPE, based on the naive "interacting effective Hamiltonian" in Eq. (7). As detailed in Appendix A, the growth rate s q associated with this instability can be computed analytically using a perturbative method [40,60]. From the knowledge of those individual rates, it is then straightforward to infer the total instability rate Γ and the most unstable mode [see Eqs. (3) and (4)]. Figure 1 (bottom panel) shows the analytical stability diagram associated with the present model, as obtained from the aforementioned analytical perturbative method; as a technical remark, we point out that the latter calculation was performed up to second order with respect to the perturbation's amplitude α q defined in Eq. A11; see Appendix A and Refs. [40,60] for details. It shows very good agreement with the numerical diagram of Fig. 1, and we attribute the small discrepancies to the perturbative nature of the method (higher order terms are generically expected to provide small corrections).
Importantly, the analytical approach explains the existence of stable regions in the vicinity of the cancellation points J eff ≈ 0, as identified in Ref. [44]. Indeed, when J eff vanishes, the time-averaged Bogoliubov dispersion E av (q) becomes trivially flat, so that no excitation mode fulfills the resonance criterion E av (q) ≈ ω associated with the existence of parametric instability.
Besides, the analytical approach also offers a simple view on the boundaries separating stable and unstable regions; in particular, it predicts the nature of the most unstable mode responsible for the onset of instabilities in the vicinity of these boundaries. In the regime where ω > 4|J eff |, no resonance can occur at g = 0, so that the system is stable; upon increasing g, a first mode can become unstable, which is the one of maximum E av (q), namely q = π; see Eq. (8). Conversely, for ω < 4|J eff |, there always exists a particular mode fulfilling the resonance condition : the system is unstable at any inter- action strength, and the most unstable mode is located at a certain momentum q < π (at lowest order, q mum is the mode at resonance; see Eq. (8)). We stress that such conclusions are similar to those found in the non-resonant case [40].

D. Extension to higher dimensional model
The present analysis can be straightforwardly extended to models featuring transverse directions (be it lattice or continuous directions).
For instance, in the case where of a continuous transverse degree of freedom is present, we find the stability diagram of Fig. 2. As observed in Ref. [40], the presence of transverse modes enhance the instabilities by opening new instability channels. More precisely, in the case of an unbounded bandwidth, there always exists Bogoliubov modes that are resonant with the drive frequency ω; as a consequence, instabilities occur at any finite interaction strength g. We note that, in contrast with the purely 1D case treated above, the cancellation of the effective tunneling does not result in a vanishing of E av ; hence, in this case, the stability regions located near the zeros of J 1 (K/ω) are rather imputable to the cancellation of the perturbation's amplitude α q that is associated with the underlying (effective) parametric oscillator [Eq. A11 in Appendix A], and which also scales as J 1 (K/ω).
Finally, in that case, simple analytical formulas are obtained for both Γ and q mum [see Appendix for details; here the index x denotes the lattice direction]: (ii) If ω < 4|J eff |(4|J eff | + 2g), we find We note a strong similarity with the non-resonant case [40]; a notable difference, though, is the fact that the rate Γ does not depend on the shaking amplitude K in the low-frequency regime [Eq. (12)]. These predictions, in particular the existence of two different regimes characterized by very specific dependences of Γ and q mum on the model parameters, could be readily tested by present-day experiments, providing a clear and unambiguous signature of parametric instabilities in resonantly-modulated systems.

IV. MOVING LATTICES WITH A SPACE-DEPENDENT PHASE: THE 1D CASE
We now consider another type of modulation scheme, which involves a main (primary) optical lattice that is perturbed by a (secondary) moving-lattice potential [26,61]. Introducing a space-dependent phase in this moving-lattice potential has been shown to generate non-trivial effective gauge fields and topological band structures in the context of 2D neutral gases [17][18][19][20][21]31]. Before addressing the case of 2D systems (Section V), where such gauge structures appear, we first investigate the properties of a simpler 1D toy model.

The model
We consider a 1D model described by the Hamiltonian whereâ n annihilates a particle at lattice site n, J > 0 denotes the tunneling amplitude of nearest-neighbor hopping, U > 0 is the repulsive on-site interaction strength, and ∆ > 0 is the energy difference between two consecutive sites. The second line of Eq. (13) captures the effects of the secondary (moving-lattice) potential; the latter is characterized by the amplitude K, the frequency ω, and a phase difference of θ = π between two consecutive sites.
In the following, we consider the resonant case where ω = ∆. Before analyzing the existence of parametric instabilities in this model, we point out that such moving-lattice potentials generically produce momentum kicks, which can be directly revealed in momentum distributions [26]. However, these effects have a zero average over one period of the drive [7]; in particular, such momentum kicks do not influence the parametric instabilities explored in the present work.
It is convenient to perform our analysis in the rotating frame, defined by the unitary transformation with α = K/ω. In this frame, the Hamiltonian readŝ so that translational invariance (with a periodicity of two lattice sites) is restored.
In the absence of interactions (U = 0), we recall that the Floquet effective Hamiltonian associated with Eq. (14) simply reads [26] where J eff = JJ −1 (2α). The Hamiltonian in Eq. (15) has a periodicity of two lattice sites, so that its spectrum [shown in Fig. 3(a)] displays two energy bands, We note that the corresponding eigenstates are labelled by the quasi-momentum p, which is defined in a reduced Brillouin zone p ∈ [−π/2; π/2]. Here, the ground state corresponds to the state with p 0 = π/2 in the " − " band; see Fig. 3(a). Thus, this state has alternate on-site amplitudes from one unit cell to the consecutive one; furthermore, depending on the sign of J eff , we find that its one-site coefficients are either equal or opposite within each cell; therefore, the ground state is found to be the same, globally, whatever the sign of J eff [see Fig. 3(c)]; this is in striking contrast with the models analyzed in the previous Section III and in Ref. [40].

Numerical results
As previously in Section III, the initial state is again taken to be the ground state of the naive "interacting effective Hamiltonian" [i.e. Eq. (7) with Eq. (15)], which is obtained numerically through imaginary-time propagation. Interestingly, we find that this state is still characterized by the quasi-momentum p 0 = π/2, even in the presence of interactions.
The stability diagram obtained from this initial state is displayed in Fig. 4, which shows the instability rate Γ as a Time-averaged Bogoliubov dispersion for the same parameters and g = 3J. The two branches are reminiscent of those of the single-particle spectrum (the global shift of π/2 comes from the fact that excitation momenta are defined with respect to the ground-state, which is a π/2 state), but these are modified by the interactions; in particular, Goldstone modes arise at q ≈ 0 in the lowest band, as generically expected. (c) Qualitative picture of the ground-state in which condensation occur; this state does not depend on the sign of J eff . function of the modulation amplitude and the interaction strength.
The stability diagram in Fig. 4 is dominated by a quasi "periodic" structure, as a function of α = K/ω, which arises from the dependence of the effective tunneling along the lattice direction, JJ −1 (2α): similarly to what was observed for other modulated band models (see Refs. [40,44] and Section III C), we note that stable re-gions are indeed privileged when the effective tunneling J eff ≈ 0, which is consistent with the parametric resonance criterion E av (q) = ω introduced in Section III C.
This important observation allows one to anticipate the presence of stable regions in other time-modulation schemes. An immediate extension is the case where the phase difference θ between consecutive sites [which we took equal to θ = π in Eq. (13)] takes another value. Effective Hamiltonians for such models can be simply obtained using the formulas of Ref. [7]; for instance, considering a moving lattice with a phase difference of θ = π/2, we find an effective tunneling J eff = JJ −1 (4α) [see [7] for details], which implies a stability diagram associated with narrower instability lobes. More generally, for mod-els of the form given in Eq. (13) with arbitrary phase θn, stable zones are obtained whenever JJ −l (pα) ≈ 0, where l is the order of the resonance (∆ = lω) and where p is the spatial periodicity of the phase nθ entering the time-modulation (i.e. the moving lattice).
Besides, our numerical calculations reveal that the onset of instability (i.e. the boundaries of the stability diagram in Fig. 4) is dominated by a most unstable mode, which in this case corresponds to the q = 0 mode; we remind that this momentum value is evaluated with respect to the ground-state, as always tacitly implied within our formalism. To understand this, let us compute the timeaveraged Bogoliubov dispersion E av (q) for the present model, which yields E av (q) = 4|J eff | 2 (1 + cos 2 q) + 4|J eff |g ± 16|J eff | 2 g 2 cos 4 q + 64|J eff | 3 cos 2 q(|J eff | + g).
As shown in Fig. 3(b), this dispersion is made of two branches, which are reminiscent of the two-branch single particle spectrum [62], but are modified by the interactions. In particular, Goldstone modes arise at q ≈ 0 in the lowest band, as generically expected as a result from the broken U (1) symmetry [63]. The absolute maximum is reached for q = 0 in the upper branch, which indeed precisely corresponds to the most unstable mode identified in our numerics: similarly to the previous model discussed in Section III, we thus find that the onset of parametric instability is governed by the mode of highest time-averaged Bogoliubov energy, which is consistent with the fact that this mode will be the first one to resonate [E av (q = 0) = ω] as one increases the interaction strength g [40]. We emphasize that this most unstable mode differs from the one identified for the resonantly-shaken Wannier-Stark ladder of Section III.
It is remarkable that the conclusions emanating from our analysis of the moving lattice are qualitatively similar to those related to other shaken-lattice models; see Sec. III and Ref. [40]. This suggests a reliable and intuitive guideline to predict stable zones of the stability diagram and to identify the most unstable mode, based on the simple knowledge of the effective band structure (i.e. the effective tunneling and the time-averaged Bogoliubov dispersion).

V. THE DRIVEN-INDUCED HARPER-HOFSTADTER MODEL
In this Section, we extend our previous analysis [Section IV] to a two-dimensional setting, which has been exploited to realize the Harper-Hofstadter Hamiltonian in cold atoms [17][18][19][20][21]31].

The model
We consider a two-dimensional extension of the previous model [Eq. (13)], which we define by the Hamiltonian whereâ m,n annihilates a particle at lattice site (m, n), J x (resp. J y ) is the tunneling along the x (resp. y) direction, and U models on-site repulsive interactions. A linear (Wannier-Stark) potential is introduced along the x-direction; besides the model features a moving lattice of amplitude K, frequency ω and a π-phase dependence along both directions (θ x,y = π). In the following, we set the resonance condition ω = ∆. This model, and other variants, have been realized in several recent ultracold-atom experiments [17][18][19][20][21]31]. Here, we focus on the "π-flux" configuration [19,20], which corresponds to setting θ y = π in Eq. (18), but we note that several of our results hold for other values of the synthetic flux (see below).
In the absence of interactions (U = 0), the effective Hamiltonian associated with Eq. (18) reads [7] with α = K/ω, J x eff = J x J −1 (2α), and J y eff = J y J 0 (2α). This corresponds to the Harper-Hofstadter model [32], with a magnetic flux Φ = π in each unit cell, and with  5. (a) Single-particle spectrum of the effective Hamiltonian in Eq. (19) for Jx,y = 1 and α = 0.5; each branch is twofold degenerate, and the two ground-states correspond to states at px = π/2, py = 0 in the " − " band; (b) Cut of the single-particle spectrum at py = 0; (c) Time-averaged Bogoliubov dispersion for the same parameters and g = 3J, plotted as a function of qx at fixed qy = 0. The four branches are reminiscent of the two branches of the single-particle spectrum, with a lifting of their twofold degeneracy through interactions; the global shift of π/2 in the x-direction comes from the fact that excitation momenta are defined with respect to the ground-state, which is a px = π/2 state; in particular, Goldstone modes arise at q ≈ 0 in the lowest band, as generically expected. different effective tunnelings along the x and y directions. As stated above, other choices for the movinglattice phase (θ y = π) lead to different fluxes per unit cell [7,18,31]. We recall that the Harper-Hofstadter model displays the well-known "Hofstadter's butterfly" spectrum [32], a rich fractal structure that hosts Dirac semimetals and Chern insulating phases [33]. In this sense, at the single-particle level, the time-dependent model under consideration [i.e. Eq. (18) with U = 0] is one of the simplest systems realizing artificial gauge fields and topological band structures for neutral atoms in 2D optical lattices [17][18][19][20][21]31]; we point out that a BEC was recently realized in the "π-flux" configuration of this model, which further motivates the present study [20].

Numerical results
In order to extract the instability properties of the 2D model in Eq. (18), we now apply the same procedure as for the 1D model of Section IV. The initial state is again chosen as the ground state of the naive "interacting effective Hamiltonian"; however, we point out that a complexity arises here from the degeneracy of this condensation state [32,64].
To see this, one observes that the single-particle effective Hamiltonian in Eq. (19) is periodic over a 2 × 2 cell, so that its eigenstates are labelled by the quasimomentum p defined in the reduced Brillouin zone p x , p y ∈ [−π/2; π/2]. The spectrum features two energy bands (labelled ±) which are both twofold degenerate [see Figs. 5(a-b)]. The ground state corresponds to the state with p = (π/2, 0) in the " − " band, and is twofold degenerate [65].
In the presence of interactions, the ground state still features this twofold degeneracy [64], as well as the momentum characteristic (p x = π/2, p y = 0). Therefore, several choices can a priori be made for the initial state of our analysis, within the whole degenerate ground space, as we now explore.
Stability diagram and sensitivity to the ground state Figure 6 shows the instability diagrams obtained for two different orthogonal ground-states in this subspace. At first sight, they are found to be very similar in the sense that the rates are of the same order of magnitude, and the diagrams feature very similar structures as a function of the model parameters; in particular, the same "periodicity" as in Fig. 4 is observed as a function of modulation amplitude α = K/ω. Yet, the exact positions of the stable and unstable regions are shifted, reflecting that the later are very sensitive to the ground-state in which the system is prepared. The later observation that stable regions cannot be simply obtained by analyzing the behavior of the effective tunneling J x,y eff [Section III C] can be related to the fact that these quantities never cancel simultaneously; indeed, in this specific 2D model, the time-averaged Bogoliubov dispersion never vanishes, meaning that the simple "flat-band" criterion of Section III C cannot be applied (a priori, there might always exist some modes fulfilling the resonance condition E av = ω). Consequently, we conjecture that the stable zones in Fig. 4 should rather be governed by the amplitude of the perturbation entering the underlying effective parametric-oscillator model [i.e. the analogue of the quantity α q that was defined in Eq. A11, in Appendix A, for the model of Sec. III]; although its full analytical derivation cannot be obtained in the present case, it is natural to believe that this amplitude is indeed ground-state dependent, since the Bogoliubov treatment leading to this effective parametric-oscillator model relies on the actual condensation state that is chosen within the degenerate ground manifold.

Most unstable mode
In contrast, we find that the most unstable mode associated with the onset of the instability is robust with respect to the ground-state. Indeed, independently of the ground state, and within the whole parameter range of the diagram, we find that this most unstable mode always corresponds to the q = 0 mode (with respect to the ground-state). This can again be accounted for by the fact that it is the mode of highest time-averaged Bogoliubov energy. To verify this, we calculate the timeaveraged Bogoliubov spectrum E av , which is obtained by numerically diagonalizing the Bogoliubov Hamiltonian derived from the GPE combined with the static effective Hamiltonian in Eq. (19); we note that analytics exist in the symmetric case where J x eff = J y eff ; see [64]. The resulting Bogoliubov spectrum E av is plotted in Fig. 5(c). Consistently with the general features reported in [64], this spectrum is made of four branches (the degeneracy observed at the single-particle level being lifted) and indeed has its absolute maximum at q = 0 in the upper branch.
Interestingly, this behavior is in fact expected to be very generic. For a model configuration leading to a flux Φ = 1/3, i.e. θ y = π/3, we also find a most unstable mode at q = 0, consistently with the fact that it is the mode of highest time-averaged Bogoliubov energy [see Ref. [64] for the Bogoliubov spectrum in that case].
We conclude this Section by noting that the stability diagram of the driven-induced Harper-Hofstadter model [ Fig. 6] displays relatively large stability regions, which extend up to significant values of the interaction strength g, hence suggesting potentially favorable regimes of operation (as far as parametric instabilities are concerned). However, one should point out that current experiments typically feature a transverse ("tube") direction, which is generically expected to trigger or increase instabilities [40]. Besides, our result that the stability diagram strongly depends on the prepared ground state appears as an important feature of these time-modulated systems, which should be taken in account in experiments.

VI. CONCLUSION
In this work, we analyzed the parametric instabilities that occur in a variety of resonantly-driven band models. While the general method used to identify and characterize these instabilities had already been applied to nonresonant models [40,44], we have hereby demonstrated its usefulness and flexibility to address the relevant case of resonant time-modulations. Ab initio predictions for instability rates and for the most unstable mode (which appears to be the most robust and directly accessible signature of parametric instabilities) were obtained; these could be directly tested in current ultracold-atom experiments [17][18][19][20][21]31], for instance, in view of optimizing their operating regimes.
Our study has confirmed the generic role played by parametric resonances in the instabilities that appear in these driven system, and which directly involve the drive frequency and the dispersion of the Floquet-Bogoliubov modes. In particular, while genuine analytical results can only be obtained in simple cases, many qualitative features of these instabilities appear to be generic, hence providing an intuitive picture that could guide experiments to stable regimes of operation.
On the one hand, instability rates exhibit a dependence on the modulation amplitude that is, in most cases, governed by the effective tunneling, and favorable regions are generally found whenever this effective tunneling is weak (based on the resonance criterion E av = ω underlying parametric instabilities and the fact that E av ≈ 0 in these regions); while this conclusion is justified in models where all effective tunneling amplitudes can simultaneously vanish (as in the simple 1D models discussed in this work), it should however be treated with care in higher dimensions, where stability regions seem to be rather governed by the amplitude of the perturbation entering the underlying effective parametric-oscillator model [i.e. the analogue of the quantity α q in Eq. A11 in Appendix A]; in particular, in the presence of degenerate ground states, as in the "Harper-Hofstadter" case treated in Section V, the instability rates were found to be very sensitive to the actual ground-state in which the system was prepared. Besides, we note that the energy scales of the (effective) system mostly rely on the effective tunneling, which indicates that a compromise must be found in actual implementations (in the context of fractional Chern insulators, it is crucial to maximize the size of topological gaps in view of generating strongly-correlated states in the presence of finite temperature).
On the other hand, for models exhibiting a finite bandwidth, the most unstable mode responsible for the onset of instability is always found to be that of highest effective Bogoliubov energy E av (for experimental values of the drive frequency [40]). When a continuum is present and the dispersion is unbounded, the most unstable mode is found (at lowest order) to fulfill the resonance condition E av (q mum ) = ω. Therefore, the simple knowledge of the effective band structure [J eff and E av ] allows one to anticipate such behaviors, which are expected to be directly relevant to experiments.
Altogether, this work constitutes a first step in view of controlling and exploiting the potentialities of interacting modulated BECs. In particular, achieving stable BECs in the context of time-modulated optical lattices, with reduced instabilities and losses, will constitute a first step towards the cold-atom realization of strongly-correlated states of matter with topological features, such as fractional Chern insulators [34][35][36][37][38]. However, many open problems remain. First of all, this study was performed it the mean-field interaction regime, which is expected to give access to the short-term dynamics. Longer-time dynamics, including the possibility for thermalization, are dominated by non-linear couplings between the excitation modes and the condensate, which go beyond the present Bogoliubov treatment. Moreover, the analysis presented in this work neglects the effects of higher bands, which are indeed expected to be weak when setting the drive frequency within a gap of the (effective) spectrum; however, higher-band effects are also expected to become important at longer times. Then, independently from the driving scheme itself and the obtained stability properties, a central question is that of the preparation of the initial state, as loading the system into a desired eigenstate with highest fidelity is far from being trivial: one solution could be to apply adiabatic perturbation theory in the presence of the periodic drive [56,57], but there might exist other alternatives and their impact on the stability properties of the prepared state remains uncharacterized. Finally, the interplay between parametric instabilities and other instability mechanisms neglected in our approach, especially inter-band transitions [39,51,52], is expected to lead to rich behaviors that still remain to be studied.

Mapping on a parametric oscillator model
Proceeding as in [40], we then introduce the two following changes of basis: first, where cosh(2θ q ) ≡ (ε av (q) + g)/E av (q) and sinh(2θ q ) ≡ g/E av (q) with E av (q) = ε av (q)(ε av (q) + 2g) (A6) the time-averaged Bogoliubov dispersion and ε av (q) = 4|J eff | sin 2 (q x /2) + ⊥ q ⊥ ; and then, This allows us to write the Bogoliubov equations under the form where E av (q) is the time-averaged Bogoliubov dispersion associated with the effective GPE, within the Bogoliubov approximation [see Eq. (A6)],Ŵ q (t) is a diagonal matrix of zero average over one driving period, As shown in [40], Eqs. A8 are formally equivalent to a set of independent parametric oscillators (one for each mode q), since they describe harmonic oscillators of eigenfrequencies E av (q) driven by a weak time-periodic modulation of frequency ω. As can be anticipated from a RWA treatment of Eq. A8, a parametric instability will thus appear in the system as soon as one of the harmonics of the modulation is close to twice the energy of any of the (effective, timeaveraged) Bogoliubov modes, 2E av (q), i.e. mω ≈ 2E av (q) (resonance condition). The later will result in a long-term explosion in the (stroboscopic) dynamics of the corresponding mode. Assuming that the harmonics can be treated independently, the instability rate of a given mode q and a given harmonics can be analytically computed following the perturbative method developed in [40,60], and the total instability rate and most unstable mode straightforwardly inferred from Eqs. (3) and (4). 3. Explicit analytics in the strict resonant case ω = ∆ As an illustration, we now consider the simple case l = 1 (i.e. ω = ∆) and φ = 0 addressed in the main text. For the typical parameters used in Fig. 1, one immediately sees that only the smaller harmonics of the modulation will allow the resonance condition to be fulfilled; this corresponds to the terms m = ±1 in Eq. (A9). We can thus disregard other harmonics, as well as the termŴ q (t) which has no long-term stroboscopic contribution [66], and we get the effective model i∂ t ũ q v q = E av (q)1 + α q E av (q) 2 0 cos(2ωt)e −2iEav(q)t − cos(2ωt)e 2iEav(q)t 0 with α q = 8JJ −1 (α) sin 2 (q x /2) g [E av (q)] 2 (A11) As shown in Refs. [40,45], Eq. A10 is indeed formally equivalent to the equations of motion describing a parametric oscillator of eigenfrequency E av (q), driven by a perturbation of amplitude α q and frequency 2ω. The quantity α q will thus be referred to as the perturbation's amplitude in the effective parametric oscillator model; it is the small parameter of the perturbative method used to compute the instability rates [40,60]. At lowest order in α q , we find in this case In the purely 1D case (no transverse direction), this yields the analytical stability diagram on the bottom panel of Fig. 1 (although the latter is actually computed using the analytical expression at second order, which is an implicit equation). In that case, the most unstable mode is q mum x = π at the onset of instability. In the case where a continuous degree of freedom is present, we find two regimes : (i) If ω > 4|J eff |(4|J eff | + 2g), one finds q mum x = π; (q mum ⊥ ) 2 /2m = g 2 + ω 2 − g − 4|J eff |. (A13)