Excitation modes of bright matter-wave solitons

TL;DR

This paper experimentally investigates the excitation modes of bright matter-wave solitons in a quasi-one-dimensional Bose-Einstein condensate, analyzing their breathing frequencies, higher-order modes, and wave packet splitting.

Contribution

It provides new experimental insights into the excitation dynamics of matter-wave solitons, including the effects of quench parameters on their modes and comparison with theoretical models.

Findings

Breathing modes are excited and their frequencies depend on atom number and confinement.

Signatures of higher-order solitons and wave packet splitting are observed.

Experimental results agree with analytical and numerical predictions.

Abstract

We experimentally study the excitation modes of bright matter-wave solitons in a quasi-one-dimensional geometry. The solitons are created by quenching the interactions of a Bose-Einstein condensate of cesium atoms from repulsive to attractive in combination with a rapid reduction of the longitudinal confinement. A deliberate mismatch of quench parameters allows for the excitation of breathing modes of the emerging soliton and for the determination of its breathing frequency as a function of atom number and confinement. In addition, we observe signatures of higher-order solitons and the splitting of the wave packet after the quench. Our experimental results are compared to analytical predictions and to numerical simulations of the one-dimensional Gross-Pitaevskii equation.

Figures

Equations

$$f(z)=\frac{1}{\sqrt{2l_z}}\mathrm{sech}\left(\frac{z}{l_z}\right),$$

$$T^{(2)}=\frac{8\pi }{\hslash }m{\left(l_z^{(1)}\right)}^2.$$

$$i\hslash \frac{\partial }{\partial t}\psi (\mathbf{r},t)=\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+V(\mathbf{r})+gN{\left|\psi (\mathbf{r},t)\right|}^2\right]\psi (\mathbf{r},t),$$

$$\psi (\mathbf{r},t)=f(z,t)\frac{1}{\sqrt{\pi }{a}_r\sigma (z,t)}\exp \left[-\frac{(x^2+y^2)}{2a_r^2\sigma (z,t{)}^2}\right],$$

$$\begin{array}{rl}i\hslash \frac{\partial }{\partial t}f(z,t)& =[-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial z^2}+V(z)\\ & +\frac{gN}{2\pi a_r^2\sigma (z,t{)}^2}{\left|f(z,t)\right|}^2\\ & +\frac{\hslash {\omega }_r}{2}\left(\sigma (z,t{)}^2+\frac{1}{\sigma (z,t{)}^2}\right)]f(z,t),\end{array}$$

$$\sigma (z,t{)}^2=\sqrt{1+2aN|f(z,t){|}^2}.$$

$$\begin{array}{rl}i\hslash \frac{\partial }{\partial t}f(z,t)& =[-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial z^2}+V(z)\\ & +\frac{gN}{2\pi a_r^2}{\left|f(z,t)\right|}^2+\hslash {\omega }_r]f(z,t).\end{array}$$

$$f(z,0)=\frac{1}{\sqrt{2l_z}}\mathrm{sech}\left(\frac{z}{l_z}\right)$$

$$z_{rms}(t)={\left(\frac{1}{N}\int n(z,t)(z-\bar{z}{)}^{2}dz\right)}^{1/2}.$$

$$\begin{gathered} E[\psi ]=\int d^3\mathbf{r}[\frac{{\hslash }^2}{2m}{|\nabla \psi (\mathbf{r})|}^2+ \\ \frac{1}{2}m({\omega }_r^2{r}^2+{\omega }_z^2{z}^2){|\psi (\mathbf{r})|}^2+\frac{gN}{2}|\psi (\mathbf{r}){|}^4] \end{gathered}$$

$$\psi (r,z)=\frac{1}{\sqrt{2l_z}}\mathrm{sech}\left(\frac{z}{l_z}\right)\cdot \frac{1}{\sqrt{\pi l_r}}\exp \left(-\frac{r^2}{2l_r^2}\right),$$

$${\epsilon }_{GP}=\frac{1}{2{\gamma }_r^2}+\frac{{\gamma }_r^2}{2}+\frac{1}{6{\gamma }_z^2}+\frac{{\pi }^2}{24}{\lambda }^2{\gamma }_z^2+\frac{\alpha }{3{\gamma }_r^2{\gamma }_z},$$

$$\frac{{\pi }^2}{4}{\lambda }^2{\gamma }_z^4+\sqrt{\zeta }{\gamma }_z-1=0,$$

