↩︎

1 Introduction↩︎

Since the first experimental realization of the Bose-Einstein condensate [1]–[3], outstanding progress has been sealed in terms of realizing, monitoring and controlling multicomponent atomic gases [4]–[6]. This is largely owed to exploiting Feshbach resonances [7] and trapping techniques [8], [9] to manipulate the interatomic interaction strength or range and the external confinement landscape respectively. Starting from single-component bosonic gases, a multitude of correlated phases beyond the weakly interacting mean-field (MF) limit were observed, such as the (Super-) Tonks-Girardeau states for strong (attractions [10], [11]) repulsions [12]–[15]. Turning to bosonic mixtures arguably enriched correlated phases, already for weak interactions, occur due to the competition between the intra- and inter-species interactions. This has lead, for instance, to phase separation [16]–[18], when the interspecies repulsion surpasses the average intraspecies one, bubble phases at the immiscibility threshold [19], or the recently realized quantum droplets [20]–[23] when the interspecies attraction slightly overcomes the intraspecies repulsions.

The emergence of the self-bound liquid-type configurations known as quantum droplets in contact interacting three-dimensional (3D) bosonic mixtures [24], [25], as well as in dipolar gases [26], [27], is a manifestation of the impact of quantum fluctuations. These may be captured by the Lee-Huang-Yang (LHY) energy [20], [28], representing the first-order correction to the Bogoliubov MF theory, which arrests the expected collapse of the ensuing attractive system at the edge of the MF stability regime [5]. Several properties of these intriguing many-body bound states of matter have already been explored. These refer, exemplarily, to their spectrum [29], [30], collective excitations [20], [29], [31], collisions [31], [32] and coexistence with nonlinear soliton [33], [34] and vortex [35]–[38] states, within the context of the LHY theory [21], [22], [39], see also Refs. [40]–[43] for the exposition of beyond-LHY effects by deploying ab-initio methods. Experimentally, quantum droplets have been observed in 3D in both homonuclear [24], [25], [44] and heteronuclear [45]–[47] mixtures. One of the main challenges faced by these experimental efforts stems from the relatively short droplet lifetimes caused (primarily) by three-body recombination [48]. These lifetimes are increased in relevant 1D settings [49], [50] due to lower densities but also in heteronuclear mixtures [45].

A particularly interesting property is the dependence of the LHY term on the dimensionality [50]–[52], in contrast to MF interactions [6]. In 1D that we focus herein, the LHY term is attractive [53] instead of being repulsive as in 3D, significantly changing the droplet parametric region and its phenomenology. Concretely, the 1D geometry causes the droplet parametric regions to i) not feature collapse, ii) host stable bright-soliton solutions [54]–[56] and iii) coincide with the respective MF stability regimes. Hence, 1D platforms share the premise of hosting, comparatively easier to control, enriched droplet phases even more in the genuine two- [57]–[61] and three-component [62] settings. Additionally, quantum droplets and more generally the impact of quantum fluctuations in 1D can be studied throughout the MF stability regime. However, the majority of the current investigations have been motivated by the correspondence to the 3D case, and thus solely focused on the attractive threshold of the MF stability regime. Hence, a general framework being able to describe droplets across the MF stability region but also heteronuclear mixtures in 1D is still lacking. Additionally, extraction and investigation of the LHY energy in higher-component setups is an open question. Indeed, previous studies on three-component mixtures in 1D have been focused on the case of a single-impurity in the third component [63]–[65] or a MF density coupling between the droplet and the impurities [66], [67]. Instead, appropriate inclusion of the LHY term among all components may be able to unveil unseen self-bound structures.

To address the above-discussed open problems, we derive the LHY energy (in a closed form whenever possible) and associated extended Gross-Pitaevskii equations (eGPEs) for two-, three- and four-component 1D bosonic mixtures based on the perturbative Bogoliubov treatment. Previous treatments restricted to two-component settings [53], can be obtained as a limiting case of our framework and are restricted to the 3D wave collapse threshold. In contrast to those, which (for convenience) we dub herein as “original", our exact approach remains valid throughout the 1D MF stability regime. Namely, it encompasses both attractive and repulsive interspecies interactions and interestingly it is able to tackle homonuclear as well as heteronuclear two-component settings. Direct comparisons of our approach to the approximate one in the attractive intercomponent interaction regime reveals quantitative deviations. These are imprinted as a reduced droplet saturation density and increased width in the exact case. Turning to the repulsive side, our results predict an early onset of phase-separation for both homonuclear and heteronuclear mixtures which is absent within the approximate framework.

Next, we provide exact closed form expressions for the LHY energies and ensuing eGPEs of a three-component mixture with either three or two identical components. It is explicated that such settings offer a broad range of possibilities for creating mixed droplet phases that are absent in lower-component setups. Additionally, a closed form LHY expression is extracted for the fully symmetric four-component mixture, constituting the higher-component mixture for which the Bogoliubov modes may be derived analytically in the general case. We present the derivations in a detailed and structured manner aiming to provide a useful resource for researchers entering the rapidly expanding field of quantum droplets and in general attractively interacting bosonic setups.

This work is structured as follows. In Section 2, we introduce the general multicomponent bosonic mixture featuring short-range contact interactions, the Bogoliubov transformation, and the original approximate two-component LHY theory being valid at the edge of the MF stability regime. Section 3 is devoted to the derivation and solution of the exact eGPEs for the genuine two-component mixture across the entire MF stability regime. This treatment is generalized to the three-component (Sec. 4) and the four-component (Sec. 5) bosonic mixture. Concluding remarks and perspectives for future investigations are offered in Sec. 6. In Appendix 7, we provide further details regarding the three-component mixture stability conditions. Appendix 8 discusses the existence of asymmetric (domain-wall like) droplet distributions in the phase-separated regime.

2 Basic settings and background↩︎

2.1 Multicomponent bosonic droplet setting↩︎

We consider a multicomponent bosonic mixture consisting of \(N_{\sigma}\) atoms of mass \(m_{\sigma}\) (where \(\sigma = A , B , C, D, \dots\)) residing in 1D free space and featuring periodic boundary conditions. Notice, however, that our results (to be presented below) persist also in the case of a box potential provided that its length is sufficiently large and does not impact the edges of the discussed self-bound droplet configurations. The mixture is assumed to be close to zero temperature, where \(s\)-wave scattering dominates [68]. Hence, interparticle interactions are represented by contact potentials. Their effective coupling strengths are intracomponent repulsive (i.e. \(g_{\sigma\sigma}\equiv g_{\sigma}>0\)) and intercomponent attractive (\(g_{\sigma\sigma^\prime}<0\)) or repulsive (\(g_{\sigma\sigma^\prime}>0\)). In a corresponding experiment, they can be adjusted either via Fano-Feshbach resonances [69], [70] through the 3D scattering length or through confinement induced resonances [68] by modifying the transverse trap frequency, \(\omega_{\perp}\). This is significantly tight such that transverse excitations are prohibited and the atomic motion is restricted across the elongated \(x\)-direction as in corresponding 1D experiments [71].

The many-body Hamiltonian in second quantized form reads \[\begin{align} \label{MB95Hamilt} \begin{aligned} H = \int dx \Bigg( \sum_{\sigma} \Big(& \Psi_{\sigma}^\dagger \hat{h}_{\sigma} \Psi_{\sigma} + \frac{g_{\sigma}}{2} \Psi_{\sigma}^\dagger \Psi_{\sigma}^\dagger \Psi_{\sigma}\Psi_{\sigma} + \sum_{\sigma\prime > \sigma} \frac{g_{\sigma\sigma^\prime}}{2} \Psi_{\sigma}^\dagger \Psi_{\sigma^\prime}^\dagger \Psi_{\sigma}\Psi_{\sigma^\prime} \Big) \Bigg), \end{aligned} \end{align}\tag{1}\] where \(\hat{h}_{\sigma}=-\frac{\hslash^2}{2m_\sigma} \left(\frac{\partial^2}{\partial {x_\sigma}^2}\right)\) is the underlying single-particle hamiltonian and \(\Psi_\sigma(x)\) refers to the field operator annihilating a \(\sigma\) species boson at position \(x\). In what follows, for computational convenience, we rescale the above Hamiltonian with respect to the energy scale set by \(m_Ag_A^2/\hbar^2\). Hence, the length, time, and interaction strengths are given in units of \(\hbar^2/(m_Ag_A)\), \(\hbar^3/(m_Ag_A^2)\), and \(g_A\) respectively.

2.2 Bogoliubov Transformation↩︎

For an arbitrary set of \(n\) bosonic operators \(A^\dagger = (\hat{\alpha}_1^\dagger,\;\hat{\alpha}_2^\dagger,\; \dots,\; \hat{\alpha}_n^\dagger\;, \hat{\alpha}_1,\;\hat{\alpha}_2,\;,\dots,\;\hat{\alpha}_n )\) the general Bogoliubov transformation is defined by \(A=U B\). Here, \(B^\dagger = (\hat{\beta}_1^\dagger,\;\hat{\beta}_2^\dagger,\; \dots,\; \hat{\beta}_n^\dagger\;, \hat{\beta}_1,\;\hat{\beta}_2,\;,\dots,\;\hat{\beta}_n )\) are the Bogoliubov operators and \(U = \begin{pmatrix} M & N \\ N^* & M^* \end{pmatrix}\) is the Bogoliubov transformation matrix, where \(M,N\) are \(n \times n\) matrices [72]. Demanding that both \(A\) and \(B\) satisfy bosonic commutation relations translates to \(\mathcal{M}_b = [A,A^\dagger] = [B,B^\dagger] = \begin{pmatrix} I_{n} & 0 \\ 0 & -I_{n} \end{pmatrix}\), with \(I_{n}\) denoting the \(n\)-dimensional identity matrix. As such, utilizing the arbitrary Bogoliubov transformation yields \[\label{BosonicMetric} \begin{align} \mathcal{M}_b = U \mathcal{M}_b U^\dagger. \end{align}\tag{2}\] Thus, demanding that Eq. 2 holds is equivalent to imposing that the Bogoliubov transformation preserves the commutation relations of the set of bosonic operators. However, in this case, it is not required to explicitly derive the constraints on the elements of the transformation matrix \(U\), as done in the usual textbook derivation [5], [6].

Another useful property of the Bogoliubov transformation is that it can be used to diagonalize a bilinear (in terms of the arbitrary set of bosonic operators \(A\)) operator, such as \(\hat{H} = A^\dagger H A\). Here, \(\hat{H}\) is an one-body operator (expressed in second quantized form) with respect to \(\hat{a}_{i}\), \(\hat{a}_{i}^\dagger\) and \(H\) is its corresponding matrix representation. Accordingly, it holds that \(\hat{H} = A^\dagger H A = B^\dagger U^\dagger H UB = B^\dagger \mathcal{M}_b U^{-1} \mathcal{M}_b H U B\). Then, by requiring that the Bogoliubov transformation¹ diagonalizes \(\hat{H}\) in terms of the bilinear form with respect to \(B\) or equivalently diagonalizes the matrix \(\mathcal{H}_M = \mathcal{M}_b H\), as \(U^{-1} \mathcal{H}_M U= D\), where \(D\) is a diagonal matrix, we obtain the eigenvalues of \(\mathcal{H}_M\) through the characteristic equation \[\label{CharacteristicEquation} \det{(\mathcal{M}_b H - \lambda I_{2n})}=0.\tag{3}\] This leads to \(\mathcal{M}_b D = {\rm diag}(\lambda_1, ..., \lambda_n, -\lambda_{n+1}, ..., -\lambda_{2n} )\) and therefore the operator \(\hat{H}\) takes the form \[\label{DiagonalizedOperator} \begin{align} \hat{H} &= B^\dagger \mathcal{M}_bD B = \sum_{i=1}^{n} \big( \lambda_{i} \hat{\beta}_{i}^\dagger \hat{\beta}_{i} - \lambda_{n+i} \hat{\beta}_{i} \hat{\beta}_{i}^\dagger \big) = \sum_{i=1}^{n} (\lambda_{i}-\lambda_{n+i}) \hat{\beta}_{i}^\dagger \hat{\beta}_{i} - \sum_{i=1}^{n} \lambda_{n+i}. \end{align}\tag{4}\] In the last step, we have used that \([\hat{b}_{i}, \hat{b}_{i}^\dagger] = 1\). In this sense, \(h_0 = - \sum_{i=1}^{n} \lambda_{n+i}\) is the eigenvalue of \(\hat{H}\) corresponding to the quasi-particle vacuum in terms of the bosonic quasi-particle operators \((\hat{b}_{i}, \;\hat{b}_{i}^\dagger\)), while \(h_i = \lambda_{i}-\lambda_{n+i}\) are the eigenvalues corresponding to the quasi-particle excitations.

