# Generating soliton trains through Floquet engineering

P. Blanco-Mas  Departamento de Física de Materiales, Universidad Complutense
de Madrid, E-28040 Madrid, Spain  C.E. Creffield  Departamento de Física de
Materiales, Universidad Complutense de Madrid, E-28040 Madrid, Spain

###### Abstract

We study a gas of interacting ultracold bosons held in a parabolic trap in the
presence of an optical lattice potential. Treating the system as a discretised
Gross-Pitaevskii model, we show how Floquet engineering, by rapidly “shaking”
the lattice, allows the ground-state of the system to be converted into a
train of bright solitons by inverting the sign of the hopping energy. We study
how the number of solitons produced depends on the system’s nonlinearity and
the curvature of the trap, show how the technique can be applied both in the
high and low driving-frequency regimes, and demonstrate the phenomenon’s
stability against noise. We conclude that the Floquet approach is a useful and
stable method of preparing solitons in cold atom systems.

##  I Introduction

A soliton is a localized excitation in a medium that preserves its shape in
time as it travels. First observed in water waves propagating in narrow
channels Russell (1844), they are ubiquitous in Nature, arising in such
diverse contexts as laser pulses in optical fibers Hasegawa (2002), the
dynamics of tsunamis Constantin and Henry (2009), and kinks moving along DNA
Davydov (1973). Their stability arises from a balance between a localized
wavepacket’s intrinsic tendency to spread in time, and a non-linear attractive
interaction which opposes this spreading.

Bose-Einstein condensates (BECs) are particularly interesting candidates to
study soliton formation, as the systems are very clean and highly
controllable, and, in particular, the interatomic interaction can be
manipulated accurately in the laboratory. When the interaction is repulsive
the BEC is stable. However, if the interaction is rapidly changed to be
attractive the BEC is rendered vulnerable to modulational instability, causing
it to break up into ripples which are then focused by the interaction to form
solitons. In this way, soliton trains have been produced in BECs of lithium-7
atoms Strecker  _et al._ (2002); Khaykovich  _et al._ (2002), rubidium-85
Cornish, Thompson, and Wieman (2006) and in cesium-133 Mez̆nars̆ic̆  _et al._
(2019); Di Carli  _et al._ (2019), using the Feshbach resonance technique to
rapidly change the interaction strength. More recently, an alternative method
of controlling the atomic interaction has been developed which involves
coupling internal atomic states with different scattering lengths with an RF
field, which has been successfully used Sanz  _et al._ (2022) to produce
solitons in BECs of potassium-39.

In this work we revisit a protocol proposed by Carr and Brand (CB) Carr and
Brand (2004a, b) for converting a trapped BEC into a soliton train. It
consists of two parts. Initially a BEC with a repulsive interaction is held in
a trapping potential, as normal. This stable situation is then perturbed by
flipping the sign of the interatomic interaction, while simultaneously
inverting the trap so that it becomes expulsive. If we write the Hamiltonian
of the system as a sum of the kinetic energy $T$, the trapping potential $V$,
and the interaction term $U$, we can schematically represent this process as

$H=T+V+U\ \Longrightarrow T-V-U\ .$ (1)

This protocol thus requires control over both the interaction strength and the
sign of the trap, which may not always be experimentally feasible. We instead
propose to use a variant of this technique by addressing a single parameter:
the sign of the kinetic energy. This process can be represented instead as

$H=T+V+U\ \Longrightarrow-T+V+U\ ,$ (2)

which is clearly equivalent to the CB protocol, but with an overall minus
sign.

