Propagation properties and stability of dark solitons in weakly interacting Bose-Bose droplets

We investigate dark solitons in two-component Bose systems with competing interactions in one dimension. Such a system hosts a liquid phase stabilized by the beyond-mean field corrections. Using the generalized Gross-Pitaevskii equation, we reveal the presence of two families of solitonic solutions. The solitons in both of them can be engineered to be arbitrarily wide. One family of solutions, however, has got an anomalous dispersion relation and our analyses show one of its branches is unstable. We find the presence of a critical velocity demarcating the stable from unstable solutions. Nonetheless, grey anomalous solitons are able to exist inside quantum droplets and can be treated as solitonic excitations thereof.


I. INTRODUCTION
Ultracold atomic systems with competing interactions have been a subject of an extensive theoretical research, powered up by various experiments.The addition of corrections for quantum fluctuations [1] to the mean-field theory led to the correct description of the emergence of self-bound objects called quantum droplets in two-component systems [2][3][4][5][6].A quite unexpected appearance of droplets in dipolar gases [7][8][9] also became explicable by this effect.All in all, the mean-field approach tumbled down in these two regimes and the use of its version with LHY corrections [10,11] -known as the generalized Gross-Pitaevskii equation (GGP), became crucial.
The GGP [12] proves its usefulness in the weakly-interacting regime [13] of Bose-Bose mixtures.It has been employed to study the properties of quantum droplets [14], including elementary excitations [15] like breathing modes [16], even in the case of unequal number of bosons of each species [17,18].
There is a variety of already known dark-soliton-like solutions of the GGP in the miscible regime, such as kink-type solitons [24], dark quantum droplets [25] and standard dark solitons [26].
Let us briefly describe the most important properties of these excitations.The kink-type soliton density profile has got two different asymptotic values with a rapidly varying region located at its origin.In this case, the phase of the wave function remains constant.
Dark quantum droplets, on the other hand, have a shape of an inverted quantum droplet.Their phase pattern is non-trivial, with −π/2 phase on one side of the depletion, π/2 on the other one with an intermediate step in the wide depleted density region of phase equal to 0. All three regions are smoothly connected.
Standard dark solitons can be divided into fully and partially depleted -black and grey correspondingly.Black solitons have a π-jump in their phase.The phase of grey solitons changes smoothly and the total phase difference is smaller than π.All in all, they share the features of dark solitons in single-component Bose gases [27,28].
Yet, there is still another type of solitons, so far found only in a beyond-LHY description called LLGPE of 1D Bose-Bose mixtures [29] and dipolar Bose gases [30], but had been retrieved as the solutions of some differential equations [31] and called "solitonlike bubbles" there.These anomalous solitons can be arbitrarily wide, are never fully depleted and have a constant phase profile.Moreover, they have a peculiar dispersion relation with an additional subbranch and a cusp.As they were not found in [25], one may think they appear due to beyond-LHY contributions to the method.
In this article, we derive the solitonic solution of the GGP and reveal the presence of anomalous solitons.We also check the stability of the solutions and show that dark solitons can exist inside quantum droplets (see figure 1).