$${\gamma }_z^{\ast }=\frac{F}{\sqrt{\zeta }},$$

$$F=-\sqrt{\frac{G}{2}}+\frac{1}{2}\sqrt{-2G+\frac{4\sqrt{2}}{{\pi }^2\sqrt{G}}{\left(\frac{\zeta }{\lambda }\right)}^2},$$

$$G=\frac{\Delta }{{\pi }^{\frac{4}{3}}}{\left(\frac{\zeta }{\lambda }\right)}^{\frac{4}{3}}-\frac{4}{3{\pi }^{\frac{2}{3}}}\frac{1}{\Delta }{\left(\frac{\zeta }{\lambda }\right)}^{\frac{2}{3}},$$

$$\Delta =\sqrt[3]{1+\sqrt{1+\frac{64{\pi }^2}{27}{\left(\frac{\lambda }{\zeta }\right)}^2}}.$$

$$\left(\frac{{\pi }^2}{12}\right){\ddot{\gamma }}_z=\frac{{\gamma }_z^{-3}}{3}-\frac{{\pi }^2}{12}{\lambda }^2{\gamma }_z+\frac{\alpha }{3}{\gamma }_z^{-2},$$

$${\omega }_z={\omega }_r\sqrt{\frac{12}{{\pi }^2}\left({{\gamma }_z^{\ast }}^{-4}+\frac{{\pi }^2}{12}{\lambda }^2+\frac{2\alpha }{3}{{\gamma }_z^{\ast }}^{-3}\right)}.$$

Full text

The dispersionless propagation of solitary waves is one of the most striking features of nonlinear dynamics, with multiple applications in hydrodynamics, nonlinear optics and broadband long-distance communications T. Dauxois (2006). In fiber optics, one-dimensional (1D) “bright” solitons, i.e. solitons presenting a local electric field maximum with one-dimensional propagation, have been observed Mollenauer et al. (1980). They exhibit a dispersionless flow and excitation modes such as breathing or higher-order modes Satsuma and Yajima (1974); Mollenauer et al. (1980); Stolen et al. (1983). Matter waves can also display solitary dispersion properties. Typically, bright matter-wave solitons are created in quasi-1D systems by quenching the particle interaction in a Bose-Einstein condensate (BEC) from repulsive to attractive Khaykovich et al. (2002). Recent experiments demonstrated the collapse Donley et al. (2001), collisions Nguyen et al. (2014), reflection from a barrier Marchant et al. (2013), and the formation of trains Strecker et al. (2002); Nguyen et al. (2017); Mežnaršič et al. (2019) of bright solitons.

In this letter, we experimentally study the excitation modes of a single bright matter-wave soliton. In previous studies, other dynamical properties have been observed, such as center-of-mass oscillations of solitons in an external trap Nguyen et al. (2014), excitations following the collapse of attractive BECs Donley et al. (2001); Cornish et al. (2006), and quadrupole oscillations of attractive BECs in three dimensions (3D) Everitt et al. (2015). Here, we probe the fundamental breathing mode of a single soliton by measuring its oscillation frequency and the time evolution of its density profile. In addition, we observe signatures of higher-order matter-wave solitons, which can be interpreted as stable excitations with periodic oscillations of the density profile and phase, or as a bound state of overlapping modes Satsuma and Yajima (1974); Carr and Castin (2002).

The shape-preserving evolution of a matter-wave soliton is due to a balancing of dispersive and attractive terms in the underlying 3D Gross-Pitaevskii equation (GPE) Dalfovo et al. (1999). For quasi-1D systems with tight radial confinement, we can approximate the matter wave in the 3D-GPE by the product of a Gaussian wave function for the radial direction and a function $f(z)$ for the longitudinal direction (see Sup ). Depending on the ansatz for the Gaussian with either constant or varying radial sizes, $f(z)$ satisfies either the 1D-GPE or the non-polynomial Schrödinger equation (NPSE) Salasnich et al. (2002). We reference to the analytical solutions of the 1D-GPE in the manuscript, but use both equations in our numerical simulations Sup .

For the 1D-GPE, an ansatz for the normalized longitudinal wave function $f(z)$ is of the form

[TABLE]