The required inversion of the kinetic energy, equivalent to endowing the atoms
with a negative effective mass, can be achieved by a technique known as
“Floquet engineering” Eckardt (2017) by applying an external driving field
which oscillates periodically in time. The dynamics of the system can then be
factored into a rapid micromotion oscillating at the same frequency as the
driving field, and an effective static Hamiltonian $H_{\mathrm{eff}}$. The
parameters of $H_{\mathrm{eff}}$ can be controlled very precisely by the
driving field, and in particular “shaking” a tight-binding lattice model by
rapidly oscillating the lattice in space, allows the intersite tunneling to be
coherently manipulated. The intersite tunneling, for example, can be tuned to
zero to produce the effect known as “coherent destruction of tunneling” (CDT)
Grossmann  _et al._ (1991), it can be rendered complex to allow the system to
simulate the effect of a synthetic magnetic field Aidelsburger  _et al._
(2013); Miyake  _et al._ (2013); Creffield  _et al._ (2016), or its sign can
be inverted Lignier  _et al._ (2007); Zenesini  _et al._ (2009), to provide a
negative effective mass.

Inverting the sign of the effective mass has been previously used in a static
system to produce gap solitons Eiermann  _et al._ (2004), by moving a trapped
condensate to the edge of the first Brillouin zone (FBZ) where the dispersion
relation has negative curvature. More recently, Ref. Mitchell  _et al._ (2021)
performed an experiment using Floquet engineering to flip the sign of the
tunneling by high-frequency shaking of a cesium BEC to obtain a solitonic
wavepacket at the center of the FBZ. By employing a similar driving of this
type we will show how Floquet engineering can be used in conjunction with
varying the trapping potential and the magnitude of the atomic interaction to
create stable soliton trains containing a specific number of solitons. We will
then go on to examine both the high and low frequency driving regimes, and
show that soliton formation occurs in both. Having the ability to vary the
driving frequency in this way gives the flexibility to avoid exciting the
atoms to a higher band at high driving frequencies, and also to evade
parametric resonances Lellouch  _et al._ (2017), which would otherwise heat
and eventually destroy the condensate. Finally we will study the effect of the
phase of the driving on the protocol, and demonstrate the scheme’s remarkably
high robustness to noise.

##  II Method

###  II.1 Model

A BEC held in a trap potential can be described well by the Gross-Pitaevskii
Hamiltonian

$H_{\mathrm{G}P}=\frac{-\hbar^{2}}{2m}\partial_{x}^{2}+V(x)+g\left|\psi(x)\right|^{2}\
,$ (3)

where $g$ is the interatomic coupling constant, proportional to the $s$-wave
scattering length, $m$ is the atomic mass, and the condensate wavefunction,
$\psi(x)$, is normalized to one. The trapping potential is denoted by $V(x)$,
and will be taken to be quadratic, $V(x)=V_{0}x^{2}$, where the distance $x$
is measured from the center of the system. If we now discretise space by
imposing an optical lattice potential, Eq. 3 can be rewritten in second-
quantized form,

$H_{\mathrm{latt}}=-J\sum_{j}\left(a_{j}a_{j+1}^{\dagger}+\mathrm{H.c.}\right)+\sum_{j}V(x_{j})n_{j}+g\sum_{j}n_{j}^{2}\
,$ (4)

where the kinetic energy is now given in terms of the tunneling $J$ between
nearest-neighbor lattice sites, $a_{j}/a_{j}^{\dagger}$ are the annihilation /
creation operators for a boson on site $j$, and $n_{j}=a_{j}^{\dagger}a_{j}$
is the standard number operator. Henceforth we will take $\hbar=1$, and
measure all energies and frequencies in units of $J$, and write all distances
in units of the lattice constant.

We will take the initial state for the simulations to be the ground state of
Hamiltonian (4). This is obtained by starting with the solution for $g=0$ (a
Gaussian) and evolving it under Eq. 4 in imaginary time, maintaining the
correct normalization, until convergence is achieved. The results for several
different values of $g$ are shown in Fig. 1, and the evolution from a narrow
Gaussian distribution to the broader inverted parabolic form predicted by the
Thomas-Fermi approximation in the limit of large $g$ is clearly visible.

