Phase separation in binary Bose mixtures at finite temperature

We investigate the magnetic behavior of finite-temperature repulsive two-component Bose mixtures by means of exact path-integral Monte-Carlo simulations. Novel algorithms are implemented for the free energy and the chemical potential of the two components. Results on the magnetic susceptibility suggest that the conditions for phase separation are not modified from the zero temperature case. This contradicts previous predictions based on approximate theories. We also determine the temperature dependence of the chemical potential and the contact parameters for experimentally relevant balanced mixtures.


Introduction
The realization of mixtures of ultracold gases in the quantum-degenerate regime has opened new interesting directions to study the simultaneous presence of superfluidity in multicomponent systems, which could not be addressed with traditional quantum fluids such as liquid 4 He or standard superconductors.The first examples of superfluid mixtures have been produced both with atoms obeying the same [1,2] or different statistics [3].In particular, two-component mixtures of bosonic species below the Bose-Einstein transition temperature provide one with the simplest set up to investigate the interplay between quantum magnetism and superfluid properties.This includes novel phenomena such as combined mass and spin superfluidity [4], non dissipative spin drag [5], and Bose-enhanced magnetic effects [6].In the case of repulsive mixtures, the zero temperature scenario is well described by mean-field theory [7]: the ground state is paramagnetic if the interspecies coupling constant is below a threshold set by the strength of interactions within each component, and is instead fully ferromagnetic, i.e. full phase separation between the two components occurs, if the coupling exceeds this critical value.This scenario has been also confirmed in a series of experiments [8][9][10][11] and quantum Monte Carlo simulations for trapped mixtures, both at zero [12] and finite temperature [13].At finite temperatures, perturbative approaches, such as Hartree-Fock (HF) and Popov theories, predict an intriguing scenario holding for mixtures below the Bose-Einstein condensation (BEC) temperature: The paramagnetic state at low temperature can turn ferromagnetic at higher temperature if the interspecies coupling is close enough to the T = 0 threshold [14][15][16].According to these theoretical schemes, the mechanism responsible for the magnetic transition are beyond mean-field effects induced by temperature, which destabilize the paramagnetic phase.Similar effects of pure quantum nature have instead a stabilizing role in attractive mixtures and lead to the formation of self-bound droplets [17][18][19].An important question, which needs to be answered, is whether the predictions of perturbative approaches are accurate enough to include the relevant role played by fluctuations around the transition temperature.
In this work we use exact path-integral Monte Carlo (PIMC) simulations to investigate the magnetic and thermodynamic properties of a repulsive two-component Bose mixture.In particular, novel algorithms are implemented to obtain precise unbiased predictions for the chemical potentials of the two separate components and for the total free energy.This provides us with crucial information on the chemical equilibrium at finite polarization and on the occurrence of stable free energy minima.We find that the magnetic susceptibility at finite temperature is well described by the simple zero temperature mean-field prediction and also that there is no evidence of a temperature-induced ferromagnetic transition.Consequently, the conditions for phase separation remain unchanged from the T = 0 case.Furthermore, for the choice of interspecies coupling corresponding to a balanced mixture of sodium atoms, we calculate chemical potential and contact parameters as a function of temperature, pointing out their deviations in the critical region from the predictions of perturbative methods.In particular, the interspecies contact parameter features a suppression at intermediate temperatures caused by statistical effects which indicates an enhanced repulsive correlation between the two components.

