Thermal suppression of demixing dynamics in a binary condensate

We investigate the demixing dynamics in a binary two-dimensional (2D) Bose superfluid using classical-field dynamics. By quenching the interspecies interaction parameter, we identify a strong and weak separation regime depending on the system temperature and the quench parameter. In the strong separation regime our results are in agreement with the inertial hydrodynamic domain growth law of binary fluids and a Porod scaling law for the structure factor at zero temperature is found. In the weak separation regime thermal fluctuations modify both the domain growth law and the Porod tail of the structure factor. Near the superfluid transition temperature the scaling dynamics approaches the diffusive growth law of a 2D conserved field. We then analyze the demixing dynamics in a box cloud. For low quench we find distinctive domain dynamics dictated by the boundary condition. Otherwise, the dynamics are qualitatively similar to those of systems with periodic boundary conditions.


I. INTRODUCTION
When two immiscible fluids such as water and oil are allowed to mix, they separate into two distinct phases [1].Such phase separation phenomenon is very well established in science, with relevant implications for important technological applications [2][3][4].In physics, phase separation occurs in a variety of condensed matter systems such as polymers, fluid mixtures, gels, ferroelectrics, membranes, superfluids, superconductors and the like.
According to classical theories of phase-ordering dynamics, the domain growth follows a characteristic power-law behavior L(t) ∼ t η , where L is the average domain size and η is the scaling exponent.The dynamics is universal such that the time evolution of an observable is solely governed by L(t).In practice, this scaling hypothesis is tested by the equal-time correlation function C(r , t) = φ(r + r , t)φ(r, t) , where φ(r) is an order parameter characterizing the dynamical evolution of the system and .. denotes the statistical average.The Fourier transform of C(r, t) is the structure factor S(k, t) = φ k (t)φ −k (t) , where φ k is the Fourier transform of φ(r).From a dimensional consideration the structure factor obeys the scaling relation S(k, t) = L d f (kL(t)) [5], where d is the spatial dimensionality and the scaling function f (q) is independent of time.Due to the presence of domain walls, f (q) exhibits a power-law tail f (q) ∼ q −(d+1) at large q, which is referred to as the Porod law [6].The scaling theory hypothesizes various power laws for domain coarsening, which describe the time dependence of characteristic length scales.The domain growth law for a two-dimensional (2D) conserved field is L ∼ t 1/3 , which characterizes the diffusive transport of the order parameter [5,7,8].In binary fluids a competition between the viscous and inertial flow leads to two growth regimes: viscous hydrodynamic (L ∼ t) [9] and inertial hydrodynamic (L ∼ t 2/3 ) [10].The viscous hydrodynamic regime has been confirmed by experiments as well as simulations [5].However, the inertial hydrodynamic regime has not been observed yet as viscous flow is non-negligible in classical fluids.
In this paper, we investigate how the thermal fluctuations influence the dynamical scaling of coarsening in a binary 2D Bose superfluid.To this end, we employ semiclassical-field simulations to address the demixing dynamics at nonzero temperature, which is triggered by quenching the interspecies interaction parameter.As a key result, we show how the characteristic scaling of domain growth is modified by temperature and quench parameter.The interplay of the quench and thermal energy results in two phase separation regimes.The first one is the strong separation regime that occurs at low temperature and high quench, for which the Porod scaling law of the structure factor S(k) ∼ k −3 holds, and the domain  coarsening follows the inertial hydrodynamic growth law of binary fluids, i.e., L(t) ∼ t 2/3 .The other regime is the weak separation regime at high temperatures, where thermal fluctuations modify both the domain growth law and the Porod tail of the structure factor.Near the superfluid critical temperature the domain coarsening approaches the diffusive growth law of a 2D conserved field (L ∼ t 1/3 ) and the Porod tail of the structure factor scales close to S(k) ∼ k −1 .Furthermore, we examine the demixing dynamics in a box cloud and find an intriguing interplay of the box symmetry and the dynamics.In particular, for low quench and small clouds, with sizes comparable to the spin healing length, demixing occurs via the creation of domains of regular patterns due to the boundary condition.For high quench and large clouds we recover the dynamics that is similar to the system with periodic boundary conditions.

