Universal quantum corrections of two-body correlation in a weakly interspecies interacting binary Bose mixture

Abstract

We investigate universal quantum corrections to two-body correlations in a zero-temperature binary Bose mixture using the Cornwall-Jackiw-Tomboulis two-particle-irreducible effective action formalism. In the weak interspecies-coupling regime, a saddle-point treatment based on Hubbard-Stratonovich transformations can be combined with a two-loop expansion and a gapless Hartree-Fock correction, thereby preserving the Goldstone theorem and reducing the coupled two-component problem to two analytically solvable single-component theories. Within this framework, we derive the ground-state energy density as a low-density expansion in the gas parameter, together with the quantum depletion and chemical potentials. The results exhibit a simple mapping to the single-component case: the universal quantum corrections of the mixture are obtained by evaluating the known single-component series at an effective scattering length $a_{\sigma\sigma}-a_{12}$ for each species, where $a_{\sigma\sigma}$ and $a_{12}$ are the intra- and interspecies $s$-wave scattering lengths. This reproduces Petrov’s equation of state at one-loop order in the weak-coupling limit and yields beyond-Lee-Huang-Yang corrections at two-loop order. We also analyze the role of mass imbalance, which enters the energy density through the exact rescaling factor $(1+m_{1}/m_{2})/2$.

I Introduction

Ultracold quantum gases offer one of the cleanest laboratories for exploring how quantum fluctuations renormalize the properties of interacting many-body systems. Even in a homogeneous, weakly interacting Bose gas at zero temperature, the ground state is not a trivial condensate: zero-point fluctuations generate a series of universal quantum corrections to the equation of state (EOS). Starting with the pioneering works of Bogoliubov 3, Lee, Huang, and Yang 12, and Wu 28, it was established that the equilibrium properties of a dilute Bose gas admit a systematic low-density expansion in powers of the gas parameter $\sqrt{na^{3}}$, where $n$ is the particle density and $a$ the $s$-wave scattering length. The leading correction to the mean-field energy density, the celebrated Lee-Huang-Yang (LHY) term proportional to $n^{2}(na^{3})^{1/2}$, together with the accompanying quantum depletion $\propto n(na^{3})^{1/2}$, constitutes the so-called two-body-correlation universal quantum effect: at this order all microscopic details of the interatomic potential beyond $a$ are irrelevant 2; 4. Beyond the LHY term, the next-order correction in the perturbative low-density expansion contains a logarithmic term of order $n^{2}(na^{3})\ln(na^{3})$ that arises from three-body correlations 28. The constant under the logarithm depends on a three-body coupling constant and was first evaluated by Braaten and Nieto 4. Beyond that order the expansion becomes nonuniversal, being sensitive to the finite range of the interatomic potential 4; 30 and, as established by Tan 26, to the three-body parameter.

A binary Bose mixture extends the physics of universal quantum corrections in a qualitatively richer direction. Following the first production of two overlapping condensates by sympathetic cooling 17, binary condensates have been realized both in hyperfine-state mixtures of the same atomic species and in heteronuclear mixtures of distinct atomic species 27. Their miscibility is governed by the celebrated criterion $g_{12}^{2}<g_{1}g_{2}$ for interspecies coupling 8. Thanks to Feshbach resonances, the interspecies $s$-wave scattering length $a_{12}$ is widely tunable, giving access to phenomena with no single-component counterpart. These include the miscible-immiscible phase transition when the interspecies repulsion exceeds the geometric mean of the intraspecies couplings 22, and most remarkably, self-bound quantum droplets: when the interspecies attraction is strong enough to render the mean-field energy negative, mechanical collapse is stabilized by the repulsive LHY quantum fluctuations 20, as confirmed experimentally in homonuclear 5 and heteronuclear 7 mixtures. The droplet regime has stimulated intense theoretical activity, including consistent many-body treatments based on bosonic pairing that remove the loophole of a complex (softened) Bogoliubov spectrum in Petrov’s theory 9. However, most existing many-body theories of mixtures either remain at the one-loop (LHY) level or focus on strongly coupled regimes, e.g., phase separation under strong repulsion 22 and droplet formation under strong attraction 20; 5; 7; 9. By contrast, the regime of weak interspecies coupling, $\absolutevalue{a_{12}}\ll\absolutevalue{a_{11}},\absolutevalue{a_{22}}$, where the mixture stays miscible and the spectrum remains well defined for both attractive and repulsive $a_{12}$, has received much less attention, especially regarding quantum corrections beyond the LHY term, and it is precisely in this regime that the universal quantum corrections of two-body correlation in a mixture should be established on a firm footing.

On the methodological side, the Cornwall-Jackiw-Tomboulis (CJT) two-particle-irreducible (2PI) effective action 6 provides a systematic framework for self-consistent quantum fluctuations, in which the Nambu-Goldstone theorem can be protected at the two-loop level by the gapless Hartree-Fock resummation scheme 11. Being a functional of the full (dressed) propagator, the 2PI effective action resums an infinite class of Feynman diagrams through the self-consistent dressing of the propagators, which is a nonperturbative partial resummation of the perturbative expansion. Its stationarity yields the Schwinger-Dyson equations. For a single-component Bose gas, this scheme has been shown to recover the full series of universal two-body-correlation corrections, i.e., the condensate fraction, energy density, chemical potential, and sound speed as expansions in the gas parameter 24, and has recently been extended to the nonuniversal EOS with finite-range interactions, yielding analytical beyond-LHY corrections that are within reach of breathing-mode spectroscopy 30. For binary mixtures, the CJT approach has been applied to finite-temperature phase transitions and symmetry-restoration phenomena 21 and to the Casimir effect of confined dual Bose-Einstein condensates (BECs) 23, yet a closed analytical two-loop EOS of a homogeneous mixture at zero temperature in the weakly coupled regime, i.e., the counterpart of the single-component low-density expansion, has remained absent. The technical bottleneck is that the interspecies interaction couples the two condensate sectors already at the level of the free propagator. The existing CJT treatments of binary mixtures 21; 23 effectively decouple the fluctuations of two species by adopting two independent $2\times 2$ propagators, thereby omitting the mixed density-fluctuation modes of a full $4\times 4$ treatment. The resulting coupled Schwinger-Dyson equations are difficult to solve analytically beyond the lowest order.

In this work, we remove this bottleneck by a different route, which, to the best of our knowledge, has not been applied to binary Bose mixtures before. We decouple the interspecies interaction via a Hubbard-Stratonovich transformation that introduces a real auxiliary field $\varphi$, with a real (imaginary) Gaussian weight for attractive (repulsive) interspecies coupling, and then treat the auxiliary field at the saddle-point level. The two species then decouple into two independent single-component problems with renormalized intraspecies couplings $\delta g_{\sigma}=g_{\sigma}-g_{12}$ and shifted chemical potentials $\mu_{\sigma}+\hbar\varphi$. Each sector can be solved separately with the full machinery of the two-loop CJT effective action combined with the gapless Hartree-Fock correction 11, after which the auxiliary field is integrated out exactly by a Gaussian integral that restores the interspecies mean-field term. This construction is valid for both signs of $g_{12}$, preserves the Goldstone theorem in each sector, and reduces the coupled self-consistent problem to tractable single-component equations. The interspecies coupling is thereby included in a mixed and partial way, entering the fluctuations of each species through the renormalized couplings $\delta g_{\sigma}$ and the chemical-potential shifts, while the dynamical mixing of the two species’ fluctuation modes is kept only at the static (saddle-point) level. The neglected dynamical part is of order $a_{12}/a$ and can be quantified by comparing with the exact one-loop results of Petrov 20. The physical content of this scheme, which is the static analog of the random-phase approximation (RPA) of condensed-matter theory, is discussed in Sec. II.

