Soliton trains after interaction quenches in Bose mixtures

André Cidrim Departamento de Física, Universidade Federal de São Carlos, 13565-905 São Carlos, Brazil. Luca Salasnich Dipartimento di Fisica e Astronomia “Galileo Galilei”, Università di Padova, via Marzolo 8, 35131 Padova, Italy INFN, Sezione di Padova, via Marzolo 8, 35131 Padova, Italy CNR-INO, via Nello Carrara, 1 - 50019 Sesto Fiorentino, Italy Padua Quantum Technologies Research Center, University of Padova, Via Gradenigo 6/b, 35131 Padova, Italy Tommaso Macrì Departamento de Física Teórica e Experimental, Universidade Federal do Rio Grande do Norte, and International Institute of Physics, Natal-RN, Brazil

Abstract

We investigate the quench dynamics of a two-component Bose mixture and study the onset of modulational instability, which leads the system far from equilibrium. Analogous to the single-component counterpart, this phenomenon results in the creation of trains of bright solitons. We provide an analytical estimate of the number of solitons at long times after the quench for each of the two components based on the most unstable mode of the Bogoliubov spectrum, which agrees well with our simulations for quenches to the weak attractive regime when the two components possess equal intraspecies interactions and loss rates. We also explain the significantly different soliton dynamics in a realistic experimental homonuclear potassium mixture in terms of different intraspecies interaction and loss rates. We investigate the quench dynamics of the particle number of each component estimating the characteristic time for the appearance of modulational instability for a variety of interaction strengths and loss rates. Finally we evaluate the influence of the beyond-mean-field contribution, which is crucial for the ground-state properties of the mixture, in the quench dynamics for both the evolution of the particle number and the radial width of the mixture. In particular, even for quenches to strongly attractive effective interactions we do not observe the dynamical formation of solitonic droplets.

Keywords: multi-component BECs, quench dynamics, modulational instability, solitons

1 Introduction

Modulational instability (MI) is a generic phenomenon that consists of the spontaneous exponential growth of perturbations resulting from the interplay between nonlinearity and anomalous dispersion. MI occurs in several areas of physics. It has been observed in classical [1, 2] and quantum [3, 4, 5, 6, 7] fluids, in waveguides [8] and in lattices [9], as well as in nonlinear optics [10] and in charged plasmas [11].

In ultracold Bose-Einstein condensates several experimental [3, 4, 5, 6, 7] and theoretical [12, 13, 14, 15, 16] works examined the conditions for the appearance of MI. In such systems the combination of dissipative nonequilibrium dynamics and the nonlinearity due the interactions results in MI that, after a variable time interval, induces the formation of a train of solitons [3, 4, 5, 17, 18, 19]. On the theoretical side, the far-from-equilibrium dynamics induced by an interaction quench from the repulsive regime to the attractive is usually well captured by mean-field approaches based on the solution of the Gross-Pitaevskii equation with the inclusion of dissipative three-body losses. When the BEC is confined in a quasi one-dimensional waveguide a description based on the non polynomial nonlinear Schrödinger equation typically predicts accurately the number of solitons for relatively weak attractive interactions [20, 5]. The relevant scale that determines the insurgence of MI is the (inverse) most unstable wavenumber $q_{MI}$ . Then, one finds that the number of solitons increases monotonically with the final value of the final scattering length. Interestingly, it is also possible to observe that the particle number of the BEC decreases with a universal power law of the holding time rescaled by the characteristic time for the creation of the modulational instability, independently of the strength of the quench [4]. Recently, the excitation spectrum of matter-wave solitons has also been measured in a quasi one-dimensional cesium condensate [21].

Whereas most works on MI in BECs focused on a single component BEC, the realization of multicomponent systems offers a natural playground to observe nonequilibrium effects in a more general framework. Restricting to two-component BECs, one notices already a rich variety of phases in the ground state. In purely repulsive mixtures, one observes a homogeneous superfluid or a phase separation when inter-species repulsion overcomes the intra-species interaction strength. In the attractive regime a series of recent experiments showed the formation of dilute self-bound droplet states in a two-component BEC both in a tight optical waveguide [22, 23] and in free space [24, 25, 26], closely following the theoretical predictions [27]. In the quasi one-dimensional geometry, upon varying the mean-field interaction from the weakly to the strongly attractive regime, one observes a smooth crossover between bright soliton states and self-bound droplets. Notably, solitons are excitations appearing genuinely in low-dimensional systems. If the one-dimensional interaction strength is attractive (focusing nonlinearity) one retrieves bright solitons, while for repulsive condensates (self-defocusing nonlinearity) one finds dark solitons. Droplets instead result from the competition between mean-field and quantum fluctuation energies with opposite sign.