II. SYSTEM AND METHODOLOGY
We consider a cloud of 87 Rb atoms in two different hyperfine states |F = 1, m F = 0 and |F = 2, m F = 0 , which is motivated by the experiments [51,52].Thus, the two species have the same masses (m 1 = m 2 = m) and the intraspecies scattering lengths are a 11 /a B = 100.86and a 22 /a B = 94.58 [53], where a B is the Bohr radius.The interspecies scattering length is a 12 /a B = 98.9 [54], resulting in the parameter α = 1.012, which is defined as Since α is slightly above 1, the two species are weakly immiscible and thermal fluctuations play a prominent role in the demixing dynamics, as we show below.We describe the system by the Hamiltonian with and where a = 1, 2 represent the two species and ψa ( ψ † a ) are the corresponding annihilation (creation) operators.The intraspecies interactions g aa and interspecies interaction g 12 are given by, respectively, z = /(mω z ) is the harmonic oscillator length of the trapping potential in the transverse direction, where ω z is the trap frequency.For a condensate with a large number of atoms we replace ψa by complex numbers ψ a .Using Eq. 2 we obtain the coupled equations of motion which govern the dynamics of binary condensates.n a = |ψ a | 2 are the densities.This system hosts two excitation branches of collective modes [55] where E a = k ( k + 2g aa n a ) are the single-component Bogoliubov spectra and k = 2 k 2 /(2m).The coupling g 12 results in hybridized branches E k,± .A direct consequence of this hybridization is that the lowmomentum part of E k,− vanishes when α = 1 (or equivalently g 12 = √ g 11 g 22 ) and becomes imaginary for α > 1.This leads to the creation of unstable modes when α is above 1, which is responsible for the demixing of the two species.We note that α c = 1 is the quantum critical point separating the miscible (α < α c ) and immiscible (α > α c ) states at zero temperature [13].The range of unstable modes is determined by setting E k,− = 0, giving where ξ a = / √ 2mg aa n a are the single-component healing lengths.k 0 vanishes when α = 1 and increases with increasing α for α > 1.The wavelength λ 0 = 2π/k 0 and the lifetime τ = /E k0,− give an estimate of length and time scale for the emergence of domains.We investigate the phase-separation dynamics using the classical-field method of Refs.[56,57].For the numerical simulations we discretize the space on a lattice of size N x × N y and a discretization length l = 1 µm.We note that l is chosen to be smaller than or comparable to the healing length and the thermal de Broglie wavelength [58].This maps the continuum Hamiltonian on the discrete Bose-Hubbard model, which introduces J = 2 /(2ml 2 ) as the tunneling energy and U aa = g aa l −2 and U 12 = g 12 l −2 as the onsite repulsive interactions.We use ω z = 2π × 4.6 kHz, leading to U 11 /J = 0.336 and U 22 /U 11 = 0.938 [51].The quench protocol is described in Fig. 1.We start with a 2D superfluid cloud of total density n = 10 µm −2 at temperature T .The initial states ψ 1 (r) of this system are sampled in a grandcanonical ensemble of chemical potential µ and temperature T via a classical Metropolis algorithm [56].We choose T in a wide range of T /T 0 = 0.1 − 1.1, where T 0 is an estimate of the critical temperature for the superfluid transition in weakly interacting 2D Bose gases [59,60].For the other species we sample the initial states with vacuum fluctuations, i.e., |ψ 2 (r i )| 2 = 1/(2l 2 ) [61], where the index i corresponds to the lattice site and .. denotes the ensemble average.At time t = 0 we use a π/2 pulse to obtain a uniform superposition of the states This results in the two cloud densities n 1 ≈ n 2 ≈ n/2 = 5 µm −2 , since g 11 and g 22 are similar.We then quench g 12 in the demixed regime (α > 1) and determine the time evolution ψ 1/2 (r, t) via Eqs.( 6) and (7).As schematically shown in Figs.1(c) and (d), the initial time evolution proceeds via nucleation of small-sized domains and the long-time evolution results in two spatially separated clouds.For our simulations we vary α in the range 1.03 ≤ α ≤ 1.55 to explore both weakly and strongly immiscible regimes, which covers a wide range of immis-  cible regime that can be realized with mixtures of other species.To analyze the demixing dynamics, as an order parameter, we calculate the local density imbalance We show m(r, t) for a periodic-boundary system in Fig. 2 and its average m(r, t) for a box system in Fig. 6, where .. denotes an average over the initial ensemble.

