Driven Bose-Hubbard Model with a Parametrically Modulated Harmonic Trap

We investigate a one-dimensional Bose-Hubbard model in a parametrically driven global harmonic trap. The delicate interplay of both the local atom interaction and the global driving allows to control the dynamical stability of the trapped quantum many-body state. The mechanism is illustrated for weak interaction by a discretized Gross-Pitaevskii equation within a Gaussian variational ansatz, yielding to a Mathieu equation for the condensate width. The parametric resonance condition can be tuned by the atom interaction strength. For stronger interaction, this mechanism is confirmed by results of the numerically exact time-evolving block decimation scheme. The global modulation also induces an effective time-independent inhomogeneous hopping strength for the atoms.

Novel concepts of driven quantum many-body systems can be studied in atomic quantum gases, see Refs. [23][24][25][26][27] for recent reviews. A trapped Bose-Einstein condensate (BEC) with weak interactions is well described by the mean-field Gross-Pitaevskii (GP) equation. In absence of any additional optical lattice, a homogeneous BEC in a time-dependent setup has been considered in different constellations for a long time.
In an early work, Castin and Dum analytically studied a homogeneous BEC in a parametrically modulated harmonic trap [28]. They showed that the driving induces a parametric instability in the global motion of the condensate which gets depleted exponentially fast and noncondensed modes become dominantly populated due to this effective "heating." The effect of a parametrically driven trap potential was also studied in Ref. [29] within the GP approach. It was shown that the dynamics of the condensate wave function is described by the classical Mathieu equation of a parametrically forced oscillator, by which one obtains stability criteria. In another mean-field study, the effect of a time-dependent scattering length on the collective motion of a BEC was studied in Ref. [30].
In the presence of an optical lattice, the BEC is described by the Bose-Hubbard model which is known to have the two distinct phases of a superfluid or a Mott insulator. Jaksch et al. [31] have investigated the case of a Bose-Hubbard model with a time-dependent lattice depth which leads to a variation of both the on-site interaction and the hopping amplitude. Starting out from the superfluid phase, the atoms are driven to the Mott insulator phase and converted there into molecules. Eventually, the melting of the molecular Mott insulating phase produces a molecular superfluid [31].
An interesting regime which is less explored is realized when a strongly interacting gas in the Mott phase is exposed to a time-dependent external driving of the global trapping potential. When a system is driven parametrically, it exchanges energy with the driving field and, in principle, can be heated to infinite temperature [49][50][51]. On the other hand, the parametric oscillator has regions of dynamical stability as well. So the natural question arises of how does a strong atomic interaction affect the stability of a globally parametrically driven quantum many-body system? Can strong short-range interaction stabilize a quantum gas in a parametric trap which would otherwise be unstable? In turn, can we obtain information on the atomic interaction by externally tuning the system to an unstable dynamical state?
In this work, we show that a global parametric modulation of the trapping potential, which does not have to be tuned to local properties, can be used to control the stability of the interacting quantum gas in an optical lattice. In particular, the global dynamics of the quantum many-body system in a parametrically modulated trap can be stabilized or destabilized by tuning the atomic interaction strength. Conversely, locating the onset of the instability can be used to determine the atom interaction strength. To illustrate the mechanism, we investigate the parametrically driven Bose-Hubbard model with repulsive interaction in two regimes. First, we consider the regime of weakly interacting atoms in the lattice in the presence of a parametrically modulated global trap. This can be treated by a mean-field Gross-Pitaevskii ansatz for the condensate wave function and is supported by a numerically exact treatment in terms of the time-evolving block-decimation (TEBD) method. Second, we aim to investigate the interplay of the strongly interacting quantum gas in the Mott regime with an additional parametrically modulated trap. To this end, we have calculated the time-dependent dynamics in this regime numerically by the TEBD approach. We find that the parametric driving leads to a breathing of the width of a local central Mott region which becomes resonant at frequencies which are shifted as compared with the noninteracting case. In the Mott regime, energy absorption is increasingly suppressed due to the strongly reduced compressibility of the Mott region.
After introducing the underlying driven Bose-Hubbard model in Sec. II, we present our mean-field analysis for parametric resonance based on a discretized GP equation in Sec. III. To go beyond the weak-interacting regime, we show our complementary results for strong interactions based on the exact numerical TEBD in Sec. IV. The connection between periodically driven harmonic trap and site-dependent hopping is clarified in Sec. V. We summarize our work in Sec. VI.