Figure 1: Probability density of the ground state of the lattice GPE (Eq. 4)
for different values of the non-linearity $g$. The curvature of the trap is
given by $V_{0}=2\times 10^{-5}$. Note how the wavefunction evolves from a
Gaussian for $g=0$, to a broader, more flattened form as $g$ increases and the
system approaches the Thomas-Fermi limit.

We now wish to introduce the time-dependent driving potential. If the optical
lattice is periodically oscillated in space, or “shaken”, an observer in the
rest frame of the lattice perceives an inertial force described by a time-
dependent lattice tilt, resulting in the Hamiltonian

$H(t)=H_{\mathrm{latt}}+K(t)\sum_{j}j\ n_{j}\ ,$ (5)

where $K(t)$ describes the form of the shaking. For the common case of
sinusoidal shaking, this function can be written as $K(t)=K\cos\left(\omega
t+\phi\right)$, where $\omega$ is the frequency of the driving, $K$ is its
amplitude, and for generality we have included a driving-phase $\phi$. All
numerical results were obtained by first preparing the initial state using the
imaginary time relaxation method described above, and then integrating this
state in time under Eq. 5, using a fourth-order Runge-Kutta routine.

###  II.2 Floquet engineering

As the time-dependent driving (5) is $T$-periodic, $H(t)=H(t+nT)$, solutions
of the time-dependent Schrödinger equation are of Floquet form

$\left[i\partial_{t}-H(t)\right]\psi_{n}(t)=\epsilon_{n}\psi_{n}(t)$ (6)

where the Floquet states $\psi_{n}(t)$ have the same $T$-periodicity as the
Hamiltonian, and provide a complete basis to describe the time evolution of
the driven system. The eigenvalues $\epsilon_{n}$ are called the
quasienergies, and govern the long-term dynamics of the system, as the time-
dependence of the Floquet states only operates over timescales within each
driving period, providing the so-called “micromotion” Eckardt (2017) of the
system.

The quasienergies can be obtained as the eigenenergies of an effective static
Hamiltonian, $H_{\mathrm{eff}}$, which depends on the parameters and form of
the driving. The process of Floquet engineering then consists of choosing the
appropriate driving function to produce the desired form of
$H_{\mathrm{eff}}$, and thus the quasienergies. In general it is difficult to
obtain closed-form solutions for the quasienergies for a given drive. It is
possible, however, to obtain expressions for $H_{\mathrm{eff}}$ as series
expansions in inverse frequency, such as the Magnus series Blanes  _et al._
(2009) and the Van Vleck series Eckardt and Anisimovas (2015), which become
exact in the limit of infinite driving frequency. To first order, it can be
shown that for a sinusoidally-driven two-level system the quasienergies are
given by Creffield (2003) $\epsilon_{\pm}=\pm J\ {\cal
J}_{0}\left(K/\omega\right)$, where ${\cal J}_{0}$ is the zeroth Bessel
function of the first kind. Fig. 2 shows the excellent agreement between this
result and the numerical results for a two-level system driven at $\omega=16$.
Good agreement continues to be obtained as long as $\omega>J$: the high-
frequency regime. However, for lower frequencies the quasienergies behave
differently, as can be seen for the case of $\omega=1$, indicating that more
terms Barata and Wreszinski (2000) must be included in the series.

The effective tunneling between the levels is proportional to the difference
between the quasienergies,
$J_{\mathrm{eff}}=\left(\epsilon_{+}-\epsilon_{-}\right)/2$, and thus in the
high-frequency limit the effective tunneling is given by

$J_{\mathrm{eff}}=J\ {\cal J}_{0}\left(K/\omega\right)\ .$ (7)

Altering the parameter $K/\omega$ therefore allows $J_{\mathrm{eff}}$ to be
tuned to a desired value. In particular, if we set $K/\omega=2.404$ (the first
zero of the Bessel function) the two quasienergies become degenerate, as can
be seen in Fig. 2a, and the effective tunneling vanishes, producing CDT.
Increasing $K/\omega$ beyond this value causes $J_{\mathrm{eff}}$ to become
negative in the interval between the first and second zeros of ${\cal J}_{0}$
(the shaded region in Fig. 2 b), which is thus the region of interest for our
soliton generation method (2).