with a single parameter, $l_{z}$, that determines both the longitudinal size and the amplitude of the soliton. Solitons form with a value of $l_{z}$ that minimizes the total energy and that provides a compromise between the kinetic and the interaction energies. This is illustrated in Fig. 1b, which shows the energy of the wave packet for varying sizes $l_{z}$ Parker et al. (2007). The kinetic energy provides a potential barrier for small $l_{z}$ that prevents the collapse of the soliton, while its spreading is inhibited by the interaction energy, which increases for large $l_{z}$.

Even without an external longitudinal potential, the soliton is stable against small perturbations of $l_{z}$. In a way, a bright matter-wave soliton creates its own trapping potential, which defines its size and excitation modes. Variational methods provide accurate predictions of its size at the energy minimum which can be calculated analytically Carr and Castin (2002); Longenecker and Mueller (2019) or numerically Parker et al. (2007). For the fundamental solution (order $n=1$) of the 1D-GPE with an atom number $N$, s-wave scattering length $a$ and radial trapping frequency $\omega_{r}$, the size $l_{z}$ corresponds to the healing length at the peak density of the soliton, i.e. $l_{z}^{(n=1)}=a_{r}^{2}/(N|a|)$ Carr and Castin (2002); Parker et al. (2007). Here, $a_{r}=\sqrt{\hbar/(m\omega_{r})}$ is the radial harmonic oscillator length. Small deviations of $l_{z}$ close to the energy minimum lead to oscillations of the soliton size. We use those oscillations resulting from an initial mismatch of $l_{z}$ to experimentally measure the self-trapping frequency of the soliton potential.

Our experimental starting point is a Bose-Einstein condensate of $500-2000$ cesium (Cs) atoms in the state $\ket{F=3,m_{F}=3}$ at scattering length of $a=+7\,a_{0}$, where $a_{0}$ is Bohr’s radius. The BEC is levitated by a magnetic field gradient, and it is confined by an optical dipole trap formed by the horizontal and vertical laser beams $L_{H}$ and $L_{V}$ (Fig. 1a). An additional magnetic offset field allows us to tune the scattering length by means of a broad magnetic Feshbach resonance Berninger et al. (2013). Details about our experimental setup, the levitation scheme and the removal of atoms can be found in references Carli et al. (2019); Sup .

Our matter-wave solitons are confined to a quasi-1D geometry with almost free propagation along the horizontal direction and strong radial confinement of $\omega_{r}=2\pi\times 95$ Hz provided by laser beam $L_{H}$. They are generated with a quench of the scattering length towards attractive interaction ($a_{i}\rightarrow a_{f}$), and by a reduction of the longitudinal trap frequency ($\omega_{z,i}\rightarrow\omega_{z,f}$). When changing $a$ and $\omega_{z}$ independently, the quenches excite inward and outward motions, respectively. Usually, it is desirable to minimize the excitations of the soliton by matching the initial Thomas-Fermi density profile of the BEC closely to the density profile of the soliton (inset Fig. 1a). However, we deliberately mismatch the quench parameters to create breathing oscillations of the soliton in order to study its self-trapping potential. Quenches with different parameters are labelled by the symbols Q1-Q7 (see Sup ). Following an evolution time $t$ in quasi-1D and after a short period of 16 ms of expansion in free space, we take absorption images to determine the density profile of the atoms (Fig. 1c). The cloud size, $l_{z}(t)$, is determined by fitting the function $A(\mathrm{sech}(z/B))^{2}$ to the integrated 1D-density profiles with fit parameters $A$ and $B$ Sup .

The response of the atomic cloud to the different steps of the quench protocol is presented in Fig. 1d. We first quench only the longitudinal confinement by 25% to $\omega_{z,f}=2\pi\times 4.3(2)$ Hz (quench Q1 in Sup ) while keeping the repulsive interaction strength constant (Fig. 1d, diamonds). The BEC starts an outwards motion with an oscillation frequency of $2\pi\times 7.5(1)\,\mathrm{Hz}\approx\sqrt{3}\omega_{z,f}$ as expected for a BEC in the Thomas-Fermi regime Menotti and Stringari (2002); Haller et al. (2009). In a second step, we additionally quench the interaction strength $a_{f}$ to $-5.4\,a_{0}$ and increase $\omega_{z,i}$ to match the initial size of the BEC to the expected size of the soliton, (Q2, Fig. 1d, squares). As a result, we observe almost dispersionless solitons with a linear increase of $l_{z}$ of $0.7(3)\,\mu$m/s (Fig. 1d, green line). Finally, we deliberately mismatch the initial size of the BEC by reducing $\omega_{z,i}$ (Q3), and generate small amplitude oscillations of the soliton with a frequency $\omega_{sol}$ of $2\pi\times 12.8(4)\,$Hz (Fig. 1d, circles). This breathing frequency of the soliton is significantly larger than any breathing frequency of a BEC or of non-interacting atoms, $2\omega_{z,f}=2\pi\times 8.6(3)\,$Hz. We observe no discernible oscillation in the radial direction after the quenches.

