Excitation modes of bright matter-wave solitons

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.

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 [1].In fiber optics, one-dimensional (1D) "bright" solitons, i.e. solitons presenting a local electric field maximum with one-dimensional propagation, have been observed [2].They exhibit a dispersionless flow and excitation modes such as breathing or higher-order modes [2][3][4].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 [5].Recent experiments demonstrated the collapse [6], collisions [7], reflection from a barrier [8], and the formation of trains [9][10][11] 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 [7], excitations following the collapse of attractive BECs [6,12], and quadrupole oscillations of attractive BECs in three dimensions (3D) [13].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 [3,14].
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) [15].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 [16]).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) [17].We reference to the analytical solutions of the 1D-GPE in the manuscript, but use both equations in our numerical simulations [16].
For the 1D-GPE, an ansatz for the normalized longitudinal wave function f (z) is of the form 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 [18].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 [14,19] or numerically [18].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 ω r , the size l z corresponds to the healing length at the peak density of the soliton, i.e. l (n=1) z = a 2 r /(N|a|) [14,18].Here, a r = h/(mω 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 |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 [20].Details about our experimental setup, the levitation scheme and the removal of atoms can be found in refer- ences [16,21].
Our matter-wave solitons are confined to a quasi-1D geometry with almost free propagation along the horizontal direction and strong radial confinement of ω r = 2π × 95 Hz provided by laser beam L H .They are generated with a quench of the scattering length towards attractive interaction (a i → a f ), and by a reduction of the longitudinal trap frequency (ω z,i → ω z, f ).When changing a and ω 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 [16]).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(sech(z/B)) 2  to the integrated 1D-density profiles with fit parameters A and B [16].
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 ω z, f = 2π × 4.3(2) Hz (quench Q1 in [16]) while keeping the repulsive interaction strength constant (Fig. 1d, diamonds).The BEC starts an outwards motion with an oscillation frequency of 2π × 7.5(1) Hz ≈ √ 3ω z, f as expected for a BEC in the Thomas-Fermi regime [22,23].In a second step, we additionally quench the interaction strength a f to −5.4 a 0 and increase ω 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) µm/s (Fig. 1d, green line).Finally, we deliberately mismatch the initial size of the BEC by reducing ω z,i (Q3), and generate small amplitude oscillations of the soliton with a frequency ω sol of 2π × 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ω z, f = 2π × 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, ω 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 ω sol without changing the quench protocol (Q4, Fig. 2a, circles).The measured values of ω sol decrease for lower N, and they approach the breathing frequency 2ω 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 [16] (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 [14,16]).We find that both models agree well with the trend of the measurements of ω 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 ω sol as we reduce the longitudinal trapping frequency ω z, f (Q5).Two regimes of ω sol can be identified in Fig. 2b for varying the values of ω z, f .For large values of ω z, f , the trap dominates the breathing of the soliton and ω sol approaches twice the trap frequency 2ω z, f .For small values of ω z, f , interactions dominate the breathing of the soliton and ω 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 [7,24,25], causing, e.g., modulations of the soliton's tails due to residual non-autonomous terms of the 1D-GPE in a harmonic potential [26].For the following experiments, however, we employ trap frequencies that are significantly smaller than the observed oscillation frequencies of the soliton (2ω z < ω sol ) to decouple the influence of the trapping potential.In summary, for small-amplitude oscillations we find good agreement of ω sol between our experimental results and analytical and numerical predictions based on the 1D-GPE (and NPSE [16]).
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) [3], and has been observed for optical solitons in silica-glass fibers [2,4].A soliton of order n can be interpreted as a bound state of n strongly overlapping solitons [14].By exploiting the equivalence of the NSE and 1D-GPE, similar effects were later proposed for bright matter-wave solitons [14,27], where it was suggested that n th -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 (n) z that is the n 2 multiple of the healing length l (1) z [16].Within the 1D-GPE theory, both creation methods result in the periodic development of multi-peaked structures for higher-order solitons [3,28], 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 [16].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 [3].For a second-order soliton, the predicted oscillation period T (2) is [14] Recently, excitation modes of higher-order have also been used as a testbed for various theoretical models beyond GPtheory.The fragmentation of solitons with an increased initial width was predicted within the multiconfigurational timedependent Hartree method for bosons [29] and critically discussed [30], and the influence of quantum effects on the dissociation process was investigated [31][32][33].
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 ω z,i before ramping a and ω z to −5.3 a 0 and 2π × 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π × 2.4(2) Hz (Fig. 3c).This frequency is significantly smaller than the expected breathing frequency of first-order solitons, 2π × 6.0 Hz, and it matches well to the expected frequency of 2π × 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 [3], 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 [14], we implement a double-quench protocol, with a first quench to generate a soliton with weak attractive interaction, a f = −0.8 a 0 , ω f = 2π × 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π × 5.6( 6) Hz (Fig. 4b), which matches the expected frequency of 2π × 5.2 Hz for second-order solitons (2π × 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 [21].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 ω sol for first-order solitions and its dependence on N and ω z .For stronger excitations and for a double-quench protocol, we observed signatures of secondorder 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 [31][32][33].

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 ω r = 2π × 95 Hz provided by laser beam L H .In quasi-1D geometry, bright matter-wave solitons collapse for large densities and interactions [6], which for our typical experimental scattering length of approximately −5 a 0 corresponds to a critical atom number of 2500 [14].As a result, we need to strongly reduce the atom number to avoid collapse, modulation instabilities [9] and three-body loss [34] 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 ±100 for 600 atoms and ±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%.Q5 Quench to determine the dependence of ω sol on the trap frequency ω z, f .We vary ω z, f from approximately 1 Hz to 9 Hz.Smaller values of ω z, f result in larger equilibrium sizes of the soliton, and we need to reduce the initial trap frequencies ω z,i to keep the oscillation amplitudes comparable during the measurements.The typical difference between ω z,i and ω z, f is approximately 3 Hz.N ≈ 1500, a f = −5.4 a 0 , ω r = 2π×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 , ω z,i = 2π × 4.9(2) Hz, ω z, f = 2π × 0.0(6) Hz, initial atom number N ≈ 3000 drops to 1100 after shedding of atoms, quench duration 13 ms, ω r = 2π×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 2: reduce interaction strength in 2 ms to a f = −4.6 a 0 , no change of other parameters.

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(sech(z/B)) 2 , with fitparameters A and B, to the integrated 1D-density profiles (dotted blue line in Fig. 5).

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, ψ(r,t), is given by, where g = 4π h2 a/m, m is the atomic mass, and a is the twobody 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 1 is satisfied [35], where n is the particle density.Our external potential V (r) is given by a 3D (anisotropic) harmonic trap.
For tight radial trapping potentials, ω r ω 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, where a r is the harmonic oscillator length in the radial direction and σ (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) [17] ih where ω r = h/ma 2 r .The condition for σ (z,t) that minimizes the action functional integrated along the trajectories in phase space is [17], For σ (z,t) = 1, we obtain the ground state of a harmonic oscillator in the radial directions, and we recover the usual 1D-GPE We have numerically integrated Eqs 5 and 7 using the splitstep Fourier transform method [36], where we exploit the fact that the kinetic and potential terms in the Hamiltonian are diagonal in momentum and real space, respectively.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,

SOLITON BREATHING FREQUENCY
where l z = a 2 r /(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 ω z = 2π × 5 Hz and atom number N = 1300, which is characteristic of the behaviour for all other ω 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.

FIG. 1 .
FIG. 1. Experimental setup and oscillation measurements.(a) Sketch of the experimental setup.Inset: density profiles for a BEC (solid red line) and for a soliton (dashed blue line).(b) Total energy of a soliton, a = −5.2 a 0 , ω r = 2π × 95 Hz, N = 2000, with an external trap, ω z = 2π × 5 Hz (dashed blue line), and without external trap, ω z = 0 Hz (solid red line).(c) Absorption images after a free expansion time of 16 ms (from data set with circles in d), integrated density profile for t = 60 ms (blue line) and fit (dashed red line).(d) Oscillations of a quantum gas after the steps of a quench procedure.Blue diamonds: quench of only ω z for a BEC (Q1), Red circles: additional interaction quench to create soliton (Q3), Green squares: optimized quench parameters to minimize breathing of the soliton (Q2).Uncertainty intervals indicate ±1 standard error.

FIG. 2 .
FIG. 2. Breathing frequency ω sol of the soliton.(a) Atom number dependence (Q4).Red circles: experimental data, the uncertainty bars for the atom number indicate the standard deviation of N over the first 100 ms of each frequency measurement.Blue triangles: simulation of the 1D-GPE [16].Red line: analytical approximation [14, 16].Dashed gray line: oscillation frequency of a non-interacting gas, 2ω z, f .(b) Dependence of ω sol on the trap frequency (Q5).Red circles: experimental data points for N ≈ 1450.Blue area: simulation of the 1D-GPE for N = 1300 to 1500.Red line: analytical approximation.Dashed gray line: 2ω z, f .

FIG. 3 .
FIG. 3. Time evolution after a strong quench of interactions and trap frequency (Q6).(a) Absorption images at time t after the quench and after 11 ms of free expansion.(b) 1D-GPE simulation of the density profiles for a second-order soliton with 1100 atoms, a f = −5.3 a 0 , and with an oscillation period T (2) of 432 ms.(c) Time evolution of the measured width l z of the central wave packet (red circles), sinusoidal fit with period 420(30) ms (dashed red line).The uncertainty intervals indicate ±1 standard deviation.

FIG. 4 .
FIG. 4. Second-order soliton and splitting after the double quench Q7.(a) Absorption images at time t after the quench and after 7 ms of free expansion.(b) Time evolution of the measured width l z of the central wave packet (red circles), sinusoidal fit with period 180(20) ms (dashed red line).The expected period from the 1D-GPE simulations is 192 ms.The histogram counts the fraction of images showing a splitting of the wave function (9 repetitions per time step).Inset: absorption image of a split matter wave for t = 210 ms.
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 ω z,i = 2π × 5.8(2) Hz to ω z, f = 2π × 4.3(2) Hz with a linear ramp of the laser power of beam L V over 4 ms.Atom number N ≈ 1800, constant interaction strength a i = +7 a 0 , ω r = 2π×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 ω z,i = 2π × 11.2(2) Hz.Atom number N ≈ 1800, ω z, f = 2π × 4.3(2) Hz, ω r = 2π×95 Hz.Q3 We mismatch the initial size of the BEC before the quench with ω z,i = 2π × 12.8(4) Hz to generate small amplitude oscillations of the width of the soliton.Atom number N ≈ 1700, ω z, f = 2π × 4.3(2) Hz, ω r = 2π×95 Hz.Q4 Quench to determine the atom-number dependence of ω sol .We vary the atom number N from 500 to 1700 for the measurement.ω z, f = 2π × 4.3(2) Hz, a i = +7 a 0 , a f = −5.4 a 0 , ramp duration 4 ms, ω r = 2π×95 Hz.

FIG. 5 .
FIG. 5. 1D-density profile of a soliton.Red line: integrated density profile of the absorption image for t = 60 ms in Fig. 1c (main article).Dotted blue line: fitted profile according to Eq. 1 in the main article.

FIG. 6 .
FIG. 6. Simulation results for the soliton breathing frequency, for comparison with Fig. 2 in the main text.(a) Frequency spectrum calculated using the 1D-GPE with a longitudinal frequency of ω z = 2π × 5 Hz and atom number N = 1300.(b) Breathing frequency (first peak in the spectrum as in a)) vs. trap frequency.N ≈ 1300 − 1500 atoms, a f = −5.4 a 0 for the NPSE (red) and the GPE (green).The simulations were evolved in time to 4000 ms.