Undamped nonequilibrium dynamics of a nondegenerate Bose gas in a 3D isotropic trap

The canonical Boltzmann equation [22, 23] describes the phase space evolution of a cloud of thermal atoms obeying classical statistics:

$$\begin{eqnarray}&&\displaystyle \frac{\partial f({\bf{r}},{{\bf{v}}}_{1},t)}{\partial t}+\{f({\bf{r}},{{\bf{v}}}_{1},t),H\}={I}_{\mathrm{coll}}[f],\end{eqnarray} \tag{ 1 }$$

where $f({\bf{r}},{{\bf{v}}}_{1},t)$ is the single particle phase space distribution that depends on the generalized coordinate ${\bf{r}}$ and velocity ${{\bf{v}}}_{1}$ of all the particles, H is the single particle Hamiltonian, and $\{\}$ denotes the classical Poisson bracket. The left-hand side describes the single-particle evolution of the thermal cloud and the right-hand side is the collision integral which drives the system toward a state of statistical equilibrium

$$\begin{eqnarray}&&{I}_{\mathrm{coll}}[f]=\displaystyle \frac{\sigma }{4\pi }\int {{\rm{d}}}^{2}{\rm{\Omega }}{{\rm{d}}}^{3}{v}_{2}| {{\bf{v}}}_{2}-{{\bf{v}}}_{1}| [f({{\bf{v}}}_{1}^{\prime })f({{\bf{v}}}_{2}^{\prime })-f({{\bf{v}}}_{1})f({{\bf{v}}}_{2})]\end{eqnarray} \tag{ 2 }$$

at a rate determined by the total scattering cross section σ. The velocities of the incoming particles are ${{\bf{v}}}_{1}$ and ${{\bf{v}}}_{2}$, and the outgoing particles leave with velocities ${{\bf{v}}}_{1}^{\prime }$ and ${{\bf{v}}}_{2}^{\prime }$. For binary elastic collisions, both the center of mass (COM) momentum and the magnitude of the relative momentum are conserved, and the solid angle integral is over the full range of possible final directions for the relative momentum.

As a direct consequence of Boltzmann's H-theorem, the necessary condition for the vanishing of the collision integral is that $\mathrm{log}(f)$ be composed of quantities that are elements in the set of collisional invariants $\{{\chi }_{1},{\chi }_{2},\ldots \}$. This defines a class of solutions of the form

$$\begin{eqnarray}&&\mathrm{log}(f({\bf{r}},{\bf{v}},t))={\chi }_{1}+{\chi }_{2}+\ldots .\end{eqnarray} \tag{ 3 }$$

For a binary elastic collision, this set is spanned by the kinetic energy, momentum, and any velocity independent function, giving

$$\begin{eqnarray}&&\mathrm{log}\left(f({\bf{r}},{\bf{v}},t)\right)=a+{\bf{b}}\cdot {\bf{v}}+c{{\bf{v}}}^{2},\end{eqnarray} \tag{ 4 }$$

where $a,{\bf{b}},{\bf{c}}$ are possibly functions of position and/or time. The familiar Maxwell–Boltzmann distribution of a stationary dilute gas follows from this expression for the choice $| {\bf{b}}| =0$, $c={-1/k}_{{\rm{b}}}T$ and $a=\mathrm{log}n{\left(m/2\pi {k}_{{\rm{b}}}T\right)}^{3/2}$.

Although the parameters that lead to the Maxwell–Boltzmann distribution are the most well-known; it should be emphasized that Boltzmann noted a curious class of parameters that corresponded to collective motion in a harmonic trap [6]. Boltzmann considered solutions where c had periodic time-dependence describing temperature oscillations, as well as a position and time-dependent ${\bf{b}}$. This class of solutions corresponds to the so-called breathing or monopole mode and describes undamped collective motion of a thermal cloud in an isotropic 3D trap, regardless of the oscillation amplitude and interaction strength. In a recent paper [24], the class of undamped nonequilibrium solutions permitted by the H-theorem was extended to time-dependent confinement. For the remainder of the paper only the case of harmonic confinement is discussed.

2.1.1. Collective modes in a harmonic trap

The collective modes of a 3D isotropic trap with trapping frequency ${\omega }_{0}$ reflect the underlying spherical symmetry and the quadratic form of the trapping potential. Any mode of the trap can be expressed in terms of irreducible spherical tensors, in this case they are multipole modes denoted by the labels (l, m). The l = 0 mode is the monopole mode of the system, and is characterized by an undamped oscillation at $2{\omega }_{0}$. The l = 1 dipole mode is a COM oscillation independent of the character of the interaction and therefore undamped. The l = 2 quadrupole mode oscillates at $2{\omega }_{0}$ in the CL regime and $\sqrt{2}{\omega }_{0}$ in the HD regime. It is collisionally damped in the HD crossover between the two extremes. For small amplitude oscillations about equilibrium, one finds a general form for the oscillation frequency $\omega ={\omega }_{r}+{\rm{i}}{\rm{\Gamma }}$ of the collective mode [25]:

$$\begin{eqnarray}&&{\omega }^{2}={\omega }_{\mathrm{CL}}^{2}+\displaystyle \frac{{\omega }_{\mathrm{HD}}^{2}-{\omega }_{\mathrm{CL}}^{2}}{1+{\rm{i}}\omega {\tilde{\tau }}_{\mathrm{coll}}},\end{eqnarray} \tag{ 5 }$$

