Parametric Instability Rates in Periodically Driven Band Systems

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

Popular Summary

At very low temperatures, quantum matter can show radically different, and often counterintuitive, behavior compared to traditional liquids, solids, and gases. Such quantum phases of matter (e.g., superfluids and superconductors) are used in practical applications such as levitating trains and medical imaging. In the field of quantum simulation, researchers are attempting to engineer phases of matter that are even more exotic. One appealing approach is to start with a standard quantum fluid and “drive it” by shaking it, rotating it, or subjecting it to an external field. However, this method faces a serious difficulty: The phase is unstable. We have developed mathematical tools that reveal the origin and signatures of these instabilities that are directly relevant to current experiments.

Our rigorous theoretical methods offer simple analytical expressions that apply to a wide variety of driven systems, such as shaken atomic gases in optical lattices as well as irradiated solids. We identify the relevant instabilities that appear over very short times—in an explosive manner—due to the interplay between the external drive and interparticle collisions. These instabilities are caused by a buildup of collective excitations of the quantum gas upon driving. We present a simple but powerful method to calculate the instability rates associated with these excitations and identify clear manifestations (e.g., momentum distributions) that can be observed in experiments.

Understanding how quantum systems respond to an external drive, and the mechanisms through which such systems destabilize, will ultimately allow for the development of robust engineered quantum setups that are immune to instabilities and heating processes. This is not only crucial for our understanding of quantum matter but also for technological applications such as quantum computing.

Figure 1

The quasiparticle momentum distribution $n_q$ of a 2D shaken BEC trapped in a one-dimensional lattice (along $x$), and allowed to move along a continuous (tube) transverse direction ($y$), develops unique resonance structures in the presence of a periodic drive. Dominating the dynamics in the early stage of evolution, they signal the onset of instabilities and provide a clear experimental signature of the parametric resonance phenomenon explained in this work. The momentum distribution is shown after 60 driving periods ($t=60T$). Here, the distribution $n_q$ is divided by the system’s volume Vol, $\xi$ denotes the BEC healing length, and $a_x$ is the lattice spacing constant. The precise model and the corresponding parameters are the same as in Fig. 17.

Figure 2

Numerical (top panel) and analytical (bottom panel) stability diagram for the model described by Eq. (2). The instability rate $\mathrm{\Gamma }$ (expressed in units of the hopping amplitude $J$, see text) is plotted as a function of the modulation amplitude $K/\omega$ [defined in Eq. (2)] and interaction strength $g$, for a driving frequency $\omega =5J$. The vertical lines and letters summarize the agreement and disagreement regions, as discussed in Sec. 4c. A: Away from the zeros of ${\mathcal{J}}_0(K/\omega )$ and ${\mathcal{J}}_2(K/\omega )$, the analytics successfully capture the instability rates extracted from numerics. B: Agreement zone near the (close) zeros of ${\mathcal{J}}_0(K/\omega )$ and ${\mathcal{J}}_2(K/\omega )$, where both analytics and numerics predict a stable behavior. C/D: Close to a zero of ${\mathcal{J}}_0(K/\omega )$, the analytical perturbative approach breaks down. E: Potential quantitative disagreement when the contribution of the second harmonic has a weak amplitude; the instability is then partly due to higher harmonics (see Sec. 4b).

Figure 3

Numerical instability rate $\mathrm{\Gamma }$ as a function of interaction strength $g=U\rho$ and modulation amplitude $K/\omega$, for two values of the driving frequency ($\omega =10J$ and $\omega =J$). At large $\omega$ (i.e., outside the “low-frequency” regime, see text), the instability rates depend mostly only on $K/\omega$ and $g/{\omega }^2$, while in the low-frequency regime, instabilities can occur at infinitesimal interaction strength.

Figure 4

Time evolution of the energy density of the fluctuations (in units of $J$) over eight driving periods for two values of $g$. In the stable regime (top panel, $g=4$), the energy displays variations due to micromotion but no long-term growth, while in the unstable regime (bottom panel, $g=7$), the energy curve features exponential growth, defining a heating rate compatible with the maximal Lyapunov exponent $\mathrm{\Gamma }$ (see Fig. 5). The parameters are $A_q=1$ and $B_q=0$ (for the initial condition, see text), $\omega =5J$, and $K/\omega =7.5$.

Figure 5

Heating rate, obtained by fitting the exponential growth of the energy in Fig. 4. The resulting diagram is very similar to the instability diagram (Fig. 2), up to an overall numerical factor $\approx 2$ (see text).

Figure 6

Time evolution of the energy density of the fluctuations (in units of $J$) over eight driving periods for the initial conditions $A_q=1$ and $B_q=0$ (see text), for $\omega =5J$, $K/\omega =7.5$, and $g=5.7$. In this regime, where $\mathrm{\Gamma }$ is small, three regimes clearly show up: At very short times $t\ll T$, the initial growth is not described by our Floquet approach, as previously explained; at intermediate times $T\ll t\ll {\mathrm{\Gamma }}^{-1}$, the growth is linear and described by Eq. (18) (red dotted line); at large times $t\gg {\mathrm{\Gamma }}^{-1}$, the growth is exponential and described by Eq. (17) (red short-dotted line).