II. MODEL
We consider a one-dimensional Bose-Hubbard model with a global harmonic potential with a time-dependent curvature V (t) = V 0 + δV sin t. The potential has a time-averaged curvature V 0 , which is parametrically modulated with the strength δV and the frequency . The model Hamiltonian reads, withh = 1, where J is the hopping amplitude and U is the on-site interaction strength. Furthermore, b (b † ) are the bosonic annihilation (creation) operators at site , and n = b † b denotes the local occupation number operator. We consider a lattice with M sites loaded with N bosonic atoms. Thus, the lattice center is located at 0 = (M − 1)/2.

III. QUANTUM MANY-BODY PARAMETRIC RESONANCE IN MEAN-FIELD REGIME
In the noninteracting limit U = 0 and for no driving, the system can be mapped exactly to a discretized quantum harmonic oscillator with frequency ω 0 = 2 √ J V 0 and unity mass with a Gaussian ground state. When the parametric driving is switched on, the parametric resonance at n = 2ω 0 produces regions of instability in the parameter space [21,22] with diverging position and momentum variances. The driven single-particle problem is still exactly solvable in terms of the Mathieu equation with its known stability diagram. For interacting particles, this is no longer possible. To elucidate the impact of quantum many-body interactions on the parametric resonance, we first consider the mean-field regime.

A. Mean-field description
A dilute, weakly interacting atom gas at zero temperature is described by the mean-field Lagrangian density obtained from the Hamiltonian (1), with the mean-field condensate wave function |ψ = ψ b † |0 . By extremizing the Lagrangian with respect to ψ (ψ * ), one arrives at a discretized version of the Gross-Pitaevskii equation. To obtain an analytic approximation for the time evolution of the boson gas, we use a Gaussian trial wave function [52] with a time-dependent width α ≡ α(t) and β ≡ β(t) and minimize L with respect to α and β.
In the following, we consider the regime J V 0 , where the condensate is extended over many sites, i.e., α 1. Then, all sums over can be approximated by continuous integrals and the Lagrangian takes the form The Euler-Lagrange equations ∂ x L = d dt ∂ẋL for x = α, β provide the equations of motionα = 4J γ αβ and with γ = 1 4α 2 + α 2 β 2 and J γ = J e −γ . Aiming at a linear stability analysis, we expand the width α(t) = α 0 + δα(t) in terms of small deviations δα(t) from its equilibrium value α 0 , i.e., we assume α 0 |δα(t)| and linearize Eq. (5) with respect to δα. Taking into account correspondingly β(t) = β 0 + δβ(t) with β 0 = 0 we find that the stationary solution α 0 follows implicitly from 043604-2 2π and that the deviation δα from equilibrium obeys (6) where both the hopping J = J e −1/4α 2 0 and the driving strength δV = δV (1 + 1/2α 2 0 ) are renormalized. Furthermore, the atomic interaction renormalizes the potential curvature such that It ranges from V 4V 0 in the noninteracting limit to V 3V 0 in the Thomas-Fermi limit, i.e., when the kinetic term can be neglected. Equation (6) is the well-known Mathieu equation with an additional time-dependent-force term. Note that the inhomogeneity does not influence the parametric resonance condition [35]. Thus, for δV V , δα(t) exhibits a resonant behavior when the parametric resonance condition n = 2ω with the resonance frequency ω = 2 √ J V is fulfilled for n = 1,2, . . ..