Methods
We consider the following Hamiltonian describing a system of N = N 1 + N 2 Bose particles belonging to two distinguishable components with equal mass m The intraspecies potentials are assumed to be the same, denoted by v(r), and v 12 (r) describes interspecies interactions.All potentials are repulsive and modeled by hard spheres, i.e. the potential is infinite inside the diameter of the sphere and zero outside.The two parameters a and a 12 define, respectively, the range of the v and v 12 potential and the corresponding value of the s-wave scattering length.In the dilute regime of interest, interaction effects only depend on the coupling strengths: g = 4πℏ 2 a m and g 12 = 4πℏ 2 a 12 m .To discuss magnetic properties we introduce the component densities n 1 + n 2 = n and the polarization parameter p = (n 1 − n 2 )/n.For such a symmetric mixture, mean-field theory at zero temperature predicts miscibility (p = 0) if g 12 < g and a fully separated state (p = 1) if g 12 > g [7].Furthermore, the same theory yields the expression χ = 2 g−g 12 for the magnetic susceptibility in the paramagnetic phase.
In a PIMC simulation, we use periodic boundary conditions in a box of volume V at fixed density n = N/V .We work with the well established worm algorithm in continuous space to efficiently sample bosonic permutations [20].Recently the method has been further developed to be fully consistent with periodic boundary conditions and was applied to the study of the single-component gas [21].The algorithm is described in details in Ref. [22].
In the present study, we implemented also the calculation of the total free energy, of the free energy differences for different polarizations, and of the chemical potentials for both components in the canonical ensemble, generalizing to mixtures the technique first proposed in Ref. [23].The details of the Monte Carlo moves added to the PIMC algorithm can be found in the Appendices A and B. In addition, we use also HF and Popov theories to compare with PIMC results.Details on the derivation of the free energy and related quantities within the HF and Popov scheme are given in Appendix C.

