# Quenched dynamics of the momentum distribution of the unitary Bose gas

###### Abstract

We study the quenched dynamics of the momentum distribution of a unitary Bose gas under isotropic harmonic confinement within a time-dependent density functional approach based on our recently calculated Monte Carlo (MC) bulk equation of state. In our calculations the inter-atomic s-wave scattering length of the trapped bosons is suddenly increased to a very large value and the real-time evolution of the system is studied. Prompted by the very recent experimental data of Rb atoms at unitarity [Nature Phys. 10, 116 (2014)] we focus on the momentum distribution as a function of time. Our results suggest that at low momenta, a quasi-stationary momentum distribution is reached after a long transient, contrary to what found experimentally for large momenta which equilibrate on a time scale shorter than the one for three body losses.

###### Keywords:

Unitary Bose Gas Time Dependent Density Functional Local Density approximation Momentum Distribution^{†}

^{†}journal: Few-Body Systems

∎

## 1 Introduction

Recent experiments claim to have achieved a metastable degenerate gas of ultracold and dilute bosonic atoms with infinite s-wave scattering length [1; 2; 3; 4]. This is the so-called unitary Bose gas, which is characterized by remarkably simple universal laws arising from scale invariance [5]. While the unitary Fermi gas has been largely investigated both experimentally and theoretically [6], its bosonic counterpart has been only marginally theoretically addressed [7; 8; 9; 10; 11; 12; 13] because generally considered as experimentally inaccessible [14] due to the dominant three-body losses at very large values of the scattering length.

By using a Monte Carlo (MC) approach [13] we have recently investigated the zero-temperature properties of a dilute homogeneous Bose gas by tuning the interaction strength of the two-body potential to achieve arbitrary positive values of the s-wave scattering length , while avoiding the formation of the self-bound clusters present in the ground-state [13]. More recently, the MC equation of state has been the key ingredient of a local density approximation (LDA) [15] to the energy density functional to be used in density functional theory (DFT). Remarkably, the density profiles of a unitary Bose gas in a harmonic trap calculated with DFT using such functional compare very well with the MC ones. In addition, by using a time-dependent formulation, we have also investigated the excitation frequencies as a function of the scattering length. Interestingly, the calculated values for the monopole breathing mode, which reproduce the expected limiting values of the ideal and unitary regimes, exhibit a non-monotonous behavior as varies [15].

Unfortunately, the experimental data on the unitary Bose gas [1; 2; 3; 4] to compare with are quite scarce, and they do not provide an easy tool for validating theoretical approaches. The most suitable data in this sense are the time-resolved measurements of the momentum distribution of a Bose-condensed gas quenched at unitarity (sudden increase of to a very large value) provided by the very recent experiment of Makotyn et al. [4]. The main result coming from Ref. [4] is that, for where is the atom number density, evolves to a quasi-steady-state distribution on a time scale shorter than the one characterizing three-body losses. We thus try to reproduce the observed phenomenon by solving a nonlinear Schrödinger equation (NLSE) obtained from time-dependent DFT (TDDFT) with our LDA functional, based on the MC equation of state and where we have introduced a dissipative term to take into account three-body losses, whose role is relevant during the quenched dynamics [16]. We must notice however, that we expect our single orbital TDDFT to be fully reliable only at low momenta, where the collective long-wavelength dynamics dominates. Therefore our analysis will be confined to low momenta, in a range not directly comparable with the one studied experimentally. The results we find, are therefore complementary to those of Makotyn et al. [4] and suggest that, as expected, while the system equilibrates locally in a short time as found by experiments, it takes a very long time to reach quasi equilibrium all over the trap.

## 2 Method

On the basis of the density functional theory [17], a reliable energy functional of the local density for an inhomogeneous system of interacting bosons at is given by

(1) |