2.3 Original two-component LHY treatment↩︎

For a two-component, homonuclear \((m_A=m_B \equiv m)\), 1D mixture the MF stability regime is given by the condition \(g_{AB}^2 \leq g_{A}g_{B}\) [5]. In the attractive limit where \(g_{AB} = - \sqrt{g_{A}g_{B}}\) holds, taking into account the first-order (LHY) quantum correction it was shown [53] to result in the following coupled set of “approximate" eGPEs describing two-component 1D quantum droplets \[\label{ApproxeGPE} \begin{align} i\hbar\frac{\partial \Psi_\sigma}{\partial t} = \Bigg(& - \frac{\hbar^2}{2m} \frac{\partial^2 }{\partial x_\sigma^2} + g_{\sigma}\left| \Psi_{\sigma} \right|^2 + g_{AB}\left| \Psi_{\sigma^\prime\neq\sigma} \right|^2 - \frac{g_{\sigma}\sqrt{m}}{\pi\hslash} \sqrt{g_A|\Psi_A|^2 +g_B|\Psi_B|^2 } \Bigg) \Psi_{\sigma}. \end{align}\tag{5}\] In these equations, the first terms represent the usual mean-field Gross-Pitaevskii equations (GPEs) and the last term provides the LHY correction. In the symmetric system, where the intracomponent interactions \(g_\sigma \equiv g\) and the species densities \(n_\sigma \equiv n\), with \(N_\sigma = \int n_\sigma dx =\\ \int \left| \Psi_\sigma \right|^2 dx\) the components behave equivalently. Accordingly, it has been demonstrated [21], [23] that these eGPEs predict a transition of the droplet density from a Gaussian to a flat-top (FT) configuration upon either increasing the particle number (\(N\)) or decreasing the intercomponent attraction (\(|g_{AB}|\)) with \(0< - g_{AB} < g\) [53]. The emergent FT structures feature a saturated peak density at [53] \(n_{0}=\frac{8mg^3}{9\hbar^2\pi^2(g+g_{AB})^2}\), with energy density \(\mathcal{E}_{0}(n=n_0) = -(g+g_{AB}) n_{0}^2\), chemical potential density \(\mu_0/L = \frac{\partial \mathcal{E}_{0}}{\partial n}(n=n_0) = -(g+g_{AB})n_0\), and healing length \(\xi_0 = \\ \sqrt{ - \frac{\hbar^2}{2m\mu_0/L} } = \sqrt{\frac{\hbar^2}{2m(g+g_{AB})n_0}}\).

It is imperative to clarify that the extraction of Eq. (5 ) is based on the assumption that we are operating at the interaction boundary \(g_{AB} = - \sqrt{g_{A}g_{B}}\). This is indeed well motivated by the respective 3D setting in which the LHY contribution is repulsive and provides stabilization outside the MF stability regime where droplets emerge [20]. However, in 1D systems the LHY correction turns out to be attractive [23], [53], and therefore the droplet configurations are hosted within the MF stability regime. Hence, as we will explicate in detail below, no assumptions (beyond the usual ones deployed for a dilute and weakly interacting gas) are needed, and we can analytically extract eGPEs which are valid throughout the MF stability regime and not only in the boundary as was currently known. This is one of the main results of our work which presents the generalized two-component eGPEs and exposes deviations from the predictions of the “original" ones.

3 Two-component mixture↩︎

