Miscibility regimes in a $^{23}$Na-$^{39}$K quantum mixture

Effects of miscibility in interacting two-component classical fluids are relevant in a broad range of daily applications. When considering quantum systems, two-component Bose-Einstein condensates provides a well controlled platform where the miscible-immiscible phase transition can be completely characterized. In homogeneous systems, this phase transition is governed only by the competition between intra- and inter-species interactions. However in more conventional experiments dealing with trapped gases, the pressure of the confinement increases the role of the kinetic energy and makes the system more miscible. In the most general case, the miscibility phase diagram of unbalanced mixtures of different atomic species is strongly modified by the atom number ratio and the different gravitational sags. Here, we numerically investigate the ground-state of a $^{23}$Na-$^{39}$K quantum mixture for different interaction strengths and atom number ratios considering realistic experimental parameters. Defining the spatial overlap between the resulting atomic clouds, we construct the phase diagram of the miscibility transition which could be directly measured in real experiments.

As for its classical counterpart, a mixture of two fluids is miscible if the fluids totally overlap forming a homogeneous solution or immiscible if the fluids remain phase-separated [23,24]. In the case of homogeneous quantum fluids, the miscible-immiscible phase transition is well defined and it is mediated by the miscibility parameter [25,26] where u 11 and u 22 are the intraspecies interaction coupling constants of species 1 and 2, respectively, and u 12 gives the interspecies interaction. This is an intuitive parameter based on the competition between intra-and interspecies interactions: if u 12 overcomes the intraspecies interaction terms (δ > 0), the fluids strongly repel each other making the system immiscible. On the contrary, if u 12 is smaller than the intraspecies interactions (δ < 0), the fluids overlap and the system is miscible. In such a picture, the miscibility regime of a two-component quantum gas can be controlled by varying the interaction coupling terms, which can be experimentally realized with the use of Feshbach resonances [27].
However, until very recent [28], homogeneous atomic BECs were not experimentally produced.
Instead, trapped atomic BECs exhibit an inhomogeneous density distribution as a result of the confinement. The increased role of the kinetic energy in such systems contributes to a more miscible mixture where phase-separation occurs for larger u 12 than the condition set by Equation 1. The shift at the miscible-immiscible critical point has been obtained in the case of mixtures composed of distinct hyperfine states of the same atomic species [29][30][31]. In the more broad scenario of unbalanced mixtures of different atomic species, the atom number ratio [32], the mass imbalance and the difference in trapping configurations between the components were also shown to affects the boundary of the miscibility phase transition [33][34][35][36]. The contribution of gravity, relevant for all real experiments due to the induced gravitational sag [37,38], is rarely taken into account in numerical simulations.
In this work, we perform numerical simulations of the ground-state of a two-component quantum mixture of 23 Na and 39 K atoms for different interaction strengths, according to the relevant Feshbach resonances for magnetic fields in the range of 95 − 117 G [14,39], in order to show the realistic miscibility regimes accessible in the experimental setup being developed in our laboratory [40] in the presence of gravity. We explore the effect of changing the number of atoms of the minority species ( 39 K), therefore changing the atom number ratio η, and calculate the spatial overlap between the atomic clouds as a quantity able to characterize the change in the miscibility regime of the system. The numerical simulations are performed at zero temperature, which satisfactory reproduces the experimental results for the case of strongly degenerate atomic mixtures [41], although theoretical works at finite temperature have shown a change of the miscibility condition of the system favoring phase separation [42][43][44].
The article is organized as follows. In Section II, we describe the two-component quantum gas at zero temperature in terms of a pair of coupled Gross-Pitaevskii equations (GPEs) (II A) and the numerical simulation method used to obtain the ground-state of the system (II B). In Section III, we first present our experimental system producing the 23 Na-39 K atomic mixtures (III A), followed by the results of the numerical simulation performed with realistic experimental parameters (III B) and the construction of the phase diagram of the miscible-immiscible transition for such a mixture (III C). Finally, in Section IV, we highlight our main findings and discuss some future perspectives for identifying the miscibility regime of a quantum mixture comparing with the results presented in this article.