where the quantum-pressure gradient term takes into account effects due to density variations [18], is the energy per atom of the homogeneous system with density equal to the local density and describes the external confinement, which we assume to be an isotropic harmonic potential . The values of have been recently fitted to the results of a MC calculation [13] for a wide range of (positive) values of the scattering length characterizing the interparticle interaction. In the weakly interacting regime (, where is the average distance between bosons) the MC results for are very close to , the universal Bogoliubov prediction [19] as corrected by Lee, Huang and Yang (LHY) [20]. In the strong-coupling regime (), instead, MC data reach a plateau and, in the unitarity limit (), a finite and positive energy per particle is found, . The equation of state of the homogeneous system[13] from such MC calculation can be well interpolated as:

(2) |

In Ref. [15] we have shown that the energy functional (1) is very accurate in reproducing the MC static density profiles of the inhomogeneous unitary Bose gas under harmonic confinement. The dynamics of the system can be obtained by generalizing the energy functional (1) into a density-phase action functional , where is the phase of a single-valued quantum-mechanical wave function representing the macroscopic wave function of the Bose-Einstein condensate of the superfluid [24], which is related to the superfluid velocity by the relation [25]. Such a density-phase action functional is written as

(3) |

where

(4) |

is the kinetic Lagrangian of Popov [26], and is given by Eq. (1) under the assumption of a time-dependent local density. It is straightforward to derive the equations of superfluid hydrodynamics with a quantum-pressure term (directly related to the von Weizsacker gradient term of Eq. (1)) by extremizing the action functional (3) [27]. Moreover, introducing the wave function

(5) |

the equations of superfluid hydrodynamics can be re-written in terms of a nonlinear Schrödinger equation

(6) |

which is the Euler-Lagrange equation obtained from the action functional (3) taking into account Eq. (5).

We stress that Eq. (6) can be alternatively obtained from time dependent density functional theory (TDDFT) [17; 28; 29] in the Local Density Approximation using a single Kohn-Sham orbital to describe the degenerate boson system [27; 15]. The computational approach based on Eq. (6) has been adopted, for instance, to describe superfluid helium [30; 31], ultracold bosonic atoms [27; 15], and also superfluid fermions in the BCS-BEC crossover [32; 33].

In the dynamics of the unitary Bose gas three-body losses play a relevant dissipative role, especially close or at unitarity. As done in previous applications of Eq. (6) [16], we model the effect of three-body losses by adding a phenomenological dissipative term (with cm/sec [4]) in equation (6). Thus, the dissipative NLSE we use for our numerical simulations is given by

(7) |

We have numerically solved Eq.(7) to obtain the real-time evolution closely simulating the experimental conditions of Makotyn et al. [4]. Therefore we consider as initial configuration a cloud of atoms, confined in a spherical harmonic trap with frequency Hz and scattering length (where is the Bohr radius) prepared in its ground state by evolving in imaginary time Eq. (6) without the dissipative term (i.e. with ). Then we switch on the dissipative term (with cm/sec), increase the scattering length to a very high, but otherwise arbitrary, value , and let the system evolve in real time for a time . At such time we analyze the momentum distribution of the particles, as described in the following. Notice that according to our calculated MC equation of state (2), the chosen value for is practically equivalent to , i.e. it gives the same energy per atom as in the truly unitary limit[13].

## 3 Numerical results

Due to the presence of the three-body recombination, whose rate increases as , atoms are continuously disappearing from the trap. At early times, s, experimental data are compatible with an exponential decay with a time constant of about s [4]. In our simulations the number of particles in the trap suffers a similar depletion, driven by the phenomenological dissipative term proportional to . The time-dependent behavior of the total number of trapped atoms can be monitored by solving Eq.(7) and by computing

(8) |

The obtained is shown as a solid line in Fig. 1 and is compared to the experimental data of Ref. [4], reported as filled squares. Our DFT results are fairly comparable with the experimental points, where these are available, but we note that for later times our decay is far from being exponential. It seems to be characterized by two different time scales: one, fast, driving the first depletion of the atoms within the trap, that covers the experimental data range; and a second one, which is much slower, dictating the emptying of the trap at very large times. We find that, in this second regime, a relevant number of atoms still populate the trap, and the depletion is much slower, leading to an apparent stationary condition, as discussed in the following, and allowing for measurements of properties that give the impression of being equilibrated.

