Finite temperature spin dynamics of a two-dimensional Bose-Bose atomic mixture

We examine the role of thermal fluctuations in uniform two-dimensional binary Bose mixtures of dilute ultracold atomic gases. We use a mean-field Hartree-Fock theory to derive analytical predictions for the miscible-immiscible transition. A nontrivial result of this theory is that a fully miscible phase at $T=0$ may become unstable at $T\neq0$, as a consequence of a divergent behaviour in the spin susceptibility. We test this prediction by performing numerical simulations with the Stochastic (Projected) Gross-Pitaevskii equation, which includes beyond mean-field effects. We calculate the equilibrium configurations at different temperatures and interaction strengths and we simulate spin oscillations produced by a weak external perturbation. Despite some qualitative agreement, the comparison between the two theories shows that the mean-field approximation is not able to properly describe the behavior of the two-dimensional mixture near the miscible-immiscible transition, as thermal fluctuations smoothen all sharp features both in the phase diagram and in spin dynamics, except for temperature well below the critical temperature for superfluidity.

The interaction in dilute gases of ultracold atoms is entirely determined by the s-wave scattering length, a, which in turn fixes the interaction coupling constant g entering all relevant equations at the mean-field level. In a two-component mixture, one has two intra-component coupling constants, g 11 and g 22 , and one inter-component coupling constant g 12 . The theoretical constraint for phase-separation at zero temperature is that these constants must satisfy the inequality g 2 12 > g 11 g 22 [1,2]. However, at finite temperature, deviations from this constraint are expected to appear. Existing theoretical studies have mainly addressed quasi-one-dimensional (1D) or 3D systems employing the mean-field treatment including the Hartree-Fock [25][26][27], Hartree-Fock-Bogoliubov-Popov [9,28,29], and Zaremba-Nikuni-Griffin formalism [30,31]. Stochastic growth dynamics for quasi-1D and quasi-2D multi-component mixtures has been investigated in Refs. [32,33]. However, there have been few finite temperature studies for 2D homogeneous Bose mixtures using beyond mean-field theory [34,35]. Recently, Ota et al. [10] predicted a temperature induced magnetic phase transition in an uniform 3D Bose-Bose mixture us- * These two authors contributed equally to this work.
ing the Popov theory. It is then natural to ask whether such a phase-transition also exists in 2D.
An important feature of 2D Bose gases is that true Bose-Einstein condensation is not possible at finite temperature due to Mermin-Wagner-Hohenberg theorem [36,37]. Nevertheless, superfluidity still exists below the critical temperature T BKT for the Berezinskii-Kosterlitz-Thouless (BKT) phase transition [38][39][40][41]; the transition from normal gas to superfluid follows from the binding and unbinding of vortex-antivortex pairs at T BKT [42]. Observation of such transition in the domain of ultracold quantum gases has been possible with quasi-uniform box traps [43,44]. One major advantage of using 2D box traps is that the dynamics of phase-separation depends only on the interplay between kinetic and interaction energy and one can encounter a novel phase-diagram which would otherwise be non-existent due to trap inhomogeneity. Finally, the 2D planar configuration is more suitable to access local densities fluctuations with respect to a 3D system, in which most of the time only column density can be measured.
This sets the stage for our current work. Our goal is to extend the investigation to the beyond mean-field level and explore the case of a uniform 2D Bose-Bose mixture occupying two different hyperfine states and satisfying the miscibility condition at zero temperature. As a first step, we calculate the phase-diagram as a function of temperature and inter-species interaction strength and identify the miscible and immiscible regions using the mean-field Hartree-Fock (HF) theory [25,27]. By treating the atoms in the two components as up and down spin states of spin-1/2 particles, the transition from a miscible to an immiscible mixture can be seen as a magnetic phase-transition. In the HF theory it occurs where the spin susceptibility exhibits a sudden discontinuity. A remarkable result of the theory is that, in the superfluid regime, temperature increase tends to favor phase separation. This result seems counter-intuitive, as thermal fluctuations usually acts against order. An example of a similar type of anomaly in the context of classical fluid is the temperature driven phase-transition in watertriethylamine mixtures [45].
As a second step, we use the Stochastic (Projected) Gross-Pitaevskii (SGPE) [46,47] theory for the same mixtures which includes effects of thermal fluctuations both in the density and spin channels, going beyond the HF description. We calculate the density profiles of the two components at equilibrium at different temperatures and interaction strengths and we simulate the response of the system to a weak external perturbation producing a spin oscillation (i.e., an out-of-phase motion of the two components). From the spin dynamics we also obtain the velocity of spin sound waves. The SGPE simulations show evidences of the miscible-immiscible phase transition in the expected range of parameters, in qualitative agreement with the HF theory, however thermal fluctuations makes the transition much smoother, thus hindering the sharp features predicted by the mean-field theory. Only for temperature T 0.5T BKT , the two theories provide consistent results.