Although the ground-state properties of binary BEC mixtures have been extensively investigated, the study of the conditions leading to MI in these setups has only recently received attention, both theoretically [28] and experimentally [23]. MI has been observed in the counterflow dynamics of two-component BECs in the miscible (purely repulsive) phase [25]. Also, a recent experiment with coherently coupled BECs rapidly quenched into the attractive regime reported the creation of bright soliton trains formed by dressed-state atoms [29].

In this work we thoroughly investigate the quench dynamics in a binary mixture of Bose-Einstein condensates from the repulsive to the attractive regime in an elongated quasi one-dimensional waveguide. We provide results for quenches from the repulsive to the weakly and strongly attractive regime, where solitonic states and quantum droplets respectively are expected in the ground state. We quantitatively characterize the resulting nonequilibrium dynamics by computing the number of particles and solitons as a function of the holding time after the quench. In section 2 we present a theoretical model for a two-component mixture which allows us to provide a quantitative estimate of the number of solitons following a quench dynamics based on the most unstable mode. In section 3 we discuss our numerical results. We perform simulations for the quench dynamics in three different regimes: repulsive to soliton, repulsive to droplet and soliton to droplet. Finally in section 4 we present our conclusions. In the appendix we provide some details about the ground-state phase diagram, the dynamics of the radial width of the two-component mixture, the effect of beyond-mean-field corrections on the nonequilibrium dynamics and the algorithm used to monitor the number of solitons in our simulations.

2 Model and simulations

We describe the nonequilibrium dissipative dynamics of a binary homonuclear mixture in a quasi one-dimensional waveguide with radial trapping frequency $\omega_{r}$ and components $i,j=1,2$ (with $i\neq j$ ) by the generalized coupled Gross-Pitaevskii equations (GPEs), which, in rescaled units, read

	$\displaystyle i\frac{\partial\psi_{i}}{\partial t}=$	$\displaystyle\left(-\frac{1}{2}\nabla^{2}+V(\mathbf{r})+4\pi\frac{a_{i}}{a_{r}}\left\|\psi_{i}\right\|^{2}+4\pi\frac{a_{ij}}{a_{r}}\left\|\psi_{j}\right\|^{2}+\right.$		(1)
		$\displaystyle\left.\frac{128\sqrt{\pi}}{3}\frac{a_{i}}{a_{r}}\left(\frac{a_{i}}{a_{r}}\left\|\psi_{i}\right\|^{2}+\frac{a_{j}}{a_{r}}\left\|\psi_{j}\right\|^{2}\right)^{\frac{3}{2}}-i\Gamma_{i}\left\|\psi_{i}\right\|^{4}\right)\psi_{i}.$		(1)

Here time, energy and space are scaled in units of the inverse radial trapping frequency $\omega_{r}^{-1}$ , the radial trapping energy $\hbar\omega_{r}$ and the associated harmonic oscillator length $a_{r}=\sqrt{\hbar/m\,\omega_{r}}$ , respectively and sum over the indices in the interaction term is assumed. The coefficients $a_{i}$ and $a_{ij}$ furnish the intra- and inter-species corresponding scattering lengths. The effective mean-field scattering length is thus given by the parameter $\delta a=a_{12}+\sqrt{a_{1}a_{2}}$ . The number of particles in each component equals $N_{1,2}$ . We consider a cigar-shaped trapping harmonic potential $V(\mathbf{r})=m\left[\omega_{r}^{2}(x^{2}+y^{2})+\omega_{z}^{2}z^{2}\right]/2$ with frequencies $\omega_{r}/2\pi=346\,\mathrm{Hz}$ and $\omega_{z}/2\pi=7.6\,\mathrm{Hz}$ , similarly to [4], thus having a radial harmonic oscillator length of $a_{r}=0.87\,\mu$ m. The generalization comes first phenomenologically including a three-body loss term $\Gamma_{i}=L^{(i)}_{3}N_{i}^{2}/(2\omega_{r}a_{r}^{6})$ . Following [24] we do not include mixed two-body or three-body inelastic loss rates. Secondly, in order to take into account corrections due to quantum fluctuations, we introduce the Lee-Huang-Yang (LHY) terms in the GPEs, with component-dependent strength. The ratio between the mean-field chemical potential and the transverse harmonic oscillator energy $(g_{1}n_{1}+g_{2}n_{2})/\hbar\omega_{r}\sim 1$ for typical parameters used in this work, which allows us to use consistently the three-dimensional LHY expression [30]. Note that for purely one-dimensional dynamics an attractive quantum fluctuation term would appear [28].