where the subscripts denote the mode frequency in the hydrodynamic (HD) and collisionless (CL) limits, and ${\tilde{\tau }}_{\mathrm{coll}}$ is related to ${\tau }_{\mathrm{coll}}$, the thermal relaxation time, by ${\tilde{\tau }}_{\mathrm{coll}}=({\omega }_{\mathrm{CL}}^{2}/{\omega }_{\mathrm{HD}}^{2}){\tau }_{\mathrm{coll}}$. The imaginary part of ω corresponds to the damping and the real part to the mode frequency. This general form is common to many kinds of relaxation phenomenon [26], and, assuming ${\rm{\Gamma }}\ll {\omega }_{r}$ can be written in an asymptotic form near the collisional and HD regimes:

$$\begin{eqnarray}{{\rm{\Gamma }}}_{\mathrm{HD}} & \approx & \displaystyle \frac{{\tilde{\tau }}_{\mathrm{coll}}}{2}({\omega }_{\mathrm{CL}}^{2}-{\omega }_{\mathrm{HD}}^{2}),\\ {{\rm{\Gamma }}}_{\mathrm{CL}} & \approx & \displaystyle \frac{1}{2{\omega }_{\mathrm{CL}}^{2}{\tilde{\tau }}_{\mathrm{coll}}}({\omega }_{\mathrm{CL}}^{2}-{\omega }_{\mathrm{HD}}^{2}).\end{eqnarray} \tag{ 6 }$$

In the transition regime, again provided that ${\rm{\Gamma }}\ll {\omega }_{r}$, the damping takes on the simple form [26]:

$$\begin{eqnarray}&&{\rm{\Gamma }}\approx \displaystyle \frac{{\tilde{\tau }}_{\mathrm{coll}}}{2}\displaystyle \frac{{\omega }_{\mathrm{CL}}^{2}-{\omega }_{\mathrm{HD}}^{2}}{1+{\omega }_{\mathrm{CL}}^{2}{\tilde{\tau }}_{\mathrm{coll}}^{2}}.\end{eqnarray} \tag{ 7 }$$

For traps where the spherical symmetry is broken, the monopole mode couples to the quadrupole modes. It is therefore collisionally damped at the HD crossover occurring between the HD and CL regimes. In section 3, the impact and strength of this coupling on the monopole mode is discussed.

Of the many schemes for magnetically trapping ultracold neutral atoms, the TOP trap [27] allows flexibility in tailoring the level of anisotropy present in the system. TOP traps utilize a fast rotating bias field oriented perpendicular to a quadrupole magnetic field. This avoids Majorana losses that would otherwise occur due to the zero point of the field being inside the cloud [28]. The quadrupole and fast rotating bias fields produce a magnetic field that has Cartesian components

$$\begin{eqnarray}&&\vec{B}=\left(\begin{array}{c}\displaystyle \frac{{B}_{z}}{2}x+{B}_{0}\mathrm{cos}({\rm{\Omega }}t)\\ \displaystyle \frac{{B}_{z}}{2}y+{B}_{0}\mathrm{sin}({\rm{\Omega }}t)\\ -{B}_{z}z\end{array}\right),\end{eqnarray} \tag{ 8 }$$

where B₀ is the field strength of the bias field, B_z is the strength of the quadrupole field. This bias field rotates at angular frequency Ω, which is typically much larger than the Larmor frequency. Near the minimum of the trap the cooled atoms experience a time averaged potential

$$\begin{eqnarray}&&{U}_{\mathrm{TOP}}(r,z)=\mu {B}_{0}+\displaystyle \frac{\mu {B}_{z}}{16{B}_{0}}({r}^{2}+8{z}^{2}),\end{eqnarray} \tag{ 9 }$$

where μ is the magnetic moment and $r=\sqrt{{x}^{2}+{y}^{2}}$.

In the high-field limit $(\mu {B}_{z}\gg {mg})$ the trap is anisotropic with aspect ratio ${\omega }_{z}/{\omega }_{r}=\sqrt{8}$. As the quadrupole field strength is decreased, the effects of gravity become important as the atoms in the trap sag away from the high-field limit trap minimum. The extent of the sag is characterized by the dimensionless quantity ${\rm{\Lambda }}={mg}/\mu {B}_{z}$ and shifted trap minimum ${z}_{0}=-\left({B}_{0}{\rm{\Lambda }}\right)/\left({B}_{z}\sqrt{1-{{\rm{\Lambda }}}^{2}}\right)$. The potential expanded about the sag position z₀ is

$$\begin{eqnarray}{U}_{\mathrm{sag}}(r,z) & = & \mu {B}_{0}\sqrt{1-{{\rm{\Lambda }}}^{2}}+\displaystyle \frac{\mu {B}_{z}^{2}}{16{B}_{0}}(1+{{\rm{\Lambda }}}^{2}){r}^{2}\sqrt{1-{{\rm{\Lambda }}}^{2}}\\ & & +\displaystyle \frac{\mu {B}_{z}^{2}}{2{B}_{0}}{z}^{2}{(1-{{\rm{\Lambda }}}^{2})}^{3/2}.\end{eqnarray} \tag{ 10 }$$

The ratio of the trap frequencies is

$$\begin{eqnarray}&&\displaystyle \frac{{\omega }_{z}}{{\omega }_{r}}=\sqrt{8\displaystyle \frac{1-{{\rm{\Lambda }}}^{2}}{1+{{\rm{\Lambda }}}^{2}}}\end{eqnarray} \tag{ 11 }$$

with ${\rm{\Lambda }}=7/9$ giving an isotropic trapping potential. For further details of the modified TOP trap used in the JILA experiment, we refer the reader to [29].