Figure 7

Top panel: Instability rate $\mathrm{\Gamma }$ as a function of the interaction strength $g$ for $\omega =5J$ and $K/\omega =7.5$, numerically computed from the Bogoliubov equations (cf. Sec. 2), taking into account only a few harmonics of the modulation in the function $h_q(t)$. As theoretically predicted from the expression of $h_q(t)$ [Eq. (20)], the odd harmonics do not contribute, while the first even ones successfully describe the full numerical calculation. Bottom panel: Comparing the analytical and numerical results for the instability rate $\mathrm{\Gamma }(g)$, using the same parameters as above. The red crosses correspond to the full numerical solution, while the blue dashed lines are numerically obtained by only taking the second harmonic of the drive into account. As generically observed in the A zones of the stability diagram (Fig. 2), the numerical instability rates are well captured by the perturbative analytical approach, implemented here up to order 1 (red dotted-dashed line) and 2 (blue dotted line).

Figure 8

Total instability rate $\mathrm{\Gamma }$ (red crosses) and instability rates $s_q$ of a few individual modes (dashed lines), numerically computed as a function of interaction strength $g$ for $\omega =5J$ and $K/\omega =7.5$. The instability rate associated with an individual mode $q$ is maximal for $E_{\mathrm{av}}(q)=\omega$ and decreases around it. The total instability rate is given by the envelope of all those curves.

Figure 9

Total instability rate and contributions of a few modes as a function of the interaction strength $g$, for $\omega =J$ (low-frequency regime) and $\omega =5J$ (outside the low-frequency regime), and for an increasing number of lattice sites $N=2$, 4, 6, 64. The total instability rate is the envelope of the curves associated with “available” modes in the system (for $N=2$, only $q=0$, $\pi$; for $N=4$, only $q=0,\pi /2,\pi ,3\pi /2\mathrm{...}$). The curve gets smoother when increasing the number of sites. Outside the low-frequency regime, the transition to instability is governed by the mode $\pi$ and already captured by a two-site system. In the low-frequency regime, the instability rate at low $g$, which is governed by a certain mode $q_{\mathrm{res}}$, is probed with more and more accuracy when gradually increasing $N$: Here, the two-site system is always stable since $q_{\mathrm{res}}>\pi$, the four-site system displays a transition since $q_{\mathrm{res}}>\pi /2$, and then this transition lowers when probing $q_{\mathrm{res}}$.

Figure 10

Numerical instability rate as a function of interaction strength $g=U\rho$ and modulation amplitude $K/\omega$, for $\omega =10J$ (i.e., outside the low-frequency regime) and for two different numbers of lattice sites in the transverse direction: $N_y=4$ (top panel) and $N_y=16$ (bottom panel).

Figure 11

Numerical instability rate as a function of interaction strength $g=U\rho$ and modulation amplitude $K/\omega$, for $\omega =5J$ (i.e., in the low-frequency regime) and for two different numbers of lattice sites in the transverse direction: $N_y=4$ (top panel) and $N_y=16$ (bottom panel).

Figure 12

Analytical instability rate as a function of interaction strength $g=U\rho$ and modulation amplitude $K/\omega$, for $\omega =10J$ (i.e., outside the low-frequency regime) and for two different numbers of lattice sites in the transverse direction: $N_y=4$ (top panel) and $N_y=16$ (bottom panel). This figure is to be compared with Fig. 10.

Figure 13

Numerical (top panel) and analytical (bottom panel) instability rates $\mathrm{\Gamma }$ as a function of the modulation amplitude $K/\omega$ and interaction strength $g$, for a high drive frequency $\omega =32J$.

Figure 14

Top panel: Instability rate as a function of the interaction strength $g$ for $\omega =32J$ and for two values of the modulation amplitude. As captured by the analytical expression Eq. (36), the instability rate increases linearly with $g$. Bottom panel: Instability rate as a function of the modulation amplitude $K/\omega$ for $\omega =32J$ and for two values of the interaction parameter $g$. As captured by the analytical expression Eq. (36), the instability rate increases quadratically with $K/\omega$ at low amplitudes.

Figure 15

Loss rate (exponential growth rate of the noncondensed density) for the 2D lattice model of Sec. 6a, for $\omega =5J$ and $N_y=16$. This rate is obtained by fitting the exponential growth of the noncondensed fraction, similarly to what was done in Fig. 5 for the energy.

Figure 16