We perform simulations in two different regimes: (i) symmetric and (ii) realistic. In (i) we set equal initial intraspecies scattering lengths $a_{1i}=a_{2i}=53.0\,a_{0}$ and variable final interactions with the constraint $a_{1f}=a_{2f}$ , where subscripts $i$ and $f$ indicate values before and after quench. We also set $a_{12}=a_{12i}=a_{12f}=-52.0\,a_{0}$ and fix the three-body losses and particle numbers to the same value for each component to $L_{3}=3.6\times 10^{-27}\mathrm{cm^{6}/s}$ . Case (ii) is motivated by current experiments on dilute quantum droplets, with homonuclear mixtures made of two hyperfine states of ${}^{39}K$ , $\left|F=1,m_{f}=0\right>$ and $\left|F=1,m_{f}=-1\right>$ [24, 23]. We consider the interval of magnetic fields $56\,\text{G}<B<57\,\text{G}$ , where for the two hyperfine states the intraspecies scattering lengths $a_{1}$ and $a_{2}$ are both positive (repulsive). The interspecies scattering length $a_{12}$ , instead, is negative. In the setup we are studying the initial intra-component scattering lengths are chosen to be $a_{1i}=75.0\,a_{0}$ and $a_{2}=a_{2i}=a_{2f}=37.5\,a_{0}$ . Component 1 presents a largely varying scattering length $a_{1}$ for a particular experimentally tunable range (via Feshbach resonances), while $a_{2}$ and the inter-species scattering length $a_{12}=a_{12i}=a_{12f}=-52.0\,a_{0}$ remain practically constant for the same range [22, 23, 24]. The three-body loss rates $L_{3}^{(1)}=3.6\times 10^{-27}\mathrm{cm^{6}/s}$ and $L_{3}^{(2)}=6.0\times 10^{-29}\mathrm{cm^{6}/s}$ [24].

Refer to caption — Figure 1: Quench protocol for a two-component BEC mixture from the repulsive to the attractive regime. At time $t=0$ the effective scattering length $\delta a$ is quenched using a linear ramp of the duration of $\delta t=1\,\text{ms}$ from $\delta a=\delta a_{12}+\sqrt{a_{11}a_{22}}>0$ to the attractive regime with $\delta a<0$ either in the soliton phase or to the strongly attractive in the droplet phase. Within the same time interval we switch off the longitudinal confinement. This allows the system to expand freely and to monitor the creation of bright soliton trains propagating along the $z$ -direction. We also consider the quench from the solitonic state to the self-bound droplet which is discussed in the text.

In the absence of a harmonic confining potential and neglecting beyond-mean-field effects one would obtain a dilute gas phase for $\delta a>0$ and a collapsing BEC for $\delta a<0$ . The addition of the beyond-mean-field contribution to the equation of state of the two-component system leads to the stabilization of self-bound quantum droplets due to the competition of attractive mean-field terms proportional to $n^{2}$ and the repulsion due to LHY-type terms proportional to $n^{5/2}$ in the corresponding equation of state [31]. The introduction of the harmonic trap leads to a rich new physical phenomenology. The detailed ground-state phase diagram for the parameters studied in this work is described in A together with the details of the numerical algorithm to obtain the ground states. The algorithm to investigate the soliton dynamics is presented in D.

In most cases, when choosing initial states for our simulations, we start from repulsive inter-species interactions (i.e. $\delta a>0$ ) before quenching to attractive values. The quench is performed by a rapid variation of the scattering length with a linear ramp of $1$ ms. Concurrently, we switch off the longitudinal trapping potential and observe the system expanding in a time-of-flight fashion. See figure 1 for a schematic representation of the quench protocol. We also notice that the small aspect ratio of our trapping potential (i.e. $\omega_{z}/\omega_{r}\ll 1$ ) implies that the most relevant dynamics will happen along the axial direction. In our simulations of the GPE in equation (1) we assume cylindrical symmetry. This choice reduces the computation to an effective two-dimensional calculation of $\psi(r,z)$ , over a grid size of $16\times 8192$ points and a domain of $5\,a_{r}\times 1290\,a_{r}$ allowing to resolve in detail the dynamics along the axial direction. In order to evaluate the derivatives involved in the kinetic term, we employ a discrete (zero-order) Hankel transform in the radial direction and the usual fast Fourier transform in the longitudinal direction [32].

3 Results

3.1 Modulational instability in a binary Bose mixture

In this section we present a theoretical model to describe the MI in a binary BEC after a quench to the attractive regime $\delta a_{f}<0$ . To characterize quantitatively the MI we introduce the Bogoliubov spectrum for a two-component BEC, which reads [33]

E_{\pm}(q)=\sqrt{\frac{\varepsilon_{1}^{2}+\varepsilon_{2}^{2}}{2}\pm\sqrt{\frac{(\varepsilon_{1}^{2}-\varepsilon_{2}^{2})^{2}}{4}+\frac{g_{12}^{2}n_{1}n_{2}q^{4}}{m_{1}m_{2}}}},