II. METHODS
A. Description of an atomic quantum mixture Consider a mixture of two different bosonic atoms, labeled 1 and 2, at T = 0 in the weakly interacting regime where interactions are treated as contact interactions. Let N 1 and N 2 be the number of particles and φ 1 (r) and φ 2 (r) be the corresponding normalized single-particle wave functions. In such a picture, and neglecting terms of the order of 1/N 1 and 1/N 2 , the energy functional of the system [25,26,45] can be written as where m i (with i = 1, 2) is the mass of atomic species i, ϑ i (r) is the corresponding external potential, u ii = 4π 2 a ii /m i are the intra-species interaction terms and u 12 = 2π 2 a 12 /m 12 is the inter-species interaction term with m 12 = m 1 m 2 /(m 1 + m 2 ), the reduced mass of the system. For all relations, a ij is the associated two-body s-wave scattering length. The wave-functions ψ 1 (r) and ψ 2 (r) are the condensate wave-function of each atomic species, defined as ψ 1 (r) = N 1 φ 1 (r) and ψ 2 (r) = N 2 φ 2 (r).
Minimizing the energy functional of Equation 2 under the constraint of fixed number of particles, N 1 and N 2 , one obtains the time-independent coupled Gross-Pitaevskii equations where µ 1 and µ 2 are the chemical potential of atomic species 1 and 2, respectively. If the interspecies interaction vanishes (u 12 = 0), Equations 4 and 5 are no longer coupled and each species behave as a single species atomic cloud. In this case, approximations such as the Thomas-Fermi approximation [25,26,45], for which the kinetic term of the GPE is neglected, can be used to find a solution for the ground-state of the system. On the other hand, when u 12 = 0, the competition between inter-and intraspecies interactions gives rise to a phase transition from a miscible to an immiscible (phase-separated) phase when increasing the positive inter-species interaction strength.
The existence of overlapping and non-overlapping regions between the atomic clouds dramatically changes the ground-state configuration of the system and it is not always possible to find analytical solutions for it, even relying on approximations [46]. A more powerful technique to obtain the ground-state of a trapped two-component BEC makes use of a numerical simulation with imaginary time evolution of the coupled GPEs.

B. Numerical simulation of the ground-state
The numerical simulation used to obtain the ground-state of the two-species BEC consists of projecting onto the minimum of the GPEs each initial trial state by propagating them in imaginary time [47]. To describe the method, let us first consider a system described by a Hamiltonian H for which the time evolution of one of its eigenstates, ψ n (r, t) with Hψ n (r, 0) = E n ψ n (r, 0), is easily obtained as: where E n is the energy associated with the n-eigenstate. The time evolution of an arbitrary trial function Ψ(r, 0), written as a linear combination of the system's eigenstates, is simply given by If ones calculates Ψ(r, t) for t = −iτ , the complex exponentials in Eq. 7 are replaced by exponential decays with decay constants given by E n / . By evaluating Ψ(r, t) at different time steps ∆τ with τ = ξ∆τ , Ψ(r, τ ) → ψ 0 (r, τ ), the ground-state of the system. The exact convergence is only obtained when τ → ∞, however, convergence methods based on the variation of the total energy of the system are used to set an upper limit for τ .
In the numerical simulations performed in this work, we define a trial function Ψ i (r, 0) for each species i with time evolution given by where H i ψ i (r) = µ i ψ i (r) from Equations 4 and 5. Considering t = −i∆τ with ∆τ infinitesimal, the resulting exponential can be expanded in a Taylor series and the time evolution of Ψ i (r, t + ∆τ ) is given by In order to achieve sufficient long times in the simulations let it be t final = ξ∆t, with ξ being an integer, Equation 9 is calculated ξ times. The resulting wave-function obtained after each time step is normalized in order to preserve the atom number.

III. RESULTS
The numerical simulations performed in this work are done following the parameters of the experimental setup being developed in our laboratory. For this reason, we first start the Results Section, subsection III A, with a description of the experimental setup and its current status in the preparation of a two-species BEC of 23 Na and 39 K. Later, the results from the numerical simulations are presented and discussed in the following two subsections.

