Out-of-equilibrium dynamics of a Bose Einstein condensate in a periodically driven band system

We report on the out-of-equilibrium dynamics of a Bose-Einstein condensate (BEC) placed in an optical lattice whose phase is suddenly modulated. The frequency and the amplitude of modulation are chosen to ensure a negative renormalized tunneling rate. Under these conditions, staggered states are nucleated by a spontaneous four wave mixing mechanism. The nucleation time is experimentally studied as a function of the renormalized tunnel rate, the atomic density and the modulation frequency. Our results are quantitatively well accounted for by a Truncated Wigner approach and reveal the nucleation of gap solitons after the quench. We discuss the role of quantum versus thermal fluctuations in the nucleation process and experimentally address the limit of the effective Hamiltonian approach.

Cold atoms provide powerful and versatile platforms for quantum simulators of many-body systems [1][2][3][4], and give access to the rich out-of-equilibrium dynamics of such systems. A remarkable progress for tunability has been achieved with the renormalization of the tunneling rate in time-dependent double-well potentials and optical lattices with the periodic shaking of the potential energy landscape [5][6][7]. It has opened many new possibilities for quantum simulations with the possibility to engineer effective Hamiltonians and simulate topological phases [8]. This includes the realization of the Hofstadter [9] and Haldane models [10] and the investigation of frustrated magnetism [11], to name a few recent examples.
This renormalization can be readily understood in the context of a one-body analysis. However, interactions play a key role in the quantum transition triggered by the modulation. An interesting situation in this respect is provided by the phase modulation of a 1D optical lattice in a regime for which the renormalized tunneling rate becomes negative. In this regime, the quantum gas experiences a dynamical instability associated with a spontaneous four wave mixing mechanism [5,12,13]. This generates a new quantum state, commonly called a staggered state, for which neighbouring sites acquire opposite phases [5,6]. This phase transition can be readily observed in momentum space after a time of flight expansion. The interference patterns [14] observed as a result of the spatial periodicity of a static optical lattice are modified: new peaks in between the former static peaks are observed [see Fig. 1(a)]. Interestingly this transition is not quantitatively accounted for with a standard meanfield approach.
In this article, we experimentally investigate the nucleation of such staggered states and compare our result with a time-dependent beyond mean-field approach. Our experiments are complementary to nucleation studies of vortices in a BEC using the rotating spoon technique [15][16][17], which provides another example of phase transition triggered by a dynamical instability [18,19]. In those lat-ter experiments, however, the kinetics of the transition could not be studied as a function of the density since the rotation weakens the transverse 2D confinement. We report hereafter a one-fold variation of the nucleation time of staggered states by varying the atomic density, and investigate experimentally and numerically the role of the renormalized tunneling rate and the modulation frequency on the transient out-of-equilibrium dynamics. We also explore the limit of the effective Hamiltonian approach in this context.
Our experiments have been realized on our rubidium-87 BEC machine that relies on a hybrid (magnetic and optical) trap [20]. It produces pure BECs of 10 5 atoms in the lowest hyperfine level F = 1, m F = −1. The 1D optical lattice is generated by superimposing two counterpropagating laser beams at 1064 nm (lattice spacing d = 532 nm) to the horizontal optical guide of the hybrid trap. The relative phase between the two arms of the lattice is controlled in time via synthesizers whose frequencies are imprinted on light using acousto-optic modulators. With a sinusoidal phase modulation of tunable amplitude ϕ 0 , the atoms experience the potential where ω ext accounts for the longitudinal confinement of the hybrid trap along the guide axis.
To obtain a theoretical understanding of the behaviour of a BEC in the presence of such a modulated potential, we first perform a gauge transformation to the comoving frame in which the lattice is periodically tilted instead of shaken. A single-band approximation can be justified in this latter representation provided the amplitude of the lattice is sufficiently strong such that V 0 > 4πϕ 0 mν 2 d 2 , which is the case in our experiments. The motion of the atoms along the lattice can then be modeled by means of a one-dimensional tight-binding Hamiltonian in which each well of the lattice is represented by one site. This  Hamiltonian is constituted by site-dependent on-site energies that are essentially given by the longitudinal confinement of the trap, as well as by an approximately siteindependent inter-site hopping matrix element J that can be tuned by varying the strength and/or wavelength of the lattice. As a result of the modulation, this inter-well tunneling rate J is renormalized by a Bessel function according tō where is the lattice characteristic energy. This result is readily derived from a one-body analysis [21] but turns out to remain valid in the presence of two-body interactions provided the latter as well as J are significantly smaller than hν. A qualitative picture of the impact of this renormalization on the physics of the system can be worked out perturbatively with the expression for the energy of the lowest band using the Peierls substitution: E 0 (k) = −2J cos(kd). ForJ > 0 the minimum of the band is located at k = 0, and the Fourier transform of the wave function consists in a comb of peaks centered about k = 0 with a spacing 2π/d. For J < 0, the minima are located on the border of the Brillouin zone at k = ±π/d. When the sign ofJ is changed through phase modulation, the system is therefore put in a metastable state. While this would not affect the mean-field dynamics of a BEC in the presence of a translational invariant lattice, any deviation from perfect homogeneity in the condensate wavefunction or the lattice will give rise to a shrinking amplitude of the periodic condensate mode and to an exponentially increasing population of staggered modes at k = ±π/d. The mechanism for the nucleation of the staggered states has been identified as a dynamical instability also referred to as a parametric instability triggered by interactions [22]. Going beyond the mean-field description of the condensate, deviations from translational invariance are naturally present in form of quantum fluctuations by means of which the actual many-body dynamics of the quantum gas at hand can be conveniently modeled and understood. As a consequence, a quantum depletion of the condensate is expected to arise in the course of time evolution, followed by the formation of an incoherent, thermal-like, cloud of atoms within staggered modes near the border of the Brillouin zone [23].
To obtain a qualitative understanding of this dynamical instability, we follow Ref. [22] and neglect the presence of the external longitudinal confinement for a moment. The system is then modeled by a homogeneous Bose-Hubbard Hamiltonian whereb l denotes the annihilation operator associated with the Wannier function centered on the l-th well of the lattice and U accounts for a site-independent twobody interaction strength. Assuming the presence of a perfectly homogeneous condensate at t = 0 exhibiting n atoms per site, we make the Bogoliubov ansatz [24] b l (t) = √ n + d 2π where µ = −2J + nU is the chemical potential of the condensate. Linearizing the resulting equation for the thereby introduced de-excitation operatorsΛ(k, t) yields whose solution evolves sinusoidally in time according tô sin Ω k t with the Bogoliubov phonon frequencies The time-dependent population of non-condensed modes can then be directly evaluated through A dynamical instability occurs whenJ becomes negative since Ω k becomes imaginary for each quasiparticle mode that satisfies the relationJ(1 − cos kd) < 0 < J(1 − cos kd) + nU , which according to Eq. (6) implies that tiny initial populations of such modes experience an exponential growth with time. For nU > −4J, this growth is most pronounced for the staggered Bogoliubov modes defined by kd = ±π, which describe an antiperiodic Bloch function within the lattice; the population of these staggered modes grows with the Lyapunov exponent λ = [−8J(2J + nU )] 1/2 / . This expression has stimulated our experiments described below, which explore how the nucleation time of the staggered states depend on the value ofJ and on the atomic density.
To allow for a more quantitative comparison with the experiment, we performed numerical simulations of the time evolution of the condensate by using the Truncated Wigner (TW) method [25]. This method approximately accounts for the effect of quantum fluctuations about the initial coherent state of the condensate. Applied to the Bose-Hubbard Hamiltonian (2), it essentially amounts to sampling the time evolution of the quantum bosonic many-body state in terms of classical trajectories that evolve according to a discrete Gross-Pitaevskii equation.
For the practical implementation of the TW method, we take into account the presence of the longitudinal confinement of the hybrid trap, which gives rise to sitedependent on-site energy terms that have to be added to the Bose-Hubbard Hamiltonian (2), as well as to a spatially inhomogeneous (Gaussian-like) condensate wavefunction which is calculated through imaginary time propagation. We furthermore incorporate the fact that the interaction energy within each site generally scales in a non-quadratic manner with the population of that site, due to the presence of a transverse Thomas-Fermi profile within the pancake-shaped lattice wells [26]. However, the presence of thermal fluctuations in the initial state (which can also be taken into account by the TW method) are neglected, i.e., we make the simplifying assumption that we start with a perfect BEC at T = 0 K.
The transition to staggered states is quantitatively characterized by the nucleation time representing the instance at which the population of a periodic condensate state within the shaken lattice becomes less significant than the population of a staggered state with antiperiodic nature. To this end, we define by P 0 the population of the central condensate peak located at k = 0, and by S ± the population of staggered states centered about ±h/(2d) in momentum space. In the numerical simulations, P 0 and S ± are determined by integrating the momentum-space density ρ(k) = n(k)a || exp(−k 2 a 2 || )/(π 1/2 d) within the intervals −π/2d < k < π/2d and −π/2d < k ∓π/d < π/2d, respectively, where n(k) is the 2π/d-periodic momentum density that results from the TW simulation of the Bose-Hubbard dynamics [see the shaded areas in Fig. 2(a)]. The nucleation time is then defined by the time beyond which S + +S − exceeds P 0 . It corresponds to the instance at which the peaks at k = ±π/2d have the same accumulated contrast as the peak at k = 0 in the experimental absorption images.
The first experimental study that we performed deals with the nucleation time of staggered states for different values of the renormalized tunneling rateJ. We specifically explored the dynamics at the border of the zone for whichJ becomes negative. The Lyapunov exponent predicts a divergence of the time required to populate the staggered states. We indeed observed a strong increase of the nucleation time when we approach an amplitude of modulation that corresponds to the zero of the Bessel function [see Fig. 1(c)]. Good agreement is obtained with TW computations of the nucleation time indicated by the solid (dashed) lines in Fig. 1(b), which were conducted according to the above prescriptions assuming the presence of 10 5 (5 × 10 4 ) atoms in the condensate. As these simulations were performed at T = 0, thermal fluctuations seem to play a minor role in this set of data for which the temperature was indeed so low that it could not be measured by time-of-flight methods.
The evolution of the atomic gas in momentum space is illustrated in Fig. 1(a) which shows a sequence of timeof-flight absorption images taken after various evolution times. We clearly see the transition from an initially coherent BEC, characterized by a sharp central peak at k = 0 and by two side peaks at k = ±2π/d, to an incoherent thermal cloud that oscillates between two maxima at k = ±π/d. These latter peaks are significantly broader than the initial condensate peaks, which is indicative of an effective increase of the temperature as a consequence of the dynamical instability mechanism [22,[27][28][29]. This is also observed in the TW simulations, as is seen in Fig. 2(a) which shows a snapshot of the numerically computed momentum distribution of the atomic gas after the nucleation process. Note that the sharp coherent BEC peaks at k = 0 and ±π/d are still present in the numerical simulations, in contrast to the experiment where they are nearly completely washed out after the nucleation. We attribute this to the fact that our TW approach is based on a single-band approximation and does not account for the effect of quantum fluctuations in the transverse degrees of freedom within the lattice, which are expected to provide further significant contributions to the quantum depletion of the condensate [30].
The TW approach allows one to obtain complementary insight into the nature of the staggered states that would not be easily accessible in the experiment. This concerns, in particular, the behaviour of lattice site occupancies. While their mean values n l = b † lb l do not display any notable feature in the course of time evolution, their rms widths ( n 2 l − n l 2 ) 1/2 dramatically increase as soon as staggered states become significantly populated. This is illustrated in Fig. 2(b) where the mean lattice site occupancies (averaged over 10000 trajectories) are plotted together with the occupancies that were obtained from a single trajectory (shown in red), at an evolution time t = 7.5 ms that exceeds the nucleation process.
The occurrence of pronounced spatial fluctuations is strongly reminiscent of gap solitons [31] and indicates that the formation of staggered states in momentum space is accompanied by the generation of solitons. This interpretation is consistent with the Bogoliubov mode analysis in a 1D shaken optical lattice of Ref. [22]. It is further confirmed by the fact that the spatial extension of the density peaks is indeed on the order of the theoretical prediction σ ≃ 2.6d[2|J|/(n max U )] 1/2 ≃ 0.7d for the full width at half maximum of a gap soliton according to Ref. [31] with n max ≃ 6000 the maximal occupancy within a lattice site. We should keep in mind, however, that in the actual (three-dimensional) experimental context, these gap solitons are expected to quickly disintegrate into vortices and vortex rings through the snake instability [32]. This effect is not accounted for in our numerical simulations which are based on a one-dimensional approximation for the optical lattice.
To explore the dependence of the nucleation time on the atomic density, we performed experiments where we drastically diminished the intensity of the vertical beam of the crossed dipole trap after the BEC production. In this manner, we could vary the atomic peak density between 5.5 × 10 13 at.cm −3 and 10 13 at.cm −3 [33] and observe a one fold increase of the nucleation time. As shown in Fig. 3(a), the experimental results are in fairly good agreement with TW simulations, except for high peak densities n peak 5 × 10 13 at.cm −3 . Indeed, to get such high densities we compress adiabatically the trap and therefore increase the temperature. Our results correspondingly suggest that thermal fluctuations are playing here a more important role in the nucleation process.
Finally, we investigate in Fig. 3(b) the scaling of the nucleation time with the driving frequency ν where we adapt the shaking amplitude ϕ 0 = (0.889π/ν) × 1.5 kHz such that the argument of the Bessel function is kept constant according to Eq. (1) yieldingJ ≃ −0.33J. As we have the same time-averaged Bose-Hubbard Hamiltonian (2) for all ν, the nucleation time is found to vary only rather weakly with the driving frequency. Note, however, that this behaviour is expected to change for ν V 0 /(4πϕ 0 νmd 2 ) ≃ 6.5 kHz where the single-band approximation can no longer be justified within the comoving frame. Indeed, a significant population of higher bands is expected to occur in this latter high-frequency regime as a consequence of near-resonant transitions, which in turn should lead to a drastic decrease of the nucleation time. This is exemplified in an experimental study performed for ν = 14 kHz and ϕ 0 = 0.028π (all other parameters were chosen as usual) for which a transition to staggered states is found after a few ms, despite the fact thatJ > 0 for this combination of parameters [26]. We attribute this instability of the condensate to a near-resonant transition from the ground band to the first excited band in the lattice. The latter exhibits a maximum at k = 0 in the Brillouin zone, which gives rise to a negative inter-site hopping matrix element.
In summary, we presented an experimental and theoretical study of the formation of staggered states within a BEC that is subjected to a periodically modulated optical lattice. Our measurements are in good agreement with numerical simulations using the TW method at T = 0 K which are based on a single-band description of the lattice, except for high densities where we have evidence that thermal fluctuations are playing an important role.
The present experimental setting can therefore be exploited in order to yield a quantitative diagnostic tool for determining interaction and fluctuation effects in BECs. Further studies in the high-frequency regime will be useful and interesting in order to explore the interplay with near-resonant inter-band transitions in more detail.
We thank J. Dalibard, N. Goldman, and M. Oberthaler for useful discussions. This work was supported by Programme Investissements d'Avenir under the program ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT, and by the Programme Hubert Curien 2016. In this supplementary material, we provide some technical details on the implementation of the Truncated Wigner method which was used to perform the numerical simulations. We furthermore provide a reminder on our lattice depth calibration technique [1]. We finally show results of phase modulation experiments performed in an extended range of parameters in both modulation frequency and lattice depth with respect to the experiments presented in the main article and discuss the nucleation of staggered states in each case.

I. IMPLEMENTATION OF THE TRUNCATED WIGNER METHOD
The starting point for the implementation of the Truncated Wigner method is the effective many-body Bose-Hubbard Hamiltonian which describes the dynamics of a Bose-Einstein condensate within a periodically shaken lattice. Here we employ a single-band approximation within the reference frame that is comoving with the shaken lattice. The shaking is incorporated by means of a periodically time-dependent Peierls phase within the inter-site hopping matrix elements of the lattice. The on-site energies account for the presence of the longitudinal harmonic confinement with the oscillation frequency ω || . For lattice strengths V 0 = sE L with s 2, we can approximately represent the Wannier function within the lth well by the Gaussian Φ l (x) = 1 where is the oscillator length associated with the frequency ω 0 = s 1/2 E L / of oscillations within each well of the lattice. The effective hopping parameter between adjacent wells can be determined from a semiclassical Wentzel-Kramers-Brillouin (WKB) ansatz [2,3]. We obtain with the dimensionless imaginary action where we introduce the parameter η = exp − 1 2s 1/2 − 1 2s 1/2 (14) and the incomplete elliptic integral of the second kind E(ϕ, k) = ϕ 0 1 − k 2 sin 2 θdθ. Note that this WKB ansatz is based on the approximate expression E 0 = (1 − η)s 1/2 ω 0 /2 for the ground state energy within each well, which is obtained from the expression (10) for the Wannier function using first-order perturbation theory.
A complication is introduced by the fact that the lattice wells in the experiment are not truly one-dimensional but rather exhibit a pancake shape, as the transverse confinement frequency ω ⊥ of the trap is comparable to the longitudinal one ω || and hence much larger than ω 0 . For the total atom numbers at hand, the BEC exhibits a parabolic Thomas-Fermi profile in its transverse density distribution within lattice wells that are located near the center of the trap. The interaction energy within such a well therefore scales in a non-quadratic manner (namely ∝ n 3/2 l ) with the occupancy n l of that well, which implies that the effective on-site interaction parameters U l appearing in the Bose-Hubbard Hamiltonian (7) scale as n −1/2 l with the occupancies of the corresponding sites. For weakly populated lattice sites that are located at the edge of the atomic cloud, on the other hand, we can justify the perturbative approximation g ≃ 2 ω ⊥ a s [4] for the effective one-dimensional interaction strength with a s ≃ 5.3 × 10 −9 m the s-wave scattering length of 87 Rb atoms, which yields U l ≃ 2 ω ⊥ a S /( √ 2πa 0 ) for those sites.
In order to simultaneously account for weakly and strongly occupied sites of the lattice, we employ a heuristic interpolation formula between the perturbative and the Thomas-Fermi regime [5]. We furthermore make the simplifying assumption that the lattice site occupancies n l remain fairly constant and are not substantially altered in the course of time evolution. While this constraint appears to be well respected on average, significant fluctuations of the lattice site occupancies about their average values may nevertheless arise if the condensate undergoes a transition to a staggered state [as seen in Fig. 2(b) of the main article].
Having thereby determined all relevant parameters of the Bose-Hubbard system (7), we can then derive the discrete nonlinear Schrödinger equation (16) by means of which the classical fields ψ l that sample the quantum many-body state of the system have to be evolved with time in the framework of the Truncated Wigner method. For the initial state we assume the presence of a perfect BEC within the optical lattice, which is in an idealized manner described by a coherent state in the classical field space. This leads to the choice ψ l (0) = φ l + χ l for the initial value of ψ l where φ l corresponds to the condensate wavefunction within the (unshaken) lattice at t = 0. In practice, φ l is determined by imaginary-time propagation of the Gross-Pitaveskii equation that describes the condensate at t = 0 [which is nearly identical with Eq. (16), except that it exhibits real hopping and the usual U l |ψ l | 2 ψ l interaction term]. This thereby yields n l = |φ l | 2 for the determination of the onsite interaction strengths according to Eq. (15). χ l is a complex random number drawn from a Gaussian probability distribution which is centered about the origin in the complex plane and yields the variance |χ l | 2 = 1/2. Effects due to the presence of finite temperature and initial quantum depletion within the atomic cloud are therefore neglected in this study. Expectation values of one-body operators are then evaluated according to the prescription where ψ * l (t)ψ l ′ (t) denotes the statistical average (over Truncated Wigner trajectories) of the expression ψ * l (t)ψ l ′ (t). This yields in particular b † lb l t = |ψ l (t)| 2 − 1/2 for the average occupancy of the site l at time t. A similar subtraction of a "half-particle" is needed in order to determine the average distribution of atoms in momentum space. We perform a precise calibration of the lattice depth using the out-of-equilibrium dynamics of a chain of BECs in an optical lattice following the method we demonstrated in [1]. To do so, we first load adiabatically a BEC in the lattice, creating a chain of BECs trapped at the bottom of the lattice sites. At t = 0, we suddently shift the phase of the lattice, which triggers the centerof-mass motion of the atomic wave packets in each well (shift of θ 0 = 90 • here). After a given holding time in the shifted lattice, we perform a 25 ms time-of-flight. The experimental absorption images (see Fig.4) show the interferences of the different wave packets located in each lattice site. The interference figure is centered on the central interference peak (0th order) corrresponding to wave packets that were released while being at rest in the shifted lattice, i.e. at the turning points of the oscillatory motion. We observe oscillations of the population in the 0th order from which we extract the period of the center-of-mass motion of the atomic wave packets. This period can be directly related to the depth of the optical lattice. Indeed, the intrasite dipole mode is coupled to a two phonon transition between the ground state band and the second excited band at k = 0:

II. LATTICE DEPTH CALIBRATION
In this way, we find a lattice of depth 1.9 E L corresponding to the measured period T 0 =153 µs. Interestingly, we have shown that this kind of measurement is robust against the value of the phase shift of excitation θ 0 , the atom-atom interaction strength, and the external confinement superimposed to the lattice. For these values of the phase shift and the lattice depth, tunneling plays an important role in the wave packet dynamics insofar as an important fraction of the wave packets tunnels to the neighboring lattice sites. It explains for instance the population of the -1 order for an evolution time in the shifted lattice of 40 µs but also more generally the asymmetry of the populations in the ± 1 orders. We have studied in details such effect in Ref. [1].

A. Frequency range
We extend the study of the effect of the modulation frequency ν on the nucleation of staggered states performed in the main article to determine the frequency range in which such states appear. We first use a lattice of depth 1.9 E L , characterized by a center-of-mass oscillation of frequency ν c = 6.55 kHz (see above) and perform phase modulation experiments at different modulation frequencies below ν c [see Fig.5 (1-4)]. We maintain the product ϕ 0 ν constant with ϕ 0 the amplitude of modulation, yielding an effective tunneling rateJ that is constant and negative. The argument of the Bessel function [see equation (1) of the main article] is chosen equal to 3.24.
At the lowest modulation frequencies (2 and 3 kHz), we clearly observe the nucleation of staggered states with the population of staggered modes at momentum p = ±h/2d, whereas for larger modulation frequencies (ν = 5 kHz), we do not significantly populate the staggered modes [see Fig.5(c4)]. More specifically, defining by ν c the centerof-mass oscillation frequency within a lattice well, we observe the population of staggered modes for modulation frequencies up to typically ν c /2, meaning for modulation frequencies that are not too close to excitations frequen-cies towards the excited bands. The experiments presented in the main article on the effect of the modulation frequency (see Fig. 3 of the main article) were performed with a lattice of depth 2.6E L , corresponding to a centerof-mass oscillation frequency ν c = 8 kHz. As those experiments were performed for modulation frequencies up to 4 kHz, we could see the population of the staggered modes for each chosen modulation frequency.
Additionally, we perform modulation experiments in a lattice of larger depth, characterized by a center-ofmass ocillation of frequency ν c = 10.6 kHz. Intriguingly, we could observe the nucleation of staggered states [see Fig.5 (5)] for a frequency of modulation (ν = 14 kHz) larger than ν c and a very low amplitude of modulation (ϕ 0 = 0.028π). For such parameters, the argument of the Bessel function is equal to 0.94 and would correspond to a positive effective tunneling rate in a single-band theory. However, the modulation frequency is close to the onephonon line from the lowest to the first excited band. Thus, the theoretical treatment used in the main article, which is based on a single-band approximation, is no longer applicable. The study of the appearance of staggered states in this frequency regime requires a more involved treatment taking into account higher bands.

B. Lattice depth range
We perform phase modulation experiments for different lattice depths ranging from 1.2E L to 3.2E L and observe the nucleation of staggered states in this range of lattice depth (see Fig.6). The product ϕ 0 × ν is kept constant and corresponds to a negative effective tunneling rateJ. The dynamics for the population of staggered modes look similar to the dynamics shown in the main article (see Fig. 1 in the main article), but for the lattice depth 1.2E L . For the lattice depth 1.2E L , the measurement is performed at a lower modulation frequency and as a consequence with a larger amplitude of modulation. For such an amplitude of modulation, the associated classical phase space starts to exhibit chaotic zones, which is probably responsible for the disrupted aspect of the dynamics observed in this case [see Fig.6(a)].