Figure 2: (a) Quasienergies for the two-level model. Black solid lines show
the first-order approximation $\epsilon=\pm{\cal J}_{0}\left(K/\omega\right)$
which becomes exact in the limit $\omega\rightarrow\infty$. Black circles show
the exact numerical results for $\omega=16$; at this high-frequency the
results are excellently approximated by the perturbative result. The red
triangles show the quasienergies for $\omega=1$, which differ considerably
from the Bessel function behaviour of the high frequency results. b)
$J_{\mathrm{eff}}$ for $\omega=16$; note how this vanishes at
$K/\omega=2.404$, and is negative in the shaded region between the first and
second zeros of ${\cal J}_{0}$. The arrow marks $K/\omega=3.80$, the value
used in the Floquet engineering method. c) $J_{\mathrm{eff}}$ for $\omega=1$.
Again, the shaded region indicates where the effective tunneling is negative.
The quantities $\epsilon$ and $J_{\mathrm{eff}}$ are measured in units of $J$
(see text).

Proceeding to an $N$-site lattice system, the quasienergies now present a
band-like structure Creffield (2009)

$\epsilon_{n}=-2J_{\mathrm{eff}}\cos k_{n}\ ,$ (8)

where $k_{n}$ are the permitted momenta in the FBZ, and $J_{\mathrm{eff}}$ is
the same intersite tunneling discussed above (Eq. 7). Accordingly we can
regulate the effective mass of a particle moving in the lattice in the same
way that we can control the tunneling in the two-level model, by appropriately
choosing the value of $K/\omega$.

##  III Results

Soliton generation –  In Fig. 3 we show the results of applying the Floquet
protocol to a system with nonlinearity parameter $g=2$. The system was
prepared in its ground state, and then at $t=0$ the driving potential was
suddenly turned on, with phase $\phi=0$, and a frequency of $\omega=16$
placing the system firmly in the high-frequency regime. As shown in Fig. 2b,
the amplitude of the driving was chosen to be $K/\omega=3.80$ that is, the
first minimum of Eq. 7, the corresponding value of $J_{\mathrm{eff}}=-0.403$
being the largest negative value possible for the effective tunneling.

Figure 3: Soliton generation using the Floquet engineering protocol (2).
Physical parameters: 512 lattices sites, $V_{0}=2\times 10^{-5}$, $g=2$, and
the driving frequency $\omega=16$. Left: contour plot of the condensate
density $|\psi(x,t)|^{2}$. The system is initialized in its ground state, and
at $t=0$ the driving potential is applied. At $t=100$ ripples begin to form in
the condensate, which gradually focus into three solitons. The central soliton
remains close to the center of the lattice, remaining in unstable equilibrium
with the lattice potential, while the other two solitons accelerate
exponentially away. Right: cross sections through the condensate at different
times. The ripples form initially at the edges of the wavepacket, indicating
they arise from self-interference, and then sharpen into the typical soliton
shape.

When the driving potential is applied, small ripples develop in the profile of
the condensate. As Carr and Brand demonstrated, the dominant source of these
ripples is self-interference of the condensate wavefunction, which causes them
to develop first at the edges of the wavepacket. Modulational instability
causes the ripples to grow in amplitude, and they are then focused by the
nonlinearity into forming sharp peaks with a typical soliton profile. For this
value of the nonlinearity the initial wavepacket divides into three solitons.
The one which forms near the center of the trap remains essentially stationary
over the remainder of the time evolution, while the other two accelerate away
at an exponentially increasing velocity. Although the trapping potential does
not change during the protocol, the solitons have a negative effective mass
due to the inversion of $J_{\mathrm{eff}}$, and so “fall uphill” against the
potential and are accelerated outwards, just like a normal particle in an
inverted potential. Note that once formed, the solitons are stable, and retain
their form throughout their trajectory.

