Dynamics-induced freezing of strongly correlated ultracold bosons

The theoretical study of non-equilibrium dynamics in closed quantum systems has seen great progress in recent years [1] mainly due to the possibility of realization of such dynamics using ultracold atoms in optical lattices [2, 3]. For bosonic atoms, such systems are well described by the Bose-Hubbard model with on-site interaction strength U and nearest-neighbor hopping amplitude J [4, 5]. Several theoretical studies have been carried out on the quench and ramp dynamics of this model [6–11]; some of them have also received support from recent experiments [3]. In contrast, studies on periodically driven closed quantum systems have been undertaken in the past mainly on driven two-level systems [12, 13] or on weakly interacting or integrable many-body systems which can be modeled by them [14, 15]. Among these, ref. [15] has predicted freezing of the time-averaged value of the order parameter (magnetization) of a periodically driven one-dimensional (1D) Ising or XY model, when the temporal average is performed over several drive cycles, at specific drive frequencies. Such a freezing occurs in the high-frequency regime and exhibits non-monotonic dependence on the drive frequency. However, to the best of our knowledge, the phenomenon of dynamics-induced freezing has never been demonstrated for dynamics involving a single drive cycle and/or for non-integrable quantum systems. Recent studies of periodic dynamics of the Bose-Hubbard model have not addressed this issue [16, 17].

In this work, we demonstrate that the periodically driven Bose-Hubbard model may exhibit dynamical freezing of the boson wave function by observing that after one cycle the wave function comes back to its original value ψ(T) = ψ(0) for specific values of the drive frequency ω = 2π/T. Our driving protocol constitutes a time-dependent hopping amplitude of the bosons J(t) = J₀ + δJcos(ωt) with J₀ and δJ chosen such that the drive takes the system from a superfluid (SF) (J = J₀ + δJ) to the Mott insulator (MI) state (J = J₀ − δJ) and back through the tip of the Mott lobe where μ = μ_tip. We provide a semi-analytical understanding of the freezing phenomenon using mean-field theory and compute the defect formation probability P = 1 − F (where F = |〈ψ(t = 0)|ψ(t = T)〉|² is the wave function overlap), the superfluid order parameter Δ(T) = 〈ψ(T)|b|ψ(T)〉 (where b denotes the boson annihilation operator), and the residual energy Q = E(t = T) − E_G (where E(t = T) is the energy of the system at the end of the drive cycle and E_G is the initial ground-state energy) as a function of ω. We also show, via inclusion of quantum fluctuation by a projection operator approach [11] and numerical mean-field study of a trapped boson system that the freezing phenomenon is qualitatively robust against quantum fluctuations and the presence of a trap. We also show that such a phenomenon persists for multiple drive cycles. We note that such a freezing behavior has two novel characteristics which distinguishes it from its counterpart in ref. [15]. First it always occurs at low frequencies, i.e., ℏω/U ≪ 1 for all ω for which the freezing occurs. Second, it occurs for a single cycle of the drive. To the best of our knowledge, the dynamics-induced freezing phenomenon at such low drive frequencies and non-integrable models has not been studied in the context of closed quantum systems; our work therefore constitutes a significant advance in our understanding of periodic dynamics of closed non-integrable quantum systems.

The Hamiltonian describing a system of ultracold bosonic atoms confined by a trap and in an optical lattice is given by

$$\begin{eqnarray} {\mathcal H} &=& \sum_{\langle {\bf r},{\bf r'}\rangle} -J b_{{\bf r}}^{\dagger} b_{{\bf r'}} + \sum_{{\bf r}} \left[-\mu_{\bf r} {\hat n}_{{\bf r}} + \frac{U}{2} {\hat n}_{{\bf r}}({\hat n}_{{\bf r}}-1)\right], \end{eqnarray}\noindent \tag{ 1 }$$

where μ_r denotes the chemical potential at site r, r' denotes one of the z nearest-neighboring sites of r, and ${\hat n}_{\mathbf { r}}= b_{\bf r}^{\dagger } b_{\bf r}$. In the absence of a trap, μ_r = μ for all sites and for zJ ≪ U, the ground state of the model is a MI state with ${\bar n}$ bosons per site with ${\bar n}=1$ for 0 ⩽ μ/U ⩽ 1. For zJ ≫ U, the bosons are delocalized and the system, for d ⩾ 2, is in a SF state. In between, at J = J_c, the system undergoes a SF-MI transition. The equilibrium phase diagram of the model constitutes the well-known Mott lobe structure [4, 5].