Magnetic behavior of binary mixtures
We first focus on the magnetic properties of the mixture, analyzing how the chemical potential and the total free energy depend on the polarization.We choose the value na 3 = 10 −4 for the gas parameter which on one side emphasizes the interesting effects due to interactions and on the other side ensures that the results are universal in terms of solely the gas parameter.However it is worth pointing out that stronger interactions could be realized in resonantly interacting gases [24].
In Fig. 1, the chemical potentials of the two components are plotted against polarization at fixed temperature T = 0.794T 0 c , where k B T 0 c = 2πℏ 2 m (n/2ζ(3/2)) 2/3 is our reference energy scale, corresponding to the BEC transition temperature of a balanced (p = 0) non-interacting mixture.The majority component 1 is Bose condensed for all values of p shown in the figure, while, at this temperature, the minority component 2 turns normal at the critical polarization p c ≃ 1 − (T /T 0 c ) 3/2 ≃ 0.292 corresponding to the maximum of the HF and Popov results for µ 2 .In panels (b) and (c), referring to g 12 > 0, we notice that HF and Popov theories predict a crossing of chemical potentials at finite polarization p > p c .This crossing corresponds to a minimum in the free energy F according to the thermodynamic relation µ  is Bose condensed and in equilibrium with the minority one in the normal phase [15].This behavior of the HF and Popov free energies is shown in Fig. 2 [see panels (b) and (c)].Note that the chosen value of g 12 /g = 0.93 corresponds to the |F = 1, m F = 1⟩ and |F = 1, m F = −1⟩ Bose-Bose mixture of 23 Na atoms investigated experimentally in Refs.[25,4].According to HF and Popov theories, this mixture should provide an example of the striking phenomenon of a paramagnetic state at low temperature which turns ferromagnetic at higher temperatures, as predicted in Ref. [15].However, the PIMC results for µ 1 and µ 2 at g 12 /g = 0.93, do not confirm this scenario.The majority component chemical potential µ 1 is in good agreement with the Popov result, but µ 2 deviates significantly in the region p > p c and does not exhibit the peak predicted by HF and Popov theories.As a result, no crossing occurs for p > p c and no minimum appears in F (p) other than at p = 0. Furthermore, from the thermodynamic relation χ holding at small polarization, we find a good agreement using the zero temperature mean-field result χ = 2/(g − g 12 ) of the magnetic susceptibility, as can be seen in panels (a) and (b) of Fig. 2, where the MF prediction, shifted to coincide with the PIMC data at p = 0, well reproduces the p 2 behavior of the PIMC data.Similar results are obtained for the fully symmetric case g 12 = g, where the chemical potentials exactly coincide for p < p c and separate without crossing for larger polarizations.As a result the free energy is flat as a function of polarization and the magnetic susceptibility diverges.Interestingly, this behavior, which is well understood at T = 0 where the ground state is degenerate with respect to polarization, remains valid at finite temperature as long as both condensates are present.The results shown in panels  interaction between the two components: µ 2 decreases with p, although without a small peak, and F monotonically increases following the mean-field magnetic susceptibility.More interesting is the case g 12 = 1.2g,where the mixture is phase separated already at T = 0. Notice that Popov theory can not be applied here if both condensates are present because spin excitations acquire an unphysical complex energy.The minority chemical potential µ 2 displays a maximum, although not as large as predicted by HF theory, and a crossing point with µ 1 .As a consequence, the free energy indicates instability at p = 0 and shows a clear minimum at p > p c , corresponding to the phase separated state with partial polarization.
We further analyze the magnetic behavior of the mixture in Fig. 3 where we show the free energy difference ∆F = F (p) − F (0) as a function of p 2 for the intermediate value g 12 = 0.5g of the interspecies coupling constant and at T = 0.794T 0 c .This choice of parameters and, in particular, the choice of temperature emphasizes thermal effects in HF and Popov theories yielding important corrections to the T = 0 magnetic susceptibility.We also note that finitesize effects in PIMC simulations of the free energy are negligible if one increases further the total number of particles.We find that F depends linearly on p 2 over a large range of values extending also beyond the critical polarization p c .Furthermore, the coefficient of the linear dependence, proportional to χ −1 , is well reproduced by the mean-field result χ −1 = (g −g 12 )/2 shown in the figure by the MF line.In contrast, HF and Popov results provide a poor account of the polarization dependence of the free energy.A possible explanation of this inadequacy involves the role of critical fluctuations which control the thermodynamics close to the transition point and, in general, can not be described using perturbative methods such as HF and Popov theories.The width of the critical region is predicted to shrink as na but for experimentally relevant values of the gas parameter (na 3 ≃ 10 −4 − 10 −6 ) it remains of the same order as the transition temperature itself.
From these results we conclude that, in contrast to HF and Popov predictions, the magnetic susceptibility depends very little on the temperature, and the conditions for phase separation seem to remain the same as at T = 0.In fact, if g 12 < g, our results indicate that the only thermodynamically stable phase is the paramagnetic state at p = 0.A ferromagnetic state forms when g 12 > g and the effect of temperature is to reduce the equilibrium polarization from the p = 1 value achieved only at zero temperature.This is found at a high temperature not far from the BEC transition point and we expect the same to be true also for lower temperatures, where thermal effects not captured by the mean-field description should play a minor role.In this respect one should also notice that higher order interaction effects at T = 0 do not change the critical value g 12 = g for the onset of ferromagnetism (see Ref. [16]).As an additional remark, we point out that our results do not exclude a non trivial interplay between ferromagnetic and critical fluctuations in the close vicinity of the transition point.To carefully investigate these effects would require a much deeper analysis of the shift of the transition point in interacting mixtures beyond the scope of this work.Furthermore, we expect the simple T = 0 scenario to hold also at densities lower than na 3 = 10 −4 .Numerical checks show that for vanishing gas parameter the free energy difference between the p = 0 state and the stable minimum at finite p predicted by Popov theory is suppressed as g 3/2 and furthermore the minimum is shifted towards higher temperatures occurring closer to the transition point.As a consequence, we expect critical fluctuations to play a major role in the magnetic response of the mixture also in the regime of extremely low densities, thereby invalidating the predictions of Popov theory.