Utilizing gravity to symmetrize the trap results in trapping frequencies on the order of 10 Hz. In such a loose trap anharmonic corrections become important and must be included in calculations. The potential used in calculations involving the anharmonic corrections has the form

$$\begin{eqnarray}{U}_{\mathrm{ah}}(r,z) & = & \displaystyle \frac{m}{2}\left({\omega }_{z}^{2}{z}^{2}+{\omega }_{r}^{2}{r}^{2}\right)+\displaystyle \frac{m\alpha }{3}{z}^{3}+\displaystyle \frac{m\beta }{3}{r}^{2}z\\ & & +\displaystyle \frac{m\kappa }{4}{z}^{4}+\displaystyle \frac{m\delta }{4}{r}^{2}{z}^{2}+\displaystyle \frac{m\epsilon }{4}{r}^{4}\end{eqnarray} \tag{ 12 }$$

for general formulation of the cubic and quartic terms (see [29]). The values of the anharmonic coefficients in SI unites used in this paper were quoted from the JILA experiment [1] as $\alpha =5.58\times {10}^{6}$ ${{\rm{s}}}^{-2}$ ${{\rm{m}}}^{-1}$, $\beta =4.17\times {10}^{6}$ ${{\rm{s}}}^{-2}\;{{\rm{m}}}^{-1}$, $\kappa =5.94\times {10}^{9}$ ${{\rm{s}}}^{-2}\;{{\rm{m}}}^{-2}$, $\delta =1.716\times {10}^{9}$ ${{\rm{s}}}^{-2}\;{{\rm{m}}}^{-2}$, and $\epsilon =2.469\times {10}^{8}$ ${{\rm{s}}}^{-2}\;{{\rm{m}}}^{-2}$.

The monopole mode is driven experimentally by applying a sinusoidal variation in the trapping frequency over four periods of monopole oscillation. The strength of the TOP field is changed along with a vertical bias which fixes the minimum of the potential (see [29]). The net result is a roughly 25% increase in the mean cloud size.

By taking moments of the Boltzmann equation (1), a set of coupled equations of motion can be derived that describe the dynamical evolution of various collective modes [22]. For an arbitrary quantity $\chi ({\bf{r}},{\bf{v}})$, the equation of motion is

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle \chi \rangle }{{\rm{d}}t}-\langle {\bf{v}}\cdot {{\rm{\nabla }}}_{{\bf{r}}}\rangle -\left\langle\displaystyle \frac{{\bf{F}}}{m}\cdot {{\rm{\nabla }}}_{{\bf{v}}}\chi \right\rangle=\langle \chi {I}_{\mathrm{coll}}\rangle ,\end{eqnarray} \tag{ 13 }$$

where the average quantity is given by

$$\begin{eqnarray}&&\langle \chi \rangle =\displaystyle \frac{1}{N}\int {{\rm{d}}}^{3}r{{\rm{d}}}^{3}{vf}({\bf{r}},{\bf{v}},t)\chi ({\bf{r}},{\bf{v}}),\end{eqnarray} \tag{ 14 }$$

and the collisional average of this quantity is given by

$$\begin{eqnarray}&&\langle \chi {I}_{\mathrm{coll}}\rangle =\displaystyle \frac{1}{4N}\int {{\rm{d}}}^{3}r{{\rm{d}}}^{3}{{\bf{v}}}_{1}{\rm{\Delta }}\chi {I}_{\mathrm{coll}}[f],\end{eqnarray} \tag{ 15 }$$

where ${\rm{\Delta }}\chi ={\chi }_{1}+{\chi }_{2}-{\chi }_{{1}^{\prime }}-{\chi }_{{2}^{\prime }}$ with ${\chi }_{i}=\chi ({\bf{r}},{{\bf{v}}}_{i})$.

From equation (13), it is straightforward to form a closed set of coupled equations for the quadratic moments, recalling that moments of the form in equation (4) have no collisional contribution. For an isotropic trap, taking the moment $\langle {{\bf{r}}}^{2}\rangle$ yields a closed set of equations describing the evolution of the monopole mode:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {{\bf{r}}}^{2}\rangle }{{\rm{d}}t}-2\langle {\bf{v}}\cdot {\bf{r}}\rangle =0,\\ &&\displaystyle \frac{{\rm{d}}\langle {\bf{v}}\cdot {\bf{r}}\rangle }{{\rm{d}}t}-\langle {{\bf{v}}}^{2}\rangle +{\omega }_{0}^{2}\langle {{\bf{r}}}^{2}\rangle =0,\\ &&\displaystyle \frac{{\rm{d}}\langle {{\bf{v}}}^{2}\rangle }{{\rm{d}}t}+2{\omega }_{0}^{2}\langle {\bf{v}}\cdot {\bf{r}}\rangle =0.\end{eqnarray} \tag{ 16 }$$

Solutions to equation (16), where the moments oscillate around their equilibrium values as ${{\rm{e}}}^{{\rm{i}}\omega t}$, confirm Boltzmann's general result as previously discussed in section 2 for an undamped monopole mode at oscillation frequency $2{\omega }_{0}$.

One can also take $(l=2,m=0)$ quadrupole moments of equation (13) to derive a closed set of six coupled equations for an isotropic trap. However, for the quadrupole mode the moment, $\langle 3{v}_{z}^{2}-{{\bf{v}}}^{2}\rangle$ cannot be expressed in the form of equation (4), which leads to a nontrivial collisional contribution $\langle (3{v}_{z}^{2}-{{\bf{v}}}^{2}){I}_{\mathrm{coll}}\rangle$ in its equation of motion. For small amplitude oscillations about equilibrium, collisional contributions of this form can be approximately rewritten in the relaxation time approximation [22, 25, 30] as