To obtain an semi-analytic insight to the freezing phenomenon, we first analyze the periodically driven Bose-Hubbard model in the absence of a trap and within mean-field approximation. The time-dependent mean-field Hamiltonian is given by

$$\begin{eqnarray} {\mathcal H}_{{\rm mf}} &=& \sum_{{\bf r}} \left[-\mu {\hat n}_{{\bf r}}+ \frac{U}{2} {\hat n}_{{\bf r}}({\hat n}_{{\bf r}}-1)\right] +\left(\Delta'_r(t) b_r^{\dagger} + {\rm h.c.} \right),\nonumber\\ \end{eqnarray}\noindent \tag{ 2 }$$

where $\Delta '_r(t) = -J(t) \sum _{\langle r'\rangle } \langle b_{r'} \rangle$ and Δ'₀ = Δ'(t = 0). Within homogeneous mean-field theory, the Gutzwiller wave function for the bosons reads $|\psi ({\mathbf { r}},t)\rangle _{{\rm mf}}= \prod _{{\bf r}} \sum _n f_n(t) |n\rangle$ [18]. The Schrödinger equation $i \partial _t |\psi ({\mathbf { r}},t)\rangle _{{\rm mf}} = {\mathcal H}_{{\rm mf}}(t) |\psi ({\bf r},t)\rangle _{{\rm mf}}$ yields the time-dependent mean-field equations for f_n(t) ≡ f_n:

$$\begin{eqnarray} &&(i \partial_t -E_n)f_n = {\tilde \Delta}(t) \sqrt{n} f_{n-1} + {\tilde \Delta}^{\ast}(t) \sqrt{n+1} f_{n+1}, \end{eqnarray} \tag{ 3 }$$

where ${\tilde \Delta }(t) = -z J(t) \sum _n \sqrt {n} f_{n-1}^{\ast } f_n$, and E_n = −μn + U(n − 1)n/2. In what follows, we shall choose J₀ and δJ such that the ground state of H_mf with J = J₀ + δJ is a SF state close to the QCP so that f_n(t = 0) ≃ 0 for n ⩾ 3 and f₁(t = 0) ≫ f₀(t = 0) ,f₂(t = 0).

Our strategy for solving these equations will be the following. First, we shall solve the mean-field equations numerically keeping states with n ⩽ 4. Such numerics yields the result that in the MI and SF states close to the QCP all f_n(t) for n ⩾ 3 remains small for all t during the dynamics. We shall then use this fact to develop semi-analytic understanding of the freezing phenomenon using equations for f₀, f₁ and f₂ (see footnote ¹).

The equations for f₀, f₁ and f₂ then read (suppressing the time dependence of f_n(t) for clarity)

$$\begin{eqnarray} i \partial_t f_0 &=& - zJ(t) [ |f_1|^2 f_0 + \sqrt{2} f_2^{\ast} f_1^2] , \nonumber \\ i \partial_t f_2 &=& E_2 f_2 - zJ(t) [2 |f_1|^2 f_2 + \sqrt{2} f_0^{\ast} f_1^2], \\ i \partial_t f_1 &=& E_1 f_1 - zJ(t) [(2 |f_2|^2 +|f_0|^2) f_1 + 2 \sqrt{2} f_1^{\ast} f_2 f_0 ]. \nonumber \end{eqnarray} \tag{ 4 }$$

Multiplying each of the equations for f_n in eqs. (4) by f*_n and subtracting the equations thus obtained from their complex conjugates, it is easy to see that |f_n|, for n ⩽ 2, obeys the relation

$$\begin{eqnarray} \partial_t |f_0|^2 &=& \partial_t |f_2|^2 = - \partial_t |f_1|^2/2. \end{eqnarray} \tag{ 5 }$$

Parameterizing $f_n=r_n(t) \exp [i \phi _n (t)]$, with the choice that ϕ_n(0) = 0 and r_n(0) = f_n(0), one can write

$$\begin{eqnarray} r_{2[0]}^2(t) &=& -(r_1^2(t) -1)/2 +[-] \eta, \end{eqnarray} \tag{ 6 }$$

where η is a time-independent parameter whose value is fixed by the initial values r_n (see footnote ²) Note that η represents the magnitude of the particle-hole asymmetry since r₀ = r₂ for η = 0. Substituting eq. (6) in eq. (4), we get

$$\begin{eqnarray} \partial_t r_1 &=& -\sqrt{2} z J(t) \sin(\phi_s) r_1 g_0(r_1), \nonumber\\ \partial_t \phi_s &=& -U +z J(t)\left[ g_1(r_1) - g_2(r_1) \cos(\phi_s) \right],\\ \partial_t \phi_d &=& -U+2\mu + z J(t) r_1^2 \left[ 1 - 4 \sqrt{2} \eta \cos(\phi_s)/g_0(r_1) \right], \nonumber \end{eqnarray} \tag{ 7 }$$

where we have suppressed the time dependence of r₁ and ϕ_s(d) for clarity, ϕ_s = ϕ₀ + ϕ₂ − 2ϕ₁ and ϕ_d = ϕ₂ − ϕ₀ are the sum and differences of the relative phases of the Gutzwiller wave function, and the functions g_i(r₁) are given by

$$\begin{eqnarray} g_0(r_1) &=& \sqrt{(1-r_1^2)^2-4\eta^2}, \quad g_1(r_1)= 6 r_1^2-3 -2 \eta, \nonumber\\ g_2(r_1) &=& 2 \sqrt{2} \left[ r_1^2(r_1^2-1)/g_0(r_1)+g_0(r_1) \right]. \end{eqnarray} \tag{ 8 }$$

We note that the first two of the equations in eqs. (7) are coupled equations describing the evolution of r₁ and ϕ_s, while the third describes the evolution of ϕ_d in terms of r₁ and ϕ_s. Furthermore, using a scaled variable t' = ωt/(2π), we find that the relation between r₁ and ϕ_s can be written as

$$\begin{eqnarray} {\rm d}r_1/{\rm d}\phi_s &=& \frac{-\sqrt{2} \sin(\phi_s) r_1 g_0(r_1)}{[g_1(r_1)-g_2(r_1)\cos(\phi_s)]- U/zJ(t')}. \end{eqnarray} \tag{ 9 }$$

This allows us to symbolically write ϕ_s = ξ(r₁,t'), where ξ is an unknown function, and thus establish an ω-independent relation between r and ϕ_s for any fixed t'.

Equations (6), (7) and (8) constitute the central result of this work. They constitute a complete description of the evolution of f₀, f₁, and f₂ in the presence of the periodic drive and provide an understanding of the freezing phenomenon as follows. First, we find that a numerical solution of eq. (7), together with eq. (6), allows us to obtain r_n, ϕ_s and ϕ_d as functions of time. The plots of 1 − r₁(t) and ϕ_s(t) as functions of t for ω/J_c = 0.52 is shown in panel (a) of fig. 1. We find that r₁ changes appreciably when J(t) is close to J_c; however, r₁(T) ≃ r₁(0) at the end of the evolution. Note that this also implies, via eq. (6), that r₂₍₀₎(0) ≃ r₂₍₀₎(T). Panel (b) of fig. 1 shows that the relation r₁(T) ≃ r₁(0) holds for a significant range ω/J_c ⩽ 0.8. Second, we note that, ϕ_s(t) undergoes rapid oscillation when J(t) ⩽ J_c; however, it also comes back close to its initial value at the end of the drive: ϕ_s(T) ≃ ϕ_s(0). Since ϕ_s and r₁ satisfies an ω-independent relation, ϕ_s = ξ(r₁,t'), we infer that ϕ_s must remain close to its initial value for the same range of ω for which r₁(T) ≃ r₁(0); this is verified numerically in panel (b) of fig. 1. Finally, we note, from the panels (c) and (d) of fig. 1, that ϕ_d(T) is a monotonic function of ω. Thus, we may define ω = ω*_m for which ϕ_d(T) ≃ 4πm (m being an integer). Together with the fact that r₁(T) ≃ r₁(0) and ϕ_s(T) ≃ ϕ_s(0) = 0, we find that at ω = ω*_m, both the relative phases satisfy ϕ₂ − ϕ₁ = −(ϕ₀ − ϕ₁) ≃ 2πm leading to |ψ_mf(T)〉 ≃ |ψ_mf(0)〉 up to a global phase. This constitutes the dynamics freezing of |ψ〉_mf.

Zoom In Zoom Out Reset image size

To obtain an accurate estimate of the degree of freezing and to check the validity of the discussion above, we compute the defect density P = 1 − F = 1 − |〈ψ_mf(T)|ψ_mf(0)〉|² via direct numerical solution of eq. (3) keeping n ⩽ 4 states per site. The plot of P as a function of ω clearly shows that P → 0 at ω = ω*_m. A plot of Log₁₀P vs. ω near ω*₉/J_c ≃ 0.47, shown in panel (a) of fig. 2, reveals that Log₁₀P ∼ − 6 indicating that the overlap, up to a global phase, is exact within our numerical accuracy. We have checked for all ω*_m ⩽ 0.8J_c, Log₁₀P < −4 which indicates a near perfect freezing. We also compute the residual energy Q(T) and the SF order parameter

$$\begin{eqnarray} \Delta &=& r_1 e^{i\phi_d/2} \left(r_0 e^{-i \phi_s/2} + \sqrt{2} r_2 e^{i\phi_s/2}\right), \end{eqnarray} \tag{ 10 }$$

at t = T as a function of ω. We find from eq. (10) that |Δ| is independent of ϕ_d. Thus |Δ(T)|/|Δ₀| and Q/U (which can also be shown to be independent of ϕ_d) remain close to unity and zero respectively over the entire range of ω/J for which r₁ and ϕ_s remain close to their initial values as shown in panels (c) and (d) of fig. 2. Such a behavior distinguishes these quantities from P which depends on ϕ_d and hence vanishes at discrete ω_m. Finally, we find that for all values of J₀ shown in panel (b) of fig. 2, there is an appreciable range of ω/J_c within which the freezing phenomenon occurs and that ω*_m decreases monotonically as a function of J₀ over this range.

Zoom In Zoom Out Reset image size

A physical picture of the freezing phenomenon described above can be obtained as follows. The dynamics of the bosons in the MI and the SF phases near the transition constitutes change in both the amplitudes and the relative phases of the different components of the boson wave function. We find that the change in amplitudes after a complete cycle of the drive (at t = T) is close to zero over a range of ω ⩽ J_c. Our treatment of the mean-field equations above brings out the crucial fact that this also means that the change in the sum of the relative phases, ϕ_s(T), will be small after a complete cycle in this regime since r₁ and ϕ_s are related by eq. (9). In contrast, the difference between the relative phases, ϕ_d(T), evolves with the drive frequency ω. Thus, at special values of ω = ω*_m, where ϕ_d(T) = 4πm, the mean-field boson wave function comes back to itself after a complete drive cycle leading to freezing.

Next, we study the effect of quantum fluctuations on the freezing phenomenon³. We incorporate such fluctuations by using a projection operator method developed in ref. [11] which provides an accurate treatment of dynamics with fluctuations for zJ(t)/U ≪ 1 (typical value of J_c/U are ∼0.034 for d = 3 and ∼0.061 for d = 2) as seen for linear dynamics in where it provides a qualitative match with experimental results [3]. The idea behind this approach, as detailed in ref. [11], is to introduce a projection operator P_ℓ = |n₀〉〈n₀|_r × |n₀〉〈n₀|_r' which lives on the link ℓ between the neighboring sites r and r' of the lattice. Using P_ℓ, one can write the boson hopping term as

$$\begin{eqnarray} T' &=& \sum_{\langle {\bf r} {\bf r'} \rangle} -J(t) b_{{\bf r}}^{\dagger} b_{\bf r'} = \sum_{\ell} T'_{\ell} \nonumber\\ &=& \sum_{\ell} [ (P_{\ell} T'_{\ell} + T'_{\ell} P_{\ell}) + P_{\ell}^{\perp} T'_{\ell} P_{\ell}^{\perp}], \end{eqnarray} \tag{ 11 }$$

where P^⊥_ℓ = (1 − P_ℓ). In the strong-coupling regime where zJ(t)/U ≪ 1, the term T^'0_ℓ[J] = (P_ℓT'_ℓ + T'_ℓP_ℓ), at any instant, represents hopping processes which takes the system out of the instantaneous low-energy manifold. Thus, one can devise a time-dependent canonical transformation via an operator S given by

$$\begin{eqnarray} &&S[J(t)]= -\frac{i}{U} \sum_{\ell} [P_{\ell},T'_{\ell}], \end{eqnarray} \tag{ 12 }$$

which eliminates T^'0_ℓ[J(t)] up to first order in J(t)/U and leads to the effective instantaneous time-dependent Hamiltonian $H^{\ast }=\exp (-iS[J(t)]) {\mathcal H} \exp (iS[J(t)])$. Such a canonical transformation is equivalent to a transformation on the system wave function |ψ〉: $|\psi '\rangle = \exp (-i S[J(t)]) |\psi \rangle$. We note that |ψ〉 and |ψ'〉 coincides for J = 0 which leads us to the natural choice $|\psi '\rangle = \prod _{\mathbf { r}} \sum _n f_{\bf r}^n(t) |n_{\bf r}\rangle$. Note that |ψ〉 is not of Gutzwiller form; it involves spatial correlation due to $\exp (iS[J(t)])$ factor. The instantaneous energy of the system is given by

$$\begin{eqnarray} E[\{f_{\bf r}^n(t)\}] &=& \langle \psi|{\mathcal H} | \psi \rangle = \langle \psi'|H^{\ast}| \psi'\rangle + {\rm O}(J(t)^3/U^3)\nonumber\\ \end{eqnarray} \tag{ 13 }$$

and includes O(J²/U²) quantum fluctuation corrections. As shown in ref. [11], this formalism allows one to describe the dynamics of the bosons by solving for the Schrödinger equation for |ψ'〉:

$$\begin{eqnarray} &&(i \hbar \partial_t + \partial S[J(t)]/\partial t) |\psi'\rangle = H^{\ast}[J(t)] |\psi' \rangle. \end{eqnarray} \tag{ 14 }$$

Using the expression of |ψ'〉 and E[{f_n}], one can convert eq. (14) to a set of equations for {fⁿ_r(t)} [11]. Defining $\varphi _n = \sqrt {n+1} f_n^{\ast } f_{n+1}$, one gets

$$\begin{eqnarray} i \hbar \partial_t f_n &=& \delta E[\{f_n(t)\};J(t)]/ \delta f_n^{\ast} + \frac{i z \hbar}{U} \frac{\partial J(t)}{\partial t}\nonumber\\ && \times \Big( \sqrt{n} f_{n-1}\Big[\delta_{n\bar{n}}\varphi_{\bar{n}}- \delta_{n,\bar{n}+1}\varphi_{\bar{n}-1}\Big] \nonumber\\ && + \sqrt{n+1} f_{n+1} \Big[ \delta_{n\bar{n}}\varphi_{\bar{n}-1}^*- \delta_{n,\bar{n}-1} \varphi_{\bar{n}}^* \Big] \Big).\quad\quad \end{eqnarray} \tag{ 15 }$$

A numerical solution of eq. (15) yields f_n(t) and hence |ψ'〉 using which one can compute $|\psi (t)\rangle = \exp [iS] |\psi '(t)\rangle$ perturbatively to O[J(t)²/U²]. Similarly, the expectation value of any operator O at any instant t can be calculated in terms of |ψ'(t)〉:

$$\begin{eqnarray} \langle O\rangle &=& \langle \psi'(t)|e^{-iS} O e^{iS}|\psi'(t)\rangle \nonumber &=& \langle \psi'(t)|O| \psi'(t)\rangle - \langle \psi'(t)|[iS,O]| \psi'(t)\rangle + \ldots,\quad\quad \end{eqnarray} \tag{ 16 }$$

where the ellipsis indicate higher-order terms in J(t)/U. Note that the second term in the expression originates from quantum fluctuation and modifies mean-field result (first term). Using the above-mentioned procedure detailed in ref. [11], we compute P(T), Q(T) and |Δ(T)| as shown in fig. 3. We find that key effects of the quantum fluctuations is to change numerical values of ω*_m and the precise range of ω over which freezing occurs; however the mean-field results hold qualitatively in the sense that P → 0 for several ω*_m with Log₁₀P ⩽ − 4 for all ω*_m. Further ω*_m also decreases monotonically with J₀ as shown in panel (b) of fig. 3 for ω* ≃ 0.6J_c.

Zoom In Zoom Out Reset image size

Next, we consider the effect of a harmonic trap on the freezing phenomenon. For this part, we numerically solve eq. (3) for d = 2 with μ_r = μ₀ + 0.01U[(r_x − 1/2)² + (r_y − 1/2)²], for N₀ = 576 sites (linear dimension 24) and with fixed total particle number N₀. We choose the trap parameters so that the ground state of the bosons in the center of the trap at t = T/2 is MI phase with ${\bar n}=1$. The evolution of the density profile of the bosons is shown in panel (a) of fig. 4 for t = 0 (left) and t = T/2 (right). Panel (b) of fig. 4 indicates value of ϕ_d(T) as a function of the position of the bosons in the trap along the line y = 0. The plot shows that for all ω ⩽ 0.3J_c, ϕ_d evolves coherently with negligible spatial variation. The plots for ϕ_s and r₁ are similar in nature; thus, we expect the boson evolution to have the same qualitative properties as that found within a homogeneous mean-field approach. A plot of P (|Δ(T)|/|Δ₀|) as a function of Log₁₀(ω/J_c) in the panels (c), (d) of fig. 4 confirms this expectation. We find that the main effect of the trap is to push the freezing phenomenon to lower frequencies leaving its qualitative nature unchanged. The largest freezing frequency occurs at ≃0.2J_c which is large compared to frequencies ≃0.05J_c where momentum-conserving boson pair production at finite momenta, which is not captured within mean-field theory, is expected to become significant [19]. Figure 4 also demonstrates that the freezing phenomenon disappears at higher drive frequencies where the trapped bosons do not evolve coherently leading to spatial variation of ϕ_d.

Zoom In Zoom Out Reset image size

Finally, we note that the freezing phenomenon persists for multiple periods of the drive. To show this explicitly, we revert to single-site mean-field theory and plot P around a freezing frequency ω* = 0.47J_c for several cycles of the drive in fig. 5. The plot shows that P(ω*) increases with number of drive cycles; however, even after t = 10T, its value remains small (⩽10⁻⁴) which indicates that the freezing persists to a good degree of accuracy for longer drive times.

Zoom In Zoom Out Reset image size

For experimental verifications of our work, we suggest interference of two bosonic condensates in the presence of an optical lattice, near the QCP which are separated after creation by a double-well potential and allowed to evolve separately for a fixed holdout time. It is well known that recombination of such separated condensates can act as a readout scheme for their relative phases [20]. We propose such a readout when one of the condensates is driven periodically with a frequency ω during the holdout for a single period T = 2π/ω. Our specific prediction is that the relative phase measured for such a drive with ω = ω*_m is going to match the phase without any drive indicating dynamics-induced freezing. For all such experiments one needs to estimate an optimal temperature T₀ at which they can be carried out. The typical value of U deep inside the Mott phase is ≃2 KHz = 200 nK leading to a melting temperature of T_m ≃ 0.2U = 40 nK for d = 3. The SF phase near the Mott tip has a coherence temperature of T_c ≃ zJ_c ≃ 35 nK [21]. Thus, a temperature of a few nano-kelvins (T₀ ≪ T_m,T_c), which is currently within the experimental reach, would be ideal for testing our prediction.

In conclusion, we have demonstrated that periodic dynamics of the ultracold bosons described by the Bose-Hubbard model leads to dynamics-induced freezing of the boson wave function at specific drive frequencies which are determined by the condition ϕ_d(T) = 4πm. The freezing phenomenon manifests itself at discrete drive frequencies ω*_m ⩽ J_c via the presence of dips in the defect density and can be detected by suitable interference experiments. We have shown that this effect, which can be qualitatively understood using mean-field theory, is robust against quantum fluctuations, as incorporated using a projection operator technique, and the presence of a trap. We note that designing assisted driving protocols for near-adiabatic drives in interacting many-body systems has been the subject of a plethora of theoretical and experimental studies lately [22]. Our work provides a rare example of near-adiabatic driving without assistance⁴ in an experimentally realizable non-integrable quantum many-body system.

SM and KS thanks K. Ray for several stimulating discussions. KS thanks DST for support through grant SR/S2/CMP-001/2009. DP acknowledges support from the Lee A. DuBridge fellowship.