Within this scheme we study the zero-temperature condensed phase of the mixture, in which both components form BECs. The universal quantum corrections of two-body correlation considered here originate from the quantum fluctuations around the two condensates. Our main results are as follows. (i) At one-loop order, the EOS reproduces Petrov’s result 20; 9 in the symmetric limit and becomes exact as $a_{12}/a\to 0$. (ii) At two-loop order, we obtain the beyond-LHY corrections as explicit series in the gas parameter, revealing a remarkably simple structure: the universal quantum corrections of the mixture are nothing but the single-component universal series evaluated at an effective scattering length $\delta a_{\sigma}=a_{\sigma\sigma}-a_{12}$ for each species, while the direct interspecies contribution is entirely contained in the mean-field term. This “single-component mapping” encodes the two-body-correlation universal effect of the mixture in a compact and predictive way. (iii) We further analyze the effect of mass imbalance in heteronuclear mixtures, and find that, for symmetric densities and intraspecies scattering lengths, the entire energy density obeys $\mathcal{E}(n;m_{1},m_{2})=(1+m_{1}/m_{2})\mathcal{E}(n;m_{1},m_{1})/2$, an exact identity at both one- and two-loop orders.

The paper is organized as follows. In Sec. II, we introduce the model and the Hubbard-Stratonovich decoupling scheme. Section III derives the effective potential and the EOS at one-loop and two-loop levels. In Sec. IV, we present and discuss our results, including comparisons with the existing one-loop (Petrov) theory and the analysis of mass imbalance. We conclude in Sec. V.

II Model

We consider a heteronuclear Bose-Bose mixture in three dimensions, in which both components form BECs. The system is described by the partition function

$$\mathcal{Z}=\int\mathcal{D}\phi_{1}^{*}\mathcal{D}\phi_{1}\mathcal{D}\phi_{2}^{*}\mathcal{D}\phi_{2}\,e^{-S[\phi_{1}^{*},\phi_{1};\phi_{2}^{*},\phi_{2}]/\hbar}$$ (1)

in coherent state path integral formalism 19; 25, where the Euclidean action

$$S[\phi_{1}^{*},\phi_{1};\phi_{2}^{*},\phi_{2}]=\int_{x}\biggl\{ \phi_{1}^{*}(x)\left[\hbar\partial_{\tau}-\frac{\hbar^{2}\nabla^{2}}{2m_{1}}-\mu_{1}\right]\phi_{1}(x)$$

$$+\phi_{2}^{*}(x)\left[\hbar\partial_{\tau}-\frac{\hbar^{2}\nabla^{2}}{2m_{2}}-\mu_{2}\right]\phi_{2}(x)$$

$$+\frac{g_{1}}{2}\absolutevalue{\phi_{1}(x)}^{4}+\frac{g_{2}}{2}\absolutevalue{\phi_{2}(x)}^{4}$$

$$+g_{12}\absolutevalue{\phi_{1}(x)}^{2}\absolutevalue{\phi_{2}(x)}^{2}\biggr\}.$$ (2)

Here, we use the standard notations $x=(\tau,{\bf\it x})$ and $\int_{x}=\int dx=\int_{0}^{\hbar\beta}d\tau\int d^{D}{\bf\it x}$, with $D=3$ and the inverse temperature $\beta\equiv 1/(k_{\mathrm{B}}T)$. The couplings $g_{1}=4\pi\hbar^{2}a_{11}/m_{1}$ and $g_{2}=4\pi\hbar^{2}a_{22}/m_{2}$, and $g_{12}=2\pi\hbar^{2}a_{12}/m_{\mathrm{red}}$, with $a_{11}$, $a_{22}$, and $a_{12}$ being the corresponding $s$-wave scattering lengths, and with the reduced mass $m_{\mathrm{red}}=m_{1}m_{2}/(m_{1}+m_{2})$.

The interspecies interaction can be decoupled with Hubbard-Stratonovich transformations that introduce real auxiliary fields $\varphi$. When there is attraction between species, i.e., $g_{12}<0$, we employ the Gaussian integral

$$\int\mathcal{D}\varphi\exp\biggl\{\int_{x}\biggl[\frac{1}{2}\varphi(x)\hbar g_{12}^{-1}\varphi(x)$$

$$+\left(\absolutevalue{\phi_{1}(x)}^{2}+\absolutevalue{\phi_{2}(x)}^{2}\right)\varphi(x)\biggr]\biggr\}$$

$$= \exp\biggl\{-\int_{x}\biggl[g_{12}\absolutevalue{\phi_{1}(x)}^{2}\absolutevalue{\phi_{2}(x)}^{2}$$

$$+\frac{g_{12}}{2}\absolutevalue{\phi_{1}(x)}^{4}+\frac{g_{12}}{2}\absolutevalue{\phi_{2}(x)}^{4}\biggr]/\hbar\biggr\},$$ (3)

the theory then becomes

$$\mathcal{Z}=\int\mathcal{D}\phi_{1}^{*}\mathcal{D}\phi_{1}\mathcal{D}\phi_{2}^{*}\mathcal{D}\phi_{2}\mathcal{D}\varphi\,e^{-S[\phi_{1}^{*},\phi_{1};\phi_{2}^{*},\phi_{2};\varphi]/\hbar},$$ (4)

where the effective action

$$S[\phi_{1}^{*},\phi_{1};\phi_{2}^{*},\phi_{2};\varphi]=\int_{x}\biggl\{ \sum_{\sigma}\phi_{\sigma}^{*}(x)\left[\hbar\partial_{\tau}-\mu_{\sigma}\right]\phi_{\sigma}(x)$$

$$-\sum_{\sigma}\frac{\hbar^{2}}{2m_{\sigma}}\phi_{\sigma}^{*}(x)\nabla^{2}\phi_{\sigma}(x)$$

$$+\frac{\delta g_{1}}{2}\absolutevalue{\phi_{1}(x)}^{4}+\frac{\delta g_{2}}{2}\absolutevalue{\phi_{2}(x)}^{4}$$

$$-\hbar\left(\absolutevalue{\phi_{1}(x)}^{2}+\absolutevalue{\phi_{2}(x)}^{2}\right)\varphi(x)$$

$$-\frac{1}{2}\varphi(x)\hbar^{2}g_{12}^{-1}\varphi(x)\biggr\},$$ (5)

with the sums running over $\sigma=1,2$, and $\delta g_{1}=g_{1}-g_{12}$, $\delta g_{2}=g_{2}-g_{12}$. When the interspecies interaction is repulsive, i.e., $g_{12}>0$, we use the Gaussian integral

$$\int\mathcal{D}\varphi\exp\biggl\{-\int_{x}\biggl[\frac{1}{2}\varphi(x)\hbar g_{12}^{-1}\varphi(x)$$

$$+\mathrm{i}\left(\absolutevalue{\phi_{1}(x)}^{2}+\absolutevalue{\phi_{2}(x)}^{2}\right)\varphi(x)\biggr]\biggr\}$$

$$= \exp\biggl\{-\int_{x}\biggl[g_{12}\absolutevalue{\phi_{1}(x)}^{2}\absolutevalue{\phi_{2}(x)}^{2}$$

$$+\frac{g_{12}}{2}\absolutevalue{\phi_{1}(x)}^{4}+\frac{g_{12}}{2}\absolutevalue{\phi_{2}(x)}^{4}\biggr]/\hbar\biggr\},$$ (6)

which leads to the effective action