In a second experiment, we demonstrate that the breathing frequency, $\omega_{sol}$, depends on the interaction term $Na$ in the 1D-GPE, a property typical of the nonlinear character of the soliton. We choose to change $N$, since the initial removal process is independent of the interaction quench, and we can study $\omega_{sol}$ without changing the quench protocol (Q4, Fig. 2a, circles). The measured values of $\omega_{sol}$ decrease for lower $N$, and they approach the breathing frequency $2\omega_{z,f}$ for non-interacting atoms in a harmonic trap (Fig. 2a, dashed line).

We compare our experimental data points to two theoretical models. In a numerical simulation of the 1D-GPE, we use the ansatz in Eq. 1 to set the starting conditions, and we determine the breathing frequency from a spectral analysis of the time evolution of the wave function Sup (Fig. 2a, triangles). In addition, we use an analytical approximation for the breathing frequency (red line) calculated with a Lagrangian variational analysis at the energy minimum of the 3D-GPE Carr and Castin (2002); Sup ). We find that both models agree well with the trend of the measurements of $\omega_{sol}$, although our experimental data points are systematically lower for large $N$ than our theoretical predictions. We speculate that this is due to non-harmonic contributions to the energy of the soliton on the breathing oscillations for finite oscillation amplitudes (Fig. 1b).

To determine the influence of the trapping potential, we measure the variation of $\omega_{sol}$ as we reduce the longitudinal trapping frequency $\omega_{z,f}$ (Q5). Two regimes of $\omega_{sol}$ can be identified in Fig. 2b for varying the values of $\omega_{z,f}$. For large values of $\omega_{z,f}$, the trap dominates the breathing of the soliton and $\omega_{sol}$ approaches twice the trap frequency $2\omega_{z,f}$. For small values of $\omega_{z,f}$, interactions dominate the breathing of the soliton and $\omega_{sol}$ reaches a constant value. This offset of the breathing frequency is a result of the “self-trapping” potential of a free soliton.

Again, we compare the experimental results with our theoretical model (Fig. 2b, red line) and the numerical simulations of the 1D-GPE. The blue band in Fig. 2b indicates the simulated frequencies for $N=1300$ to $N=1500$. The simulation predicts a lower breathing frequency for the free soliton than the analytical approximation, but all curves are within the uncertainly range of the experimental data.

External trapping potentials can in principle alter the soliton dynamics Serkin et al. (2015); Serkin and Belyaeva (2019); Nguyen et al. (2014), causing, e.g., modulations of the soliton’s tails due to residual non-autonomous terms of the 1D-GPE in a harmonic potential Tenorio et al. (2005). For the following experiments, however, we employ trap frequencies that are significantly smaller than the observed oscillation frequencies of the soliton ($2\omega_{z}<\omega_{sol}$) to decouple the influence of the trapping potential. In summary, for small-amplitude oscillations we find good agreement of $\omega_{sol}$ between our experimental results and analytical and numerical predictions based on the 1D-GPE (and NPSE Sup ).

Breathing oscillations of $l_{z}$ close to the equilibrium size are not the only possible excitation modes of solitons. The existence of higher-order solitons has been predicted in the nonlinear Schrödinger equation (NSE) Satsuma and Yajima (1974), and has been observed for optical solitons in silica-glass fibers Mollenauer et al. (1980); Stolen et al. (1983). A soliton of order $n$ can be interpreted as a bound state of $n$ strongly overlapping solitons Carr and Castin (2002). By exploiting the equivalence of the NSE and 1D-GPE, similar effects were later proposed for bright matter-wave solitons Carr and Castin (2002); Golde et al. (2018), where it was suggested that nth-order solitons can be generated by a rapid increase of the attractive interaction strength by a factor $n^{2}$. Similarly, our simulations of the 1D-GPE show that higher-order solitons can be created for an increased initial size of the wave packet. An $n^{th}$-order soliton forms for a sech-shaped wave function with an initial size $l_{z}^{(n)}$ that is the $n^{2}$ multiple of the healing length $l_{z}^{(1)}$, i.e. $l_{z}^{(n)}=n^{2}\;l_{z}^{(1)}$ Sup .