Particle-position snapshots
Visualizing instantaneous particle positions during PIMC simulations allows us to shed some light on the ferromagnetic transition.To minimize the effects due to inter-domain interfaces, we consider large scale simulations comprising N = 8000 particles, with N 1 = N 2 .The gas parameter is na 3 = 10 −4 .Fig. 4 shows the position snapshots observed after thermalization is reached.Two initial particle configurations are considered.They feature either vertically separated or mixed components.In the separated configuration, the first component is uni-  formly randomly distributed in the lower half of the 3D simulation box, while the second component is in the upper half.In the mixed initial configuration, both components are uniformly distributed in the whole box.In panels (a) and (b), the interspecies coupling strength is g 12 /g = 0.93, i.e., below the T = 0 MF critical point.In the first panel, the temperature is relatively low, namely, T /T 0 c ∼ = 0.238.Here, even Popov theory would predict a paramagnetic state.In the second, it is closer to the BEC transition temperature, where Popov theory would predict a ferromagnetic state.The two simulations start in the separated configuration.Despite of being initially separated, the two components rapidly mix, indicating a paramagnetic state, both at low temperature and closer to the BEC transition.In panel (c), the inter-species interaction strength (g 12 /g = 1.2) is beyond the critical point predicted by the MF theory.In this case, the two components keep the initial spatial separation along the vertical direction, with only minor mixing close to the interface separating the two domains.Interestingly, even when they start from a mixed configuration [panel (d)], they form two well defined ferromagnetic domains.This indicates that large-scale PIMC simulations are able to simulate phase separated states.Chiefly, these observations further corroborate the claim that the finite temperature transition corresponds to the T = 0 MF scenario, in contrast to the HF and Popov predictions.When the temperature is raised [panel (e)], the interface is less regular and it looses memory of the initial position.One also notices a larger impurity density, corresponding to a ferromagnetic state with partial polarization.Moving even closer to the BEC transition temperature T 0 c [panel (f)], the two domains can be hardly identified by naked eye.However, we argue that in the thermodynamic limit one would still observe a (partially) ferromagnetic state, meaning that the Curie critical temperature where melting occurs is even higher.

Thermodynamic properties of balanced mixtures
We now turn our attention to the study of thermodynamic quantities, focusing on the g 12 = 0.93g sodium mixture in the balanced state p = 0.In this case we have chosen the value na 3 = 10 −6 for the gas parameter which is closer to experimentally relevant conditions in the absence of Feshbach resonances.The PIMC results for the thermodynamic quantities shown below are the extrapolations to the thermodynamic limit of the data computed with up to 512 total particles.In Fig. 5 we show the chemical potential µ = µ 1 = µ 2 of the mixture as a function of temperature below and above the transition point and we compare it with the results of perturbative approaches.The results are in good agreement with both HF and Popov predictions when the temperature is not too close to the critical point.In the critical region around T 0 c , deviations are sizable.They tend to suppress the maximum, similarly to the results of Fig. 1 for the minority component.We notice that a maximum in the temperature dependence of the chemical potential should be expected on general grounds from the theory of superfluids and has been recently observed in a single-component dilute Bose gas [27].The PIMC results for µ in a single-component gas are discussed in the appendix as a test study of the chemical potential algorithm.
In Fig. 6 we show the results for the contact parameters, important thermodynamic quantities sensitive to short-range correlations.In a symmetric unpolarized mixture one defines two contact parameters C 11 = C 22 = C and C 12 associated to correlations within each component and between the two components respectively The contact parameter C has been measured as a function of temperature in a single-  component gas [28] and in a mixture of a Bose gas with impurities [29].In our PIMC simulations we have computed C and C 12 from the short-range behavior of the pair correlation function for particles belonging to the same and to different components [21].The results for C are in good agreement with both HF and Popov predictions, showing deviations only in the vicinity of T 0 c .For C 12 , instead, the HF prediction does not depend on the temperature, while the Popov prediction yields a small minimum.Our PIMC findings show a small minimum around T ≃ 0.7T 0 c which reproduces this.This enhanced repulsive correlation between the two components at intermediate temperatures has been already discussed in repulsive mixtures [30,31] and deserves further investigation.