The size of trapped cloud predicted on the basis of our DFT calculation is also in reasonable agreement with the experimental data, given their uncertainty. In Fig. 2 we report the average radius of the bosonic cloud as a function of time . The solid line refers to the results obtained using the dissipative NLSE (7) with the MC equation of state (2) and . They are fully compatible with the experimental findings [4] reported as horizontal dashed lines in Fig. 2. For the sake of comparison, in Fig. 2 we also plot with a dotted line the results for the same average radius, computed using the dissipative NLSE Eq. 6 with the same finite scattering length and the same three-body loss coefficient as above, but with the energy density from the Gross-Pitaevskii equation of state instead than from MC. The figure clearly shows that the radius obtained within the GPE theory (dotted line) rapidly exits the experimental range as measured in the first 500s [4] and displays local oscillations due to interference effects in the outer part of the atomic cloud interacting with the trap walls.

Given the solution of the dissipative NLSE (7), the time dependent momentum distribution is easily calculated as:

(9) |

In Fig. 3 we plot obtained at different times, up to a maximum value of ms, with , where and is the density at the center of the trap. Notice that in the figure we report at low momenta (), because our single-orbital TDDFT is supposed to be fully reliable only at long wavelengths. During the first s we observe an expansion of the cloud both in real and in momentum spaces, and then a quasi-stationary momentum distribution seems to develop for larger times.

However, we must point out that at longer times (not shown in the figure) is still changing. This could be ascribed to the fact that, due to the sudden quench of to so high values, the bosonic cloud expands until the wave-front interacts with the steeper part of the trap walls and is reflected back, letting the cloud to shrink again at much larger times. Thus the momentum distribution first decreases (close to the turning point) and then it increases again as the contraction of the cloud occurs.

To better clarify the time scales involved in the momentum changes, we plot in Fig. 4 the time evolution of the number of atoms with a given momentum for different values of during the first after the sudden quench. We note that the evolution of the distribution for these small momenta is not uniform, being oscillatory first, with the amplitude of the oscillations extinguishing on a time scale of .

This is to contrasted with what found in the experiments of Ref.[4] for higher () components in the momentum distribution. There, a shorter time scale governs the saturation of the k-component, increasing from about to about as decreases from to . Our results complement these observations, in finding a still larger time-scale when . Our results can be compared also with recent studies on the dynamics of after a sudden quench to very large value of [34]. In particular, calculations based on a bath approach [34] lead to a multistep equilibration process, where larger modes equilibrate faster (having a shorter relaxation time), the oscillations in are damped on intermediate times, leading to an apparent equilibrated state, and the full equilibration is reached only on a much larger time scale, when also the condensate attains its final time-independent state.

## 4 Conclusions

We have numerically investigated the quenched dynamics of a Bose gas of Rb atoms under isotropic harmonic confinement after the sudden increase of the s-wave scattering length from to . To take into account three body losses we have introduced a dissipative term into the time dependent nonlinear Schrödinger equation (TDNLSE) obtained from a local approximation of a TD energy functional. We have shown that, while the equation obtained by using the Gross-Pitaevskii equation of state in the TD energy functional is unable to fully capture the experimental quenched dynamics of the atomic cloud, the corresponding equation derived from a functional containing the energy density fitted to our recently calculated Monte Carlo bulk equation of state [13] gives values of the average radius of the cloud in fairly good agreement with very recent experimental data [4]. We take this result as supporting the use of our MC equation of state in the energy density functional.

The solution of our dissipative TDNLSE shows a fast depletion of the condensate at short times, s due to three-body losses, accompanied by the relative increase in the population of large momenta. Also this feature is in agreement with the experimental findings [4]. As shown in Fig. 4 though, the momentum densities at various momenta reach a ’quasi equilibrium’ distribution, at very large times of about s. One must notice however that even the larger momenta reported in our figures are much smaller than those studied experimentally, corresponding to values for which our TDDFT is questionable. We reconcile our results with those reported in Ref. [4] by observing that they refer to different length scales: while from the experiment a fast equilibration on a short length scale has been found, reminiscent of local equilibrium, our calculation shows that a much longer time is needed to get quasi-equilibrium on a large length scale. The authors of Ref. [4] interpret their results (see their fig. 4) as a saturation of the value of at times of the order of s, shorter than the time scale for three body losses ( s), thus indicating a different mechanism as responsible for the equilibration of the momentum distribution. Such a mechanism could be the interaction between the condensed and non condensed fractions of boson atoms which is not present in our equation, since the dissipative term in our equation lets the total number of atoms to decrease while leaving the system in a degenerate state.

