Dynamically generated correlations in a trapped bosonic gas via frequency quenches

TL;DR

This paper investigates how stochastic frequency quenches in a trapped bosonic gas induce nontrivial correlations and a nonequilibrium steady state, revealing universal scaling laws and complex statistical structures.

Contribution

It introduces a novel protocol of repeated quantum quenches that drives a bosonic gas into a unique NESS with emergent long-range correlations and provides analytical and numerical characterization of this state.

Findings

Discovery of a CIID structure in the joint probability density function

Identification of long-range dynamical correlations in the NESS

Universal scaling laws in order and gap statistics

Abstract

We study a system of $N$ noninteracting bosons in a harmonic trap subjected to repeated quantum quenches, where the trap frequency is switched from one value to another after a random time duration drawn from an exponential distribution. Each cycle contains two steps: (i) changing the trap frequency to enable unitary evolution under a Hamiltonian, and (ii) reapplying the original trap at stochastic times to cool the gas back to its initial state. This protocol effectively makes it an open quantum system and drives it into a unique nonequilibrium steady state (NESS). We analytically and numerically characterize the NESS, uncovering a conditionally independent and identically distributed (CIID) structure in the joint probability density function (JPDF) of the positions. The JPDF in the CIID structure is a product of Gaussians with a common random variance, which is then averaged with respect to its distribution, making the JPDF non-factorizable, giving rise to long-range emergent dynamical correlations. The average density profile of the gas shows significant deviations from the initial Gaussian shape. We further compute the order and the gap statistics, revealing universal scaling in both bulk and edge regimes. We also analyze the full counting statistics, exposing rich parameter-dependent structure. Our results demonstrate how stochastic quenches can generate nontrivial correlations in quantum many-body systems.

Equations

$$H({\omega }_1)=\sum _{j=1}^N\left(\frac{p_j^2}{2m}+\frac{1}{2}m{\omega }_1^2{x}_j^2\right),$$

$$\psi (x_1,x_2,\dots ,x_N)=\prod _{j=1}^N{\psi }_0(x_j),\text{where}{\psi }_0(x)={\left(\frac{m{\omega }_1}{\hslash \pi }\right)}^{1/4}\exp \left(-\frac{m{\omega }_1{x}^2}{2\hslash }\right)$$

$$Prob.(x_1,x_2,\dots x_N,t)=e^{-rt}\prod _{j=1}^{N}p(x_j,t)+r{\int }_0^{t}d\tau e^{-r\tau }\prod _{j=1}^{N}p(x_j,\tau ).$$

$$K(x,y,\tau )\equiv ⟨x|e^{-i\tau H({\omega }_2)/\hslash }|y⟩=\frac{m{\omega }_2}{2i\pi \hslash \sin ({\omega }_2\tau )}\exp \left(\frac{im{\omega }_2}{2\hslash \sin ({\omega }_2\tau )}\left[(x^2+y^2)\cos ({\omega }_2\tau )-2xy\right]\right).$$

$$p(x,\tau )\equiv p(x|V(\tau ))=\frac{1}{\sqrt{2\pi {\sigma }_1^{2}V(\tau )}}\exp \left(-\frac{x^2}{2{\sigma }_1^{2}V(\tau )}\right),\text{where}{\sigma }_1=\sqrt{\frac{\hslash }{m{\omega }_1}},$$

$$V(\tau )=\frac{1}{4}\left[(1+f^2)+(1-f^2)\cos (2{\omega }_2\tau )\right],\text{with f=ω2ω1. }$$

$$V_1=\frac{1}{2}\min [f^2,1],\text{and }{V}_2=\frac{1}{2}\max [f^2,1].$$

$$\mathrm{Prob}.(x_1,x_2,\dots x_N)={\int }_{V_1}^{V_2}dVh(V)\prod _{j=1}^{N}p(x_j|V),\text{where}h(V)=r{\int }_0^{\infty }d\tau e^{-r\tau }\delta [V-V(\tau )]$$

