Collective Dynamics and Atom Loss in Bright Soliton Matter Waves

TL;DR

This paper models the dynamics of bright solitons in cold gases, revealing how excitations affect atom loss and enabling precise calculations of the critical atom number for stability.

Contribution

It introduces a variational approach that decouples radial and axial modes, providing new analytical tools for understanding soliton stability and atom loss.

Findings

Atom loss is insensitive to excitation degree.

The variational method accurately predicts the critical atom number.

Decoupling of modes simplifies the analysis of soliton dynamics.

Abstract

Motivated by recent experiments, we model the dynamics of bright solitons formed by cold gases in quasi-1D traps. A dynamical variational ansatz captures the far-from equilibrium excitations of these solitons. Due to a separation of scales, the radial and axial modes decouple, allowing for closed-form approximations for the dynamics. We explore how soliton dynamics influence atom loss, and find that the time-averaged loss is largely insensitive to the degree of excitation. The variational approach enables us to perform high precision calculations of the critical atom number (ie. the maximum number of atoms that can exist in a single soliton before the attractive forces overwhelm quantum pressure, leading to collapse).

Equations

$$\begin{gathered} i\hslash {\partial }_{\tilde{t}}\tilde{\psi }=-\frac{{\hslash }^2}{2m}{\tilde{\nabla }}^2\tilde{\psi }+\frac{m}{2}({{\tilde{\omega }}_r}^2{\tilde{r}}^2+{{\tilde{\omega }}_z}^2{\tilde{z}}^2)\tilde{\psi } \\ +\tilde{g}|\tilde{\psi }{|}^2\tilde{\psi }-i{\tilde{\Gamma }}_3|\tilde{\psi }{|}^4\tilde{\psi }, \end{gathered}$$

$$E_{\mathrm{gauss}}=\frac{N}{2}\left(\frac{1}{{\ell }_r^2}+\frac{1}{2{\ell }_z^2}+{\ell }_r^2+\frac{gN}{{\ell }_r^2{\ell }_z}\right)$$

$$\begin{gathered} \psi (r,z,t)=\frac{A(r,z,t)}{{\pi }^{\frac{3}{4}}{\ell }_z(t{)}^{\frac{1}{2}}{\ell }_r(t)}\exp [-\frac{r^2}{2{\ell }_r(t{)}^2}-\frac{z^2}{{\ell }_z(t{)}^2}] \\ \exp [i{r}^2{\varphi }_r+i{z}^2{\varphi }_z], \end{gathered}$$

$$A$$

$$N=\int d^{3}x|\psi {|}^2=\sum _{l=0}^n\sum _{m=-\frac{l}{2}}^{\frac{l}{2}}|c_{lm}{|}^2.$$

$$\frac{\delta S}{\delta {\psi }^{\ast }}=-i{\Gamma }_3|\psi {|}^4\psi$$

$$S$$

$$\mathcal{H}$$

$$\frac{\delta S}{\delta c_{lm}^{\ast }}=\int d^{3}x\frac{\delta S}{\delta {\psi }^{\ast }}\frac{\partial {\psi }^{\ast }}{\partial c_{lm}^{\ast }}.$$

$$\frac{\delta S}{\delta c_{lm}^{\ast }}=\int d^{3}x(-i{\Gamma }_3|\psi {|}^4\psi )\frac{\partial {\psi }^{\ast }}{\partial c_{lm}^{\ast }}.$$

$$\frac{\delta S}{\delta {\ell }_{r,z}}=\int d^{3}x(-i{\Gamma }_3|\psi {|}^4\psi )\frac{\partial {\psi }^{\ast }}{\partial {\ell }_{r,z}}+h.c.$$

$$S$$

$$+\frac{N{\omega }_z^2{\ell }_z^2}{4}+\frac{1}{2}N{\ell }_z^2(2{\varphi }_z^2+\dot{{\varphi }_z})+\frac{g{N}^2}{{\ell }_r^2{\ell }_z}+N\dot{\chi })$$

$${\varphi }_r-\frac{\dot{{\ell }_r}}{2{\ell }_r}$$

$${\varphi }_z-\frac{\dot{{\ell }_z}}{2{\ell }_z}$$

