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

#### **Thudiyangal Mithun 1,2,\* ,†, Aleksandra Maluckov 3,\* ,†, Kenichi Kasamatsu 4,† , Boris A. Malomed 5,† and Avinash Khare 6,†**


Received: 9 December 2019; Accepted: 14 January 2020; Published: 18 January 2020

**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.

**Keywords:** quantum droplet; binary Bose–Einstein condensate; modulational instability

#### **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–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 <sup>39</sup>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 <sup>41</sup>K-87Rb

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–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–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 <sup>39</sup>K-39K and <sup>23</sup>Na-87Rb 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*<sup>1</sup> > 0 and *g*<sup>2</sup> > 0, and inter-component attraction, *g*<sup>12</sup> < 0. In the weak-interaction limit, the corresponding energy density, which includes the MF terms and LHY correction, was derived in [5]:

$$\mathcal{E}\_{\rm ID} = \frac{\left(\sqrt{\mathcal{g}\_1}\rho\_1 - \sqrt{\mathcal{g}\_2}\rho\_2\right)^2}{2} + \frac{\mathcal{g}\delta\mathfrak{g}\left(\sqrt{\mathcal{g}\_2}\rho\_1 + \sqrt{\mathcal{g}\_1}\rho\_2\right)^2}{\left(\mathcal{g}\_1 + \mathcal{g}\_2\right)^2} - \frac{2\sqrt{\mathcal{m}}\left(\mathcal{g}\_1\rho\_1 + \mathcal{g}\_2\rho\_2\right)^{3/2}}{3\pi\hbar},\tag{1}$$

where *m* is the atomic mass (the same for both components), *ρ<sup>j</sup>* = |Ψ*<sup>j</sup>* | 2 (*j* = 1, 2) is the density of the *j*-th component, represented by the MF wave function Ψ*<sup>j</sup>* , and