(2)

where we defined the Bogoliubov energies of each component

\varepsilon_{i}(q)=\sqrt{\frac{q^{2}}{2m_{i}}\left(\frac{q^{2}}{2m_{i}}+g_{ii}n_{i}\right)}.

(3)

In this work we focus on the equal mass case $m_{1}=m_{2}$ . We define the total density of the system as $n=n_{1}+n_{2}$ . In our simulations we fix the ratio of the densities (and particle numbers) at the initial time to $\frac{n_{1}}{n_{2}}=\sqrt{\frac{a_{2}}{a_{1}}}$ , i.e. the equilibrium configuration that minimizes the energy functional in the attractive regime [34]. Under the condition $\delta a=a_{12}+\sqrt{a_{1}a_{2}}<0$ the lower branch $E_{-}(q)$ becomes unstable. The most unstable mode corresponds to the wavenumber $q_{\text{min}}$ that minimizes the argument of $E_{-}(q)$ .

Upon solving $dE_{-}(q)/dq=0$ we obtain

\displaystyle q_{\text{min}}^{2}

\displaystyle=

\displaystyle 2n\sqrt{a_{1}a_{2}}\left(\sqrt{1-\frac{8\,\delta a}{(\sqrt{a_{1}}+\sqrt{a_{2}})^{2}}+\frac{4\,\delta a^{2}}{\sqrt{a_{1}a_{2}}(\sqrt{a_{1}}+\sqrt{a_{2}})^{2}}}-1\right),

(4)

where we scaled the scattering lengths $a_{i}$ and consequently also $\delta a$ by a factor $2\pi$ to account for the effective soliton dynamics along the longitudinal direction. When $|\delta a|\ll(\sqrt{a_{1}}+\sqrt{a_{2}})^{2}$ then we can expand $q_{\text{min}}$ to a reduced expression to first order in $\delta a$ which reads

\displaystyle q_{\text{min}}^{2}\approx(q_{\text{min}}^{\text{red}})^{2}

\displaystyle=

\displaystyle 8n\,|\delta a|\frac{\sqrt{a_{1}a_{2}}}{(\sqrt{a_{1}}+\sqrt{a_{2}})^{2}}.

(5)

We can now define the associated wavelength $\lambda_{0}=2\pi/q_{\text{min}}$ . Starting with a quasi-1D binary Bose mixture of length $L$ (ex. the longitudinal Thomas-Fermi radius), a sudden quench to the attractive regime induces the formation of a train of solitons. The number of solitons $N_{s}$ can then be estimated as

N_{s}=\frac{L}{\lambda_{0}}.

(6)

3.2 Number of solitons

We now discuss the creation of soliton trains induced by the MI. In figure 2 we show snapshots of the density for the symmetric regime (i) at two different times: one right before the quench from the repulsive mean-field regime to the attractive one and another at $t=30$ ms. Figure 2(a) corresponds to a quench to the weakly attractive regime $\delta a=-3\,a_{0}$ where the ground state of the system is an extended soliton (see A for a thorough discussion of the phase diagram). Figure 2(b) corresponds to a quench to the strongly attractive regime $\delta a=-9\,a_{0}$ where the ground state is instead a self-bound droplet. For both cases the initial state is the ground state of the system for repulsive interactions. Notice that for the sake of visualization the density is normalized at its peak value, therefore the color code is the same for both components to emphasize the density modulations, even if the number of particles changes with time. We observe that, as soon as MI sets in, density peaks are formed, creating a soliton train. The number of solitons increases with the strength of the attractive interaction.

We compute the number of solitons numerically from the peaks of the density distribution for each of the two components, employing an algorithm that we describe in detail in D. The results of this analysis are shown in figure 3.

In figure 3 the average number of solitons computed from the numerics (points) is compared to our prediction from equations (4) and (6) (dashed black line). We notice that the expression for $q_{\text{min}}^{\mathrm{red}}$ in equation (5) reproduces to an excellent approximation $q_{\text{min}}$ for the parameters used in this work. For weakly attractive interactions the agreement between the numerics and the analytical result is good for $|\delta a|\lesssim 6\,a_{0}$ . For stronger attractive interactions we observe a larger deviation of the analytical prediction from the numerics, likely due to the far from equilibrium dynamics involved in the creation of the soliton train.

In figure 2(c) we show the snapshot of the density starting from a solitonic, weakly attractive configuration to the regime of strong attraction. We observe that the initial condensate splits into just two bright solitons. This has to be compared to figure 2(b), where, due to the larger initial longitudinal length, the quench dynamics produces a soliton train with several density peaks.