$$2{\ell }_r(\dot{{\varphi }_r}+2{\varphi }_r^2)+{\ell }_r-\frac{1}{{\ell }_r^3}-\frac{2gN}{{\ell }_r^3{\ell }_z}$$

$$2{\ell }_z(\dot{{\varphi }_z}+2{\varphi }_z^2)-\frac{1}{{\ell }_z^3}-\frac{2gN}{{\ell }_r^2{\ell }_z^2}$$

$$\dot{N}+\frac{2{\Gamma }_3{N}^3}{3\sqrt{3}{\pi }^3{\ell }_r^4{\ell }_z^2}$$

$$\dot{\chi }+\frac{1}{{\ell }_r^2}+\frac{1}{2{\ell }_z^2}+\frac{7gN}{2{\ell }_r^2{\ell }_z}$$

$$\ddot{{\ell }_r}+{\ell }_r-\left(1+\frac{2gN}{{\ell }_z}\right)\frac{1}{{\ell }_r^3}=0.$$

$$\ddot{{\ell }_z}-\frac{2gN}{{\ell }_z^2{\ell }_r^2}-\frac{1}{{\ell }_z^3}=0.$$

$${\ell }_r^2=$$

$${\ell }_z(\eta )$$

$$t(\eta )$$

$$\ddot{r}+\frac{1}{m}\frac{dV(r)}{dr}-\frac{L^2}{m^2{r}^3}=0$$

$$\vec{r}(t)=a\hat{x}\cos t+\frac{L}{a}\hat{y}\sin t.$$

$${\ell }_z(\eta )$$

$$t(\eta )$$

$$⟨{\ell }_r^{-2}⟩=\frac{1}{T_r}{\int }_0^{T_r}\frac{dt}{{\ell }_r^2}=\frac{1}{\sqrt{1+\frac{2gN}{{\ell }_z}}}\approx 1,$$

$$T_z⟨{\ell }_z^{-2}⟩={\int }_0^{T_z}\frac{dt}{{\ell }_z^2}=2\pi ,$$

$$\dot{N}\approx N^5{E}_r(1-{ϵ}^2{)}^{\frac{3}{2}}8{\Gamma }_3{g}^2/3\sqrt{3}{\pi }^3$$

$$N(t)\approx \frac{N(0)}{{\left(1+\frac{4N(0{)}^4(2g{)}^2(1-{ϵ}_z^2{)}^{\frac{3}{2}}{E}_r}{3\sqrt{3}{\pi }^3}2{\Gamma }_{3}t\right)}^{\frac{1}{4}}}.$$

Full text

I Introduction

Solitons are localized excitations of nonlinear media which robustly maintain their shape kasman2010 ; drazin1989 . Due to their robustness they play key roles in settings as disparate as weather patterns, optical fibers, and quantum degenerate atoms remoissenet1999 . The latter has provided a particularly large degree of control, enabling the quantitative investigation of a wide range of soliton phenomena abdullaev2008 . For example, cold atom solitons can collapse and explode in a “Bosenova” cornish2006 , form mesoscopic Bell states gertjerenken2013 , and improve the noise floor of a matter-wave interferometer polo2013 . Here we model the far from equilibrium dynamics of a single bright soliton formed from an atomic gas in a highly anisotropic trap.

This study is motivated by experiments at Rice University where a quantum degenerate gas of bosonic lithium atoms is confined by a quasi-1D trap hulet2017 . The interactions are suddenly changed from repulsive to attractive, which makes the cloud unstable. This modulational instability causes the cloud to break up into a train of solitons. The solitons bounce off one another and undergo large shape deformations. Here we model the dynamics of the shape deformations of a single soliton in this train. We use a time dependent variational method to derive equations of motion for parameters which describe the soliton’s shape. We furthermore make a series of approximations which yield closed-form expressions for the length and thickness of the soliton. For much of the experimentally relevant parameter range, we find that the closed-form expressions provide a good approximation to the dynamics.