Raising the value of $g$ reveals that the number of solitons produced in a
trap with a given curvature increases weakly as a function of the
nonlinearity, as observed previously in Refs. Nguyen, Luo, and Hulet (2017);
Everitt  _et al._ (2017). The trap curvature can also be used as a parameter
to control the soliton number. We show the combined effect of these factors in
Fig. 4a. From this plot it is clear that high soliton numbers are favored by a
high value of the nonlinearity and a low trap curvature. This can be
understood qualitatively by a simple scaling argument. The growth of the
ripples is governed by the most unstable Bogoliubov mode, the wavelength of
which Nguyen, Luo, and Hulet (2017) is related to the healing length of the
condensate, $\xi$. The other length scale of the problem is the harmonic
oscillator length $\ell$ which gives an estimate for the effective width of
the initial wavepacket. On dimensional grounds, the number of solitons will
vary approximately as the ratio of these lengths

$n\sim\xi/\ell\sim\left(g/\sqrt{V_{0}}\right)^{1/2}\ ,$ (9)

in agreement with the trends observed.

Figure 4: a) Number of solitons produced as the trap curvature $V_{0}$ and the
nonlinearity $g$, both measured in units of $J$, are varied. The number of
solitons increases for high $g$ and low $V_{0}$, as predicted by Eq. 9. In all
cases a lattice of 256 sites was used. b) Plot of the critical interaction
strength, $g_{c}$, at which the system changes from producing two solitons to
producing three. The dashed line is a least squares fit to
$g_{c}=a_{1}+a_{2}\sqrt{V_{0}}$, where $a_{1}$ and $a_{2}$ are fitting
parameters, which describes the behaviour with reasonable accuracy.

To examine this behaviour more quantitatively, we plot in Fig. 4b the critical
value of the nonlinearity parameter, $g_{c}$, at which the number of solitons
produced in the system changes from two to three – that is, the lowest
boundary curve plotted in Fig. 4a. According to Eq. 9, this quantity should
vary as $g_{c}\sim\sqrt{V_{0}}$. As can be seen, the data can indeed be fitted
with reasonable accuracy by this expression.

Low-frequency regime –  The results presented so far have been for
$\omega=16$, which is well within the high-frequency regime. The results
obtained are essentially identical to those obtained by performing the
simulations without the time-dependent driving, but setting the value of $J$
to $J_{\mathrm{eff}}=-0.403$ by hand at $t=0$, which we shall term the
“switched protocol”. The excellent agreement between the results indicates how
well the Floquet protocol duplicates the switched protocol for this frequency.

As $\omega$ is lowered we should expect this agreement to reduce, as more
terms must be included in the Magnus expansion, and $H_{\mathrm{eff}}$ can no
longer be approximated as a nearest-neighbor hopping Hamiltonian with a single
effective tunneling. Perhaps surprisingly, although the amplitude and position
of the generated solitons alter slightly for smaller values of $\omega$, the
process of soliton formation itself remains robust. The important point is
just that the nearest-neighbor tunneling must become negative. Although the
positions of quasienergy degeneracies will drift away from the zeros of the
Bessel function as the frequency is reduced Creffield (2003), there will
nonetheless still be some intervals over which $J_{\mathrm{eff}}$ will be
negative.

As an example, in Fig. 2c, we plot the effective tunneling for a system at a
low driving frequency of $\omega=1$. As we noted previously,
$J_{\mathrm{eff}}$ is related to the difference of the quasienergies, but
unlike the high-frequency case, in this instance we do not know which
quasienergy corresponds to $\epsilon_{+}$ and which to $\epsilon_{-}$. As a
result, even knowing the values of the quasienergies we are uncertain of the
sign of $J_{\mathrm{eff}}$. Simulating the driven system in a lattice with an
additional static tilt reveals that an initial wavepacket moves up the
potential for $K/\omega<4.8$, indicating that its effective mass is negative,
while it moves down the potential for higher driving amplitudes, corresponding
to its effective mass being positive. Accordingly, driving the system at
$K/\omega=3.8$, as done in the high-frequency case, should also produce
solitons.