We also performed simulations in the realistic case (ii) for the experimental parameters of section 2. The dynamics is significantly more complex than in the symmetric case. First the asymmetry in the number of particles (see section 3.3 and inset of figure 4) is such that the particles in the second component is almost constant during the expansion dynamics after the quench, whereas $N_{1}(t)$ is greatly reduced after tens of milliseconds, similarly to the symmetric case. The effect is that the second component is only weakly affected by the attractive dynamics due to the limited overlap with the first component. Therefore the soliton trains observed after the quench are poorly described by the theory described in section 3. We provide further details in D.

3.3 Number of atom loss

In this section we discuss the evolution of the number of particles as a function of time after the quench. In figure 4 we show the results for both quantities from GPE simulations of the two coupled components $N_{1}(t)$ and $N_{2}(t)$ for (a) the case of symmetric interactions and losses (b) and for the realistic experimental parameters of section 2. In the symmetric case $N_{1}(t)=N_{2}(t)$ for all times, whereas in the realistic case the number of particles in components $1$ and $2$ decrease with time at different rates because of the different three-body losses coefficients $L_{3}^{(1)},L_{3}^{(2)}$ . The first component has a much larger three-body decay, resulting in a more complex dynamics. We observe that for all the final attractive mean-field interactions considered in figure 4, the number of particles for short times slowly decreases before establishing the MI at $t\approx t_{MI}$

\displaystyle t_{MI}\,\omega_{r}

\displaystyle=

\displaystyle\frac{(\sqrt{a_{1}}+\sqrt{a_{2}})^{2}}{\sqrt{a_{1}a_{2}}}\frac{1}{4n_{0}|\delta a|a_{r}^{2}}.

(7)

where $n_{0}$ is the peak density of the initial configuration [20]. Consistently with the prediction of equation (7), in our simulations for the symmetric (a) and the realistic case (b), quenching to the strongly attractive regime leads to faster decrease of the number of particles.

4 Conclusions

In this work we studied the nonequilibrium dynamics of a two-component Bose-Einstein condensate after a quantum quench to the attractive interspecies interactions. We specialized to the experimentally relevant case of potassium binary mixture. Quenching the effective mean-field scattering length from repulsive to attractive values in a wide interval we observed a MI and the creation of soliton trains. We characterized quantitatively the number of solitons via numerical simulations of the coupled Gross-Pitaevskii equations. In the stationary, long-time limit we observed that an analytical model based on the calculation of the most unstable Bogoliubov mode is in reasonable quantitative agreement with the number of solitons for both components in the symmetric configuration. The experimentally relevant case, with asymmetric intraspecies interactions and different loss rates, leads to a more intricate dynamics which is only qualitatively captured by our model. The related time scale for the rise of the instability however does not translate into a universal scaling for the losses of both components, in contrast to what was recently observed for a single component lithium BEC with small final scattering length $|a_{f}|$ [4].

We emphasize that this work focuses on a far-from-equilibrium dynamical regime. For the fast magnetic field ramps considered here, even in the strongly attractive regime $|\delta a_{f}|\gtrsim 7\,a_{0}$ , the solitonic bumps in the density are not self-bound droplets, as their width equals the transverse harmonic oscillator length and the atom number decay is different from what has been observed in the formation of self-bound droplets [34] (see also B).

The MI analysis can be used to study also other Bose-condensed systems. A sudden quench of the s-wave scattering length can be applied not only to atomic gases in the same or different hyperfine states, but also to heteronuclear bosonic mixtures. Moreover, in bosonic systems with spin-orbit and Rabi couplings the MI can be induced by varying these one-body couplings. However, our work strongly suggests that generically one cannot trust only the analytical calculations based on the most unstable mode of the elementary excitations: a comparison with numerical simulation is needed to obtain reliable predictions. Extensions of this work may include a systematic study of the effects of the coherent coupling of a two-components mixture [35, 29, 36], the inclusion of long-range dipolar interactions [37], or the investigation of finite-temperature effects [38] across the normal-to-BEC transition in the attractive regime and its connection to the Kibble-Zurek mechanism [39, 40].

We thank D. Luo, R. Hulet, and S. Wüster for useful discussions. This research was developed with the help of XMDS2 software [41]. We thank the High Performance Computing Center (NPAD) at UFRN for providing computational resources. T.M. acknowledges CNPq for support through Bolsa de produtividade em Pesquisa n.311079/2015-6. This work was supported by the Serrapilheira Institute (grant number Serra-1812-27802), CAPES-NUFFIC project number 88887.156521/2017-00. T.M. thanks the Physics Department of the University of L’Aquila for the hospitality where part of the work was done. A.C. is supported by FAPESP through Grant No. 2017/09390-7. LS acknowledges the BIRD project “Superfluid properties of Fermi gases in optical potentials” of the University of Padova for financial support.

Appendix A Ground-state phase diagram

