Miscibility regimes in a ²³Na-³⁹K 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 $\phi_{1}(\mathbb{r})$ and $\phi_{2}(\mathbb{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 [45, 25, 26] can be written as

where $m_{i}$ (with $i=1,2$) is the mass of atomic species $i$, $\vartheta_{i}(\mathbb{r})$ is the corresponding external potential, $u_{ii}=4\pi\hbar^{2}a_{ii}/m_{i}$ are the intra-species interaction terms and $u_{12}=2\pi\hbar^{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 $\psi_{1}(\mathbb{r})$ and $\psi_{2}(\mathbb{r})$ are the condensate wave-function of each atomic species, defined as

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 $\mu_{1}$ and $\mu_{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  [45, 25, 26], 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}\neq 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.

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, $\psi_{n}(\textbf{r},t)$ with $H\psi_{n}(\textbf{r},0)=E_{n}\psi_{n}(\textbf{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 $\Psi(\textbf{r},0)$, written as a linear combination of the system’s eigenstates, is simply given by

If ones calculates $\Psi(\textbf{r},t)$ for $t=-i\tau$, the complex exponentials in Eq. 7 are replaced by exponential decays with decay constants given by $E_{n}/\hbar$. By evaluating $\Psi(\textbf{r},t)$ at different time steps $\Delta\tau$ with $\tau=\xi\Delta\tau$, $\Psi(\textbf{r},\tau)\rightarrow\psi_{0}(\textbf{r},\tau)$, the ground-state of the system. The exact convergence is only obtained when $\tau\rightarrow\infty$, however, convergence methods based on the variation of the total energy of the system are used to set an upper limit for $\tau$.

In the numerical simulations performed in this work, we define a trial function $\Psi_{i}(\textbf{r},0)$ for each species $i$ with time evolution given by

where $\widehat{H}_{i}\psi_{i}(\textbf{r})=\mu_{i}\psi_{i}(\textbf{r})$ from Equations 4 and 5. Considering $t=-i\Delta\tau$ with $\Delta\tau$ infinitesimal, the resulting exponential can be expanded in a Taylor series and the time evolution of $\Psi_{i}(\mathbb{r},t+\Delta\tau)$ is given by

In order to achieve sufficient long times in the simulations let it be $t_{\text{final}}=\xi\Delta t$, with $\xi$ being an integer, Equation 9 is calculated $\xi$ times. The resulting wave-function obtained after each time step is normalized in order to preserve the atom number.

A complete description of the experimental setup and experimental sequence for producing a Bose-Einstein condensate of ²³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 magneto-optical 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 [50, 40], 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 ³⁹K until it reaches the saturation value ($\sim 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_{\text{Na}}^{0}/N_{\text{K}}^{0}$ 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_{\text{Na}}\sim 1\times 10^{9}~{}$atoms and $N_{\text{K}}\sim 1\times 10^{6}~{}$atoms both at $T=220~{}\mu$K trapped in the $|F=1,m_{F}=-1>$ hyperfine ground-state.

Evaporative cooling [52] of sodium is done with microwave radiation at $\sim 1.7~{}$GHz while potassium atoms are sympathetic cooled [53, 54] decreasing its temperature without significant atom loss. At $T\sim 6-7~{}\mu$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_{\text{Na}}=5\times 10^{6}$ and $N_{\text{K}}=8\times 10^{5}$ at the beginning of the ODT for maximum atoms number of ³⁹K. In single-species operation for sodium under the same conditions we obtain an almost pure BEC (with BEC fraction $>80\%$) with $N=1\times 10^{6}~{}$atoms at $T\sim 80~{}\mu$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 $\omega_{\text{x,y}}=2\pi\times 107(137)~{}$Hz and $\omega_{\text{z}}=2\pi\times 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.2 which results are presented in the following.

The ground-state of  ²³Na - ³⁹K mixtures was obtained with the numerical simulation method described in II.2. In the simulations, we discretize the space with a three-dimensional grid of $69\times 69\times 69$. The grid step size was chosen to be equal to $0.6~{}\mu$m resulting in a total volume of $41\times 41\times 41~{}\mu$m³. The time interval for the simulations were $\Delta t=50\times 10^{-6}$ in units of $1/\overline{\omega}_{1}$, where $\overline{\omega}_{1}=(\omega_{x}\omega_{y1}\omega_{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_{\text{final}}\sim 3000$.

The number of sodium atoms was chosen $N_{\text{Na}}=5\times 10^{5}~{}$atoms in agreement with the numbers obtained in the experiment. The number of potassium atoms was varied with $N_{\text{K}}=1\times 10^{4}-5\times 10^{5}~{}$atoms setting $\eta=N_{\text{Na}}/N_{\text{K}}=50-1$. The trapping frequencies were also set from the experimental values with $\omega_{\text{x,y}}=2\pi\times 107(137)~{}$Hz and $f_{\text{z}}=2\pi\times 148(193)~{}$Hz for Na(K), respectively. The sodium scattering length was fixed to $a_{\text{Na}}=52~{}a_{0}$, with $a_{0}$ being the Born radius, while the scattering length of ³⁹K, $a_{\text{K}}$, and the interspecies scattering length, $a_{\text{NaK}}$, was varied according to the Feshbach resonances occuring at magnetic fields smaller than $300~{}$G [39, 14]. In Figure 1, we show the values of the scattering lengths ($a_{\text{Na}}$, $a_{\text{K}}$ and $a_{\text{NaK}}$) as a function of the magnetic field in the region with $B=95-117.2~{}$G. In this region, both $a_{\text{K}}$ and $a_{\text{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 $\delta=0$) occurs at $B_{0}=109.1~{}$G [14]. The potassium scattering length was obtained with the simple relation:

where $a_{\text{bg}}=-19a_{0}$ is the background scattering length, $B_{01}=32.6~{}$G and $B_{02}=162.8~{}$G are the position of the first and second resonances for ³⁹K at the $|F=1,m_{F}=-1\rangle$ hyperfine state and $\Delta_{1}=55~{}$G and $\Delta_{2}=-37~{}$G are the corresponding resonance widths [39]. The $a_{\text{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 Figure 2, we show the simulated density profiles along the $z$-axis of ²³Na (in blue) and ³⁹K (in red) for Feshbach fields $B=100~{}$G in (a) with $\delta=4.32$, $B=108~{}$G in (b) with $\delta=0.22$ and $B=111~{}$G in (c) with $\delta=-0.50$ and different atom number ratio $\eta=50,10,5$ in solid, dashed and dotted lines, respectively.

For a fixed miscibility parameter we observe different behaviours of the system when changing $\eta$. 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 $\eta$. In Figure 2 (b), the system is expected to be immiscible according to the miscibility parameter ($\delta=0.54>0$), however for $\eta=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 $\delta=-0.30<0$ but in the case of $\eta=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 of the system, in most real experiments where the condensates are confined by a harmonic trap, additional information is necessary to establish the critical point for the miscible-immiscible phase transition.

The ground-state configurations obtained in the previous section show the flexibility of the ²³Na - ³⁹K mixture in achieving different miscibility regimes with the change of the Feshbach magnetic field $B$ and the atom number ratio $\eta$.

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}(\mathbb{r})$ and $n_{2}(\mathbb{r})$ are the atomic densities of species $1$ and $2$, respectively.

In Fig. 3, we show the spatial overlap normalized by the overlap, $O_{0}=O(a_{12}=0)$, at each case as a function of the Feshbach field for different atom number ratio, $\eta$. The normalization was done considering that the case of vanishing inter-species interaction exhibit the maximum spatial overlap between the atomic clouds possible in each configuration. In the immiscible region ($B<B_{0}=109.1~{}$G), $O_{\text{norm}}$ exhibit small values with $O_{\text{norm}}\rightarrow 0$ when reducing the magnetic field for all $\eta$. Approaching $B_{0}$, the spatial overlap increases differently for each $\eta$ with larger $\eta$ showing an earlier increase on the spatial overlap. The limit case of $\eta=50$ shows a significant increase of $O_{\text{norm}}$ for $B>104~{}$G, increasing over a broad range of magnetic fields. This increase occurs before the critical point estimated with the miscibility parameter, in accordance with was already observed in the density profiles of Fig. 2. Reducing $\eta$ increases the magnetic field for which the spatial overlap significantly increases. For the other limit with $\eta=1$, the spatial overlap starts 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 $\delta=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_{\text{peak}}$) for which $O_{\text{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 $\eta$ are displayed in Figure. 4. The maximum value of the curves drifts to larger magnetic fields as $\eta$ decreases. In Fig. 4 (b), we show $B_{\text{peak}}$ as a function of $\eta$. The almost linear behaviour of the points in the semilog scale suggests a dependence of $B_{\text{peak}}(\eta)=B_{\text{peak}}^{0}-\alpha\ln{\eta}$. We find $B_{\text{peak}}^{0}=113.5~{}$G and $\alpha=1.70$.

The miscible-immiscible phase diagram for the ²³Na-³⁹K mixture under our experimental conditions is shown in Figure 5. The colormap represents the value of the normalized overlap for each combination of $\eta$ and $B$ ranging from zero to unit. The transition point for each $\eta$ obtained from Figure 4 (b) is displayed by the solid light gray curve with the shaded area covering its uncertainty. The black dashed line represents the transition point for the homogeneous case setting $\delta>0$ to the left side of the curve and $\delta<0$ to the right.

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 $\eta=1$ the system is more immiscible and the miscible phase only occurs for $\delta\sim-0.5$. A deeper analysis of the $\eta=1$ case is presented in Fig. 6, where we show the normalized overlap for $\eta=1$ under different trapping conditions: real experimental conditions (star), without gravity (asterisk), considering equal trapping potentials with $\vartheta_{1}=\vartheta_{2}$ (plus sign) and for the homogeneous case (gray solid curve) obtained setting the external potentials $\vartheta_{1}(\mathbb{r})=\vartheta_{2}(\mathbb{r})=0$. The dashed vertical line in black represents the point for $\delta=0$ which indeed matches the abrupt transition of $O_{\text{norm}}$ from $0$ to unity observed in the homogeneous case. The role of gravity and different trapping configurations for each species is clear in the data of Figure 6: setting $g=0$ and $\vartheta_{1}=\vartheta_{2}$ drifts the transition point to smaller magnetic fields approaching $B_{0}$. However, due to the large difference between the intraspecies scattering lengths for sodium and potassium ($a_{\text{Na}}=52a_{0}$ and $a_{\text{K}}\sim 7.59-8.73a_{0}$), the system still behaves more immiscible than the homogeneous case.

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 $\delta<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]), $\delta>0$ is not always sufficient, specially in the case of large atom number imbalance between the atomic species.