Due to three-body processes, atoms are continually scattered out of the solitons and lost from the trap. Interestingly, the experimental loss rate appears to be orders-of-magnitude faster than would be expected from models of static solitons. One possible hypothesis is that the extra loss is due to the large-scale oscillations of the solitons. We model this process, exploring how atom loss is influenced by the dynamics. In the regime of large oscillations, collisions predominantly occur when the cloud is most compressed. We quantify these episodic losses and characterize their dependence on the amplitude of the oscillations. Radial oscillations slightly increase the loss rate, while axial oscillations slightly decrease it. Since the size of the effect is relatively small, we conclude that this hypothesis is not able to account for the greatly enhanced loss observed in the experiment.

The nonlinear dynamics literature contains some discussion of dynamical solitons – particularly “breathers,” which are a class of time dependent solutions to equations such as the nonlinear Schrodinger equation. Breathers are important for understanding phenomena such as rogue waves, and have some connections to the physics explored here.

We use two distinct approaches to modeling the soliton dynamics. The first, and most conventional, is to discretize space and numerically integrate the Gross-Pitaevskii equation. The second is to write down an ansatz wavefunction and derive a set of coupled ordinary differential equations for the parameters from a time-dependent variational principle. The second approach is much more numerically efficient than the first, as it identifies the most important degrees of freedom. We systematically increase the number of parameters in the ansatz. The six-parameter Gaussian ansatz captures the qualitative features, while the sixteen-parameter ansatz is more accurate than our spatial discretization approach (using a grid-spacing appropriate for calculations on a laptop computer). Using a separation-of-scales, we simplify the equations for the Gaussian ansatz, and make a formal mapping onto two classic problems, simple harmonic and planetary motion. We thereby produce approximate closed-form solutions to the equations and find invariants which have physical significance.

Throughout this paper we model a single soliton. By contrast, the most recent experiments at Rice University study a train of 10-20 solitons. Although generalizing our results to a soliton train is straightforward, it is beyond the scope of this study. Qualitative features, such as the dependence of loss-rate on excitation amplitude, are nonetheless directly relevant to these experiments hulet2017 .

II Model

The cold gas experiments of hulet2017 are well-modeled by the three-dimensional, axially symmetric Gross-Pitaevskii equation with a harmonic trap and three-body loss pethick2008 ,

[TABLE]

where the tilde symbols are used to denote quantities with physical units. Planck’s constant is $\hbar$, and the atomic mass is $m$. The axial trap frequency $\tilde{\omega}_{z}$ is much smaller than the radial frequency $\tilde{\omega}_{r}$: in the experiments of hulet2017 , $\frac{\tilde{\omega}_{r}}{2\pi}=346$ Hz, and $\frac{\tilde{\omega}_{z}}{2\pi}=7.4$ Hz. For the soliton dynamics that we study, the axial trap plays a very small role. Thus, we set $\tilde{\omega}_{z}=0$ for most of the discussion. The atoms interact via a short-range potential with scattering length $a_{s}$, yielding a coupling constant $\tilde{g}=4\pi\hbar^{2}a_{s}/m$. In the experiments, $a_{s}$ is tuned via a Feshbach resonance chin2010 from approximately $3a_{0}$ to various negative values of order $-a_{0}$, where $a_{0}$ is the Bohr radius. Atom losses from inelastic three-body collisions are modeled by the term proportional to $\Gamma_{3}$, where $\Gamma_{3}=L_{3}/12\sim 10^{-29}$ cm6/s. Physically, the atomic density is $n(r)=|\psi(r)|^{2}$. We adimensionalize this equation, scaling length as $x=\sqrt{\frac{m\tilde{\omega}_{r}}{\hbar}}\tilde{x},$ time as $t=\tilde{\omega}_{r}\tilde{t},$ frequency as $\omega=\frac{\tilde{\omega}}{\tilde{\omega}_{r}},$ and the wavefunction as $\psi=(\frac{\hbar}{m\tilde{\omega}_{r}})^{\frac{3}{2}}\tilde{\psi}$. The dimensionless coupling constants are $g=4\sqrt{2\pi^{3}}\hbar\tilde{\omega}_{r}(\frac{\hbar}{m\tilde{\omega}_{r}})^{\frac{3}{2}}\tilde{g}$ and $\Gamma_{3}=m^{3}\tilde{\omega}_{r}^{2}\hbar^{-3}\tilde{\Gamma}_{3}$