II. MISCIBLE-IMMISCIBLE PHASE TRANSITION
A. Hartree-Fock theory We start from a mean-field description of the uniform 2D interacting Bose mixtures, provided by the Hartree-Fock theory. This theory has been widely used in three dimensions, to investigate the equilibrium [48,49] properties of the mixtures at finite-temperature. At equilibrium, the chemical potential of a given component i = {1, 2} below the superfluid-normal phase transition is given by [2,25,27] where n i = n i0 + n iT is the density of atoms of each component. The coupling constants for the intra-species interaction, g ii , and inter-species interaction, g 12 , are related to the s-wave scattering lengths a ij and mass m 1 = m 2 = m as g ij = 4π 2 a ij /m. In our work, we consider the case of symmetric mixtures with g 11 = g 22 = g and N 1 = N 2 = N/2 as the number of atoms. Although the above expression Eq. (1) is similar in three or two dimensions, the physical meaning behind the quantity n i0 is formally different; while in 3D the quantity n i0 corresponds to the condensate density, this identification no longer holds in 2D, where Bose-Einstein condensation (BEC) is ruled out at finite temperature [36,37] and n i0 corresponds to the quasi -condensate density, which characterizes the suppression of density fluctuations. According to [50,51], one can define whereΨ i is the bosonic field operator and · · · denotes the statistical average. For sufficiently low temperatures, the quasi-condensate is known to play the same role as the true condensate [52]. As for the thermal component, it can be defined through the Bose distribution func- where S is the area and β = (k B T ) −1 . This expression takes into account the mean-field effects at the level of the single-particle energy through the expression z i = exp[β(µ i − 2g ii n i − g 12 n 3−i )] for the fugacity. The thermal atoms density thus takes the explicit form where λ T is the thermal de Broglie wavelength. It is worth noticing that Eq. (3) yields a divergent behavior for µ i = 2g ii n i + g 12 n 3−i , in accordance with the fact that BEC does not exist in 2D at finite temperature. However, in the BKT superfluid phase one can safely use expression (1) for the evaluation of thermodynamic quantities.