The Hamiltonian density for a short-range (contact) interacting two-component bosonic mixture in a box of length \(L\), can be easily deduced from the general multicomponent Hamiltonian given by Eq. (1 ) with \(\sigma=A, B\) and \(\sigma'=A,B \neq \sigma\). A corresponding experimentally relevant system of choice here could be, for instance, two different hyperfine states of \(^{39}\)K as in the 3D case [24]. Next, we express the field operators in the plane-wave basis as \(\hat{\Psi}_A(x) = \frac{1}{\sqrt{L}} \hat{a}_0 + \sum_{k\neq 0} \frac{e^{ikx}}{\sqrt{L}} \hat{a}_k\) and \(\hat{\Psi}_B(x) = \frac{1}{\sqrt{L}} \hat{b}_0 + \sum_{k\neq 0} \frac{e^{ikx}}{\sqrt{L}} \hat{b}_k\), where we have implicitly assumed periodic boundary conditions. Substitution of these operators into the aforementioned two-component Hamiltonian and integration over the box length, leads to the Hamiltonian \[\label{Hamilt2Comp} \begin{align} \hat{\mathcal{H}}& = \frac{g_{A}}{2L} \hat{a}_0^\dagger \hat{a}_0^\dagger \hat{a}_0 \hat{a}_0 + \frac{g_{B}}{2L} \hat{b}_0^\dagger \hat{b}_0^\dagger \hat{b}_0 \hat{b}_0 + \frac{g_{AB}}{L} \hat{a}_0^\dagger \hat{b}_0^\dagger \hat{a}_0 \hat{b}_0 + \frac{\hbar^2}{2m_A} \sum_{k\neq 0} k^2 \hat{a}_k^\dagger \hat{a}_k + \frac{\hbar^2}{2m_B} \sum_{k\neq 0} k^2 \hat{b}_k^\dagger \hat{b}_k \\ & + \frac{g_{A}}{2L} \sum_{k\neq 0} \big[ \hat{a}_0^\dagger \hat{a}_0^\dagger \hat{a}_k \hat{a}_{-k} + \hat{a}_0 \hat{a}_0 \hat{a}_k^\dagger \hat{a}_{-k}^\dagger + 4\hat{a}_0^\dagger \hat{a}_0 \hat{a}_k^\dagger \hat{a}_{k} \big] \\ & + \frac{g_{B}}{2L} \sum_{k\neq 0} \big[ \hat{b}_0^\dagger \hat{b}_0^\dagger \hat{b}_k \hat{b}_{-k} + \hat{b}_0 \hat{b}_0 \hat{b}_k^\dagger \hat{b}_{-k}^\dagger + 4\hat{b}_0^\dagger \hat{b}_0 \hat{b}_k^\dagger \hat{b}_{k} \big] \\ & +\frac{g_{AB}}{L} \sum_{k\neq 0} \big[ \hat{a}_0^\dagger \hat{b}_0^\dagger \hat{a}_k \hat{b}_{-k} + \hat{a}_0^\dagger \hat{a}_0 \hat{b}_k^\dagger \hat{b}_{k} + \hat{a}_0^\dagger \hat{b}_0 \hat{b}_k^\dagger \hat{a}_{k} + \hat{b}_0^\dagger \hat{a}_0 \hat{a}_k^\dagger \hat{b}_{k} + \hat{b}_0^\dagger \hat{b}_0 \hat{a}_k^\dagger \hat{a}_{k} + \hat{a}_0\hat{b}_0 \hat{a}_k^\dagger \hat{b}_{-k}^\dagger \big] \\ & + \mathcal{O}(\hat{a}_k^3). \end{align}\tag{6}\] Here, terms higher than bilinear in the bosonic operators \(\hat{a}_k\), \(\hat{b}_k\) are neglected. From this expression, it is possible to readily retrieve the zeroth order (MF) approximation by setting \(\hat{a}_0=\sqrt{N_A}\), \(\hat{b}_0=\sqrt{N_B}\) and neglecting all other terms with \((\hat{a}_{k\neq0},\hat{b}_{k\neq0})\). Such a process results in the known MF energy density of the two-component setting \[\label{MFenergy2Components} \mathcal{E}_{\rm MF} = \frac{E_{\rm MF}}{L} =\frac{1}{2} g_{A}n_A^2 + \frac{1}{2} g_{B}n_B^2 + g_{AB}n_An_B .\tag{7}\] To obtain the first-order correction to the MF energy, one needs to account for the depletion of the condensate due to quantum fluctuations. This is achieved by expressing the particle number operators as \(N_A = \hat{a}_0^\dagger \hat{a}_0 + \sum_{k\neq 0} \hat{a}_k^\dagger \hat{a}_{k}\) and \(N_B = \hat{b}_0^\dagger \hat{b}_0 + \sum_{k\neq 0} \hat{b}_k^\dagger \hat{b}_{k}\). Inserting those into the Hamiltonian of Eq. 6 and defining the set of operators \(\Phi ^\dagger=(\hat{a}_k^\dagger,\hat{b}_k^\dagger,\hat{a}_{-k},\hat{b}_{-k})\), we arrive at the bilinear Hamiltonian \[\label{bilinearTwoComp} \begin{align} \hat{H}=& \sum_{k>0} \Phi^\dagger H \Phi + \frac{g_{A}N_A^2}{2L} + \frac{g_{B}N_B^2}{2L}+ \frac{g_{AB}N_B N_A}{L} - \sum_{k>0} \Big( \frac{\hbar^2 k^2}{2m_A} +\frac{\hbar^2 k^2}{2m_B} + g_{A}n_A + g_{B}n_B \Big). \end{align}\tag{8}\] The operators, \(\Phi\), satisfy the bosonic commutation relations, i.e. \(\mathcal{M}_b = [\Phi,\Phi^\dagger] = \begin{pmatrix} I_{2} & 0 \\ 0 & -I_{2} \end{pmatrix}\), where \(I_{2}\) is the \(2 \times 2\) identity matrix. Also, the Hamiltonian matrix \[\label{HamiltMatrixTwoComp} H = \begin{pmatrix} h_A(k) & h_{AB} & g_{A}n_A & h_{AB} \\ h_{AB} & h_B(k) & h_{AB} & g_{B}n_B\\ g_{A}n_A & h_{AB} & h_A(k) & h_{AB} \\ h_{AB} & g_{B}n_B & h_{AB} & h_B(k) \end{pmatrix},\tag{9}\] where \(h_\sigma(k) = \frac{\hbar^2 k^2}{2m_\sigma} + g_{\sigma}n_\sigma\) and \(h_{AB}= g_{AB}\sqrt{n_A n_B}\). Note that the sums in Eq. 8 are now over positive momenta only. As argued in Sec. 2.2, we can diagonalize the Hamiltonian of Eq. 8 in terms of an appropriate Bogoliubov transformation, by solving the characteristic equation (Eq. 3 ) \(\det{(\mathcal{M}_b H - \omega_k I_{4})}=0\) which in this case takes the form \(\omega^4 + \beta \omega^2 + \gamma=0\). The parameters are \(\beta = - \epsilon_A^2(k) - \epsilon_B^2(k)\) and \(\gamma = \epsilon_A^2(k) \epsilon_B^2(k) - g_{AB}^2n_An_B \frac{\hbar^4k^4}{m_Am_B}\), with \(\epsilon_\sigma(k) = \\\sqrt{\frac{\hbar^4k^4}{4m_\sigma^2} + \frac{\hbar^2k^2}{m_\sigma}g_{\sigma}n_{\sigma} }\) representing the corresponding single-component Bogoliubov dispersion relations [6]. The solutions of the characteristic polynomial equation yield the underlying dispersion relation for the two-component system in the presence of the first-order quantum correction \[\label{SolutionsTwoComponent} \begin{align} \omega^2_{\pm } = &\frac{\epsilon_A^2 + \epsilon_B^2 \pm \sqrt{(\epsilon_A^2 - \epsilon_B^2)^2 + 4 g_{AB}^2n_An_B \frac{\hbar^4k^4}{m_Am_B}}}{2}. \end{align}\tag{10}\]

Finally, utilizing the diagonal form of a bilinear operator within the Bogoliubov transformation prescribed by Eq. 4 and collecting the constant terms in Eq. 8 , we find the respective ground state energy \[\label{LHYTwoComponent} \begin{align} &E_0 =\frac{g_{A}N_A^2}{2L} + \frac{g_{B}N_B^2}{2L}+ \frac{g_{AB}N_B N_A}{L} + \sum_{k>0} \Big( \omega_{+} + \omega_{-} - \frac{\hbar^2 k^2}{2m_A} -\frac{\hbar^2 k^2}{2m_B} - g_{A}n_A - g_{B}n_B \Big). \end{align}\tag{11}\] Apparently, the first terms in Eq. 11 corresponds to the standard MF energy, while the sum is the first-order LHY quantum correction for the two-component mixture. To calculate the LHY contribution, we substitute the sum in Eq. 11 with an integral over \(k\) according to \(\sum_{k>0} \xrightarrow{} (L/2\pi) \int_0^{+\infty} dk\). We note in passing that in the 3D case, additional terms (\(\propto +\frac{mg_{\sigma\sigma^\prime} n_\sigma}{\hbar^2 k^2}\)) enter the corresponding sum over \(k\), while the latter has to be substituted with the corresponding 3D integral [6].

3.1 Exact closed form solution for the homonuclear mixture↩︎

As already discussed above, unlike the 3D case (see Ref. [20]), in 1D it is not necessary to impose the approximation \(g_{AB}^2\approx g_{A}g_{B}\) (which leads to Eq. 5 ) in the dispersion relation of Eq. 10 . In particular, it is possible to calculate exactly the ground state energy for a homonuclear (\(m_A=m_B \equiv m\)) mixture. Indeed, setting \(p:=\frac{4(g_{AB}^2-g_{A}g_{B})n_An_B}{(g_{A}n_A+g_{B}n_B)^2}\), and by employing the change of variables \(\hbar k = \sqrt{m (g_{A}n_A + g_{B}n_B)} y\) to calculate the integral in Eq. 11 , the ground state energy density becomes \[\label{LHYGSEnergyMassBalance2Comp} \mathcal{E} = \frac{E_0}{L} = \mathcal{E}_{\rm MF} + \underbrace{\frac{\sqrt{m}(g_{A}n_A+g_{B}n_B)^\frac{3}{2}}{2\hbar\pi} \mathcal{I}(p)}_{\mathcal{E}_{\rm LHY}},\tag{12}\] where the function of the beyond MF energy contribution \[\label{Ifactor} \mathcal{I}(p) = - \frac{\sqrt{2}}{3} \Bigg( (1-\sqrt{p+1})^\frac{3}{2} + (1+\sqrt{p+1})^\frac{3}{2} \Bigg ).\tag{13}\] Importantly, the generalized two-component energy density described by Eq. 12 is valid for any interaction strength and particle number for which the MF treatment is valid. Moreover, the LHY energy is real-valued and negative, as long as \(-1\leq p\leq 0\), i.e. \(\left| g_{AB} \right| \leq \sqrt{g_{A}g_{B}}\), so that \(p\leq 0\) (while \(p \geq -1\) is trivially satisfied). Notice that this interaction interval corresponds to the entire MF stability region, ranging from the usual droplet threshold \(( -g_{AB} \approx \sqrt{g_{A}g_{B}})\) all the way to the standard miscibility threshold \((g_{AB} \approx \sqrt{g_{A}g_{B}})\) [16], [23]. It is also worth mentioning explicitly that, in the limit of \(p=0\), we obtain \(\mathcal{I}(p=0)=-4/3\), i.e. we recover the original droplet energy which results in Eq. 5 . Moreover, at the edge of the MF stability regime where \(g_{AB}^2 \approx g_{A}g_{B}\), it holds that \(\left| p \right| \ll 1\), leading to \(I(p)\approx -\frac{4}{3} - \frac{p}{2}\). Finally, in the opposite limit of an uncoupled mixture (\(g_{AB}=0\)) we find \(I(p=-1)=-2\sqrt{2}/3\). Hence, in this case the LHY term is reduced by a factor of \(\sqrt{2}\). A similar expression was recently derived in the context of a mobile impurity immersed in a droplet host [66].

Having determined the ground state energy of the two-component bosonic mixture in the presence of the LHY correction, we can subsequently derive the exact eGPEs. This is done by evaluating the Euler-Lagrange equations with respect to \(n_{\sigma}\) yielding: \[\label{eGPETwoComponent} \begin{align} &i\hbar\frac{\partial \Psi_\sigma}{\partial t} = (- \frac{\hbar^2}{2m} \frac{\partial^2 }{\partial x^2} + \frac{\partial \mathcal{E}}{\partial n_\sigma})\Psi_\sigma = \Bigg(- \frac{\hbar^2}{2m} \frac{\partial^2 }{\partial x^2} + g_{\sigma}\left| \Psi_{\sigma} \right|^2 + g_{AB}\left| \Psi_{\sigma^\prime\neq\sigma} \right|^2 + \frac{\partial \mathcal{E}_{\rm LHY}}{\partial n_\sigma}\Bigg)\Psi_\sigma . \end{align}\tag{14}\] In general, the full eGPEs have a quite complicated form, since

\[\label{LHYterm2Comp} \begin{align} \frac{\partial \mathcal{E}_{\rm LHY}}{\partial n_{\sigma}} &= \frac{3g_{\sigma}\sqrt{m(g_{A}n_A + g_{B}n_B)}}{4\hbar\pi}\mathcal{I}(p) +\frac{\sqrt{m(g_{A}n_A + g_{B}n_B)^3}}{2\hbar\pi}\frac{\partial \mathcal{I}}{\partial p} \frac{\partial p}{\partial n_{\sigma}}, \end{align}\tag{15}\] where \[\begin{align} \frac{\partial \mathcal{I}}{\partial p} = \frac{\sqrt{2} ( \sqrt{1-\sqrt{p+1}} - \sqrt{1 +\sqrt{p+1}}) }{4\sqrt{p+1}},\tag{16}\\ \frac{\partial p}{\partial n_{\sigma}} = \frac{4 n_{\sigma\prime} (g_{A}g_{B} - g_{AB}^2)( g_{\sigma}n_{\sigma} - g_{\sigma\prime}n_{\sigma\prime} )}{(g_{A}n_A + g_{B}n_B)^3}, \tag{17} \end{align}\] with \(\sigma \neq \sigma\prime\). Notably, the form and subsequent impact of the LHY correction away from the droplet limit \(( - g_{AB} \approx \sqrt{g_{A}g_{B}})\) and even for repulsive intercomponent interactions (where e.g. bubble configurations have been predicted [19]) is largely unknown. Considering stronger intercomponent interactions outside the MF stability, i.e. \(\left| g_{AB} \right|>\sqrt{g_{A}g_{B}}\) where \(p>0\), results in a complex LHY correction. A common practice in the literature is to simply neglect the imaginary part of \(\mathcal{E}\) [20], [62], assuming that it is sufficiently small and hence the growth rate of the resulting instability is slower than e.g. the decay due to three-body recombination. Examining the validity of this treatment by comparing its predictions e.g. with MF computations, ab-initio calculations [40]–[43] or experimental data would be a particularly interesting direction for future studies, especially in 1D. This is because, aside from being more amenable to ab-initio calculations, 1D systems do not experience MF collapse in this limit, unlike their 3D counterparts where stable solutions exist only in the presence of the LHY term [20].

An insightful application of our generalized framework is to consider the limiting case of a fully balanced mixture characterized by \(g_{A}=g_{B} \equiv g\) and \(n_A=n_B \equiv n\), and therefore \(p = (\frac{g_{AB}}{g})^2 -1\). Accordingly, the second term in Eq. 15 vanishes and the LHY term becomes \(\sim (g^3n)^{1/2}I(p)\). Thus, in this case, the exact LHY term has the same density dependence as in the “approximate" eGPEs [53] given by Eq. 5 , but they are quantitatively different since the exact LHY term depends also on \(g_{AB}\). Using the exact energy of the full two-component system given by Eq. 12 , we find the minimum of the energy per particle for the aforementioned balanced mixture to occur at the saturation density \[\label{SaturationDensity} \begin{align} n_s &= \frac{mg^3I^2(p)}{2\hbar^2\pi^2(g+g_{AB})^2} = \frac{2mg^3 \big[ 1+3(g_{AB}/g)^2 +(1-(g_{AB}/g)^2)^{3/2} \big]}{9\hbar^2\pi^2(g+g_{AB})^2}. \end{align}\tag{18}\] Here, the energy density \(\mathcal{E}(n=n_s) = \frac{-m^2g^6I^4(p)}{8\hbar^4\pi^4(g+g_{AB})^3} = - (g+g_{AB}) n_s^2\), the chemical potential density \(\mu_s/L = \frac{-mg^3I^2(p)}{2\hbar^2\pi^2(g+g_{AB})}= - (g+g_{AB}) n_s\), and the healing length \(\xi_s = \sqrt{ - \frac{\hbar^2}{2m\mu_s/L} } = \sqrt{\frac{\hbar^2}{2m(g+g_{AB})n_s}}\).

Apparently, these exact results exhibit quantitative differences from the predictions of the original LHY approach described by Eq. (5 ). This can be readily deduced, for instance, by comparing them with the approximate saturation density (\(n_0\)) as well as the chemical potential and healing length discussed in Sec. 2.3 for the balanced mixture. Note that the chemical potentials and energies for all 1D droplet configurations discussed throughout this work, are negative. This supports their self-bound nature in the entire MF stability regime.

To further elucidate the deviations between the original (Eq. 5 ) and the exact (Eq. 14 ) eGPEs predictions we numerically compute, via imaginary time-propagation, the corresponding ground state droplet densities for the symmetric mixture with \(g_{A}=g_{B} \equiv g\) and \(N_A=N_B\equiv \\ N=10\). Since, the exact eGPEs (Eq. 14 ) are valid throughout the MF stability regime, we present results for the standard droplet regime \(g_{AB}<0\) [Fig. 1 (a), (b)], but also for the nearly intercomponent decoupled (\(g_{AB}\approx 0\)) system [Fig. 1 (c)] and even for repulsive interactions within the miscible regime quantified by \(0<g_{AB} \leq g\) [Fig. 1 (d)]. The underlying LHY energy, \(\mathcal{E}(n=n_s) = - (g+g_{AB}) n_s^2\), attains its maximum (or minimum absolute) value for the decoupled system, whilst it decreases as \(g_{AB}\) decreases or increases. Interestingly, it acquires a relatively small imaginary component for \(\left| g_{AB} \right|> g\) as shown in Fig. 1 (e).

A close inspection of Fig. 1 reveals that the original eGPEs (Eq. 5 ) are only accurate when \(\left| g_{AB} \right| = g\). In contrast, deviating from this condition they significantly overestimate the droplet’s saturation density and localization in all other interaction regimes, see Fig. 1 (a)-(d). In particular, the difference in the saturation density among the two approaches becomes larger for intermediate intercomponent couplings lying in the interval \(-g<g_{AB}<g\) as illustrated in Fig. 1 (f). Even more, for \(g_{AB}= - 0.6\) the original eGPEs predict a Gaussian profile for the ground state density and they fail to capture the FT density profile of the droplet as predicted by the exact eGPEs, see Fig. 1 (b). Additionally, as mentioned before, it is clear that the droplet solution exists throughout the MF stability regime within both approaches, exhibiting FT configurations (for sufficiently large particle number) irrespectively of the intercomponent interaction, \(g_{AB}\), value [Fig. 1 (a)-(d)]. Specifically, the droplet becomes progressively more delocalized and its saturation density decreases as \(g_{AB}\) varies from \(g_{AB}=-g\) to \(g_{AB} = g\). Namely, the droplet saturation density (width) reduces (increases) sharply between the two limits, reaching a minimum value of \(n_s(g_{AB}=g) =\frac{2mg}{9\hbar^2\pi^2}\approx 0.0225\) for the parameter values considered here, see also Fig. 1 (f). As such, it could be that in practice, after a certain value of \(g_{AB}\) the droplet would be dilute and extended enough appearing to have effectively disperse rendering challenging its experimental resolution by the corresponding imaging apparatus.

As a next step, we turn our attention to an intracomponent interaction imbalanced homonuclear (\(m_A=m_B \equiv m\)) mixture, with \(g_{A}\neq g_{B}\) and fixed \(g_B=0.5g_{A}\), as well as \(N_A=N_B \equiv N=10\). The corresponding two-component droplet configurations within the exact [Eq. 14 ] and original [Eq. 5 ] eGPEs, but also the solutions obtained by calculating the LHY energy through direct numerical integration of the sum in Eq. 11 are visualized in Fig. 2 upon different variations of the \(g_{AB}\) attraction². For clarity, the ensuing solutions are marked by “Exact",”Original", and “Num" respectively. In general, it appears that as long as the intercomponent coupling is attractive (i.e. \(g_{AB}<0\)) the two components prefer to maximize their spatial overlap [30], [61], while exhibiting similar density profiles as can be seen in Fig. 2 (a). Moreover, tuning \(g_{AB}\) towards weaker attractions, both components become less localized, until a transition from a Gaussian-type to a FT density profile takes place, compare in particular Fig. 2 (a) and (b). Once more, as in the fully symmetric mixture, the original eGPEs predict appreciably more localized droplet densities alongside a higher FT value.

Additionally, approaching the weakly coupled intercomponent regime (\(g_{AB}\approx0\)) the individual component droplet separation becomes more prominent, with the more strongly repulsive component exhibiting a more localized density profile and a higher saturation density [Fig. 2 (b)]. Strikingly, this effect is absent within the original eGPEs, predicting instead a largely overlapping behavior for all intercomponent interactions shown in Fig. 2. Turning to repulsive intercomponent couplings, e.g. \(g_{AB} = 0.4g_A\), we observe a phase separation of the droplet density between the two components, even though we are still operating well within the MF stability regime, see Fig. 2 (c). A behavior that is completely absent within the original eGPEs. This phase separation can be understood from energetic arguments. Indeed, by ignoring the spatial overlap between the two components, it is possible to approximate the system shown in Fig. 2 (c) by three distinct single-component droplets. Estimating their energy, in the homogeneous FT case, through minimization of the energy per particle (see also Eq. 12 ) for each component in the absence of the other, results in \(E_\sigma^{\rm 1comp}=-\frac{8m_\sigma^3g_\sigma^3}{81 \hbar^4 \pi^4}L_\sigma\). Here, \(L_\sigma\) designates the spatial extent of each individual droplet. With this approximation (setting \(L_\sigma\) as the full-width-at-half-maximum for each structure) we find the total energy \(E =E^{\rm1comp}_A + 2E^{\rm1comp}_B \approx-0.2832\) which is in excellent agreement with the numerically obtained energy \(E=-0.2816\). The relatively small deviations stem from the (positive) contributions of the finite overlap and the kinetic terms. In particular, we can approximate the kinetic energy of each droplet as \(K.E._\sigma = \int \frac{\hbar^2}{2m_\sigma} |\frac{\partial \Psi_\sigma}{\partial x}|^2 dx \approx 2 \frac{\hbar^2}{2m_\sigma} |\frac{ \sqrt{n_s}_\sigma}{ \xi_\sigma}|^2 \xi_\sigma = \frac{8m_\sigma g_\sigma^2}{27\hbar^2 \pi^2} = \frac{2g_\sigma}{3\pi}n_s\), where \(\xi_\sigma\) is the healing length of the droplet and the factor of two stems from the two edges of the droplet. Interestingly, this suggests that the kinetic energy of the droplet scales as \(\propto m_\sigma g_\sigma^2\) explaining why the more weakly interacting (and as we shall see below the lighter component) splits preferentially in the phase separated configurations, in order to reduce the systems kinetic energy. Clearly, an overlapping arrangement, with component B in a single FT configuration centered at \(x=0\) and length \(L_B^\prime\approx 2L_B\) would be associated with higher energy due to the repulsive intercomponent coupling. The energy contribution of the latter is \(E_{AB}\approx g_{AB} n_A n_B L_{overlap}\), and it would be of order \(\sim +0.1\) assuming similar saturation densities to those in Fig. 2 (c) and \(L_{overlap} \approx L_A\leq L_B^\prime\). Hence, this hypothetical overlapping configuration would have a total energy of order \(E^\prime \approx E^{\rm1comp}_A + E^{\rm1comp}_B(L_\sigma= L_B^\prime) + E_{AB} = E + E_{AB} \sim -0.18\), which is significantly larger that the phase separated configuration.

Note that throughout this work, we focus on parity symmetric³ density configurations, an assumption that pertains only to phase-separated structures. The droplet equations, however, support also parity asymmetric (or domain-wall like) configurations which, in the case of phase-separation, possess slightly lower energies (typically of the order of \(0.3 \%\)) than the parity symmetric splited configurations e.g. depicted in Fig. 2 (c). Such parity asymmetric distributions while existing within the classical field approximation are not, in general, permitted within a many-body wavefunction perspective which in the absence of an explicitly symmetry breaking contribution in the Hamiltonian should retain the symmetry. In this context, the lowest energy parity symmetric configurations correspond to a good approximation to the reduced one-body density of the underlying atomic gas in the ground state, see also Appendix 8 for further details and the respective asymmetric density distributions. It is important to emphasize that these asymmetric single domain-wall like configurations are intriguing on their own right and deserve further investigation, while they are promising for being exploited for unprecedented nucleation of topological excitations and hydrodynamic instabilities in the droplet realm. Finally, it is worth noting that in the presence of e.g. an external harmonic confinement, the components phase-separate at the standard miscibility threshold [23], [73].

3.2 Heteronuclear two-component mixtures↩︎

Besides homonuclear mixtures, droplets can also be hosted in genuinely heteronuclear settings consisting of two different isotopes. Such systems have already been experimentally investigated in 3D deploying the isotopes \(^{87}\)Rb and \(^{41}\)K [45], [74] or \(^{87}\)Rb and \(^{23}\)Na [75] but also in 1D with the former composition [47]. It is a known fact that the ensuing 3D eGPEs are far more complicated compared to the homonuclear ones due to the arguably complex form of the LHY contribution. However, thus far, in 1D (to the best of our knowledge) the extraction of the eGPEs remains elusive. Below, we elaborate on their construction and provide corresponding numerical results within the eGPEs framework.

Specifically, in the most general case of a heteronuclear mixture (\(m_A\neq m_B\)) the integral encompassing the LHY contribution represented by Eq. 11 cannot be calculated in a closed analytical form. An interesting limit arises for \(\frac{g_{B}n_B}{g_{A}n_A} = \frac{m_B}{m_A}\) and \(g_{AB} \approx 0\). Recall that the first of these conditions when \(m_A=m_B\) reduces to the usual density fixing relation where the two-component eGPE framework reduces to a single-component one [20], [53]. Retaining both (interaction) conditions, the LHY energy can be calculated analytically acquiring the closed form \[\label{MassImbApprox} \mathcal{E}_{LHY} = \frac{E_{LHY}}{L} = - \frac{m_A^2 + m_B^2}{2\hbar\pi} \Big( \frac{g_An_A}{m_A} + \frac{g_Bn_B}{m_B} \Big)^{\frac{3}{2}}.\tag{19}\] However, this scenario is somewhat specialized. It requires the presence of appreciable experimental tunability in terms of Feshbach resonances that is not yet reached for these recently realized systems. As such, in what follows, we focus on general heteronuclear setups that do not satisfy the above-discussed conditions.

More concretely, we first and foremost numerically calculate the LHY energy density given by Eq. (11 ) and hence \(\frac{\partial \mathcal{E}_{LHY}}{\partial n_\sigma}\) to be incorporated in the corresponding eGPEs featuring intercomponent mass imbalance, \(m_A \neq m_B\). For the sake of comparison, we will also contrast the above-discussed ground state droplet densities with the ones obtained from either the original (Eq. 5 ) or exact (Eq. 14 ) eGPEs (exact only in the mass-balanced case, \(m_A=m_B \equiv m\)), by directly substituting the different masses of each component in the respective equation such that we can “artificially" emulate the heteronuclear setting. We mainly focus on moderate mass-imbalances with \(m_B=2m_A\) to unveil the basic droplet structures and afterwards comment on the impact of increasing imbalance.

For sufficiently strong intercomponent attraction, see Fig. 3 (a), namely within the Gaussian droplet regime, the exact homonuclear eGPEs prediction is in excellent agreement with the fully general numerical approach. Deviations among the two approaches arise and progressively increase for decreasing intercomponent attraction, where the system transitions towards the FT region, see Fig. 3 (b) where \(g_{AB}=-0.4g_A\). In both cases intercomponent mixing is minor, i.e. the droplet profiles are close to one another, and becomes gradually more pronounced as we move to the noninteracting limit presented in Fig. 3 (c) where the droplets maintain their FT structure. In contrast, for repulsive \(g_{AB}\), the heteronuclear system prefers to phase separate into a central droplet occupied by the heavier component, while the lighter species fragments into two FT droplets on either side [Fig. 3 (d)]. This preference of the lighter component to allocate into two different fragments is attributed to its reduced kinetic energy \(\sim m_\sigma g_\sigma^2\) as compared to the corresponding one of the heavier component (see also the discussion at the end of Sec. 3.1). In this case we observe a relatively small disagreement (hardly visible) with the fully general numerical approach, primarily at the edges of the droplet. As in the case of homonuclear mixtures we enforce parity symmetry, while asymmetric configurations are addressed in Appendix 8.

On the other hand, the original eGPEs [Eq. 5 ] completely miss this phase-separation behavior, predicting miscible droplets in the repulsive limit, while overestimating (underestimating) as before the droplet peak (width). For completeness, note that by employing a lighter species for component B (e.g. \(m_B=0.5m_A\)) results in the same to the above-discussed phenomenology with the droplet configurations reversed while occupying overall larger length scales, compare for instance Fig. 3 (b) and (e) with \(m_B=2m_A\) and \(m_B=0.5m_A\) respectively. In general, we find that both the exact (for \(m_A=m_B \equiv m\)) eGPEs [Eq. 14 ] and the full numerical treatment of the heteronuclear mixture do not unveil fundamentally altered droplet configurations with respect to the mass-imbalance. As an example, for \(m_B=0.2m_A,\; 0.5m_A, \; 2m_A\), or \(5m_A\) variations mainly the droplet localization are affected (not shown for brevity). Instead, the original eGPEs [Eq. 5 ] exhibit high sensitivity to the mass ratio predicting significantly different behavior, see for instance Fig. 3 (d), (f) with \(m_B=0.5m_A\) and \(m_B=5m_A\) where structural deformations are at play.

Concluding our investigation on the characterization of 1D two-component droplet setups, it is worth mentioning explicitly that neglecting the LHY correction and hence reducing our description to the standard GPEs the system displays the properties of a uniform gas. This in part stems from the absence of wave collapse in 1D. Indeed, all the droplet structures that we revealed as well as the ones to be discussed in the following Sections are inherently owed to beyond MF effects. This means that they represent a direct manifestation of the involvement of quantum fluctuations as accounted for by the perturbative LHY correction.

4 Three-component mixture↩︎

Here, we extend our considerations on the construction of the relevant eGPEs to three-component 1D droplet settings. For their realization, three different states of \(^{39}\)K could be potentially employed. The understanding and properties of such states are largely unexplored in all spatial dimensions [62], [63], [65], [67], while their 1D description through eGPEs taking the appropriate LHY correction is, to the best of our knowledge, missing. Following the standard assumptions, we only consider two-body contact interparticle interactions and restrict the Hamiltonian to include up to bilinear combinations of operators. It is then straightforward to obtain the Hamiltonian for the three-component mixture, by employing the corresponding field operators (\(\hat{\Psi}_A(x)\), \(\hat{\Psi}_B(x)\), \(\hat{\Psi}_C(x)\)) as was done in Eq. 6 . Along these lines, the three-component Hamiltonian is readily found to be (see also Ref. [62] for the 3D case) \[\label{LHYHamilt3Comp} \begin{align} &\hat{\mathcal{H}} = \frac{g_{A} N_A^2}{2L} + \frac{g_{B} N_B^2}{2L} + \frac{g_{C} N_C^2}{2L} + \frac{g_{AB} N_A N_B}{L} + \frac{g_{AC} N_A N_C}{L} + \frac{g_{BC} N_B N_C}{L} \\ &+ \sum_{k\neq 0} \Big[ \big( \frac{\hbar^2k^2}{2m_A} + g_{A} n_A \big) \hat{a}_k^\dagger \hat{a}_k + \big( \frac{\hbar^2 k^2}{2m_B} + g_{B} n_B \big) \hat{b}_k^\dagger \hat{b}_k + \big( \frac{\hbar^2 k^2}{2m_C} + g_{C} n_C \big) \hat{c}_k^\dagger \hat{c}_k \Big ] \\ & + \frac{1}{2} \sum_{k\neq 0} \Big[ g_{A} n_A \big( \hat{a}_k \hat{a}_{-k} + \hat{a}_k^\dagger \hat{a}_{-k}^\dagger \big) +g_{B} n_B \big( \hat{b}_k \hat{b}_{-k} + \hat{b}_k^\dagger \hat{b}_{-k}^\dagger \big)+ g_{C} n_C \big( \hat{c}_k \hat{c}_{-k} + \hat{c}_k^\dagger \hat{c}_{-k}^\dagger \big) \Big] \\ & +g_{AB}\sqrt{n_A n_B} \sum_{k\neq 0} \big( \hat{a}_k \hat{b}_{-k} + \hat{a}_{k}^\dagger\hat{b}_{-k}^\dagger + \hat{a}_{k} \hat{b}_k^\dagger+ \hat{a}_k^\dagger \hat{b}_{k} \big) \\ &+g_{AC}\sqrt{n_A n_C} \sum_{k\neq 0} \big( \hat{a}_k \hat{c}_{-k} + \hat{a}_{k}^\dagger\hat{c}_{-k}^\dagger + \hat{a}_{k} \hat{c}_k^\dagger+ \hat{a}_k^\dagger \hat{c}_{k} \big) \\ & +g_{BC}\sqrt{n_B n_C} \sum_{k\neq 0} \big( \hat{b}_k \hat{c}_{-k} + \hat{b}_{k}^\dagger\hat{c}_{-k}^\dagger + \hat{b}_{k} \hat{c}_k^\dagger+ \hat{b}_k^\dagger \hat{c}_{k} \big) + \mathcal{O}(\hat{a}_k^3). \end{align}\tag{20}\] Next, by defining the set of operators \(\Phi ^\dagger=(\hat{a}_k^\dagger,\hat{b}_k^\dagger,\hat{c}_k^\dagger,\hat{a}_{-k},\hat{b}_{-k},\hat{c}_{-k})\), which satisfy the bosonic commutation condition \(\mathcal{M}_b = [\Phi,\Phi^\dagger] = \begin{pmatrix} I_{3} & 0 \\ 0 & -I_{3} \end{pmatrix}\), where \(I_{3}\) is the \(3 \times 3\) identity matrix, we can express the Hamiltonian of Eq. 20 in the following bilinear form \[\label{bilinearThreeComp} \begin{align} \hat{H}= & \sum_{k>0} \Phi^\dagger H \Phi + \frac{g_{A}N_A^2}{2L} + \frac{g_{B}N_B^2}{2L} + \frac{g_{C} N_C^2}{2L} + \frac{g_{AB} N_A N_B}{L} + \frac{g_{AC} N_A N_C}{L} + \frac{g_{BC} N_B N_C}{L} \\ &-\sum_{k>0} \sum_{\sigma} \Big( \frac{\hbar^2 k^2}{2m_\sigma} + g_{\sigma}n_\sigma \Big). \end{align}\tag{21}\] The respective Hamiltonian matrix (\(H\)), due to its involved form is provided in Appendix 7. As it was shown in Sec. 2.2, this Hamiltonian can be diagonalized by solving its characteristic equation (see also Eq. 3 ) \(\det{(\mathcal{M}_b H - \omega_k I_{6})}=0\), which yields the characteristic polynomial \[\label{CharactericticPolynomial3Comp} \begin{align} & \omega^6 - (\epsilon_A^2 +\epsilon_B^2 +\epsilon_C^2 )\omega^4 + \big[\epsilon_A^2\epsilon_B^2 +\epsilon_A^2\epsilon_C^2 +\epsilon_B^2\epsilon_C^2 -4(J_{AB} + J_{AC}+ J_{BC}) \big]\omega^2 + \\ &+ \big[ - \epsilon_A^2 \epsilon_B^2 \epsilon_C^2 +4 ( J_{AB}\epsilon_C^2 + J_{AC}\epsilon_B^2+ J_{BC}\epsilon_A^2) - 16 \mathop{\mathrm{sgn}}{(g_{AB}g_{AC}g_{BC})}\sqrt{\left| J_{AB} J_{AC} J_{BC} \right|}) \big] = 0. \end{align}\tag{22}\] As expected, it contains the single-component Bogoliubov energies \(\epsilon_\sigma^2 = (\frac{\hbar^2 k^2}{2m_\sigma})^2 + \frac{\hbar^2 k^2}{m_\sigma}g_{\sigma\sigma}n_{\sigma}\) and \(J_{\sigma\sigma^\prime} = \frac{\hbar^2 k^2}{2m_\sigma}\frac{\hbar^2 k^2}{2m_\sigma^\prime} g_{\sigma\sigma^\prime}^2n_{\sigma}n_{\sigma^\prime}\) which represents the interspecies coupling contribution.