We study systems with attractive interactions: $a_{s}<0$. In the absence of loss, Eq. (1) will have a metastable soliton solution as long as $N<N_{c}=k_{c}\sqrt{\frac{\hbar}{m\tilde{\omega}_{r}}}\frac{1}{|a_{s}|}$ perezgarcia1998 ; parker2007 . In Sec. III, we calculate that $k_{c}=0.677986(2)$ when $\omega_{z}=0$. As explained in billam2012 , this equilibrium physics can be understood through a Gaussian ansatz $\psi\propto\exp(-r^{2}/(2\ell_{r}^{2})-z^{2}/(2\ell_{z}^{2}))$. The energy as a function of $\ell_{r}$ and $\ell_{z}$ is

[TABLE]

An example energy landscape is shown in Fig. 2 for a particular $N<N_{c}$. It consists of two basins, separated by a saddle. In the region where both $\ell_{z}$ and $\ell_{r}$ are small, the attractive interactions dominate, causing the cloud to shrink rapidly. Consequently, the density becomes very high and rapid three-body recombination leads to large atom losses wieman2001 ; hulet2013 . This process of shrinking and evaporating is often referred to as a collapse. When $N<N_{c}$ the energy landscape also contains a local minimum at finite $\ell_{z}$ and $\ell_{r}$ where the soliton is metastable. The energy landscape is highly anisotropic around this minimum: the curvature of the energy landscape in the $\ell_{z}$ direction is much smaller than the curvature in the $\ell_{r}$ direction. In the limit $N\ll N_{c}$, the low energy dynamics become purely axial, and Eq. (1) reduces to a 1D nonlinear Schrodinger equation billam2012 .

When $N>N_{c}$ the energy landscape is qualitatively different, with no stationary points at nonzero $\ell_{z}$ and $\ell_{r}$. In this limit a single soliton will always collapse. As is illustrated by the experiment, the only way to have $N>N_{c}$ is to have multiple solitons.

To model the dynamics of a soliton, we extend the ansatz to a completely general time-dependent form,

[TABLE]

where again $\ell_{r}(t)$ is the radial size, and $\ell_{z}(t)$ is the axial size. The phases proportional to $r^{2}$ and $z^{2}$ are related to currents.

The factor $A(r,z,t)$ determines the number of particles in the soliton and accounts for all deviations from the Gaussian shape. We impose cylindrical and inversion symmetry. Consequently, we expand $A$ in a set of even-degree orthogonal polynomials, up to total degree $2n$:

[TABLE]

where $H_{k}(z)$ is the $k^{th}$ Hermite polynomial and $L_{k}(r)$ is the $k^{th}$ Laguerre polynomial. The factors are chosen so that

[TABLE]

In the limit $n\to\infty$, the exact solution to Eq. (1) can be written in the form of the ansatz in Eq. (4). By increasing $n$, we can increase the accuracy of our approximation, at the expense of increasing the number of coupled equations which need to be solved.

To derive the equations of motion for our variational parameters we write the Gross-Pitaevskii equation, Eq. (1), as

[TABLE]

where the action $S$ is a functional of $\psi$,

[TABLE]

Note that the dissipative term cannot be incorporated into the action.

Taking the variational derivative of Eq. (7), we see that

[TABLE]

We then use the equations of motion, Eq. (6), to arrive at our equations for the variational parameters,

[TABLE]

Similarly we have two equations associated with $\ell_{r}$ and $\ell_{z}$,

[TABLE]

We use a computer algebra system to perform the integrals in Eqs. (10)-(11), to arrive at a coupled set of equations for the variational parameters. Section IV.1 gives explicit examples in the Gaussian case.

In addition to the time-dependent variational approach, we analyze Eq. (1) by discretizing time and space. We then use a split-step algorithm to numerically calculate the dynamics of a soliton javanainen2006 . The results of the discretization approach agree with the results of the variational approach. The comparison is discussed in section IV.2.

III Critical Atom Number