A. Demixing dynamics
In Fig. 2 we show the time evolution of m(r, t) of a single trajectory at T /T 0 = 0.21 for α = 1.03 and 1.55.We employ a fixed system size of 256 × 256 µm 2 for all periodic-boundary simulations.The quench to the demixed state at t = 0 triggers the nucleation of domains of each component, where smaller domains are present for α = 1.55 than that for α = 1.03.We estimate the initial domain size L 0 by the momentum range k 0 of unstable modes in Eq. 9. We obtain L 0 ∼ λ 0 = 19 and 4.7 µm and the nucleation time t ∼ τ = 26 ms and 1.5 ms for α = 1.03 and 1.55, respectively.The intermediate time evolution manifests the coarsening process, where small domains shrink and large ones grow.There are small patches of other component in the domains, which are due to the initial fluctuations that suppress the dynamics in the weak separation regime.This is the scenario for α = 1.03, whereas for α = 1.55 the dynamics is weakly affected by these fluctuations as the system is in the strong separation regime.The weak versus strong separation regime occurs as an interplay between the thermal and quench energy.At high temperatures thermal fluctuations dominate the dynamics and no phase separation occurs as we show in Appendix A.
To identify the interplay of collective modes in the demixing dynamics we calculate the dynamic structure factor of the density where n 1 (k, ω) is the Fourier transform of the density n 1 (r, t) of component 1 in space and time: T s = 0.55 s is the sampling time for the numerical Fourier transform and N l = N x N y is the number of lattice sites.In Fig. 3 we show S 1 (k, ω) as a function of the wave vector k = kê x and frequency ω for α = 1.03, 1.13, 1.34 and 1.55.We observe both dynamically stable and unstable modes, where stable modes appear as two excitation branches and unstable ones as a broad spectrum of low-energy excitation.We compare these results with the Bogoliubov spectra E k,± of Eq. 8.For our discretized system, the free-particle dispersion takes the form k = 2J[1 − cos(k x l)], where J = 2 /(2ml 2 ) is the tunneling energy.We show the real-valued predictions of E k,± as the continuous lines in Fig. 3, which capture the excitation branches of stable modes for low and intermediate α and show deviations for high α.We also show the imaginary solution of E k,− , which qualitatively captures the broad spectrum of unstable modes.The momentum range of unstable modes increases with increasing α and is close to the predictions of Eq. 9.There is a peaklike excitation at small k corresponding to the structure of a macroscopic domain that the system forms at time t = T s .This excitation peak shifts to a lower k for large α, implying a rapid growth of domains at large α, which is consistent with the dynamics shown in Fig. 2. Furthermore, we compare the two excitation branches of stable modes with the spectra of phase separated clouds.In this case, the cloud density is twice the initial density, i.e., n i,f = 2n i and the Bogoliubov spectrum reads ).This result agrees with the upper branch of S 1 (k, ω) for all α in Fig. 3.The other component being spatially separated from component 1 acts as a thermal cloud whose free-particle dispersion captures the lower branch for all α in Fig. 3.