$$\cos (2{\omega }_2\tau )=\alpha (V)=\frac{4V-f^2-1}{1-f^2},$$

$$h(V)=\frac{q}{2\sinh (q\pi /2)}\frac{\cosh \left[\frac{q}{2}{\cos }^{-1}(-\alpha )\right]}{\sqrt{(V_2-V)(V-V_1)}}\text{with}\alpha (V)=\frac{4V-1-f^2}{1-f^2},$$

$$\rho (x)={\int }_{V_1}^{V_2}dVh(V)\frac{e^{-x^2/(2V)}}{\sqrt{2\pi V}}.$$

$$h(V_{2}-\epsilon)\simeq\frac{A_{f}}{\sqrt{\epsilon}}\quad\text{where}\quad A_{f}=\frac{q\cosh\big{[}(q/2)\,\cos^{-1}(-\alpha(V_{2}))\big{]}}{\sqrt{2|1-f^{2}|}\,\sinh(q\pi/2)}.$$

$$\rho (x)\sim A_f\frac{e^{-x^2/(2V_2)}}{\sqrt{2\pi V_2}}{\int }_0^{\infty }\frac{dϵ}{\sqrt{ϵ}}{e}^{-ϵx^2/(2V_2^2)}=A_f\frac{\sqrt{V_2}}{|x|}{e}^{-x^2/(2V_2)},$$

$$\mathrm{Prob}.(M_k=w)=\frac{w}{{\theta }_{\alpha }^2}h\left(\frac{w^2}{2{\theta }_{\alpha }^2}\right),$$

$$\mathrm{Prob}.(M_k-M_{k+1}=d_k)=\frac{1}{\lambda }F\left(\frac{d_k}{\lambda }\right),F(z)={\int }_{V_1}^{V_2}dVh(V)\frac{e^{-z/\sqrt{V}}}{\sqrt{V}}.$$

$$F(z)\sim \frac{A_f}{\sqrt{V_2}}{e}^{-z/\sqrt{V}}{\int }_0^{\infty }\frac{dϵ}{\sqrt{ϵ}}{e}^{-ϵz/(2V_2^{3/2})}=\sqrt{2\pi }{A}_f\frac{V^{1/4}}{\sqrt{z}}{e}^{-z/\sqrt{V}},$$

$$\mathrm{Prob}.(\kappa ,L)=\frac{\sqrt{\pi }{L}^2}{2[{\mathrm{erf}}^{-1}(\kappa ){]}^3}\exp \left([{\mathrm{erf}}^{-1}(\kappa ){]}^2\right)h\left(\frac{L^2}{2[{\mathrm{erf}}^{-1}(\kappa ){]}^2}\right).$$

Full text

1 Introduction

Quantum gases have emerged as a promising ground for exploring nonequilibrium quantum dynamics [1], with cold atom experiments providing exceptional control over interactions [2, 3], trapping geometries [4, 5, 6], and driving protocols [7, 8]. Among the many ways to drive such systems, quantum quenches, which are abrupt changes in a system parameter, play a central role in probing relaxation, thermalization, correlations, and phase transitions [8, 9, 7, 10, 11, 12, 13]. A particularly promising class of driving protocols, which has been implemented in classical systems, involves stochastic sequences of quenches, in which the time between parameter changes is drawn from a probability distribution [14, 15]. Such protocols are experimentally feasible and can fundamentally alter the long-time dynamics, producing steady states that differ from both conventional thermal equilibrium and those of periodically driven systems. The repeated random quenches typically act as an effective driven-dissipative mechanism, rendering the system an open quantum system with an effectively non-Hermitian description. Even though the microscopic evolution between quenches is unitary, the stochastic timing of the quenches, together with averaging over realizations, leads to a breakdown of purely Hamiltonian dynamics. As a result, the system can settle into a nonequilibrium stationary state (NESS) characterized by dynamically induced correlations even among otherwise independent (noninteracting) bosons.