In this appendix we discuss the ground-state properties of a two-component mixture in the attractive regime in a cigar-shaped harmonic potential. In figure 5 we show the numerical and the variational phase diagram obtained by imaginary-time evolution of the generalized GPE of equation (1) (numerical) and by minimizing the corresponding energy functional (variational). For the variational approach we employ a gaussian wavefunction

\psi(\textbf{r})=\sqrt{\frac{N}{\pi^{\frac{3}{2}}\sigma_{r}^{2}\sigma_{z}}}\exp\left(-\sum_{r_{i}=x,y,z}\frac{r_{i}^{2}}{2\sigma_{r_{i}}^{2}}\right)

(8)

with variational parameters $\sigma_{r}$ and $\sigma_{z}$ . The wavefunction is normalized to the total number of particles $||\psi||^{2}=N$ [31].

We observe a smooth crossover from the droplet to the soliton phase. First, for low particle number or, equivalently, small values of $|\delta a|$ , the ground state of the system corresponds to a soliton, whose shape depends on the external trapping, for which $\sigma_{r}\sim a_{r}$ , while $\sigma_{z}\gg a_{r}$ . Specifically, beyond-mean-field corrections are not necessary for the stability of this state [20]. Reducing $\delta a$ the ground state is (almost) isotropic $\sigma_{r}\sim\sigma_{z}<a_{r}$ , independent of the confinement aspect ratios. The existence of this self-bound state is enabled by taking into account the contribution of Gaussian quantum fluctuations in the variational energy equation (1).

We compare the ground states obtained from the variational analysis with the numerical simulation of the GP equation. After relaxation in imaginary time, we find that the ground state of the system for the strongly attractive inter-species scenario is a self-trapped, spherical droplet state.

Appendix B Evolution of the radial width

The results of the dynamics of the radial widths of each of the two components of the mixture after the quench are shown in figure 6. The radii of both components fluctuate closely to the initial radial harmonic oscillator length $a_{r}$ for both quenches to the weakly attractive and strongly attractive regimes. In particular, for the cases where $\delta a_{f}=-9.0\,a_{0}$ (symmetric) $\delta a_{f}=-10.92\,a_{0}$ (realistic) and which have a self-bound droplet phase as a ground state (see figure 5) we observe no signature of such a state throughout the dynamics.

Appendix C Effect of quantum fluctuations

The effect of the quantum fluctuations for two-component Bose gas is described by the term $g_{i}^{\mathrm{LHY}}=\frac{128\sqrt{\pi}}{3}\frac{a_{i}}{a_{r}}\left(\frac{a_{i}}{a_{r}}\left|\psi_{i}\right|^{2}+\frac{a_{j}}{a_{r}}\left|\psi_{j}\right|^{2}\right)^{\frac{3}{2}}$ in the extended GPE equation (1). We tested the effect of the removal of this term in our simulations for quenches into the weak and strong attractive regime. The results are shown in figure 7 for symmetric interactions and loss terms and for the realistic experimental configuration. Whereas the results for the second component in the realistic case are almost identical, the first component displays a much faster particle number decay in the absence of beyond-mean-field terms. This can be explained by the attractive nature of the mean-field interaction which is not balanced by the additional repulsive LHY correction leading to higher densities and therefore to higher losses. At the same time, MI takes place on a shorter time scale, leading to a faster stabilization of the particle number. Therefore, at longer times, in the absence of beyond-mean-field effects we observe a larger particle number.

Appendix D Details on the soliton number algorithm and comparison with theory

We briefly describe the algorithm to compute the number of solitons from the dynamics of the density $n(z,t)$ integrated along the transverse directions as a function of time. Our method identifies local maxima in $n(z,t)$ as bright solitons. We start by disregarding, at each instant, peaks in low-density regions (with amplitudes $\lesssim 10\%$ of the mean initial density). A weak Gaussian filter is applied to smooth out most numerical effects at distances much smaller than the typical healing length $\xi$ of the system. Finally, we avoid overcounting solitons which are undergoing a probable splitting process by treating as one visible peaks that are apart by distances smaller than $\xi$ at a particular instant.

In figure 8 we provide further numerical results about the comparison between the realistic and the symmetric case with our theoretical model. In the realistic case, as discussed in the main text, we only plot the number of solitons in the first component $N_{s}^{1}$ . The dynamics in the realistic case is significantly more complex than in the symmetric one. Therefore the soliton trains observed after the quench are inadequately described by the theory in this context. Due to this poor agreement between the numerics and the estimate from equation (6), we only show the effect of three-body losses on the solitons number in the more controllable symmetric case. We observe that in the absence of losses equation (6) is still able to provide a reasonable estimate of the number of solitons for quenches to intermediate values of $|\delta a_{f}|$ . However, for larger $|\delta a_{f}|$ the deviation can be as much as $100\%$ . In the realistic case it can be argued that a similar increase in the number of solitons might be observed in the absence of dissipation, therefore partially compensating for the mismatch of the theoretical (dashed line) and the numerics. However a more refined theory is needed to explain the behavior of both components.