Conclusion
We have investigated the magnetic and thermodynamic properties of repulsive Bose mixture using exact numerical methods.For the values of the parameters considered in the simulations we do not find the ferromagnetic transition predicted to occur at finite temperature by perturbative approaches and we find good agreement with the magnetic susceptibility from simple mean-field theory at zero temperature.We further argue that a similar conclusion is expected to hold for lower values of the gas parameter.This claim is further corroborated by the analysis of particle-positions snapshots.Thermodynamic quantities reveal the role of critical fluctuations close to the BEC transition point and the behavior of the contact parameters contains important information on short-range correlations in the mixture that can be measured in future experiments.Our findings indicate the importance of unbiased simulations for atomic mixtures, in contrast to previous perturbative treatments of repulsive and attractive two-component Bose gases.

A PIMC computation of chemical potential and free energy
In this appendix, we present the details of the PIMC algorithm we employ for the computation of the chemical potential of a Bose gas.The basic idea is to recognize that the chemical potential can be derived from the ratio of the partition functions for the systems with N + 1 and N particles (at fixed volume and temperature) as As noted in Ref. [23] the above formula can be leveraged in a canonical PIMC calculation by enlarging the configurational space to include the sector with one additional particle.The ratio Z N +1 /Z N is then evaluated as the relative time spent by the simulation in the two sectors.The simulation resembles a grand canonical one, with the difference that it is restricted to states with either N or N + 1 particles, thus providing higher statistics for the computation of µ(N, T ).Combining the chemical potential with the pressure, we can obtain the free energy where Ω = −P V is the grand canonical potential.

A.1 Details of the algorithm
In order to extend the algorithm described in Ref. [22] and enable the computation of the chemical potential we work with N + 1 polymers and implement a boolean variable for each polymer to activate or deactivate it.The configurational space is now composed by four sectors: the original Z N and G N together with the corresponding sectors with one additional particle Z N +1 and G N +1 and one needs to introduce appropriate Monte Carlo moves to allow the Markov-Chain to visit all the configurations within these sectors.The four sectors together with the sector-changing moves are summarized in Fig. 7.In general one can introduce a grand canonical chemical potential µ gc as a simulation parameter to be used to increase the sampling efficiency.In particular it The complementary move consists in removing a one-polymer ring with zero winding from the Z N +1 sector by deactivating it.The move is accepted with probability

A.2 Benchmarks
Following the strategy of Ref. [22] we carefully check our implementation by running a number of tests.First of all, we verify that we correctly recover the values of the chemical potential for the ideal Bose gas for each system size N .In Fig. 8 we show the PIMC results at the temperature T = 1.5T 0 c compared with the exact values (obtained via the recursion formulas as in Refs.[32,33] and reviewed in Ref. [22]) and with the result in the thermodynamic limit given by where z is an effective fugacity that determines the total density of the gas via the equation nλ 3 T = g 3/2 (z) with g ν (z) the usual special Bose functions.The agreement between the PIMC data and the expected values is perfect at any size and does not depend on the number of imaginary-time slices used in the simulation.Moreover we verify that below T 0 c the PIMC results are compatible with a zero chemical potential.
We then benchmark the interacting gas, where the repulsive interaction is modeled by a hard sphere potential.As in Ref. [21], we use the pair-product ansatz [34] for the computation of the potential energy ∆U .In Fig. 9 we compare the PIMC results for the chemical potential and the free energy with the perturbative predictions.The PIMC data are extrapolated to the thermodynamic limit using a linear fit in 1/N of the results for four sizes N = 128, 256, 384, 512.The number of imaginary-time slices is 16 for all sizes.For the chemical potential, we also compare our results with the predictions from the universal relations of Ref. [26], using their data for the density shift λ(X) ∝ n − n c to extract the reduced temperature t = T /T 0 c and then mapping it to the corresponding chemical potential using the data for the chemical potential shift X ∝ µ − µ c .In particular, expressing the relations in our units, we find that each value of λ(X) can be mapped to a value of t solving the equation The PIMC values, extrapolated to the thermodynamic limit (black crosses), are compared with the perturbative results of Hartree-Fock (green dashed line) and Popov (red solid line) theories.Left panel: Results for the chemical potential, also compared with the predictions from the universal relations of Ref. [26] (blue dots).Right panel: difference in freee energy with the ideal Bose gas result.
in the region where the universal relations can be applied, namely t ∼ 1.The numerical constant C is determined as C = 0.0142(4).From the corresponding values of X we then determine the chemical potential shift as Finally, to get the sought-after values of µ, we need to add the values of µ c as obtained from Ref. [35] where K = 0.673(1) is a numerical constant1 .The data for µ obtained from the universal are represented by the blue dots in the left panel of Fig. 9 and show a good agreement with the PIMC data for T < T 0 c , while, for T > T 0 c , a discrepancy builds up for increasing temperatures, since the universal relations are valid only in the regime of large occupation numbers for singleparticle modes.In that regime, the PIMC data nicely reproduce the HF predictions.With the benchmarks shown so far, we are now confident that the PIMC algorithm correctly reproduces the physics of a single-component Bose gas both in the non-interacting and in the interacting case.In the following section we show how to extend the algorithm for Bose mixtures.