II. FRAMEWORK
We look into a weakly interacting two-component gas of Ñ bosons in 1D.The intraspecies interaction is repulsive, whereas the intercomponent ones are attractive.Their masses are the same, i.e. m ↑ = m ↓ = m with σ = {↑, ↓} denoting the component.The energy density functional for this system is given by [12]: where ρσ is the σ-component density and g σσ ′ for σ = σ ′ is the intracomponent interaction strength and the intercomponent interaction strength otherwise, i.e. when σ ̸ = σ ′ .In the miscible system without any external confinement the densities are tied up via the condition ρ↑ /ρ ↓ = g ↓↓ /g ↑↑ [1].This leads to a simplification of the energy density functional E int .An approach based on the local density approximation enables us to use E int and write the generalized Gross-Pitaevskii equation for the wave function Φ(x, t) as follows [12]: where Φ(x, t) is related to single-component wave functions Φσ via Φσ (x, t) = g 1/4 σ σ Φ(x, t) √ g ↑↑ + √ g ↓↓ , where ↑ =↓ and ↓ =↑.Moreover, we have defined δg = g ↑↓ + √ g ↑↑ g ↓↓ , which is assumed to be positive throughout the work.
We now introduce the units of length , energy E 0 = ℏ/t 0 and normalization factor of the wave function Φ 0 = [15].With t = t 0 t, x = x 0 x, Φ = Φ 0 Φ and Ẽ = E 0 E, we can rewrite (2) in the dimensionless form [14]: The normalization condition of the wave function Φ is related to the real number of atom in the system Ñ via , assuming a box of size L with periodic boundary conditions.

A. Speed of sound and equilibrium density
The stable and unstable liquid phases in Bose-Bose mixtures are demarcated by the point where the speed of sound becomes imaginary [13].We use the following time-dependent Ansatz to linearize (3): where Φ 0 = N/L, µ 0 = N L − N L is the chemical potential corresponding to the constant density profile and δΦ is a small perturbation.Then, with ω ρ = ρ 0 − 1 2 √ ρ 0 , we obtain a set of Bogoliubov-de Gennes (BdG) equations Assuming that δΦ = ue (ikx−iωt) + v * e −(ikx−iωt) [32], we diagonalize the Hamiltonian matrix an obtain the following eigenfrequencies The spectrum is linear in the limit of low momenta, well approximated by the phonon dispersion law ω = ck, where is the speed of sound.On the other hand, when k ≫ 1, ( 6) reduces to the free particle spectrum ω = k 2 /2.From ( 7), we see the speed of sound is imaginary when ρ 0 < 1/4.It indicates that the stable-to-unstable liquid transition happens when ρ 0 = 1/4 ≡ ρ ins and this is a phonon instability.
The liquid phase is characterized by the presence of a minimum in the energy per particle function [13,14] This minimum occurs when ρ 0 = 4/9 and we will further refer to this value as the equilibrium density ρ eq .When N/L < ρ eq , the system prefers to form a quantum droplet with a bulk density equal to ρ eq .This interpretation is valid if N ≫ 1.In such a case, the surface energy needed to form the droplet density profile is small in comparison to the bulk energy and therefore negligible.