The roots (\(\omega_1^2, \; \omega_2^2, \; \omega_3^2\)) of the characteristic polynomial can be calculated analytically, albeit they have a rather involved form in the general case. Accordingly, the ground state energy becomes \[\label{LHYThreeComponent} \begin{align} E_0 =&\frac{g_{A}N_A^2}{2L} + \frac{g_{B}N_B^2}{2L}+ \frac{g_{C}N_C^2}{2L}+ \frac{g_{AB}N_A N_B}{L} + \frac{g_{AC}N_A N_C}{L}+ \frac{g_{BC}N_B N_C}{L} \\ &+\sum_{k>0} \Big( \omega_{1} + \omega_{2} +\omega_{3} - \frac{\hbar^2 k^2}{2m_A} -\frac{\hbar^2 k^2}{2m_B} - \frac{\hbar^2 k^2}{2m_C} - g_{A}n_A - g_{B}n_B - g_{C}n_C\Big). \end{align}\tag{23}\] Apparently, already at the MF level, the system is characterized by 12 independent variables spanning the parameter space (\(g_{\sigma\sigma\prime}\), \(n_\sigma\), \(m_\sigma\)). This renders rather challenging the systematic exploration of ensuing phases and development of a unified understanding in the general case. Hence, it is more convenient and instructive (at least for the purpose of this work) to address only certain characteristic cases, both from the MF and even more so from the LHY perspective.

