Phase Separation and Dynamics of two-component Bose-Einstein condensates

The miscibility of two interacting quantum systems is an important testing ground for the understanding of complex quantum systems. Two-component Bose-Einstein condensates enable the investigation of this scenario in a particularly well controlled setting. In a homogeneous system, the transition between mixed and separated phases is fully characterised by a `miscibility parameter', based on the ratio of intra- to inter-species interaction strengths. Here we show, however, that this parameter is no longer the optimal one for trapped gases, for which the location of the phase boundary depends critically on atom numbers. We demonstrate how monitoring of damping rates and frequencies of dipole oscillations enables the experimental mapping of the phase diagram by numerical implementation of a fully self-consistent finite-temperature kinetic theory for binary condensates. The change in damping rate is explained in terms of surface oscillation in the immiscible regime, and counterflow instability in the miscible regime, with collisions becoming only important in the long time evolution.

An important property of a binary mixture is its miscibility. For a homogeneous system, the miscible-immiscible transition is uniquely characterized by the miscibility parameter = (g 11 g 22 /g 2 12 ) − 1. Depending on the strength of the intraspecies (g 11 , g 22 ) and interspecies (g 12 ) interaction, the two components can either overlap in space ( > 0) or phase separate ( < 0). This spatial overlap has practical consequences on, e.g., rethermalization rate [24], coarse-graining dynamics [25][26][27], structures of vortex lattice [28], or instabilities in fluid dynamics [29]. Based on the assumption of overlapping trap centers of the two components, numerical studies [30][31][32][33][34][35] have  39 K mixture at temperature T = 0 for various total numbers of 87 Rb (N Rb , species 1) and 39 K (N K , species 2) atoms. The symmetrical phase is characterized by a normalized trap-center density n norm [see Eq. (4)] while the asymmetrical phase is shown in white. Typical density profiles along the z axis are shown in (e) for asymmetrical phase and (f)-(i) for symmetrical phase. White circles in (a) indicate points accessible experimentally using the Feshbach resonances reported in Ref. [14]. also seen experimentally [11,36], a full characterization of the boundaries separating the three phases is still missing.
The aim of this article is fourfold: (i) to demonstrate that for a trapped binary mixture, = 0 is generally no longer the optimal criterion for the transition boundary; (ii) to characterize the full phase diagram (see Fig. 1) based on the identification of a new parameter; (iii) to propose measurements of the frequency and damping rates of induced dipole oscillations as a universal experimental tool for mapping out the phase diagram; and (iv) to demonstrate the importance of thermal effects on the dynamics, by providing a numerical implementation of a fully self-consistent finite-temperature model for binary mixtures [37]. This extends the successful Zaremba-Nikuni-Griffin (ZNG) model [38] to two components. The validity of the ZNG model has previously been verified in studies of collective modes of single-component BECs [39,40] and macroscopic excitations [41,42].
Typical studies to date have probed the phases by varying a single parameter (e.g., g 12 [32,36] or the number of atoms [11,31]), equivalent to a line scan in a multidimensional parameter space. In these investigations it was thus impossible to obtain a full picture of the problem. Here we present a ground-state phase diagram of the density profiles in a four-dimensional space, treating the ratios of the interaction strengths (g 11 /g 12 , g 12 /g 22 ) and the numbers of atoms (N Rb , N K ) as independent variables [ Fig. 1(a)-1(d)]. Not only do we find a clear deviation of the boundaries from the uniform case (g 11 /g 12 = g 12 /g 22 ), which depends on the atom numbers, but we can also predict whether certain phases (e.g., the asymmetrically separated phase) are accessible to a particular mixture or experiment.
The mapping of the transition boundary in experiments is a nontrivial task. A typical procedure [8,14,16] involves free expansion of the BECs and a measurement of the separation of the COMs of the expanded clouds as a function of . The boundary is then identified as the point where the separation starts to grow. While this method yields quantitative agreement with the transition boundary in the uniform case, we have previously [14] shown that the measured separation can be influenced by the repulsion developed during the expansion dynamics rather than the in-trap phase separation. In the second part of this work, we therefore propose to map the transition boundary of a trapped mixture by monitoring the damping rate and the frequency of dipole oscillations, which we numerically find to be sensitive to the in-trap density profiles. Based on experimentally relevant parameters, our numerical simulations of the oscillation dynamics at both zero and finite temperature indicate an abrupt increase in both the damping rate and the frequency when crossing from the immiscible to the miscible regime, which generally occurs at = 0.

II. MODEL
To address the role of temperature in such dynamics, we use the two-component generalization of the ZNG model, which has been previously demonstrated to pass the stringent test of the undamped Kohn mode, essential for a correct modeling of collective modes [43]. Our kinetic model [37] describes the self-consistent coupling of two BECs, each coupled to their own thermal cloud, and additionally includes coupling between the thermal clouds. This approach enables us to consider the relative importance of damping arising from mean-field coupling (U j c , U j n ) and thermal-condensate (C .. 12 , C kj 12 ) or thermal-thermal (C .. 22 ) collisions. Each condensate wave function, φ j (r) obeys a dissipative Gross-Pitaevskii equation while the Wigner distribution function of the thermal atoms f j (p,r,t) obeys a quantum Boltzmann equation The effective potentials are related to the condensate density n c,j (r) = |φ j (r)| 2 , and thermal atom densityñ j (r) = dp/(2π ) 3 f j (p,r,t), as Here g kj = 2π 2 a kj /m kj denotes the effective interaction strength, where a kj defines the s-wave scattering length between atoms in components j and k, and m −1 defines the reduced mass. The source terms in Eq. (1) are related to the first three collision integrals in Eq. 12 . These collision integrals are sampled in dynamical simulation through the direct simulation Monte Carlo [40] method.

III. SIMULATION PARAMETERS
For our simulations, we consider the 87 Rb -39 K mixture [14], but our analysis is general and can be extended to other mixtures. The atoms are confined in a harmonic potential V j (r) = 1 2 m j [ω 2 r,j r 2 + ω 2 z,j (z − z j ) 2 ] with radial and axial angular frequencies, ω r,j = 2πν r,j and ω z,j = 2πν z,j , respectively, with ν r = 119 Hz (178 Hz) and ν z = 166 Hz (248 Hz) for 87 Rb ( 39 K) atoms-labeled as species 1 (2)-in such a way that m 1 ω 2 r,1 = m 2 ω 2 r,2 (and similarly for the axial frequencies). For atomic species with different masses, there is a gravitational sag between the two trap centers, i.e., |z 1 − z 2 | > 0. Since this sag can be eliminated, we treat the sag z 1 − z 2 as a variable in our modeling. The scattering length of 87 Rb is fixed at a 11 = 99 a 0 .

IV. EQUILIBRIUM PHASE DIAGRAM
At equilibrium, both the source terms and the collision integrals vanish, hence we can set φ j (r,t) = φ j (r)e −iμ j t/ with chemical potential μ j in Eq. (1) and use a semiclassical Hartree-Fock approximation for the thermal cloudñ j (r) = g 3/2 ( j )/λ 3 j with thermal wavelength λ j = 2π 2 /(m j k B T ) and local fugacity j = exp[(μ j − U j n )/(k B T )] at temperature T . We also set the sag to be zero. We then obtain the density profiles by solving Eq. (1) (using imaginary-time propagation) andñ j (r) self-consistently.
Focussing initially on the zero-temperature case, we consider four different values (3 × 10 4 , 5 × 10 4 , 10 5 , and 5 × 10 5 ) for the total number of 87 Rb atoms (N Rb ) and 39 K atoms (N K ), resulting in a total of 16 combinations. The ground states, shown in Figs. 1(e)-1(i), are then given by the final density profiles with the lowest total energy, probed numerically by means of different initial states [44].
The ambiguity in classifying a trapped mixture into a mixed or demixed phase becomes evident in Fig. 1(h). Based on the usual prescription (motivated by the homogeneous system), this would typically be classified as a mixed phase, since = 0.5. To avoid such an evident inconsistency, we instead propose characterizing the mixture (with symmetrical distribution) by the difference in the normalized trap-center density, If the two components overlap strongly [e.g., Fig. 1  To the lowest order we neglect the effect of 87 Rb on 39 K, hence we can adopt the single-component Thomas-Fermi approximation, n c,2 ≈ [μ 2 − m 2 2 (ω 2 r,2 r 2 + ω 2 z,2 z 2 )]/g 22 . Taking into account the mean field g 12 n c,2 , 87 Rb atoms experience an effective trap potential (1 − g 12 /g 22 )m 1 (ω 2 r,1 r 2 + ω 2 z,1 z 2 )/2 at the trap center. If g 12 /g 22 < 1, the effective trap remains harmonic, hence 87 Rb has a peak at the center. However, if g 12 /g 22 > 1, the effective trap turns into a potential barrier and 87 Rb develops a shell structure. A complete phase diagram showing data for all 16 combinations of N Rb and N K can be found in the Appendix (Fig. 6) to further support our reasoning. We also show there (Fig. 7) a similar observation for a phase diagram of an isotropic trap which rules out trap anisotropy as the underlying origin of our observations. We next investigate how this boundary is changed at finite temperature. Consider a line scan in the g 12 /g 22 -g 11 /g 12 space that corresponds to tuning the scattering lengths a 12 and a 22 in the experimental setup [14] by using Feshbach resonances. This is shown as the white circles in Fig. 1(a). We simulate the density profiles at T = 100 nK 0.4 T c and 150 nK 0.6 T c for N Rb = (6-8) × 10 4 , and N K = (3-4) × 10 4 , chosen such that the number of BEC atoms does not deviate significantly from 5 × 10 4 ( 87 Rb) and 3 × 10 4 ( 39 K) respectively [Figs. 2(a) and 2(b)], used at T = 0 nK in Fig. 1(a). We find the equilibrium condensate densities to be very similar for the three different temperatures for the same miscibility parameter [Figs. 2(c)-2(h)], even though the thermal clouds have different magnitudes. As a result, n norm will be approximately the same even at different temperatures [ Fig. 3(a)] as long as the condensate atom numbers are comparable.

V. DETECTION THROUGH DIPOLE OSCILLATION
The question arises whether n norm is a physically meaningful quantity. To investigate this, we dynamically simulate the dipole oscillation of the mixture at different miscibilities and temperatures. Starting with a mixture at equilibrium, we excite the dipole mode in our simulation by increasing the trap-center separation from 0 μm to 0.2 μm in 5 ms linearly in time, followed by a decrease in the separation back to 0 μm in 1 ms linearly in time. The separation is chosen to be small compared to the Thomas-Fermi radii.
In Fig. 4, we show the simulated displacements of the condensates. The top panels display the undamped oscillation at zero temperature for an immiscible mixture (left) and a rapidly damped oscillation for a miscible mixture (right). At finite temperature (middle panels), numerical solution of our full binary kinetic model [37] (including all possible collisional processes) reveals damping even for the immiscible mixtures, due to the interaction with the thermal cloud. These numerical results are compared for different conditions by fitting the COM of the 39 K atoms with a damped sinusoidal function for time t > 10 ms [45], as z = A exp(−γ t) cos(2πνt + ϕ). Both the damping rate γ and the frequency ν, shown in Figs. 3(b) and 3(c), respectively, increase markedly at ≈ 0.2 (vertical dashes), at which n norm also starts to change in Fig. 3(a). We have checked the results for other atom numbers and arrived at the same conclusion (e.g., critical ≈ 0.5 if N Rb = 5 × 10 4 and N K = 10 5 , or critical ≈ 0 if N Rb = 5 × 10 5 and N K = 3 × 10 4 ).
Interestingly, we find that the mean-field interactions between the BECs and the thermal clouds are by far the dominant damping mechanisms; nonetheless, inclusion of collisions appears essential at long time scales ( γ −1 ) to eliminate residual center-of-mass oscillations, as evident by the 39 K COM oscillations shown in the bottom panels of Fig. 4  (and associated inset).
It is noteworthy that the sharp changes in both quantities remain detectable at finite temperature (green squares and red diamonds in Fig. 3), albeit with a decrease in the difference between the damping rates across the transition. This indicates the relevance of these changes in a realistic experimental setup. In addition, the much lower oscillation frequency compared to any of the axial trap frequencies (166 and 248 Hz) highlights the impact of interspecies repulsion [46].
Overall, our results in Fig. 3 demonstrate that n norm is physically meaningful, and we can detect the transition by monitoring both the damping rates and the frequency. Interestingly, n norm also reflects the penetrability of the mean-field barrier. Consequently, we can reinterpret our results as related to the distinct dynamical behavior of the BEC in the presence of a penetrable or impenetrable barrier, as seen also in recent experiments on vortex generation [47] and Josephson junction [48] in single-component Bose gas.
It is interesting to note how counterintuitive the underdamped oscillation of an immiscible mixture is. Since an immiscible mixture is associated with strong interspecies repulsion, we would naively expect the dipole oscillation to be quickly damped by the strong interaction. Our numerical results reveal the opposite.
This can be explained based on the different dynamical excitations associated with the topology of the equilibrium densities. In the immiscible case, there is a large interfacial tension between the two condensates [49]. Since a large amount of energy is needed to break the surface tension, the lower energy excitation involves undamped surface shape oscillations [50]. In the top panels of Fig. 5, we show the column-integrated density of the 87 Rb condensate, n int (y,z) = dx n c,1 ( x 2 + y 2 ,z), at different times. These snapshots clearly reveal the surface oscillation, where the 87 Rb cloud encircles the 39 K cloud and periodically displays a gap at the interface.
For the miscible case, large spatial overlap of the condensates favor the counterflow instability [29]. Sound, seen as density modulation in the bottom panels of Fig. 5, can be generated if the relative speed v of the two condensates is greater than c 1 + c 2 [51], where c i = g ii n c,i (0,0)/m i is the speed of sound of each individual condensate. This has been probed experimentally with a double superfluid Bose-Fermi mixture executing dipole oscillations in an elongated trap [17,53] via a different excitation scheme. For our parameters, v ≈ (ω z,1 + ω z,2 ) × 0.1 μm ≈ 0.3 mm/s, damping emerges due to density inhomogeneities even though v is much smaller than the trap-center speeds of sound, c 1 ≈ 2.2 mm/s and c 2 ≈ 3.4 mm/s. In addition, the increase in damping rate with decreasing for > 0.5 [ Fig. 3(b)] is consistent with the counterflow instability [54].

VI. EXPERIMENTAL FEASIBILITY
The simulations presented so far apply to many existing multispecies experiments, e.g., the mixture of 87 Rb and 39 K [14]. Typically, however, the gravitational sag between the two components of different mass is not compensated and experiments are conducted with relatively small spatial overlap [10,13,14,16]. Depending on the trap parameters the differential sag can exceed the typical Thomas-Fermi radius and therefore needs to be compensated to conduct the experiments proposed above.
The gravitational sag for each of the species is given by g/ω 2 z,j , where g is the acceleration due to gravity. The differential sag between both species can be canceled by employing an optical potential with a carefully selected wavelength. One option is to select a dipole trap which leads to identical frequencies ω z,j for both species [55]. Another option is to add an additional optical potential to an existing trap which introduces an upward force proportional to the atomic mass, thus canceling the sag.
In both cases the necessary detuning δ j of the dipole potential can be found as follows. In the simplified case of a two-level system the dipole force is proportional to j /ν 3 0,j δ j , where j is the linewidth and ν 0,j is the transition frequency [2]. Hence the same upward acceleration for both species can be obtained if the criterion m 2 /m 1 = 2 ν 3 0,1 δ 1 / 1 ν 3 0,2 δ 2 is met. The detunings δ j then provide the wavelength necessary to cancel gravity. A full calculation for the 87 Rb -39 K mixture shows that an additional beam at 806 nm should hence be used to realize the second cancellation method outlined above.

VII. CONCLUSIONS
We have provided numerical evidence that for a trapped condensate mixture with overlapping trap centers, the miscible-immiscible transition depends critically on the condensate numbers, deviating from the simple homogeneous prediction used to date. We demonstrate that this transition can be mapped out experimentally by measuring the damping rate and the frequency of the dipole oscillations, the predominant contribution to which stems from mean-field coupling. We relate this change in damping rate across the transition to the different dynamical excitations due to the topology of the initial density distribution. The successful implementation of the fully self-consistent dynamical kinetic model for binary mixtures opens up a multitude of possibilities for studying coupled binary dynamics, and we hope that our work will inspire future experiments on systematic studies of the dynamical behavior of trapped multicomponent condensates.
Data supporting this publication is openly available under an "Open Data Commons Open Database License" [56]. range of realistic parameters, the transition from miscible to immiscible does not happen precisely at the assumed = 0 line.