We show the results in Fig. 5 for three different driving frequencies. At the
sample time used, $t=300$, the switched protocol produces in three solitons.
For a driving frequency of $\omega=16$, the result is essentially
indistinguishable for the switched protocol. Lowering the driving frequency to
$\omega=4$ again gives a very similar result: although the three solitons
produced have slightly different positions and amplitudes, the differences
from the $\omega=16$ result are very minor. Taking now the value $\omega=1$,
well outside the high-frequency regime, again gives a very similar result.
Using a high driving frequency is thus not a necessary requirement for this
technique, as long as $J_{\mathrm{eff}}$ changes sign.

Figure 5: Soliton development at $t=300$ for three different driving
frequencies. In the high frequency regime, $\omega=16$ (solid black line), the
results obtained are indistinguishable from the switched protocol. For
$\omega=4$ (dash-dotted red line) the solitons produced strongly resemble the
high frequency result, although small differences are visible. At low
frequency, $\omega=1$ (dashed blue line), the differences are larger, but
qualitatively a very similar result is still obtained. System parameters: 512
lattice sites, $V_{0}=0.00002$ and $g=2$.

Effect of phase –  The amplitude of the driving is set through the requirement
$K/\omega=3.8$, and we have shown that the effect is robust to varying the
driving frequency, $\omega$. This leaves one last parameter in the driving to
consider: the phase, $\phi$. We have so far used cosinusoidal driving, that
is, $\phi=0$. If we instead use sinusoidal driving ($\phi=\pi/2$) we see no
sign of any soliton creation. The wavepacket instead sloshes from side to side
of the trapping potential. This effect is produced by the so-called kick-
operator Eckardt (2017). As shown in Ref. Creffield and Sols (2008), when
$\phi$ is not zero, the driving imprints a phase onto the condensate which
excites it into motion. For the type of driving we consider the full
expression for the effective tunneling is given in the high-frequency limit by

$J_{\mathrm{eff}}=Je^{-i\left[K/\omega\right]\sin\phi}{\cal
J}_{0}\left(K/\omega\right)\ ,$ (10)

which clearly reduces to Eq. 7 for $\phi=0$. While the same Bessel function
renormalization of the modulus of the tunneling still happens for sinusoidal
driving, the additional phase factors mean that $J_{\mathrm{eff}}$ does not
flip sign when $K/\omega=3.80$. The driving instead delivers a kick to the
condensate, giving it an initial velocity $v_{0}=K/\omega$, and thus causing
it to make harmonic oscillations about the center of the trap.

Noise –  We finally consider how stable the soliton creation process is to
noise on the initial state, since in experiment it is clearly impossible to
prepare the desired initial state with complete fidelity. We have so far used
the ground state, $\psi^{0}(x)$, as the initial state, and we will now add
random noise to it

$\psi(x,t=0)=\psi^{0}(x)+r(x)\gamma\ ,$ (11)

where $r$ is a random variable uniformly distributed over $(-1,1)$, and
$\gamma$ sets the amplitude of the noise. After the noise has been added to
both the real and imaginary components of $\psi^{0}(x)$, the resulting state
is normalized to unity as usual, and used as the initial state of the
simulation.

We show the results for two different noise levels in Fig. 6. For the given
system parameters, $V=2\times 10^{-5}$ and $g=2$, we expect to produce three
solitons, as seen previously in Fig. 3. For $\gamma=0.01$ (Fig. 6a) the
initial state already appears notably irregular, but the process indeed gives
rise to three well defined solitons. Unlike in the clean system, the seeds for
soliton formation are now dominated by the imposed random fluctuations in the
condensate wavefunction, rather than by self-interference of the condensate.
As a result the positions of the solitons are randomly shifted in position
with respect to the clean system, and the solitons are more equal in size, as
they began forming at essentially the same time. Further increasing the noise
level to $\gamma=0.05$ produces two solitons instead of the three expected.
Nonetheless it is striking that even in the presence of such a high level of
noise, the process still produces clearly identifiable solitons that are
easily detectable above the background noise level