B. Static trap
In Fig. 1(a), we show the stationary condensate width α 0 as a function of the interaction strength for different potential curvatures. Below UN J , the condensate width α 0 (J /V 0 ) 1/4 is mainly determined by the potential curvature V 0 and only gradually increases with UN. For UN > J , the U term becomes comparable in size to the J term which leads to a steeper growth of α 0 . In fact, the condensate width behaves as α 0 ∼ U 1/3 in the Thomas-Fermi limit. Moreover, in Fig. 1(b), we show the condensate width α 0 as a function of the static potential curvature V 0 . In all cases, we find an algebraic decrease of α 0 ∼ V −1/η 0 with increasing V 0 , where 3 η 4. In the noninteracting limit we find η = 4, whereas in the strongly interacting limit we obtain η = 3.
The resonance frequency ω is also modified by the condensate interaction. In Fig. 2 (main), we depict the dependence of ω on the interaction strength for various potential curvatures. For ease of comparison, we show the resonance frequency scaled to √ J V 0 = ω 0 /2. For strong interaction U J /N, the resonance frequency turns out to be independent of the interaction U and approaches ω = 2 √ 3J V 0 , where J J is satisfied. In the opposite limit of noninteracting bosons, J = J is not satisfied per se. This is especially the case for a steep potential when V 0 is so large that α 0 is only of the order of several lattice sites. In this regime, the resonance frequency is given by In the inset of Fig. 2, we show the resonance frequency as a function of the potential curvature for different interaction strengths. For a small (large) enough interaction strength U the resonance frequency ω weakly decreases (increases) with the potential steepness V 0 .

C. Parametrically driven trap
Next, we address the condensate stability in the parametrically modulated trap. To get information about the onset of the parametric instability on a general basis, we write Eq. (6) in terms of two coupled first-order linear differential equations. By using the Floquet theorem [53], we determine whether the solutions for a given set of parameters are stable or not (see appendix for more details). In Fig. 3 driving amplitude is sufficient to destabilize the condensate entirely. A finite particle interaction may cause a transition from a stable to an unstable behavior. Hence, atom-atom interaction may be explored to stabilize a condensate in a parametrically modulated trap by modifying the resonance condition. Moreover, the onset of the parametric instability can be used to measure the interaction strength.

IV. STRONGLY INTERACTING ATOMS
Next, we determine the quantum many-body dynamics of a strongly interacting gas in a lattice and in a parametrically driven trapping potential. To this end, we use the numerically exact time-dependent TEBD method, which is a variant of the time-dependent density matrix renormalization group [54][55][56]. We determine the numerically exact transient dynamics and are able to investigate the onset of Mott physics and its interplay with the parametric resonance.

A. Static trap
First, we consider the static case δV = 0 and calculate the ground state of the Hamiltonian (1). Then, the condensate width is extracted as the full width at half maximum (FWHM) of the distribution of the local occupation number n [57]. We use a lattice with M = 32 sites filled with N = 16 bosons in a trap with V 0 = 0.0922J . In Fig. 4, we show the condensate width in dependence of U calculated numerically exactly in comparison with the mean-field result α 0 . The inset depicts the corresponding distribution n . A perfect agreement is found for U = 0, where both approaches yield coinciding Gaussians. Deviations between the two results increase for growing U , as expected. For U = 7J , a plateau-like Mott region with a nearly integer occupation number starts to appear in the numerical result; see inset of Fig. 4. This wedding-cake-like structure cannot be reproduced by the variational mean-field approach. The reasons for the deviations are twofold. Beyond finite-size effects in the numerically exact result, for increasing U , the condensate starts to locally form a Mott insulating state. The quantum fluctuations at larger U become more important and lead to a broadening of the population density. This is not taken into account by the GPE mean-field approach, but is captured by the TEBD. In fact, for all U , the condensate width is systematically larger by TEBD than predicted by the mean-field approach. Interesting is the step-like increase of the condensate FWHM at U ≈ 7J which accompanies the formation of the Mott plateau. In fact, this is a signature of the Berezinsky-Kosterlitz-Thouless quantum phase transition which has been shown to occur between U = 7.5J and U = 8J for the n = 1 Mott lobe [58,59].

