Universal Dynamics of a Degenerate Bose Gas Quenched to Unitarity

TL;DR

This paper investigates the universal dynamics of a degenerate Bose gas quenched to unitarity, revealing exponential momentum distribution behavior that is independent of interaction details, consistent with experimental observations.

Contribution

The study connects two-body and many-body dynamics, demonstrating universal exponential behavior in the momentum distribution after a quench to unitarity.

Findings

Long-time momentum distribution exhibits universal exponential decay.

The exponential behavior is linked to a sub-leading term in the initial wave function.

Results are consistent with experimental observations from Cambridge.

Abstract

Motivated by an unexpected experimental observation from the Cambridge group, [Eigen {\it et al.,} Nature {\bf563}, 221 (2018)], we study the evolution of the momentum distribution of a degenerate Bose gas quenched from the weakly interacting to the unitarity regime. For the two-body problem, we establish a relation that connects the momentum distribution at a long time to a sub-leading term in the initial wave function. For the many-body problem, we employ the time-dependent Bogoliubov variational wave function and find that, in certain momentum regimes, the momentum distribution at long times displays the same exponential behavior found by the experiment. Moreover, we find that this behavior is universal and independent of the short-range details of the interaction potential. Consistent with the relation found in the two-body problem, we also numerically show that this exponential form is hidden in the same sub-leading term of the Bogoliubov wave function in the initial stages. Our results establish a consistent picture to understand the universal dynamics observed in the Cambridge experiment.

Equations

$$n_{\mathbf{k}}\sim e^{-\Lambda k/k_n},$$

$$n_{\mathbf{k}}(t)=|⟨\mathbf{k}|e^{-\frac{i}{\hslash }\hat{H}t}|{\varphi }^i⟩{|}^2=|⟨\mathbf{k}|e^{\frac{i}{\hslash }{\hat{H}}_{0}t}{e}^{-\frac{i}{\hslash }\hat{H}t}|{\varphi }^i⟩{|}^2,$$

$${\hat{\Omega }}^{(-)}=\underset{t\to +\infty }{\lim }{e}^{\frac{i}{\hslash }\hat{H}t}{e}^{-\frac{i}{\hslash }{\hat{H}}_{0}t}$$

$$n_{\mathbf{k}}(t\to +\infty )=|⟨\mathbf{k}|{\hat{\Omega }}^{(-)†}|{\varphi }^i⟩{|}^2=|⟨{\mathbf{k}}^{(-)}|{\varphi }^i⟩{|}^2.$$

$$|{\mathbf{k}}^{(-)}⟩=|\mathbf{k}⟩+\frac{1}{{ϵ}_{\mathbf{k}}+i{0}^{-}-\hat{H}}\hat{V}|\mathbf{k}⟩,$$

$$⟨\mathbf{r}|{\mathbf{k}}^{(-)}⟩=\frac{1}{(2\pi {)}^{3/2}}\left(e^{i\mathbf{k}\cdot \mathbf{r}}+\frac{1}{ik}\frac{e^{-ikr}}{r}\right).$$

$$\psi (\mathbf{k})=e^{i\theta (\mathbf{k})}\sqrt{n_{\mathbf{k}}(t\to +\infty )},$$

$$|{\varphi }^i⟩=\int d^3\mathbf{p}\psi (\mathbf{p})|{\mathbf{p}}^{(-)}⟩,$$

$${\varphi }^i(\mathbf{k})=⟨\mathbf{k}|{\varphi }^i⟩=\frac{1}{(2\pi {)}^{\frac{3}{2}}}\int d^3\mathbf{r}{d}^3\mathbf{p}{e}^{-i\mathbf{k}\cdot \mathbf{r}}\psi (\mathbf{p})⟨\mathbf{r}|{\mathbf{p}}^{(-)}⟩.$$

$${\varphi }^i(k)=\frac{1}{2\pi }\left(-\frac{i}{k}\right)\underset{ϵ\to 0^{+}}{\lim }{\int }_0^{+\infty }dpp\psi (p)$$

$$\times \left[\sum _{\sigma =\pm }\frac{\sigma }{p+\sigma (k+iϵ)}+\sum _{{\sigma }^{\prime }=\pm }\frac{{\sigma }^{\prime }}{p+{\sigma }^{\prime }(k-iϵ)}\right].$$

$${\varphi }^i(k)=-\psi (k)-\frac{1}{k}\sum _j\frac{\text{Res}{\left[2z\Psi (z)\right]}_{z=z_j}}{z_j-k},$$