In this work, we investigate one such scenario: a gas of $N$ noninteracting bosons in a harmonic trap whose frequency is repeatedly switched between two values. The switching events occur after random durations sampled from an exponential distribution. Although the bosons do not interact directly, the simultaneous quench protocol induces effective all-to-all correlations. We would like to remark that the center-shift quench protocol studied in Ref. [16] is typically experimentally more demanding in cold atom setups than a frequency-switching protocol discussed in this work. This is because center-shift protocols can involve the physical transport of a lens or a magnetic coil, making it prone to excess heating due to the transfer of vibration from the transport stage [17]. On the other hand, controlled changes of the trap frequency can be implemented with significantly reduced heating. Trap frequencies can be altered with high precision using optical or magnetic confinement techniques. In a typical trap, frequencies over a wide range (e.g., $10-500$Hz) [13] can be tuned on millisecond timescales. The exponential waiting-time distributions with mean intervals from a few milliseconds to several seconds can be programmed using acousto-optic modulators, making the stochastic driving readily achievable. For ultracold bosonic gases such as 87Rb and 23Na, these parameter ranges lie well within current experimental capabilities, ensuring that the proposed frequency-quench protocol is experimentally feasible.

Remarkably, the stationary many-body joint probability density function (JPDF) exhibits a conditionally independent and identically distributed (CIID) structure, a property known from certain exactly solvable classical models with stochastic resetting [18, 19, 20, 21, 22, 16, 14]. This structure, where a set of random variables are independent and identically distributed, given a certain condition, makes it possible to compute a wide range of observables exactly, including density profiles, correlation functions, gap and extreme value statistics, and full counting statistics. Many of these quantities are directly measurable in ultracold atom experiments via absorption imaging [23, 24, 25, 26] or quantum gas microscopy [27, 28, 29], suggesting that stochastic frequency quenches offer a realistic and analytically tractable route to exploring strongly correlated nonequilibrium quantum states without direct interparticle interactions.

We briefly summarize our main findings: (i) We find a simple quenching protocol that generates strong, all-to-all correlations in a bosonic gas, with the JPDF having a CIID structure in the steady state. (ii) We find exactly the average density profile and the statistics of the gap between successive particles. We show that the tail of the average density profile remains Gaussian, and that the tail of the gap distribution remains exponential. (iii) We show that there is a strict upper and lower bound for the number of particles that can be contained within any finite volume in the large-$N$ limit.

The paper is organized as follows. In Sec. 2, we discuss the details of our setup and the protocol, along with details of numerical simulations employed. We obtain an interesting $N$-particle JPDF that turns out to be of CIID form. In Sec. 3, we exploit the rich CIID structure to derive various observables. Perfect agreement is demonstrated between our analytical results and direct numerical simulations. In Sec. 4, we summarize our results along with an outlook.

2 Setup and Protocol

Our setup consists of a gas of noninteracting bosons trapped in a harmonic potential. The protocol we employ is as follows:

The gas is allowed to cool to the ground state of the system, described by the Hamiltonian

[TABLE]

where $\omega_{1}$ is the trap frequency, $m$ is the mass of the particles, and $N$ are the number of particles. This ground state is described by the product state

[TABLE]

is the ground state of the quantum harmonic oscillator. 2. 2.

The frequency is then quenched to $\omega_{2}$, and the wavefunction evolves unitarily under the Hamiltonian $H(\omega_{2})$, with $H(\omega)$ given in Eq. (1). This unitary evolution occurs until a stochastic time $\tau$, picked from the exponential distribution $r\,e^{-r\tau}$, where $r$ is the resetting rate.