B. Parametrically driven trap
Having studied the static trap, we next address the transient dynamics of the parametrically driven trap. At initial time t = 0, we prepare the system in the ground state of the Hamiltonian H (0). In Fig. 5(a), we show the time evolution of the total energy E(t) = E(t) − E(0) as a function of time t and driving frequency with E(t) = H (t) . In Fig. 5(b), we show two cuts along the lines = 1.31ω and = 2.34ω . The resonance frequency ω can be estimated according to Eq. (7). Resonant driving in Fig. 5 is expected to occur close to /ω = 2/n. We find clear evidence of the two-photon resonance n = 2 which is manifest in a growth of the total energy after each period. The resonance is slightly shifted to 1.3ω , due to the strong interaction between the atoms which is not taken into account in the variational mean-field description. However, the n = 1 resonance is not observed in the transient dynamics, but is expected to occur at longer times. For off-resonant driving, the energy oscillates with the frequency and with small variations in the amplitudes. 043604-4 The result is shown in Fig. 6. Note that for smaller values of , the time needed to complete the mth period is longer than for larger driving frequencies. Therefore, data points for smaller are averaged over fewer periods. Yet, this does not change the result substantially. A rather broad peak shows up around the second resonance ω which is slightly shifted to higher frequencies. The slight oscillations at the flank as well as the rise of the first resonance around 2ω can be expected to vanish when longer times are taken into account for the averaging.  As can be seen in the inset of Fig. 4, when the interaction is below or still in the vicinity of the Berezinsky-Kosterlitz-Thouless phase transition, the mean-field result for particle density agrees well with that of TEBD. However, for stronger interaction, this is in general no longer true and the features of the phase transition are not caught by a mean-field ansatz. In a global harmonic potential, zones of different quantum phases may coexist, and sites that locally realize a Mott insulator state hinder the expansion and contraction of the condensate. This can be seen in Fig. 7 where we show a comparison of the time evolution of the total time-dependent energy E(t) for the cases of noninteracting bosons [ Fig. 7(b)], interacting bosons with U = J [ Fig. 7(c)], and strongly interacting bosons with U = 10J [ Fig. 7(d)]. While the ground state (with the static potential curvature V 0 = 0.0922J ) of the cases (b) and (c) occupies a superfluid state on all sites, the central about 14 sites of the ground state in case (d) occupy a Mott state. According to the mean-field result of Eq. (7), the resonance frequencies can be estimated for the three cases to be ω 1.18J [ Fig. 7(b)], ω 1.06J [ Fig. 7(c)], and ω 1.05J [ Fig. 7(d)]. We find that the resonance frequencies for the cases (b) and (c), determined by the mean-field argument, fit well to the transient behavior that we observe. This is also reflected in the time-averaged energy shown in Fig. 7(a). The slight discrepancies can be traced back to the choice of a rather steep potential such that α 0 1. Then, the mean-field ansatz becomes not very reliable and further corrections become necessary. However, for the case of strong interaction shown 043604-5 in Fig. 7(d), no significant energy absorption occurs in the vicinity of the resonance frequency and the energy returns to almost its initial value after each period. In contrast, in the cases (b) and (c) a clear growth can be seen.
Furthermore, we show in Fig. 8(a) the time evolution of the time-dependent width δ (t) = n t | − 0 | of the particle density for different values of the driving frequency. We choose a parametric drive of the form V (t) = V 0 − δV cos( t), such that V (0) V (t) is satisfied for all times, and a rather large driving strength δV = 0.5V 0 in order to enhance the impact of the parametric resonance. For each frequency, we choose the initial state to be the ground state of the initial Hamiltonian H (0). The transient behavior shows a growth of δ in the vicinity of 2ω and ω . To illustrate that the initial state locally occupies a Mott state, we show in Fig. 8(b) the density profile of the initial state n and the local compressibility κ = n 2 − n 2 . In the central about 25 sites, the system initially occupies a Mott insulating state, which is characterized by a reduced compressibility. The particle motion on these sites is strongly suppressed because an energy of the order of U is needed to move particles inside this region. In that sense, the Mott region acts as a local barrier that suppresses the contraction process during the time evolution. Consequently, particle currents are only observed inside the superfluid regions during the time evolution.