$$S[\phi_{1}^{*},\phi_{1};\phi_{2}^{*},\phi_{2};\varphi]=\int_{x}\biggl\{ \sum_{\sigma}\phi_{\sigma}^{*}(x)\left[\hbar\partial_{\tau}-\mu_{\sigma}\right]\phi_{\sigma}(x)$$

$$-\sum_{\sigma}\frac{\hbar^{2}}{2m_{\sigma}}\phi_{\sigma}^{*}(x)\nabla^{2}\phi_{\sigma}(x)$$

$$+\frac{\delta g_{1}}{2}\absolutevalue{\phi_{1}(x)}^{4}+\frac{\delta g_{2}}{2}\absolutevalue{\phi_{2}(x)}^{4}$$

$$+\mathrm{i}\hbar\left(\absolutevalue{\phi_{1}(x)}^{2}+\absolutevalue{\phi_{2}(x)}^{2}\right)\varphi(x)$$

$$+\frac{1}{2}\varphi(x)\hbar^{2}g_{12}^{-1}\varphi(x)\biggr\}.$$ (7)

Next, we adopt the saddle-point approximation and take a uniform solution $\varphi(x)=\varphi$. In this way, we can first perform integration on $\{\phi_{1}^{*},\phi_{1}\}$ and $\{\phi_{2}^{*},\phi_{2}\}$ separately. We expect this approximate treatment to be reliable, at least for weak interspecies coupling, i.e., $\absolutevalue{a_{12}}\ll\absolutevalue{a_{11}},\absolutevalue{a_{22}}$, and it makes the inclusion of higher-order quantum fluctuations more tractable.

To appreciate the advantage of this route, recall that a direct treatment of the mixture starts with the $4\times 4$ free inverse propagator $G_{0}^{-1}(k)$, whose off-diagonal blocks, proportional to $g_{12}v_{1}v_{2}$ (where $v_{1},v_{2}$ denote the condensate amplitudes of the two components), couple the density fluctuations of the two components. Existing CJT calculations for binary mixtures 21; 23 instead take two independent $2\times 2$ propagators and thereby omit these mixed fluctuation modes already at the free-propagator level. In our scheme, the interspecies coupling is absorbed into the auxiliary field: the saddle point implements its static (Hartree) level, and the dynamical part of the mixed modes enters only at order $a_{12}/a$, being systematically neglected in the weakly coupled regime.

Conceptually, this construction is the static analog of the random-phase approximation (RPA) familiar from condensed-matter many-body theory. A complete RPA treatment 18 would integrate the auxiliary field as a momentum-dependent Gaussian fluctuation, thereby resumming the density-density (ring) diagrams of the interspecies channel. Here, the saddle-point approximation freezes the auxiliary field at its uniform classical value, implementing the Hartree mean field of the interspecies channel, and the residual uniform Gaussian fluctuations, which we integrate out exactly, restore the mean field interspecies energy. This static level is justified in the weakly coupled regime, where the dynamical, momentum-dependent corrections are suppressed by the powers of $a_{12}/a$, the small parameter controlling our expansion. The scheme is non-trivial in two respects. First, it converts the coupled two-species self-consistent problem into two independent single-component problems, each of which admits an analytical two-loop solution with the Goldstone theorem protected, so that the full low-density expansion of the mixture EOS is obtained in closed form. Second, the final expressions apply uniformly to attractive and repulsive interspecies interactions: the sign of $g_{12}$ enters the Hubbard-Stratonovich decoupling only through the choice of a real or imaginary Gaussian weight, and the subsequent calculation, including the saddle-point treatment and the two-loop expansion, proceeds identically in both cases.

III Loop expansion, energy density, and quantum depletion

Within the saddle-point approximation for weak interspecies interactions, we perform self-consistent calculations by a loop expansion of the 2PI effective action, and obtain the energy density of the system when macroscopic condensation occurs at zero temperature.

III.1 One-loop level

For the case of interspecies attraction, according to Eq. (5) (the repulsive case $g_{12}>0$ follows analogously from Eq. (7), with an imaginary Gaussian weight for the auxiliary field, and yields identical final expressions), in the saddle-point approximation,

$$S\simeq \hbar\beta L^{D}\biggl[-\left(\mu_{1}+\hbar\varphi\right)v_{1}^{2}+\frac{\delta g_{1}}{2}v_{1}^{4}-\left(\mu_{2}+\hbar\varphi\right)v_{2}^{2}+\frac{\delta g_{2}}{2}v_{2}^{4}$$

$$-\frac{1}{2}\varphi\hbar^{2}g_{12}^{-1}\varphi\biggr]$$

$$+\frac{1}{2}\sideset{}{{}^{\prime}}{\sum}_{k}{\bf\it\eta}^{\top}(-k)G_{0}^{-1}(k){\bf\it\eta}(k)$$

$$+\int_{x}\biggl\{\frac{\delta g_{1}}{\sqrt{2}}v_{1}\left(\eta_{1}^{3}+\eta_{1}\eta_{2}^{2}\right)+\frac{\delta g_{1}}{8}\left(\eta_{1}^{2}+\eta_{2}^{2}\right)^{2}$$

$$+\frac{\delta g_{2}}{\sqrt{2}}v_{2}\left(\eta_{3}^{3}+\eta_{3}\eta_{4}^{2}\right)+\frac{\delta g_{2}}{8}\left(\eta_{3}^{2}+\eta_{4}^{2}\right)^{2}\biggr\}.$$ (8)

Here, we have separated the field into its coherent and fluctuating components in real field formalism, i.e., $\phi_{1}(x)=v_{1}+\left[\eta_{1}(x)+\mathrm{i}\eta_{2}(x)\right]/\sqrt{2}$, $\phi_{2}(x)=v_{2}+\left[\eta_{3}(x)+\mathrm{i}\eta_{4}(x)\right]/\sqrt{2}$, and ${\bf\it\eta}=\matrixquantity(\lx@physics@matrix\eta_{1}&\eta_{2}&\eta_{3}&\eta_{4}\endlx@physics@matrix)^{\top}$. The notation $\sideset{}{{}^{\prime}}{\sum}_{k}=\frac{1}{\hbar\beta}\sum_{\omega_{n}}\frac{1}{L^{D}}\sum_{{\bf\it k}}=\frac{1}{\hbar\beta L^{D}}\sum_{k}$ arises from the Fourier transforms, and $\omega_{n}=2\pi n/(\hbar\beta)$ are the bosonic Matsubara frequencies, $k=(\omega_{n},{\bf\it k})$. The inverse propagator

$$G_{0}^{-1}(k)=\matrixquantity(\lx@physics@matrix A_{1}&\hbar\omega_{n}&0&0\\ -\hbar\omega_{n}&A_{2}&0&0\\ 0&0&A_{3}&\hbar\omega_{n}\\ 0&0&-\hbar\omega_{n}&A_{4}\endlx@physics@matrix),$$ (9)

where

$$A_{1}= \frac{\hbar^{2}{\bf\it k}^{2}}{2m_{1}}-\mu_{1}-\hbar\varphi+3\delta g_{1}v_{1}^{2},$$ (10)

$$A_{2}= \frac{\hbar^{2}{\bf\it k}^{2}}{2m_{1}}-\mu_{1}-\hbar\varphi+\delta g_{1}v_{1}^{2},$$ (11)

$$A_{3}= \frac{\hbar^{2}{\bf\it k}^{2}}{2m_{2}}-\mu_{2}-\hbar\varphi+3\delta g_{2}v_{2}^{2},$$ (12)

$$A_{4}= \frac{\hbar^{2}{\bf\it k}^{2}}{2m_{2}}-\mu_{2}-\hbar\varphi+\delta g_{2}v_{2}^{2}.$$ (13)

The partial one-loop effective potential

$$V_{\mathrm{eff}}[v]= -\left(\mu_{1}+\hbar\varphi\right)v_{1}^{2}+\frac{\delta g_{1}}{2}v_{1}^{4}-\left(\mu_{2}+\hbar\varphi\right)v_{2}^{2}+\frac{\delta g_{2}}{2}v_{2}^{4}$$

