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

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.


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]).

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]: where m is the atomic mass (the same for both components), ρ j = |Ψ j | 2 (j = 1, 2) is the density of the j-th component, represented by the MF wave function Ψ j , and g ≡ √ g 1 g 2 , δg ≡ g 12 + g.
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][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 ∼ −ρ 3/2 (for ρ 1 = ρ 2 ≡ ρ) 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 ∼ ρ 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, +∞ −∞ E 1D dZ, gives rise to the system of GP equations, which include the LHY correction, where T and Z are the time and coordinate measured in physical units, and parameter 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: Further, rescaling mg 2 casts Equation (3) in the normalized form, where parameter determines the asymmetry of the system, in the case of P = 1. Note that, as concerns stationary solutions with chemical potentials µ 1,2 , sought for as states with mutually proportional components, φ 1 (z) = Kφ 2 (z), are only possible in the fully symmetric case with P = 1, µ 1 = µ 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 ψ 1 = ψ 2 .

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.

The Single-Component GP Model
Under the single-component reduction of the binary system, with g 1 = g 2 ≡ g and ψ 1 = ψ 2 ≡ ψ, Equation (1) simplifies to [5] ε 1D ≡h with the single dimensionless density, n = |ψ| 2 ≡ h 2 /mg ρ. Assuming a spatially uniform state, the equilibrium density and the corresponding chemical potential are given by Density n 0 corresponds to the minimum of the energy per particle, ∂ n n −1 ε 1D (n) = 0, and µ 0 is negative for δg/g > 0. The corresponding single GP equation is with normalization condition +∞ −∞ |ψ(z)| 2 dz = N, where N ≡ N 1 = N 2 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,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.

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: This solution exists in a finite range of negative values of the chemical potential µ 0 < µ < 0, featuring the FT shape at 0 [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 µ 0 towards µ = 0, the maximum density of the localized mode, monotonously decreases from n 0 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/ √ −µ 0 at µ/µ 0 ≈ 0.776. Further increase of µ towards µ = 0 makes the QD broader, its width  (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 0 = 36.025 and µ 0 = −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 , n c ). 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 µ = µ 0 + 0.00001, very close to the delocalization limit (the transition to PW).
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, µ 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 µ < 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 ∼ sech 2 −µ/2z and the nonlinear-Schrödinger one, ∼ sech −2µz . For δg/g < 0, the N(µ) dependence for the soliton family carries over into the following form, which is an analytical continuation of Equation (15). This dependence also satisfies the VK criterion.

The Plane-Wave Solution
The PW solution of Equation (12 )can be presented in a form ψ(z, t) = √ n exp (iK PW z − iµt) with wavenumber K PW and constant density n, which determine the corresponding chemical potential: The Galilean invariance of Equation (12) implies that any quiescent solution ψ 0 (z, t) generates a family of moving ones, with arbitrary velocity c. Therefore, K PW may be canceled by means of transformation For given µ, Equation (19) produces two different branches of the density as a function of µ (here, 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 µ: µ c = −(2π 2 δg/g) −1 = (9/8)µ 0 , the respective density being Values µ = µ c and n = n c correspond to the spinodal point [5], and n + (µ 0 ) = n 0 (see Equation (13)). Note that the above-mentioned existence region of the soliton solution in terms of the chemical potential, µ 0 < µ < 0, lies completely inside that of the PW state, which is µ c ≤ µ. Thus, the soliton always coexists with the PW (this fact is also obvious in Figure 1).

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 The substitution of this expression in Equation (12) and linearization with respect to perturbation δψ leads to the corresponding Bogoliubov-de Gennes equation, By looking for perturbation eigenmodes with wavenumber k and frequency Ω, and real infinitesimal amplitudes ζ and η, Equation (23) yields a dispersion relation for the eigenfrequencies: 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 ] −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 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 with k 0 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].  (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 n c , 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][22][23]48]. Another option, which is specific to the 1D setting, is sudden decrease of density n by relaxing the transverse trapping.

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 = 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 = g 2 (i.e., P = 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 φ 1 = 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 1 = 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, +∞ −∞ dz|ψ 1 (z)| 2 = N 1 , while chemical potential µ 2 is fixed in the other one, allowing its norm N 2 to vary. Figure 4 displays essential features of weakly asymmetric droplets for δg/g = 0.05 and fixed N 1 = 100. The symmetric (completely overlapping) solution with N 1 = N 2 is found at µ 1 = µ 2 = −0.88878. When µ 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 µ 2 , but Figure 4b demonstrates that, at µ 2 → −0, ψ 2 develops small-amplitude extended tails, which are absent in ψ 1 . Due to the contribution of the tails, the approach of µ 2 < 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 µ 2 → −0 is opposite to the decay of the QD's norm in the single-component model at µ → −0, cf. Equation (15). At µ 2 ≥ 0, the ψ 2 component undergoes delocalization, with its tails developing a nonzero background at |z| → ∞, as seen in the density profile displayed in Figure 4b for µ 2 = 0, and norm N 2 (µ 2 ) diverging at µ 2 → −0 in Figure 4c.   (28) (the red dashed line pertaining to the right vertical axis), on µ 2 for fixed N 1 = 100. The parameters are P = 1 (g 1 = g 2 ) and δg/g = 0.05. The symmetric point with N 1 = N 2 = 100 and δ 21 = 0 corresponds to µ 1 = µ 2 = −0.88878.
In Figure 4c, we also plot the parameter of the asymmetry between the two components, defined as 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.

Asymmetric
QDs in the System with P = 1 (g 1 = g 2 ) Next, we consider the QDs for P = 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 δg/g = 0.05, N 1 = 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 2 and parameter δ 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 E 1D − µ 1 ρ 1 − µ 2 ρ 2 [5,14], which leads to the zero-pressure condition, From this, we obtain relation which can be rewritten in the scaled form as 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 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 = 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 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, where N 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 ∆ 21 . The results are shown in Figure 7. For the system with P = 1 , the curve demonstrates an empirical dependence ∆ 21 ∝ N −α T with exponent α ≈ 0.58. Accordingly, the asymmetry tends to vanish asymptotically for very "heavy" droplets, at N T → ∞. 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 T → ∞. 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.

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 j e −iµ j t , (j = 1, 2). The substitution of this in Equation (7) yields Accordingly, in the symmetric system with P = 1, densities of the asymmetric PW state are expressed in terms of the chemical potentials as We introduce the perturbation around the PW states as with infinitesimal amplitudes ζ j and η j , cf. Equation (24). The substitution of this in Equation (7) and the linearization with respect to ζ j and η j yields the dispersion equation for the perturbation: where P 1 = (P + GP −1 )n 1 , P 2 = (P −1 + GP)n 2 , For P = 1 and n 1 = n 2 , these results reproduce Equation (25) for the Ω − branch. A parameter region in which at least one squared eigenfrequency Ω 2 ± is negative gives rise to the MI of the two-component state.

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 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 / √ 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)).  (36) in the plane of wave number |k| and density ratio n 12 = n 1 /n 2 , are displayed for (a) P = 1 and (b) P = 1.25 with fixed δg/g = 0.05 and (n 1 + n 2 )/2 = 10. The solid and dashed white curves represent the MI boundary k = k 0 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 12 (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 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 σ, 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 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 = 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 = 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 max 12 = 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, σ(k max ), along with the respective value of the density ratio, n max 12 , as a function of P. Value n max 12 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 = 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. (e) (f) (g) 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 1 + 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 12 = n 1 /n 2 . (b,c) Snapshots of density profiles for the cases of (b) n 12 ≡ n 1 /n 2 = 1 at t = 200 and (c) n 12 = 9 at t = 400. Panels (d,e) and (f,g) are top views of the spatiotemporal evolution of the densities, n 1 (z, t) and n 2 (z, t), for n 12 = 1 and n 12 = 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.

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][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. 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 with the elliptic modulus taking all values 0 < q < 1, other parameters being The second solution is expressed in terms of the elliptic cosine, with q 2 > 1/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]): with parameters

Appendix B. Analytical Solutions for Strongly Asymmetric Fundamental and Dipole States
Here we consider analytical solutions of Equation (7) with strong asymmetry, N 1 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|ψ 1 | 2 + P −1 |ψ 2 | 2 ≈ P −1/2 |ψ 2 | is used to simplify Equation (7) to the following equations for stationary states in Equation (9): 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, 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 φ exists with and eigenvalue 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 1 > g 2 and g 1 < 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. 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: where γ is the same as in Equation (A10), the respective eigenvalue being 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 µ 2 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 γ = 1, is a delocalized one, ∼ tanh −µ 2 /2z . 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 where ν = 0 and 1 correspond to Equations (A11) and (A13), respectively, the ν-th spate existing at P 2 > ν (ν + 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.
In the limit of q → 1, solutions I and II go over into the solution Equation (A9) with γ = 1 and µ 1 = µ 2 /4.

Appendix C.2.3. Solution III
is an exact solution to Equation (A6), provided that In the limit of q → 1, solution III goes over into the solution Equation (A12) with γ = 1 and µ 1 = 0. Solutions for P 2 = 1.