As a first application of our variational wavefunction, we calculate $N_{c}$, the maximum number of atoms which can be found in a static condensate. A soliton with more atoms than $N_{c}$ is unstable and will collapse. As already explained, the critical atom number can be expressed as $N_{c}=k_{c}\sqrt{\frac{\hbar}{m\tilde{\omega}_{r}}}\frac{1}{|a_{s}|}$, where $k_{c}$ is a dimensionless quantity that only depends on the trap geometry perezgarcia1998 ; parker2007 . Much theoretical perezgarcia1998 ; parker2007 ; carr2002 ; gammal2002 ; salasnich2002 ; ruprecht1995 and experimental roberts2001 ; claussen2003 effort has been devoted to calculating $k_{c}$. In particular, Parker et al. parker2007 calculated that $k_{c}=0.675\pm 0.005$ for a cylindrical harmonic trap ($\omega_{z}=0$). They calculated this result with imaginary time evolution of the discretized Gross-Pittaevskii equation, a numerically intensive method. Their result agrees with similar analysis by Gammal et al. and Salasnich et al.gammal2002 ; salasnich2002 .

We calculate the same result using our variational wavefunction, which is a polynomial of degree $2n$ times a Gaussian. Since this is an equilibrium calculation, we can take the coefficients of the polynomials $c_{kj}$ to be real. The ansatz will have $N_{p}=(n+1)(n+2)/2+2$ parameters. We define $\bar{k}_{n}$ as the largest value of $k=N|a_{s}|\sqrt{m\tilde{\omega}_{r}/\hbar}$ for which our ansatz yields a soliton solution. Thus, the best estimate for $k_{c}$ from a given ansatz is $\bar{k}_{n}$.

To find $\bar{k}_{n}$, we first choose a modest $k$, which we know is below $\bar{k}_{n}$. We then optimize the parameters of our ansatz by minimizing the energy, using a conjugated gradient method. Subsequently, we slightly increase $k$ and repeat the minimization, using our previous solution as the starting point. We continue in this manner until a local minimum cannot be found – at that point we know we have $k\approx\bar{k}$. We further refine our estimate of $\bar{k}_{n}$ by using sequentially smaller step sizes.

The resulting $\bar{k}_{n}$ values are shown in Fig. 3. Remarkably, we find nearly identical results for $\bar{k}_{n}$ when $n=2m$ and $n=2m+1$. The even-$n$ data is well approximated by a power law, $\bar{k}_{n}=k_{c}+a(n-b)^{-c}$. We fit the parameters $a,b$, and $c$ in log space, by minimizing $\chi^{2}=\sum_{n}\log[\bar{k}_{n-2}-\bar{k}_{n}]-\log[a(n-2-b)^{-c}-a(n-b)^{-c}]$. The fit is illustrated by the inset of Fig. 3. We then extract $k_{c}$, by requiring that our ansatz pass through the final data point. The uncertainty is estimated by repeating the fit with one fewer data point, which gives an uncertainty in $k_{c}$ of $1.6\times 10^{-6}$. We further confirmed the extrapolation by fitting the odd-$n$ data, and by considering a number of alternative functional forms. In all cases the results differed by numbers of order $10^{-7}$ or smaller. We conclude that $k_{c}=0.677986(2)$, corresponding to the horizontal gray line in Fig. 3.

A comparison of the value of $N_{c}$ for different values of $n$ gives a good estimate for the accuracy of each ansatz. From the data, the Gaussian ansatz has a 15% error, the n=4 ansatz has a 0.4% error, and the n=8 ansatz has a 0.02% error.

IV Dynamics

In Sec. IV.1 we analyze the equations of motion for the Gaussian ansatz. In particular, in Sec. IV.1.1 through IV.1.4 we derive approximate analytic solutions to the equations of motion, and in Sec. IV.1.5 we describe their limitations. In section IV.2, we go beyond the gaussian approximation, calculating dynamics for the ansatz with larger $n$, and compare those results to a numerical integration of the Gross Pittaevskii equation.

IV.1 Gaussian Ansatz