$$+\frac{\hbar}{2}\sideset{}{{}^{\prime}}{\sum}_{k}\tr\ln(G_{0}^{-1}[v](k)/\hbar),$$ (14)

after summing over Matsubara frequencies,

$$\frac{\hbar}{2}\sideset{}{{}^{\prime}}{\sum}_{k}\tr\ln(G_{0}^{-1}/\hbar)=\frac{1}{2}\frac{1}{L^{D}}\sum_{{\bf\it k}}\left[\hbar\omega_{1{\bf\it k}}+\hbar\omega_{2{\bf\it k}}\right]$$ (15)

at zero temperature, where $\hbar\omega_{1{\bf\it k}}=\sqrt{A_{1}A_{2}}$, $\hbar\omega_{2{\bf\it k}}=\sqrt{A_{3}A_{4}}$. In the tree approximation, i.e., neglecting the quantum fluctuations, by taking the derivative of the condensate thermodynamic potential with respect to $v_{1}$ and $v_{2}$, we get

$$-\mu_{1}-\hbar\varphi+\delta g_{1}v_{1}^{2}= 0,$$ (16)

$$-\mu_{2}-\hbar\varphi+\delta g_{2}v_{2}^{2}= 0,$$ (17)

and the well-known gapless Bogoliubov spectra

$$\hbar\omega_{1{\bf\it k}}= \sqrt{\frac{\hbar^{2}{\bf\it k}^{2}}{2m_{1}}\left(\frac{\hbar^{2}{\bf\it k}^{2}}{2m_{1}}+M_{1}\right)},$$ (18)

$$\hbar\omega_{2{\bf\it k}}= \sqrt{\frac{\hbar^{2}{\bf\it k}^{2}}{2m_{2}}\left(\frac{\hbar^{2}{\bf\it k}^{2}}{2m_{2}}+M_{2}\right)},$$ (19)

where $M_{1}=2\delta g_{1}v_{1}^{2}$, $M_{2}=2\delta g_{2}v_{2}^{2}$. In the thermodynamic limit (particle number $N\to\infty$, system size $L\to\infty$, particle density $n=N/V\to\text{finite}$), $\frac{1}{L^{D}}\sum_{{\bf\it k}}\to\int\frac{d^{D}{\bf\it k}}{(2\pi)^{D}}$, after integration and regularization, we can obtain

$$V_{\mathrm{eff}}[v]= -\left(\mu_{1}+\hbar\varphi\right)v_{1}^{2}+\frac{\delta g_{1}}{2}v_{1}^{4}-\left(\mu_{2}+\hbar\varphi\right)v_{2}^{2}+\frac{\delta g_{2}}{2}v_{2}^{4}$$

$$+\frac{\sqrt{2}m_{1}^{3/2}M_{1}^{5/2}}{15\pi^{2}\hbar^{3}}+\frac{\sqrt{2}m_{2}^{3/2}M_{2}^{5/2}}{15\pi^{2}\hbar^{3}}.$$ (20)

Thus,

$$\mathcal{Z}=\int d\varphi\exp\left\{\beta L^{D}\left(\frac{1}{2}\varphi\hbar^{2}g_{12}^{-1}\varphi-V_{\mathrm{eff}}[v]\right)\right\},$$ (21)

and we can perform a Gaussian integral on the auxiliary field $\varphi$, which yields

$$\mathcal{Z}=\exp\biggl\{ -\beta L^{D}\biggl(-\mu_{1}v_{1}^{2}-\mu_{2}v_{2}^{2}+\frac{g_{1}}{2}v_{1}^{4}+\frac{g_{2}}{2}v_{2}^{4}+g_{12}v_{1}^{2}v_{2}^{2}$$

$$+\frac{\sqrt{2}m_{1}^{3/2}M_{1}^{5/2}}{15\pi^{2}\hbar^{3}}+\frac{\sqrt{2}m_{2}^{3/2}M_{2}^{5/2}}{15\pi^{2}\hbar^{3}}\biggr)\biggr\}.$$ (22)

According to the standard thermodynamic definitions, $\mathcal{Z}=\exp(-\beta L^D\veff)$ and $V_{\mathrm{eff}}=-P(\mu)$, which is related to energy density by the Legendre transformation $\mathcal{E}(n)=\mu n-P(\mu)$, $n=\partial P/\partial\mu$. This yields

$$\mathcal{E}= \frac{g_{1}}{2}n_{1}^{2}+\frac{g_{2}}{2}n_{2}^{2}+g_{12}n_{1}n_{2}$$

$$+\frac{8m_{1}^{3/2}\delta g_{1}^{5/2}}{15\pi^{2}\hbar^{3}}n_{1}^{5/2}+\frac{8m_{2}^{3/2}\delta g_{2}^{5/2}}{15\pi^{2}\hbar^{3}}n_{2}^{5/2},$$ (23)

with which we take $n_{1}=v_{1}^{2}$, $n_{2}=v_{2}^{2}$. When we focus on the idealized case of $m_{1}=m_{2}=m$, $a_{11}=a_{22}=a$ and $n_{1}=n_{2}=n/2$, we find that

$$\mathcal{E}= \frac{\pi\hbar^{2}}{m}(a+a_{12})n^{2}$$

$$+\frac{64\sqrt{2\pi}\hbar^{2}}{15m}a^{5/2}\left(1-\frac{a_{12}}{a}\right)^{5/2}n^{5/2}.$$ (24)

For comparison, Petrov’s one-loop result 9; 20 reads

$$\mathcal{E}=\frac{\pi\hbar^{2}}{m}(a+a_{12})n^{2}+\frac{32\sqrt{2\pi}\hbar^{2}}{15m}a^{5/2}\mathcal{F}(\frac{a_{12}}{a})n^{5/2},$$ (25)

where $\mathcal{F}(\alpha)=(1+\alpha)^{5/2}+(1-\alpha)^{5/2}$ becomes complex in the droplet phase $a+a_{12}<0$. In the weak-coupling regime $\absolutevalue{a_{12}/a}\ll 1$ the two expressions coincide up to corrections of order $a_{12}/a$.

III.2 Two-loop approximation

Let us now further consider the calculation of effective potential in the two-loop approximation. Since we neglect the setting sun diagrams and only include the double bubble diagrams in two-loop skeleton diagrams, the corresponding contribution, known as the Luttinger-Ward functional in the condensed-matter literature 14 (equivalently, the $\Phi$ functional of the $\Phi$-derivable scheme 2, i.e., of the CJT 2PI effective action 6), reads

$$V_{2}[G]= \frac{3\delta g_{1}}{8}\left(Q_{11}^{2}+Q_{22}^{2}\right)+\frac{\delta g_{1}}{4}Q_{11}Q_{22}$$

$$+\frac{3\delta g_{2}}{8}\left(Q_{33}^{2}+Q_{44}^{2}\right)+\frac{\delta g_{2}}{4}Q_{33}Q_{44},$$ (26)

where $Q_{ij}=\sideset{}{{}^{\prime}}{\sum}_{k}\hbar G_{ij}$, $i,j=1,2,3,4$. To protect the Nambu-Goldstone theorem when taking into account the field fluctuations, a gapless Hartree-Fock approximation introduces a phenomenological symmetry-restoring correction 10; 24; 11

$$\Delta V= -\frac{\delta g_{1}}{4}\left(Q_{11}^{2}+Q_{22}^{2}\right)+\frac{\delta g_{1}}{2}Q_{11}Q_{22}$$

$$-\frac{\delta g_{2}}{4}\left(Q_{33}^{2}+Q_{44}^{2}\right)+\frac{\delta g_{2}}{2}Q_{33}Q_{44}$$ (27)

to the $\Phi$-derivable scheme. Therefore, the effective potential in the two-loop approximation takes the form