Figure 6: Effect of noise on soliton generation. (a) Above: initial state with
$\gamma=0.01$. Below: the particle density of the system at $t=300$. Just as
in the clean system (Fig. 3), three solitons form. (b) Above: initial state
with $\gamma=0.05$. Below: at $t=300$ the system has evolved to present just
two solitons, the random noise having affected details of the soliton
generation process. Physical parameters: 512 lattice sites, $V=2\times
10^{-5}$, $g=2$.

##  IV Conclusions

In summary, we have shown how Floquet engineering may be used to implement a
protocol based on the CB method for converting a trapped condensate into a
train of solitons. It has the advantage that only one control parameter needs
to be altered, and as it does not involve controlling the interparticle
interaction, it is applicable to atomic species for which this control is not
easily available, such as those like rubidium-87 which lack a convenient
Feshbach resonance. The method consists of using Floquet physics to invert the
sign of the intersite tunneling, or equivalently, to give the condensate atoms
a negative effective mass. Modulational instability then causes the initial
state to break up into a train of spatially-localized pulses, which then self-
focus into solitons under the influence of the nonlinear interaction. We have
demonstrated how the curvature of the trap and the magnitude of the atomic
interaction can be used to deterministically prepare trains of a given number
of solitons. Unlike many applications of Floquet engineering, the method is
not restricted to the high-frequency regime, but works equally well for lower
driving frequencies, which have the advantage of avoiding driving the atoms
into higher bands. In experiment, soliton generation in the low-frequency
regime may be harder to attain due to the increased rate of heating Lellouch
_et al._ (2017) destroying the coherence of the condensate. Nonetheless,
experiments Lignier  _et al._ (2007) have reported coherent control of the
tunneling down to frequencies of $\omega=0.3J$, well into the low frequency
regime, and once formed the solitons are able to self-cool Carr and Brand
(2004a) by emitting small bursts of atoms. Finally we would like to emphasise
the remarkable stability of the method with respect to noise. The modulational
instability is seeded by the most unstable Bogoliubov mode, and as this
typically has a fairly long wavelength, noise on shorter length scales has
relatively little effect. The repeatability and excellent control afforded by
this method make it an excellent tool to investigate soliton dynamics and
collisions, and hold out the prospect of using such soliton trains for
precision measurement applications such as atom interferometry.

###### Acknowledgements.

This work was supported by the Universidad Complutense de Madrid through grant
no. FEI-EU-19-12.