B PIMC algorithm for a binary Bose mixture
Extending the PIMC algorithm to the case of multicomponent gases is pretty straightforward: one just needs to restrict the Swap move to involve only particles of the same species and, in the interacting case, to take into account the inter-species interaction described by the s-wave scattering length a 12 .The computation of the chemical potentials for the two species proceeds as before via the free energy difference, this time making sure the number of particles of the other species is kept fixed.For a two-component mixtures we have: where the differences are numerically computed as the ratios of Monte Carlo times spent in the different sectors.The simulation now lives in a configurational space made by 4 × 4 = 16 sectors.Several consistency checks where made on the algorithm.In Fig. 10 we show the results for two non-interacting ideal Bose gases, where we recover the known exact result as a function of the polarization.The chemical potential µ 1 for the majority component is consistent with zero, while the chemical potential µ 2 for the minority component becomes non-zero above the critical polarization, where it becomes normal.
Using the above method for computing the chemical potentials, one can obtain the value of the free energy of the mixture via the thermodynamic relation Notice that, while this quantity contains valuable information and it represents our main test-bench for the perturbative predictions, it comes at the cost of a cancellation between the pressure term and the chemical potential terms, that amplifies its final statistical error.However, when we focus on the magnetic properties of a binary mixture, we are only interested in the free energy difference among mixtures at different values of the polarization p = (N 1 − N 2 )/N , where N = N 1 +N 2 is the total number of particles.Such a difference can be evaluated more efficiently by devising an algorithm that directly samples configurations with different values of the polarization, while keeping N fixed.Denoting with Z N,p the partition function with N total particles and polarization p, the free energy difference ∆F (N, p) between the state at polarization p and the unpolarized state with p = 0 can be computed as where the ratio t(Z N,P )/t(Z N,0 ) is the ratio between the time spent in the sector with polarization p and the time spent in the sector with zero polarization.There are many possible ways to implement such an algorithm: One possibility is, for example, to combine the moves of Sec.A.1 for the two species in such a way that each time a particle of one species is created a particle of the other species is removed.In the following we mention another possibility, which is slightly more sophisticated.