$$V_{\mathrm{eff}}[v,G]= -\left(\mu_{1}+\hbar\varphi\right)v_{1}^{2}+\frac{\delta g_{1}}{2}v_{1}^{4}-\left(\mu_{2}+\hbar\varphi\right)v_{2}^{2}+\frac{\delta g_{2}}{2}v_{2}^{4}$$

$$+\frac{1}{2}\sideset{}{{}^{\prime}}{\sum}_{k}\tr\left[\ln G^{-1}(k)+G_{0}^{-1}[v]G(k)\right]$$

$$+V_{2}[v,G]+\Delta V,$$ (28)

which gives

$$-(\mu_{1}+\hbar\varphi)+\delta g_{1}v_{1}^{2}+\frac{3\delta g_{1}}{2}Q_{11}+\frac{\delta g_{1}}{2}Q_{22}=0,$$ (29)

$$-(\mu_{2}+\hbar\varphi)+\delta g_{2}v_{2}^{2}+\frac{3\delta g_{2}}{2}Q_{33}+\frac{\delta g_{2}}{2}Q_{44}=0,$$ (30)

by taking the variation of the effective potential with respect to $v_{1}$ and $v_{2}$. Within the 2PI formalism, the Schwinger-Dyson equations follow equivalently from the stationarity of the effective potential with respect to the full propagator, $\delta V_{\mathrm{eff}}[v,G]/\delta G=0$, which yields $G^{-1}=G_{0}^{-1}-\Sigma$ with

$$-\Sigma_{11}= \frac{\delta g_{1}}{2}Q_{11}+\frac{3\delta g_{1}}{2}Q_{22},$$ (31)

$$-\Sigma_{22}= \frac{\delta g_{1}}{2}Q_{22}+\frac{3\delta g_{1}}{2}Q_{11},$$ (32)

$$-\Sigma_{33}= \frac{\delta g_{2}}{2}Q_{33}+\frac{3\delta g_{2}}{2}Q_{44},$$ (33)

$$-\Sigma_{44}= \frac{\delta g_{2}}{2}Q_{44}+\frac{3\delta g_{2}}{2}Q_{33},$$ (34)

we can obtain the gap equations

$$-(\mu_{1}+\hbar\varphi)+\delta g_{1}v_{1}^{2}-\Sigma_{22}=0,$$ (35)

$$-(\mu_{2}+\hbar\varphi)+\delta g_{2}v_{2}^{2}-\Sigma_{44}=0,$$ (36)

and the full propagator

$$G^{-1}(k)=\matrixquantity(\lx@physics@matrix\frac{\hbar^2\vb*{k}^2}{2m_{1}}+M_{1}&\hbar\omega_{n}&0&0\\ -\hbar\omega_{n}&\frac{\hbar^2\vb*{k}^2}{2m_{1}}&0&0\\ 0&0&\frac{\hbar^2\vb*{k}^2}{2m_{2}}+M_{2}&\hbar\omega_{n}\\ 0&0&-\hbar\omega_{n}&\frac{\hbar^2\vb*{k}^2}{2m_{2}}\endlx@physics@matrix),$$ (37)

where

$$M_{1}=-(\mu_{1}+\hbar\varphi)+3\delta g_{1}v_{1}^{2}-\Sigma_{11},$$ (38)

$$M_{2}=-(\mu_{2}+\hbar\varphi)+3\delta g_{2}v_{2}^{2}-\Sigma_{33}.$$ (39)

By solving the poles of the Green’s function, we obtain the two Bogoliubov spectra

$$\hbar\omega_{1{\bf\it k}}= \sqrt{\frac{\hbar^{2}{\bf\it k}^{2}}{2m_{1}}\left(\frac{\hbar^{2}{\bf\it k}^{2}}{2m_{1}}+M_{1}\right)},$$ (40)

$$\hbar\omega_{2{\bf\it k}}= \sqrt{\frac{\hbar^{2}{\bf\it k}^{2}}{2m_{2}}\left(\frac{\hbar^{2}{\bf\it k}^{2}}{2m_{2}}+M_{2}\right)}.$$ (41)

At zero temperature, after integrating and regularizing in the continuum limit, we get

$$Q_{11}= \frac{\sqrt{2}m_{1}^{3/2}M_{1}^{3/2}}{3\pi^{2}\hbar^{3}},$$ (42)

$$Q_{22}= -\frac{\sqrt{2}m_{1}^{3/2}M_{1}^{3/2}}{6\pi^{2}\hbar^{3}},$$ (43)

$$Q_{33}= \frac{\sqrt{2}m_{2}^{3/2}M_{2}^{3/2}}{3\pi^{2}\hbar^{3}},$$ (44)

$$Q_{44}= -\frac{\sqrt{2}m_{2}^{3/2}M_{2}^{3/2}}{6\pi^{2}\hbar^{3}},$$ (45)

and

$$\frac{\hbar}{2}\sideset{}{{}^{\prime}}{\sum}_{k}\tr(G_{0}^{-1}[v]G)= \frac{1}{2}\left[\left(-(\mu_{1}+\hbar\varphi)+3\delta g_{1}v_{1}^{2}-M_{1}\right)Q_{11}+\left(-(\mu_{1}+\hbar\varphi)+\delta g_{1}v_{1}^{2}\right)Q_{22}\right]$$

$$+\frac{1}{2}\left[\left(-(\mu_{2}+\hbar\varphi)+3\delta g_{2}v_{2}^{2}-M_{2}\right)Q_{33}+\left(-(\mu_{2}+\hbar\varphi)+\delta g_{2}v_{2}^{2}\right)Q_{44}\right].$$ (46)

Starting from the partition function

$$\mathcal{Z}=\int d\varphi\exp\left\{\beta L^{D}\left(\frac{1}{2}\varphi\hbar^{2}g_{12}^{-1}\varphi-V_{\mathrm{eff}}[v,G]\right)\right\}$$ (47)

and integrating over $\varphi$, we find

$$-P(\mu)=V_{\mathrm{eff}}= -\mu_{1}v_{1}^{2}-\mu_{2}v_{2}^{2}+\frac{\delta g_{1}}{2}v_{1}^{4}+\frac{\delta g_{2}}{2}v_{2}^{4}+\frac{\sqrt{2}m_{1}^{3/2}M_{1}^{5/2}}{15\pi^{2}\hbar^{3}}+\frac{\sqrt{2}m_{2}^{3/2}M_{2}^{5/2}}{15\pi^{2}\hbar^{3}}$$

$$+\frac{\delta g_{1}}{8}\left(Q_{11}^{2}+Q_{22}^{2}\right)+\frac{3\delta g_{1}}{4}Q_{11}Q_{22}+\frac{\delta g_{2}}{8}\left(Q_{33}^{2}+Q_{44}^{2}\right)+\frac{3\delta g_{2}}{4}Q_{33}Q_{44}$$

$$+\frac{1}{2}\left[\left(-\mu_{1}+3\delta g_{1}v_{1}^{2}-M_{1}\right)Q_{11}+\left(-\mu_{1}+\delta g_{1}v_{1}^{2}\right)Q_{22}\right]$$

$$+\frac{1}{2}\left[\left(-\mu_{2}+3\delta g_{2}v_{2}^{2}-M_{2}\right)Q_{33}+\left(-\mu_{2}+\delta g_{2}v_{2}^{2}\right)Q_{44}\right]$$

$$+\frac{g_{12}}{2}\left[v_{1}^{2}+\frac{1}{2}\left(Q_{11}+Q_{22}\right)+v_{2}^{2}+\frac{1}{2}\left(Q_{33}+Q_{44}\right)\right]^{2}.$$ (48)

The particle densities are determined by

$$n_{1}= \partialderivative{P}{\mu_{1}}=v_{1}^{2}+\frac{1}{2}\left(Q_{11}+Q_{22}\right),$$ (49)