The most important parameters describing a soliton are the dimensionless width $\ell_{r}$, length $\ell_{z}$, and particle number $|c_{00}|^{2}=N$. In this section we analyze the Gaussian ansatz, given by Eq. (3) with $n=0$. The Gaussian ansatz has six parameters: the above three variables, the global phase of the soliton, $\arg(c_{00})=\chi$, and the variables $\phi_{r}$ and $\phi_{z}$ which account for the probability currents. As already explained, we work in units set by the radial harmonic confinement. In this case the action can be written as

[TABLE]

Taking the variational derivatives as in Eqs. (10,11) then yields

[TABLE]

IV.1.1 Decoupling of timescales

Because of the quasi-1D geometry of the soliton, the timescale of the oscillations in the radial direction is much smaller than that of the oscillations in the axial direction. Both oscillation timescales are much smaller than the timescale of the atom losses. The separation of these timescales enables an adiabatic decoupling approach hinch1991 .

First, the terms on the right side of Eqs. (13)-(18) are insignificant and will be neglected. With this approximation in place, The phases, $\phi_{r},\phi_{z}$, are related to the lengths $\ell_{r},\ell_{z}$ by $\phi=\dot{\ell}/(2\ell)$, and consequently $2\ell(\dot{\phi}+2\phi^{2})=\ddot{\ell}$. This observation allows us to combine Eq. (13) and Eq. (15) and combine Eq. (14) and Eq. (16) into second order differential equations with no first-derivative terms:

[TABLE]

[TABLE]

Now we can take advantage of the different time-scales in the problem. Since $\ell_{z}$ and $N$ are slowly varying compared to $\ell_{r}$, we first solve Eq. (19) with $\ell_{z}$ and $N$ fixed. We then solve Eq. (20), replacing $\ell_{r}^{-2}$ with its time average, $\left\langle\ell_{r}^{-2}\right\rangle$. As we will soon see, $\left\langle\ell_{r}^{-2}\right\rangle=1/\sqrt{1+2gN/\ell_{z}}\approx 1$, is not only slowly varying, but also nearly constant. We will then solve Eq. (20) when $gN/\ell_{z}\ll 1$, and substitute those solutions into Eq. (17) to find $N(t)$.

As we explain in section IV.1.2, one can map these equations onto two iconic problems in classical mechanics to arrive at analytic solutions. For $\ell_{r}$, the solution is the distance from the center of an ellipse to its edge,

[TABLE]

where $L=\sqrt{1+2gN/\ell_{z}}$ and $a$ is a parameter which determines the amplitude of motion. Similarly, $\ell_{z}(t)$ is the distance from a focus of an ellipse to its edge. The graph is a cycloid, defined parametrically by

[TABLE]

Here $\epsilon^{2}=1+{E_{z}}/({2g^{2}N^{2}\left\langle\ell_{r}^{-2}\right\rangle^{2}}),A=(2|E_{z}|(1-\epsilon^{2}))^{-1/2},B=(4E_{z}^{2}(1-\epsilon^{2}))^{-1/2}$ and $E_{z}$ is a parameter determining the amplitude of the motion.

IV.1.2 Solving equations of motion

To find the analytic results in Eqs. (21) and (26) we note that Eqs. (20) and (19) are of the same form as the central force equation:

[TABLE]

with $V\propto r^{2}$ and $V(r)\propto-{1}/{r}$. These correspond to the harmonic and gravitational potential respectively. We will find $\vec{r}(t)$, then our desired solutions are $\ell(t)=r(t)=|\vec{r}(t)|$.

For the harmonic oscillator, Eq. (19), we can characterize the elliptical motion by the semimajor axis of the orbit $a$, and the angular momentum $L=\sqrt{1+2gN/\ell_{z}}$. The energy is $E_{r}={a^{2}}/{2}+{L^{2}}/({2a^{2}})$. The particle’s trajectory can be written as

[TABLE]

Thus we see $\ell_{r}=|\vec{r}|$ is given by Eq. (21).

For Eq. (20), we can characterize the motion of the particle by the energy $E_{z}$. The equations of motion for $r(t)=\ell_{z}(t)$ can be integrated with a parametric substitution, as found in Landau’s Mechanics landau1982 . The solution is a cycloid:

[TABLE]