4.1 Fully symmetric mixture↩︎

The simplest reduction of the three-component system corresponds to that of a fully symmetric mixture, which can be retrieved by requiring \(m_\sigma = m\), \(g_{\sigma\sigma} = g\), \(g_{\sigma\sigma\prime}=G\), and \(n_{\sigma} = n\) for every \(\sigma\neq\sigma\prime\). It is important to first examine this reduction since, despite its simplicity, it offers a useful benchmark with the respective single-component droplet system [31], [53], while deviating from it allows to build-up a systematic understanding of the more complex three-component setting. Employing the above-mentioned assumptions from Eq. (23 ) we can readily obtain the underlying ground state energy density \[\label{LHYThreeComponentSymetric} \begin{align} \mathcal{E} = \frac{E_0}{L} &= \frac{3g n^2}{2} + 3G n^2 +\frac{1}{L}\sum_{k>0} \Big( \omega_{1} + \omega_{2} +\omega_{3} - \frac{3\hbar^2 k^2}{2m} - 3gn \Big) = \\ & 3(\frac{g}{2}+G)n^2 \underbrace{- \frac{2(mg^3n^3)^{\frac{1}{2}}}{3\hbar\pi} \Big[ (1+ \frac{2G}{g})^{\frac{3}{2}} + 2(1- \frac{G}{g})^{\frac{3}{2}} \Big]}_{\mathcal{E}_{\rm LHY}^{\rm 3comp}}, \end{align}\tag{24}\] where \(\omega_{1}= \omega_{2} = \sqrt{\epsilon^2-2J}\) and \(\omega_{3} = \sqrt{\epsilon^2 + 4J}\)⁴ with Bogoliubov energies \(\epsilon^2 = (\frac{\hbar^2 k^2}{2m})^2 + \frac{\hbar^2 k^2}{m}gn\) and interspecies coupling contribution \(J = \frac{\hbar^2 k^2}{2m}Gn\). Finally, for convenience we define the dimensionless density independent parameter \[F(G/g) = \mathcal{E}_{\textrm{LHY}}^{\textrm{3comp}}/\sqrt{mg^3n^3\hbar^{-2}}=-\frac{2}{3\pi} \Big[ (1+ 2G/g)^{3/2} + 2(1- G/g)^{3/2} \Big] <0.\] Note that as long as \(-g/2 \leq G \leq g\) (which corresponds to the MF stability regime as discussed in Appendix 7) the eigenvalues \(\omega_{i}\) are real. Also, in line with the two-component case, the symmetric three-component mixture reduces to an effectively single-component one. This effectively single-component system, however, depends on two interaction parameters, namely \(g\) and \(G\). It is inherently different from the corresponding single field reduction stemming from the two-component setup exclusively due to the respective LHY corrections (see the discussion below).

The LHY energy of this reduced system, again remains \(\sim n^{3/2}\) throughout the MF stability regime \(G\in[-g/2, g]\). Moreover, the density independent proportionality factor ranges from \(F(G/g)= -\sqrt{6}/\pi\) for \(G=-g/2\), to \(F(G/g)= -2\sqrt{3}/\pi\) for \(G=g\), while it takes its maximum \(F(G/g)= - 2/\pi\) value at \(G=0\), as shown in Fig. 4 (a). Interestingly, the strength of quantum fluctuations appears to be larger close to the immiscibility threshold \(G=g\) as compared to the one close to the droplet regime \(G=-g/2\), as evidenced by the more strongly attractive LHY energy (see e.g. Fig. 4 (a) and Eq. 24 ). This is in contrast to the two-component mixture where quantum fluctuations contribute equally at the two limits, see in particular Fig. 1 (e) and Eq. 13 . We attribute this behavior of the three-component system to the aforementioned “asymmetric” MF stability regime (see also Appendix 7 for a more elaborated discussion on its asymmetric nature).

On the other hand, the MF interaction term vanishes at the droplet threshold \(G\approx-g/2\), while it takes its maximum value at the immiscibility threshold \(g=G\). Hence, it is not a-priori obvious how important the role of the LHY correction becomes at each limit (e.g. from Eq. 24 ) and one has to solve the corresponding eGPE, as we shall do below. Finally, it is interesting to note that in the fully symmetric case the corresponding LHY energy of the original three-component mixture is always smaller than the one of the reduced two-component system, namely it holds that \(\left| \mathcal{E}_{\textrm{LHY}}^{\textrm{3comp}} \right| \geq \frac{3}{2} \left| \mathcal{E}_{\textrm{LHY}}^{\textrm{2comp}} \right|\) with the equality being valid for \(G=0\). As such, it appears that the effect of quantum fluctuations is enhanced in the case of a three-component mixture, even beyond the additive effect owing to the additional component, due to the interspecies coupling. Instead, for the MF energy it holds that \(\mathcal{E}_{\textrm{MF}}^{\textrm{3comp}}/3 = \mathcal{E}_{\textrm{MF}}^{\textrm{2comp}}/2\), i.e. the ratio of the MF energy divided by the number of components remains constant for any number of components in the symmetric case.

To obtain the effective single-component eGPE equation emanating from the fully symmetric three-component mixture we evaluate the Euler-Lagrange equations (from Eq. 24 ) with respect to \(n\). This standard procedure leads to the symmetric eGPE:

\[\label{eGPEThreeComponentSymmetric} \begin{align} 3i\hbar\frac{\partial \Psi}{\partial t} =& - 3\frac{\hbar^2}{2m} \frac{\partial^2 \Psi}{\partial x^2} + 3(g+2G) \left| \Psi \right|^2\Psi + F(G/g) \frac{3\sqrt{mg^3}}{2 \hbar} \left| \Psi \right| \Psi . \end{align}\tag{25}\]

Similarly to the case of a symmetric two-component mixture, this symmetric eGPE accepts a droplet solution with saturation density \(n_s = F^2(G/g)\frac{mg^3}{9\hbar^2(2G + g)^2}\), energy density \(\mathcal{E}^{3comp}(n=n_s) = \\- \frac{3}{2} (g + 2G) n_s^2\), chemical potential density \(\mu_s/L = - \frac{3}{2} (g+2G) n_s\), and healing length \(\xi = \sqrt{ - \frac{3\hbar^2}{2m\mu_s/L} } \\= \sqrt{\frac{9\hbar^4(g+2G)}{m^2g^3F^2(G/g)}} = \sqrt{\frac{\hbar^2}{m(g+2G)n_s}}\). Again, in accordance with the two-component symmetric mixture, the ground state density of the symmetric three-component system transits from a Gaussian-type profile (for \(G\approx -1/2g\)) to a FT one for decreasing intercomponent attraction as captured by the \(G\) coefficient. This behavior is explicitly illustrated in Fig. 4 (b) for a fixed atom number⁵. Additionally, the ground state density maintains its FT character throughout the repulsive interaction regime, \(G>0\), see Fig. 4 (c) and (d). The saturation density [width] decreases [increases] by over an order of magnitude as \(G \approx g\), eventually acquiring the minimum saturation density \(n_s(G=g) = \frac{4mg}{27\pi^2\hbar^2}\approx 0.015\) deep in the repulsive interaction regime as can be seen in Fig. 4 (d).

4.2 Two identical components↩︎

As it was argued above, already in the limiting case of a fully symmetric three-component mixture evident differences arise (concerning, for instance, the droplet saturation density, and the LHY correction) when compared to the corresponding reduced two-component system. However, by construction three-component mixtures offer additional possibilities for intercomponent asymmetry, that could lead to intriguing phenomena and phases that can not be captured within the two-component setting. A step forward to delve more into the three-component properties is to consider a system where two of the components are identical (i.e. symmetric), while the third one is different. To be concrete, this situation is described by \(m_A=m_B \equiv m\), \(g_{A}=g_{B} \equiv g\), and \(g_{AC}=g_{BC} \equiv G_C\), while \(g_{AB}\) and \(g_{C}\) are unrestricted parameters, see also Ref. [62] for the 3D case.

In this scenario, it can be found that the MF stability requires the following conditions to be fulfilled: i) \(g,g_C > 0\), ii) \(gg_{C} > G_C^2\), and \(g > \left| g_{AB} \right|\) as well as iii) \(G_C^2< \frac{g_{C}}{2}(g+g_{AB})\), see also Appendix 7 for further details. We remark that the same conditions hold for the respective 3D system as reported in Ref. [62], which is to be expected because the homogeneous MF energy is independent of the dimensionality. It is also important to emphasize that the last (iii) condition is more restrictive than the second (ii) one since \(g_{C}(g+g_{AB})/2 \leq g g_{C}\) if \(\left| g_{AB} \right|<g\). As a consequence, it is possible for each pair of subsystems to be stable, while the overall system is outside the MF stability region, a behavior that equally holds in 3D [62]. For simplicity, here, we focus on interaction intervals where each subsystem and the overall mixture reside within the MF stability region.

Further assuming that \(m_C=m\) (for convenience), we find that the LHY energy density for the general three-component system in Eq. 23 takes the form \[\begin{align} \label{ELHY2Plus1} \mathcal{E}_{\textrm{LHY}} =& - \frac{\sqrt{2 m}}{6\hbar\pi} \Big( (2 (g-g_{AB})n)^{\frac{3}{2}} + (W + Q)^{\frac{3}{2}} + (W - Q)^{\frac{3}{2}} \Big), \end{align}\tag{26}\] with the parameters \[\begin{align} \label{ELHY2Plus1Params} Q =& \Big( (g+g_{AB})^2n^{2} - 2 (g + g_{AB}) g_{C} n n_C + 8 G_C^{2} n n_C + g_{C}^{2} n_C^{2} \Big)^{\frac{1}{2}},\\ W =& g n + g_{AB} n + g_{C} n_C. \end{align}\tag{27}\] The LHY energy density is real and negative as long as \(g>\left| g_{AB} \right|\) and \(W - Q > 0 \xrightarrow{} 2 G_C^2 < g_{C}(g + g_{AB})\), which is exactly the MF stability regime, while it becomes complex otherwise. Moreover, by removing the third component namely setting \(G_C=g_{C}=n_C=0\), it holds that \(Q=W=(g+g_{AB})n\) and the LHY energy density becomes \(\mathcal{E}_{\textrm{LHY}} = \frac{-2\sqrt{m}(gn)^{3/2}}{3\hbar\pi}[(1-g_{AB}/g)^{3/2} + (1+g_{AB}/g)^{3/2}]\). This is exactly the LHY energy density of the symmetric two-component mixture in the MF stability regime (see Eq. 12 for \(g_A=g_B \equiv g\) and \(n_A=n_B \equiv n\)).

As in Sec. 3.1, we can derive the LHY terms of the effective two-component eGPEs describing the system by calculating their partial derivatives \[\begin{align} \label{LHYTerms3Comp} \frac{\partial\mathcal{E}_{\textrm{LHY}}}{\partial n} =& - \frac{\sqrt{2m}}{4\pi h } \Bigg( \Big( g + g_{AB} + \frac{Z}{Q} \Big) \sqrt{W - Q} + \Big( g + g_{AB} - \frac{Z}{Q} \Big) \sqrt{W + Q} + 2 \sqrt{2} (g - g_{AB})^{\frac{3}{2}} \sqrt{n} \Bigg), \\ \frac{\partial\mathcal{E}_{\textrm{LHY}}}{\partial n_C} =& - \frac{\sqrt{2m}}{4\pi h } \Bigg( \Big( g_C + \frac{P}{Q} \Big) \sqrt{W - Q} + \Big( g_C - \frac{P}{Q} \Big) \sqrt{W + Q} \Bigg), \end{align}\tag{28}\] where \[\begin{align} \label{2Plus1DerParams} Z = & -\frac{1}{2}\frac{\partial Q^2}{\partial n} = - g^{2} n - 2 g g_{AB} n + g g_C n_C - g_{AB}^{2} n + g_{AB} g_C n_C - 4 G_C^{2} n_C, \\ P = & -\frac{1}{2}\frac{\partial Q^2}{\partial n_C} = g g_{C} n + g_{AB} g_{C} n - 4 G_C^{2} n - g_{C}^{2} n_C. \end{align}\tag{29}\] These LHY terms are incorporated in the corresponding eGPEs, which by additionally requiring \(\Psi_A= \Psi_B \equiv \Psi\), read \[\begin{align} \label{eGPETwoPlusOneComponent} 2i\hbar\frac{\partial \Psi}{\partial t}& = - \frac{\hbar^2}{m} \frac{\partial^2 \Psi}{\partial x^2} +\Big( 2(g+g_{AB}) \left| \Psi \right|^2 + 2G_C \left| \Psi_C \right|^2 +\frac{\partial \mathcal{E}_{\textrm{LHY}}}{\partial n} \Big)\Psi, \\ i\hbar\frac{\partial \Psi_C}{\partial t}& = - \frac{\hbar^2}{2m} \frac{\partial^2 \Psi_C}{\partial x^2} + \Big( g_{C}\left| \Psi_{C} \right|^2 + 2G_C\left| \Psi \right|^2 + \frac{\partial \mathcal{E}_{\textrm{LHY}}}{\partial n_C} \Big)\Psi_C. \end{align}\tag{30}\] A few characteristic examples of the resulting ground state droplet configurations are presented in Fig. 5, with fixed \(g_{C}=g\), and \(g_{AB}=-0.5g\). Specifically, the case of \(G_C=-0.5g\) refers to the respective LHY-fluid limit (attained in general here for \(G_C=-g/2\) due to MF stability requirements) of the fully symmetric mixture where MF interactions cancel out. Here, the individual components become identical and possess a Gaussian type distribution, see Fig. 5 (a). As we decrease the intercomponent attraction, the system transitions to a mixed state where the symmetric components (\(A\), \(B\)) exhibit a FT configuration and the third one (\(C\)) becomes significantly distorted and delocalizes. Particularly, there is a FT droplet segment lying within the symmetric components and the remaining atoms of the third component reside outside this region in a self-bound state as shown in Fig. 5 (b), (c). However, by switching-off the intercomponent interactions among the symmetric components and the third one, i.e. reaching the decoupling limit with \(G_C=0\), the subsystems feature two independent droplet structures [Fig. 5 (d)].

On the other hand, turning to the repulsive side of the MF stability regime (e.g. for \(G_C = 0.1g\), or \(G_C = 0.4g\)), it turns out that the third component progressively separates from the symmetric two-component droplet and accommodates two smaller sized droplets on either side of the symmetric two-component droplet which remains localized at the center [Fig. 5 (e), (f)] (see also Appendix 8 for the parity asymmetric configurations). Finally, it is important to emphasize that the above-described investigation is far from being exhaustive for the three-component system. A more detailed examination of the emerging droplet states and mixed phases of matter e.g. by varying \(N_C\), \(g_C\), \(g_{AB}\), and \(G_C\) across the MF stability regime is certainly an interesting direction to be pursued in future studies.

5 Four-component mixture↩︎

As can be deduced from the above analysis, in order to obtain the ground state energy of the quasi-particle vacuum one needs to diagonalize the Hamiltonian matrix which has \(2\times n\) dimensions, with \(n\) representing the number of bosonic species. Due to the symmetries of the Hamiltonian (e.g. being a real Hermitian matrix) the resulting characteristic equation is a polynomial of degree \(n\), with respect to the eigenvalues \(\omega^2\). Therefore, the four-component mixture is the most complex system for which the characteristic equation can be analytically solved and hence extract its LHY energy. For higher-component settings the characteristic equation can only be solved numerically in the general case. As such, for completeness, we present below the LHY energy for the four-component mixture, even though such a setup would be arguably more challenging to be experimentally realized.

The process for deriving the energy of the four-component mixture is identical to the ones of the two- (Sec. 3) and three- (Sec. 4) component ones described above. For brevity, we provide directly the resulting characteristic polynomial \[\label{CharactericticPolynomial4Comp} \omega^8 + \alpha \omega^6 + \beta \omega^4 +\gamma \omega^2 + \delta = 0,\tag{31}\] containing the rather involved coefficients