At the end of the duration $\tau$, the system is once again cooled to the ground state of $H(\omega_{1})$. Thus, the system is reset to its initial state in Eq. (2). The entire procedure 1-2 is then repeated. Step 1 of our protocol is assumed to occur instantaneously. In practice, relaxation to the ground state takes a finite time, but we emulate this by disregarding the elapsed time (and making no measurements) until the system reaches the ground state with frequency $\omega_{1}$, as done in recent colloidal experiments [15, 30]. In the long-time limit, the system approaches a unique NESS, described by a CIID structure. We discuss this structure below.

Under this protocol, the JPDF of the positions of the particles at any time $t$ can be written using the renewal equation [31]

[TABLE]

Here $\tau$ is the random time from the last resetting event, and $p(x_{j},\tau)=|\psi(x_{j},\tau)|^{2}$, with $\psi(x,\tau)=\int_{-\infty}^{\infty}\,\langle x|e^{-i\tau{H}\,(\omega_{2})/\hbar}|y\rangle\,\psi_{0}(y)\,dy$. We are interested in the stationary state observables, in the limit $t\rightarrow\infty$. In this limit, the first term of the renewal equation, Eq. (3) vanishes and one finds the emergence of a CIID structure, $\mathrm{Prob}.(x_{1},x_{2},\ldots x_{N})=r\int_{0}^{\infty}d\tau\,e^{-r\tau}\,\prod_{j=1}^{N}\,p(x_{j},\tau)\,.$ The strong, all-to-all correlation is clear from this unique non-factorized form of the JPDF. The explicit form of the evolved wavefunction $\psi(x,\tau)$ can be obtained by using the quantum propagator,

[TABLE]

Since both $K(x,y,\tau)$ and $\psi_{0}(y)$ are Gaussian, $\psi(x,\tau)=\int_{-\infty}^{\infty}K(x,y,\tau)\,\,\psi_{0}(y)\,dy$ remains Gaussian at all times. Consequently, the single particle conditional PDF $p(x,\tau)=|\psi(x,\tau)|^{2}$ is again a Gaussian,

[TABLE]

and $\tau$ appears through the variable $V(\tau)$ as

[TABLE]

It is evident from Eq. (6) that $V\in(V_{1},V_{2})$ with

[TABLE]

Therefore, NESS JPDF $\mathrm{Prob}.(x_{1},x_{2},\ldots x_{N})$ can be expressed in terms of the random variance $V$ as

[TABLE]

is the distribution of $V$, and $p(x|V)$ is given in Eq. (5). To proceed further, we must obtain $h(V)$ in Eq. (8), with $V(\tau)$ given in Eq. (6) and $\tau$ drawn from the exponential distribution $r\,e^{-r\tau}$. The evaluation of $h(V)$ requires us to express $\tau$ as a function of $V$ from Eq. (6). To do this, we must find the roots of

[TABLE]

for $\tau>0$. The smallest root, $\tau_{1}$ is given by, $\tau_{1}=\cos^{-1}(\alpha)/(2\,\omega_{2})$. Since $\cos(2\omega_{2}\tau)$ has a period of $\pi/\omega_{2}$ there are roots to Eq. (9) at $\tau_{2s+1}={s\pi}/{\omega_{2}}\,+\tau_{1},\,\,\text{with}\,\,s=1,2,\dotsc$ There is an additional root, denoted by $\tau_{2s}$, within each period. The two roots $\tau_{2s-1}$ and $\tau_{2s}$ within each period are symmetrically placed with respect to $(2s-1)\pi/(2\omega_{2})$, where $s=1,2,\dotsc$. This yields $\tau_{2s}={s\pi}/{\omega_{2}}\,-\tau_{1},\,\,\text{with}\,\,s=1,2,\dotsc\,$ Summing over all the roots, the PDF $h(V)$ in Eq. (8) becomes

[TABLE]

and $q=r/\omega_{2}$, $f=\omega_{1}/\omega_{2}$, and $V_{1}$ and $V_{2}$ are given in Eq. (7). Clearly, $h(V)$ has square-root divergences at both supports. In Fig. 1, we plot the PDF $h(V)$ in Eq. (10) for different values of $f$ and $q$ and compare with direct numerical simulations. We find excellent agreement between the two.

