Parametric Instabilities in ResonantlyDriven BoseEinstein Condensates
Abstract
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 HarperHofstadter model, which combined optical superlattices and resonant timemodulations. Adding interparticle interactions to these engineered band systems is expected to lead to stronglycorrelated states with topological features, such as fractional Chern insulators. However, the interplay between interactions and external timeperiodic 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 earlystage parametric instabilities that occur in systems of resonantlydriven BoseEinstein condensates in optical lattices. We apply and extend an approach based on Bogoliubov theory [PRX 7, 021015 (2017)] to a variety of resonantlydriven band models, from a simple shaken WannierStark ladder to the more intriguing driveninduced HarperHofstadter 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 Oka and Aoki (2009); Kitagawa et al. (2010); Lindner et al. (2011); Cayssol et al. (2013); Mitrano et al. (2016); Aidelsburger et al. (2017). In the context of ultracold quantum gases, it was shown that subjecting neutral atoms to an external timeperiodic drive could be used to design artificial gauge fields in these systems Goldman et al. (2014); Goldman and Dalibard (2014); Bukov et al. (2015a); Eckardt (2016); Aidelsburger et al. (2017), hence opening promising perspectives in the quantum simulation of topological states of matter Aidelsburger et al. (2015a); Jotzu et al. (2014); Goldman et al. (2016) and quantum magnetism Eckardt et al. (2010); Struck et al. (2011); Eckardt (2016); see also the recent work Görg et al. (2017). The underlying concept of Floquet engineering Cayssol et al. (2013); Goldman and Dalibard (2014); Bukov et al. (2015a); Goldman et al. (2016); Eckardt (2016) 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 Aidelsburger et al. (2011, 2013); Miyake et al. (2013); Aidelsburger et al. (2015a); Jotzu et al. (2014); Kennedy et al. (2015); Tai et al. (2017); Tarnowski et al. (2017), some of which lead to nontrivial geometrical and topological properties Goldman et al. (2016).
A particularly interesting class of periodicallydriven setups is that featuring resonant timemodulations Kolovsky (2011); Bermudez et al. (2011); Hauke et al. (2012); Goldman et al. (2015); Creffield et al. (2016); Price et al. (2017), 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 WannierStark ladder; see Refs. Sias et al. (2008); Mukherjee et al. (2015) for experimental realizations and Ref. Eckardt (2016) for a review. This socalled photonassisted tunneling effect constitutes a natural ingredient for the generation of artificial fluxes within coldatom systems Kolovsky (2011); Hauke et al. (2012); Goldman et al. (2015); Creffield et al. (2016), as exemplified by the recent experimental realizations of the HarperHofstadter model Aidelsburger et al. (2011, 2013); Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2015b); Tai et al. (2017), a 2D lattice penetrated by a uniform magnetic flux Hofstadter (1976); the latter band model is of particular interest, due to its rich topological band structure Kohmoto (1989); Möller and Cooper (2015).
However, a major challenge remains in this context, namely, the addition of interparticle interactions in view of creating novel stronglycorrelated states of matter, e.g. fractional Chern insulators Regnault and Bernevig (2011); Sterdyniak et al. (2013); Bergholtz and Liu (2013); Grushin et al. (2014); Möller and Cooper (2015). Indeed, the interplay between interparticle interactions and an external drive, such as lattice shaking, has been shown to lead to significant heating and losses in experiments involving BoseEinstein condensates (BEC) Aidelsburger et al. (2015b); Reitter et al. (2017). 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, timemodulated 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 Lellouch et al. (2017), where instability rates, stability diagrams and robust physical signatures of such processes were obtained for simple nonresonant shaken systems. Other approaches have been adopted to characterize dynamical instabilities in shaken BECs Modugno et al. (2004); Kramer et al. (2005); Tozzo et al. (2005); Creffield (2009); Bukov et al. (2015b), and a recent study even pointed out the possibility of dynamically stabilizing a modulated BEC, inspired by the Kapitza pendulum Martin et al. (2017); the experimental evidence of staggeredstates in timemodulated BECs, whose formation also stems from an instability involving the external drive and collective excitations, was recently reported in Ref. Michon et al. (2017); we note that timemodulating the trapping potential can also be exploited to create correlated excitations in BECs Jaskula et al. (2012).
On the other hand, at longer times, dissipative processes are expected to be dominated by scattering events, which are typically associated with twobody processes and are captured by the socalled FloquetFermi golden rule Bilitewski and Cooper (2015a, b); Aidelsburger et al. (2015a). Finally, the role of (resonant) interband transitions Weinberg et al. (2015); StrÃ¤ter and Eckardt (2016); Reitter et al. (2017); Quelle and Morais Smith (2017), and the formation of collective emissions of matterwave jets upon driving Logan W. Clark (2017), have been analyzed in very recent ultracoldatom experiments.
The aim of this paper is to analyze the onset of parametric instabilities in resonantlymodulated BEC’s. Building on the tools developed in Ref. Lellouch et al. (2017), we identify the general features of those instabilities that occur in the resonantmodulations context, and provide useful results that could be readily applied to topical systems, such as the driveninduced HarperHofstadter model Aidelsburger et al. (2011, 2013); Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2015b); Tai et al. (2017).
The rest of the paper is organized as follows: In Sec. II, we recall the general method of Ref. Lellouch et al. (2017) to treat parametric instabilities. In Sec. III, we apply this method to a simple resonantlyshaken WannierStark ladder, highlighting the main features of these driveninduced instabilities. We then address the topical case of an optical lattice that is modulated by a secondary moving lattice, and which includes a spacedependent phase Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2013); we first study the 1D case in Sec. IV, and then discuss the 2D configuration (leading to the HarperHofstadter Hamiltonian Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2013)) in Sec. V. Final remarks are provided in the concluding section VI.
Ii General method
We first briefly summarize the general method of Ref. Lellouch et al. (2017) 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. Modugno et al. (2004); Kramer et al. (2005); Tozzo et al. (2005); Creffield (2009); Bukov et al. (2015b); Salerno et al. (2016).
Linear stability analysis
Consider a weaklyinteracting Bose gas, in a generic periodicallydriven system of period ; 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 fullycondensed in some welldefined state and to perform a linear stability analysis around this specific state. Let us denote by the condensate wavefunction at initial time , where 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 Weinberg et al. (2016); Novicenko et al. (2016).
Given this state as an input, we are interested in the stability of the full timedependent solution . This property can be evaluated by following the guideline below:
(i) In the weaklyinteracting regime, this solution is governed by the timedependent GrossPitaevskii Equation (tGPE), whose exact form depends on the precise model under consideration. Thus, we first determine the timeevolution of the condensate wavefunction , by solving the (tGPE) with the initial condition . This can be performed numerically using direct realtime propagation in real space. We stress that this calculation is based on the full timedependent equations, so that the dynamics of the BEC is exactly computed (including all micromotion effects Goldman and Dalibard (2014)).
(ii) Given the timedependent solution for the condensate wave function , we analyse its stability by considering a small perturbation
(1) 
and linearizing the (tGPE) in . This yields the timedependent Bogoliubovde Gennes equations, which take the general form
(2) 
where is a periodic operator (we set here and in the following). Based on this timeperiodicity, it is convenient to exploit the Floquet theorem and to focus the stability analysis on the “timeevolution” (propagator) matrix , which is obtained by timeevolving Eq. (2) over a single period . From the knowledge of , we extract the “Lyapunov” exponents , which are related to the eigenvalues of through the relation ; here we explicitly introduced the momentum , which will be used to index the corresponding excitation modes. The appearance of Lyapunov exponents with positive imaginary parts indicates a dynamical instability Tozzo et al. (2005); Creffield (2009), i.e. an exponential growth of the corresponding modes, given by the rate .
Observables
As previously discussed in Ref. Lellouch et al. (2017), a first quantitative indicator of the instability is the maximum growth rate of the spectrum,
(3) 
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 () of physical observables in the system. For instance, indicates that physical observables (e.g. the total energy, the depleted fraction, …) will exponentially grow up with the rate ^{1}^{1}1The factor of stems from the fact that is defined in amplitude while physical observables involve squared moduli of wavefunctions..
Another relevant indicator is the most unstable mode
(4) 
The excitations momentum distribution indeed shows a pronounced contribution around and structures of momentum 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 (and hence, the most unstable mode) are always defined relatively to the groundstate (which may not always be associated with a vanishing momentum).
As pointed out in Ref. Lellouch et al. (2017), 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 shorttime dynamics. While the instability rate must therefore be treated as a dynamical quantity (affected by saturation effects and other, possibly incoherent, mechanisms Bilitewski and Cooper (2015a, b)), the most unstable mode 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 resonantlyshaken WannierStark ladder
We first apply this method to a simple onedimensional (1D) resonantlydriven WannierStark ladder; we note that this toymodel is a direct extension of the shaken 1D lattice studied in Ref. Lellouch et al. (2017). 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 resonantlydriven models.
iii.1 The model
We consider a system of weaklyinteracting bosons, trapped in a shaken 1D WannierStark ladder, as described by the periodicallydriven BoseHubbard Hamiltonian Eckardt (2016)
(5)  
where annihilates a particle at lattice site , denotes the tunneling amplitude of nearestneighbor hopping, and is the repulsive onsite interaction strength. The second line in Eq. (5) captures the onsite potential term, which contains two effects: the WannierStarkladder potential, which introduces an energy shift between consecutive sites, and a timeperiodic modulation of amplitude and frequency ; this timemodulation simply corresponds to an external shaking of the 1D optical lattice, as viewed from the moving frame Eckardt (2016); Creffield et al. (2016). The frequency modulation is chosen to be resonant with the offset , i.e. with denoting some integer; see Refs. Sias et al. (2008); Mukherjee et al. (2015).
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 (noninteracting) Floquet effective Hamiltonian associated with Eq. (5), and which reads Sias et al. (2008); Choudhury and Mueller (2015); Goldman et al. (2015); Eckardt (2016)
(6) 
where the effective tunneling is , with the lth order Bessel function, and where we set .
iii.2 Numerical results
To compute the instability rates of the system, we first numerically implement the procedure detailed in Sec. II. The initial state is chosen to be the ground state of the naive “interacting effective Hamiltonian”
(7) 
which is a reasonable approximate for an experimentallyprepared groundstate; in general, we note that this groundstate can be numerically determined using imaginary time propagation. For the specific model under consideration, we find that is the Bloch state of momentum if (homogeneous condensate), and for : see analytical details in Appendix A. This fact is reminiscent of the nonresonant case; see Refs. Creffield (2009); Lellouch et al. (2017).
Figure 1 displays the behavior of the instability rate as a function of the interaction strength (with the condensate density, which enters the normalization of the initial state) and modulation amplitude , in the purely resonant case () and for a driving frequency .
We note that the general features are very similar to what was observed in Refs. Creffield (2009); Lellouch et al. (2017) in the nonresonant case : the stability diagram displays lobes, which are separated by stable regions corresponding to cancellation points of the effective tunneling (here the zeros of the function ). For large enough driving frequencies (such as in Fig. 1), the system is stable at low values of , and a transition to instability appears at finite ; close to the transition, these instabilities are found to be dominated by the Bogoliubov mode , which is the most unstable one; for smaller values of (smaller than the effective freeparticle bandwidth, i.e. here), the system would be unstable for any nonzero interaction strength, with the most unstable mode corresponding to Lellouch et al. (2017).
We find that this situation is very general: For other values of , we find similar stability diagrams, except that the positions of the lobes are now governed by the Bessel function (instead of ); besides, similar conclusions hold for the nature of the most unstable mode.
iii.3 Analytical approach
The numerical results described in the previous Section can be understood using the analytical method developed in Ref. Lellouch et al. (2017), 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 Landau and Lifshits (1969); Bukov et al. (2015b), a seminal model of periodicallydriven 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 will display a dynamical instability (characterized by an exponential growth of its population) whenever the drive frequency approaches its timeaveraged Bogoliubov energy, namely, , with
(8) 
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 associated with this instability can be computed analytically using a perturbative method Landau and Lifshits (1969); Lellouch et al. (2017). 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 defined in Eq. 31; see Appendix A and Refs. Landau and Lifshits (1969); Lellouch et al. (2017) 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 , as identified in Ref. Creffield (2009). Indeed, when vanishes, the timeaveraged Bogoliubov dispersion becomes trivially flat, so that no excitation mode fulfills the resonance criterion 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 , no resonance can occur at , so that the system is stable; upon increasing , a first mode can become unstable, which is the one of maximum , namely ; see Eq. (8). Conversely, for , there always exists a particular mode fulfilling the resonance condition : the system is unstable at any interaction strength, and the most unstable mode is located at a certain momentum (at lowest order, is the mode at resonance; see Eq. (8)). We stress that such conclusions are similar to those found in the nonresonant case Lellouch et al. (2017).
iii.4 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. Lellouch et al. (2017), 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 . 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 ; hence, in this case, the stability regions located near the zeros of are rather imputable to the cancellation of the perturbation’s amplitude that is associated with the underlying (effective) parametric oscillator [Eq. 31 in Appendix A], and which also scales as .
Finally, in that case, simple analytical formulas are obtained for both and [see Appendix for details; here the index denotes the lattice direction]:

If , one finds
(9) (10) 
If , we find
(11) (12)
We note a strong similarity with the nonresonant case Lellouch et al. (2017); a notable difference, though, is the fact that the rate does not depend on the shaking amplitude in the lowfrequency regime [Eq. (12)]. These predictions, in particular the existence of two different regimes characterized by very specific dependences of and on the model parameters, could be readily tested by presentday experiments, providing a clear and unambiguous signature of parametric instabilities in resonantlymodulated systems.
Iv Moving lattices with a spacedependent 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) movinglattice potential StamperKurn and Ketterle (2001); Goldman et al. (2015). Introducing a spacedependent phase in this movinglattice potential has been shown to generate nontrivial effective gauge fields and topological band structures in the context of 2D neutral gases Aidelsburger et al. (2011, 2013); Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2015b); Tai et al. (2017). 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
(13)  
where annihilates a particle at lattice site , denotes the tunneling amplitude of nearestneighbor hopping, is the repulsive onsite interaction strength, and is the energy difference between two consecutive sites. The second line of Eq. (13) captures the effects of the secondary (movinglattice) potential; the latter is characterized by the amplitude , 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 movinglattice potentials generically produce momentum kicks, which can be directly revealed in momentum distributions Goldman et al. (2015). However, these effects have a zero average over one period of the drive Goldman et al. (2014); 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 . In this frame, the Hamiltonian reads
(14)  
so that translational invariance (with a periodicity of two lattice sites) is restored.
In the absence of interactions (), we recall that the Floquet effective Hamiltonian associated with Eq. (14) simply reads Goldman et al. (2015)
(15) 
where . 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,
(16) 
We note that the corresponding eigenstates are labelled by the quasimomentum , which is defined in a reduced Brillouin zone . Here, the ground state corresponds to the state with in the band; see Fig. 3(a). Thus, this state has alternate onsite amplitudes from one unit cell to the consecutive one; furthermore, depending on the sign of , we find that its onesite coefficients are either equal or opposite within each cell; therefore, the ground state is found to be the same, globally, whatever the sign of [see Fig. 3(c)]; this is in striking contrast with the models analyzed in the previous Section III and in Ref. Lellouch et al. (2017).
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 imaginarytime propagation. Interestingly, we find that this state is still characterized by the quasimomentum , 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 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 , which arises from the dependence of the effective tunneling along the lattice direction, : similarly to what was observed for other modulated band models (see Refs. Creffield (2009); Lellouch et al. (2017) and Section III.3), we note that stable regions are indeed privileged when the effective tunneling , which is consistent with the parametric resonance criterion introduced in Section III.3.
This important observation allows one to anticipate the presence of stable regions in other timemodulation 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. Goldman et al. (2014); for instance, considering a moving lattice with a phase difference of , we find an effective tunneling [see Goldman et al. (2014) for details], which implies a stability diagram associated with narrower instability lobes. More generally, for models of the form given in Eq. (13) with arbitrary phase , stable zones are obtained whenever , where is the order of the resonance () and where is the spatial periodicity of the phase entering the timemodulation (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 mode; we remind that this momentum value is evaluated with respect to the groundstate, as always tacitly implied within our formalism. To understand this, let us compute the timeaveraged Bogoliubov dispersion for the present model, which yields
(17) 
As shown in Fig. 3(b), this dispersion is made of two branches, which are reminiscent of the twobranch single particle spectrum ^{2}^{2}2The global shift of comes from the fact that excitation momenta are defined with respect to the groundstate, which is a state in the present case., but are modified by the interactions. In particular, Goldstone modes arise at in the lowest band, as generically expected as a result from the broken symmetry Goldstone et al. (1962). The absolute maximum is reached for 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 timeaveraged Bogoliubov energy, which is consistent with the fact that this mode will be the first one to resonate [] as one increases the interaction strength Lellouch et al. (2017). We emphasize that this most unstable mode differs from the one identified for the resonantlyshaken WannierStark 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 shakenlattice models; see Sec. III and Ref. Lellouch et al. (2017). 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 timeaveraged Bogoliubov dispersion).
V The driveninduced HarperHofstadter model
In this Section, we extend our previous analysis [Section IV] to a twodimensional setting, which has been exploited to realize the HarperHofstadter Hamiltonian in cold atoms Aidelsburger et al. (2011, 2013); Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2015b); Tai et al. (2017).
The model
We consider a twodimensional extension of the previous model [Eq. (13)], which we define by the Hamiltonian
(18)  
where annihilates a particle at lattice site , (resp. ) is the tunneling along the (resp. ) direction, and models onsite repulsive interactions. A linear (WannierStark) potential is introduced along the direction; besides the model features a moving lattice of amplitude , frequency and a phase dependence along both directions (). In the following, we set the resonance condition . This model, and other variants, have been realized in several recent ultracoldatom experiments Aidelsburger et al. (2011, 2013); Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2015b); Tai et al. (2017). Here, we focus on the “flux” configuration Miyake et al. (2013); Kennedy et al. (2015), which corresponds to setting 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 (), the effective Hamiltonian associated with Eq. (18) reads Goldman et al. (2014)
(19)  
with , , and . This corresponds to the HarperHofstadter model Hofstadter (1976), with a magnetic flux in each unit cell, and with different effective tunnelings along the and directions. As stated above, other choices for the movinglattice phase () lead to different fluxes per unit cell Aidelsburger et al. (2013, 2015b); Goldman et al. (2014). We recall that the HarperHofstadter model displays the wellknown “Hofstadter’s butterfly” spectrum Hofstadter (1976), a rich fractal structure that hosts Dirac semimetals and Chern insulating phases Kohmoto (1989). In this sense, at the singleparticle level, the timedependent model under consideration [i.e. Eq. (18) with ] is one of the simplest systems realizing artificial gauge fields and topological band structures for neutral atoms in 2D optical lattices Aidelsburger et al. (2011, 2013); Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2015b); Tai et al. (2017); we point out that a BEC was recently realized in the “flux” configuration of this model, which further motivates the present study Kennedy et al. (2015).
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 Hofstadter (1976); Powell et al. (2011).
To see this, one observes that the singleparticle effective Hamiltonian in Eq. (19) is periodic over a cell, so that its eigenstates are labelled by the quasimomentum defined in the reduced Brillouin zone . The spectrum features two energy bands (labelled )
(20) 
which are both twofold degenerate [see Figs. 5(ab)]. The ground state corresponds to the state with in the band, and is twofold degenerate ^{3}^{3}3Let us stress that the effective Hamiltonian in Eq. (19) obtained from Eq. (18), does not correspond to the Landau gauge, hence the difference compared to the more common version of the HarperHofstadter model whose spectrum is and groundstate in ..
In the presence of interactions, the ground state still features this twofold degeneracy Powell et al. (2011), as well as the momentum characteristic . 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 groundstates 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 . Yet, the exact positions of the stable and unstable regions are shifted, reflecting that the later are very sensitive to the groundstate in which the system is prepared.
The later observation that stable regions cannot be simply obtained by analyzing the behavior of the effective tunneling [Section III.3] can be related to the fact that these quantities never cancel simultaneously; indeed, in this specific 2D model, the timeaveraged Bogoliubov dispersion never vanishes, meaning that the simple “flatband” criterion of Section III.3 cannot be applied (a priori, there might always exist some modes fulfilling the resonance condition ). Consequently, we conjecture that the stable zones in Fig. 4 should rather be governed by the amplitude of the perturbation entering the underlying effective parametricoscillator model [i.e. the analogue of the quantity that was defined in Eq. 31, 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 groundstate dependent, since the Bogoliubov treatment leading to this effective parametricoscillator 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 groundstate. 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 mode (with respect to the groundstate). This can again be accounted for by the fact that it is the mode of highest timeaveraged Bogoliubov energy. To verify this, we calculate the timeaveraged Bogoliubov spectrum , 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 ; see Powell et al. (2011). The resulting Bogoliubov spectrum is plotted in Fig. 5(c). Consistently with the general features reported in Powell et al. (2011), this spectrum is made of four branches (the degeneracy observed at the singleparticle level being lifted) and indeed has its absolute maximum at in the upper branch.
Interestingly, this behavior is in fact expected to be very generic. For a model configuration leading to a flux , i.e. , we also find a most unstable mode at , consistently with the fact that it is the mode of highest timeaveraged Bogoliubov energy [see Ref. Powell et al. (2011) for the Bogoliubov spectrum in that case].
We conclude this Section by noting that the stability diagram of the driveninduced HarperHofstadter model [Fig. 6] displays relatively large stability regions, which extend up to significant values of the interaction strength , 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 Lellouch et al. (2017). Besides, our result that the stability diagram strongly depends on the prepared ground state appears as an important feature of these timemodulated 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 resonantlydriven band models. While the general method used to identify and characterize these instabilities had already been applied to nonresonant models Lellouch et al. (2017); Creffield (2009), we have hereby demonstrated its usefulness and flexibility to address the relevant case of resonant timemodulations. 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 ultracoldatom experiments Aidelsburger et al. (2011, 2013); Miyake et al. (2013); Kennedy et al. (2015); Aidelsburger et al. (2015b); Tai et al. (2017), 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 FloquetBogoliubov 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 underlying parametric instabilities and the fact that 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 parametricoscillator model [i.e. the analogue of the quantity in Eq. 31 in Appendix A]; in particular, in the presence of degenerate ground states, as in the “HarperHofstadter” case treated in Section V, the instability rates were found to be very sensitive to the actual groundstate 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 stronglycorrelated 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 (for experimental values of the drive frequency Lellouch et al. (2017)). When a continuum is present and the dispersion is unbounded, the most unstable mode is found (at lowest order) to fulfill the resonance condition . Therefore, the simple knowledge of the effective band structure [ and ] 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 timemodulated optical lattices, with reduced instabilities and losses, will constitute a first step towards the coldatom realization of stronglycorrelated states of matter with topological features, such as fractional Chern insulators Regnault and Bernevig (2011); Sterdyniak et al. (2013); Bergholtz and Liu (2013); Grushin et al. (2014); Möller and Cooper (2015). However, many open problems remain. First of all, this study was performed it the meanfield interaction regime, which is expected to give access to the shortterm dynamics. Longertime dynamics, including the possibility for thermalization, are dominated by nonlinear 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, higherband 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 Weinberg et al. (2016); Novicenko et al. (2016), 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 interband transitions Weinberg et al. (2015); StrÃ¤ter and Eckardt (2016); Reitter et al. (2017), is expected to lead to rich behaviors that still remain to be studied.
Acknowledgments
We thank M. Aidelsburger, M. Bukov and M. Cheneau for insightful discussions.
This research was financially supported by the Belgian Science Policy Office under the Interuniversity Attraction Pole project P7/18 “DYGEST”, the FRSFNRS (Belgium) and the ERC TopoCold Starting Grant.
References
 Oka and Aoki (2009) T. Oka and H. Aoki, “Photovoltaic hall effect in graphene,” Phys. Rev. B 79, 081406 (2009).
 Kitagawa et al. (2010) T. Kitagawa, E. Berg, M. Rudner, and E. Demler, “Topological characterization of periodicallydriven quantum systems,” Phys. Rev. B 82, 235114 (2010).
 Lindner et al. (2011) N. H. Lindner, G. Refael, and V. Galitski, “Floquet topological insulator in semiconductor quantum wells,” Nature Phys. 7, 490–495 (2011).
 Cayssol et al. (2013) J. Cayssol, B. Dora, F. Simon, and R. Moessner, “Floquet topological insulators,” Phys. Status Solidi RRL 7, 101–108 (2013).
 Mitrano et al. (2016) Matteo Mitrano, Alice Cantaluppi, Daniele Nicoletti, Stefan Kaiser, A Perucchi, S Lupi, P Di Pietro, D Pontiroli, M Riccò, Stephen R Clark, et al., “Possible lightinduced superconductivity in k3c60 at high temperature,” Nature 530, 461–464 (2016).
 Aidelsburger et al. (2017) M Aidelsburger, S Nascimbene, and N Goldman, “Artificial gauge fields in materials and engineered systems,” arXiv preprint arXiv:1710.00851 (2017).
 Goldman et al. (2014) N. Goldman, G. Juzeliunas, P. Ohberg, and I. B. Spielman, “Lightinduced gauge fields for ultracold atoms,” Rep. Prog. Phys. 77 (2014).
 Goldman and Dalibard (2014) N. Goldman and J. Dalibard, “Periodicallydriven quantum systems: Effective hamiltonians and engineered gauge fields,” Phys. Rev. X 4 (2014).
 Bukov et al. (2015a) M. Bukov, L. D’Alessio, and A. Polkovnikov, “Universal highfrequency behaviour of periodicallydriven quantum systems : from dynamical stabilization to floquet engineering,” Adv. Phys. 64, 139–226 (2015a).
 Eckardt (2016) A Eckardt, ‘‘Atomic quantum gases in periodically driven optical lattices,” (2016), arxiv:1606.08041v1.
 Aidelsburger et al. (2015a) M. Aidelsburger, M. Lohse, C. Schweizer, M. Atala, J. Barreiro, S. Nascimbene, N. R. Cooper, I. Bloch, and N. Goldman, “Measuring the chern number of hofstadter bands with ultracold bosonic atoms,” Nature Phys. 11 (2015a).
 Jotzu et al. (2014) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, “Experimental realisation of the topological haldane model,” Nature (London) 515, 237–240 (2014).
 Goldman et al. (2016) N. Goldman, J.C. Budich, and P. Zoller, “Topological quantum matter with ultracold gases in optical lattices,” Nature Physics 12, 639–645 (2016).
 Eckardt et al. (2010) A. Eckardt, P. Hauke, P. SoltanPanahi, C. Becker, K. Sengtock, and M. Lewenstein, “Frustrated quantum antiferromagnetism with ultracold bosons in a triangular lattice,” Europhys. Lett. 89, 10010 (2010).
 Struck et al. (2011) J. Struck, C. Olschlager, R. Le Targat, P. SoltanPanahi, A. Eckardt, M. Lewenstein, P. Windpassinger, and K. Sengstock, “Quantum simulation of frustrated classical magnetism in triangular optical lattices,” Science 333, 996–999 (2011).
 Görg et al. (2017) F. Görg, M. Messer, K. Sandholzer, G. Jotzu, R. Desbuquois, and T. Esslinger, “Enhancement and sign reversal of magnetic correlations in a driven quantum manybody system,” ArXiv eprints (2017), arXiv:1708.06751 [condmat.quantgas] .
 Aidelsburger et al. (2011) M. Aidelsburger, M. Atala, S. Nascimbene, S. Trotzky, Y. Chen, and I. Bloch, “Experimental realization of strong effective magnetic fields in an optical lattice,” Phys. Rev. Lett. 107 (2011).
 Aidelsburger et al. (2013) M. Aidelsburger, M. Atala, M. Lohse, J. Barreiro, B. Paredes, and I. Bloch, ‘‘Realization of the hofstadter hamiltonian with ultracold atoms in optical lattices,” Phys. Rev. Lett. 111, 185301 (2013).
 Miyake et al. (2013) H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton, and W. Ketterle, “Realizing the harper hamiltonian with laserassisted tunneling in optical lattices,” Phys. Rev. Lett. 111, 185302 (2013).
 Kennedy et al. (2015) C. J. Kennedy, W. C. Burton, W. C. Chung, and W. Ketterle, ‘‘Observation of boseeinstein condensation in a strong synthetic magnetic field, arxiv:150308243,” Nature Physics 11, 859–864 (2015).
 Tai et al. (2017) M Eric Tai, Alexander Lukin, Matthew Rispoli, Robert Schittko, Tim Menke, Dan Borgnia, Philipp M Preiss, Fabian Grusdt, Adam M Kaufman, and Markus Greiner, “Microscopy of the interacting harper–hofstadter model in the twobody limit,” Nature 546, 519–523 (2017).
 Tarnowski et al. (2017) Matthias Tarnowski, F Nur Ünal, Nick Fläschner, Benno S Rem, André Eckardt, Klaus Sengstock, and Christof Weitenberg, “Characterizing topology by dynamics: Chern number from linking number,” arXiv preprint arXiv:1709.01046 (2017).
 Kolovsky (2011) A. Kolovsky, EPL 93, 20003 (2011).
 Bermudez et al. (2011) Alejandro Bermudez, Tobias Schaetz, and Diego Porras, “Synthetic gauge fields for vibrational excitations of trapped ions,” Physical review letters 107, 150501 (2011).
 Hauke et al. (2012) Philipp Hauke, Olivier Tieleman, Alessio Celi, Christoph Ölschläger, Juliette Simonet, Julian Struck, Malte Weinberg, Patrick Windpassinger, Klaus Sengstock, Maciej Lewenstein, and André Eckardt, Phys. Rev. Lett. 109, 145301 (2012).
 Goldman et al. (2015) N. Goldman, J. Dalibard, M. Aidelsburger, and N. Cooper, “Periodicallydriven quantum matter: the case of resonant modulations,” Phys. Rev. A 91 (2015).
 Creffield et al. (2016) C. E. Creffield, G. Pieplow, F. Sols, and N. Goldman, “Realization of uniform synthetic magnetic fields by periodically shaking an optical square lattice,” New J. Phys. 18, 093013 (2016).
 Price et al. (2017) Hannah M Price, Tomoki Ozawa, and Nathan Goldman, “Synthetic dimensions for cold atoms from shaking a harmonic trap,” Physical Review A 95, 023607 (2017).
 Sias et al. (2008) C. Sias, H. Lignier, Y. Singh, A. Zenesini, D. Ciampini, O. Morsch, and E. Arimondo, “Observation of photonassisted tunneling in optical lattices,” Phys. Rev. Lett. 100, 040404 (2008).
 Mukherjee et al. (2015) Sebabrata Mukherjee, Alexander Spracklen, Debaditya Choudhury, Nathan Goldman, Patrik Öhberg, Erika Andersson, and Robert R Thomson, “Modulationassisted tunneling in laserfabricated photonic wannier–stark ladders,” New Journal of Physics 17, 115002 (2015).
 Aidelsburger et al. (2015b) M. Aidelsburger, M. Lohse, C. Schweizer, M Atala, J. T. Barreiro, S. Nascimbène, N. R. Cooper, I. Bloch, and N. Goldman, Nature Physics 11, 162–166 (2015b).
 Hofstadter (1976) Douglas R. Hofstadter, “Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields,” Phys. Rev. B 14, 2239–2249 (1976).
 Kohmoto (1989) Mahito Kohmoto, “Zero modes and the quantized hall conductance of the twodimensional lattice in a magnetic field,” Physical Review B 39, 11943 (1989).
 Möller and Cooper (2015) Gunnar Möller and Nigel R Cooper, “Fractional chern insulators in harperhofstadter bands with higher chern number,” Physical review letters 115, 126401 (2015).
 Regnault and Bernevig (2011) Nicolas Regnault and B Andrei Bernevig, ‘‘Fractional chern insulator,” Physical Review X 1, 021014 (2011).
 Sterdyniak et al. (2013) Antoine Sterdyniak, Cecile Repellin, B Andrei Bernevig, and Nicolas Regnault, “Series of abelian and nonabelian states in c¿ 1 fractional chern insulators,” Physical Review B 87, 205137 (2013).
 Bergholtz and Liu (2013) Emil J Bergholtz and Zhao Liu, “Topological flat band models and fractional chern insulators,” International Journal of Modern Physics B 27, 1330017 (2013).
 Grushin et al. (2014) A. G. Grushin, A. GómezLeón, and T. Neupert, “Floquet fractional chern insulators,” Phys. Rev. Lett. 112, 156801 (2014).
 Reitter et al. (2017) M. Reitter, J. NÃ¤ger, K. Wintersperger, C. StrÃ¤ter, I. Bloch, A. Eckardt, and U. Schneider, ‘‘Interaction dependent heating and atom loss in a periodically driven optical lattice,” (2017), arxiv:1706.04819.
 Lellouch et al. (2017) S. Lellouch, M. Bukov, E. Demler, and N. Goldman, “Parametric instability rates in periodically driven band systems,” Phys. Rev. X 7, 021015 (2017).
 Modugno et al. (2004) M. Modugno, C. Tozzo, and Dalfovo F., “Role of transverse excitations in the instability of boseeinstein condensates moving in optical lattices,” Phys. Rev. A 70, 043625 (2004).
 Kramer et al. (2005) M Kramer, C Tozzo, and F Dalfovo, “Parametric excitation of a BoseEinstein condensate in a onedimensional optical lattice,” Physical Review A 71, 061602 (2005).
 Tozzo et al. (2005) C. Tozzo, M. Kramer, and F. Dalfovo, “Stability diagram and growth rate of parametric resonances in boseeinstein condensates in onedimensional optical lattices,” Phys. Rev. A 72, 023613 (2005).
 Creffield (2009) C. E. Creffield, “Instability and control of a periodicallydriven boseeinstein condensate,” Phys. Rev. A 79, 063612 (2009).
 Bukov et al. (2015b) M. Bukov, S. Gopalakrishnan, M. Knap, and Demler E., “Prethermal floquet steady states and instabilities in the periodically driven, weakly interacting bosehubbard model,” Phys. Rev. Lett. 115, 205301 (2015b).
 Martin et al. (2017) J. Martin, B. Georgeot, D. GuÃ©ryOdelin, and D. L. Shepelyansky, “Kapitza stabilization of a repulsive boseeinstein condensate in an oscillating optical lattice,” (2017), arXiv:1709.07792.
 Michon et al. (2017) E. Michon, C. CabreraGutierrez, A. Fortun, M. Berger, M. Arnal, V. Brunaud, J. Billy, C. Petitjean, P. Schlagheck, and D. GueryOdelin, ‘‘Outofequilibrium dynamics of a bose einstein condensate in a periodically driven band system,” (2017), arXiv:1707.06092.
 Jaskula et al. (2012) JC Jaskula, Guthrie B Partridge, Marie Bonneau, Raphaël Lopes, Josselin Ruaudel, Denis Boiron, and Christoph I Westbrook, “Acoustic analog to the dynamical casimir effect in a boseeinstein condensate,” Physical Review Letters 109, 220401 (2012).
 Bilitewski and Cooper (2015a) T. Bilitewski and N. R. Cooper, ‘‘Population dynamics in a floquet realization of the harperhofstadter hamiltonian,” Phys. Rev. A 91, 063611 (2015a).
 Bilitewski and Cooper (2015b) T. Bilitewski and N. R. Cooper, “Scattering theory for floquetbloch states,” Phys. Rev. A 91, 033601 (2015b).
 Weinberg et al. (2015) M. Weinberg, C. Ölschläger, C. Sträter, S. Prelle, A. Eckardt, K. Sengstock, and J. Simonet, “Multiphoton interband excitations of quantum gases in driven optical lattices,” Phys. Rev. A 92, 043621 (2015).
 StrÃ¤ter and Eckardt (2016) C. StrÃ¤ter and A. Eckardt, “Interband heating processes in a periodically driven optical lattice,” (2016), arXiv:1604.00850.
 Quelle and Morais Smith (2017) A. Quelle and C. Morais Smith, “Resonances in a periodically driven bosonic system,” ArXiv eprints (2017), arXiv:1710.09680 [condmat.statmech] .
 Logan W. Clark (2017) Lei Feng Cheng Chin Logan W. Clark, Anita Gaj, “Collective emission of matterwave jets from driven boseeinstein condensates,” (2017), arXiv:1706.05560.
 Salerno et al. (2016) G. Salerno, T. Ozawa, H. M. Price, and I. Carusotto, “Floquet topological system based on frequencymodulated classical coupled harmonic oscillators,” Phys. Rev. B 93, 085105 (2016).
 Weinberg et al. (2016) P. Weinberg, M. Bukov, S. Vajna, L. D’Alessio, A. Polkovnikov, and M. Kolodrubetz, “Adiabatic perturbation theory and geometry of periodicallydriven systems,” (2016), arxiv:1606.02229.
 Novicenko et al. (2016) V. Novicenko, E. Anisimovas, and G. Juzeliunas, “Floquet analysis of a quantum system with modulated periodic driving,” arxiv:1608.08420 (2016).
 (58) The factor of stems from the fact that is defined in amplitude while physical observables involve squared moduli of wavefunctions.
 Choudhury and Mueller (2015) S. Choudhury and E. Mueller, “Transverse collisional instabilities of a boseeinstein condensate in a driven onedimensional lattice,” Phys. Rev. A 91, 023624 (2015).
 Landau and Lifshits (1969) L. Landau and E. Lifshits, Theoretical Physics  Vol. 1 (mechanics) (1969).
 StamperKurn and Ketterle (2001) Dan M StamperKurn and Wolfgang Ketterle, “Spinor condensates and light scattering from boseeinstein condensates,” in Coherent atomic matter waves (Springer, 2001) pp. 139–217.
 (62) The global shift of comes from the fact that excitation momenta are defined with respect to the groundstate, which is a state in the present case.
 Goldstone et al. (1962) J. Goldstone, A. Salam, and S. Weinberg, “Broken symmetries,” Phys. Rev. , 965â970 (1962).
 Powell et al. (2011) S. Powell, R. Barnett, R. Sensarma, and S. Das Sarma, “Bogoliubov theory of interacting bosons on a lattice in a synthetic magnetic field,” Phys. Rev. A 83, 013612 (2011).
 (65) Let us stress that the effective Hamiltonian in Eq. (19) obtained from Eq. (18), does not correspond to the Landau gauge, hence the difference compared to the more common version of the HarperHofstadter model whose spectrum is and groundstate in .
 (66) Interestingly however, we find that this term can play a role for higher values of , with non trivial effects arising from the noncommutation between the diagonal and the offdiagonal matrix in Eqs. 28.
Appendix A Analytical treatment of the resonantlydriven WannierStark ladder
We present the analytical derivation of instability rates for the resonantlydriven WannierStark ladder of Sec. III. For the sake of generalness, we consider an even more general situation where the modulation has a nontrivial uniform phase, and with possible additional transverse directions in the model. We thus deal with the Hamiltonian
(21) 
where annihilates a particle at lattice site and transverse position , denotes the tunneling amplitude of nearestneighbor hopping along the xdirection, describes a kinetic part along transverse directions (be it a lattice or a continuous one), and is the repulsive onsite interaction strength. The onsite potential term is made of a WannierStark ladder along the xdirection introducing an energy shift between consecutive sites, and a timeperiodic modulation along the xdirection of amplitude , phase , and frequency . The frequency modulation is chosen to be resonant with the offset () with denoting some integer.
We will perform all the theoretical analysis in the rotating frame (where translational invariance is restored), which is defined by the unitary transformation (with the position operator on the lattice).
a.1 Condensate dynamics and Bogoliubov equations
The singleparticle quasienergy spectrum of this model reads
(22) 
where is the effective tunneling (with ) and is the dispersion associated with transverse directions. The effective groundstate, which is taken as the initial state in which the system is supposed to be prepared, is a Bloch state of momentum , and if or if .
Taking this state, normalized to the condensate density , as our initial condition for the propagation of the timedependent GrossPitaevskii equation (tGPE), we find that the solution of the (tGPE) reads
(23) 
where we have introduced the functions
and
Plugging it into the BogoliubovDe Gennes equations, and taking advantage of the translational invariance to rewrite these in momentum space, we find that the dynamics of a momentum mode is governed by the equation
(24) 
where
a.2 Mapping on a parametric oscillator model
Proceeding as in Lellouch et al. (2017), we then introduce the two following changes of basis: first,
(25) 
where and with
(26) 
the timeaveraged Bogoliubov dispersion and
and then,
(27) 
This allows us to write the Bogoliubov equations under the form
(28) 
where is the timeaveraged Bogoliubov dispersion associated with the effective GPE, within the Bogoliubov approximation [see Eq. (26)], is a diagonal matrix of zero average over one driving period,
and is a (realvalued) function,
(29) 
As shown in Lellouch et al. (2017), Eqs. 28 are formally equivalent to a set of independent parametric oscillators (one for each mode ), since they describe harmonic oscillators of eigenfrequencies driven by a weak timeperiodic modulation of frequency . As can be anticipated from a RWA treatment of Eq. 28, 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, , i.e. (resonance condition). The later will result in a longterm explosion in the (stroboscopic) dynamics of the corresponding mode. Assuming that the harmonics can be treated independently, the instability rate of a given mode and a given harmonics can be analytically computed following the perturbative method developed in Landau and Lifshits (1969); Lellouch et al. (2017), and the total instability rate and most unstable mode straightforwardly inferred from Eqs. (3) and (4).
a.3 Explicit analytics in the strict resonant case
As an illustration, we now consider the simple case (i.e. ) and 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 in Eq. (29). We can thus disregard other harmonics, as well as the term which has no longterm stroboscopic contribution ^{4}^{4}4Interestingly however, we find that this term can play a role for higher values of , with non trivial effects arising from the noncommutation between the diagonal and the offdiagonal matrix in Eqs. 28., and we get the effective model
(30) 
with
(31) 
As shown in Refs. Lellouch et al. (2017); Bukov et al. (2015b), Eq. 30 is indeed formally equivalent to the equations of motion describing a parametric oscillator of eigenfrequency , driven by a perturbation of amplitude and frequency . The quantity 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 Landau and Lifshits (1969); Lellouch et al. (2017). At lowest order in , we find in this case
(32) 
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 at the onset of instability.
In the case where a continuous degree of freedom is present, we find two regimes :

If , one finds
(33) (34) 
If , we find
(35) (36)