V. EFFECTIVE SITE-DEPENDENT HOPPING
The parametric driving of the global trap can also be used to create a spatially varying hopping strength. By a time-dependent unitary transformation the time dependence of the potential can be converted to a time-and site-dependent hopping amplitude. With this, the bosonic annihilation operator transforms according to The exponential can be absorbed into the hopping amplitude to define a site-and time-dependent hopping amplitude Making use of the Jacobi-Anger identity to expand the exponential in terms of the mth Bessel function of the first kind J m (x), we obtain Thus, for a large enough driving frequency , the time average yields an effective local hopping J eff = J J 0 (δV [2( − 0 ) + 1]/ ), where the spatial dependence is imprinted by the Bessel function J 0 (x).

VI. CONCLUSIONS
We have studied a driven one-dimensional Bose-Hubbard model with a periodically modulated harmonic trap. It describes an interacting gas of bosonic atoms in a lattice which is placed in a parametrically modulated trapping potential. This model allows for the investigation of the interplay of strong atomic interactions in the Mott state of the lattice and the external parametric drive. We have analyzed the parametric resonance condition first in the mean-field regime of weak atomic interactions. We find that also in the presence of the lattice, the condensate width is governed by the Mathieu equation. The resonance frequency of the condensate is shifted to lower values for increasing interaction. Moreover, the phase diagram of stable and unstable dynamics is inherited from the Mathieu equation. Furthermore, for stronger interaction, the mean-field approach becomes invalid and the numerically exact TEBD technique has to be invoked. For strong interaction, we observe the formation of a local Mott-insulator region in which the movement of atoms also in the presence of the driving becomes suppressed. Finally, we have demonstrated that the global parametric modulation yields site-dependent hopping amplitudes which can be controlled. Interestingly, locating the onset of the instability allows, in principle, to determine the atom interaction strength. Thus, dynamically probing a quantum many-body system with a periodic modulation of the harmonic confinement provides a diagnostic tool, which warrants an experimental realization in the realm of ultracold Bose gases.

ACKNOWLEDGMENT
This work was supported by the German Research Foundation (DFG), the DFG Collaborative Research Centers SFB 925 and SFB/TR185, and by the Academy of Finland (Contract No. 275245).

APPENDIX: STABILITY ANALYSIS
In this appendix, we provide more details on the stability analysis of Eq. (5). To study the appearance of the instability in the parametrically modulated trap potential, we make use of the well-known Floquet theorem. The theorem is directly applicable to systems whose dynamics are described by a set of d coupled linear differential equationṡ where A(t + T ) = A(t) is a d-dimensional matrix with periodicity T . The main idea is to evolve the fundamental solution (t) over a single period T in time. Whether the system is stable [i.e., all solutions of (A1) are stable] or unstable [i.e., there exists at least one solution of (A1) that is unstable] can be characterized by the eigenvalues λ i=1,...,d of the monodromy matrix B = −1 (0) (T ). The eigenvalues λ i are directly connected to the Floquet exponents ν i via λ i = e T ν i . When all eigenvalues satisfy |λ i | 1, the real part of each Floquet exponent is less than or equal to zero and the solution is stable. In turn, in the case that |λ i 0 | > 1 for any i 0 , there exists an unstable solution of Eq. (A1) with Re(ν i 0 ) > 0.
To apply the Floquet theorem to Eq. (5), we define the variable τ = sin t, which allows us to rewrite Eq. (5) in terms of a linear differential equation like Eq. (A1) with the vector x = (δα,δ α,τ,τ ) t and the matrix with the period T = 2π/ . To construct the fundamental solution, we have to calculate the time evolution x j =1,...,4 (t) of four linearly independent initial vectors x j (0). Then, the fundamental solution is given by (t) = (x 1 (t) x 2 (t) x 3 (t) x 4 (t)).
In practice, we choose the initial conditions such that (0) = 1 is the identity matrix. From this, we numerically calculate (t) and analyze the spectrum of the monodromy matrix B. In the eigenvalue spectrum of B, we always find two eigenvalue λ j 1 = λ j 2 = 1 which corresponds to the solutions of the differential equationτ = − 2 τ . The remaining two eigenvalues are used to characterize the stability behavior of Eq. (5).