Henceforth, we set $\sigma_{1}=1$ and compute various macroscopic and microscopic observables for this model. These observables have been validated through numerics obtained by the method described below.

*Numerical method employed.— * The numerical data is obtained as follows: (i) We first pick a random number $\tau$ from an exponential distribution $re^{-r\tau}$. (ii) Next we obtain $V(\tau)$ defined in Eq. (6). To find the distribution of $V(\tau)$ in Fig 1, $10^{7}$ realizations of $V$ are obtained, and the data are binned and averaged. For the observables: (a) We pick $\tau$ from the exponential distribution and find the corresponding $V(\tau)$. (b) For each $V$, we generate $N$ Gaussian random numbers, taking $V$ as the variance. This procedure (a)-(b) is repeated for $10^{5}-10^{6}$ number of realizations, and the observables are computed, binned, and averaged over all the realizations.

3 Observables

3.1 Average Density Profile

The average density profile $\rho(x)=N^{-1}\langle\sum_{i=1}^{N}\delta(x_{i}-x)\rangle$ is nothing but the marginal single particle distribution. From the CIID structure in Eq. (8), it follows that

[TABLE]

Although $h(V)$ has a finite support, the average density profile in Eq. (11) is supported over the entire real line. The comparison of numerical data and the analytical average density profile (with the integral evaluated numerically) in Fig. 2 shows an excellent agreement.

The asymptotic behavior of $\rho(x)$ can be obtained by noting that due to the essential singularity of $V$ in $e^{-x^{2}/(2V)}$, for $x\gg\sqrt{V_{2}}$, the most dominant contribution to the integral in Eq. (11) comes from near the upper limit $V_{2}$. Therefore, we set $V=V_{2}-\epsilon$ and keep only the leading order terms in $\epsilon$ in the integrand in Eq. (11). From Eq. (10), for $\epsilon\to 0$, we get

[TABLE]

Hence, for $x\gg\sqrt{V_{2}}$, Eq. (11) yields

[TABLE]

where we have extended the upper limit of the integral to $\infty$ (since the integrand decays very fast). Hence, the tail remains Gaussian, albeit with a different variance and a logarithmic correction. The form of the tail in Eq. (13) is validated by simulation in Fig. 2.

3.2 Order Statistics

The order statistics $\mathrm{Prob}.(M_{k}=w)$ represent the position of the $k$-th furthest particle along the line, where $k=1$ corresponds to the rightmost particle, $k=2$ to the second rightmost particle, and so on. The order statistics for this class of models are obtained by first computing these for a fixed $V$. Since the conditional single particle distributions are Gaussian, and IID, an asymptotic analysis leads to the conditional distribution being $\mathrm{Prob}.(M_{k}|V)=\delta(M_{k}-\sqrt{2V}\mathrm{erfc}^{-1}(2\alpha))$, with $\alpha=k/N$ [22, 18]. Thus, averaging this over $V$ leads to

[TABLE]

with $h(z)$ defined in Eq. (10), and $\theta_{\alpha}=\mathrm{erfc}^{-1}(2\alpha)$. This scaling function is valid for both the bulk as well as the edge, and has finite support $M_{k}/\theta_{\alpha}\in\left(\min(f,1),\,\max(f,1)\right)/\sqrt{2}$. Note that for $\alpha>1/2$, $\theta_{\alpha}$ is negative, and the corresponding $M_{k}$’s have support over negative values, but are still described by the scaling function Eq. (14). Moreover, we unravel the interesting feature that although the density profile is supported over the real line, the order statistics has finite support [21, 16]. These results are validated in Fig. 3.

3.3 Gap Statistics

We compute the statistics of the gap between successive particles, $d_{k}=M_{k}-M_{k+1}$. This is done by fixing $V$, and analyzing the IID random variables [18]. After averaging over $h(V)$, this leads to