A. Experimental setup
A complete description of the experimental setup and experimental sequence for producing a Bose-Einstein condensate of 23 Na atoms is described in [40]. Here, we present a short description of the system giving the experimental parameters relevant for the simulations performed later in this Section.
Briefly, sodium and potassium atoms coming from independent two-dimensional magnetooptical traps (2D-MOTs) [48,49] are combined in a common vacuum chamber where they will be trapped and further cooled in a three-dimensional MOT (3D-MOT). Due to the strong interspecies losses present in the Na-K mixture [40,50], the operation of an intial two-color MOT is not the best alternative in our experiment. Instead, we chose to favor the minority species (potassium) during the MOT phase, starting the MOT sequence with the loading of a single species MOT of 39 K until it reaches the saturation value (∼ 20 s). Next, we operate the two-color MOT by switching on the lights responsible for trapping and cooling sodium atoms. We control the initial atom number ratio N 0 Na /N 0 K by changing the time duration of the loading of the sodium atoms in the two-color MOT operation.
Once the two species are loaded, we perform subsequent cooling procedures followed by a fine pumping stage which transfer both species to the F = 1 ground-state before turning on an optically plugged Quadrupole trap [51]. At the beginning of the magnetic trap the atomic clouds have N Na ∼ 1 × 10 9 atoms and N K ∼ 1 × 10 6 atoms both at T = 220 µK trapped in the Evaporative cooling [52] of sodium is done with microwave radiation at ∼ 1.7 GHz while potassium atoms are sympathetic cooled [53,54] decreasing its temperature without significant atom loss. At T ∼ 6 − 7 µK, the atomic clouds are transferred to a pure optical dipole trap (ODT) [55] where the interspecies interaction can be tuned with the use of Feshbach resonances [27] by applying a uniform magnetic field. We have atomic clouds with N Na = 5 × 10 6 and N K = 8 × 10 5 at the beginning of the ODT for maximum atoms number of 39 K. In single-species operation for sodium under the same conditions we obtain an almost pure BEC (with BEC fraction > 80%) with N = 1 × 10 6 atoms at T ∼ 80 µK after applying an optical evaporation which reduces the initial ODT potential height by a factor of five in 4.2 s. The final ODT configuration exhibit a planar geometry with equal frequencies in the xy-plane perpendicular to the gravity direction. The final frequencies are ω x,y = 2π × 107(137) Hz and ω z = 2π × 148(193) Hz for Na(K), respectively. This is the actual situation of our experimental system and, following the initial atom number difference in the ODT, we estimate to be able to obtain a two-species BEC once implemented the Feshbach field. Following these experimental numbers we performed the simulations described in Section II B which results are presented in the following. where ω 1 = (ω x ω y1 ω z1 ) 1/3 is the geometric mean of the trapping frequencies for species 1. We considered species 1 (2) as the potassium (sodium) atoms. We apply convergence methods based on the difference between the wave-functions of subsequent time intervals and monitor the total energy evolution in order to ensure the achievement of the ground-state configuration for both species. With these methods, typical integration times gave t final ∼ 3000.
The number of sodium atoms was chosen N Na = 5 × 10 5 atoms in agreement with the numbers obtained in the experiment. The number of potassium atoms was varied with N K = 1 × 10 4 − 5 × 10 5 atoms setting η = N Na /N K = 50 − 1. The trapping frequencies were also set from the experimental values with ω x,y = 2π×107(137) Hz and f z = 2π×148(193) Hz for Na(K), respectively.
The sodium scattering length was fixed to a Na = 52 a 0 , with a 0 being the Born radius, while the scattering length of 39 K, a K , and the interspecies scattering length, a NaK , was varied according to the Feshbach resonances occuring at magnetic fields smaller than 300 G [14,39]. In Figure 1, we show the values of the scattering lengths (a Na , a K and a NaK ) as a function of the magnetic field in the region with B = 95 − 117.2 G. In this region, both a K and a NaK are positive and the system changes its behaviour from immiscible to miscible with increasing the magnetic field. The predicted phase transition point for a homogeneous system (with δ = 0) occurs at B 0 = 109.1 G [14]. The potassium scattering length was obtained with the simple relation: and ∆ 1 = 55 G and ∆ 2 = −37 G are the corresponding resonance widths [39]. The a NaK curve displayed in Figure 1 was obtained from [14] via a private communication.
Due to the presence of the gravitational force, each species suffers a different gravitation sag and the phase-separation at the immiscible phase occurs along the vertical direction (z-axis). In For a fixed miscibility parameter we observe different behaviours of the system when changing η. In Figure 2 (a), the system is always immiscible, i.e., the sodium and potassium atoms do not share the same position in the trap for any value of η. In Figure 2 (b), the system is expected to be immiscible according to the miscibility parameter (δ = 0.54 > 0), however for η = 50 the phase-separated region disappears and the potassium atoms always share its position in the trap with sodium atoms, which is characteristic of a miscible system. Finally, in Figure 2 (c), the system is expected to be miscible with δ = −0.30 < 0 but in the case of η = 5 there is still a region of the potassium cloud that do not shares the trap with sodium atoms remaining immiscible.
We see that, for inhomogeneous systems, the atom number ratio has a strong influence in the miscibility regime of a two-species Bose-Einstein condensate. While in the homogeneous case, the miscibility parameter is enough to set the regime The construction of a phase diagram miscible-immiscible needs a more quantitative way of defining the miscibility region of a given set of parameters for the two-species system. Proposals to characterize the regime of such system include the calculation of the binder cumulant of the system's magnetization [56], the difference between the centers of mass of each atomic cloud [32], the study of the entropy of the mixture as defined in [44,57] and the monitor of dipole oscillations of the atomic clouds in a harmonic trap [35]. Here, similar to the works presented in [29,58], we follow the definition of the miscible and immiscible phases and propose the calculation of the spatial overlap between the atomic clouds to be an indicator of the phase transition.
We define the spatial overlap of the atomic clouds as: where n 1 (r) and n 2 (r) are the atomic densities of species 1 and 2, respectively.
In Fig. 3 to increase for B > 111 G, a magnetic field larger than B 0 . The dashed vertical line in Figure 3 represents the critical point for the transition at B 0 with δ = 0.
To define the transition from immiscible to miscible from the normalized overlap we associate a threshold-like behaviour and identify the Feshbach field (B peak ) for which O norm varies the most. This is done performing the numerical second derivative of the normalized overlap and identifying its maximum value. The second derivatives as a function of the Feshbach field for all η are displayed in Figure. 4. The maximum value of the curves drifts to larger magnetic fields as η decreases. In Fig. 4 (b), we show B peak as a function of η. The almost linear behaviour of the points in the semilog scale suggests a dependence of B peak (η) = B 0 peak − α ln η. We find B 0 peak = 113.5 G and α = 1.70. Differently from earlier works performed with a balanced mixture of two distinct hyperfine states of a single species [29][30][31], in our case for η = 1 the system is more immiscible and the miscible phase only occurs for δ ∼ −0.5. A deeper analysis of the η = 1 case is presented in Fig. 6, where we show the normalized overlap for η = 1 under different trapping conditions: real experimental conditions (star), without gravity (asterisk), considering equal trapping potentials with ϑ 1 = ϑ 2 (plus sign) and for the homogeneous case (gray solid curve) obtained setting the external potentials However, due to the large difference between the intraspecies scattering lengths for sodium and potassium (a Na = 52a 0 and a K ∼ 7.59 − 8.73a 0 ), the system still behaves more immiscible than The identification of the miscibility regime of the Na-K mixture under realistic experimental conditions is important when defining the best parameters for studying different physical phenomena. In studies which the spatial overlap between the components of the mixture is important (i.e. coupled vortex dynamics [59][60][61], binary quantum turbulence [62], coupled superfluidity and excitations [63,64], etc.) it is not always sufficient to have δ < 0. The contrary is also true, when the immiscible nature of the system is relevant (i.e. in studies of dynamical instabilities [65][66][67][68]), δ > 0 is not always sufficient, specially in the case of large atom number imbalance between the atomic species.