$$\begin{eqnarray}&&\langle \chi {I}_{\mathrm{coll}}\rangle =-{\displaystyle \frac{\chi }{\tau }}_{\mathrm{coll}},\end{eqnarray} \tag{ 17 }$$

where the relaxation time ${\tau }_{\mathrm{coll}}$ is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{coll}}=\displaystyle \frac{5}{4{\gamma }_{\mathrm{coll}}},\end{eqnarray} \tag{ 18 }$$

and ${\gamma }_{\mathrm{coll}}=n(0){v}_{\mathrm{th}}\sigma /2$ is the average collision rate with ${v}_{\mathrm{th}}=\sqrt{8{k}_{{\rm{b}}}T/\pi m}$. As with the monopole mode, the set of coupled equations for the quadrupole mode can also be solved algebraically, giving frequencies of $2{\omega }_{0}$ and $\sqrt{2}{\omega }_{0}$ in the CL and HD regimes respectively. However, in the HD crossover the quadrupole mode damps as in equation (5) due to the nonvanishing collisional contribution.

The scaling ansatz method has been used for Bose gases above and below the transition temperature [31–34]. In the nondegenerate regime, the ansatz method has been used to evaluate the effects of the mean field interaction [32, 33] and to estimate the effects of anharmonic corrections to the trapping potential on the collective mode frequencies and damping rates [34]. Here this method is outlined.

When the thermal cloud is in statistical equilibrium under harmonic confinement, the Boltzmann equation for the equilibrium distribution ${f}_{0}({\bf{r}},{\bf{v}})$ is just the Poisson bracket:

$$\begin{eqnarray}&&\left(\displaystyle \sum _{i}^{3}{v}_{i}\displaystyle \frac{\partial {f}_{0}}{\partial {r}_{i}}-{\omega }_{i}^{2}{r}_{i}\displaystyle \frac{\partial {f}_{0}}{\partial {v}_{i}}\right)=0,\end{eqnarray} \tag{ 19 }$$

where the subscript denotes components of a Cartesian vector in 3D and the possibly differing frequency components in each direction. Moments of the equilibrium distribution are the same as those obtained in the previous subsection when the temporal derivatives are set to zero.

To describe time-dependent collective oscillations, a scaling ansatz on the form of the distribution function $f({\bf{r}},{\bf{v}},t)$ can be made:

$$\begin{eqnarray}&&f({\bf{r}},{\bf{v}},t)={\rm{\Gamma }}{f}_{0}({\bf{R}}({\bf{r}},{\bf{t}}),{\bf{V}}({\bf{r}},{\bf{v}},{\bf{t}})).\end{eqnarray} \tag{ 20 }$$

This ansatz utilizes the symmetry of the problem and the nature of the collective modes in the forms of the renormalized position, ${\bf{R}}({\bf{r}},{\bf{t}})$, and velocity, ${\bf{V}}({\bf{r}},{\bf{v}},{\bf{t}})$, with an additional factor Γ to enforce the normalization.

For the dipole mode, the form of the ansatz is ${R}_{i}={r}_{i}-{\eta }_{i}(t)$ and ${V}_{i}={v}_{i}-{\dot{\eta }}_{i}(t)$ with Γ equal to unity. The details of the motion are contained in the time-dependent vector of free parameters, $\vec{\eta }$. Such a scaling clearly describes a translation of the COM of the equilibrium cloud.

For the monopole and quadrupole mode, the ansatz mimics the form of Boltzmann's solution as discussed in section 2:

$$\begin{eqnarray}{R}_{i} & = & \displaystyle \frac{{r}_{i}}{{b}_{i}(t)},\\ {V}_{i} & = & \displaystyle \frac{1}{{\theta }_{i}^{1/2}}\left({v}_{i}-\displaystyle \frac{{\dot{b}}_{i}(t)}{{b}_{i}(t)}{r}_{i}\right),\\ {\rm{\Gamma }} & = & \displaystyle \frac{1}{\displaystyle \prod _{j}^{3}({b}_{j}(t){\theta }_{j}{(t)}^{1/2}),}\end{eqnarray} \tag{ 21 }$$

where there are two time-dependent vectors of free parameters $\vec{\theta }$ and $\vec{b}$. Such a scaling describes the stretching and compression of the equilibrium cloud in phase space along with a space-dependent and time-dependent translation in the local velocity, while maintaining a stationary COM.

Substituting equation (20) into equation (1) and following [32, 33], a Newton-like set of equations of motion for the free parameters can be derived

$$\begin{eqnarray}&&{\ddot{b}}_{i}+{\omega }_{i}^{2}{b}_{i}-{\omega }_{i}^{2}\displaystyle \frac{{\theta }_{i}}{{b}_{i}}=0\\ &&{\dot{\theta }}_{i}+2\displaystyle \frac{{\dot{b}}_{i}}{{b}_{i}}{\theta }_{i}=-\displaystyle \frac{1}{{\tau }_{\mathrm{coll}}}\left[{\theta }_{i}-\bar{\theta }\right],\end{eqnarray} \tag{ 22 }$$

where the quantity $\bar{\theta }={\displaystyle \sum }_{i}{\theta }_{i}/3$ is the average temperature. To obtain information about the collective oscillations, the free parameters are linearized about the equilibrium position for small amplitude oscillation, ${b}_{i}\approx 1$ and ${\theta }_{i}\approx 1$, and assumed to have time-dependence of the form ${{\rm{e}}}^{{\rm{i}}\omega t}$. For a harmonic trap, the results obtained agree with the moment method.