[TABLE]

with $\lambda=N^{-1}\sqrt{2\pi}\,e^{\theta_{\alpha}^{2}}$, $\theta_{\alpha}=\mathrm{erfc}^{-1}(2\alpha)$. The gap statistics results are validated in Fig. 4. The integrand in the scaling function $F(z)$ in Eq. (15) has a similar essential singularity as that in the average density profile in Eq. (11). Therefore, for $z\gg\sqrt{V_{2}}$, a similar analysis as in Sec. 3, yields

[TABLE]

where $A_{f}$ is given in Eq. (12). Thus, the tail of the scaling function of the gap statistics remains exponential, with the quenching protocol introducing a subleading logarithmic correction.

3.4 Full Counting Statistics

Exploiting the CIID structure, given a $V$, the probability of having $N_{L}$ particles within $[-L,L]$ is a binomial distribution, $\mathrm{Prob.}(N_{L},N|V)=\binom{N}{N_{L}}\,\phi(L|V)^{N_{L}}[1-\phi(L|V)]^{N-N_{L}}$. Here $\phi(L|V)=\int_{-L}^{L}dx\,p(x|V)$, where $p(x|V)$ is defined in Eq. (5). In the $N\rightarrow\infty$ limit, we know that the binomial distribution converges to a Gaussian distribution, and consequently the fraction of particles contained between $[-L,L]$ is a Dirac delta function about $\phi(L|V)=\mathrm{erf}(L/\sqrt{2V})$ [22]. Averaging this over $V$ gives us the required FCS, which reads

[TABLE]

This is supported over $[\kappa_{\min},\kappa_{\max}]$, with $\kappa_{\max}=\max[\mathrm{erf(L)},\,\mathrm{erf}(Lf^{-1})]$ and $\kappa_{\min}=\min[\mathrm{erf(L)},\,\mathrm{erf}(Lf^{-1})]$. Hence, in any volume, no more than $\kappa_{\max}$ and no less than $\kappa_{\min}$ fraction of particles can be contained in the thermodynamic limit of $N\to\infty$. These findings are validated in Fig. 5.

4 Conclusion and Outlook

We have studied a system of $N$ noninteracting bosons in a harmonic trap subjected to stochastic frequency quenches, in which the trap frequency alternates between two chosen values after random exponentially distributed durations. Despite the absence of direct interactions, the protocol induces strong dynamical correlations, driving the system to a unique NESS whose JPDF displays a CIID structure. The system and the driving protocol yield an effectively non-Hermitian setup. This stands out as a rare example of an open quantum system that is both analytically tractable and implementable in cold-atom experiments. The CIID mathematical structure we obtained enables exact computation of several observables, including density profiles, gap and order statistics, and full counting statistics, revealing both universal and parameter-dependent features.

Some future questions include incorporating inherent interactions, either perturbatively or within mean-field frameworks, to probe the consequences of departure from a CIID structure [32]. The intricate interplay between inherent interactions and dynamically generated correlations is expected to yield highly non-trivial steady states. It will also be interesting to explore non-exponential waiting-time distributions [33], which could lead to qualitatively different steady states. Finally, given that our setup and protocol in particularly suitable for ultracold atom experiments, combining state-of-the-art absorption imaging techniques [23, 24, 25, 26] also with high-resolution detection methods, such as quantum gas microscopy [27, 28, 29], could enable direct observation of the predicted correlations and statistics, paving the way towards controlled engineering of correlated open quantum systems.

5 Acknowledgements

MK acknowledges the support of the Department of Atomic Energy, Government of India, under project no. RTI4001. SNM and SS acknowledge the support from the Science and Engineering Research Board (SERB, Government of India), under the VAJRA faculty scheme (No. VJR/2017/000110). SNM acknowledges support from ANR Grant No. ANR-23-CE30-0020-01 EDIPS.