\[\begin{align} \label{CharactericticPolynomial4CompCoeff} \alpha = & -\sum_{\sigma} \epsilon_\sigma^2, \\ \beta = & \sum_{\sigma = A}^D \sum_{\sigma^\prime > \sigma}^D \Big( \epsilon_\sigma^2\epsilon_{\sigma^\prime}^2 - 4J_{\sigma\sigma^\prime} \Big),\\ \gamma = & \sum_{\sigma }^D \Bigg[- \epsilon_\sigma^2 \Big[ \sum_{\sigma^\prime > \sigma}^D \sum_{\sigma^{\prime\prime} > \sigma^\prime}^D \epsilon_{\sigma^\prime}^2\epsilon_{\sigma^{\prime\prime}}^2 + 4\sum_{\sigma^\prime \neq \sigma}^D \sum_{\sigma^{\prime\prime} > \sigma^\prime, \sigma^{\prime\prime} \neq \sigma }^D J_{\sigma^\prime\sigma^{\prime\prime}} \Big] \\ &+16 \prod_{ {\sigma^\prime \neq \sigma}, \sigma^\prime<\sigma^{\prime\prime} } \mathop{\mathrm{sgn}}{(g_{\sigma^\prime\sigma^{\prime\prime}})} \sqrt{J_{\sigma^\prime\sigma^{\prime\prime}}} \Bigg] \notag, \\ \delta = & \prod_{\sigma} \epsilon_\sigma^2 -4 \sum_{\mathcal{P}} \Big[ \big( \epsilon_{\sigma_1}^2\epsilon_{\sigma_2}^2 -4 J_{\sigma_1\sigma_2} \big)J_{\sigma_3\sigma_4} \Big] + 16 \sum_\sigma \epsilon_\sigma^2 \prod_{ \sigma_{i,j} \neq \sigma, \sigma_i < \sigma_j} \frac{\sqrt{|J_{\sigma_i\sigma_j} }}{\mathop{\mathrm{sgn}}{(g_{\sigma_i\sigma_j}) }} \\ & -32 \frac{\sqrt{|J_{AC}J_{AD}J_{BC}J_{BD} |} }{\mathop{\mathrm{sgn}}{(g_{AC}g_{AD}g_{BC}g_{BD} )}} -32 \frac{\sqrt{|J_{AB}J_{AD}J_{BC}J_{CD} |} }{\mathop{\mathrm{sgn}}{(g_{AB}g_{AD}g_{BC}g_{CD} )}} -32 \frac{\sqrt{|J_{AB}J_{AC}J_{BD}J_{CD} |} }{\mathop{\mathrm{sgn}}{(g_{AB}g_{AC}g_{BD}g_{CD} )}} \notag, \end{align}\tag{32}\] where \(\epsilon_\sigma^2 = (\frac{\hbar^2 k^2}{2m_\sigma})^2 + \frac{\hbar^2 k^2}{m_\sigma}g_{\sigma\sigma}n_{\sigma}\) are the single component Bogoliubov energies and \(J_{\sigma\sigma^\prime} =\\ \frac{\hbar^2 k^2}{2m_\sigma}\frac{\hbar^2 k^2}{2m_\sigma^\prime} g_{\sigma\sigma^\prime}^2n_{\sigma}n_{\sigma^\prime}\) the interspecies coupling contributions. The notation \(\mathcal{P}\) in the sums refers to all possible permutations of \(\sigma_1, \sigma_2, \sigma_3, \sigma_4\) such that \(\sigma_1<\sigma_2\), and \(\sigma_3<\sigma_4\) and thus index duplication is avoided.

The solutions (\(\omega_i\)) of the characteristic polynomial equation can be calculated exactly and the resulting ground state energy for the four-component quasi-particle vacuum reads \[\label{LHY4Component} \begin{align} E_0 =&\sum_{\sigma} \frac{g_{\sigma\sigma}N_\sigma^2}{2L} + \frac{1}{2}\sum_{\sigma\neq\sigma^{\prime}} \frac{g_{\sigma\sigma^\prime}N_\sigma N_\sigma^\prime}{L}+ \sum_{k>0} \Bigg[\sum_{i=1}^{4} \omega_{i} - \sum_\sigma \Big( \frac{\hbar^2 k^2}{2m_\sigma} + g_{\sigma\sigma}n_\sigma \Big) \Bigg]. \end{align}\tag{33}\] It is apparent that for this setup the parametric space as defined by the different intra- and inter-component coupling combinations, the atom number and mass of each component becomes exceedingly large. While a systematic study of this system is interesting on its own right, in what follows we restrict ourselves to the relatively simpler case of the symmetric mixture.

5.1 Symmetric four-component mixture↩︎

Let us consider a fully symmetric four-component mixture characterized by \(N_\sigma \equiv N\), \(g_{\sigma\sigma} \equiv g\), and \(g_{\sigma\sigma\prime} \equiv G_4\). Here, the MF energy (first two sums of Eq. (33 )) reduces to \(\mathcal{E}^{\rm 4comp}_{\textrm{MF}}=(2g+6G_4)n^2\), whilst the MF stability regime is given by \(-g/3<G_4<g\), ensuring the stability of each subsystem but also of the full four component mixture. Accordingly, the roots of the characteristic polynomial of Eq. 31 simplify to \(\omega_{1,2,3}^2 = \epsilon^2 - G_4\frac{\hbar^2k^2}{m}n\), and \(\omega_{4}^2 = \epsilon^2 + 3G_4\frac{\hbar^2k^2}{m}n\) and the LHY energy density becomes \[\label{LHY4CompSymmetric} \begin{align} \mathcal{E}_{\textrm{LHY}}^{\rm 4comp} &= - \frac{2\sqrt{mg^3n^3}}{\pi\hbar} \Bigg( \Big(1 - \frac{G_4}{g}\Big)^\frac{3}{2} + \frac{1}{3}\Big(1 + 3\frac{G_4}{g}\Big)^\frac{3}{2} \Bigg). \end{align}\tag{34}\]

As in the case of the symmetric two-component mixture [Sec. 3.1], the emergent droplet solution can be shown to have saturation density \(n_s = \frac{mg^3}{4\pi^2\hbar^2(3G_4 + g)^2}\Big( (1-\frac{G_4}{g})^{3/2} + \frac{1}{3}(1+3\frac{G_4}{g})^{3/2}\Big)\), energy density \(\mathcal{E}^{\rm 4comp}(n=n_s) = - 2 (g + 3G_4) n_s^2\), chemical potential density \(\mu_s/L = -2 (g+3G_4) n_s\), and healing length \(\xi = \sqrt{ - \frac{\hbar^2}{2m\mu_s/L} } = \sqrt{\frac{\hbar^2}{4m(g+3G_4)n_s}}\). Moreover, similarly to the three-component setup, the LHY energy of the fully symmetric four-component mixture is always smaller than twice the LHY energy of the corresponding symmetric two-component system, namely \(\mathcal{E}_{\textrm{LHY}}^{\textrm{4comp}} \\ \leq 2 \mathcal{E}_{\textrm{LHY}}^{\textrm{2comp}}\) with the equality being valid for \(G_4=0\). This implies the presence of enhanced quantum fluctuations in the symmetric four-component system, see also Fig. 4 (a), beyond an additive contribution. On the other hand, for the MF energies it holds that \(\mathcal{E}_{\textrm{MF}}^{\textrm{4comp}} = 2 \mathcal{E}_{\textrm{MF}}^{\textrm{2comp}}\). Representative ground state droplet profiles of the reduced four-component setting are illustrated in Fig. 6 for different intercomponent couplings (\(G_4\)) and fixed intracomponent ones. In line with the behavior of the droplets arising in lower-component mixtures also here a deformation from a Gaussian to a FT density profile occurs for decreasing intercomponent attractions as depicted in Fig. 6 (a). Finally, at the repulsive side, increasing the intercomponent interaction across the MF stability regime results in quantum droplet structures of significantly larger spatial width, as shown in Fig. 6 (b).

6 Conclusions and future directions↩︎

We derived, in a systematic manner, the first-order beyond mean-field LHY quantum correction and associated eGPEs for two-, three- and four-component 1D bosonic mixtures featuring short-range contact interactions. The four-component system represents the highest multicomponent setting for which the LHY energy can be analytically derived in the general case. In all cases, whenever the complexity of the system allows corresponding closed form expressions are provided. Importantly, our treatment relies only on the assumption of a dilute gas such that the perturbative Bogoliubov treatment is valid. Hence, the LHY energies and the ensuing eGPEs obtained are reliable throughout the MF stability regime, namely from the attractive to the repulsive intercomponent interaction regimes in contrast to previous treatments.

Focusing on the two-component mixture, we extracted closed form solutions for the LHY energy and constructed the corresponding exact coupled system of eGPEs for homonuclear settings and provide numerical results for heteronuclear 1D systems. Specifically, by tuning the intercomponent interactions towards smaller attractions, we retrieve the structural deformation of the droplets from Gaussian to FT density configurations. However, it is explicated that the previously discussed original two-component LHY treatment reveals quantitative deviations as compared to our exact eGPE results and can be obtained as a limiting case from the exact expressions presented herein. These deviations are reflected by an increased saturation droplet density and decreasing width. Turning to repulsive interactions, we identify even more significant qualitative differences with the original LHY treatment, especially for heteronuclear settings featuring intercomponent mass-imbalance. In particular, self-bound structures with decreasing saturation densities are found for larger intercomponent repulsions, while an early onset of phase-separation occurs for both homonuclear and heteronuclear mixtures. The numerically computed heteronuclear droplet states are found to be in good agreement with the results of the exact eGPEs for homonuclear settings upon direct substitution of the different masses. Meanwhile, the original eGPEs fail to adequately capture the droplet densities.

Additionally, exact closed form expressions of the LHY energies and the eGPEs are discussed for three-component mixtures especially in the cases of i) a fully-symmetric system and ii) two identical components coupled to a third independent one. It is argued that the LHY energy is not additive to the number of components and it is in general larger for three-component setups as compared to two-component ones. A plethora of droplet phases exist, ranging from the underlying LHY-fluid limit, to Gaussian miscible droplet components as well as to mixed states for decreasing attraction. These mixed states refer to FT configurations for the symmetric components and a delocalized droplet distribution for the third one. The latter assembles in a FT droplet segment within the symmetric components and the remaining atoms lying beyond the intercomponent overlap region are in a self-bound state. Instead, for repulsive interactions the third component phase separates from the symmetric two-component droplet and exhibits two smaller sized droplets on either side. Finally, a closed LHY form is also provided for the fully symmetric four-component mixture. Here, the transition from Gaussian to FT droplets for decreasing attraction and afterwards to large spatial FT structures upon crossing to the repulsive intercomponent regime is observed.

There are several intriguing pathways to be followed in the future based on our work which provides the first steps towards exploring ground state and dynamical multicomponent droplet configurations. A straightforward extension is to systematically study the transition towards phase-separation in heteronuclear mass-imbalanced two-component mixtures and unveil the origin of this mechanism. Similarly, the three- and four-component mixtures offer completely unexplored possibilities to create exotic droplet configurations. These include, for instance, miscible spatially deformed droplet structures, and importantly mixed droplet phases thereof as well as droplet-soliton [76] or droplet-vortex [36] coexistence. The interpolation of three-component droplet states to the few-body regime using sophisticated numerical schemes is also of particular interest [77]–[79]. On the other hand, dynamical crossing of these droplet phases is expected to reveal insights into their collective excitations and the spontaneous generation of nonlinear defects or delocalized structures such as shock-waves [80], [81]. Along these lines, emulating the time-of-flight measurement process of multicomponent droplet states in order to testify their self-bound nature is certainly worth pursuing. Finally, extending our treatment to 3D geometries in order to investigate the impact of the LHY term beyond the MF stability edge as well as in the presence of arbitrary external potentials are very relevant open questions to address.

Acknowledgements↩︎

The authors thank P. G. Kevrekidis, G. Bougas and G.C. Katsimiga for extensive discussions on the topic of droplets.

6.0.0.1 Funding information

S.I.M acknowledges support from the University of Missouri Science and Technology, Department of Physics, Startup fund.

7 Details on the stability conditions of the three-component mixture↩︎

The Hamiltonian matrix for the three-component mixture used in Eq. 21 of the main text takes the form

\[H = \begin{pmatrix} \frac{\hbar^2 k^2}{2m_A} + g_{A}n_A & g_{AB}\sqrt{n_A n_B} & g_{AC}\sqrt{n_A n_C} & g_{A}n_A & g_{AB}\sqrt{n_A n_B} & g_{AC}\sqrt{n_A n_C} \\ g_{AB}\sqrt{n_A n_B} & \frac{\hbar^2 k^2}{2m_B} + g_{B}n_B & g_{BC}\sqrt{n_B n_C}& g_{AB}\sqrt{n_A n_B} & g_{B}n_B & g_{BC}\sqrt{n_B n_C}\\ g_{AC}\sqrt{n_A n_C} & g_{BC}\sqrt{n_B n_C} & \frac{\hbar^2 k^2}{2m_C} + g_{C}n_A & g_{AC}\sqrt{n_A n_C} & g_{BC}\sqrt{n_B n_C} & g_{C}n_C \\ g_{A}n_A & g_{AB}\sqrt{n_A n_B} & g_{AC}\sqrt{n_A n_C} & \frac{\hbar^2 k^2}{2m_A} + g_{A}n_A & g_{AB}\sqrt{n_A n_B} & g_{AC}\sqrt{n_A n_C}\\ g_{AB}\sqrt{n_A n_B} & g_{B}n_B & g_{BC}\sqrt{n_B n_C} & g_{AB}\sqrt{n_A n_B} & \frac{\hbar^2 k^2}{2m_B} + g_{B}n_B & g_{BC}\sqrt{n_B n_C} \\ g_{AC}\sqrt{n_A n_C} & g_{BC}\sqrt{n_B n_C} & g_{C}n_C & g_{AC}\sqrt{n_A n_C} & g_{BC}\sqrt{n_B n_C} & \frac{\hbar^2 k^2}{2m_C} + g_{C}n_A \end{pmatrix}.\]