Within the 1D-GPE theory, both creation methods result in the periodic development of multi-peaked structures for higher-order solitons Satsuma and Yajima (1974); Gamayun and Semenyakin (2016), e.g. they create a sharp central peak with side wings for a second-order soliton (Fig. 3a) and a double peak for a third-order soliton Sup . Sizes and interaction quenches that do not fulfil the previous conditions, lead to a “shedding” of the atomic density in the $z-$direction. The wave packet oscillates and loses particles until its size and shape match the next (lower $n$) higher-order soliton Satsuma and Yajima (1974). For a second-order soliton, the predicted oscillation period $T^{(2)}$ is Carr and Castin (2002)

[TABLE]

Recently, excitation modes of higher-order have also been used as a testbed for various theoretical models beyond GP-theory. The fragmentation of solitons with an increased initial width was predicted within the multiconfigurational time-dependent Hartree method for bosons Streltsov et al. (2008) and critically discussed Cosme et al. (2016), and the influence of quantum effects on the dissociation process was investigated Weiss and Carr (2016); Yurovsky et al. (2017); Ng et al. (2019).

Here, we apply two different quench protocols to study the evolution of strongly excited solitons. Depending on the initial size and the quench parameters, we observe shedding and fragmentation of the wave packet, and we measure oscillation frequencies that indicate the creation of higher-order solitons. To demonstrate the effect of a strong quench of an elongated BEC, we increase $a_{i}$ and reduce $\omega_{z,i}$ before ramping $a$ and $\omega_{z}$ to $-5.3\,a_{0}$ and $2\pi\times 0.0(6)\,$Hz in $13\,$ms (Q6). Our quench induces an initial spreading of the wave packet, followed by a strong shedding of atoms, and finally, in the formation of a soliton that contains approximately 1/3 of the initial atom number (Fig. 3b). We determine the soliton width and find a slow oscillation of $l_{z}(t)$ with a frequency of $2\pi\times 2.4(2)$ Hz (Fig. 3c). This frequency is significantly smaller than the expected breathing frequency of first-order solitons, $2\pi\times 6.0$ Hz, and it matches well to the expected frequency of $2\pi\times 2.3$ Hz for second-order solitons in Eq. 2.

Observing shedding and oscillations agrees well with the predictions for higher-order solitons within the 1D-GPE Satsuma and Yajima (1974), however, we find a strong dependence on details of the quench protocol and on the dynamical evolution during the quench. For a closer match to theoretical works Carr and Castin (2002), we implement a double-quench protocol, with a first quench to generate a soliton with weak attractive interaction, $a_{f}=-0.8\,a_{0}$, $\omega_{f}=2\pi\times 1.4(2)\,$Hz, and, after a settling time of $25\,$ms, a second quench of only the interaction strength, $a_{f}=-4.6\,a_{0}$ (Q7). Starting with approximately 2200 atoms, we observe no shedding but a small loss of 300 atoms during the first 60 ms. The density distributions (Fig. 4a) resemble the expected profiles of a second-order soliton (Fig. 3a), and the width of the wave packet oscillates with a frequency of $2\pi\times 5.6(6)$ Hz (Fig. 4b), which matches the expected frequency of $2\pi\times 5.2$ Hz for second-order solitons ($2\pi\times 13.2$ Hz for the first-order).

For both measurements (Fig. 3c, 4b), a small percentage of absorption images show the splitting of the soliton into two fragments (inset Fig. 4b). Due to the destructive nature of our absorption images it is difficult to conclude on the evolution and on the cause of the splitting process. A double-peak structure in the density profile can indicate the generation of a third-order soliton, fragmentation due to quantum effects, or simply an insufficient technical control of our quench parameters. For our setup, the control of horizontal magnetic field gradients to avoid longitudinal accelerations is especially challenging Carli et al. (2019). The percentage of images that show a splitting of the wave packet increases for longer evolution times, and we indicate their fraction in Fig. 4b with a histogram.