B.1 Details of the algorithm for free energy differences in a mixture
An efficient algorithm that spans the configurations with different polarizations, while keeping N fixed can be devised by taking close inspiration from the original grand canonical implementation of Refs.[20,36].The Monte Carlo moves have been adapted in such a way that both the total number of polymers and the total number of beads are kept constant throughout the simulation.Within this algorithm, the worms for the two species might be present simultaneously, also in a configuration where the beads of one polymer are shared between the two worms.For example, the polymer i 0 might be filled by the species 1 up to the imaginary-time slice j 0 (corresponding to the head of the worm 1), while the rest of the slices are filled by the worm of the species 2 (that has its tail at the imaginary-time slice j 0 + 1 of the polymer i 0 ).We briefly describe below a minimal pair of moves, called Advance/Recede, that allows the simulation to span the configurations at different values of polarization (except the case at p = 1).Other moves can be included in order to improve the ergodicity of the Markov chain across the sectors, for example by combining Advance/Recede with Open/Close.The details of these combined moves will not be given here; instead we just outline the aforementioned minimal addition that can be used for small values of the polarization.
The Advance and Recede moves change the relative lengths of the worms and can only be performed when both worms are present.In the Advance move the head of the worm of species s is advanced in imaginary time from the slice j to the slice j + ∆j, by sampling the new ∆j beads with the staging algorithm as in Move Head.The tail of the worm of the other species s ′ is advanced as well by deleting all the beads between the slice j and the slice j + ∆j.Note that we must reject the move if the worm of species s ′ is completely deleted by the proposed update.The move is then accepted with probability A advance = min 1, e β∆µ ∆j/M −∆U , (B.17) where ∆µ = µ s gc − µ s ′ gc is the difference between the grand canonical chemical potentials of the two species.The complementary Recede move is completely symmetric and can be obtained as the Advance move with negative values of ∆j.It consists in receding the head of the species s by deleting ∆j beads, while simultaneously creating ∆j new beads for the species s ′ .The move is accepted with probability This pair of moves can change the species of whole polymers, thus allowing the algorithm to sample configurations with different polarizations.We have checked that the free energy differences computed directly through Eq. (B.16) reproduce those obtained from the full computation of the free energy, but deliver smaller statistical errors.

C Hartree-Fock and Popov theories
The Hartree-Fock and Popov theories of repulsive binary Bose mixtures at finite temperature are described in details in Refs.[15,16].We note that Popov's theory is also known as the finite temperature extension of Beliaev's approach and includes the important contribution of anomalous fluctuations to thermodynamic quantities [37,38].Here we report the results for the Helmholtz free energy obtained in the two approaches from which all thermodynamic quantities discussed in the main text can be derived.
Within the HF approximation one finds where n 0 = n − 2n 0 T is the condensate density calculated to lowest order in the interaction strength.In the regime of high polarization (p > p c ) the above expression reduces to

Figure 1 :
Figure 1: Chemical potentials µ 1 (blue) and µ 2 (red) as a function of the polarization p = (N 1 −N 2 )/N for a system with a total of N = N 1 + N 2 = 128 particles at temperature T = 0.794T 0 c and with gas parameter na 3 = 10 −4 , for four values of the couplings ratio g 12 /g.The dashed lines are the HF predictions, the solid lines are the Popov predictions.For the minority (red) component, the two coincide for p > p c .In panel (d) only the HF lines are shown.The vertical lines indicate the critical polarization p c .
(a) and (d) of Figs. 1 and 2 are instead in better qualitative agreement with approximate perturbative approaches.The case g 12 = 0 corresponds to no

Figure 2 :
Figure 2: Free energy as a function of the polarization p = (N 1 − N 2 )/N for a system with a total of N = N 1 + N 2 = 128 particles at temperature T = 0.794T 0 c and with gas parameter na 3 = 10 −4 , for four values of the couplings ratio g 12 /g.The dotted blue lines in the panels (a) and (b) are the parabolas obtained from the mean-field prediction of the magnetic susceptibility χ = 2/(g − g 12 ), the green dashed lines are the HF predictions, the red solid lines are the Popov predictions.In panel (d) only the HF line is shown.The vertical lines indicate the critical polarization p c .

Figure 3 :
Figure 3: Free energy difference for a mixture with na 3 = 10 −4 , and g 12 /g = 0.5 and T = 0.794T 0 c as a function of the polarization squared.The PIMC results are compared with the T = 0 mean-field (MF -blue dotted line), HF (green dashed line) and Popov (red solid line) predictions.The vertical line indicates the critical polarization p c .