Time dependence of the quasiparticle (phonon) density (top panel) and the condensate density (bottom panel) using the linearized (red line) and the conserving (blue line) approximation, for the 2D band model from Sec. 6b (with a lattice along the $x$ direction and the continuum along the $y$ direction). There are two discernible regimes: I and II (see text), delineated by a vertical dashed line. The vertical dotted line marks the time up to which the linearized mean-field dynamics and the WCCA agree well. The model parameters are $\rho /\mathrm{Vol}=17.34$, $U/J=1.011$ ($g/J=17.53$), $m=1/(20J)$, and healing length to $x$-lattice constant ratio $\xi /a_x=1.47$, $\omega /J=32.37$, and $K/\omega =1.0$. We use a total of $N=N_x{N}_y=5\times 10^4$ momentum points with $N_x=50$ and $N_y=1000$ to reach the thermodynamic limit. As an initial state, we choose the Bogoliubov ground state.

Figure 17

Quasiparticle momentum distribution function $n_q$ after a duration of $t=1000T$, for the 2D band model from Sec. 6b (with a lattice along the $x$ direction and continuum along the $y$ direction). The top panel shows the linearized mean-field model, while the bottom panel shows the conserving approximation (WCCA) to order $U$. The model parameters are the same as in Fig. 16. Here, we set the lattice constant $a_x=1$ to unity, while $\xi =1.47a_x$ denotes the BEC healing length.

Figure 18

Quasiparticle momentum distribution function $n_q$ after a duration of $t=10T$, for the 2D band model from Sec. 6b (with a lattice along the $x$ direction and a continuum along the $y$ direction), computed using the linearized dynamics, for different values of the relevant parameter $\sqrt{4|J_{\mathrm{eff}}|(4|J_{\mathrm{eff}}|+2g)}/\omega$; see Table 2. The suggested ellipses correspond to on-resonance modes [$E_{\mathrm{av}}(q)=\pi$], which are unstable, while the most unstable modes, which are revealed by the main peaks in the momentum distribution, are accurately described by our theoretical predictions (see Table 2).

Figure 19

Quasiparticle momentum distribution function $n_q$ for the 2D band model from Sec. 6b (with a lattice along the $x$ direction and a continuum along the $y$ direction), as a function of the parameter $\sqrt{4|J_{\mathrm{eff}}|(4|J_{\mathrm{eff}}|+2g)}/\omega$ and of $q_x$ (upper line) [resp. $q_{\perp }$ (lower line)], after integration over $q_{\perp }$ (resp. $q_x$), and computed as follows: (left panels) after a time $t=10T$ using the linearized dynamics, (middle panels) same time but using the WCCA, and (right panels) after a longer time $t=101T$ using the WCCA. The dotted blue line indicates our theoretical prediction for the most unstable mode (see Table 2), which is verified in all regimes: For $\sqrt{4|J_{\mathrm{eff}}|(4|J_{\mathrm{eff}}|+2g)}/\omega >1$, we find $q_{\perp }^{\mathrm{mum}}=0$ and $q_x^{\mathrm{mum}}<\pi$, as predicted by Eq. (37); for $\sqrt{4|J_{\mathrm{eff}}|(4|J_{\mathrm{eff}}|+2g)}/\omega <1$, we find $q_x^{\mathrm{mum}}=\pi$ and $q_{\perp }^{\mathrm{mum}}>0$, as given by Eq. (35). Note that, in the latter regime, the peaks are only apparent after a long time (here $101T$, right panel) because of the small values of the corresponding instability rates; see also Fig. 18, where these peaks are made visible. As already observed in Fig. 16, both linearized and WCCA agree very well on the position of the peaks, although the WCCA predicts a broadening of the instability zone due to nonlinear couplings between the modes. Note that at longer times (right panel), higher-order parametric resonances are also visible. The apparent “quantized” character of the resonant peaks for $t=101T$ is a numerical artifact due to the discretization of the variable $\sqrt{4|J_{\mathrm{eff}}|(4|J_{\mathrm{eff}}|+2g)}/\omega$ used in the numerical simulation.

Figure 20

Comparison between the numerically and analytically computed instability rate $\mathrm{\Gamma }$ in the case $K/\omega =2.4$ [zone C on the stability diagram of Fig. 2, where ${\mathcal{J}}_0(K/\omega )\approx 0$ and ${\mathcal{J}}_2(K/\omega )$ is large]. Since the effective Bogoliubov dispersion is flat, the analytical perturbative approach breaks down.

Figure 21

Comparison between the numerically and analytically computed instability rate $\mathrm{\Gamma }$ in the cases $K/\omega =5.52$ [(left panel), zone D of the stability diagram of Fig. 2, where ${\mathcal{J}}_0(K/\omega )\approx 0$] and $K/\omega =5$ [(right panel), zone E of the stability diagram of Fig. 2, where ${\mathcal{J}}_2(K/\omega )\approx 0$]. The disagreement in those cases can be due to a breakdown of the perturbative approach [when ${\mathcal{J}}_0(K/\omega )$ vanishes, (left panel)] and/or a vanishing of the second harmonic [when ${\mathcal{J}}_2(K/\omega )$ vanishes, (right panel)].