$$\hat{\mathcal{H}}=\sum _{\mathbf{k}}{ϵ}_{\mathbf{k}}{\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_{\mathbf{k}}+\frac{1}{2L^3}\sum _{\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}}{\hat{a}}_{\mathbf{k}+\mathbf{q}}^{†}{\hat{a}}_{{\mathbf{k}}^{\prime }-\mathbf{q}}^{†}V(\mathbf{q}){\hat{a}}_{{\mathbf{k}}^{\prime }}{\hat{a}}_{\mathbf{k}}.$$

$$|\Phi (t)⟩=\mathcal{A}(t)\exp \left[g_0(t){\hat{a}}_0^{†}+\sum _{\mathbf{k}\cdot \hat{z}>0}{g}_{\mathbf{k}}(t){\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_{-\mathbf{k}}^{†}\right]|0⟩.$$

$$i\hslash {\dot{g}}_0$$

$$i\hslash {\dot{g}}_{\mathbf{p}}$$

$$+\frac{1}{L^3}\sum _{\mathbf{k}\ne 0}V(\mathbf{p}-\mathbf{k})\frac{2|g_{\mathbf{k}}|{}^2{g}_{\mathbf{p}}+g_{\mathbf{k}}+g_{\mathbf{k}}^{\ast }{g}_{\mathbf{p}}^2}{1-|g_{\mathbf{k}}{|}^2}.$$

$${\tilde{n}}_{\mathbf{k}}(t)=\frac{N_{\mathbf{k}}(t)}{n},$$

Full text

Because the interaction in cold atomic systems can be controlled by optical and magnetic fields, it can be tuned in a time scale that is much shorter than the relaxation time. Cold atomic systems are also very clean and the microscopic interaction between atoms can be understood very well in terms of universal low-energy effective interactions. Because of these two reasons, cold atoms are ideal for studying far-from-equilibrium dynamics in the many-body sector from a microscopic point of view. In equilibrium, there exists a lot of phenomena that are universally applicable, independently of the details of interactions at the microscopic scale. A major question for the non-equilibrium physics is whether such universal phenomena can also be found in far-from-equilibrium situations.

A recent experiment on strongly interacting Bose gases reveals a great surprise exp . The system is initially prepared as a nearly pure Bose condensate at very low temperature and with weak interactions. Then the interaction is changed abruptly and the system is quenched to unitarity with the $s$-wave scattering length being infinite. The subsequent many-body dynamics was monitored by observing the evolution of the momentum distribution $n_{{\bf k}}$. A prethermalization stage was found where $n_{{\bf k}}$ remains a constant for a long time. The most surprising finding in the experiment is that $n_{{\bf k}}$ has a functional form

[TABLE]

where $k=|{\bf k}|$, $k_{n}=(6\pi^{2}n)^{1/3}$ is a momentum scale related to the total density $n$, and $\Lambda=3.62$ is obtained from fitting the experimental data exp . This functional form is seen to be valid for $k$ ranging from $\sim k_{n}$ to a few times $k_{n}$. There are a number of previous theoretical works that have studied weakly interacting Bose gases quenched to the strongly interacting regime theory1 ; theory2 ; theory3 ; theory4 ; theory5 ; theory6 ; theory7 ; theory8 ; theory9 ; theory10 ; theory11 , with either finite or infinite scattering lengths. However, this phenomenon has not been predicted by any theory before.

Here we should note that at unitarity, because there is no other length scale, both the two-body collisional rate and the three-body loss rate are proportional to $E_{n}$, where $E_{n}$ is given by $\hbar^{2}k^{2}_{n}/(2m)$. However, it has been shown previously that the coefficient for the two-body collisional rate is usually larger than that for the three-body loss rate Rem13 ; Fletcher13 ; Jin ; Eigen17 , such that the many-body dynamics is governed by two-body collisions for a reasonably long time before the three-body loss takes over and heats the system up. Therefore, it is very reasonable to view both the prethermalization and processes occurring before that as caused by the two-body collisions, while temperature increasing at later times is due to the three-body losses. This separation of time scales allows us to safely ignore the three-body loss and only focus on two-body collisional effects when analyzing the prethermal dynamics. In another words, we can view the prethermalization regime as the long time limit of the dynamics governed by the two-body collisions. Moreover, we have also ignored the coherent three-body effect like the Efimov effect, because usually when two-body interaction is dominated, the three-body effect only provides a small correction.

In this letter we focus on understanding of the origin of the emergent exponential behavior of $n_{{\bf k}}$ and answering whether this functional form is universal or not, and we address this issue from both two-body and many-body perspectives. The main results can be summarized as follows:

I. For the two-body problem, we prove a relation between the long time behavior of the momentum distribution and the properties of the initial wave function. This relation works for arbitrary short-range potentials. With this theorem, we can determine which property of the initial wave function is responsible for the exponential form of Eq. (1) in the momentum distribution at a long time.

II. For the many-body problem, we employ a variational time-dependent Bogoliubov wave function and by solving the time-dependent equation, we find that for a certain range of momentum, the averaged $n_{{\bf k}}$ for long time evolution indeed obeys the form of Eq. (1), although the coefficient $\Lambda$ is quantitatively different from the experimental value due to the mean-field nature of our ansatz. We use three different potentials tuned near the vicinity of a scattering resonance, that are the square well, the Gaussian potential and the Yukawa potential, and find that this behavior is independent of the short-range details.

At the end, we also discuss the connection between the two-body and the many-body results.

Two-Body Problem. Let us first start with the two-body problem whose Hamiltonian can be written in terms of the relative coordinate ${\bf r}$ as ${\hat{H}}={\hat{H}}_{0}+{\hat{V}}({\bf r})$, where $\hat{H}_{0}=-\hbar^{2}\nabla^{2}/m$ is the kinetic energy with $m$ being the mass, and $\hat{V}({\bf r})$ is a short-range potential and we only consider the $s$-wave interaction. Here we choose $t=0$ as the time right after the quench of interactions and we denote the initial wave function as $|\phi^{i}\rangle$. The momentum distribution $n_{{\bf k}}(t)$ at momentum ${\bf k}$ and time $t$ is given by

[TABLE]

where $|{\bf k}\rangle$ is a plane wave state. The second equality follows from the fact that $|{\bf k}\rangle$ is an eigenstate of $\hat{H}_{0}$ and $e^{\frac{i}{\hbar}\hat{H}_{0}t}$ only gives rise to a phase factor that does not change $n_{{\bf k}}$. Furthermore, making use of the properties of the Møller operator Taylor

[TABLE]

in scattering theory, we can derive the following relation supple

[TABLE]

Here $|{\bf k}^{(-)}\rangle$ is the inward scattering wave function defined as Taylor

[TABLE]

where $\epsilon_{{\bf k}}=\hbar^{2}{\bf k}^{2}/m$. For a short-range potential, it is straightforward to show that outside the range of interaction, $\langle{\bf r}|{\bf k}^{(-)}\rangle$ behaves as Taylor

[TABLE]

where $r=|{\bf r}|$ and $k=|{\bf k}|$. Here we have explicitly used the fact that the system is quenched to unitarity with $a_{\text{s}}=\infty$.

With the relation Eq. (4), we can determine the requirement for the initial wave function $|\phi^{i}\rangle$ that can lead to the long-time behavior of Eq. (1) in $n_{{\bf k}}$. It is important to know that {$|{\bf k}^{(-)}\rangle$} also form a complete and orthogonal basis, and we can expand the wave function in terms of this basis. Let us introduce

[TABLE]

then the exponential form of $n_{{\bf k}}$ will translate to the same kind of exponential dependence for $\psi({\bf k})$ up to a phase factor. We can then write the initial wave function $|\phi^{i}\rangle$ as

[TABLE]

and in the momentum space

[TABLE]

Considering the situation where $\psi({\bf p})$ is isotropic, i.e. it can be written as $\psi(p)$ with $p=|{\bf p}|$, we can substitute Eq. (6) into Eq. (9) and integrate out the azimuthal degrees of freedom, we find supple

[TABLE]

Let us introduce an auxiliary function $\Psi(z)$ in the complex plane, such that it satisfies the requirement at the positive side of the real axis $\Psi(z=p>0)=\psi(p)$ and at its negative side $\Psi(z=p<0)=-\psi(-p)$ note . With the help of this auxiliary function, it can be shown that

[TABLE]

where $\text{Res}[f(z)]_{z=z_{j}}$ denotes the residue of the function $f(z)$ at its pole $z_{j}$.

Eq. (11) is a very interesting result. Here we should note that the amplitude of $\psi(k)$ obeys this universal exponential form only at the momentum $\gtrsim k_{n}$, however, the auxiliary function $\Psi(z)$ will certainly depend on the small momentum behavior of $\psi(k)$. Therefore, the residues of $\Psi(z)$, as well as the coefficient for this second term in the r.h.s. of Eq. (11), are non-universal. When $\psi(k)$ is a regular function as in Eq. (1) and Eq. (7), the second term recovers the well-known $1/k^{4}$ behavior of the momentum distribution at large $k$. Hence, Eq. (11) tells us that, one can subtract the leading order $1/k^{2}$ term from fitting the large momentum, and the remaining regular sub-leading term reveals the momentum distribution at a long time. That is to say that, in order for the long time behavior of $n_{{\bf k}}$ to obey Eq. (1), the sub-leading term of the initial wave function has to obey the form given by Eq. (1) and Eq. (7).

Many-Body Problem. Now we turn into the many-body problem whose Hamiltonian can be written in second quantized form as

[TABLE]

Here $\hat{a}_{{\bf k}}^{\dagger}$($\hat{a}_{{\bf k}}$) is the creation (annihilation) operator for bosons with momentum ${\bf k}$. $L^{3}$ is the system’s volume. $V({\bf q})=\int d^{3}{\bf r}e^{i{\bf q}{\bf r}}V({\bf r})$ is the Fourier transform of the interaction potential $V({\bf r})$. We implement the Bogoliubov-type variational ansatz

[TABLE]

Here $\mathcal{A}(t)$ is a normalization factor, $|0\rangle$ is the particle vacuum; $g_{0}$ and $g_{{\bf k}}$ are variational parameters, which can determine $N_{0}(t)=|g_{0}|^{2}$ and $N_{{\bf k}}(t)=|g_{{\bf k}}|^{2}/(1-|g_{{\bf k}}|^{2})$ as the particle number at zero-momentum and finite momentum ${\bf k}$ modes. The Bogoliubov ansatz assumes that the system remains as a Bose condensate during the entire dynamics, which is indeed the case for this experiment.

The dynamical equations for $g_{0}(t)$ and $g_{{\bf k}}(t)$ can be obtained from the Euler-Lagrange equation for the Lagrangian $\mathcal{L}=\frac{i\hbar}{2}[\langle\Phi|\dot{\Phi}\rangle-\langle\dot{\Phi}|\Phi\rangle]-\langle\Phi|\hat{H}|\Phi\rangle$. It results in coupled equations for $g_{0}$ and $g_{{\bf k}}$ as theory2

[TABLE]

The total number $N=N_{0}+\sum_{{\bf k}\neq 0}N_{{\bf k}}(t)$ is a conserved quantity, and $n=N/L^{3}$ is the total density. Making use of the spherical symmetry of this system, we consider that $g_{{\bf k}}$ only depends on $k$ and can be simplified as $g_{k}$, and we can further simplify this equation by performing the azimuthal integration first supple . Without loss of generality, we take the initial state to be a pure Bose-Einstein Condensate (BEC), i.e. $g_{0}(0)=\sqrt{N}$ and all $g_{{\bf k}}(0)=0$. As in the two-body case, we start the time evolution right after the interaction quench, and therefore we set the interaction at scattering resonance.

Here, to verify whether the dynamics is universal, that is to say, whether it depends on the short-range details, we consider three different short-range potentials:

(i) The square well potential:

$V_{\text{SW}}({\bf r})=-\frac{\hbar^{2}\gamma_{\text{s}}}{mr_{0}^{2}}\Theta\left(r_{0}-r\right)$, where $\Theta$ is the Heaviside step function. The $s$-wave resonance occurs at $\gamma_{\text{s}}=(\pi/2)^{2}$.

(ii) The Gaussian potential:

$V_{\text{GW}}({\bf r})=-\frac{\hbar^{2}\gamma_{\text{g}}}{mr^{2}_{0}}e^{-r^{2}/r_{0}^{2}}$, and the $s$-wave resonance occurs at $\gamma_{\text{g}}\approx 2.68$.

(iii) The Yukawa potential:

$V_{\text{YW}}({\bf r})=-\frac{\hbar^{2}\gamma_{\text{y}}}{mr_{0}}\frac{e^{-r\left/r_{0}\right.}}{4\pi r}$, and the $s$-wave resonance occurs at $\gamma_{\text{y}}\approx 21.1$.

We numerically solve the coupled equations Eq. (14) with these three potentials by discretizing both the radial momentum and the time, from which we can obtain $g_{{\bf k}}(t)$ and $N_{{\bf k}}(t)$. Following Ref. exp , we introduce a normalized momentum distribution

[TABLE]

such that $\frac{1}{L^{3}}\sum_{{\bf k}}\tilde{n}_{{\bf k}}(t)=1$. In Fig. 1 we plot $\tilde{n}_{{\bf k}}(t)$ as a function of $t$. One can see that following a growth at the initial stage, $\tilde{n}_{{\bf k}}(t)$ exhibits an oscillatory behavior for $t\gg t_{n}$, where $t_{n}$ is a typical time scale defined as $t_{n}=\hbar/E_{n}$. We also find that this oscillatory solution is stable against noises. Though it looks surprising that the highly non-linear equations can display stable oscillatory solution, it can be understood analytically in term of a simplified version of these coupled equations supple . In reality, the interactions between quasi-particles cause Baliaev-Landau damping, which will eventually smear out the oscillation and lead to a saturation result. Here, we take a long time average of $\tilde{n}_{{\bf k}}(t)$ starting from the second peak in the oscillation, as indicated in Fig. 1. The average is denoted by $\bar{n}_{{\bf k}}$, which is taken as the long time saturation value of the momentum distribution.

In Fig. 2, we plot the dimensionless quantity $k_{n}^{3}\bar{n}_{{\bf k}}$ as a function of $k/k_{n}$. The fit shows a regime around $k\sim k_{n}$, where $\bar{n}_{{\bf k}}$ behaves as Eq. (1), consistent with the experimental observation in Ref. exp . This fitting yields a coefficient $\Lambda=5.1-5.2$. For large $k$, the fitting yields a $1/k^{4}$ behavior. Most importantly, we note that the curves obtained using the three different potentials defined above collapse onto one another, which shows that this emergent exponential behavior of the momentum distribution is independent of the short-range details of the interaction potentials.

Connection between the Two- and Many-Body Problems. To summarize the results above, on one hand, our discussion on the two-body problem has established a relation between $n_{{\bf k}}$ at the long time and the wave function behavior at the initial time; and on the other hand, our Bogoliubov calculation for the many-body problem has discovered the exponential form for $n_{{\bf k}}$ at the long time as observed in the experiment reported in Ref. exp . Now a natural question is whether the same relation also holds in the Bogoliubov wave function, namely, whether the exponential form is also hidden in the sub-leading term of the Bogoliubov wave function at the early time. To check this conjecture, we look into the wave function $g_{k}$ at times $t<t_{n}$, far before the saturation of the momentum distribution, as indicated by arrows in Fig. 1. This early stage wave function is reminiscent of the initial wave function in the two-body case. In Fig. 3(a), we plot $|g_{k}|$ as a function of $k/k_{n}$, which does not show any exponential behavior, and the large $k$ part can be well fitted by a $1/k^{2}$ tail. In Fig. 3(b), following the same spirt of Eq. (11) discovered in the two-body problem, we subtract the $1/k^{2}$ part in $|g_{k}|$, and plot the sub-leading term as a function of $k/k_{n}$. Interestingly, in the same momentum range where the long time $\bar{n}_{{\bf k}}$ plotted in Fig. 2 shows an exponential behavior, this sub-leading term in the early time wave function plotted in Fig. 3(b) also shows an exponential behavior. Notice that we start from an initial state with all atoms in the zero-momentum state, the exponential behavior at the early stage wave function may originate from the pair production process, as discussed in Ref. Cheng .

Comments on Comparison with Experiment. Aside from momentum distribution, we find that the growth time also follows the scaling law as discovered by the experiment, and we find that the condensed fraction at long time is not vanishing and is about $10\%$. However, we should emphasize that the agreement between our Bogoliubov theory and the experiment is only qualitative. Experimentally, this exponential behavior of $n_{{\bf k}}$ is valid up to $\sim 3k_{n}$ and they do not find $1/k^{4}$ behavior, but in our case it is only up to $\sim k_{n}$ and is followed by a $1/k^{4}$ tail at higher momenta. The value of $\Lambda$ is also somewhat different between our calculation and the experimental result. However, since the system is strongly interacting, we do not expect the mean-field type Bogoliubov theory to be quantitatively accurate anyway. Moreover, our calculation leads to a very fast oscillation of $n_{{\bf k}}$ at long times and its mean value saturates. In experiment, the momentum distribution eventually takes off again after a prethermalization plateau, and this is caused by the heating due to the three-body loss which we do not include in our theory. Our results offer valuable insight for understanding this observation but more involved theories are required for a more quantitative comparison with experiment.

Note Added. When finishing this paper, we became aware of another preprint, Ref. Meera , which also did the many-body calculation with the Bogoliubov wave function.

Acknowledgment. We thank Wei Zheng, Manuel Valiente, Xin Chen, Zeng-Qiang Yu, Shizhong Zhang, Ran Qi for inspiring discussion. We are particularly grateful to Manuel Valiente for carefully reading our manuscript and helpful suggestions. This work is supported by NSFC Grant No. 11734010, 11604300, 11835011, 11774315, Beijing Outstanding Young Scholar Program and MOST under Grant No. 2016YFA0301600.