IV. DISCUSSION
We have shown that the miscible-immiscible phase transition in a trapped two-component Bose-Einstein condensation of different atomic species under realistic experimental parameters (considering the effect of gravity and different trapping potentials) suffers strong influence of the atom number ratio η. In the case of large η, the system behaves more miscible than the homogeneous case with the transition occurring at δ > 0, while for η = 1, the system is more immiscible with the transition occurring at δ < 0. We have defined the miscibility regime of the system by identifying the magnetic field B peak for which the normalized spatial overlap between the atomic clouds changes the most. This value was obtained from the magnetic field for which the numerical second derivative of the normalized overlap exhibits a maximum. The behaviour of B peak with η could be easily associated with a logarithm dependence from the graph of Fig. 4 (b) making it possible to draw the critical curve in the phase diagram of the miscible-immiscible phase transition for the simulated 23 Na-39 K quantum mixture (see Fig. 5). The use of the spatial overlap to identify the miscibility regime of the system could be directly implemented in real experiments by performing high resolution in situ images of each atomic species. Further characterizations both on the experimental and theoretical sides could be performed using dynamical properties of the atomic mixture, such as the dipole oscillations proposed in [35], and considering finite temperature effects as realized in recent works [42][43][44].
V. ACKNOWLEDGEMENTS