$$\mathbf{g} \equiv \sqrt{\mathbf{g}\_1 \mathbf{g}\_2} \qquad \delta \mathbf{g} \equiv \mathbf{g}\_{12} + \mathbf{g}. \tag{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 |*δg*| *g* . In experiments, *δg* may be tuned to be both positive and negative [21–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 ∼ −*ρ* 3/2 (for *<sup>ρ</sup>*<sup>1</sup> <sup>=</sup> *<sup>ρ</sup>*<sup>2</sup> <sup>≡</sup> *<sup>ρ</sup>*) is replaced by one ∼ +*ρ* 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 *ρ* 3/2 term which was derived in [5], a term <sup>∼</sup> *<sup>ρ</sup>* 2 , which can be absorbed into the mean-field energy density, and a higher-order term ∼ *ρ* 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 ≤ 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 ≤ 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 *δ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, R <sup>+</sup><sup>∞</sup> −∞ E1D*dZ*, gives rise to the system of GP equations, which include the LHY correction,

$$\begin{split} i\hbar \frac{\partial \Psi\_{1}}{\partial T} &= -\frac{\hbar^{2}}{2m} \frac{\partial^{2} \Psi\_{1}}{\partial Z^{2}} + (\mathcal{g}\_{1} + \mathbb{G}\mathcal{g}\_{2})|\Psi\_{1}|^{2}\Psi\_{1} - (1 - \mathbb{G})\mathcal{g}|\Psi\_{2}|^{2}\Psi\_{1} - \frac{\mathcal{g}\_{1}\sqrt{m}}{\pi\hbar} \sqrt{\mathcal{g}\_{1}|\Psi\_{1}|^{2} + \mathcal{g}\_{2}|\Psi\_{2}|^{2}\Psi\_{1}} \\ i\hbar \frac{\partial \Psi\_{2}}{\partial T} &= -\frac{\hbar^{2}}{2m} \frac{\partial^{2} \Psi\_{2}}{\partial Z^{2}} + (\mathcal{g}\_{2} + \mathbb{G}\mathcal{g}\_{1})|\Psi\_{2}|^{2}\Psi\_{2} - (1 - \mathbb{G})\mathcal{g}|\Psi\_{1}|^{2}\Psi\_{2} - \frac{\mathcal{g}\_{2}\sqrt{m}}{\pi\hbar} \sqrt{\mathcal{g}\_{1}|\Psi\_{1}|^{2} + \mathcal{g}\_{2}|\Psi\_{2}|^{2}\Psi\_{2}} \end{split} (3)$$

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

$$\mathbf{G} = \frac{\mathbf{2g\delta g}}{(\mathbf{g\_1} + \mathbf{g\_2})^2} \tag{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\_{\hat{j}} = \int\_{-\infty}^{+\infty} |\Psi\_{\hat{j}}|^2 dZ. \tag{5}$$

Further, rescaling

$$\left(\frac{\text{mg}^2}{\hbar^3}\right)T \equiv t, \quad \left(\frac{\text{mg}}{\hbar^2}\right)Z \equiv z, \quad \left(\frac{\hbar}{\sqrt{\text{mg}}}\right)\Psi\_{1,2} \equiv \psi\_{1,2} \tag{6}$$

casts Equation (3) in the normalized form,

$$\begin{aligned} i\frac{\partial\psi\_1}{\partial t} &= -\frac{1}{2}\frac{\partial^2\psi\_1}{\partial z^2} + (P + GP^{-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, \\ 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. \end{aligned} \tag{7}$$

where parameter

$$P \equiv \sqrt{\frac{\mathcal{G}\_1}{\mathcal{G}\_2}} = \frac{\mathcal{G}\_1}{\mathcal{g}} \tag{8}$$

determines the asymmetry of the system, in the case of *P* 6= 1. Note that, as concerns stationary solutions with chemical potentials *µ*1,2, sought for as

$$
\psi\_{1,2}(z,t) = \exp(-i\mu\_{1,2}t)\phi\_{1,2}(z),
\tag{9}
$$

states with mutually proportional components, *φ*1(*z*) = *Kφ*2(*z*), are only possible in the fully symmetric case with *P* = 1, *µ*<sup>1</sup> = *µ*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 *ψ*<sup>1</sup> = *ψ*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*<sup>1</sup> = *g*<sup>2</sup> ≡ *g* and *ψ*<sup>1</sup> = *ψ*<sup>2</sup> ≡ *ψ*, Equation (1) simplifies to [5]

$$
\varepsilon\_{\rm{1D}} \equiv \frac{\hbar^4}{m^2 \mathcal{g}^3} \mathcal{E}\_{\rm{1D}} = \frac{\delta \mathcal{g}}{\mathcal{g}} n^2 - \frac{2^{5/2}}{3\pi} n^{3/2},\tag{10}
$$

with the single dimensionless density, *n* = |*ψ*| <sup>2</sup> <sup>≡</sup> *h*¯ <sup>2</sup>/*mg ρ*. Assuming a spatially uniform state, the equilibrium density and the corresponding chemical potential are given by

$$m\_0 = \frac{8}{9\pi^2} \left(\frac{\mathcal{g}}{\delta \mathcal{g}}\right)^2, \; \mu\_0 = -\frac{4}{9\pi^2} \frac{\mathcal{g}}{\delta \mathcal{g}}.\tag{11}$$

Density *n*<sup>0</sup> corresponds to the minimum of the energy per particle, *∂<sup>n</sup>* - *n* −1 *ε*1D(*n*) = 0, and *µ*<sup>0</sup> is negative for *δ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\tag{12}$$

with normalization condition R <sup>+</sup><sup>∞</sup> −∞ |*ψ*(*z*)| <sup>2</sup>*dz* <sup>=</sup> *<sup>N</sup>*, where *<sup>N</sup>* <sup>≡</sup> *<sup>N</sup>*<sup>1</sup> <sup>=</sup> *<sup>N</sup>*<sup>2</sup> is the number of atoms in each component.

Although coefficient *δ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 *δ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] and [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 *δ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{Ae^{-i\mu t}}{1 + B\cosh(\sqrt{-2\mu}z)}, \quad A = \sqrt{n\_0}\frac{\mu}{\mu\_0}, \quad B = \sqrt{1 - \frac{\mu}{\mu\_0}}.\tag{13}$$

This solution exists in a finite range of negative values of the chemical potential *µ*<sup>0</sup> < *µ* < 0, featuring the FT shape at 0 < *µ* − *µ*<sup>0</sup> |*µ*0|, with size *L* ≈ (−2*µ*0) <sup>−</sup>1/2 ln <sup>h</sup> (1 − *µ*/*µ*0) −1 i [5,20]. A typical density profile of the FT solution is displayed in the inset of Figure 1. At *µ* = *µ*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 *µ* from *µ*<sup>0</sup> towards *µ* = 0, the maximum density of the localized mode,

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

monotonously decreases from *n*<sup>0</sup> to 0. The QD's FWHM size, defined by condition *n* (*z* = *L*FWHM/2) = *n* (*z* = 0) /2, also shrinks at first with increasing *µ*, attaining a minimum value (*L*FWHM)min ≈ 2.36/ √ −*µ*<sup>0</sup> at *µ*/*µ*<sup>0</sup> ≈ 0.776. Further increase of *µ* towards *µ* = 0 makes the QD broader, its width diverging as *<sup>L</sup>*FWHM <sup>≈</sup> 1.71/<sup>√</sup> −*µ* at *µ* → −0.

**Figure 1.** The maximum density *n*max ≡ *n*(*z* = 0) in the FT (flat-top) state, as per Equation (14), and the PW (plane-wave)/density are displayed as functions of *µ* by the red-solid and blue-dashed curves, respectively, for *δg*/*g* = 0.05. In this case, Equation (11) yields *n*<sup>0</sup> = 36.025 and *µ*<sup>0</sup> = −0.900633. The PW solution includes upper and lower branches corresponding to *n* ±, as given by Equation (20), the lower one (marked by circles) being subject to the MI (modulational instability). The spinodal point is one with coordinates (*µc*, *nc*). For other values of *δg*/*g*, the plot can be generated from the present one by rescaling. The inset shows the density profile of the FT solution for *δg*/*g* = 0.05 and *µ* = *µ*<sup>0</sup> + 0.00001, very close to the delocalization limit (the transition to PW).

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

$$N(\mu) = 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]. \tag{15}$$

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

$$\frac{dN(\mu)}{d\mu} = -\frac{n\_0}{\mu\_0^2} \sqrt{-\frac{\mu}{2}} \frac{1}{1 - \mu/\mu\_0} < 0,\tag{16}$$

due to *µ*<sup>0</sup> < 0 and

$$0 < \mu/\mu\_0 < 1.\tag{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 *δg*/*g* < 0, when the cubic term in Equation (12) is self-attractive, like the quadratic one. In that case, *µ*<sup>0</sup> 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 *µ* < 0. Then, it follows from Equation (13) that the soliton-like mode exists for all values of *µ* < 0 (unlike the finite interval Equation (17), in which the solution exists for *δg*/*g* > 0), and it does not feature the FT shape. Rather, with the increase of −*µ*, it demonstrates a crossover between the KdV-soliton shape <sup>∼</sup> sech<sup>2</sup> p −*µ*/2*z* and the nonlinear-Schrödinger one, <sup>∼</sup> sech p −2*µz* . For *δg*/*g* < 0, the *N*(*µ*) dependence for the soliton family carries over into the following form,

$$N(\mu)\Big|\_{\delta\mathcal{G}<0} = n\_0 \sqrt{\frac{2}{\mu\_0}} \left[ \sqrt{-\frac{\mu}{\mu\_0}} - \arctan\left(\sqrt{-\frac{\mu}{\mu\_0}}\right) \right],\tag{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 *<sup>ψ</sup>*(*z*, *<sup>t</sup>*) = <sup>√</sup> *n* exp (*iK*PW*z* − *iµt*) with wavenumber *K*PW and constant density *n*, which determine the corresponding chemical potential:

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

The Galilean invariance of Equation (12) implies that any quiescent solution *ψ*<sup>0</sup> (*z*, *t*) generates a family of moving ones, with arbitrary velocity *c*. Therefore, *K*PW may be canceled by means of transformation *ψ<sup>c</sup>* (*z*, *t*) = exp *icz* <sup>−</sup> *ic*<sup>2</sup> *t*/2 *ψ*<sup>0</sup> (*z* − *ct*, *t*) with *c* = −*K*PW.

For given *µ*, Equation (19) produces two different branches of the density as a function of *µ* (here, *K*PW = 0 is set):

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

For *δg*/*g* = 0.05, these branches are shown in Figure 1. As follows from Equation (20), they exist (for *δg*/*g* > 0) above a minimum value of *µ*: *µ<sup>c</sup>* = −(2*π* 2 *δg*/*g*) <sup>−</sup><sup>1</sup> = (9/8)*µ*0, the respective density being

$$n\_{\varepsilon} = n^{\pm}(\mu\_{\varepsilon}) = \frac{1}{2\pi^2} \left(\frac{\text{g}}{\delta g}\right)^2 = \frac{9}{16} n\_0. \tag{21}$$

Values *µ* = *µ<sup>c</sup>* and *n* = *n<sup>c</sup>* correspond to the spinodal point [5], and *n* <sup>+</sup>(*µ*0) = *n*<sup>0</sup> (see Equation (13)). Note that the above-mentioned existence region of the soliton solution in terms of the chemical

potential, *µ*<sup>0</sup> < *µ* < 0, lies completely inside that of the PW state, which is *µ<sup>c</sup>* ≤ *µ*. 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*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). \tag{22}$$

The substitution of this expression in Equation (12) and linearization with respect to perturbation *δψ* 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^\*). \tag{23}$$

By looking for perturbation eigenmodes with wavenumber *k* and frequency Ω,

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

and real infinitesimal amplitudes *ζ* and *η*, Equation (23) yields a dispersion relation for the eigenfrequencies:

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

The MI takes place when Ω acquires an imaginary part. As follows from Equation (25), this occurs when the density satisfies condition *n* < [2*π* 2 (*δg*/*g*) 2 ] <sup>−</sup><sup>1</sup> = *n<sup>c</sup>* 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\mathfrak{g}}{\mathfrak{g}}n\right) \equiv k\_0^2. \tag{26}$$

The MI gain *σ* ≡ |ImΩ| is plotted in Figure 2 versus |*k*| and *δg*/*g*, for given density *n* = 40 in panel (a), and versus |*k*| and *n*, for given *δg*/*g* = 0.05 in (b). It is easy to find from Equation (25) that the largest gain is attained at wavenumber

$$k\_{\text{max}} = \frac{k\_0}{\sqrt{2}},\tag{27}$$

with *k*<sup>0</sup> defined as per Equation (26). Note that Figure 2a includes the case of the self-attractive cubic nonlinearity, with *δ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 *δ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 *δ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* ≥ 100. For these parameters, we get *k*max = 0.6508 and *σ* (*k*max) = 0.2118, which determines the wavelength of the fastest growing modulation, *λ* = 2*π*/*k*max ≈ 9.66, and the growth-time scale, *τ* = 2*π*/*σ* (*k*max) ≈ 30. The number of the generated droplets in Figure 3 is

consistent with estimate *L*/*λ* ' 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*/*λ* 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].