B. Solitonic solutions
We now look for a solution of a dispersionless wave moving with velocity v through an infinite constant background density.Such an object is represented by the following wave function: where and ζ is the comoving coordinate ζ = x − vt.
Figure 2 shows the density minimum min ρ(x, t = 0) and root mean square (RMS) width ⟨x 2 ⟩ (assuming ⟨x⟩ = 0) of the solitonic excitations as functions of the background density.One can distinguish 3 regions there: (i) unstable liquid, (ii) anomalous, and (iii) the standard one.
In the unstable region, when ρ ∞ < ρ ins , even a small perturbation to a uniform density can cause a violent dynamics in the system.Thus, it is no surprise that there are no solitonic solutions there.
Otherwise, when ρ ins < ρ ∞ < ρ eq , we are in the anomalous regime.The inequality fulfils the necessary condition for the existence of an anomalous soliton [31,33], namely 0 < ρ min < ρ ∞ < ρ eq such that U (ρ min ) = U (ρ ∞ ) = 0 and U (ρ) < 0 ∀ ρ min ⩽ ρ ⩽ ρ ∞ as well as the stability of a uniform system, i.e. c 2 > 0. A motionless anomalous soliton is shown in the inset of figure 2(a).Its characteristic feature is a partial depletion of density even if the soliton is motionless (β ≡ v/c = 0).Another important property of the motionless soliton in the anomalous regime is a constant phase profile ϕ(ζ; β = 0) = const.b,d) Soliton stability criterion.Solitons are stable when dP/dβ < 0. In the standard regime, the solitons are stable everywhere.In the anomalous one -only above the critical velocity, which in this case of ρ∞ = 0.36 is numerically evaluated to be βcr = 0.2550 (29).[Red circles in panels (c,d) mark the place, where the velocity of solitons equals to the critical velocity βcr.] In the last regime, ρ ∞ > ρ eq , the solitons have standard properties in terms of their density minimum and phase profile.Namely the density of a black (β = 0) soliton reaches zero and there is a typical π-jump in the phase.
As one can see in figure 2(b), the motionless solitons -both anomalous and standard ones -become ultrawide when ρ ∞ → ρ eq .It gives an opportunity to steer their size with the number of atoms in the system and interaction strengths g σσ ′ .The presence of arbitrarily wide solitons with a substantial density depletion in weakly interacting Bose-Bose mixtures is an opportunity for taking an in situ absorption image of a dark soliton in such systems.
Another distinct property of the anomalous solitons, already described in [30], is the appearance of a subbranch in the dispersion relation.In figure 3(a), (c), we present the relation of the renormalized2 energy [25]: and regularized momentum [34]: Also in this case, there is a violent change in the energy spectrum when we cross ρ ∞ = 4/9 from above and enter the anomalous region.Namely, an additional subbranch altogether with a cusp appears.
It makes the effective mass m eff = (d 2 E/dP 2 ) −1 not properly defined due to the lack of derivative.The stability analysis reveals another meaning of the two subbranches.The soliton is stable when dP/dβ < 0 [35].In our case, this condition is fulfilled only for the lower branch of anomalous solitons.It means there is a critical velocity β cr below which the anomalous soliton is unstable.
In figure 3(d) we show the stability criterion and see the solitons are stable only above a certain velocity β cr , marked with a red circle.If we compare this picture with figure 3(c), the upper subbranch corresponds to solitons moving with β ⩽ β cr (also marked with a red circles).Thus, the upper branch is the unstable one.
On the other hand, standard solitons with ρ ∞ > ρ eq , are always stable, which we can see in figure 3(b), namely β cr = 0. On the anomalous side, when we approach the standard regime (ρ ∞ → ρ − eq ), we have β cr → 0, whereas while getting closer to the unstable one (ρ ∞ → ρ + ins ), β cr → 1/2, which we show in figure 4. To give further evidence on the matter of the soliton stability, we numerically perform a Bogoliubov-de Gennes analysis of solitonic wave functions (see Appendix B for technical details).The lowest-state eigenenergy ω indicates whether or not the soliton is stable.We extract the critical velocity by finding the soliton velocity for which the BdG eigenvalue becomes real.The results are qualitatively consistent with the analysis based on the stability criterion from [35] (cf.figure 4).

IV. DARK SOLITON IN A QUANTUM DROPLET
The analysis of Bogoliubov excitations in quantum Bose-Bose droplets have revealed the presence of phonon modes [15].We now check whether solitons are present in quantum droplets as well.As the droplet is an inhomogeneous object, we use methods previously applied to trapped gases [36].We describe this approach below.
We take the following solution for a quantum droplet [12] ψ QD (x) = √ ρ eq µ/µ eq where µ eq = ρ eq − √ ρ eq is the chemical potential corresponding to a homogeneous density profile with density ρ = ρ eq and the number of particles in the droplet N as a function of the chemical potential µ is given by [15]: Next, we modify the droplet wave function in order to get a dark soliton inside it.We achieve it by simply multiplying (15) by the solitonic wave function ψ(x) = ρ(x)e iϕ(x) , where ρ(x) and ϕ(x) are the solutions of Eqs.(11a) and (11b) for ρ ∞ = max |ψ QD (x)| 2 and a given non-zero velocity β ̸ = 0. We will comment on the case β = 0 later.Then, we normalize the overall wave function ψ(x)ψ QD (x), impose periodic boundary conditions on the phase and numerically evolve in real time (see Appendix A for more details).
Figure 5 shows us the results of these simulations.The initial density and phase profiles of a moving solitons are shown in figure 5 The space-time diagram in figure 5(c) show us how the anomalous motionless and moving solitons behave in quantum droplets.The dynamics is very stable.We do not see any phonons or shock waves appearing.The grey soliton in figure 5(c) indeed travels at β = 0.5.When the droplet bulk is much larger than the soliton width, one should be able to use (16) along with the solution of Eqs.(11a) and (11b) to construct the wave function of a dark soliton-quantum droplet system.
Let us now return to the condition we made earlier about the non-zero soliton velocity.Are the motionless solitons not able to exist inside quantum droplets?According to our stability analysis, the motionless solitons are bound to be unstable.We emphasize, however, that one still can achieve very wide stable solitons inside quantum droplets having a small, but finite velocity β, as lim ρ∞→ρ − eq β cr = 0.