References

[1] T. Brooke Benjamin and J. E. Feir. The disintegration of wave trains on deep water part 1. theory. Journal of Fluid Mechanics, 27(3):417–430, 1967.
[2] H C Yuen and B M Lake. Instabilities of waves on deep water. Annual Review of Fluid Mechanics, 12(1):303–334, 1980.
[3] Kevin E. Strecker, Guthrie B. Partridge, Andrew G. Truscott, and Randall G. Hulet. Formation and propagation of matter-wave soliton trains. Nature, 417(6885):150–153, 2002.
[4] Jason H. V. Nguyen, De Luo, and Randall G. Hulet. Formation of matter-wave soliton trains by modulational instability. Science, 356(6336):422–426, April 2017.
[5] 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. Observation of a modulational instability in bose-einstein condensates. Phys. Rev. A, 96:041601, Oct 2017.
[6] Tadej Mežnaršič, Tina Arh, Jure Brence, Jaka Pišljar, Katja Gosar, Žiga Gosar, Rok Žitko, Erik Zupanič, and Peter Jeglič. Cesium bright matter-wave solitons and soliton trains. Phys. Rev. A, 99:033625, Mar 2019.
[7] L. D. Carr and Y. Castin. Dynamics of a matter-wave bright soliton in an expulsive potential. Phys. Rev. A, 66:063602, Dec 2002.
[8] K. Tai, A. Hasegawa, and A. Tomita. Observation of modulational instability in optical fibers. Phys. Rev. Lett., 56:135–138, Jan 1986.
[9] C. Fort, F. S. Cataliotti, L. Fallani, F. Ferlaino, P. Maddaloni, and M. Inguscio. Collective excitations of a trapped bose-einstein condensate in the presence of a 1d optical lattice. Phys. Rev. Lett., 90:140405, Apr 2003.
[10] Govind Agrawal. Nonlinear Fiber Optics. Academic Press, 2019.
[11] S.G. Thornhill and D. [ter Haar]. Langmuir turbulence and modulational instability. Physics Reports, 43(2):43 – 99, 1978.
[12] L. Salasnich, A. Parola, and L. Reatto. Modulational instability and complex dynamics of confined matter-wave solitons. Phys. Rev. Lett., 91:080405, Aug 2003.
[13] A. Smerzi, A. Trombettoni, P. G. Kevrekidis, and A. R. Bishop. Dynamical superfluid-insulator transition in a chain of weakly coupled bose-einstein condensates. Phys. Rev. Lett., 89:170402, Oct 2002.
[14] L. D. Carr and J. Brand. Pulsed atomic soliton laser. Phys. Rev. A, 70:033607, Sep 2004.
[15] T. Karpiuk, M. Brewczyk, S. Ospelkaus-Schwarzer, K. Bongs, M. Gajda, and K. Rza¸żewski. Soliton trains in bose-fermi mixtures. Phys. Rev. Lett., 93:100401, Sep 2004.
[16] Giacomo Gori, Tommaso Macrì, and Andrea Trombettoni. Modulational instabilities in lattices with power-law hoppings and interactions. Phys. Rev. E, 87:032905, Mar 2013.
[17] U. Al Khawaja, H. T. C. Stoof, R. G. Hulet, K. E. Strecker, and G. B. Partridge. Bright soliton trains of trapped bose-einstein condensates. Phys. Rev. Lett., 89:200404, Oct 2002.
[18] H. Kiehn, S. I. Mistakidis, G. C. Katsimiga, and P. Schmelcher. Spontaneous generation of dark-bright and dark-antidark solitons upon quenching a particle-imbalanced bosonic mixture. Physical Review A, 100(2), August 2019.
[19] S I Mistakidis, G C Katsimiga, P G Kevrekidis, and P Schmelcher. Correlation effects in the quench-induced phase separation dynamics of a two species ultracold quantum gas. New Journal of Physics, 20(4):043052, April 2018.
[20] L. Salasnich, A. Parola, and L. Reatto. Condensate bright solitons under transverse confinement. Phys. Rev. A, 66:043603, Oct 2002.
[21] Andrea Di Carli, Craig D. Colquhoun, Grant Henderson, Stuart Flannigan, Gian-Luca Oppo, Andrew J. Daley, Stefan Kuhr, and Elmar Haller. Excitation modes of bright matter-wave solitons. Physical Review Letters, 123(12), Sep 2019.
[22] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell. Quantum liquid droplets in a mixture of bose-einstein condensates. Science, 359(6373):301–304, 2018.
[23] P Cheiney, CR Cabrera, J Sanz, B Naylor, L Tanzi, and L Tarruell. Bright soliton to quantum droplet transition in a mixture of bose-einstein condensates. Physical review letters, 120(13):135301, 2018.
[24] G Semeghini, G Ferioli, L Masi, C Mazzinghi, L Wolswijk, F Minardi, M Modugno, G Modugno, M Inguscio, and M Fattori. Self-bound quantum droplets of atomic mixtures in free space. Physical review letters, 120(23):235301, 2018.
[25] C. Hamner, J. J. Chang, P. Engels, and M. A. Hoefer. Generation of dark-bright soliton trains in superfluid-superfluid counterflow. Phys. Rev. Lett., 106:065302, Feb 2011.
[26] C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich, F. Ancilotto, M. Modugno, F. Minardi, and C. Fort. Observation of quantum droplets in a heteronuclear bosonic mixture. Phys. Rev. Research, 1:033155, Dec 2019.
[27] DS Petrov. Quantum mechanical stabilization of a collapsing bose-bose mixture. Physical review letters, 115(15):155302, 2015.
[28] Thudiyangal Mithun, Aleksandra Maluckov, Kenichi Kasamatsu, Boris A. Malomed, and Avinash Khare. Modulational instability, inter-component asymmetry, and formation of quantum droplets in one-dimensional binary bose gases. Symmetry, 12(1), 2020.
[29] J. Sanz, A. Frölian, C. S. Chisholm, C. R. Cabrera, and L. Tarruell. Interaction control and bright solitons in coherently-coupled bose-einstein condensates, 2019.
[30] Paweł Zin, Maciej Pylak, Tomasz Wasak, Mariusz Gajda, and Zbigniew Idziaszek. Quantum bose-bose droplets at a dimensional crossover. Phys. Rev. A, 98:051603, Nov 2018.
[31] Alberto Cappellaro, Tommaso Macrì, and Luca Salasnich. Collective modes across the soliton-droplet crossover in binary bose mixtures. Phys. Rev. A, 97:053623, May 2018.
[32] The zero-order Hankel transform of a function $f(r,z)$ is given by $\tilde{f}(k_{r},k_{\phi},z)=-2i\pi\int_{0}^{\infty}f(r,z)J_{0}(k_{r}r)\,r\ dr$ . The full momentum-space representation of $f(r,z)$ is then $\tilde{f}({\bf k})=\mathscr{F}_{z}[\tilde{f}(k_{r},k_{\phi},z)]$ , with $\mathscr{F}_{z}[\cdot]$ the usual Fourier transform along the $z-$ direction.
[33] David M Larsen. Binary mixtures of dilute bose gases with repulsive interactions at low temperature. Annals of Physics, 24:89 – 101, 1963.
[34] G. Ferioli, G. Semeghini, S. Terradas-Briansó, L. Masi, M. Fattori, and M. Modugno. Dynamical formation of quantum droplets in a ${}^{39}\mathrm{K}$ mixture. Phys. Rev. Research, 2:013269, Mar 2020.
[35] Alberto Cappellaro, Tommaso Macrì, Giovanni F. Bertacco, and Luca Salasnich. Equation of state and self-bound droplet in rabi-coupled bose mixtures. Scientific Reports, 7(1), oct 2017.
[36] Ishfaq Ahmad Bhat, T. Mithun, B. A. Malomed, and K. Porsezian. Modulational instability in binary spin-orbit-coupled bose-einstein condensates. Phys. Rev. A, 92:063606, Dec 2015.
[37] Igor Ferrier-Barbut, Matthias Wenzel, Matthias Schmitt, Fabian Böttcher, and Tilman Pfau. Onset of a modulational instability in trapped dipolar bose-einstein condensates. Phys. Rev. A, 97:011604, Jan 2018.
[38] K L Lee, N B Jørgensen, L J Wacker, M G Skou, K T Skalmstang, J J Arlt, and N P Proukakis. Time-of-flight expansion of binary bose–einstein condensates at finite temperature. New Journal of Physics, 20(5):053004, may 2018.
[39] P. Comaron, F. Larcher, F. Dalfovo, and N. P. Proukakis. Quench dynamics of an ultracold two-dimensional bose gas. Phys. Rev. A, 100:033618, Sep 2019.
[40] Yiping Chen, Munekazu Horikoshi, Kosuke Yoshioka, and Makoto Kuwata-Gonokami. Dynamical critical behavior of an attractive bose-einstein condensate phase transition. Phys. Rev. Lett., 122:040406, Feb 2019.
[41] Graham R. Dennis, Joseph J. Hope, and Mattias T. Johnsson. XMDS2: Fast, scalable simulation of coupled stochastic partial differential equations. Computer Physics Communications, 184(1):201–208, January 2013.