## References

  * Russell (1844) J. S. Russell, Report of the British Association for the Advancement of Science , 311 (1844).
  * Hasegawa (2002) A. Hasegawa, Opt. Photon. News 13, 33 (2002).
  * Constantin and Henry (2009) A. Constantin and D. Henry, Zeitschrift für Naturforschung A 64, 65 (2009).
  * Davydov (1973) A. Davydov, Journal of Theoretical Biology 38, 559 (1973).
  * Strecker  _et al._ (2002) K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. Hulet, Nature 417, 150 (2002).
  * Khaykovich  _et al._ (2002) L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, Science 296, 1290 (2002).
  * Cornish, Thompson, and Wieman (2006) S. L. Cornish, S. T. Thompson, and C. E. Wieman, Phys. Rev. Lett. 96, 170401 (2006).
  * Mez̆nars̆ic̆  _et al._ (2019) T. Mez̆nars̆ic̆, T. Arh, J. Brence, J. Pis̆ljar, K. Gosar, Z̆. Gosar, R. Z̆itko, E. Zupanic̆, and P. Jeglic̆, Phys. Rev. A 99, 033625 (2019).
  * Di Carli  _et al._ (2019) A. Di Carli, C. D. Colquhoun, G. Henderson, S. Flannigan, G.-L. Oppo, A. J. Daley, S. Kuhr, and E. Haller, Phys. Rev. Lett. 123, 123602 (2019).
  * Sanz  _et al._ (2022) J. Sanz, A. Frölian, C. S. Chisholm, C. R. Cabrera, and L. Tarruell, Phys. Rev. Lett. 128, 013201 (2022).
  * Carr and Brand (2004a) L. D. Carr and J. Brand, Phys. Rev. A 70, 033607 (2004a).
  * Carr and Brand (2004b) L. D. Carr and J. Brand, Phys. Rev. Lett. 92, 040401 (2004b).
  * Eckardt (2017) A. Eckardt, Rev. Mod. Phys. 89, 011004 (2017).
  * Grossmann  _et al._ (1991) F. Grossmann, T. Dittrich, P. Jung, and P. Hänggi, Phys. Rev. Lett. 67, 516 (1991).
  * Aidelsburger  _et al._ (2013) M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Phys. Rev. Lett. 111, 185301 (2013).
  * Miyake  _et al._ (2013) H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton, and W. Ketterle, Phys. Rev. Lett. 111, 185302 (2013).
  * Creffield  _et al._ (2016) C. E. Creffield, G. Pieplow, F. Sols, and N. Goldman, New Journal of Physics 18, 093013 (2016).
  * Lignier  _et al._ (2007) H. Lignier, C. Sias, D. Ciampini, Y. Singh, A. Zenesini, O. Morsch, and E. Arimondo, Phys. Rev. Lett. 99, 220403 (2007).
  * Zenesini  _et al._ (2009) A. Zenesini, H. Lignier, D. Ciampini, O. Morsch, and E. Arimondo, Phys. Rev. Lett. 102, 100403 (2009).
  * Eiermann  _et al._ (2004) B. Eiermann, T. Anker, M. Albiez, M. Taglieber, P. Treutlein, K.-P. Marzlin, and M. K. Oberthaler, Phys. Rev. Lett. 92, 230401 (2004).
  * Mitchell  _et al._ (2021) M. Mitchell, A. Di Carli, G. Sinuco-León, A. La Rooij, S. Kuhr, and E. Haller, Phys. Rev. Lett. 127, 243603 (2021).
  * Lellouch  _et al._ (2017) S. Lellouch, M. Bukov, E. Demler, and N. Goldman, Phys. Rev. X 7, 021015 (2017).
  * Blanes  _et al._ (2009) S. Blanes, F. Casas, J. Oteo, and J. Ros, Physics Reports 470, 151 (2009).
  * Eckardt and Anisimovas (2015) A. Eckardt and E. Anisimovas, New Journal of Physics 17, 093039 (2015).
  * Creffield (2003) C. Creffield, Phys. Rev. B 67, 165301 (2003).
  * Barata and Wreszinski (2000) J. C. A. Barata and W. F. Wreszinski, Phys. Rev. Lett. 84, 2112 (2000).
  * Creffield (2009) C. E. Creffield, Phys. Rev. A 79, 063612 (2009).
  * Nguyen, Luo, and Hulet (2017) J. H. V. Nguyen, D. Luo, and R. G. Hulet, Science 356, 422 (2017).
  * Everitt  _et al._ (2017) P. J. Everitt, M. A. Sooriyabandara, M. Guasoni, P. B. Wigley, C. H. Wei, G. D. McDonald, K. S. Hardman, P. Manju, J. D. Close, C. C. N. Kuhn, S. S. Szigeti, Y. S. Kivshar, and N. P. Robins, Phys. Rev. A 96, 041601 (2017).
  * Creffield and Sols (2008) C. E. Creffield and F. Sols, Phys. Rev. Lett. 100, 250402 (2008).