Slight anisotropies in the harmonic confinement and anharmonic corrections to the trapping potential equation (12) present possible sources of damping. Modeling these imperfections is discussed in this subsection using trapping data from the JILA experiment.

3.3.1. Damping due to trap anisotropies

Anisotropies in the confining potential lead to coupling between the monopole and quadrupole modes, which are damped in the HD crossover regime. The moment method can account for this coupling in a straightforward way given a potential of the form

$$\begin{eqnarray}&&U(x,y,z)=\displaystyle \frac{m{\omega }_{0}^{2}}{2}\left((1+\epsilon ){x}^{2}+(1-\epsilon ){y}^{2}+{\lambda }^{2}{z}^{2}\right),\end{eqnarray} \tag{ 23 }$$

where $\lambda ={\omega }_{z}/{\omega }_{0}$, ${\omega }_{0}\sqrt{1+\epsilon }={\omega }_{x}$, and ${\omega }_{0}\sqrt{1-\epsilon }={\omega }_{y}$.

Starting with $\langle {{\bf{r}}}^{2}\rangle$, the moment method yields a set of nine coupled equations for the monopole and quadrupole modes (see [35]). In figure 1, frequency measurements from the JILA experiment with various anisotropies were used to estimate the damping rate for small oscillations. For nonzero and λ different from unity, the monopole couples to multiple quadrupole modes and the shape of the curve departs from that expected for the typical transition region shape when $\epsilon =0$.

Zoom In Zoom Out Reset image size

The monopole damping measured in the JILA experiment is on the order of 0.1 s⁻¹ or larger. Therefore, the trap can roughly be treated as isotropic with $\lambda =1$ and $\epsilon =0$ for the lowest four curves in figure 1. Thus, the observed monopole damping must then be due to some other mechanism that is present even when the harmonic piece of the trapping potential is virtually isotropic.

3.3.2. Damping due to anharmonic corrections

Computing the monopole damping using the moment method for nonzero anharmonic coefficients generates an infinite set of coupled equations for successively higher order moments and so is not a viable method. However, the scaling ansatz method, which takes as an input moments of f₀, can provide an estimate for the effects of the anharmonic coefficients which enter through higher order moments of f₀. The usefulness of such an approach hinges on being able to neglect deformations of the equilibrium distribution due to the anharmonic corrections.

Using the scaling ansatz, equation (22) can be rewritten in a more general form for an arbitrary trapping force, F_i:

$$\begin{eqnarray}&&{\ddot{b}}_{i}-\displaystyle \frac{{\theta }_{i}}{{b}_{i}}\langle {v}_{i}^{2}{\rangle }_{0}-\displaystyle \frac{1}{m}\langle {F}_{i}({b}_{j}{r}_{j}){r}_{i}{\rangle }_{0}=0,\\ &&\dot{{\theta }_{i}}+2\displaystyle \frac{{\dot{b}}_{i}}{{b}_{i}}{\theta }_{i}=-\displaystyle \frac{1}{{\tau }_{\mathrm{coll}}}\left[{\theta }_{i}-\bar{\theta }\right],\end{eqnarray} \tag{ 24 }$$

where the 0 subscript indicates a moment as in equation (14) over the equilibrium distribution f₀, and ${b}_{j}{r}_{j}$ is a short hand notation for the set $\{{b}_{x}{r}_{x},{b}_{y}{r}_{y},{b}_{z}{r}_{z}\}$. From equation (24), the procedure for including even order corrections to the trapping potential is straightforward. However, to include odd order corrections, the equilibrium distribution, f₀, must be deformed. A simple deformation is

$$\begin{eqnarray}&&{f}_{0}^{\mathrm{ah}}({\bf{r}},{\bf{v}})=\left(1-\displaystyle \frac{m\beta }{3{k}_{{\rm{b}}}T}\left({x}^{2}+{y}^{2}\right)z-\displaystyle \frac{m\alpha }{3{k}_{{\rm{b}}}T}{z}^{3}\right){f}_{0}({\bf{r}},{\bf{v}}),\end{eqnarray} \tag{ 25 }$$

which is a perturbative expansion to first order in the odd order anharmonic corrections from equation (12). The deformation does not change the overall normalization nor the value of the even-order moments, and the added terms are collisional invariants of the form equation (4).

In section 5, the validity of this deformation compared to numerical and experimental results for the damping of the monopole mode is addressed. However, in the remainder of this section, the general solution of equation (24) is discussed for small amplitude oscillations, following the derivation in [34] and noting that the derivation is the same for ${f}_{0}^{\mathrm{ah}}$.

Equation (24) is solved in both the CL and HD regimes. Making the assumption that the monopole damping is much less than the collective mode frequency $({{\rm{\Gamma }}}_{{\rm{M}}}\ll {\omega }_{{\rm{M}}})$ leads to a generally applicable formula equation (7) for the damping. In the JILA experiment, ${{\rm{\Gamma }}}_{{\rm{M}}}/{\omega }_{{\rm{M}}}$ is on the order of 10⁻³ to 10⁻².

In the CL regime, ${\tau }_{\mathrm{coll}}\to \infty$, which gives a power law relation between the two scaling parameters: ${\theta }_{i}=\frac{1}{{b}_{i}^{2}}$. Using this power law relation, equation (24) collapses into three differential equations for each component of the parameter b_i. For small amplitude oscillations ${b}_{i}\approx 1$, and the substitution ${\eta }_{i}={b}_{i}-1$ linearizes the three differential equations for each component of ${\eta }_{i}$:

$$\begin{eqnarray}&&{\ddot{\eta }}_{i}+\displaystyle \frac{3}{m}\displaystyle \frac{\langle {U}_{i}{r}_{i}{\rangle }_{0}}{\langle {r}_{i}^{2}{\rangle }_{0}}{\eta }_{i}+\displaystyle \frac{1}{m}\displaystyle \sum _{j}\displaystyle \frac{\langle {U}_{{ij}}{r}_{i}{r}_{j}{\rangle }_{0}}{{\langle r}_{i}^{2}{\rangle }_{0}}{\eta }_{j}=0,\end{eqnarray} \tag{ 26 }$$

where ${U}_{i}=\partial U/\partial {r}_{i}$ and ${U}_{{ij}}={\partial }^{2}U/\partial {r}_{i}\partial {r}_{j}$. Equation (26) can be treated with matrix methods as in [34] to extract the collective mode frequencies.

In the HD regime, ${\tau }_{\mathrm{coll}}\to 0$, and thus the system is in local equilibrium everywhere. Furthermore, the components of the ${\theta }_{i}$'s must be equal to their average, ${\theta }_{i}=\bar{\theta }$. From the normalization of the ansatz then follows a power law relation between the parameters: ${\theta }_{i}=1/{\displaystyle \prod }_{j}{b}_{j}^{2/3}$. This implies that equation (24) collapses to a set of three differential equations for each component of the parameter b_i that couple when the power law relation for ${\theta }_{i}$ is substituted. Making the small amplitude assumption leads to the following equation:

$$\begin{eqnarray}&&\ddot{{\eta }_{i}}+\displaystyle \frac{1}{m}\displaystyle \frac{\langle {U}_{i}{r}_{i}{\rangle }_{0}}{\langle {r}_{i}^{2}{\rangle }_{0}}\left(\displaystyle \frac{5}{3}{\eta }_{i}+\displaystyle \frac{2}{3}\displaystyle \sum _{i\ne j}{\eta }_{j}\right)+\displaystyle \frac{1}{m}\displaystyle \sum _{j}\displaystyle \frac{\langle {U}_{{ij}}{r}_{i}{r}_{j}{\rangle }_{0}}{\langle {r}_{i}^{2}{\rangle }_{0}}=0.\end{eqnarray} \tag{ 27 }$$

As in the CL case, matrix methods can be used to extract the collective mode frequencies for the monopole mode and the (l = 2, m = 0) and (l = 2, m = 2) quadrupole modes.

The calculated mode frequencies have temperature dependence through higher-order moments coming from the anharmonic corrections. For a given temperature, the scaling ansatz theory gives the shift in the frequency ${\rm{\Delta }}{\omega }_{{\rm{M}}}$ and damping Γ from that anticipated by the harmonic limit ${\omega }_{{\rm{M}}}^{0}=2{\omega }_{0}$. The damping Γ is from collisional relaxation in the system.

The shift in the frequency ${\rm{\Delta }}{\omega }_{{\rm{M}}}={\omega }_{{\rm{M}}}^{0}-{\omega }_{{\rm{M}}}^{{\rm{P}}}$, measures the gap between the result in the harmonic limit ${\omega }_{{\rm{M}}}^{0}$ and the peak frequency component ${\omega }_{{\rm{M}}}^{{\rm{P}}}$ present in the cloud. Centered about the peak component is a spread of frequencies whose width $\delta {\omega }_{{\rm{M}}}$ is not predicted from the scaling ansatz method. To obtain full information about the width, one can perform a numerical simulation as discussed in section 4. The width provides a measure of the dephasing-induced damping ${{\rm{\Gamma }}}_{\mathrm{dephase}}$, which is responsible for the appearance of actual relaxation in the system and is separate from the collisional damping estimate from the scaling ansatz theory. In [34] it was argued that for a given width $\delta {\omega }_{{\rm{M}}}$, the dephasing time $\delta t$ is given approximately by $\delta {\omega }_{{\rm{M}}}\delta t\approx 2\pi$, and therefore ${{\rm{\Gamma }}}_{\mathrm{dephase}}\approx 2\pi /\delta t\approx \delta {\omega }_{{\rm{M}}}$. For an anisotropic trap, the frequency shift was used in [34] in lieu of the actual width to obtain dephasing-induced damping rates which agreed with experimental results within an order of magnitude. This approximation is reasonable only in the CL regime where the spread of frequencies in the cloud begins with the collective oscillation of the lowest energy atoms at ${\omega }_{{\rm{M}}}^{0}$. As will be made precise in section 5, when the collision rate increases, $\delta {\omega }_{{\rm{M}}}$ narrows about the peak value ${\omega }_{{\rm{M}}}^{{\rm{P}}}$ and decreases proportionately with ${{\rm{\Gamma }}}_{\mathrm{dephase}}$. The shift ${\rm{\Delta }}{\omega }_{{\rm{M}}}$ becomes then a poor estimate for the width $\delta {\omega }_{{\rm{M}}}$. As ${{\rm{\Gamma }}}_{\mathrm{dephase}}$ decreases a transition regime exists where dephasing effects in the cloud are destroyed through collisions, and collisional damping becomes the dominant effect. Such an effect can account for anomalous damping of the monopole mode in regimes where ${\omega }_{{\rm{M}}}{\tau }_{\mathrm{coll}}\gg 1$, whereas collisional damping dominates when ${\omega }_{{\rm{M}}}{\tau }_{\mathrm{coll}}\approx 1$. We refer to the regime where $\delta {\omega }_{{\rm{M}}}{\tau }_{\mathrm{coll}}\approx 1$ as the dephasing crossover, and in section 5 this crossover is studied in detail with the aid of numerical simulation.