B. Thermodynamic quantities
In order to calculate the phase diagram of the mixture, one can investigate the dynamical stability of the system. The mixture is found to be stable against density and spin fluctuations if the compressibility κ T and spin susceptibility κ M are positive [53,54]. The two quantities are defined through the chemical potential Eq. (1) following the thermodynamic relation κ T (κ M ) = [∂(µ 1 ± µ 2 )/∂(n 1 ± n 2 )] −1 . In the particular case of the symmetric mixtures considered in this work, the compressibility as well the susceptibility assume the simple form where we have introduced the isothermal compressibility for the single-component Bose gas: where n 0 = n 10 + n 20 is the total quasi-condensate atom density. At T = 0, Eq. (5) yields the well known result κ sc T = 1/g, and the spin susceptibility (4) reduces to κ M = 2/(g − g 12 ), revealing its large increase near the miscible-immiscible phase transition occurring for g 12 = g [55,56]. At finite temperature instead, the miscible mixture is stable under the condition 1 − g 12 κ sc T > 0. Remarkably, the isothermal compressibility in Eq. (5) increases at finite temperature, compared to the value predicted at T = 0. This behavior is the direct consequence of exchange effects [1,2], which is responsible for both the factor 2 in the denominator of Eq. (5), as well as the temperature dependence of the chemical potential in Eq. (1).
In Fig. 1(a), we show the single-component compressibility κ sc T calculated from Eq. (5), using a typical interaction parameter mg/ 2 = 0.095 [57]. For comparison, we also show the results obtained from the universal relation. The latter predicts that the equation of state of a single-component Bose gas of density n sc in the vicinity of the critical density, n BKT , depends on a single variable X related to the interaction g and to the reduced chemical potential βµ according to n sc − n BKT = f (X). The universal function f has been evaluated in Ref. [52] using classical Monte-Carlo simulations. In addition, universal relations also provide a prediction for the BKT superfluid transition temperature of a single-component  7), respectively. The red dashed line is an estimate of the location of the second magnetic instability, above T 0 BKT , obtained from universal relations, namely, by interpolating the data points for the equation of state given in Ref. [52]. According to the HF theory, the mixture is miscible above the blue and red lines, while it is dynamically unstable against phase separation in the gray area. The array of markers correspond to the parameters used in SGPE numerical simulation and the letters serve as a guide to read Figures 3 and 7. uniform Bose gas: T 0 BKT = 2π 2 n sc /[mk B ln 380 2 /mg ] [42]. It is worth noticing that the Hartree-Fock theory is able to capture the qualitative behavior of the isothermal compressibility below T 0 BKT , with its characteristic increase when T approaches the BKT transition temperature, although it does not predict the typical peak just above the critical point.
Our result for the spin susceptibility of the mixture, Eq. (4), is reported in Fig. 1(b) for g 12 /g = 0.9. Temperature is normalized to T 0 BKT which for the mixture we define according to Since the inter-species interaction is expected to play a little role for the critical point of the BKT superfluid phase transition [34,35], the value of T 0 BKT for a symmetric miscible mixture is close to the actual value of BKT critical temperature; it significantly underestimates it only in the case of full phase separation. Remarkably, the instability condition κ sc T = 1/g 12 , which is graphically represented in panel (a) by the intersection between the compressibility curve and the red horizontal line, takes place at a temperature well below the BKT transition, and the spin susceptibility diverges at The pole of the spin susceptibility 1 − g 12 κ sc T = 0 identifies the temperature T M at which a magnetic instability occurs. We can calculate its location for different values of the interaction strengths. The result is shown as a blue line in Fig. 2. It is worth mentioning that the decrease of κ sc T above T 0 BKT predicted by the universal relations in Fig. 1(a) suggests that the instability condition characterizing the magnetic phase transition should also be satisfied above the BKT transition. We estimate its location by interpolating the data points for the equation of state given in Ref. [52] and we plot the corresponding T M as a red dashed line. In the interaction vs. temperature phase diagram of Fig. 2, the mixture is expected to be dynamically unstable against phase separation in the gray area below the blue or red lines.
We can also derive an analytical result for the temperature T M by using the small-g expression for the single component compressibility: This simple expression allows one to see how the isothermal compressibility reflects the universal behavior of a 2D Bose gas. Indeed, κ sc T /κ sc T (T = 0) does not explicitly depend on the value of the coupling constant g and, in the regime T T 0 BKT , one has D 0 λ 2 T n 1 with n = n 1 + n 2 , thus obtaining the simple estimate valid for gn k B T 0 BKT and 1 − g 12 /g 1. This expression is shown as a black dotted line in Fig. 2, which nicely agrees with the behavior of the pole of the spin susceptibility (blue line) at low T .
The above results for the magnetic instability of the binary mixture suggest the occurrence of a first order phase transition, the value of T M corresponding to the spinodal temperature above which the unpolarized uniform configuration of the mixture is dynamically unstable. The actual transition to a demixed configuration is then expected to take place at smaller values of the temperature and could be identified by comparing the free energy of the uniform unpolarized configuration with the one of the phase-separated configuration. In a similar way, the phase diagram for the 3D Bose mixtures has been obtained in Ref. [10], by means of the Popov theory. The equilibrium configuration in the new phaseseparated phase was found to be characterized by a full space separation of the Bose-Einstein condensed components of the two atomic species, their thermal components remaining instead mixed, but with a finite magnetization. However in 2D, the calculation of the free energy, as well as the characterization of the new phase can not be assessed within the actual HF theory. In fact, any mean-field framework based on the notion of the quasi-condensate is expected to fail above the BKT transition point, where vortex proliferation destroys quasilong range order, and the quasi-condensate becomes ill defined. We therefore need a reliable theoretical framework, which allows for the description of the 2D binary mixtures in both the superfluid and normal regimes.
The c-fields represent the coherent region of the energy spectrum, which includes a large but finite number of low-lying modes up to an energy cutoff i cut . The energy cut-off is chosen as [46,47,[69][70][71][72] i This choice guarantees that the mean occupation of the modes below i cut is of order ∼ 1 or larger. The projec-torP maintains the c-fields within the coherent region at each step of the numerical solution. The modes above the cut-off represent the incoherent region of the energy spectrum; it is the source of a stochastic Gaussian random noise which satisfies the following fluctuation-dissipation theorem where · · · denotes the averaging over different noise realizations. The amount of coupling between the coherent and incoherent regions is fixed by the parameter γ, which accounts for the thermal equilibration rate. In this work, we choose γ = 0.01, which is the same of Ref. [68]. Similar values were also used in [70] and [73]; in the latter case, the parameter γ was optimized to reproduce typical experimental growth rates of single component condensates in 3D. It is worth mentioning that in SGPE individual results obtained with independent noise realizations are equivalent to the individual results obtained from independent experimental runs; due to the random nature of the noise, the outcomes of each noise realizations will differ from one another as is the case in experiments. Furthermore, since SGPE describe a grand-canonical ensemble, the number of atoms is not a conserved quantity and one has to use the chemical potential in order to normalize the number density to the desired value. The density of the atoms in the c-fields is n i (x, t) = |ψ i (x, t)| 2 .
In order to perform simulations which are meaningful for feasible experiments, we choose our parameters compatible to those of the ongoing experiments [43,57], with confinement in rectangular box potential in the x-y plane and harmonic trap in the z-direction. The confinement along z is sufficiently strong to freeze all degrees of freedom in that direction. In this work, we consider a uniform 2D hard-wall square box of dimensions L x × L y = (50 × 50)µm. The frequency of the harmonic confinement in the experimental geometry, ω z , can be used to relate the actual 3D s-wave scattering length a ij of the atoms in different hyperfine levels to the 2D coupling constants g and g 12 used in this work, related by g ij = √ 8πa ij /a z such as to ensure that our simulations correspond to realistic configurations. Here, a z = /mω z is the oscillator length along z. A typical atom number in our simulation is 1.8 × 10 4 in each component, which corresponds to µ 0 /k B = 4.8 nK and T 0 BKT = 33.25 nK, if ω z = 2π × (1500 Hz).

IV. EQUILIBRIUM DENSITY PROFILES
We obtain the equilibrium configurations at a given temperature T by numerically propagating Eq. (8) in real-time starting from purely random c-fields until equilibrium is reached, with the constraint µ 1 = µ 2 = µ; for a finite-size uniform 2D Bose-Bose mixture one can choose the value of the chemical potential from the relationship µ ≈ (1 + g 12 /g)µ 0 obtained from Eq. (1) at zero temperature, µ 0 being the chemical potential of a single component Bose gas. We have verified, using the universal relations of Ref. [52], that for a single-component 2D Bose gas, the temperature dependence of the chemical potential is weak enough, so that one can use its T = 0 value in the whole region of temperature of our interest.
In the simulations one has also to consider that, depending on the interaction parameters, the initialization of chemical potential µ can have dramatic effects on the composite system at equilibrium. In particular, with equal µ i and g 12 /g > 1 and for certain noise realizations, the number of atoms (N 1 or N 2 ) of one of the components may decrease considerably in order to lower the total energy. Furthermore, in order to probe the physics of the binary mixture through the equilibrium density profiles, one has to rely on individual simulations. In fact, ensemble averaging over different stochastic realizations, as mentioned before, would wash away signatures of phase-separation because of the different possible spatial orientation of the density patterns. These consequences stem from the random nature of the model. For the present study, since µ 1 = µ 2 , we ought to use the thermally equilibrated profiles obtained from individual runs with N 1 ≈ N 2 for further analysis.
The density profiles obtained from SGPE simulations with different values of g 12 /g and T /T 0 BKT are shown in Fig. 3. Each panel represents the density n 1 (x, t) obtained in the evolution of a mixture starting from a single noise realization. The density n 2 (x, t) is similar, except that it is larger when n 1 is smaller, in such a way that the total density is almost uniform in the box, apart from thermal fluctuations. The first two rows correspond to values of the interaction parameters for which HF theory predicts full miscibility at all temperatures. Conversely, the lowest row falls into the gray area of Fig. 2 where the mixture is expected to be immiscible. The SGPE density profiles are consistent with these predictions, with phase separation clearly visible in the last row, especially at low T . The third and fourth rows of Fig. 3 are those where signatures of the miscible-immiscible phase-transition are expected to appear. According to HF theory, the parameters of panels (o), (s) and (t) are such that the mixture should be phase-separated, though it is miscible at T = 0 for the same interaction strength. Actually, SGPE shows random patterns, with archipelagos of atoms, albeit marred by fluctuations in the equilibrium density profiles, having a close resemblance to the density profiles for g 12 /g = 1.1, in panels (x) and (y), where the immiscibility tends to get suppressed because of thermal fluctuations.
From the equilibrium density profiles we can also calculate the mixing entropy and the overlap integral. To obtain the mixing entropy we divide the 2D box in N cell cells and compute the number of atoms n 1,j and n 2,j in the cell j, with j = 1, 2, ..., N cell for an individual SGPE realization. For the present work N cell equals the number of grid points used for numerical discretization. Then we use the definition [74][75][76] S mix = − 1 2 N cell j=1 n 1,j n 1,j + n 2,j ln n 1,j n 1,j + n 2,j + n 2,j n 1,j + n 2,j ln n 2,j n 1,j + n 2,j .
The result is given in Fig. 4(a). Where S mix has been averaged over N = 40 realizations. In the low temperature limit the mixing entropy is large for g 12 /g < 1 (miscible mixture) and small for g 12 /g > 1 (immiscible), while at large temperature all the curves tends to the same value. A similar behavior is also found in the temperature dependence of the overlap integral defined as shown in panel Fig. 4(b). Here too, Λ has been calculated for equilibrium density profiles obtained from a single SGPE simulation, and then averaged over N .
In both cases, the purple and green line corresponding to g 12 /g = 0.8, 0.9, respectively tend to bend more toward lower entropy and lower overlap than the curves for smaller values of g 12 /g, the difference being larger at finite T than at T = 0. This mild tendency to demixing seems to qualitatively agree with the expectations of HF theory; however, the disagreement at the quantitative level remains large, as there is no sign of sharp transitions. The convergence of all curves to a constant value at large T is consistent with the presence of the Gaussian white noise source Eq. (10) in the SGPE. Ideally, in SGPE formalism, the ensemble averaging of individual c-field realizations is first performed to suppress the effects of random noise, followed by the computation of different physical quantities [47]. However, in the cur-rent work, we first calculate the pertinent observables corresponding to the equilibrium c-field density profiles obtained from a single SGPE run, followed by its averaging. This is essentially done to retain the signatures of immiscibility within the equilibrium solution, which, as mentioned earlier can get washed away because of ensemble averaging at the beginning. On the contrary, the inherent presence of noise being dominant at high temperature within the individual realization of equilibrium density profiles makes the values of S mix and Λ saturated.

A. Center of mass drift
The miscibility of binary Bose mixtures can also be assessed within the SGPE theory by evaluating the spin dynamics. For this purpose, we follow the approach adopted in the experimental work of Ref. [55], and evaluate the spin center of mass oscillation by applying an external potential. We first generate the equilibrium density profiles using Eq. (8) at given interaction strengths and temperature as described in the previous Section. The mixture is then suddenly subjected to a weak potential tilt, acting on the two components in opposite directions. We choose a linear potential of the form V ext = V 0 xσ z θ(t), where V 0 is the strength of the potential which is of the order 10 −2 µ 0 and σ z the Pauli matrix (see Fig. 5). The system is then let to evolve in absence of dissipation, i.e., with γ = 0. A similar experimental protocol has been employed to measure the speed of sound in a box for a single-component Bose gas, both in 2D [57] and 3D [77].
The application of the external potential renders each species of the mixture adapt itself to a displacement along the direction of the potential minima. These bring about marked changes in the density and sound waves are emitted. From the evolution, we extract the center of mass drift along the x and y direction through the definition The dynamics of the center of mass drift through Eq. (13) is then averaged over N = 20 ∼ 40. In what follows, we investigate in particular the spin dipole moment, given by the center of mass drift according to We notice that the method of evolving in time the states stochastically generated at t = 0 − , using the projected Gross-Pitaevskii equation for the classical field, has also been earlier used to model the growth of quasicondensate on an atom chip [59], and is similar in essence to the truncated Wigner method for Bose condensed gases [78].

B. Linear Response Theory
Before presenting the numerical results for the dynamics of the spin dipole moment we address the problem by means of linear response theory [2]. In particular we evaluate the dynamic spin response function of the 2D Bose mixtures using the Random Phase Approximation (RPA): [79] χ M (k, ω) = 2χ(k, ω) 1 − g 12 χ(k, ω) , with χ(k, ω) as the density response function for a singlecomponent Bose gas, The reference response functions for the quasi-condensate χ 0 0 and thermal atoms χ 0 T are evaluated consistently within the HF description of Sec. II (we consider the external perturbation along the x-direction k = ke x and the long wavelength limit k → 0): where f (p) = [e (p 2 /(2m)+gn0/2)/(k B T ) − 1] −1 is the Bose distribution function of the 2D ideal Bose gas. We note on passing that the density response function used in Ref. [68] is retrieved from Eq. (16) by explicitly putting n 0 = 0 and replacing 2g by g, thus considering the gas as a normal gas with suppressed density fluctuations. Although such procedure is found to give a good description of the 2D Bose gas in the normal phase, in the present work one needs to properly take into account the effects of density fluctuations through the temperature dependence of the quasi-condensate, in order to describe the magnetic phase transition. The present RPA is compatible with the quasi-condensate HF theory of Sec. II, and one can immediately verify that the compressibility sum-rule χ M (k → 0, ω = 0) = κ M yields the HF result for the spin susceptibility Eq. (4).
In order to evaluate the linear response of the system under the perturbation described in Fig. 5, let us simplify the problem and consider a uniform mixture subjected to a sudden application of a spatially periodic potential, producing a sinusoidal modulation of wavevector q 0 = π/L x . This essentially amounts to neglect finite-size effects coming from the hard-walls. Under such perturbation, the magnetization density starts to oscillate with a time-dependent amplitude given by [80,81]: corresponding to the x-component of the spin dipole moment, and where χ M is the imaginary part of the spin response function given by Eq. (15).

C. Spin Oscillation Dynamics
In Fig. 6, we first show the spin oscillation dynamics M(t) obtained from SGPE simulation at low temperature T /T 0 BKT = 0.1, for two different values of coupling constant g 12 /g = 0.5 (upper panel) and = 1.1 (lower panel). As one could expect from the zero-temperature miscibility criterion g 12 < g, we find that for g 12 /g = 0.5, the two components are fully mixed, going back and forth in the direction where the external potential is induced (blue solid line) until reaching eventually the new equilibrium position given by the spin susceptibility (see Eq. (18)). On the other hand for g 12 /g = 1.1, the two components are already phase-separated at t = 0 and the effects of external potential is to induce a density oscillation in the respective component. Another signature of phaseseparation is provided by the oscillation in the direction perpendicular to the external potential, given by the red dotted line in Fig. 6. While for the miscible case this oscillation is relatively small, in the g 12 /g = 1.1 case phase-separation is responsible for the coupling between the x and y components, leading to a well-defined oscillation in the perpendicular direction too. For the miscible case in the upper panel of Fig. 6, we also show the results obtained from RPA approach Eq. (18). In RPA, the oscillation frequency at T 0 is given by the spin mode of the Bogoliubov sound as ω = c s q 0 , with c s (T = 0) = (g − g 12 )n/(2m) where n = n 1 + n 2 is the total density of atoms, and we find a good agreement with SGPE simulation. As for the phase-separated case we have verified that both components oscillate at the frequency of the Bogoliubov phonon mode for a singlecomponent Bose gas occupying half area of the box.
A clear distinction between the behavior of a miscible mixture and a phase-separated gas, as in Fig. 6, can only be expected at low temperature. In particular, a well defined phase-separated dynamics can be seen only in configurations like in panels (u) and (v) of Fig. 3. The dynamical response of the system at higher T and closer to the region of magnetic instability is more complex.
In Fig. 7, we show the spin center of mass oscillation obtained with SGPE simulations starting from each of the configurations in panels (a)-(t) of Fig. 3, followed by averaging over N . In order to extract the frequency of the oscillation we use the expression where A 0 ∝ κ M is related to the spin susceptibility, whilẽ ω = c/q 0 and Γ correspond to the frequency and damping rate of the spin sound mode, respectively. Equation (19) has been obtained by assuming a damped harmonic oscillator model for the response function χ M in Eq. (18).
The fitted functions are shown as green dotted lines in Fig. 7. While we obtain good fit at low and intermediate temperature, the fitting seems to become less accurate at high temperature, where the oscillations are strongly damped and thermal fluctuations become dominant. The values of the spin sound velocity c extracted from the fit to the SGPE oscillations are shown in Figure 8, where they are normalized by the single-component Bogoliubov sound c 0 = gn/(2m). The solid lines correspond to the RPA predictions for the same interaction strengths. In the superfluid phase T < T 0 BKT , the spin sound velocity is found to decrease by increasing both the temperature and the ratio g 12 /g. Remarkably, RPA predicts a vanishing sound velocity at finite temperature, associated to the magnetic instability discussed in Sec. II. This can be understood by recalling that the ratio between the energy weighted and inverse energy weighted moments of the dynamic response function provides an estimate for the mean excitation energy [2]: Therefore, the divergence of the spin susceptibility is responsible for the vanishing of spin sound mode. On the other hand, the speed of sound extracted from SGPE simulation do not show this kind of behavior. This discrepancy points out that, contrary to the single-component 2D Bose gas where mean-field approach provides a qual-itatively correct picture of the system below T 0 BKT (see Fig. 1(a) and Ref. [68]), the HF theory seems inadequate in describing 2D Bose mixtures, except at low T where the sound velocity is indeed close to the SGPE results.
As regards the damping rate, in general, we observe that the oscillation becomes broad and damped as one increases the temperature (that is, from left to right in each row). This is likely due to Landau damping, arising from the coupling between the collective sound mode and the thermally populated single-particle excited states, which is properly included in the SGPE theory (see Ref. [68] for a discussion in the case of a single-component 2D gas). However, the quality of the fit, especially at high T , is not good enough to permit a quantitative analysis. In addition, in the T → 0 limit, the SGPE theory, while correctly approaching the T = 0 results for the speed of sound, is not expected to be reliable in determining the damping rate, which is more sensitive to the way in which fluctuations are included in the theory; this may also be the origin of the anomalous high damping in panels (k) and (p) in Fig. 7.  FIG. 8. Spin sound velocity normalized to the Bogoliubov speed of sound for a single component Bose-gas calculated for mg/ 2 = 0.095. From top to bottom: g12/g = 0.2, 0.5, 0.8, 0.9. Points refer to the spin sound velocity obtained from SGPE simulation extracted using the fitting function (19), while lines are the RPA prediction.

VI. CONCLUSIONS
We have carried out a systematic investigation of the role of thermal fluctuations in two-dimensional homogeneous bosonic mixtures of dilute atomic gases. To this end, we have employed the Hartree-Fock mean field framework complemented with beyond mean-field effects through Stochastic (Projected) Gross-Pitaevskii formalism to demonstrate the role of finite temperature on the miscibility-immiscibility phase transition. A remarkable result predicted by the mean-field theory is the existence of phase-separation induced by thermal fluctuations, occurring even for mixtures which are miscible at zero temperature (g 12 < g). Analytical predictions based on mean-field HF theory are then compared with the equilibrium density solutions obtained by numerically solv-ing SGPE at different coupling strengths and temperature. The spin dynamics of the mixture brought about by a weak external perturbation is also simulated, which shows some qualitative agreement with the mean-field predictions for temperatures well below the critical temperature for superfluidity. However, with the inclusion of beyond mean-field effects through SGPE, both the density profiles as well as spin dynamics are devoid of any sharp distinctive features near the phase transition point predicted by the HF theory. This, we believe, would also be the scenario in actual experiments. The recent availability of box like traps would be of value to corroborate these features.
Important open issues concern the extension of the HF theory to take into account effects of thermal fluctuations in the spin channel, to clarify the observed discrepancy between the mean-field and stochastic approaches. Indeed, while HF theory is widely used in describing singlecomponent Bose gas as well as 3D mixtures in the density channel, our work points out its failure for the investigation of spin physics. Development of beyond mean-field theory for the 2D Bose mixtures along the line of those in Refs. [82,83] could provide further insight on magnetic properties of binary mixtures at finite temperatures.

ACKNOWLEDGMENTS
We are grateful to Stefano Giorgini and Sandro Stringari for many fruitful discussions concerning the results reported in Sec. II. A. Roy thanks S. Gautam for useful discussions. This work is supported by Provincia Autonoma di Trento and from INFN-TIFPA under the project FIS . A.R. acknowledges support from the Italian MIUR under the PRIN2017 project CEnTraL. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. This material is based upon work supported with GCP research credits by Google Cloud.