Soliton trains after interaction quenches in Bose mixtures

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.


Introduction
Modulational instability 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 [4,17,3,5]. On the theoretical side, the far-fromequilibrium 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 [18,5]. The relevant scale that determines the insurgence of modulational instability is the (inverse) most unstable wavenumber q M I . 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 [19].
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 [20] and in free space [21,22,23], closely following the theoretical predictions [24]. In the quasi onedimensional 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 [25] and experimentally. MI has been observed in the counterflow dynamics of two-component BECs in the miscible (purely repulsive) phase [22]. 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 [26].
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. We then continue presenting 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.

Model and simulations
We describe the nonequilibrium dissipative dynamics of a binary homonuclear mixture in a quasi one-dimensional waveguide with radial trapping frequency ω r and components i, j = 1, 2 (with i = j) by the generalized coupled Gross-Pitaevskii equations (GPEs), which, in rescaled units, read Here time, energy and space are scaled in units of ω −1 r , ω r and a 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 δa = a 12 + √ a 1 a 2 . The number of particles in each component equals N 1,2 . We consider a cigar-shaped trapping harmonic potential V (r) = m [ω 2 r (x 2 + y 2 ) + ω 2 z z 2 ] /2 with frequencies ω r /2π = 346 Hz and ω z /2π = 7.6 Hz, similarly to [4], having thus a radial harmonic oscillator length of a r = /m ω r = 0.87 µm. The generalization comes first phenomenologically including a three-body loss term Γ i = L (i) 3 N 2 i /(2ω r a 6 r ). 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 [21]. The ratio between the mean-field chemical potential and the transverse harmonic oscillator energy (g 1 n 1 + g 2 n 2 )/ ω r ∼ 1 for typical parameters used in this work, which allows us to use consistently the threedimensional LHY expression [27]. Note that for purely one-dimensional dynamics an attractive quantum fluctuation term would appear [25].
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 fix the three-body losses and particle numbers to the same value for each component. Case (ii) is motivated by current experiments on dilute quantum droplets, with homonuclear mixtures made of two hyperfine states of 39 K, |F = 1, m f = 0 and |F = 1, m f = −1 [21,20]. We consider the interval of magnetic fields 56 G < B < 57 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 [28,21]. The three-body loss rates L In the absence of a harmonic confining potential and neglecting beyond-mean-field effects one would obtain a dilute gas phase for δa > 0 and a collapsing BEC for δa < 0. The addition of the beyond-mean-field contribution to the equation of state of the twocomponent 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 [29]. 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 Appendix A.
In most cases, when choosing initial states for our simulations, we start from repulsive inter-species interactions (i.e. δa > 0) before quenching to attractive values. The quench is performed by a rapid variation of the scattering length with a linear ramp of 1ms. 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. ω z /ω r 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 twodimensional calculation of ψ(r, z), over a grid size of 16 × 8192 points and a domain of  5 a r × 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 (zeroorder) Hankel transform in the radial direction and the usual fast Fourier transform in the longitudinal direction [30].

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 δa f < 0. To characterize quantitatively the MI we introduce the Bogoliubov spectrum for a two-component BEC, which reads [31] where we defined the Bogoliubov energies of each component Figure 2. Snapshots of the time-evolution of the density at the initial time t = 0 and at t = 30ms after a quench with a 1ms-ramp for different initial and final scattering lengths. We consider the initial number of particles N 1 = N 2 = 2.5×10 4 and symmetric intraspecies interactions a 1 = a 2 which are quenched simultaneously. The interspecies scattering length a 12 is left unaltered during the quench. The system is initialized in the ground state for the corresponding value of δa i . The three-body loss parameters are chosen equal to L 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 n 1 n 2 = a 2 a 1 , i.e. the equilibrium configuration that minimizes the energy functional in the attractive regime [32]. Under the condition δa = a 12 + √ a 1 a 2 < 0 the lower branch E − (q) becomes unstable. The most unstable mode corresponds to the wavenumber q min that minimizes the argument of E − (q).
Upon solving dE − (q)/dq = 0 we obtain where we scaled the scattering lengths a i and consequently also δa by a factor 2π to account for the effective soliton dynamics along the longitudinal direction. When |δa| ( √ a 1 + √ a 2 ) 2 then we can expand q min to a reduced expression to first order in δa which reads We notice that the expression for q red min in equation (5) reproduces to an excellent approximation q min for the parameters used in this work.
We can now define the associated wavelength λ 0 = 2π/q 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

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 = 30ms. Figure 2(a) corresponds to a quench to the weakly attractive regime δa = −3 a 0 where the ground state of the system is an extended soliton (see Appendix A for a thorough discussion of the phase diagram). Figure 2(b) corresponds to a quench to the strongly attractive regime δa = −9 a 0 where the ground state is instead a selfbound 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 for t > t M I modulational instability sets in, creating peaks in the density in the form of 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 Appendix 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). For weakly attractive interactions the agreement between the numerics and the analytical result is good for |δa| 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 and therefore not shown.

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 . 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 modulational instability at t ≈ t M I where n 0 is the peak density of the initial configuration [18]. 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. We observe that, contrarily to what was verified for a single component BEC, the dynamics of the number of particles for quenches into the strongly attractive regime does not collapse to a universal curve [4].
3 , therefore the second component (inset) practically maintains its particle number throughout the evolution. In (a) the dissipation acts on the two components equally, therefore N 1 (t) = N 2 (t) for all times. Different curves in (a) and (b) correspond to different final effective interaction strength as in the color scale bar.

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 modulational instability 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 [4].
We emphasize that this work focuses on a far-from-equilibrium dynamical regime. Therefore the connection to the physics of the ground state which displays (extended) solitons and droplet state, as described in the appendix, is not evident. For example, upon quenching to the strongly attractive regime |δa f | 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 [32].
Extensions of this work may include a systematic study of the effects of the coherent coupling of a two-components mixture [33,34,35], the inclusion of long-range dipolar interactions [36], or the investigation of finite-temperature effects [37] across the normalto-BEC transition in the attractive regime and its connection to the Kibble-Zurek mechanism [38,39]. For large negative δa the system is in the droplet phase σ r ∼ σ z a r . For small negative δa the system is in the soliton phase σ r ∼ a r < σ z . Numerical simulation of the ground state wavefunction obtained by imaginary-time evolution of the generalized GPE of equation (1) in the absence of dissipation are shown with orange (σ r ) and blue (σ z ) points.
In this apendix we discuss the ground-state properties of a two-component mixture in the attractive regime in a cigar-shaped harmonic potential. In figure A1 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 with variational parameters σ r and σ z . The wavefunction is normalized to the total number of particles ||ψ|| 2 = N [29]. We observe a smooth crossover from the droplet to the soliton phase. First, for low particle number or, equivalently, small values of |δa|, the ground state of the system corresponds to a soliton, whose shape depends on the external trapping, for which σ r ∼ a r , while σ z a r . Specifically, beyond-mean-field corrections are not necessary for the stability of this state [18]. Reducing δa the ground state is (almost) isotropic σ r ∼ σ z < a r , independent of the confinement aspect ratios. The existence of this selfbound 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 B1. 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 δa f = −9.0 a 0 (symmetric) δa f = −10.92 a 0 (realistic) and which have a self-bound droplet phase as a ground state (see figure A1) we observe no signature of such a state throughout the dynamics.  Figure C1. Effect of beyond-mean-field corrections on the dynamics of the particle number decay after the quench for (a) the symmetric case with equal interactions a 1 = a 2 and loss rates L  (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 C1 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, modulational instability 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. Soliton number algorithm
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 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 ξ 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 ξ at a particular instant.