In conclusion, we experimentally studied the creation and the excitation of breathing modes of bright matter-waves solitons in a quasi-one-dimensional geometry after a quench of interaction and longitudinal confinement. We measured the “self-trapping” frequency $\omega_{sol}$ for first-order solitions and its dependence on $N$ and $\omega_{z}$. For stronger excitations and for a double-quench protocol, we observed signatures of second-order solitons and the shedding and splitting of the wave function. Further measurements of the splitting process and the damping of the oscillations due to shedding are necessary to distinguish technical fluctuations from higher-order solitons and fragmentation due to quantum effects Weiss and Carr (2016); Yurovsky et al. (2017); Ng et al. (2019).

We thank L. D. Carr for helpful initial discussions during his visit. We acknowledge support by the EU through “QuProCS” (GA 641277) and the ETN “ColOpt” (GA 721465), and by the EPSRC Programme Grant “DesOEQ” (EP/P009565/1). AdC acknowledges financial support by EPSRC and SFC via the IMPP.

Supplemental Material

I Experimental methods

I.1 Controlling the atom number in the BEC

The solitons are confined to a quasi-1D geometry with almost free propagation along the horizontal direction and strong radial confinement of $\omega_{r}=2\pi\times 95$ Hz provided by laser beam $L_{H}$. In quasi-1D geometry, bright matter-wave solitons collapse for large densities and interactions Donley et al. (2001), which for our typical experimental scattering length of approximately $-5\,a_{0}$ corresponds to a critical atom number of 2500 Carr and Castin (2002). As a result, we need to strongly reduce the atom number to avoid collapse, modulation instabilities Strecker et al. (2002) and three-body loss Everitt et al. (2017) for a deterministic and reproducible creation of the soliton. We remove atoms with a small additional magnetic field gradient, which pushes the atoms over the edge of the optical dipole trap. Our precise control of magnetic field strengths allows us to reduce the atom number down to $200$ atoms, with a reproducibility of $\pm 100$ for 600 atoms and $\pm 350$ for 4500 atoms, measured as the standard deviation of the atom number in 50 consecutive runs. A removal period of 4 s and smooth ramps of the magnetic field strength are necessary to minimize excitations of the BEC. Following the removal procedure we measure residual fluctuations of the width of the BEC below 3.5%.

I.2 Quench parameters

Several different quench protocols are employed for the measurements. The quenches are labeled by the symbols Q1-Q7 in the main article:

• Q1

We quench only the trap frequency from $\omega_{z,i}=2\pi\times 5.8(2)$ Hz to $\omega_{z,f}=2\pi\times 4.3(2)$ Hz with a linear ramp of the laser power of beam $L_{V}$ over $4\,$ms. Atom number $N\approx 1800$, constant interaction strength $a_{i}=+7\,a_{0}$, $\omega_{r}=2\pi\times$95 Hz.

• Q2

In addition to the quench Q1 of the trap frequency, we also quench the interaction strength from $a_{i}=+7\,a_{0}$ to $a_{f}=-5.4\,a_{0}$ in 4 ms. We minimize oscillations of the width of the soliton by reducing the initial size of the BEC with $\omega_{z,i}=2\pi\times 11.2(2)\,$Hz. Atom number $N\approx 1800$, $\omega_{z,f}=2\pi\times 4.3(2)$ Hz, $\omega_{r}=2\pi\times$95 Hz.

• Q3

We mismatch the initial size of the BEC before the quench with $\omega_{z,i}=2\pi\times 12.8(4)$ Hz to generate small amplitude oscillations of the width of the soliton. Atom number $N\approx 1700$, $\omega_{z,f}=2\pi\times 4.3(2)$ Hz, $\omega_{r}=2\pi\times$95 Hz.

• Q4

Quench to determine the atom-number dependence of $\omega_{sol}$. We vary the atom number $N$ from 500 to 1700 for the measurement. $\omega_{z,f}=2\pi\times 4.3(2)$ Hz, $a_{i}=+7\,a_{0}$, $a_{f}=-5.4\,a_{0}$, ramp duration 4 ms, $\omega_{r}=2\pi\times$95 Hz.

• Q5

Quench to determine the dependence of $\omega_{sol}$ on the trap frequency $\omega_{z,f}$. We vary $\omega_{z,f}$ from approximately 1 Hz to 9 Hz. Smaller values of $\omega_{z,f}$ result in larger equilibrium sizes of the soliton, and we need to reduce the initial trap frequencies $\omega_{z,i}$ to keep the oscillation amplitudes comparable during the measurements. The typical difference between $\omega_{z,i}$ and $\omega_{z,f}$ is approximately $3\,$Hz. $N\approx 1500$, $a_{f}=-5.4\,a_{0}$, $\omega_{r}=2\pi\times$95 Hz.