$$n_{2}= \partialderivative{P}{\mu_{2}}=v_{2}^{2}+\frac{1}{2}\left(Q_{33}+Q_{44}\right),$$ (50)

which can be combined with eqs. 31, 32, 33, 34, 35 and 36 and (38), (39), and leads to the self-consistent equations

$$M_{1}= 2\delta g_{1}n_{1}-2\delta g_{1}Q_{11},$$ (51)

$$M_{2}= 2\delta g_{2}n_{2}-2\delta g_{2}Q_{33}.$$ (52)

We can solve Eqs. (51), (52) and expand the solutions in the gas parameters $\sqrt{n_{1}(a_{11}-a_{12})^{3}}$ and $\sqrt{n_{2}(a_{22}-a_{12})^{3}}$, respectively, for a dilute (weakly interacting) mixture, which yields

$$x_{1}\equiv \sqrt{M_{1}}=\sqrt{2\delta g_{1}n_{1}}\left[1-\frac{16\sqrt{n_{1}(a_{11}-a_{12})^{3}}}{3\sqrt{\pi}}\left(1-\frac{40\sqrt{n_{1}(a_{11}-a_{12})^{3}}}{3\sqrt{\pi}}\right)+\mathcal{O}\left(\left(n_{1}(a_{11}-a_{12})^{3}\right)^{3/2}\right)\right],$$ (53)

$$x_{2}\equiv \sqrt{M_{2}}=\sqrt{2\delta g_{2}n_{2}}\left[1-\frac{16\sqrt{n_{2}(a_{22}-a_{12})^{3}}}{3\sqrt{\pi}}\left(1-\frac{40\sqrt{n_{2}(a_{22}-a_{12})^{3}}}{3\sqrt{\pi}}\right)+\mathcal{O}\left(\left(n_{2}(a_{22}-a_{12})^{3}\right)^{3/2}\right)\right].$$ (54)

Using eqs. 48, 49 and 50, we get the EOS

$$\mathcal{E}(n)= \mu n-P(\mu)$$

$$= \frac{\delta g_{1}}{2}v_{1}^{4}+\frac{\delta g_{2}}{2}v_{2}^{4}$$

$$+\frac{\sqrt{2}m_{1}^{3/2}M_{1}^{5/2}}{15\pi^{2}\hbar^{3}}+\frac{\sqrt{2}m_{2}^{3/2}M_{2}^{5/2}}{15\pi^{2}\hbar^{3}}$$

$$+\frac{\delta g_{1}}{8}\left(Q_{11}^{2}+Q_{22}^{2}\right)+\frac{3\delta g_{1}}{4}Q_{11}Q_{22}$$

$$+\frac{\delta g_{2}}{8}\left(Q_{33}^{2}+Q_{44}^{2}\right)+\frac{3\delta g_{2}}{4}Q_{33}Q_{44}$$

$$+\frac{1}{2}\left[\left(3\delta g_{1}v_{1}^{2}-M_{1}\right)Q_{11}+\left(\delta g_{1}v_{1}^{2}\right)Q_{22}\right]$$

$$+\frac{1}{2}\left[\left(3\delta g_{2}v_{2}^{2}-M_{2}\right)Q_{33}+\left(\delta g_{2}v_{2}^{2}\right)Q_{44}\right]$$

$$+\frac{g_{12}}{2}\left(n_{1}+n_{2}\right)^{2},$$ (55)

then we take $v_{1}^{2}=n_{1}-(Q_{11}+Q_{22})/2$ and $v_{2}^{2}=n_{2}-(Q_{33}+Q_{44})/2$ from Eqs. (49) and (50), so that a self-consistent universal energy density expressed in terms of the particle densities and the microscopic interaction parameters is obtained. In addition, the analytical expression of the quantum depletion

$$n_{\mathrm{ex}} =n_{1\mathrm{ex}}+n_{2\mathrm{ex}}=(n_{1}-v_{1}^{2})+(n_{2}-v_{2}^{2})$$

$$=\frac{1}{2}\left(Q_{11}+Q_{22}+Q_{33}+Q_{44}\right)$$ (56)

can be derived by substituting the solutions for $x_{1}$ and $x_{2}$.

IV Results

We first consider the balanced case, i.e., $m_{1}=m_{2}=m$, $a_{11}=a_{22}=a$, and $n_{1}=n_{2}=n/2$. In this situation the self-consistent equations derived in Sec. III admit closed analytical solutions, and the ground-state energy density takes the compact form

$$\mathcal{E}=\frac{g_{12}}{2}n^{2}+\frac{\delta gn^{2}}{4}\left[1+\frac{128}{15\sqrt{\pi}}\sqrt{\frac{n}{2}(a-a_{12})^{3}}-\frac{1024}{9\pi}\frac{n}{2}(a-a_{12})^{3}+\mathcal{O}\left(\left(\frac{n}{2}(a-a_{12})^{3}\right)^{3/2}\right)\right],$$ (57)

where $\delta g\equiv g-g_{12}=4\pi\hbar^{2}(a-a_{12})/m$ and $n=n_{1}+n_{2}$ is the total density. Equation (57) is one of the central results of this work, and its structure admits a transparent physical interpretation. The first term is the interspecies mean-field energy. The second term is precisely the sum of the two-loop EOSs of two independent single-component Bose gases, each with density $n/2$ and with the intraspecies scattering length replaced by the effective value $\delta a\equiv a-a_{12}$. Indeed, for a single-component gas the two-loop EOS reads $\mathcal{E}_{\mathrm{s}}=(g_{\mathrm{s}}\rho^{2}/2)[1+\frac{128}{15\sqrt{\pi}}\sqrt{\rho a_{\mathrm{s}}^{3}}-\frac{1024}{9\pi}\rho a_{\mathrm{s}}^{3}+\cdots]$ 24, which Eq. (57) reproduces species by species under the mapping $g_{\mathrm{s}}\to\delta g$, $\rho\to n/2$, $a_{\mathrm{s}}\to a-a_{12}$. In the limit $a_{12}/a\to 0$ Eq. (57) therefore reduces exactly to the two-loop EOS of two independent single-component gases, a nontrivial consistency check of our Hubbard-Stratonovich saddle-point scheme. The interspecies interaction enters the universal quantum corrections exclusively through the combination $\delta a=a-a_{12}$: attraction ($a_{12}<0$) strengthens the effective scattering length of each species and enhances quantum fluctuations, whereas repulsion ($a_{12}>0$) weakens them. We term this the “single-component mapping” of the two-body-correlation universal quantum effect in a weakly coupled mixture. Physically, the mapping originates from the saddle-point treatment, which converts the interspecies interaction into a static (Hartree) background for each species: it renormalizes the effective intraspecies coupling to $\delta g_{\sigma}$ but leaves the fluctuation structure of each sector unchanged. We stress that the next-to-LHY correction in Eq. (57) is a pure two-body-correlation effect: the logarithmic three-body correction $\propto n^{2}(na^{3})\ln(na^{3})$ of a single-component gas 28; 4 arises from the setting-sun skeleton diagrams, which are neglected in our two-loop approximation (only the double-bubble diagrams are retained), so that the mapping is exact within this truncation. The same conclusions hold equally for repulsive interspecies interactions ($a_{12}>0$), as demonstrated by the (b) panels of Figs. 1 and 3 and by the corresponding curve in Fig. 2.

The corresponding chemical potential is obtained from $\mu=\partial\mathcal{E}/\partial n$ evaluated along the symmetric configuration $n_{1}=n_{2}=n/2$ (equivalently, $\mu=\mu_{1}=\mu_{2}$), which yields

$$\mu=g_{12}n+\frac{\delta gn}{2}\left[1+\frac{32}{3\sqrt{\pi}}\sqrt{\frac{n}{2}(a-a_{12})^{3}}-\frac{512}{3\pi}\frac{n}{2}(a-a_{12})^{3}+\mathcal{O}\left(\left(\frac{n}{2}(a-a_{12})^{3}\right)^{3/2}\right)\right],$$ (58)