The collision algorithm, which is a procedure for Monte Carlo sampling the collision integral [41], follows that described in [36]. After binning the particles in space, random pairs in each bin are selected and their collision probability is calculated. If a collision is successful, the particle velocities are updated according to the differential cross-section for an s-wave collision.

For each set of simulation parameters we check that the equilibrium collision rate averaged over the entire cloud matches with the analytic result:

$$\begin{eqnarray}&&{\gamma }_{\mathrm{eq}}={N}_{\mathrm{th}}\displaystyle \frac{\sigma {\omega }_{0}^{3}\;m}{2{\pi }^{2}{k}_{{\rm{b}}}T},\end{eqnarray} \tag{ 29 }$$

to $\lt 2\%$. The velocity independent cross-section used is that for bosons $\sigma =8\pi {a}^{2}$.

To mimic the experimental monopole drive scheme, the trap frequency, ${\omega }_{0}$, is modulated over four periods of monopole oscillation

$$\begin{eqnarray}{\omega }_{0}(t)=\left\{\begin{array}{lr}{\omega }_{0}\left(1+A\mathrm{sin}(2{\omega }_{0}t)\right) & t\leqslant 4\pi /{\omega }_{0}\\ {\omega }_{0} & t\gt 4\pi /{\omega }_{0}\end{array}\right\},\end{eqnarray} \tag{ 30 }$$

where A is a unitless measure of the strength of the drive. During the drive, the anharmonic corrections are neglected, and are switched on when the next oscillation minimum occurs after the drive is turned off. Figure 2 illustrates the increase of the monopole amplitude during the drive and the effect of turning on the anharmonic shifts on the total cloud energy. For A = 0.15 and a starting temperature of 152 nK, the percent increase in the total energy when the anharmonic corrections are switched on is on the order of $0.001\%$. Neglecting the anharmonic terms during the drive effectively captures the experimental scenario without including the modulation of the bias, quadrupole, and all of the shimming fields. Results are presented then over a range of drive strengths that produce an increase in the mean cloud size at the end of the drive phase on the order with that experimentally measured. Driving the cloud ideally emphasizes the impact of the anharmonic shifts on the monopole mode and eliminates any residual driving of the quadrupole and dipole modes which are not the subject of the current study. For the full details of the drive, we refer readers to [1].

Zoom In Zoom Out Reset image size

During the drive there is also a shift in the average temperature and full width at half maximum (FWHM) of the oscillation due to an increase in the energy of the cloud from the work done on it by pumping of the trap.

In the CL regime, the many-body dynamics of the thermal cloud are dictated by the single-particle trajectories. Here, the system behaves in analogy with a simple pendulum. For small amplitude oscillations the pendulum executes simple-harmonic motion. However, as the amplitude grows, the small-angle approximation breaks down and anharmonic corrections become important. The pendulum traces out a repeating trajectory but with an energy-dependent period.

The procedure for replacing anharmonic corrections by an energy-dependent harmonic potential in 1D is well-known [26]. The renormalized trapping potential for small amplitude motion along the z-axis is:

$$\begin{eqnarray}&&U(0,0,z)\approx \displaystyle \frac{m{\tilde{\omega }}_{z}^{2}(E)}{2}{z}^{2},\\ &&{\tilde{\omega }}_{z}={\omega }_{z}(1+{\xi }_{z}E),\\ &&{\xi }_{z}=-\displaystyle \frac{5{\alpha }^{2}}{6\;{m}^{3}{\omega }_{0}^{6}}+\displaystyle \frac{3\kappa }{4\;{m}^{2}{\omega }_{0}^{4}}.\end{eqnarray} \tag{ 31 }$$

In figure 3(a), atoms in the anharmonic trap have been binned according to their energy and the collective motion in each bin averaged to obtain the monopole period in each bin. For the first several bins the period of collective oscillation increases linearly with the bin energy. For higher energy bins the scaling is still monotonic; however, the deviation increases as the difference in the linear correction between the different axes becomes more apparent. At sufficiently high binning energy, the linear correction equation (31) breaks down as higher-order terms become important, and the atom number in each bin decreases, spoiling the appearance of an undamped collective mode.

Zoom In Zoom Out Reset image size

In the absence of collisions, the population of atoms in each bin is static over the entire simulation. The statistical weight of each frequency component, as shown in the normalized histogram of figure 3(b), is then also static. The determining characteristics of the frequency distribution are its width and shift of the peak from ${\omega }_{{\rm{M}}}^{0}=2{\omega }_{0}$. The theoretical prediction of ${\rm{\Delta }}{\omega }_{{\rm{M}}}$ from the scaling ansatz using ${f}_{0}^{\mathrm{ah}}$ is for the shift of the peak (grey region of figure 3(b)) and agrees with the numerical result. However, the width, denoted $\delta {\omega }_{{\rm{M}}}$, of the frequency distribution determines the relative importance of dephasing effects. If the frequency distribution is shifted and sharply peaked, dephasing induced damping is minimal, mimicking the large amplitude oscillation of a single pendulum. Otherwise, if the distribution is broad, as in figure 3(a), the net sum of many 'pendulums' oscillating at shifted frequencies is dephasing induced damping through interference. The scaling ansatz theory also gives a prediction for the damping due to collisional relaxation, and in the CL regime, this prediction is vanishingly small as expected (see figure 5).