V. CONCLUSIONS
We have shown the basic properties of a single dark soliton in a weakly interacting Bose-Bose mixture.There are two types of solitons: standard and anomalous ones.The standard motionless solitons are black and they have a π-jump in the phase, whereas the motionless do not have a fully depleted density nor a phase jump.In the anomalous regime, the solitons have also a peculiar dispersion relation with a cusp and subbranch.We found there is a critical velocity of the solitons β cr , below which the solitons are unstable.The stability analysis plays a crucial role in the possible experimental realization of solitons in quantum droplets.First, the unstable solitons will be suppressed in a potential experiment.Second, it is an opportunity to study the instability mechanism.
Our study disfavours the hypothesis that anomalous solitons are a beyond-LHY feature.Our observations [29,30] are consistent with the condition for the presence of anomalous solitons [31,33], where a preferred equilibrium density ρ eq plays the key role in the existence of these objects.
Most importantly, we show that this uncanny type of solitons can exist inside a quantum droplet.According to our numerical simulations, the stable anomalous grey solitons might be experimentally observable in the droplets.
A theoretical study of multiple-soliton systems is a natural way to extend our research, especially if it goes for two-soliton interactions.
of soliton velocity for a series of different densities in order to estimate the transition from an unstable profile to a stable one; a negative value of ω 2 indicates the solitonic profile is unstable.It is worth noting the values obtained for β = 0 are consistent with [40].
The lack of quantitative agreement in figure 4 might be caused by the numerical imperfections, especially related to the numerical grid size.Here, we show results for a numerical grid with 2048 nodes, i.e. resulting in 4096 x 4096 matrix to diagonalize.However, it is noteworthy to mention that the results obtained with grid sizes of 512 and 1024 nodes exhibited comparatively poorer performance.This suggests that the choice of a coarser grid significantly impacted the accuracy of the outcomes, reinforcing the importance of an appropriately refined numerical grid for obtaining more reliable numerical results.As the grid size increases, so does the dimension of the matrix that needs to be diagonalized.This relationship highlights that larger grids entail more complex computational processes due to the higher-dimensional nature of the underlying matrices.This seems to be the main limitation of our BdG analysis.

FIG. 3 .
FIG. 3. (a,c) Dispersion relation of solitons [standard -top row, anomalous -bottom row].The red dashed line shows the upper subbranch.(b,d) Soliton stability criterion.Solitons are stable when dP/dβ < 0. In the standard regime, the solitons are stable everywhere.In the anomalous one -only above the critical velocity, which in this case of ρ∞ = 0.36 is numerically evaluated to be βcr = 0.2550(29).[Red circles in panels (c,d) mark the place, where the velocity of solitons equals to the critical velocity βcr.]

FIG. 4 .
FIG.4.Critical velocity based on the stability criterion dP/dβ < 0 (solid green line) and its estimation based on the Bogoliubovde Gennes analysis of solitonic wave function.