• Q6

Strong quench starting from an elongated BEC to excite higher-order oscillations and shedding. The ratio between the calculated initial Thomas-Fermi radius of the BEC and the expected width $l_{z}$ of the soliton is 24. $a_{i}=56\,a_{0}$, $a_{f}=-5.3\,a_{0}$, $\omega_{z,i}=2\pi\times 4.9(2)$ Hz, $\omega_{z,f}=2\pi\times 0.0(6)$ Hz, initial atom number $N\approx 3000$ drops to 1100 after shedding of atoms, quench duration $13\,$ms, $\omega_{r}=2\pi\times$86 Hz.

• Q7

Double quench to create a stable soliton in step 1 and quench the scattering length by approximately a factor of 4 in step 2. Step 1: $\omega_{z,i}=2\pi\times 4.9(2)$ Hz, $\omega_{z,f}=2\pi\times 1.4(2)$ Hz, $a_{i}=29\,a_{0}$, $a_{f}=-0.8\,a_{0}$, $\omega_{r}=2\pi\times$86 Hz, quench duration 15 ms, $N\approx 2200$. Settling delay of 25 ms between quenches. Step 2: reduce interaction strength in 2 ms to $a_{f}=-4.6\,a_{0}$, no change of other parameters.

I.3 Fit of density profiles

We employ absorption imaging to measure the 2D-density profile of the soliton, and we integrate over one radial axis to determine the 1D-density profile (red line in Fig. 5). The width $l_{z}$ of the soliton in Eq. 1 of the main article, is determined by fitting the function $A(\mathrm{sech}(z/B))^{2}$, with fit-parameters $A$ and $B$, to the integrated 1D-density profiles (dotted blue line in Fig. 5).

II Theoretical methods

II.1 The Model

The time evolution of the collective wave function of $N$ atoms in an external potential with the 3D Gross-Pitaevski equation (GPE) for a time and space-dependent collective atomic wave-function, $\psi(\mathbf{r},t)$, is given by,

[TABLE]

where $g=4\pi\hbar^{2}a/m$, $m$ is the atomic mass, and $a$ is the two-body s-wave scattering length. This semi-classical field equation can be seen as a mean-field computation, and describes the dynamics of many weakly interacting particles at low temperatures when the condition $n|a|^{3}\ll 1$ is satisfied Cowell et al. (2002), where $n$ is the particle density. Our external potential $V(\mathbf{r})$ is given by a 3D (anisotropic) harmonic trap.

For tight radial trapping potentials, $\omega_{r}\gg\omega_{z}$, we can approximate the 3D wave function with a Gaussian solution in the radial directions and an arbitrary component, $f(z,t)$, in the longitudinal direction,

[TABLE]

where $a_{r}$ is the harmonic oscillator length in the radial direction and $\sigma(z,t)$ is a free parameter that dictates the width of the radial wavefunction. Substituting this ansatz into the 3D-GPE, and integrating over the radial directions, we arrive at the so called non-polynomial Schödinger equation (NPSE) Salasnich et al. (2002)

[TABLE]

where $\omega_{r}={\hbar/ma_{r}^{2}}$. The condition for $\sigma(z,t)$ that minimizes the action functional integrated along the trajectories in phase space is Salasnich et al. (2002),

[TABLE]

For $\sigma(z,t)=1$, we obtain the ground state of a harmonic oscillator in the radial directions, and we recover the usual 1D-GPE

[TABLE]

We have numerically integrated Eqs 5 and 7 using the split-step Fourier transform method Suarez (2015), where we exploit the fact that the kinetic and potential terms in the Hamiltonian are diagonal in momentum and real space, respectively.

III Soliton Breathing Frequency

In this section we explain how the numerical calculations of the soliton breathing frequencies shown in Fig. 2 of the main text were carried out. We begin with the order 1 soliton solution,

[TABLE]