We now compare directly to the small cloud data from the JILA experiment. Starting with the temperature and atom number quoted in the experiment, the drive strength is varied and the amplitude of the resulting oscillation is fit to a simple exponential decay to compare directly to the fitting function used in the experiment. As shown in figure 4(a), a drive strength of A = 0.15 matches the experimental data, suggesting that the anomalous damping seen in the experiment for the small cloud data is mainly due to dephasing effects. Such effects depend on the settled temperature of the cloud as shown in figure 4(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although the width of the frequency distribution determines the strength of the dephasing induced damping, in the CL regime where the distribution is broad, the frequency shift ${\rm{\Delta }}{\omega }_{{\rm{M}}}$ and the width scale proportionately (see discussion at the end of section 3.3.2) and can be used as a fitting function, looking then for a scaling fit of the form

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{dephase}}(T)={\rm{\Delta }}{\omega }_{{\rm{M}}}(\zeta T),\end{eqnarray} \tag{ 32 }$$

where ζ is a free parameter independent of the temperature. For the experimental trap parameters, $\zeta =0.4$ from a fit to the data shown in figure 4(b).

The scaling ansatz result using f₀ instead of ${f}_{0}^{\mathrm{ah}}$ predicts a decrease of the period with increasing bin energy, which disagrees directly with the numerical result and equation (31). Additionally, from the red and blue lines in figure 4(b), as the settled temperature increases the deformed gaussian predicts that the peak frequency decreases (${\rm{\Delta }}{\omega }_{{\rm{M}}}(T)\gt 0$).

The CL regime picture of a collection of uncoupled monopole modes oscillating at shifted frequencies begins to break down when the dephasing period ${\tau }_{\mathrm{dph}}=2\pi /\delta {\omega }_{{\rm{M}}}$ exceeds ${\tau }_{\mathrm{coll}}$. This is typically in a different collisional regime than the HD crossover where ${\omega }_{0}{\tau }_{\mathrm{coll}}\approx 1$, and an average atom suffers a collision on a timescale faster than the oscillation period. Therefore, as the collision rate increases there are two important regimes: the dephasing crossover and the HD crossover

$$\begin{eqnarray}&&\delta {\omega }_{{\rm{M}}}{\tau }_{\mathrm{coll}}\approx 1(\mathrm{Dephasing}\ \mathrm{crossover}),\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\omega }_{0}{\tau }_{\mathrm{coll}}\approx 1(\mathrm{Hydrodynamic}\ \mathrm{crossover}).\end{eqnarray} \tag{ 34 }$$

In the remainder of this subsection an individual experimental run from the large cloud data set at 60.6 nK with $3.463\times {10}^{5}$ atoms with measured trap frequencies $({f}_{x}=9.036\;\mathrm{Hz},{f}_{y}=9.034\;\mathrm{Hz},{f}_{z}=9.034\;\mathrm{Hz})$ is analyzed. From the criteria equation (34) can be obtained crude estimates for the dephasing crossover ${\gamma }_{\mathrm{coll}}\approx {10}^{-1}-1$ s⁻¹ and the HD crossover ${\gamma }_{\mathrm{coll}}\approx {10}^{2}$ s⁻¹. This data set corresponds to the green curve in figure 1 for a quoted collision rate ${\gamma }_{\mathrm{coll}}=8.88$ s⁻¹ which lies between the two regimes. The maximum damping rate due to anisotropies from figure 1 is on the order of 10⁻³ s⁻¹ compared to the experimentally measured rate 0.14 ± 0.02 s⁻¹. Therefore, the harmonic part of the trapping potential is treated as isotropic.

To quantify the effect of the anharmonic corrections in the crossover regime, numerical simulations with A = 0.05 over a range of collision rates using the trapping frequency ${\nu }_{0}=9.035\;\mathrm{Hz}$ and initial temperature 60 nK were performed. Figure 5 contains the results along with the prediction for the collisional damping from the scaling ansatz theory using ${f}_{0}^{\mathrm{ah}}$. The scaling ansatz prediction is plotted for several temperatures, illustrating the effect of the increase in temperature of the cloud post-drive on the damping.

The estimates for the HD and dephasing crossover regimes also agree with the qualitative structure of figure 5. The dephasing crossover is marked by a decrease in the dephasing induced damping; whereas the HD crossover is marked by an increase in the collisional damping. The region in between the two crossovers is characterized by a local minimum in the damping rate. Comparing the 65 nK damping rate from figure 4 with the 60.6 nK damping rate from this subsection draws experimental support for this result. That the CL regime is characterized by higher damping than the dephasing crossover is a common feature of all of the large cloud data compared for similar temperatures to the small data [1]. The intermediate regime between dephasing and HD crossovers then provides an experimental window where damping from trap anharmonicities can be minimized.

Whereas the physics of dephasing in the CL regime is analogous to the uncoupled oscillation of a collection of pendulums with different oscillation energies, in the dephasing crossover the pendulums begin to couple, and, as illustrated in figure 6, the width of the spread in frequency components narrows. The coupling effectively synchronizes the different pendulums, which leads to a local minima in the damping rate as seen in figure 5. As the coupling increases, the system nears the HD crossover, and the spread in the frequencies is minimal. The damping is then due to the appearance of a temperature dependent anisotropic trap, where the monopole mode damps naturally through coupling with the quadrupole modes whose moments are not collisionally invariant.

Zoom In Zoom Out Reset image size