where $\epsilon^{2}=1+{E_{z}}/({2g^{2}N^{2}\left\langle\ell_{r}^{-2}\right\rangle^{2}})$ is the eccentricity of the orbit. Bound gravitational orbits, which correspond to self-trapped solitons, have $E_{z}<0$. When $E_{z}$ approaches zero from below, $\ell_{z,max}$ diverges as $1/|E_{z}|$ while $\ell_{z,min}$ approaches a constant value. Physically, $E_{z}>0$ means that a soliton is no longer self-trapped and can expand to infinite size. However, any nonzero trapping in the axial direction ($\omega_{z}\neq 0$) would restrain this unlimited expansion.

IV.1.3 Invariants

One insightful result from the decoupling approximation is the existence of certain invariants, which are found by considering the angular momentum $L=mr^{2}\dot{\theta}$ of the classical system. The integral of $\dot{\theta}$ over one period of radial oscillation is $\pi$ or $2\pi$ in the harmonic and graviational problems. For the radial motion this yields

[TABLE]

where we have used that the period of oscillation is $T_{r}=\pi$. For the axial motion we instead find

[TABLE]

where $T_{z}=2\pi/[4g^{2}N^{2}(1-\epsilon^{2})^{3/2}\langle\ell_{r}^{-2}\rangle^{2}]$ is the period of the axial motion.

We have three applications of these invariants. First, we used Eq. (28) in our adiabatic decoupling. Second, these invariants are useful in analyzing the evolution of the phase $\chi$. Though $\chi$ is a physically irrelevant global phase for the case of the single soliton, in the context of a soliton train, phase differences between adjacent solitons affect the strength and sign of the soliton-soliton interaction. The phase evolution in Eq. (18) is dominated by terms proportional to $\langle\ell_{r}^{-2}\rangle$, which are independent of the amplitude of radial oscillation. Consequently the phase relationship between neighboring solitons in a train will be robust against the oscillation of individual solitons. There is some experimental evidence for this phenomenon hulet2017 .

Finally, these invariants can be used to analyze the atom loss in Eq. (17). The time averaged loss is proportional to $\langle\ell_{z}^{-2}\rangle\approx 4g^{2}N^{2}(1-\epsilon^{2})^{3/2}$.

IV.1.4 Atom Loss

As explained in the introduction, the experimentally observed atom loss rate is orders-of-magnitude larger than one would expect for static solitons. We investigated whether or not the oscillations of the solitons were responsible for this discrepency. Using Eqs. (28) and (29), the average of $\dot{N}$ over these oscillations is

[TABLE]

The atom loss rate depends on the amplitudes of radial and axial oscillation through $E_{r}$ and $\epsilon$, respectively. In the absence of oscillations, $\epsilon=0$ and $E_{r}=1+{\cal O}(gN/\ell_{z})$. Axial oscillations increase $\epsilon$, and hence decrease $\dot{N}$. Radial oscillations increase $E_{r}$ and hence increase $\dot{N}$. This dependence, however, is weak. If $(\ell_{z})_{\rm max}/(\ell_{z})_{\rm min}=2$ one has $\epsilon^{2}=1/9$, while if $(\ell_{r})_{\rm max}/(\ell_{r})_{\rm min}=2$ one has $E_{r}=1.25$. We conclude that the oscillations cannot explain the orders-of-magnitude enhancement of the loss in hulet2017 .

We can integrate Eq. (30) over time-scales large compared to the oscillation period:

[TABLE]

Fig. 4c shows the atom loss in a soliton undergoing a large amplitude axial oscillation with $\epsilon=0.87$. Thus, the dynamic soliton loses atoms at a lower rate: the green curve in the plot, given by Eq. (31), decreases more slowly than the red curve, which represents the atom number of a non-oscillating soliton.

IV.1.5 Beyond Decoupling

To study corrections to the decoupling approximation, we numerically integrate Eqs. (13) through (18). In Fig. 4 the blue lines show the numerical solution within the Gaussian approximation, while the green lines show the closed-form predictions from Eqs. (21), (26), and (31).