again in exact correspondence with the single-component universal series 24; 30 under the same mapping: both the leading LHY coefficient $32/(3\sqrt{\pi})$ and the next-to-LHY coefficient $512/(3\pi)$ are reproduced species by species. Likewise, substituting the solutions of the self-consistent equations into Eqs. (49) and (50), the quantum depletion fraction is found as

$$\frac{n_{\mathrm{ex}}}{n}=\frac{8}{3\sqrt{\pi}}\sqrt{\frac{n}{2}(a-a_{12})^{3}}\left[1-\frac{16}{\sqrt{\pi}}\sqrt{\frac{n}{2}(a-a_{12})^{3}}+\mathcal{O}\left(\frac{n}{2}(a-a_{12})^{3}\right)\right],$$ (59)

with $n_{\mathrm{ex}}=n_{1\mathrm{ex}}+n_{2\mathrm{ex}}$. Under the same mapping, Eq. (59) reproduces the single-component depletion $\rho_{\mathrm{ex}}/\rho=(8/(3\sqrt{\pi}))\sqrt{\rho a^{3}}[1-(16/\sqrt{\pi})\sqrt{\rho a^{3}}+\cdots]$ 24; 30, whose leading term has been verified experimentally in a homogeneous single-component gas 13. Equations (57)–(59) form a complete, parameter-free analytical description of the weakly interacting binary Bose mixture at the two-loop level: for given intraspecies scattering lengths and a given (weak) interspecies scattering length, all equilibrium thermodynamic quantities are determined to beyond-LHY order.

In Fig. 1 we compare the particle-density dependence of the energy density obtained from our one-loop and two-loop results with the prediction of Petrov’s theory 20; 9 for $a_{12}=-0.1a$ [panel (a)] and $a_{12}=0.1a$ [panel (b)]. Several observations are in order. First, the two-loop curves lie systematically below the one-loop curves, reflecting the negative beyond-LHY coefficient $-1024/(9\pi)$ in Eq. (57). At the largest density shown, $na^{3}\sim 8\times 10^{-3}$, the two-loop correction reaches about $55\%$ of the one-loop (LHY) correction for $a_{12}=-0.1a$, i.e., it is as important as in the single-component case. Second, our one-loop result agrees with Petrov’s theory in the limit $a_{12}/a\to 0$ and deviates from it by corrections of order $a_{12}/a$ for finite interspecies coupling, here $|a_{12}|/a=0.1$. This deviation is the price of the saddle-point treatment, which treats the fluctuations of each species as independent single-component fluctuations and neglects the mixed interspecies fluctuation modes. This neglect is systematically controlled by the small parameter $|a_{12}|/a$ and is precisely the regime our approach is designed for. Third, the deviation between the one-loop and two-loop results grows with density and is asymmetric between attraction and repulsion: because the effective gas parameter $\sqrt{(n/2)(a-a_{12})^{3}}$ is larger (smaller) than $\sqrt{(n/2)a^{3}}$ for $a_{12}<0$ ($a_{12}>0$), quantum corrections are enhanced by interspecies attraction and suppressed by interspecies repulsion. This asymmetry has no counterpart in a single-component gas and is a genuine mixture effect, directly traceable to the single-component mapping.

Figure 1: Energy density $\mathcal{E}$ as a function of the total particle density $n$ for (a) $a_{12}=-0.1a$ and (b) $a_{12}=0.1a$, with $m_{1}=m_{2}=m$, $a_{11}=a_{22}=a$, and $n_{1}=n_{2}=n/2$. Our one-loop and two-loop results are compared with Petrov’s one-loop theory 20. The curves are labeled in the figure. $\mathcal{E}$ is in units of $10^{-4}\hbar^{2}/(2ma^{5})$ and $n$ in units of $10^{-3}a^{-3}$.

The quantum depletion fraction is shown in Fig. 2 for the same parameters. It grows monotonically with density, consistent with Eq. (59). The next-order correction in Eq. (59) is non-negligible at the highest densities shown, as in the single-component case 24. In line with the single-component mapping, the depletion is enhanced by interspecies attraction ($a_{12}=-0.1a$, for which $\delta a=1.1a$) and suppressed by interspecies repulsion ($a_{12}=0.1a$, $\delta a=0.9a$). Since the condensate fraction is $n_{c}/n=1-n_{\mathrm{ex}}/n$, the measurement of the condensate fraction provides a direct probe of the effective-scattering-length shift $\delta a-a=-a_{12}$ induced by the interspecies interaction, without any additional fitting parameter.

Figure 2: Quantum depletion fraction $n_{\mathrm{ex}}/n$ as a function of the total density $n$ for $a_{12}=-0.1a$ and $a_{12}=0.1a$. Parameters and units are the same as in Fig. 1.

Finally, we turn to mass-imbalance effects, which are unique to heteronuclear mixtures and constitute a distinctive feature of this work. For $m_{1}\neq m_{2}$, keeping $a_{11}=a_{22}=a$ and $n_{1}=n_{2}=n/2$, the mean-field energy depends on the masses only through $g_{12}\propto 1/m_{\mathrm{red}}=1/m_{1}+1/m_{2}$ and $g_{\sigma}\propto 1/m_{\sigma}$. The fluctuation contributions of species $\sigma$ likewise scale as $m_{\sigma}^{3/2}\delta g_{\sigma}^{5/2}\propto 1/m_{\sigma}$ at one-loop order, and at two-loop order the beyond-LHY coefficients depend on the gas parameter $\sqrt{n_{\sigma}(a_{\sigma\sigma}-a_{12})^{3}}$, which is manifestly mass independent. The entire energy density therefore obeys the exact scaling identity

$$\mathcal{E}(n;m_{1},m_{2})=\frac{m}{2}\left(\frac{1}{m_{1}}+\frac{1}{m_{2}}\right)\mathcal{E}(n;m,m),$$ (60)

where $m$ is an arbitrary reference mass and $\mathcal{E}(n;m,m)$ is the mass-balanced EOS of Eq. (57). The emergence of the reduced-mass combination $1/m_{1}+1/m_{2}$ mirrors its ubiquitous role in two-body physics, where the relative motion depends on the masses only through the reduced mass $m_{\mathrm{red}}$. Equivalently, in terms of the mass ratio $r\equiv m_{2}/m_{1}$,

$$\mathcal{E}(n;r)=\frac{1+1/r}{2}\,\mathcal{E}(n;1),$$ (61)

an identity valid already at the one-loop level and persisting at the two-loop level. Although this identity is purely algebraic at the mean-field level, its persistence through the quantum fluctuations is a nontrivial structural property. At one-loop order it relies on the cancellation $m_{\sigma}^{3/2}\delta g_{\sigma}^{5/2}\propto 1/m_{\sigma}$ inside the zero-point integrals. At two-loop order, the density-fluctuation integrals $Q_{\sigma}\propto m_{\sigma}^{3/2}M_{\sigma}^{3/2}$ are themselves mass independent, because the self-consistent solution $M_{\sigma}\propto\delta g_{\sigma}n_{\sigma}$ carries the same inverse-mass factor, and every two-loop term inherits the $1/m_{\sigma}$ scaling separately. No mixed mass-ratio combination such as $\sqrt{m_{1}m_{2}}$ or $\ln(m_{1}/m_{2})$ appears, because the saddle-point decoupling renders each species an independent single-component problem whose gas parameter $\sqrt{n_{\sigma}(a_{\sigma\sigma}-a_{12})^{3}}$ is mass independent. A direct treatment based on the mixed $4\times 4$ modes would instead introduce mass-imbalance combinations already at the one-loop level. For arbitrary mass ratios, mass imbalance therefore enters the weakly coupled EOS only as a global amplitude, leaving its density dependence unchanged. Moreover, the self-consistently closed nature of the calculation, in which all quantities follow from the same solution of the Schwinger-Dyson equations and thermodynamic consistency is guaranteed by the $\Phi$-derivable structure, makes the scaling law an exact property of the scheme itself rather than an accidental feature of a particular perturbative order. Equations (60) and (61) show that, in the symmetric-density setup, mass imbalance amounts to a global rescaling of the EOS, with the lighter species dominating both the mean-field energy and the quantum fluctuations, as both the coupling $g_{\sigma}$ and the zero-point energy of the Bogoliubov modes scale as $1/m_{\sigma}$ at fixed scattering lengths. In Fig. 3 we verify this scaling for the mass ratios $r=2,1,1/2,1/3$: the numerical curves collapse onto Eq. (61) for both signs of $a_{12}$.

