Modulational Instability, Inter-Component Asymmetry, and Formation of Quantum Droplets in One-Dimensional Binary Bose Gases

Abstract

Quantum droplets are ultradilute liquid states that emerge from the competitive interplay of two Hamiltonian terms, the mean-field energy and beyond-mean-field correction, in a weakly interacting binary Bose gas. We relate the formation of droplets in symmetric and asymmetric two-component one-dimensional boson systems to the modulational instability of a spatially uniform state driven by the beyond-mean-field term. Asymmetry between the components may be caused by their unequal populations or unequal intra-component interaction strengths. Stability of both symmetric and asymmetric droplets is investigated. Robustness of the symmetric solutions against symmetry-breaking perturbations is confirmed.

1. Introduction

The mean-field (MF) theory of weakly interacting dilute atomic gases rules out formation of a liquid state [1,2]. However, it has been recently shown that a liquid phase arises if one takes into account beyond-MF effects originating from quantum fluctuations around the MF ground state of weakly interacting binary (two-component) Bose gases [3]. A fundamental property that allows one to interpret this phase as a fluid is incompressibility: It maintains a limit density which cannot be made larger (see details below), hence adding more atoms leads to spatial expansion of the state. Another fundamental feature of this quantum-fluid phase is that it facilitates self-trapping of quantum droplets (QDs), which are stabilized by the interplay between the contact MF interaction and the beyond-MF Lee–Huang–Yang (LHY) correction [4]. Binary Bose–Einstein condensates (BECs) with competing intra- and inter-component MF interactions of opposite signs offer a remarkable possibility for the generation of QDs, as proposed by Petrov [3]. This possibility was further elaborated in various settings, including different effective dimensions [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. In particular, the dynamics of QDs with the flat-top (FT) or Gaussian shape, which correspond to large or relatively small numbers of particles, respectively, was addressed in the framework of the one-dimensional (1D) reduction of the model [20]. The theoretical prediction was followed by experimental creation of QDs in mixtures of two different atomic states of ${}^{39}$K, with quasi-2D [21,22] and fully 3D [23,24] shapes (see also recent reviews [25,26]). Very recently, the creation of especially long-lived QDs was reported in a heteronuclear ${}^{41}$K-${}^{87}$Rb system [27]. Another theoretically predicted and experimentally realized option for the creation of QDs makes use of the single-component condensate with dipole–dipole interactions [28,29,30,31,32,33,34,35]. It is relevant to mention that the formation of multiple droplets was also predicted and experimentally observed as an MF effect in strongly nonequilibrium (turbulent) states of BECs [36].

Collective modes of QDs are a subject of special interest, as they reveal internal dynamics of the droplets [20,24,32,37,38]. In particular, the stable existence of the QDs is secured if the particle-emission threshold lies below all excitation modes, hence a perturbation in the form of such modes will not cause decay of the droplet.

We aim to address issues that are related to the creation of QDs in the 1D setting and were not addressed in previous works. First, we consider modulational instability (MI) of spatially uniform plane-wave (PW) states, in the framework of the coupled system of Gross–Pitaevskii (GP) equations with the LHY corrections, for the two-component MF wave function of the binary condensate. This is the system which was originally derived in [5]. Recently, MI has been experimentally demonstrated in BECs with attractive interactions [39,40,41]. Other examples of the MI are provided by the binary BEC with the linear Rabi coupling or the spin-orbit coupling [41,42], and by a system combining the MF and LHY terms [43]. The linear-stability analysis, followed by direct simulations of the corresponding GP equations, shows that the lower branch of the PW states exhibits MI, the instability splitting the PW into a chain of localized droplet-like structures. Secondly, we address properties of the QDs in the binary condensate in the framework of the two-component GP system, without assuming effective inter-component symmetry, which reduces the system to a single-component GP equation. The asymmetry implies different MF self-repulsion coefficients in the two components, and/or unequal norms in them. Although properties of QDs have been studied by using the two-component GP system in some papers [6,11,14,15,17,18], the explicit asymmetry of the system parameters has not been addressed, except for [14] in which the situation for ${}^{39}$K-${}^{39}$K and ${}^{23}$Na-${}^{87}$Rb atomic mixtures have been considered. We conclude that the population difference between the components does not significantly affect density profiles of QDs in the system with equal MF self-repulsion strengths in the two components. On the other hand, we find that profiles of the QD solutions are essentially asymmetric when the self-repulsion coefficients are different in the components. Generally, the numerical findings corroborate stability of the known symmetric states against symmetry-breaking perturbations. We also address the MI of the two-component system, and demonstrate that chains of asymmetric QDs can be generated by the MI-induced nonlinear evolution.

The paper is organized as follows. In Section 2 we introduce the model and discuss conditions necessary for the formation of the droplets. Section 3.1 deals with the single-component version of the symmetric system. We consider various solutions admitted by it (PW, FT, periodic, etc.), and apply the linear-stability analysis of the PW solution to assess the MI, in a combination with direct simulations. In Section 3.2, we address the stability of asymmetric droplets, as well as the formation of droplets in the two-component asymmetric system via the MI. The paper is concluded by Section 4. Additional symmetric and asymmetric exact and approximate analytical solutions are presented in Appendices.

2. Model and Methods

We consider the 1D model of the two-component condensate with coefficients of the intra-component repulsion, $g_1>0$ and $g_2>0$, and inter-component attraction, $g_{12}<0$. In the weak-interaction limit, the corresponding energy density, which includes the MF terms and LHY correction, was derived in [5]:

$${\mathcal{E}}_{1\mathrm{D}}=\frac{{\left(\sqrt{g_1}{\rho }_1-\sqrt{g_2}{\rho }_2\right)}^2}{2}+\frac{g\delta g{\left(\sqrt{g_2}{\rho }_1+\sqrt{g_1}{\rho }_2\right)}^2}{{(g_1+g_2)}^2}-\frac{2\sqrt{m}{\left(g_1{\rho }_1+g_2{\rho }_2\right)}^{3/2}}{3\pi \hslash },$$

(1)

where m is the atomic mass (the same for both components), ${\rho }_j={|{\mathsf{\Psi }}_j|}^2$$(j=1,2)$ is the density of the j-th component, represented by the MF wave function ${\mathsf{\Psi }}_j$, and

$$g\equiv \sqrt{g_1{g}_2},\delta g\equiv g_{12}+g.$$

(2)

The last term in Equation (1) represents the LHY correction. Derivation of Equation (1) assumes that the binary BEC is close to the point of the MF repulsion-attraction balance, with $\left|\delta g\right|\ll g$. In experiments, $\delta g$ may be tuned to be both positive and negative [21,22,23].

Equation (1) is valid in the case of tight confinement applied in the transverse dimensions, which makes the setting effectively one-dimensional. In the 3D case, the LHY term $\sim -{\rho }^{3/2}$ (for ${\rho }_1={\rho }_2\equiv \rho$) is replaced by one $\sim +{\rho }^{5/2}$. A detailed consideration of the crossover from 3D to 1D [12,44,45] in the two-component system is a problem which may be a subject of a separate work. Here, it is relevant to compare the symmetric version of Equation (1) for the energy density with that recently presented in [12]. It demonstrates that an accurately derived LHY contribution to the energy density of the 1D system contains, in addition to the ${\rho }^{3/2}$ term which was derived in [5], a term $\sim {\rho }^2$, which can be absorbed into the mean-field energy density, and a higher-order term $\sim {\rho }^3$, which was omitted in the analysis reported in [12]. A conclusion formulated in that work is that the energy density originally derived in [5] is literally valid if the ratio of the mean-field energy to that of the transverse confinement takes values $\le 0.03$. For typical experimental parameters, this implies that the difference between absolute values of scattering lengths of the mean-field intra-component repulsion and inter-component attraction should be $\le 1$ nm, which may be achieved in the experiment. The 1D QDs originate from the balance of the second term in Equation (1), corresponding to the weakly repulsive MF interaction, with $\delta g>0$, and the LHY term, which introduces effective attraction in the 1D setting, on the contrary to the repulsion in the 3D setting [5,20].

The energy functional, ${\int }_{-\infty }^{+\infty }{\mathcal{E}}_{1\mathrm{D}}dZ$, gives rise to the system of GP equations, which include the LHY correction,

$$\begin{array}{cc}\hfill i\hslash \frac{\partial {\mathsf{\Psi }}_1}{\partial T}& =-\frac{{\hslash }^2}{2m}\frac{{\partial }^2{\mathsf{\Psi }}_1}{\partial Z^2}+(g_1+G{g}_2)|{\mathsf{\Psi }}_1{|}^2{\mathsf{\Psi }}_1-(1-G)g{|{\mathsf{\Psi }}_2|}^2{\mathsf{\Psi }}_1-\frac{g_1\sqrt{m}}{\pi \hslash }\sqrt{g_1|{\mathsf{\Psi }}_1{|}^2+g_2{|{\mathsf{\Psi }}_2|}^2}{\mathsf{\Psi }}_1,\hfill \\ \hfill i\hslash \frac{\partial {\mathsf{\Psi }}_2}{\partial T}& =-\frac{{\hslash }^2}{2m}\frac{{\partial }^2{\mathsf{\Psi }}_2}{\partial Z^2}+(g_2+G{g}_1)|{\mathsf{\Psi }}_2{|}^2{\mathsf{\Psi }}_2-(1-G)g{|{\mathsf{\Psi }}_1|}^2{\mathsf{\Psi }}_2-\frac{g_2\sqrt{m}}{\pi \hslash }\sqrt{g_1|{\mathsf{\Psi }}_1{|}^2+g_2{|{\mathsf{\Psi }}_2|}^2}{\mathsf{\Psi }}_2,\hfill \end{array}$$

(3)

where T and Z are the time and coordinate measured in physical units, and parameter

$$G=\frac{2g\delta g}{{(g_1+g_2)}^2},$$

(4)

measures the deviation from the MF repulsion–attraction balance point, see Equation (2). The normalization of the components of the wave function is determined by numbers of bosons in each component:

$$N_j={\int }_{-\infty }^{+\infty }{|{\mathsf{\Psi }}_j|}^{2}dZ.$$

(5)

Further, rescaling

$$\left(\frac{m{g}^2}{{\hslash }^3}\right)T\equiv t,\left(\frac{mg}{{\hslash }^2}\right)Z\equiv z,\left(\frac{\hslash }{\sqrt{mg}}\right){\mathsf{\Psi }}_{1,2}\equiv {\psi }_{1,2}$$

(6)

casts Equation (3) in the normalized form,

$$\begin{array}{cc}\hfill i\frac{\partial {\psi }_1}{\partial t}& =-\frac{1}{2}\frac{{\partial }^2{\psi }_1}{\partial z^2}+(P+G{P}^{-1})|{\psi }_1{|}^2{\psi }_1-(1-G){|{\psi }_2|}^2{\psi }_1-\frac{P}{\pi }\sqrt{P|{\psi }_1{|}^2+P^{-1}{|{\psi }_2|}^2}{\psi }_1,\hfill \\ \hfill i\frac{\partial {\psi }_2}{\partial t}& =-\frac{1}{2}\frac{{\partial }^2{\psi }_2}{\partial z^2}+(P^{-1}+GP)|{\psi }_2{|}^2{\psi }_2-(1-G){|{\psi }_1|}^2{\psi }_2-\frac{1}{\pi P}\sqrt{P^{-1}|{\psi }_2{|}^2+P{|{\psi }_1|}^2}{\psi }_2,\hfill \end{array}$$

(7)

where parameter

$$P\equiv \sqrt{\frac{g_1}{g_2}}=\frac{g_1}{g}$$

(8)

determines the asymmetry of the system, in the case of $P\ne 1$. Note that, as concerns stationary solutions with chemical potentials ${\mu }_{1,2}$, sought for as

$${\psi }_{1,2}\left(z,t\right)=exp(-i{\mu }_{1,2}t){\varphi }_{1,2}\left(z\right),$$

(9)

states with mutually proportional components, ${\varphi }_1\left(z\right)=K{\varphi }_2\left(z\right)$, are only possible in the fully symmetric case with $P=1$, ${\mu }_1={\mu }_2$, and $K=1$. In previous works [5,20], 1D solutions for QDs were considered only in the framework of the single GP equation which corresponds to symmetric system Equation (7) with $P=1$ and ${\psi }_1={\psi }_2$.

3. Modulation Instability Versus QDs

In this section, we address MI of PWs in both symmetric and asymmetric GP systems, and relate it to formation of the QDs in the binary bosonic gas. To the best of our knowledge, this is the first work aiming to associate the MI with the formation of the 1D droplets in the system with unequal components. We first consider MI in the framework of the single-component reduction of the symmetric version of Equation (7), after briefly reviewing stationary solutions of the GP equation. Next, we extend the analysis for the two-component GP system, which makes it possible to produce asymmetric QDs, starting from the MI of asymmetric PW states.

3.1. The Single-Component GP Model

Under the single-component reduction of the binary system, with $g_1=g_2\equiv g$ and ${\psi }_1={\psi }_2\equiv \psi$, Equation (1) simplifies to [5]

$${\epsilon }_{1\mathrm{D}}\equiv \frac{{\hslash }^4}{m^2{g}^3}{\mathcal{E}}_{1\mathrm{D}}=\frac{\delta g}{g}{n}^2-\frac{2^{5/2}}{3\pi }{n}^{3/2},$$

(10)

with the single dimensionless density, $n={|\psi |}^2\equiv \left({\hslash }^2/mg\right)\rho$. Assuming a spatially uniform state, the equilibrium density and the corresponding chemical potential are given by

$$n_0=\frac{8}{9{\pi }^2}{\left(\frac{g}{\delta g}\right)}^2,{\mu }_0=-\frac{4}{9{\pi }^2}\frac{g}{\delta g}.$$

(11)

Density $n_0$ corresponds to the minimum of the energy per particle, ${\partial }_n\left[n^{-1}{\epsilon }_{1\mathrm{D}}\left(n\right)\right]=0$, and ${\mu }_0$ is negative for $\delta g/g>0$. The corresponding single GP equation is

$$i\frac{\partial \psi }{\partial t}=-\frac{1}{2}\frac{{\partial }^2\psi }{\partial z^2}+\frac{\delta g}{g}{|\psi |}^2\psi -\frac{\sqrt{2}}{\pi }|\psi |\psi ,$$

(12)

with normalization condition ${\int }_{-\infty }^{+\infty }{|\psi \left(z\right)|}^{2}dz=N$, where $N\equiv N_1=N_2$ is the number of atoms in each component.

Although coefficient $\delta g/g$ can be scaled out in Equation (12), as done in [20], we keep it here as a free parameter. This option is convenient for the subsequent consideration of the MI, treating $\delta g/g$ and density n as independent constants, which may be matched to experimentally relevant parameters.

Below, we address two stationary solutions of Equation (12). One is the QD bound state of a finite size, which was studied in detail in [5,20]. The other solution is the PW with uniform density. Here, we briefly recapitulated basic properties of these solutions for the completeness of the presentation. In Section 3.1.3 we address the MI of the PWs and associate it with the spontaneous generation of chains of localized modes. Additional families of exact analytical solutions of Equation (12) are given in Appendix A.

3.1.1. The Droplet Solution

As shown in [5,20,46], at $\delta g/g>0$ Equation (12) gives rise to an exact soliton-like solution representing a QD, maintained by the balance between the effective cubic self-repulsion and quadratic attraction:

$$\psi (z,t)=\frac{A{e}^{-i\mu t}}{1+Bcosh\left(\sqrt{-2\mu }z\right)},A=\sqrt{n_0}\frac{\mu }{{\mu }_0},B=\sqrt{1-\frac{\mu }{{\mu }_0}}.$$

(13)

This solution exists in a finite range of negative values of the chemical potential ${\mu }_0<\mu <0$, featuring the FT shape at $0<\mu -{\mu }_0\ll \left|{\mu }_0\right|$, with size $L\approx {\left(-2{\mu }_0\right)}^{-1/2}ln\left[{\left(1-\mu /{\mu }_0\right)}^{-1}\right]$ [5,20]. A typical density profile of the FT solution is displayed in the inset of Figure 1. At $\mu ={\mu }_0$, the size of the droplet diverges, and the solution carries over into the delocalized PW with uniform density, $n=n_0$. The fact that the density of the condensate filling the FT state cannot exceed the largest value, $n_0$, implies its incompressibility. For this reason, the condensate may be considered as a fluid, as mentioned above. With the increase of $\mu$ from ${\mu }_0$ towards $\mu =0$, the maximum density of the localized mode,

$$n_{max}\equiv n(z=0)=n_0{\left(\frac{\mu }{{\mu }_0}\right)}^2{\left(1+\sqrt{1-\frac{\mu }{{\mu }_0}}\right)}^{-2},$$

(14)

monotonously decreases from $n_0$ to 0. The QD’s FWHM size, defined by condition $n\left(z=L_{\mathrm{FWHM}}/2\right)=n\left(z=0\right)/2$, also shrinks at first with increasing $\mu$, attaining a minimum value ${\left(L_{\mathrm{FWHM}}\right)}_{min}\approx 2.36/\sqrt{-{\mu }_0}$ at $\mu /{\mu }_0\approx 0.776$. Further increase of $\mu$ towards $\mu =0$ makes the QD broader, its width diverging as $L_{\mathrm{FWHM}}\approx 1.71/\sqrt{-\mu }$ at $\mu \to -0$.

The norm of the exact QD solutions given by Equation (13) is

$$N\left(\mu \right)=n_0\sqrt{-\frac{2}{{\mu }_0}}\left[ln\left(\frac{1+\sqrt{\mu /{\mu }_0}}{\sqrt{1-\mu /{\mu }_0}}\right)-\sqrt{\frac{\mu }{{\mu }_0}}\right].$$

(15)

It satisfies the well-known Vakhitov–Kolokolov (VK) necessary stability criterion [47],

$$\frac{dN\left(\mu \right)}{d\mu }=-\frac{n_0}{{\mu }_0^2}\sqrt{-\frac{\mu }{2}}\frac{1}{1-\mu /{\mu }_0}<0,$$

(16)

due to ${\mu }_0<0$ and

$$0<\mu /{\mu }_0<1.$$

(17)

Full stability of the QD family has been verified by direct simulations of the evolution of perturbed QDs in the framework of Equation (12).

It is relevant to mention that the exact solution of Equation (13) is valid too at $\delta g/g<0$, when the cubic term in Equation (12) is self-attractive, like the quadratic one. In that case, ${\mu }_0$ is positive, as per Equation (11), while the chemical potential of the self-trapped state remains negative, as the solution of Equation (13) may exist only at $\mu <0$. Then, it follows from Equation (13) that the soliton-like mode exists for all values of $\mu <0$ (unlike the finite interval Equation (17), in which the solution exists for $\delta g/g>0$), and it does not feature the FT shape. Rather, with the increase of $-\mu$, it demonstrates a crossover between the KdV-soliton shape $\sim {\mathrm{sech}}^2\left(\sqrt{-\mu /2}z\right)$ and the nonlinear-Schrödinger one, $\sim \mathrm{sech}\left(\sqrt{-2\mu }z\right)$. For $\delta g/g<0$, the $N\left(\mu \right)$ dependence for the soliton family carries over into the following form,

$$N\left(\mu \right){|}_{\delta g<0}=n_0\sqrt{\frac{2}{{\mu }_0}}\left[\sqrt{-\frac{\mu }{{\mu }_0}}-arctan\left(\sqrt{-\frac{\mu }{{\mu }_0}}\right)\right],$$

(18)

which is an analytical continuation of Equation (15). This dependence also satisfies the VK criterion.

3.1.2. The Plane-Wave Solution

The PW solution of Equation (12) can be presented in a form $\psi (z,t)=\sqrt{n}exp\left(i{K}_{\mathrm{PW}}z-i\mu t\right)$ with wavenumber $K_{\mathrm{PW}}$ and constant density n, which determine the corresponding chemical potential:

$${\mu }_{\mathrm{PW}}=\frac{\delta g}{g}n-\frac{\sqrt{2}}{\pi }\sqrt{n}+\frac{1}{2}{K}_{\mathrm{PW}}^2.$$

(19)

The Galilean invariance of Equation (12) implies that any quiescent solution ${\psi }_0\left(z,t\right)$ generates a family of moving ones, with arbitrary velocity c. Therefore, $K_{\mathrm{PW}}$ may be canceled by means of transformation ${\psi }_c\left(z,t\right)=exp\left(icz-i{c}^{2}t/2\right){\psi }_0\left(z-ct,t\right)$ with $c=-K_{\mathrm{PW}}$.

For given $\mu$, Equation (19) produces two different branches of the density as a function of $\mu$ (here, $K_{\mathrm{PW}}=0$ is set):

$$\sqrt{n^{\pm }\left(\mu \right)}=\frac{1}{\sqrt{2}\pi }\frac{g}{\delta g}\pm \sqrt{\frac{1}{2{\pi }^2}{\left(\frac{g}{\delta g}\right)}^2+\frac{g}{\delta g}\mu }.$$

(20)

For $\delta g/g=0.05$, these branches are shown in Figure 1. As follows from Equation (20), they exist (for $\delta g/g>0$) above a minimum value of $\mu$: ${\mu }_c=-{(2{\pi }^2\delta g/g)}^{-1}=(9/8){\mu }_0$, the respective density being

$$n_c=n^{\pm }\left({\mu }_c\right)=\frac{1}{2{\pi }^2}{\left(\frac{g}{\delta g}\right)}^2=\frac{9}{16}{n}_0.$$

(21)

Values $\mu ={\mu }_c$ and $n=n_c$ correspond to the spinodal point [5], and $n^{+}\left({\mu }_0\right)=n_0$ (see Equation (13)). Note that the above-mentioned existence region of the soliton solution in terms of the chemical potential, ${\mu }_0<\mu <0$, lies completely inside that of the PW state, which is ${\mu }_c\le \mu$. Thus, the soliton always coexists with the PW (this fact is also obvious in Figure 1).

3.1.3. Modulational Instability of the Plane Waves

Here, we aim to analyze the MI of PW solutions in the framework of the single-component GP Equation (12) and demonstrate how the development of the MI can help to generate QDs. We perform the analysis for the PWs with zero wavenumber $K_{\mathrm{PW}}=0$, which is sufficient due to the aforementioned Galilean invariance of the underlying equation.

A small perturbation is added to the stationary PW state as

$$\psi (z,t)=\left[\sqrt{n}+\delta \psi (z,t)\right]exp\left(-i\mu t\right).$$

(22)

The substitution of this expression in Equation (12) and linearization with respect to perturbation $\delta \psi$ leads to the corresponding Bogoliubov–de Gennes equation,

$$i\frac{\partial }{\partial t}\delta \psi =-\frac{1}{2}\frac{{\partial }^2}{\partial z^2}\delta \psi +\frac{\delta g}{g}n(\delta \psi +\delta {\psi }^{*})-\frac{\sqrt{n}}{\sqrt{2}\pi }(\delta \psi +\delta {\psi }^{*}).$$

(23)

By looking for perturbation eigenmodes with wavenumber k and frequency $\mathsf{\Omega }$,

$$\delta \psi =\zeta cos(kz-\mathsf{\Omega }t)+i\eta sin(kz-\mathsf{\Omega }t),$$

(24)

and real infinitesimal amplitudes $\zeta$ and $\eta$, Equation (23) yields a dispersion relation for the eigenfrequencies:

$${\mathsf{\Omega }}^2=\frac{k^4}{4}+\left(\frac{\delta g}{g}n-\frac{\sqrt{n}}{\sqrt{2}\pi }\right)k^2.$$

(25)

The MI takes place when $\mathsf{\Omega }$ acquires an imaginary part. As follows from Equation (25), this occurs when the density satisfies condition $n<{\left[2{\pi }^2{(\delta g/g)}^2\right]}^{-1}=n_c$ see (Equation (21)), which corresponds to branch $n^{-}$ of the PW state. The instability region in terms of k is given by

$$k^2<4\left(\frac{\sqrt{n}}{\sqrt{2}\pi }-\frac{\delta g}{g}n\right)\equiv k_0^2.$$

(26)

The MI gain $\sigma \equiv \left|\mathrm{Im}\mathsf{\Omega }\right|$ is plotted in Figure 2 versus $|k|$ and $\delta g/g$, for given density $n=40$ in panel (a), and versus $|k|$ and n, for given $\delta g/g=0.05$ in (b). It is easy to find from Equation (25) that the largest gain is attained at wavenumber

$$k_{max}=\frac{k_0}{\sqrt{2}},$$

(27)

with $k_0$ defined as per Equation (26). Note that Figure 2a includes the case of the self-attractive cubic nonlinearity, with $\delta g/g<0$, which naturally displays much stronger MI, as in this case it is driven by both the quadratic and cubic nonlinear terms. In fact, the extension of the MI chart to $\delta g/g<0$ makes it possible to compare the MI in the present system and its well-known counterpart in the setting with the fully attractive nonlinearity.

Comparing parameter values at which the QD solutions are predicted to appear, and the MI condition for the PW with the corresponding density, the MI is expected to provide a mechanism for the creation of the QDs. This is confirmed by direct simulations of the GP Equation (12), as shown in Figure 3. The PW with $n=10$ is taken as the input, so that it is subject to the MI for $\delta g/g=0.05$, as seen in Figure 2b. As shown in Figure 3, small initial perturbations trigger the emergence of multiple-QD patterns (chains) at $t\ge 100$. For these parameters, we get $k_{\max }=0.6508$ and $\sigma \left(k_{\max }\right)=0.2118$, which determines the wavelength of the fastest growing modulation, $\lambda =2\pi /k_{\max }\approx 9.66$, and the growth-time scale, $\tau =2\pi /\sigma \left(k_{\max }\right)\approx 30$. The number of the generated droplets in Figure 3 is consistent with estimate $L/\lambda \simeq 10$, where $L=100$ is the size of the simulation domain. We have checked that the number of generated droplets is approximately given by $L/\lambda$ for other values of parameters as well. This dynamical scenario is similar to those observed in other models in the course of the formation of soliton chains by MI of PWs [39,40]. The long-time evolution in Figure 3a shows that the number of the droplets becomes smaller due to merger of colliding droplets into a single one, which agrees with dynamical properties of 1D QDs reported in [20].

To implement this mechanism of the generation of a chain of solitons in the experiment, i.e., make the density smaller than the critical value $n_c$, one may either apply interaction quench (by means of the Feshbach resonance), suddenly decreasing $\delta g/g$, as was done in recent experimental works for different purposes [21,22,23,48]. Another option, which is specific to the 1D setting, is sudden decrease of density n by relaxing the transverse trapping.

3.2. The Two-Component Gross–Pitaevskii Model

In this section, we revert to the full two-component GP system Equation (7), aiming to explore the formation of QD states in it. The two-component setting may include parameter imbalance between the two components, as indicated theoretically [3] and observed experimentally [21,22,23,27]. Here, we present the analysis of asymmetric QDs in two cases: (i) the two-component GP system with different populations, $N_1/N_2\ne 1$, and equal intra-component coupling strength, $g_1=g_2$ (i.e., $P=1$, see Equation (8)), and (ii) the system with different intra-component coupling strengths, $g_1\ne g_2$ (i.e., $P\ne 1$). These options suggest a possibility to check the stability of the solutions of the symmetric system, reduced to the single-component form, against symmetry-breaking perturbations. That objective is relevant because, in the real experiment, scattering lengths of the self-interaction in the two components are never exactly equal [21,22,23,24]. We address, first, an asymmetric single-droplet solution, and, subsequently, MI of the PW states in the two-component system.

Because, as said above, solutions with mutually proportional components (written as ${\varphi }_1=K{\varphi }_2$) are possible solely in the strictly symmetric setting, asymmetric QDs cannot be found in an exact analytical form. As shown in Appendix B (see Equations (A6)–(A12)), asymptotic analytical solutions can be obtained for strongly asymmetric states, with one equation replaced by its linearized version. In this section, we chiefly rely on numerical solution of Equation (7).

3.2.1. Asymmetric QDs with Unequal Populations ($N_1\ne N_2$) for $g_1=g_2$ ($P=1$)

In the system with $P=1$ (see Equation (8)), we calculated the droplet states as stationary solutions of Equation (7) by means of the imaginary-time-evolution method with the Neumann’s boundary conditions, under the constraint that the norm is fixed in the first component, ${\int }_{-\infty }^{+\infty }dz{|{\psi }_1\left(z\right)|}^2=N_1$, while chemical potential ${\mu }_2$ is fixed in the other one, allowing its norm $N_2$ to vary.

Figure 4 displays essential features of weakly asymmetric droplets for $\delta g/g=0.05$ and fixed $N_1=100$. The symmetric (completely overlapping) solution with $N_1=N_2$ is found at ${\mu }_1={\mu }_2=-0.88878$. When ${\mu }_2$ deviates from this value, profiles of the two components become slightly different, as shown in Figure 4a. The profiles of the droplet solution hardly change for different values of ${\mu }_2$, but Figure 4b demonstrates that, at ${\mu }_2\to -0$, ${\psi }_2$ develops small-amplitude extended tails, which are absent in ${\psi }_1$. Due to the contribution of the tails, the approach of ${\mu }_2<0$ towards zero leads to the increase of norm $N_2$, as seen in Figure 4c. Note that the growth of $N_2\left({\mu }_2\right)$ at ${\mu }_2\to -0$ is opposite to the decay of the QD’s norm in the single-component model at $\mu \to -0$, cf. Equation (15). At ${\mu }_2\ge 0$, the ${\psi }_2$ component undergoes delocalization, with its tails developing a nonzero background at $|z|\to \infty$, as seen in the density profile displayed in Figure 4b for ${\mu }_2=0$, and norm $N_2\left({\mu }_2\right)$ diverging at ${\mu }_2\to -0$ in Figure 4c.

In Figure 4c, we also plot the parameter of the asymmetry between the two components, defined as

$${\delta }_{21}=\frac{n_2(z=0)-n_1(z=0)}{n_2(z=0)+n_1(z=0)}.$$

(28)

It increases almost linearly with ${\mu }_2$, although its absolute value does not exceed $0.02$. Thus, the droplet tends to keep a nearly symmetric profile, with respect to the two components, in the symmetric system, even if the population imbalance is admitted. In fact, this circumstance makes the analysis self-consistent, as the use of the GP system with the LHY correction implies that the MF intra- and inter-component interactions nearly cancel each other, which is possible only if shapes of the two components are nearly identical.

3.2.2. Asymmetric QDs in the System with $P\ne 1$ ($g_1\ne g_2$)

Next, we consider the QDs for $P\ne 1$, setting $P>1$ without loss of generality. Then, the MF energy is minimized for $n_2>n_1$; the situation with $n_1>n_2$ can be considered too, replacing P by $P^{-1}$.

Following the procedure similar to that employed in Section 3.2.1, we produce QD solutions for $\delta g/g=0.05$, $N_1=100$, and several different values of P, varying ${\mu }_2$. In Figure 5a, we plot density profiles for three different values of P. Naturally, the difference of the two components increases with the increase of P. In Figure 5b we display $N_2$ and parameter ${\delta }_{21}$ (see Equation (28)) of the asymmetric QDs for $P=1.25$ and $1.67$. All these states have been checked to be stable in time-dependent simulations.

The density difference at the center of the droplet can be determined by the condition of the existence of the liquid phase in the free space. This condition is obtained by minimizing the grand-potential density ${\mathcal{E}}_{1\mathrm{D}}-{\mu }_1{\rho }_1-{\mu }_2{\rho }_2$ [5,14], which leads to the zero-pressure condition,

$$\begin{array}{cc}\hfill p({\rho }_1,{\rho }_2)& =-{\mathcal{E}}_{1\mathrm{D}}+\sum _{j=1,2}\left(\frac{\partial {\mathcal{E}}_{1\mathrm{D}}}{\partial {\rho }_j}\right){\rho }_j\hfill \\ \hfill & \equiv -{\mathcal{E}}_{1\mathrm{D}}+{\mu }_1{\rho }_1+{\mu }_2{\rho }_2=0.\hfill \end{array}$$

(29)

From this, we obtain relation

$$\frac{{(\sqrt{g_1}{\rho }_1-\sqrt{g_2}{\rho }_2)}^2}{2}+\frac{g\delta g{(\sqrt{g_2}{\rho }_1+\sqrt{g_1}{\rho }_2)}^2}{{(g_1+g_2)}^2}-\frac{\sqrt{m}}{3\pi \hslash }{(g_1{\rho }_1+g_2{\rho }_2)}^{3/2}=0,$$

(30)

which can be rewritten in the scaled form as

$$\frac{P+G{P}^{-1}}{2}{n}_1^2+\frac{P^{-1}+GP}{2}{n}_2^2+(G-1)n_1{n}_2=\frac{1}{3\pi }{\left(P{n}_1+\frac{n_2}{P}\right)}^{3/2}.$$

(31)

For given $n_1$, we solved Equation (31) to find the respective value of $n_2$, which is shown in Figure 6 for $\delta g/g=0.05$ and several values of P. There are two branches of the solutions, that enclose the negative-pressure region, in which QDs may exist. The maximum value of $n_j$ at the tip of the negative-pressure region corresponds to the density in the droplet’s FT segment. The ascending negative-pressure region for each P nearly follows relation $n_2=P{n}_1$, which is derived by the minimization condition for the dominant first term in Equation (30) It is seen that a larger difference in the profiles of the two components occurs for larger P, as expected. Also, for given $n_1$, the negative-pressure region becomes wider with respect to $n_2$ for larger P (note that the figure displays a log–log plot).

As the QDs have a finite norm, it is relevant to characterize the asymmetry in terms of the norm, rather than density. Here, we aim to find a largest value of the norm difference,

$${\Delta }_{21}=(N_2-N_1)/N_{\mathrm{T}},$$

(32)

where $N_{\mathrm{T}}=N_1+N_2$ is the total norm, which admits the existence of the QDs. For given $N_1$, we obtain the upper bound for $N_2$ above which the solution becomes delocalized, and calculate the corresponding critical value of ${\Delta }_{21}$. The results are shown in Figure 7. For the system with $P=1$, the curve demonstrates an empirical dependence ${\Delta }_{21}\propto N_{\mathrm{T}}^{-\alpha }$ with exponent $\alpha \approx 0.58$. Accordingly, the asymmetry tends to vanish asymptotically for very “heavy” droplets, at $N_{\mathrm{T}}\to \infty$. As the system becomes slightly asymmetric, with $P=1.25$, exponent $\alpha$ is significantly reduced for small $N_{\mathrm{T}}$, and converges to a certain finite value at $N_{\mathrm{T}}\to \infty$. Thus, it is again confirmed that values $P>1$ maintain conspicuous asymmetry between the QD’s components. Finally, strongly asymmetric non-FT (Gaussian-shaped [20]) solutions can be obtained in an approximate analytical form for any value of P, as shown in Appendix C.

3.2.3. The MI of the Asymmetric PW States

The MI of two-component asymmetric PWs is a relevant subject too. Such solutions are written as ${\psi }_j(z,t)=\sqrt{n_j}{e}^{-i{\mu }_{j}t},(j=1,2)$. The substitution of this in Equation (7) yields

$$\begin{array}{ccc}\hfill {\mu }_1& =& (P+G{P}^{-1})n_1+(-1+G)n_2-\frac{P}{\pi }\sqrt{P{n}_1+\frac{n_2}{P}},\hfill \\ \hfill {\mu }_2& =& (P^{-1}+GP)n_2+(-1+G)n_1-\frac{1}{\pi P}\sqrt{P{n}_1+\frac{n_2}{P}}.\hfill \end{array}$$

(33)

Accordingly, in the symmetric system with $P=1$, densities of the asymmetric PW state are expressed in terms of the chemical potentials as

$$n_j=\frac{1}{4}\left[\frac{1}{{\pi }^2{G}^2}+\frac{{\mu }_1+{\mu }_2}{G}+{(-1)}^{j+1}({\mu }_1-{\mu }_2)\right]\pm \frac{\sqrt{1+2{\pi }^{2}G({\mu }_1+{\mu }_2)}}{4{\pi }^2{G}^2}.$$

We introduce the perturbation around the PW states as

$${\psi }_j(z,t)=\left[\sqrt{n_j}+\delta {\psi }_j(z,t)\right]e^{-i{\mu }_{j}t},$$

(34)

$$\delta {\psi }_j={\zeta }_{j}cos(kz-\mathsf{\Omega }t)+i{\eta }_{j}sin(kz-\mathsf{\Omega }t),$$

(35)

with infinitesimal amplitudes ${\zeta }_j$ and ${\eta }_j$, cf. Equation (24). The substitution of this in Equation (7) and the linearization with respect to ${\zeta }_j$ and ${\eta }_j$ yields the dispersion equation for the perturbation:

$${\mathsf{\Omega }}_{\pm }^2=\frac{k^2}{4}\left[k^2+2(P_1+P_2-Q_1-Q_2)\right]\pm \frac{k^2}{2}\sqrt{{(P_1-P_2-Q_1+Q_2)}^2+4{(R-S)}^2},$$

(36)

where

$$\begin{array}{c}{P}_1=(P+G{P}^{-1})n_1,P_2=(P^{-1}+GP)n_2,\\ Q_1=\frac{P^2{n}_1}{2\pi \sqrt{P{n}_1+P^{-1}{n}_2}},Q_2=\frac{P^{-2}{n}_2}{2\pi \sqrt{P{n}_1+P^{-1}{n}_2}}.\\ R=(-1+G)\sqrt{n_1{n}_2},S=\frac{\sqrt{n_1{n}_2}}{2\pi \sqrt{P{n}_1+P^{-1}{n}_2}},\end{array}$$

(37)

For $P=1$ and $n_1=n_2$, these results reproduce Equation (25) for the ${\mathsf{\Omega }}_{-}$ branch. A parameter region in which at least one squared eigenfrequency ${\mathsf{\Omega }}_{\pm }^2$ is negative gives rise to the MI of the two-component state.

3.2.4. The MI for $P=1$

In Figure 8, we plot the gain spectrum $\sigma =\mathrm{Im}\left(\mathsf{\Omega }\right)$ for the asymmetric PWs in the symmetric system with $P=1$ and $\delta g/g=0.05$, in the plane of wavenumber k and density ratio $n_{12}=n_2/n_2$. For the consistency with the single-component situation displayed in Figure 3, we here fix the total density as $(n_1+n_2)/2=10$. For given $n_{12}$, the MI occurs at $|k|<k_0$, and the gain attains its maximum at $k=k_{\max }=k_0/\sqrt{2}$. The largest gain is obtained at equal densities, $n_{12}=1$. Both the k-band of the instability and magnitude of the gain slowly decrease as the deviation of $n_{12}$ from unity increases. This means that the MI occurs in the PW states with a large density difference, thus giving rise to the formation of solitons with large asymmetry even for equal intra-component MF interaction strengths, $P=1$ (see Equation (8)).

In Figure 9 we display typical examples of the numerically simulated development of the MI in the symmetric two-component system with $P=1$ and population imbalance. Figure 9a shows the evolution of central-point values of the density of the first component, $n_1(z=0)$, for different values of the density ratio, $n_{12}=n_1/n_2$. Time required for the actual onset of the instability increases with the increase in $n_{12}$, as is clearly shown by the density-plot evolution in Figure 9d,e for $n_{12}=1$ and Figure 9f,g for $n_{12}=9$. This observation can be understood in terms of the MI gain $\sigma$, as shown in Figure 8c, where $\sigma$ at $k=k_{\max }$ becomes smaller with increasing $n_{12}$.

Spatial profiles at fixed time, which are plotted in Figure 9b,c for these two cases, show fragmentation of the profiles into sets of localized structures. The decrease in the number of fragments with the increase of $n_{12}$ is explained by the decrease of $k_{max}$, see Figure 8a. For $n_{12}=1$, the results are the same as in the single-component case, as coinciding profiles in the two components of the symmetric system are stable against spontaneous symmetry breaking. On the other hand, when $n_{12}\ne 1$ an in-phase two-component localized structure appears, keeping the initial density imbalance. Since one can select an arbitrary ratio of densities of the two components for the initial PW state, a highly asymmetric structure, like the one displayed in Figure 9c, may emerge even for $P=1$, as a result of the MI-induced nonlinear evolution.

3.2.5. The MI for $P\ne 1$

Figure 8b represents the MI gain for $P=1.25$ and a fixed total density, $(n_1+n_2)/2=10$, in the case of slightly different strengths of the intra-component repulsion. The peak value of the MI gain is attained at $n_{12}=n_{12}^{\max }=0.577$, below the equal-densities point $n_{12}=1$. This is consistent with the fact that, at $P>1$, unequal values $n_1<n_2$ are suitable to the formation of an asymmetric soliton structure, as seen in Figure 5a. In Figure 8c, we plot the peak MI gain, $\sigma \left(k_{\max }\right)$, along with the respective value of the density ratio, $n_{12}^{\max }$, as a function of P. Value $n_{12}^{\max }$ monotonously decreases as a function of P, while the peak gain attains a minimum at $P=1$.

In Figure 10, we present the development of the MI in the two-component system for $P=1.25$ and a fixed total density, $(n_1+n_2)/2=10$. Figure 10a displays the evolution of the central-point density of the first component, $n_1(z=0)$, for different values of the density ratio, $n_{12}=n_1/n_2$. It shows that time required for the development of the MI increases with the increase in the asymmetry of the density. This is also made evident by the density plots of the temporal evolution of the first component, shown in Figure 10e–g. This result is consistent with Equation (36), which shows a decrease of the MI gain with the increase of the asymmetry even for $P\ne 1$. Spatial profiles at fixed time, displayed in Figure 10b–d, show fragmentation of the profiles. Figure 10c clearly indicates that, even for $n_{12}=1$, the MI generates asymmetric droplet-like structures similar to Figure 5a, where the complete overlapping of the two densities does not occur.

4. Conclusions

The main purpose of this work is to associate the MI (modulation instability) of plane waves (PWs) to the mechanism of the creation of QDs (quantum droplets) in the system described by the coupled GP (Gross–Pitaevskii) equations including the LHY (Lee–Huang–Yang) terms in the 1D setting. This system is the model of weakly interacting binary Bose gases with approximately balanced interactions between the intra-component self-repulsion and the inter-component attraction. We have investigated, analytically and numerically, the MI of the lower branch of PW states in both symmetric (effectively single-component) and asymmetric (two-component) GP systems, and ensuing formation of a chain of droplet-like states. In particular, numerical solution for QDs which are asymmetric with respect to the two components are obtained, both in the system with equal repulsion strengths but unequal populations in the two components, and in the one with different self-repulsion strengths. The results corroborate that the previously known symmetric states are robust against symmetry-breaking disturbances.

These predictions can be tested experimentally by preparing uniform binary Bose gases with equal or different densities of two components, and suddenly reducing the strength of the effective MF (mean-field) interaction by means of the Feshbach-resonance quench, in order to enhance the relative strength of the LHY terms. In particular, for typical values of physical parameters, an estimate of the characteristic time of the modulation instability growth for typical values of the physical parameters is ∼1 $\mathsf{\mu }$s. This time is much smaller than a typical lifetime of the droplet, which is ≳100 $\mathsf{\mu }$s [21,22,23,27], thus making the observation of the MI feasible. The present analysis being restricted to the 1D setting, effects of the tight transverse confinement and crossover to the 3D configuration [12,44,45] deserves further consideration.