Figure 5 :
Figure5: Chemical potential of an unpolarized mixture with na 3 = 10 −6 , and g 12 /g = 0.93 as a function of the temperature.The PIMC results in the thermodynamic limit (black points) are compared with the HF (green dashed line) and the Popov (red solid line) predictions.

Figure 6 :
Figure 6: Intraspecies (panel (a)) and interspecies (panel (b)) contact parameters of an unpolarized mixture with na 3 = 10 −6 , and g 12 /g = 0.93 as a function of the temperature.The PIMC results in the thermodynamic limit (black points) are compared with the HF (green dashed line) and the Popov (red solid line) predictions.

Figure 8 :
Figure 8: Chemical potential of an ideal Bose gas with N particles at temperature T = 1.5T 0c .The PIMC values (black diamonds) are compared with the exact results at fixed N (connected by the blue dotted line) and with the result in the thermodynamic limit (green horizontal line).Inset: the difference between the PIMC and the exact values.

11 )Figure 9 :
Figure9: Results for an interacting Bose gas with gas parameter na 3 = 10 −6 as a function of the temperature.The PIMC values, extrapolated to the thermodynamic limit (black crosses), are compared with the perturbative results of Hartree-Fock (green dashed line) and Popov (red solid line) theories.Left panel: Results for the chemical potential, also compared with the predictions from the universal relations of Ref.[26] (blue dots).Right panel: difference in freee energy with the ideal Bose gas result.

Figure 10 :
Figure10: Chemical potentials µ 1 (blue) and µ 2 (red) for a non-interacting mixture of ideal Bose gases as a function of the polarization.The temperature is kept fixed at T = 0.5T 0 c and the total number of particles is N = 128.The PIMC points are compared to the exact results connected by the dotted lines.The critical polarization, at which the minority component becomes normal, is signaled by a gray vertical line.Inset: the difference between the PIMC data and the exact results.

12 n 1 n 2 + gn 0 T 2 + 1 βV k ln 1 −T 2 + 2 + 1 βV k ln 1 2 )/λ 3 T= 1 2 gn 0 ± (g 2 − g 2 12 )n 2 p 2 + g 2 12 n 2 0 ,
e −β(ϵ k +gn 1,0 ) + ln 1 − e −β(ϵ k +gn 2,0 ) , (C.19) holding when both condensates are present, i.e. in the polarization range p < p c set by the critical polarization p c = 1 − (T /T 0 c ) 3/2 at which the minority component 2 turns normal.Here ϵ k = ℏ 2 k 2 /(2m) is the single-particle kinetic energy and n 0 T = ζ(3/2)/λ 3 T is the noninteracting thermal density written in terms of the thermal wavelength λ T = 2πℏ 2 /mk B T and ζ(3/2) ≃ 2.612.Furthermore, n i,0 (i = 1, 2) correspond to the condensate density of the two components calculated to lowest order in the interaction strength: n i,0 = n i − n 0 T .When p > p c and the density n 2 of the minority component does not exceed the thermal density n 0 T , the above expression for free energy becomesg 12 n 1 n 2 + µ IBG 2 n − e −β(ϵ k +g(n 1 −n 0 T )) + ln 1 − e −β(ϵ k −µ IBG 2 ) , (C.20)where the effective chemical potential µ IBG 2 is fixed by the normalization condition of the minority component n 2 = g 3/2 (e βµ IBG , with g 3/2 (z) the usual special Bose function.Notice that expressions (C.19) and (C.20) coincide at p = p c where n 2 = n 0 T and µ IBG 2 = 0.The Popov theory includes the contribution from collective excitations (density and spin waves) into the thermodynamics of the mixture yielding the following expression for the free energy: components are in the condensed state (p < p c ).The first term in the second line collects the thermal contribution from the excitation spectrum in the density and spin channel E ± k = ϵ 2 k + 2Λ ± ϵ k whereas the last term survives also at T = 0 yielding the Lee-Huang-Yang beyond mean-field corrections to the ground-state energy.Both terms involve the effective chemical potentials Λ ±