B. Dynamical scaling
To characterize the scaling behavior we calculate the structure factor of the imbalance with where m(k) is the Fourier transform of m(r).In Figs.
4(a1-e1) we show S(k, t) as a function of the wave vector k = kê x and time t for α = 1.55 and various T /T 0 .The nucleation of domains is indicated by the spectral peak at finite k, which gradually moves to smaller k as the domains coarsen.The location of the peak describes an average size of the domain, whereas the peak broadening reflects the influence of thermal fluctuations on the dynamics.The thermal effect is strong at high temperature, resulting in a decreasing peak amplitude in Figs.4(a1-e1).
For T /T 0 = 0.86, only after a short time evolution, the spectral peak vanishes due to strong diffusion induced by thermal fluctuations.We fit the structure factor with the Gaussian distribution g(k where A 0 , k d , and σ are the fitting parameters.From dependence of the initial superfluid fraction ns/n, which we determine using the method described in Ref. [62].The results are obtained for the system size 256 × 256 µm 2 . k d we determine the average domain size L = 2π/k d .In Figs.4(a2-e2) we show the determined values of L(t) on a log-log scale.The growth of L(t) demonstrates a powerlaw behavior that is typical for coarsening of macroscopic domains L ξ s , where the spin healing length is defined as ξ s = / √ 2mng s , with g s = (2g 12 −g 11 −g 22 )/2.ξ s is a length scale on which the two species interact to nucleate domains.We fit L(t) with the function f (t) = c 0 t η , where c 0 is the fitting parameter.This way, we determine the scaling exponent η, see caption of Fig. 4. We note that the value of η decreases with increasing temperature.In Figs.4(a3-e3) we show the scaled structure factor S(k, t)/L(t) 2 as a function of the scaled wave vector kL(t) for various t and T /T 0 .The different-time results collapse on one single curve, confirming the scaling hypothesis.The momentum tail follows a power-law behavior and its decay varies with temperature.For T /T 0 = 0.21 we find a dependence S(k) ∼ k −3 , which is consistent with the Porod law [5].This occurs due to the presence of domain walls that lead to the dependence S(k) ∼ k −3 at large k in 2D.At high temperature, the momentum tail decays slowly as the process of phase separation is suppressed by dominant thermal fluctuations.For T /T 0 = 0.86 we find a momentum-tail behavior S(k) ∼ k −1. 25 , where no phase separation is visible in the time evolution, see Appendix A. and α influence the value of η, exemplifying the interplay between strong and weak separation regimes.At low T and large α, the system is in the strong separation regime in the sense that the dynamics is marginally affected by initial fluctuations.Here we obtain η close to the zerotemperature prediction η 0 = 0.68, which suggests that the dominant process for domain growth is the inertial hydrodynamic transport of superfluid from low-density to high-density regions.This scenario is consistent with the simulations of binary condensates at zero temperature [50].At high T and small α, initial fluctuations suppress the dynamics of phase separation and result in a renormalization of η and the high-momentum tail of the structure factor.We refer to this regime as the weak separation regime.To quantify these regimes we calculate the time evolution of the average squared imbalance m 2 .In Fig. 5(c) we show m 2 at t = 5.5 s as a function of T /T 0 for α = 1.03, 1.13 and 1.55.It decreases with increasing T /T 0 and decreasing α.As shown in the inset of Fig. 5(c), m 2 increases during the time evolution and reaches the steady state in the long-time evolution.We find that the dynamics is in the strong separation regime for m 2 0.64, where we recover both the zerotemperature η 0 and the Porod tail of the structure factor.The weak separation regime sets in when m 2 0.64, where thermal fluctuations renormalize the scaling parameters.Here the dynamics is influenced by thermal fluctuations that suppress the initial superfluid order of the system [Fig.5(d)].For m 2 0.4 we observe no phase separation.

C. Demixing dynamics in a square box
We now turn to the demixing dynamics of a homogeneous 2D cloud confined in a square-box geometry, which is motivated by the experiments [51,52].Compared to a periodic boundary system, where domain locations are spontaneous, finite boundaries break the translational invariance and act as a pinning potential for the formation of domains [26].We choose the same density and the same quench protocol as above.We first analyze the demixing dynamics in a box cloud of size 64 × 64 µm 2 , which is comparable to the experiments [51,52].In Fig. 6 we show the time evolution of the average imbalance m(x, y) at T /T 0 = 0.21 for α = 1.03, 1.13, and 1.55.Indeed, the nucleation of domains is pinned by the box boundaries, which stems from a density difference at the edges since a 11 and a 22 are different, serving as a seed for the creation of domains.On the contrary, in the case of periodic-boundary systems domain nucleation is seeded from the fluctuations of the field.The box symmetry results in a qualitatively different average dynamics than in infinite systems.For α = 1.03, the time evolution proceeds via formation of regular patterns that undergo a continuous transformation to create structures of striking, geometric shape.At t ∼ 0.5 s the time evolution shows the creation of one macroscopic domain of component 1 and 2, which then equilibrates after t ∼ 5 s.These results are close to the measurements that show the creation of similar structures for time evolution up to 100 ms [51,52].The pinning effect is suppressed when α is high, see the dynamics for α = 1.13 and 1.55 in Figs.6(b) and  (c).The reason for this is the smaller spin healing length ξ s at higher α, which supports the creation of small-sized domains.ξ s is 4.6, 2.2 and 1.1 µm for α = 1.03, 1.13, and 1.55, respectively.For the cloud size considered, we obtain L x /ξ s ≈ 14 and 58 for α = 1.03 and 1.55, respectively, where the former supports the creation of regularshaped domains due to the boundary condition.So, for L x /ξ s 1, the dynamics approaches the one obtained for a system with periodic boundary conditions.
Next, we analyze the growth laws for the domains in a box cloud of sizes between 64×64 µm 2 and 512×512 µm 2 .We calculate the structure factor S(k, t) to determine the average domain size using the procedure described above.From the power-law growth of domains we ascertain the scaling exponent η, which is analogous to Fig. 4. In Fig. 7(a) we show η as a function of 1/L x at T /T 0 = 0.21 for α = 1.03 and 1.55.L x is the linear dimension of the box.The variation of system size by a factor of 64 allows us to perform a reliable finite-size scaling, which gives access to the scaling exponent η ∞ in the thermodynamic limit.We obtain η ∞ = 0.59 and 0.68 for α = 1.03 and 1.55, respectively.The results of η ∞ are close to the values obtained for periodic-boundary conditions.In Fig. 7(b) we show η as a function of 1/L x at T /T 0 = 0.65 for α = 1.55.For this system we find η ∞ = 0.53, confirming the renormalization of the scaling exponent at high temperature.

IV. CONCLUSION
We have studied the demixing dynamics of a binary 2D Bose superfluid using classical-field simulations.By quenching the interspecies interaction parameter we have analyzed the coarsening dynamics at various values of temperature and the quench parameter.We have demonstrated that the dynamical scaling of domain growth in-terpolates between the inertial hydrodynamic growth law of binary fluids and the diffusive growth law of a 2D conserved field.Specifically, for low temperature and high quench we have found the inertial hydrodynamic growth law L(t) ∼ t 2/3 and the Porod scaling law of the structure factor S(k) ∼ k −3 , where L is the average domain size and k is the wave vector.We have pointed out that at high temperature thermal fluctuations suppress the demixing dynamics and modify both the domain growth law and the Porod tail of the structure factor.We have shown that near the superfluid transition temperature the scaling dynamics approaches the diffusive growth law of a 2D conserved field, L(t) ∼ t 1/3 , and the Porod tail scales similar to S(k) ∼ k −1 .We have then studied the demixing dynamics in a box cloud.We have shown that for low quench and small clouds of sizes comparable to the spin healing length the box symmetry gives rise to distinctive dynamics, which is characterized by domains of geometric shapes.By varying the system size we have determined the scaling exponents of the growth law and found them to be consistent with the results of systems with periodic boundary conditions.
Our results highlight the fundamental interplay of the quench and thermal energy in phase separation, which modifies the underlying scaling laws of coarsening dynamics.The experimental realization of our results provides a quantum simulation of scaling laws of binary fluids.Furthermore, Bose mixtures in a ring trap offer the capability to study the solid-body rotation and persistent currents in multicomponent quantum mixtures.

FIG. 2 . 55 FIG. 3 .
FIG. 2. Nucleation of domains and the coarsening dynamics.Time evolution of the two-species density imbalance m(x, y) of a single trajectory for α = 1.03 (upper row) and 1.55 (lower row), displaying nucleation of domains of two components (red and blue) and their coarsening dynamics.The spatial dimensions for each panel are 256 × 256 µm 2 .

Figs. 5 FIG. 6 . 55 FIG. 7 .
FIG.6.Domain formation in a square box of size 64 × 64 µm 2 .(a-c) show the time evolution of the average imbalance m(x, y) at T /T0 = 0.21 for α = 1.03 (upper row), 1.13 (middle row) and 1.55 (lower row).We present the full time evolution as videos in the Supplementary material, which displays a continuous transformation between the domains of different shapes.

FIG. 8 .
FIG. 8. Influence of temperature on the demixing dynamics in a periodic-boundary system of size 256 × 256 µm 2 for α = 1.55.(a-d) show the time evolution of the imbalance m(x, y) of a single trajectory after the quench deep in the demixed state at temperatures T /T0 = 0.43, 0.65, 0.75 and 0.86.The red and blue colors denote the two components.