where $l_{z}=a_{r}^{2}/(N|a_{i}|)$, and we have used $a_{i}=7\,a_{0}$. We then evolve this initial state either with the 1D-GPE or NPSE to a simulation time of $4000$ ms and evaluate the frequency spectrum of the oscillation of the soliton’s centre ($z=0$). In Fig. 6 we present the frequency spectrum for the GPE and a longitudinal frequency of $\omega_{z}=2\pi\times 5$ Hz and atom number $N=1300$, which is characteristic of the behaviour for all other $\omega_{z}$ data points. We observe several prominent frequency modes in the signal, but we select the lowest frequency peak to compare to the experimental measurements, because the resolution in the experiment is restricted to low frequency components.

Fig. 6b also shows the results of the simulation using both the 1D-GPE and the NPSE (compare with Fig. 2 of the main text). We can see that for these atom numbers there are differences between the predictions of the 1D-GPE and NPSE. However these differences are small compared to the uncertainty in the experimental results.

IV Higher order solitons

Figure 7 shows numerical simulations of the 1D-GPE for the time evolution of second- and third-order solitons with initial sizes $l_{z}^{(2)}$ and $l_{z}^{(3)}$. Large initial soliton sizes lead to the periodic formation of local maxima and minima of the density profile. Striking characteristics of the time evolution are the periodic development of a sharp central peak with side wings for the second-order soliton (Fig. 7a,b), and the periodic formation of a broad double-peak structure for the third-order soliton (Fig. 7c,d).

We also simulate the time evolution of solitons with the same start conditions using the NPSE and analyse the results using the root mean square width of the wave packet for a quantitative comparison (Fig. 8)

[TABLE]

Here, $\bar{z}$ is the mean position of the wave packet and $n(z,t)$ is the 1D-density. We observe small quantitative differences between the two equations but the overall behaviour is very similar.

V Variational approach for the BREATHING FREQUENCY

In this section, we show how the longitudinal breathing frequency plotted in Fig. 2 of the main article can be determined from the variational ansatz for the soliton. For a cylindrical cigar-shaped potential the energy functional of Eq. 3 is given by Carr and Castin (2002); Parker et al. (2007)

[TABLE]

The energy of a soliton can be determined with a variational method using the following ansatz for the wave function

[TABLE]

where the transverse width $l_{r}$ and longitudinal width $l_{z}$ are the variational parameters Carr and Castin (2002); Parker et al. (2007). Combining Eqs. 10 and 11, and rescaling the variables by the transverse frequency $\omega_{r}$, provides an equation for the normalized energy of the soliton Carr and Castin (2002)

[TABLE]

with $\varepsilon_{GP}=E/\hbar\omega_{r}$, $\gamma_{r}=l_{r}/\sigma_{r}$, $\gamma_{z}=l_{z}/\sigma_{r}$, $\lambda=\omega_{z}/\omega_{r}$, $\alpha=Na/\sigma_{r}$, and $\sigma_{r}=\sqrt{\hbar/m\omega_{r}}$. We can simplify Eq. 12 for our system with weak interactions and strong transverse confinement by neglecting variations of the radial soliton size, i.e. $\gamma_{r}=1$. The energy minimum is found by calculating the zero-crossing of the first derivative of Eq. 12 with respect to $\gamma_{z}$

[TABLE]

where $\alpha=-|\alpha|=-\sqrt{\zeta}$. Eq. 13 has been solved for an expulsive potential with $\omega^{2}_{z}<0$ Carr and Castin (2002). Here, we provide the solution for a trapping potential with $\omega^{2}_{z}>0$. The longitudinal size of the soliton $\gamma_{z}^{*}$ at the energy minimum is

[TABLE]

with

[TABLE]

where

[TABLE]

with

[TABLE]

In order to find the oscillation frequency $\omega_{z}$ of the soliton, the equations of motion for the variational parameters are determined with a Lagrangian variational analysis Carr and Castin (2002)

[TABLE]

where the time derivative is calculated with respect to the normalised time $\tau=\omega_{r}t$. Again, we have assumed that the radial size of the soliton is constant, i.e. $\gamma_{r}=1$.

For small deviations of the soliton size from its equilibrium value, we can write the solution as $\gamma_{z}=\gamma^{*}_{z}+\delta\gamma_{z}$, where $\gamma^{*}_{z}$ is the minimum given by Eq. 14 and $\delta\gamma_{z}$ is a small deviation. A linear expansion of Eq. 18 leads to the expression for the longitudinal breathing frequency

[TABLE]

We compare our experimental measurements of the breathing frequency of the soliton to the predictions of Eq. 19 in Fig. 2 of the main article (red line).