**Figure 2.** Color-coded values of the MI gain, *σ* = Im(Ω), are displayed for fixed *n* = 40 in (**a**), and for fixed *δg*/*g* = 0.05 in (**b**). Note that panel (**a**) covers both signs of the cubic nonlinearity, *δg* > 0 and *δg* < 0. Solid and dashed white curves represent the MI boundary (Equation (26)) and the peak value of the MI gain (Equation (27)), respectively.

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 *nc*, one may either apply interaction quench (by means of the Feshbach resonance), suddenly decreasing *δg*/*g*, as was done in recent experimental works for different purposes [21–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–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*<sup>2</sup> 6= 1, and equal intra-component coupling strength, *g*<sup>1</sup> = *g*<sup>2</sup> (i.e., *P* = 1, see Equation (8)), and (ii) the system with different intra-component coupling strengths, *g*<sup>1</sup> 6= *g*<sup>2</sup> (i.e., *P* 6= 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–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 *φ*<sup>1</sup> = *Kφ*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).

**Figure 3.** A typical example of the MI development, starting from an unstable PW state, with density *n* = 10 and *δg*/*g* = 0.05, which is subject to the MI, pursuant to Figure 2. In (**a**), the spatiotemporal pattern of the evolution of the condensate density is shown. In the right-hand panels, cross sections of the density profiles are displayed at *t* = 100 (**b**), *t* = 120 (**c**), and *t* = 140 (**d**). The simulations were performed in domain −50 < *z* < +50 with 2500 grid points and periodic boundary conditions.

3.2.1. Asymmetric QDs with Unequal Populations (*N*<sup>1</sup> 6= *N*2) for *g*<sup>1</sup> = *g*<sup>2</sup> (*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, R <sup>+</sup><sup>∞</sup> −∞ *dz*|*ψ*1(*z*)| <sup>2</sup> = *N*1, while chemical potential *µ*<sup>2</sup> is fixed in the other one, allowing its norm *N*<sup>2</sup> to vary.

Figure 4 displays essential features of weakly asymmetric droplets for *δg*/*g* = 0.05 and fixed *N*<sup>1</sup> = 100. The symmetric (completely overlapping) solution with *N*<sup>1</sup> = *N*<sup>2</sup> is found at *µ*<sup>1</sup> = *µ*<sup>2</sup> = −0.88878. When *µ*<sup>2</sup> 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 *µ*2, but Figure 4b demonstrates that, at *µ*<sup>2</sup> → −0, *ψ*<sup>2</sup> develops small-amplitude extended tails, which are absent in *ψ*1. Due to the contribution of the tails, the approach of *µ*<sup>2</sup> < 0 towards zero leads to the increase of norm *N*2, as seen in Figure 4c. Note that the growth of *N*2(*µ*2) at *µ*<sup>2</sup> → −0 is opposite to the decay of the QD's norm in the single-component model at *µ* → −0, cf. Equation (15). At *µ*<sup>2</sup> ≥ 0, the *ψ*<sup>2</sup>

component undergoes delocalization, with its tails developing a nonzero background at |*z*| → ∞, as seen in the density profile displayed in Figure 4b for *µ*<sup>2</sup> = 0, and norm *N*2(*µ*2) diverging at *µ*<sup>2</sup> → −0 in Figure 4c.

**Figure 4.** (**a**) Stationary weakly asymmetric (with respect to the two components) solutions of Equation (7), obtained for *µ*<sup>2</sup> = −0.4 with fixed *N*<sup>1</sup> = 100. Dashed and solid curves display density profiles of the first (*n*<sup>1</sup> ) and second (*n*2) components, respectively. (**b**) The semi-log plot of the density profiles of *n*<sup>2</sup> for *µ*<sup>2</sup> = −0.4, −0.04, and 0 at *z* > 0. (**c**) Dependences of *N*<sup>2</sup> (black dots: the left vertical axis) and asymmetry parameter *δ*21, defined as per Equation (28) (the red dashed line pertaining to the right vertical axis), on *µ*<sup>2</sup> for fixed *N*<sup>1</sup> = 100. The parameters are *P* = 1 (*g*<sup>1</sup> = *g*2) and *δg*/*g* = 0.05. The symmetric point with *N*<sup>1</sup> = *N*<sup>2</sup> = 100 and *δ*<sup>21</sup> = 0 corresponds to *µ*<sup>1</sup> = *µ*<sup>2</sup> = −0.88878.

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)}. \tag{28}$$

It increases almost linearly with *µ*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* 6= 1 (*g*<sup>1</sup> 6= *g*2)

Next, we consider the QDs for *P* 6= 1, setting *P* > 1 without loss of generality. Then, the MF energy is minimized for *n*<sup>2</sup> > *n*1; the situation with *n*<sup>1</sup> > *n*<sup>2</sup> 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 *δg*/*g* = 0.05, *N*<sup>1</sup> = 100, and several different values of *P*, varying *µ*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*<sup>2</sup> and parameter *δ*<sup>21</sup> (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 E1D − *µ*1*ρ*<sup>1</sup> − *µ*2*ρ*<sup>2</sup> [5,14], which leads to the zero-pressure condition,

$$p(\rho\_1, \rho\_2) = -\mathcal{E}\_{\rm ID} + \sum\_{j=1,2} \left(\frac{\partial \mathcal{E}\_{\rm ID}}{\partial \rho\_j}\right) \rho\_j$$

$$\equiv -\mathcal{E}\_{\rm ID} + \mu\_1 \rho\_1 + \mu\_2 \rho\_2 = 0. \tag{29}$$

*Symmetry* **2020**, *12*, 174

From this, we obtain relation

$$\frac{\left(\sqrt{\mathcal{G}\_1}\rho\_1 - \sqrt{\mathcal{G}\_2}\rho\_2\right)^2}{2} + \frac{\mathcal{g}\delta\mathcal{g}\left(\sqrt{\mathcal{G}\_2}\rho\_1 + \sqrt{\mathcal{G}\_1}\rho\_2\right)^2}{(\mathcal{g}\_1 + \mathcal{g}\_2)^2} - \frac{\sqrt{m}}{3\pi\hbar}(\mathcal{g}\_1\rho\_1 + \mathcal{g}\_2\rho\_2)^{3/2} = 0,\tag{30}$$

which can be rewritten in the scaled form as

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

For given *n*1, we solved Equation (31) to find the respective value of *n*2, which is shown in Figure 6 for *δ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<sup>j</sup>* 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*<sup>2</sup> = *Pn*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*<sup>2</sup> for larger *P* (note that the figure displays a log–log plot).

**Figure 5.** (**a**) Stationary solutions of Equation (7), obtained for *δg*/*g* = 0.05 and *N*<sup>1</sup> = 100. From the left panel to the right one, the parameter in Equation (8) is *P* = 1.25, 1.67, and 2.5, and the chemical potential for the second component is *µ*<sup>2</sup> = −0.018, −0.011, and −0.006, respectively, just below the threshold above which the tails of *ψ*<sup>2</sup> extend to infinity. Dashed and solid curves represent the density of the first (*n*<sup>1</sup> ) and second (*n*2) components. (**b**) Dependences of *N*<sup>2</sup> (black dots: the left vertical axis) and asymmetry parameter *δ*21, defined as per Equation (28) (the red dashed line pertaining to the right vertical axis), on *µ*<sup>2</sup> for fixed *N*<sup>1</sup> = 100 and *P* = 1.25 or *P* = 1.67.

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\_{\rm T} \tag{32}
$$

where *N*<sup>T</sup> = *N*<sup>1</sup> + *N*<sup>2</sup> is the total norm, which admits the existence of the QDs. For given *N*1, we obtain the upper bound for *N*<sup>2</sup> above which the solution becomes delocalized, and calculate the corresponding critical value of ∆21. The results are shown in Figure 7. For the system with *P* = 1 , the curve demonstrates an empirical dependence ∆<sup>21</sup> ∝ *N* −*α* <sup>T</sup> with exponent *α* ≈ 0.58. Accordingly, the asymmetry tends to vanish asymptotically for very "heavy" droplets, at *N*<sup>T</sup> → ∞. As the system becomes slightly asymmetric, with *P* = 1.25, exponent *α* is significantly reduced for small *N*T, and converges to a certain finite value at *N*<sup>T</sup> → ∞. 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.

**Figure 6.** The negative-pressure region in the (*n*<sup>1</sup> , *n*2) plane for *δg*/*g* = 0.05 and values of asymmetry parameter in Equation ( 8) *P* = 1 (the solid curve), 1.25 (dashed), 2.5 (dashed-dotted), and 10 (dotted). Boundaries are determined by the zero-pressure condition, as given by Equation (31). The negative pressure, at which localized states may exist, occurs inside the boundaries. Thin lines represent relation *n*<sup>2</sup> = *Pn*<sup>1</sup> .

**Figure 7.** The inverse of the largest relative norm difference ∆21, up to which the asymmetric droplets exist (see Equation (32)), shown as a function of the total number, *N*T, at different values of asymmetry parameter of Equation (8). Here we set *δg*/*g* = 0.05.

#### 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 *ψj*(*z*, *t*) = √*n<sup>j</sup> e* −*iµ<sup>j</sup> t* , (*j* = 1, 2). The substitution of this in Equation (7) yields

$$\mu\_1 \quad = (P + GP^{-1})n\_1 + (-1 + G)n\_2 - \frac{P}{\pi} \sqrt{P n\_1 + \frac{n\_2}{P}},$$

$$\mu\_2 \quad = (P^{-1} + GP)n\_2 + (-1 + G)n\_1 - \frac{1}{\pi P} \sqrt{P n\_1 + \frac{n\_2}{P}}.\tag{33}$$

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

$$m\_{\overline{l}} = \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\_{\dot{\jmath}}(z,t) = \left[\sqrt{n\_{\dot{\jmath}}} + \delta\psi\_{\dot{\jmath}}(z,t)\right]e^{-i\mu\_{\dot{\jmath}}t},\tag{34}
$$

$$\delta\psi\_{\rangle} = \mathbb{Z}\_{\rangle}\cos(kz - \Omega t) + i\eta\_{\rangle}\sin(kz - \Omega t),\tag{35}$$

with infinitesimal amplitudes *ζ<sup>j</sup>* and *η<sup>j</sup>* , cf. Equation (24). The substitution of this in Equation (7) and the linearization with respect to *ζ<sup>j</sup>* and *η<sup>j</sup>* yields the dispersion equation for the perturbation:

$$
\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}} \tag{36}
$$

where

$$\begin{aligned} P\_1 &= (P + GP^{-1})n\_1, \quad P\_2 = (P^{-1} + GP)n\_2, \\ Q\_1 &= \frac{P^2 n\_1}{2\pi\sqrt{P n\_1 + P^{-1} n\_2}}, \quad Q\_2 = \frac{P^{-2} n\_2}{2\pi\sqrt{P n\_1 + P^{-1} n\_2}}.\end{aligned} \tag{37}$$

$$\mathcal{R} = (-1 + \mathcal{G})\sqrt{n\_1 n\_2}, \quad \mathcal{S} = \frac{\sqrt{n\_1 n\_2}}{2\pi\sqrt{P n\_1 + P^{-1} n\_2}}.$$

For *P* = 1 and *n*<sup>1</sup> = *n*2, these results reproduce Equation (25) for the Ω<sup>−</sup> branch. A parameter region in which at least one squared eigenfrequency Ω<sup>2</sup> <sup>±</sup> 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 *σ* = Im(Ω) for the asymmetric PWs in the symmetric system with *P* = 1 and *δg*/*g* = 0.05, in the plane of wavenumber *k* and density ratio *n*<sup>12</sup> = *n*2/*n*2. For the consistency with the single-component situation displayed in Figure 3, we here fix the total density as (*n*<sup>1</sup> + *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/ √ 2. The largest gain is obtained at equal densities, *n*<sup>12</sup> = 1. Both the *k*-band of the instability and magnitude of the gain slowly decrease as the deviation of *n*<sup>12</sup> 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)).

**Figure 8.** Color-coded values of the MI gain, *σ* = Im(Ω), for asymmetric PWs, as calculated from Equation (36) in the plane of wave number |*k*| and density ratio *n*<sup>12</sup> = *n*1/*n*2, are displayed for (**a**) *P* = 1 and (**b**) *P* = 1.25 with fixed *δg*/*g* = 0.05 and (*n*<sup>1</sup> + *n*2)/2 = 10. The solid and dashed white curves represent the MI boundary *k* = *k*<sup>0</sup> and the peak value of the MI gain at *k* = *k*max = *k*0/ √ 2, respectively. In (**c**), we plot *σ*(*k*max) (circles) and *n* max <sup>12</sup> (triangles) versus *P*.

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*<sup>12</sup> = *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*<sup>12</sup> = 1 and Figure 9f,g for *n*<sup>12</sup> = 9. This observation can be understood in terms of the MI gain *σ*, as shown in Figure 8c, where *σ* 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*<sup>12</sup> is explained by the decrease of *k*max, see Figure 8a. For *n*<sup>12</sup> = 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*<sup>12</sup> 6= 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* 6= 1

Figure 8b represents the MI gain for *P* = 1.25 and a fixed total density, (*n*<sup>1</sup> + *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*<sup>12</sup> = *n* max <sup>12</sup> = 0.577, below the equal-densities point *n*<sup>12</sup> = 1. This is consistent with the fact that, at *P* > 1, unequal values *n*<sup>1</sup> < *n*<sup>2</sup> are suitable to the formation of an asymmetric soliton structure, as seen in Figure 5a. In Figure 8c, we plot the peak MI gain, *σ*(*k*max), along with the respective value of the density ratio, *n* max <sup>12</sup> , as a function of *P*. Value *n* max <sup>12</sup> 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*<sup>1</sup> + *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*<sup>12</sup> = *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* 6= 1. Spatial profiles at fixed time, displayed in Figure 10b–d, show fragmentation of the profiles. Figure 10c clearly indicates that,

even for *n*<sup>12</sup> = 1, the MI generates asymmetric droplet-like structures similar to Figure 5a, where the complete overlapping of the two densities does not occur.

**Figure 9.** Numerically simulated development of the MI of asymmetric PW states in the two-component system, with *P* = 1 and *δg*/*g* = 0.05 . The initial PW states are taken with fixed total density, (*n*<sup>1</sup> + *n*2)/2 = 10. (**a**) The evolution of the central density of the first component, *n*1(*z* = 0), for different density ratios in the two components, *n*<sup>12</sup> = *n*1/*n*2. (**b**,**c**) Snapshots of density profiles for the cases of (**b**) *n*<sup>12</sup> ≡ *n*1/*n*<sup>2</sup> = 1 at *t* = 200 and (**c**) *n*<sup>12</sup> = 9 at *t* = 400. Panels (**d**,**e**) and (**f**,**g**) are top views of the spatiotemporal evolution of the densities, *n*<sup>1</sup> (*z*, *t*) and *n*2(*z*, *t*), for *n*<sup>12</sup> = 1 and *n*<sup>12</sup> = 9, respectively. Simulations were performed in the domain −50 ≤ *z* ≤ +50 with 2048 grid points, subject to periodic boundary conditions. In this figure and in Figure 10, the scaled time unit corresponds to ∼ 1 µs in physical units.

**Figure 10.** Numerically simulated development of the modulational instability in the two-component system with *δg*/*g* = 0.05 and *P* = 1.25. The initial PW states are taken with a fixed total density, (*n*<sup>1</sup> + *n*2)/2 = 10. (**a**) The evolution of the central density of the first component, *n*1(*z* = 0), for different density ratios in the two components, *n*<sup>12</sup> = *n*1/*n*2. (**b**–**d**) Snapshots of density profiles for the cases of (**b**) *n*<sup>12</sup> ≡ *n*1/*n*<sup>2</sup> ∼ 0.1 at *t* = 300, (**c**) *n*<sup>12</sup> = 1 at *t* = 200 and (**d**) *n*<sup>12</sup> = 9 at *t* = 600. Panels (**e**–**g**) represent the top view of the spatiotemporal evolution of the densities, *n*<sup>1</sup> (*z*, *t*), corresponding to (**b**–**d**), respectively (the evolution of *n*<sup>2</sup> (*z*, *t*) shows similar patterns). Simulations were performed in the domain −50 ≤ *z* ≤ +50 with 2048 grid points, subject to periodic boundary conditions.

#### **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 µs. This time is much smaller than a typical lifetime of the droplet, which is &100 µs [21–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.

**Author Contributions:** Conceptualization and investigation, T.M., A.M., K.K., B.A.M., and A.K.; analytics, A.K. and B.A.M.; numerics, T.M., K.K., and A.M.; writing—original draft preparation, T.M.; writing—review and editing, A.M., K.K., B.A.M., and A.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** We appreciate valuable comments received from M. Modugno. T.M. acknowledges support from IBS (Project Code IBS-R024-D1). A.M. acknowledges support from the Ministry of Education, Science, and Technological Development of Republic of Serbia (project III45010) and the COST Action CA 16221. The work of K.K. is partly supported by the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for Scientific Research (KAKENHI Grant No. 18K03472). B.A.M. appreciates support from the Israel Science Foundation through grant No. 1287/17. A.K. thanks the Indian National Science Academy for the grant of INSA Scientist Position at Physics Department, Savitribai Phule Pune University).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Other Exact Solutions for the Single-Component GP Equation**

Here we briefly list other types of exact solutions of the single-component Equation (12), in addition to the FT and PW solutions in Equations (13) and (22) which were considered in detail above (solutions to Equation (12) in the form of dark and anti-dark solitons were reported in [46]).The stability of a majority of these solutions is not addressed here, as it should be a subject for a separate work.

#### *Appendix A.1. δg*/*g* > 0

In the case of comparable quadratic self-attraction and cubic repulsion in Equation (12) with *δg* > 0, exact spatially-periodic solutions with odd parity can be expressed in terms of the Jacobi's elliptic sine, whose modulus *q* is an intrinsic parameter of the family:

$$\psi(z,t) = \exp\left(-i\mu\_{\text{sn}}t\right) \left[A\,\text{sn}(\beta z, q) + B\right],\tag{A1}$$

where

$$B = \frac{\sqrt{2}}{3\pi} \frac{\mathcal{g}}{\delta g} > 0, \quad A = \sqrt{\frac{2}{1+q^2}} B > 0, \quad \mu\_{\text{sn}} = -2\frac{\delta g}{\mathcal{g}} \mathcal{B}^2 < 0, \quad \mathcal{\beta}^2 = \frac{2}{(1+q^2)} \frac{\delta g}{\mathcal{g}} \mathcal{B}^2.$$

In the limit of *q* → 1, Equation (A1) goes over into the kink (the same as found in [46]),

$$\psi(z,t) = \exp\left(-i\mu\_{\text{kink}}t\right) \left[A\tanh(\beta z) + B\right],\tag{A2}$$

with parameters

$$A = B = \frac{\sqrt{2}}{3\pi} \frac{\text{g}}{\delta g} > 0, \quad \mu\_{\text{kink}} = -2 \frac{\delta \text{g}}{\mathcal{S}} B^2, \quad \mathcal{S}^2 = \frac{\delta \text{g}}{\mathcal{S}} B^2.$$

#### *Appendix A.2. δg*/*g* < 0

In the case when the inter-species MF attraction is stronger than the intra-species repulsion, resulting in *δg* < 0, spatially-periodic solutions are expressed in terms of even Jacobi's elliptic functions, dn(*x*, *q*) and cn(*x*, *q*). First, it is

$$\psi(z,t) = \exp\left(-i\mu\_{\mathrm{dn}}t\right) \left[A\,\mathrm{dn}(\beta z, q) + B\right],\tag{A3}$$

with the elliptic modulus taking all values 0 < *q* < 1, other parameters being

$$B = \frac{\sqrt{2}}{3\pi} \frac{\mathcal{g}}{\delta g} < 0, \quad A = -\sqrt{\frac{2}{2 - q^2}} B > 0, \quad \mu\_{\text{dn}} = -2B^2 \frac{\delta g}{\mathcal{g}} > 0, \quad \mathcal{J}^2 = -\frac{2}{(2 - q^2)} \frac{\delta g}{\mathcal{g}} B^2.$$

The second solution is expressed in terms of the elliptic cosine, with *q* <sup>2</sup> > 1/2:

$$\psi(z,t) = \exp\left(-i\mu\_{\rm cn}t\right) \left[A \operatorname{cn}(\beta z, q) + B\right],\tag{A4}$$

$$B = \frac{\sqrt{2}}{3\pi} \frac{\mathcal{g}}{\delta g} < 0, \quad A = -\sqrt{\frac{2}{2q^2 - 1}} B > 0, \quad \mu\_{cn} = -2\frac{\delta g}{\mathcal{g}} B^2 > 0, \quad \mathcal{J}^2 = -\frac{2}{(2q^2 - 1)} \frac{\delta g}{\mathcal{g}} B^2.$$

In the limit of *q* → 1, both solutions in Equations (A3) and (A4) carry over into a state of the "bubble" type [49], which changes the sign at two points (the same solution was reported as an "W-shaped soliton" in [46]):

$$\psi(z,t) = \exp\left(-i\mu\_{\text{bulk}}t\right) \left[A\text{sech}(\beta z) + B\right],\tag{A5}$$

with parameters

$$B = \frac{\sqrt{2}}{3\pi} \frac{\text{g}}{\delta \text{g}} < 0, \quad A = -\sqrt{2}B > 0, \quad \mu \text{ (bubble } = \beta^2 = -2\frac{\delta \text{g}}{\mathcal{S}}B^2 > 0.1)$$

#### **Appendix B. Analytical Solutions for Strongly Asymmetric Fundamental and Dipole States**

Here we consider analytical solutions of Equation (7) with strong asymmetry, *N*<sup>1</sup> *N*2, which can be found under small-amplitude conditions, *n*1(*z* = 0) *n*2(*z* = 0) *n*0. Then, cubic terms may be neglected in Equation (7), and approximation p *P*|*ψ*1| <sup>2</sup> + *P*−<sup>1</sup> |*ψ*2| <sup>2</sup> ≈ *<sup>P</sup>* <sup>−</sup>1/2 <sup>|</sup>*ψ*2<sup>|</sup> is used to simplify Equation (7) to the following equations for stationary states in Equation (9):

$$
\mu\_1 \phi\_1 = -\frac{1}{2} \frac{d^2 \phi\_1}{dz^2} - \frac{\sqrt{P}}{\pi} \phi\_2 \phi\_1 \tag{A6}
$$

$$
\mu\_2 \phi\_2 = -\frac{1}{2} \frac{d^2 \phi\_2}{dz^2} - \frac{1}{\pi P^{3/2}} \phi\_2^2. \tag{A7}
$$

Although this case is somewhat formal, in terms of the underlying concept of the quantum droplets, which is essentially based on the competition of residual MF and LHY terms, it is interesting to consider it too.

The soliton solution of Equation (A7) is obvious,

$$\phi\_2(z) = \frac{3\pi}{2} \left( -\mu\_2 \right) \frac{P^{3/2}}{\cosh^2 \left( \sqrt{-\mu\_2/2} z \right)} \tag{A8}$$

where the solution in Equation (13) takes essentially the same form in the limit of |*µ*| *µ*0. Then, the substitution of Equation (A8) in Equation (A6) makes it tantamount to the linear Schrödinger

equation with the Pöschl-Teller potential [50]. The ground-state (GS) solution of Equation (A6) for *φ*1, with arbitrary amplitude *φ* (0) 1 ,

$$\left(\phi\_1(z)\right)\_{\rm GS} = \frac{\phi\_1^{(0)}}{\left[\cosh\left(\sqrt{-\mu\_2/2}z\right)\right]^{\gamma}}\tag{A9}$$

exists with

$$\gamma = \frac{1}{2} \left( \sqrt{24P^2 + 1} - 1 \right),\tag{A10}$$

and eigenvalue

$$(\mu\_1)\_{\rm GS} = \left(\sqrt{24P^2 + 1} - 1\right)^2 \frac{\mu\_2}{16}.\tag{A11}$$

In this case, the QD solutions are quasi-Gaussian objects [20]. Note that, in the symmetric system with *P* = 1, Equations (A10) and (A11) yield *γ* = 2 and (*µ*1)GS = *µ*2, i.e., the eigenmode and eigenvalue coincide with their counterparts in the soliton solution in Equation (A8), while they are different in the asymmetric system, the GS level lying below or above the chemical potential of soliton in Equation (A8) at *g*<sup>1</sup> > *g*<sup>2</sup> and *g*<sup>1</sup> < *g*2, respectively.

In Figure A1 we compare a typical asymptotic solution given by Equations (A8) and (A9) with a numerically obtained GS solution for the same values of the parameters. It is seen that the analytical and numerical results match well.

**Figure A1.** Comparison of the asymptotic analytical solutions, given by Equations (A8) and (A9), with their numerically obtained counterparts. The density of the first (*n*<sup>1</sup> ) and second (*n*2) components are displayed in top and bottom panels, respectively. Solid blue lines represent the numerical results, while dashed red lines represent the analytical solution. Here, parameters are *δg*/*g* = 0.05, *N*<sup>1</sup> = 0.0001067, *N*<sup>2</sup> = 0.0148044 and (*µ*2)GS = *µ*<sup>2</sup> = −0.005.

Further, it is also possible to produce the first excited state of Equation (A6) in the form of the dipole (antisymmetric) mode with an arbitrary amplitude:

$$\Phi\left(\phi\_1(z)\right)\_{\text{dip}} = \frac{\phi\_1^{(0)} \sinh\left(\sqrt{-\mu\_2/2}z\right)}{\left[\cosh\left(\sqrt{-\mu\_2/2}z\right)\right]^\gamma},\tag{A12}$$

where *γ* is the same as in Equation (A10), the respective eigenvalue being

$$\left(\left(\mu\_1\right)\_{\text{dip}} = \left(\sqrt{24P^2 + 1} - 3\right)^2 \frac{\mu\_2}{16}\right) \tag{A13}$$

which is obviously higher than its GS counterpart in Equation (A11) (at *P* = 1, Equation (A13) yields (*µ*1)dip = *µ*2/4, and (*µ*1)dip falls below *µ*<sup>2</sup> for *P* > √ 2). Unlike the GS, the dipole mode exists not at all values of *P*, but only for *P* > √ 1/3. Exactly at *P* = √ 1/3, one has (*µ*1)dip = 0, and the dipole mode in Equation (A12), with *<sup>γ</sup>* <sup>=</sup> 1, is a delocalized one, <sup>∼</sup> tanh p −*µ*2/2*z* .

Linear Schrödinger Equation (A6) with the Pöschl-Teller potential may give rise to higher bound states of integer order *ν* as well, with eigenvalues

$$(\mu\_1)\_\nu = \left(\sqrt{24P^2 + 1} - (1 + 2\nu)\right)^2 \frac{\mu\_2}{16} \tag{A14}$$

where *ν* = 0 and 1 correspond to Equations (A11) and (A13), respectively, the *ν*-th spate existing at *P* <sup>2</sup> > *ν* (*ν* + 1) /6. The number of such solutions is always finite.

Unlike solutions considered in Appendices A and C, the stability of the solutions given by Equations (A8)–(A14) is obvious.

#### **Appendix C. Other Exact Solutions in the Case of** *N*<sup>1</sup> *N*<sup>2</sup>

Here we provide periodic solutions to the semi-linear system of Equations (A6) and (A7) in terms of Jacobi elliptic functions. In the limit of *q* → 1, they go over into solutions given in the main text, in the form of Equations (A8), (A9), and (A12).

#### *Appendix C.1. Solution of Equation (A7)*

An exact periodic solution of Equation (A7) with the quadratic nonlinearity is

$$\phi\_2 = A \left[ \text{dn}^2(\beta z, q) + p \right], \tag{A15}$$

with

$$\beta^2 = -\frac{\mu\_2}{2\sqrt{1-q+q^2}}, \quad A = -\frac{3\pi\mu\_2 P^{3/2}}{2\sqrt{1-q+q^2}}, \quad p = \frac{-(2-q) + \sqrt{1-q+q^2}}{3}.\tag{A16}$$

In the limit of *q* → 1, the solution in Equation (A15) goes over into the solution in Equation (A8). Note that *p* is vanishing in this limit, according to Equation (A16).

#### *Appendix C.2. Solutions of Equation (A6)*

We now show that, with *φ*<sup>2</sup> given by Equation (A15), linear Equation (A6) *φ*<sup>1</sup> has several particular solutions depending on the value of *P*.

**Solutions For** *P* <sup>2</sup> = 1/3

Appendix C.2.1. Solution I

It is easy to check that

$$
\phi\_1 = \phi\_1^{(0)} \text{dn}(\beta z, q) \tag{A17}
$$

is an exact solution to Equation (A6), provided that

$$P^2 = \frac{1}{3}, \quad \mu\_1 = \left(\frac{\mu\_2}{12}\right) \frac{2 - q + 2\sqrt{1 - q + q^2}}{\sqrt{1 - q + q^2}}.$$

Appendix C.2.2. Solution II

$$\boldsymbol{\phi}\_{1} = \boldsymbol{\phi}\_{1}^{(0)} \mathbf{c} \mathbf{n}(\boldsymbol{\beta} \boldsymbol{z}, \boldsymbol{q}) \tag{A18}$$

is an exact solution to Equation (A6), provided that

$$P^2 = \frac{1}{3}, \quad \mu\_1 = \left(\frac{\mu\_2}{12}\right) \frac{2q - 1 + 2\sqrt{1 - q + q^2}}{\sqrt{1 - q + q^2}}.$$

In the limit of *q* → 1, solutions I and II go over into the solution Equation (A9) with *γ* = 1 and *µ*<sup>1</sup> = *µ*2/4.

Appendix C.2.3. Solution III

$$\phi\_1 = \phi\_1^{(0)} \text{sn}(\beta z, q) \tag{A19}$$

is an exact solution to Equation (A6), provided that

*P* <sup>2</sup> = 1 3 , *µ*<sup>1</sup> = *µ*2 12 2 p 1 − *q* + *q* <sup>2</sup> − (<sup>1</sup> + *<sup>q</sup>*) p 1 − *q* + *q* 2 .

In the limit of *q* → 1, solution III goes over into the solution Equation (A12) with *γ* = 1 and *µ*<sup>1</sup> = 0. **Solutions for** *P* <sup>2</sup> = 1**.**

Appendix C.2.4. Solution IV

It is easy to check that

$$\phi\_1 = \phi\_1^{(0)} [\text{dn}^2(\beta z, q) + p] \tag{A20}$$

is an exact solution to Equation (A6), provided that

$$P^2 = 1, \quad \mu\_1 = \mu\_2 \dots$$

Appendix C.2.5. Solution V

$$\phi\_1 = \phi\_1^{(0)} \text{cn}(\beta z, q) \text{dn}(\beta z, q) \tag{A21}$$

is an exact solution to Equation (A6), provided that

$$P^2 = 1, \quad \mu\_1 = \left(\frac{\mu\_2}{2}\right) \frac{q + \sqrt{1 - q + q^2}}{\sqrt{1 - q + q^2}}.$$

In the limit *q* = 1, solutions IV and V go over into solution Equation (A9) with *γ* = 2 and *µ*<sup>1</sup> = *µ*2. Appendix C.2.6. Solution VI

$$\phi\_1 = \phi\_1^{(0)} \text{sn}(\beta z, q) \text{dn}(\beta z, q) \tag{A22}$$

is an exact solution to Equation (A6), provided that

$$P^2 = 1, \quad \mu\_1 = \left(\frac{\mu\_2}{4}\right) \frac{\Im(1-q) + \sqrt{1-q+q^2}}{\sqrt{1-q+q^2}}.$$

Appendix C.2.7. Solution VII

$$\boldsymbol{\phi}\_{1} = \boldsymbol{\phi}\_{1}^{(0)} \operatorname{sn}(\beta z, \boldsymbol{\eta}) \operatorname{cn}(\beta z, \boldsymbol{\eta}) \tag{A23}$$

is an exact solution to Equation (A6), provided that

$$P^2 = 1, \quad \mu\_1 = \left(\frac{\mu\_2}{4}\right) \frac{2\sqrt{1-q+q^2} - (2-q)}{\sqrt{1-q+q^2}} \frac{1}{4}$$

In the limit of *q* → 1, solutions VI and VII go over into Equation (A12), with *γ* = 2 and *µ*<sup>1</sup> = *µ*2/4.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