For realistic heteronuclear mixtures, e.g., ${}^{41}\mathrm{K}$–${}^{87}\mathrm{Rb}$ with $r\simeq 87/41\simeq 2.12$, the scaling factor is $\simeq 0.74$, i.e., a $\sim 26\%$ reduction of the energy density with respect to the mass-balanced EOS with the same scattering lengths and densities, a sizable effect that provides a clear benchmark for future numerical studies of heteronuclear mixtures.

Figure 3: Energy density as a function of the total density $n$ for different mass ratios $r=m_{2}/m_{1}=2,1,1/2$, and $1/3$ (as labeled in the figure), with (a) $a_{12}=-0.1a$ and (b) $a_{12}=0.1a$. Other parameters are as in Fig. 1. The curves confirm the scaling identity (61): $\mathcal{E}(r)=(1+1/r)\mathcal{E}(1)/2$. $\mathcal{E}$ is in units of $10^{-4}\hbar^{2}/(2m_{1}a^{5})$ and $n$ in units of $10^{-3}a^{-3}$.

Before concluding, we comment on the range of applicability of our results. The parameters of Figs. 1–3, $a_{12}=\pm 0.1a$ and $na^{3}\lesssim 8\times 10^{-3}$, lie deep inside the weak-coupling regime $\absolutevalue{a_{12}}\ll a$ and the dilute regime $na^{3}\ll 1$, where both the saddle-point treatment and the low-density expansion are quantitatively controlled. Two remarks are in order. First, the dilute limit enters only through the analytic expansion of the self-consistent solutions: the two-loop CJT scheme itself is a nonperturbative $\Phi$-derivable resummation, with dressed propagators and the Goldstone theorem preserved by the gapless Hartree-Fock correction, so that the same self-consistent equations can be solved numerically at substantially larger gas parameters. Second, the restriction $\absolutevalue{a_{12}}\ll a$ is tied to the saddle-point approximation rather than to the $\Phi$-derivable framework itself, and stronger interspecies couplings can be addressed by including the fluctuations of the auxiliary field. Within the analytic expansion, the expressions remain formally well defined for all $\absolutevalue{a_{12}}<a$ (i.e., $\delta g_{\sigma}>0$), away from the droplet boundary $a_{12}\to-a$ and the miscibility boundary $a_{12}\to a$. The accuracy of the saddle-point treatment is quantified at leading order by comparing our one-loop EOS with Petrov’s exact one-loop result: the relative deviation of the LHY coefficient, $|2(1-\alpha)^{5/2}-\mathcal{F}(\alpha)|/\mathcal{F}(\alpha)$ with $\alpha=a_{12}/a$, grows from about $25\%$ at $|\alpha|=0.1$ to about $65\%$ at $|\alpha|=0.3$. Since the LHY term contributes only a fraction of the total energy in the dilute regime, this translates into an error of a few percent in the total energy density at $|\alpha|=0.1$ and of order ten percent at $|\alpha|=0.3$. We therefore expect quantitative control for $\absolutevalue{a_{12}}/a\lesssim 0.1$, as used in the figures, and semiquantitative usefulness up to $\absolutevalue{a_{12}}/a\sim 0.3$. Within this range the predicted two-loop correction, which reaches about $55\%$ of the LHY term at $na^{3}\sim 8\times 10^{-3}$, exceeds the estimated systematic error of the saddle-point treatment, so that the beyond-LHY signal is not masked by the approximation error.

V Conclusion

In summary, we have developed a unified and fully analytical theory of the universal quantum corrections to two-body correlations in a weakly interspecies interacting binary Bose mixture at zero temperature. The central strategy is to decouple the interspecies interaction via a Hubbard–Stratonovich transformation and to treat the resulting auxiliary field at the saddle-point level. This procedure reduces the binary mixture to two independent single-component problems with renormalized intraspecies couplings $\delta g_{\sigma}=g_{\sigma}-g_{12}$, which are solved within the two-loop Cornwall–Jackiw–Tomboulis effective action supplemented by the gapless Hartree–Fock correction. The resulting $\Phi$-derivable scheme, formulated in terms of dressed propagators, is self-consistent beyond the plain Bogoliubov expansion and preserves a gapless excitation spectrum. Integrating out the auxiliary field then restores the interspecies mean-field contribution. The approach applies equally to attractive and repulsive interspecies interactions, protects the Goldstone theorem in each sector, and renders the coupled self-consistent problem analytically tractable.

Within this framework, we derived the ground-state energy density, the chemical potential, and the quantum depletion as explicit low-density expansions in powers of the gas parameter. The results reveal a remarkably simple structure, which we call the single-component mapping: the universal quantum corrections of the mixture are obtained by evaluating the single-component universal series at an effective scattering length $\delta a_{\sigma}=a_{\sigma\sigma}-a_{12}$ for each species, while the direct interspecies contribution is entirely contained in the mean-field term. At one-loop order, the equation of state reproduces Petrov’s result 20; 9 in the limit $a_{12}/a\to 0$, and at two-loop order it yields the beyond-LHY corrections, which are enhanced by interspecies attraction and suppressed by interspecies repulsion. For heteronuclear mixtures, we further established an exact scaling identity for the energy density under mass imbalance at the one-loop level, which survives the two-loop truncation and captures the mass-ratio dependence of the equation of state through the simple factor $(1+m_{1}/m_{2})/2$. This result is directly relevant to realistic mixtures such as ⁴¹K–⁸⁷Rb.

Our results are universal in that they depend only on the $s$-wave scattering lengths, and they clarify how interspecies coupling renormalizes the two-body-correlation effects of each species. The analytical equation of state derived here also provides a controlled starting point for addressing nonuniversal finite-range effects and three-body correlations. More broadly, multi-component Bose gases continue to attract significant interest, as exemplified by quantum droplets in three-component mixtures 16; 15. For a fully symmetric $N$-component mixture, in which all species share the same mass, intraspecies coupling, and interspecies coupling, the present scheme extends directly: the single-component mapping remains valid, and the energy density is obtained from Eq. (57) by replacing the per-species density $n/2$ with $n/N$, namely, $\mathcal{E}=g_{12}n^{2}/2+(\delta g\,n^{2}/2N)[1+\frac{128}{15\sqrt{\pi}}\sqrt{(n/N)(a-a_{12})^{3}}-\frac{1024}{9\pi}(n/N)(a-a_{12})^{3}+\cdots]$. The chemical potential and the depletion follow from the same substitution. Such closed analytical results for the weakly coupled multi-component case, which follow directly from the single-component mapping, provide useful benchmarks for future studies of multi-component Bose gases. For non-identical species, by contrast, the Hubbard–Stratonovich decoupling requires a matrix of auxiliary fields and the analysis becomes considerably more involved. We expect that the single-component mapping and the mass-ratio scaling identity derived here will serve as useful benchmarks for experiments and numerical simulations of weakly coupled binary Bose mixtures.

*