The effect of the interaction between the condensate and excited atoms on the dynamics of momentum distribution of a Bose gas after a rapid quench to a very large value of has been explicitly addressed at in some recent papers [34; 35; 36; 37; 38; 39] and supports the existence of different times domains in the dynamics of momentum distribution. In Ref.[34], in particular, they correspond to a fast depletion of the condensate followed or accompanied by a slow re-adjustment of the momentum distribution and, only at very large times, the reaching of a metastable state. We will explore the possibility of coupling the wave function obeying our TDNLSE to a bath [34] describing the small noncondensed fraction with the goal of better reproducing the experimental data at large momenta and short times.

The authors are grateful to J. P. Corson for a critical reading of a previous version of this report and for useful remarks. The authors acknowledge for partial support the Ministero Istruzione Università Ricerca (PRIN project 2010LLKJBX).

## References

- [1] Li w. and Ho T.-L. (2012) Bose Gases near Unitarity. Phys. Rev. Lett. 108, 195301.
- [2] Rem B.S., Grier A.T., Ferrier-Barbut I., Eismann U., Langen T., Navon N., Khaykovich L., Werner F., Petrov D.S., Chevy F. and Salomon C. (2013) Lifetime of the Bose Gas with Resonant Interactions. Phys. Rev. Lett. 110, 163202.
- [3] Fletcher R.J., Gaunt A.L., Navon N., Smith R.P. and Hadzibabic Z. (2013) Stability of a Unitary Bose Gas. Phys. Rev. Lett. 111, 125303.
- [4] Makotyn P., Klauss C.E., Goldberger D.L., Cornell E.A. and Jin D.S. (2014) Universal dynamics of a degenerate unitary Bose gas. Nature Phys. 10, 116.
- [5] Castin Y. and Werner F. (2012) The Unitary Gas and its Symmetry Properties. Lecture Notes in Physics Vol 836, Ed. Wilhelm Zwerger, Springer, Berlin, pp 127-191.
- [6] Zwerger W. Ed. (2012), The BCS-BEC Crossover and the Unitary Fermi Gas. Springer, Berlin.
- [7] Cowell S., Heiselberg H., Mazets I.E., Morales J., Pandharipande V.R. and Pethick C.J. (2002) Cold Bose Gases with Large Scattering Lengths. Phys. Rev. Lett. 88, 210403.
- [8] Adhikari S.K. and Salasnich L. (2008) Nonlinear Schrödinger equation for a superfluid Bose gas from weak coupling to unitarity: Study of vortices. Phys. Rev. A 77, 033618.
- [9] Song J.L. and Zhou F. (2009) Ground State Properties of Cold Bosonic Atoms at Large Scattering Lengths. Phys. Rev. Lett. 103, 025302.
- [10] Lee Y.-L. and Lee Y.-W. (2010) Universality and stability for a dilute Bose gas with a Feshbach resonance. Phys. Rev. A 81, 063613.
- [11] Diederix J.M., van Heijst T.C.F. and Stoof H.T.C. (2011) Ground state of a resonantly interacting Bose gas. Phys. Rev. A 84, 033618.
- [12] van Heugten J.J.R.M. and Stoof H.T.C. (2013) Theory of unitary Bose gases. arXiv:1302.1792.
- [13] Rossi M., Salasnich L., Ancilotto F. and Toigo F. (2014) Monte Carlo simulations of the unitary Bose gas. Phys. Rev. A 89, 041602(R).
- [14] Ho T.-L. (2004) Universal Thermodynamics of Degenerate Quantum Gases in the Unitarity Limit. Phys. Rev. Lett. 92, 090402.
- [15] Rossi M., Ancilotto F., Salasnich L. and Toigo F. Density functional theory of a trapped Bose gas with tunable scattering length: from weak-coupling to unitarity. arXiv:1408.3945 to be published in Eur. Phys. J..
- [16] Lahaye T., Metz J., Frolich B., Koch T., Meister M., Griesmaier A., Pfau T., Saito H., Kawaguchi Y. and Ueda M. (2008) d-Wave Collapse and Explosion of a Dipolar Bose-Einstein Condensate. Phys. Rev. Lett. 101, 080401.
- [17] Hohenberg P. and Kohn W. (1964) Inhomogeneous Electron Gas. Phys. Rev. 136, B864.
- [18] von Weizsäcker C.F. (1935) Zur theorie der kernmassen. Z. Phys. 96 431.
- [19] Bogoliubov N.N. (1947) On the Theory of Superfluidity. J. Phys. (USSR) 11, 23.
- [20] Lee T.D., Huang K. and Yang C.N. (1957) Eigenvalues and Eigenfunctions of a Bose System of Hard Spheres and Its Low-Temperature Properties. Phys. Rev. 106, 1135.
- [21] The coefficients ’s have been determined by ensuring continuity of the interpolating function and of its first three derivatives at and once the and the ’s have been fitted to the calculated values. This procedures automatically well reproduces the MC values in the interval .
- [22] Gross E.P. (1961) Structure of a quantized vortex in boson systems. Nuovo Cimento 20 454;
- [23] Pitaevskii L.P. (1961) Vortex lines in an imperfect Bose gas. Soviet Physics JETP 13, 451.
- [24] Leggett A.J. (2006) Quantum Liquids. Oxford Univ. Press, Oxford.
- [25] Landau L.D. and Lifshitz E.M. (1987) Fluid Mechanics, Course of Theoretical Physics, vol. 6. Pergamon Press, London.
- [26] Popov V.N. (1983) Functional Integrals in Quantum Field Theory and Statistical Physics. Reidel, Dordrecht.
- [27] Kim Y.E. and Zubarev A.L. (2003) Density-functional theory of bosons in a trap. Phys. Rev. A 67, 015602.
- [28] Kohn W. and Sham L.J. (1965) Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 140, A1133.
- [29] Runge E. and Gross E.K.U. (1984) Density-Functional Theory for Time-Dependent Systems. Phys. Rev. Lett. 52, 997.
- [30] Dalfovo F., Lastri A., Pricaupenko L., Stringari S. and Treiner J. (1995) Structural and dynamical properties of superfluid helium: A density-functional approach. Phys. Rev. B 52, 1193.
- [31] Ancilotto F., Dalfovo F., Pitaevskij L. and Toigo F. (2005) Density pattern in supercritical flow of liquid He-4. Phys. Rev. B 71, 104530.
- [32] Kim Y.E. and Zubarev A.L. (2004) Time-dependent density-functional theory for trapped strongly interacting fermionic atoms. Phys. Rev. A 70, 033612.
- [33] Manini N. and Salasnich L. (2005) Bulk and collective properties of a dilute Fermi gas in the BCS-BEC crossover. Phys. Rev. A 71, 033625.
- [34] Rançon A. and Levin K. (2014) Equilibrating dynamics in quenched Bose gases: Characterizing multiple time regimes. Phys. Rev. A 90, 021602(R).
- [35] Yin X. and Radzihovsky L. (2013) Quench dynamics of a strongly interacting resonant Bose gas. Phys. Rev. A 88, 063611.
- [36] Kain B. and Ling H.Y. (2014) Nonequilibrium states of a quenched Bose gas. Phys. Rev. A 90, 063626.
- [37] Sykes A.G., Corson J.P., D’Incao J.P., Koller A.P., Greene C.H., Rey A.M., Hazzard K.R.A. and Bohn J.L. (2014) Quenching to unitarity: Quantum dynamics in a three-dimensional Bose gas. Phys. Rev. A 89, 021601(R).
- [38] Kira M. (2014) Excitation picture of an interacting Bose gas, Ann. Phys. 351, 200.
- [39] Corson J.P. and Bohn J.L. (2015) Bound-state signatures in quenched Bose-Einstein condensates, Phys Rev A 91, 013616.