Though these closed-form solutions are fairly accurate for large ranges of the parameters, they do not explain certain features of the dynamics. First, the soliton loses atoms at the greatest rate when the cloud is smallest during each oscillation. This short timescale behavior is not captured in Eq. (31), but when averaged over one period, the closed-form prediction is quite accurate.

Second, the analytic solutions do not capture the slow equilibration between the radial and axial motions. The interaction term can slowly shift energy between the radial and axial modes, until their amplitudes have reached equilibrium. With the current experimental conditions, though, the soliton degrades long before its motion approaches equilibrium. Thus, this effect is not currently experimentally observable.

It is possible to find perturbative analytic solutions which capture these two behaviors. However, for this level of detail, it is more illuminating to investigate the features of the full numerical results.

IV.2 Beyond Gaussian Approximation

To more accurately model the dynamics we integrate the coupled ordinary differential equations produced by Eqs. (10) and (11) for more complicated ansatzes with larger values of $n$. Fig. 4 shows an example simulation comparing an ansatz with $n=2$ (with 6 complex variational parameters controlling the functional form) to the Gaussian ansatz (with 1 complex variational parameter controlling the functional form) and to the case where the soliton is not oscillating.

The radial dynamics in panel (a) of Fig. 4 show quantitative agreement between the analytic solution and the full simulations. This agreement is also underscored in Fig. 5, where we compare our variational calculations to a numerical method based on discretizing Eq. (1). This discretization approach follows javanainen2006 and uses a split-step integrator, taking an adimensional spatial grid spacing of 0.05 and a timestep of 0.001. We ran the simulation over only one radial oscillation because of the computational intensity of the discretization method. Fig. 5 validates the variational results, since they are nearly indistinguishable from the results of spatial discretization.

For the axial dynamics in panel (b) of Fig. 4, quantitative predictive power is lost after a single axial oscillation. The decoupling approximations break down when the soliton is sufficiently compressed. To illustrate this breakdown, in Fig. 6a, we show a longer simulation using the n=2 ansatz (16 real variational parameters) with the same initial conditions as in Fig. 4. We use a density-plot to show the ”axial density”, $\rho(z)=\int\!d^{3}rn(r,z)$ as a function of time. For an excitation of this amplitude, the dynamics are quite rich. Energy moves between the various degrees of freedom, and it is likely that the behavior is chaotic: the variations in the amplitude cause large variations in the period of oscillation, which in turn affect the amplitude. Fig. 6b shows the results of using the Gaussian ansatz, which does not contain enough degrees of freedom to capture the richness of the dynamics. The Gaussian approximation is, however, reliable over shorter timescales or for small amplitudes.

V Summary

To summarize, we studied a systematic set of variational ansatzes for the static and dynamic properties of bright solitons in an attractive Bose gas in a cylindrically symmetric trap. We presented a precise variational calculation of $k_{c}=N_{c}|a_{s}|\sqrt{m\tilde{\omega}_{r}/\hbar}$, where $N_{c}$ is the critical atom number, above which a single soliton will collapse. We found that for a cylindrical harmonic potential, $k_{c}=0.677986(2)$. We studied the dynamics, and found closed-form solutions to the equations from the simplest ansatz by a formal mapping to two familiar classical systems. This simple solution captures all of the significant qualitative features and some of the quantitative features of the far-from-equilibrium dynamics of a bright soliton.

The oscillations of the soliton affect the atom loss, causing it to occur episodically. When time-averaged, however, we find that the loss rate only depends weakly on the oscillation amplitude. We conclude that the surprisingly large loss rate seen in the experiments [8] cannot be explained by this dynamical effect.

Although we solely considered the behavior of a single soliton, the generalization to a soliton train is straightforward. Similar techniques could also be used to explore the soliton generation process, where a monolithic atomic cloud breaks up into an array of solitons.

Acknowledgements

We would like to acknowledge Randy Hulet for discussions about the experiments and the physics of cold atom solitons, Junkai Dong for discussions of both techniques and phenomena, and Matthiew Reichl for providing some of the computer code which formed the basis for our spatial discretization approaches. This material is based upon work supported by the National Science Foundation under Grant No. PHY-1806357 and the ARO-MURI Non-equilibrium Many-body Dynamics Grant No. W9111NF-14-1-0003.