Before exploring the LHY term in more detail, it is important to determine the MF stability regime. In general, the stability condition is given by demanding that the Hessian matrix (\(Hess\)) in terms of derivatives with respect to the densities (\(Hess_{\sigma,\sigma\prime} = \frac{\partial^2 \mathcal{E}_{\textrm{MF}} }{\partial n_\sigma\partial n_\sigma\prime}\)) has only positive eigenvalues [5]. Namely, the Hessian matrix in our case reads \[\label{Hessian3Comp} Hess = \begin{pmatrix} g_{A} & g_{AB} & g_{AC} \\ g_{AB} & g_{B} & g_{BC} \\ g_{AC} & g_{BC} & g_{C} \\ \end{pmatrix},\tag{35}\] with eigenvalues dictated by the roots of the characteristic polynomial \[\label{HessCharPol3Comp} f(x) = x^3 - x^2 \sum_\sigma g_{\sigma\sigma} + x \sum_{\sigma<\sigma\prime} \Delta g_{\sigma\sigma\prime} + C.\tag{36}\] In this expression, \(\Delta g_{\sigma\sigma\prime}= g_{\sigma} g_{\sigma\prime} - g_{\sigma\sigma\prime}^2\) and \(C = g_{A}g_{BC}^2 + g_{B}g_{AC}^2 + g_{C}g_{AB}^2 - g_{A}g_{B}g_{C} - 2 g_{AB}g_{AC}g_{BC}\). Note that the stability of the individual components requires \(\Delta g_{\sigma\sigma\prime} >0\) (and \(g_{\sigma\sigma}>0\)) for all components. Then, MF stability requires \[\label{3CompMFstability} \begin{align} &[1.]\; g_{\sigma\sigma} > 0, \quad [2.]\; \Delta g_{\sigma\sigma\prime} > 0, \quad [3.] \; C<0 \quad \\ &[4.] \; \frac{1 }{3}\sum_\sigma g_{\sigma\sigma} - \sqrt{\Delta} >0,\; {\rm where} \; \Delta = \sum_{\sigma} g_{\sigma\sigma}^2 + \sum_{\sigma\neq\sigma\prime} (g_{\sigma\sigma\prime}^2 - \Delta g_{\sigma\sigma\prime}/2 ) ~~\textrm{and}~~\\ &[5.]\; f\Big(\frac{\sum_\sigma g_{\sigma\sigma} - \sqrt{\Delta } }{3} \Big) \geq 0 \; {\rm and } \; f\Big(\frac{ \sum_\sigma g_{\sigma\sigma} + \sqrt{\Delta } }{3} \Big) \leq 0 . \end{align}\tag{37}\] In particular, the last two conditions guarantee the existence of three real solutions. However, the Hessian matrix (Eq. 35 ) is a real and symmetric (normal) matrix, and hence it is guaranteed to be diagonalizable. As such, the characteristic polynomial (Eq. 36 ) is already guaranteed to have three real roots by construction and we do not need to check the (rather complicated) conditions (5.) explicitly. Finally, condition (4.) is automatically satisfied as long as condition (2.) holds.

In the case of two symmetric components (i.e., \(g_A=g_B=g\)) coupled to the third component with the same interaction strength (namely \(g_{AC}=g_{BC}=G_C\)) the characteristic polynomial of the Hessian matrix becomes \[\label{HessCharPolTwoPlusOneComp} \begin{align} f(x) =& x^3 - x^2 (2g + g_{C}) + x (g^2 - 2G_C^2 - 2g g_{C} - g_{AB}^2) + [2G_C^2g - 2G_C^2g_{AB} - g^2 g_{C} + g_{AB}^2 g_{C}]. \end{align}\tag{38}\] Therefore, MF stability for these three-component mixtures requires \[\label{MFcondition2plus1} [1.] \; g > 0 \; {\rm and} \; g_{C}>0 \; [2.] \; \left| g_{AB} \right|> g \; {\rm and} \; gg_{C} > G_C^2 \; [3.] \; G_C^2< \frac{g_{C}}{2}(g+g_{AB}).\tag{39}\]

For the case of the fully symmetric mixture i.e. \(g_A=g_B=g_C=g\) and \(g_{AB}=g_{AC}=g_{BC}\), the system is significantly simplified and the MF stability requires \(-g/2 < G < g\). Evidently, the MF stability window becomes reduced compared to the respective two-component one, given by \(-g<g_{AB}<g\) as discussed in Sec. 3.1, and asymmetric with respect to the sign of the inter-species interaction strength. An intuitive explanation for this phenomenology is the following. Consider a mixture consisting of \(\mathcal{N}\) identical components with \(N\) particles per component, featuring intracomponent interactions \(g_\sigma =g\) and intercomponent ones \(g_{\sigma\sigma^\prime} =G\), for all \(\sigma\), \(\sigma^\prime\). Ignoring any complications stemming from the exchange statistics and distinguishability, a particle of species \(\sigma\) may scatter from \(N-1\) \(\sigma\)-species particles with strength \(g\) and \((\mathcal{N}-1)N\) atoms of other components with strength \(G\). Hence we expect that, on average, the intracomponent interaction strength will be suppressed by a factor of \(\mathcal{N}-1\) (assuming \(N \gg 1\)) compared to the intercomponent one, yielding an average interaction strength \(\propto (G+ \frac{g}{\mathcal{N}-1})\). We anticipate the stability of this effective single-component mixture to be ensured as long as its average interaction strength is repulsive (i.e. \(G \geq -\frac{g}{\mathcal{N}-1}\)), since the effective single-component picture does not impose any limitations at the repulsive side. Indeed this pattern appears to hold for \(\mathcal{N}=2,\) (see Sec. 3.1) \(\mathcal{N}=3,\) (see Sec. 4.1) and \(\mathcal{N}=4\) (see Sec. 5.1), where the lower stability bound corresponds to \(G\geq -g\), \(G\geq -g/2\), and \(G\geq -g/3\) respectively. On the other hand, the only upper bound imposed to the intercomponent interaction strength \(G\), stems from the stability of its sub-systems which requires \(G<g\) regardless of the number of components. This means that at the repulsive side all components will either be stable and fully overlap or all components will be mutually unstable and phase separate from all other components. A hypothetical partially phase separated configuration, with \(\mathcal{M}\) overlapping components in one cluster and \(\mathcal{N} - \mathcal{M}\) in the other, would feature average intra-cluster interaction strengths \(\propto G+ \frac{g}{\mathcal{M}-1}\) and \(\propto G+ \frac{g}{\mathcal{N} - \mathcal{M}-1}\) respectively, both of which are larger than the inter-cluster coupling strength \(G\) (considering only repulsive interactions). As such, there is no reason for the clusters to phase separate, which explains why such configurations do not seem to appear.

8 Asymmetric phase-separated configurations↩︎

The eGPEs derived in the main text support also ground state solutions with spontaneously broken parity symmetry in the case of mass-, or interaction-imbalanced mixtures with inter-component repulsion. Namely, these solutions are present in the cases where the phase separated configurations discussed in the main text appear as illustrated in Fig. 2 (c), Fig. 3 (d), (f), and Fig. 5 (e), (f). In Fig. 7 (a)-(e) we present the corresponding parity asymmetric solutions for all cases mentioned above. In all cases, these solutions feature slightly lower energies (e.g. the relative energy difference is \(\approx 0.3\%\) [\(\approx 0.25 \%\)] for the homonuclear [heteronuclear with \(m_B=5m_A\)] two-component mixture presented in Fig. 7 (a) [Fig. 7 (c)]) than the parity symmetric configurations shown in the main text. However, as mentioned in the main text, we consider such asymmetric phase-separated configurations to be associated with the classical-field approximation since they are not in principle permitted under a many-body ansatz. As such, within the main text we present and discuss the respective parity symmetric solutions as the ground state of the system in-spite of the slightly larger energy owing to the larger kinetic energy. We have confirmed that both types of solutions (i.e., symmetric and asymmetric) are robust when exposed to small amplitude perturbations, e.g., emulated by a random noise distribution across the spatial profile of the solution, and they afterwards let to evolve at least for total evolution times \(T= 10^4\;\hbar^3/(m_Ag_A^2)\) that we have checked. It is also important to remark that the existence of these configurations should not be under-appreciated. They can naturally emerge in corresponding experiments e.g. even in the presence of relatively small magnetic field gradients or other symmetry breaking contributions due to un-avoidable imperfection sources. Moreover, these domain-wall like distributions are interesting on their own right e.g. for dynamical applications in order to create topological excitations or study hydrodynamic instabilities.

To clarify our argument about the many-body wavefunction structure, consider as an example the homonuclear two-component mixture. Let \(\Phi_A(x)\), \(\Phi_B(x)\), be the single particle wavefunctions corresponding to the densities depicted in Fig. 7 (a). We note in passing that all the configurations shown in Fig. 7 (a)-(e) are, of course, doubly degenerate and their parity transformed counterpart is also a solution with the same energy, due to the parity symmetry of the ensuing Hamiltonian. Hence, neglecting any higher order correlations, the simplest form of the many-body wavefuntion for our homonuclear mixture, consistent with the results from the eGPEs would be⁶: \[\label{ManyBodyWfn} |\Psi\rangle = \sqrt{\frac{N_AN_B}{2}}\Bigg( |\Phi_A(x)\rangle |\Phi_B(x)\rangle + |\Phi_A(-x)\rangle |\Phi_B(-x)\rangle\Bigg).\tag{40}\] Therefore, the corresponding diagonal of the reduced single-particle density matrix (namely, \(\rho^1_{\sigma}(x) = Tr_\sigma \langle \Psi|\Psi\rangle\)) of the underlying many-body system is not the one presented in Fig. 7 (a), but it has the significantly more involved spatial profile illustrated in Fig. 7 (f). Obviously, this structure can be associated with a significantly higher energy (at the level of the classical-field approximation) than either the parity symmetric solution shown in Fig. 2 (c), or the parity asymmetric solution depicted in Fig. 7 (a). Finally, let us mention explicitly that the parity symmetric solutions discussed throughout the text can be directly associated with the reduced single-particle density of a many-body wavefunction (at the same level of approximation) since they are not degenerate. Namely, if \(\Phi_A^S(x)\), \(\Phi_B^S(x)\) are the parity symmetric single particle wavefunctions corresponding to the densities illustrated in Fig. 2 (c), then \(|\Psi\rangle = \sqrt{N_AN_B} |\Phi_A^S(x)\rangle |\Phi_B^S(x)\rangle\) is a valid (parity symmetric) many-body wavefunction, with the diagonal of its reduced single-particle density matrix being exactly the structure shown in Fig. 2 (c). Thus, the densities of the parity symmetric solutions of the eGPEs can be interpreted as the reduced single-particle densities of the ground-state many-body wavefunction of the atomic gas and their energy can be interpreted as the energy of the many-body configuration (within the accuracy of the perturbative calculation). In contrast, parity asymmetric configurations do not permit this interpretation even if they correspond to a lower energy within the classical-field approximation, although upon proper symmetrization they could correspond to a valid (in general excited) many-body configuration (as the one shown in Fig. 7 (f)). Of course as mentioned already, the situation changes in the presence of external parity breaking contributions which would lift the degeneracy of the parity asymmetric configurations and void the need for the symmetrization, thus making the corresponding parity asymmetric configurations valid candidates for the ground-state of the atomic gas.

References↩︎