Universality of miscible-immiscible phase separation dynamics in two-component Bose-Einstein condensates

We investigate the non-equilibrium dynamics across the miscible-immiscible phase separation in a binary mixture of Bose-Einstein condensates. The excitation spectra reveal that the Landau critical velocity vanishes at the critical point, where the superfluidity spontaneously breaks down. We analytically extract the dynamical critical exponent $z=2$ from the Landau critical velocity. Moreover, by simulating the real-time dynamics across the critical point, we find the average domain number and the average bifurcation delay show universal scaling laws with respect to the quench time. We then numerically extract the static correlation length critical exponent $v=1/2$ and the dynamical critical exponent $z=2$ according to Kibble-Zurek mechanism. The scaling exponents $(v=1/2, z=2)$ in the phase separation driven by quenching the atom-atom interaction are different from the ones $(v=1/2, z=1)$ in the phase separation driven by quenching the Rabi coupling strength [PRL \textbf{102}, 070401 (2009); PRL \textbf{107}, 230402 (2011)]. Our study explores the connections between the spontaneous superfluidity breakdown and the spontaneous defect formation in the phase separation dynamics.

In this paper, we investigate the non-equilibrium dynamics of phase separation in a binary mixture of atomic BECs. When the system is driven across the critical point at a finite rate, the critical dynamics across a miscibleimmiscible phase transition is studied. Through calculating the Bogoliubov excitation spectrum, we find that the Landau critical velocity vanishes and the superfluidity breaks down at the critical point, and we analytically extract the dynamical critical exponent and static correlation length critical exponent from the Landau critical velocity. To show how the non-equilibrium dynamics appears, we numerically simulate the real-time dynamics in quench process, in which the intra-component interaction strength is linearly swept through the critical point. From the non-equilibrium dynamics far from the critical point, we numerically extract two universal scalings for the average domain number and the average bifurcation delay with respect to the quench rate, and we find the critical exponents derived from the numerical results consist well with the analytical exponents. These scaling exponents in the phase separation induced by tuning the atom-atom interaction are different from the ones [5][6][7][8], in which the phase separation is induced by tuning the Rabi coupling strength.
The paper is organized as follows. In section 2, we describe the model and discuss its ground states. In section 3, we implement the Bogoliubov-de Gennes (BdG) analysis and obtain the Landau critical velocity. In section 4, we analytically extract the Kibble-Zurek (KZ) scalings and numerically simulate the real-time dynamics of the phase transition and extract the universal scalings. Finally, we give a brief summary and discussion in section 5.

Model
We consider a mixture of two weakly interacting BECs in a ring trap. When the transverse frequency ω ⊥ is sufficiently strongly, we can integrate out the transverse degrees of freedom and obtain the one-dimension ) ( ) characterize the intra-and intercomponent interactions, with m j being the mass of a atom in component j, a jj and a 12 respectively denoting the intra-component and inter-component s-wave scattering lengths. In experiments, g jj and g 12 can be adjusted by using Feshbach resonance [52][53][54]. For simplicity, we only consider the homogeneous case, namely the external trapping potential V(x)=0.
In the mean-field (MF) theory, the system obeys the Gross-Pitaevskii equations (GPEs) The nature of the ground states is determined by the competition between the intra-component and intercomponent interactions [44][45][46]. If the intra-component interaction dominates, i.e. g g g 11 22 12 2 > , the energy is minimized when the two components occupy all available volume. In such a miscible phase, the two BEC wavefunctions coexist at all positions. However, if the inter-component interaction dominates, i.e. g g g 11 22 12 2 < , the energy of the system is minimized when the two BECs are separated in space. In such an immiscible phase, the two BEC wave-functions occupy different spatial regions. For the ground state in the immiscible phase, it is not easy to obtain the exact solution analytically. However, by utilizing the imaginary time propagation method, we numerically solve equation (3) with a split-step method [55] with the Wick rotation t=−iτ and obtain the ground state under the period boundary condition. In figure 1, we show a typical ground state in the immiscible phase. The two BEC wave-functions indeed occupy different space. We observe an increase of the interface with the intra-component interaction strength g jj . If the intra-component interaction strength g jj is stronger than a threshold g jj c (or the inter-component interaction g 12 is smaller than a threshold g c 12 ), the interface will cover the whole space, which is referred as miscible-immiscible phase transition.
To understand the critical dynamics near the phase transition, we introduce a dimensionless distance = is the critical point, and is a linearly quenched parameter with τ Q being the quench time. Below, we quench the intra-component interaction strength g 22 with different quench time τ Q over several orders of magnitude. Then, from our numerical simulation, we extract the relation between the number of topological defects, the bifurcation delay and the quench time τ Q .

Bogoliubov excitation and spontaneous superfluidity breakdown near the critical point
In this section, we analyse the Bogoliubov excitations near the critical point of phase separation and investigate the spontaneous breakdown of superfluidity. According to the Landau criterion, if the superfluid velocity is smaller than the Landau critical velocity, elementary excitations are prohibited due to the conservation of energy and momentum. However, around the critical point, the Landau critical velocity vanishes and so that elementary excitations appear spontaneously. We implement a Bogoliubov analysis to obtain the excitation modes over the ground states. In the miscible phase, the nonlinear Schrödinger equations have an obvious homogenous solution , with f j the ground state wave-function of the jth component and L the length of the system. The chemical potentials are given as μ 1 =g 11 ρ 1 +g 12 ρ 2 and μ 2 =g 22 ρ 2 +g 12 ρ 1 , and g 11 =g 12 in our system. To derive the Bogoliubov excitation spectrum, we consider the perturbed ground state Inserting the state(6) into the equation (3) and keeping the first-order terms, in units of , the linearized equations for the perturbations are given as, The perturbations δf 1,2 can be written as in which q is the excitation quasimomentum, ω is the excitation frequency, u j,q and v j,q , ( j=1, 2) are the Bogoliubov amplitudes. Inserting equation (9) into the linearized equation equation (7), (8) and comparing the coefficients for the terms of e iqx−iωt and e iqx+iωt , one obtains the BdG equations Assuming that the two BECs have the same particle number, N 1 =N 2 =N/2, we obtain the excitation spectrum by diagonalizing the matrix q ( ) for each q 2 , 12 with ò 0 =q 2 /2 and > , η − is a negative value, ω − becomes imaginary in the long-wavelength limit. This means that the phase transition corresponds to the dynamical instability of Bogoliubov excitations. In [5,6], the KZ mechanism originates from the gradually vanishing of the Higgs modes when the system approaches to the critical point, that is, the gapless Higgs modes at the critical point. However, in our work, the long-wavelength excitations are always gapless, and the KZ mechanism originates from the spontaneous superfluidity breakdown induced by the dynamical instability.
On the other hand, at the critical point, g g g g

KZ scalings
In this section, we extract two universal scalings from the Bogoliubov excitation in the miscible phase. As the softening of the phonon mode near the critical point, we can extract the dynamical critical exponent z from the Landau critical velocity. Furthermore, we can make use of the Landau critical velocity to define correlation length, and extract another universal critical exponent v. At the same time, we perform numerical simulations of non-equilibrium dynamics in the GPE(3). Starting from the ground state in the miscible phase, according to equation (5), we linearly sweep the intra-component interaction strength g 22 to drive the system across the critical point g c 22 . When the evolution time t increases from t<0 to t>0, the system goes from the miscible phase to the immiscible phase. Through analysing the domains and bifurcation delay, we numerically obtain two universal critical exponents, which are consist well with the analytical exponents.
The dynamical critical exponent z can be derived from the low-energy long-wave excitations. At the critical point, the low-energy excitation in the long-wavelength limit behave as 20,56], we have the dynamical critical exponent z=2. The static correlation length critical exponent v can be derived from the divergence of the correlation length ξ near the critical point. The Landau critical velocity v L provide a definition of the correlation length ξ according to mv L  x = ( ) [57]. Therefore, the Landau critical velocity v L should have a power-law scaling behavior around the critical point as In the vicinity of the critical point, we can simply expand η − to the leading order of ò, that is g g As v L h = -, thus we have the static correlation length critical exponent v 1 2 = .
In figure 2, we show our numerical results of v L for different  | |near the critical point in a log-log coordinate. Through linear fitting, we find a power-law v L v  µ | | with the static correlation length critical exponent v=0.4996±0.0001, which agrees well with the analytical one v=1/2.
In addition, one can also derive the critical exponent from the sound velocity [44]. In the long-wavelength limit (q → 0), the sound velocity reads We can expand c − to the leading order of ò, that is, c > , c 2 -becomes negative and so that the system becomes unstable. Therefore, the fastest growing modes with the wave number q c 2 f = -| |are given by minimizing 2 wwith respect to q, and these modes grow initially at a rate E c q 2 = -| | . We expect that these modes will set the length scale of the pattern that is formed, so that the condensates will separate into domain structure, and the length scale ξ of domain is given as At the same time, we can estimate the time scale τ for the formation process of domains Combining the equations (19) and (20), one can analytically obtain the two critical exponents v z 1 2, 2 = = ( ). In our simulations, we choose quench times τ Q over four orders of magnitude and perform 100 runs of simulation for each τ Q . Since the quantum fluctuations that trigger the phase transition are ignored in the MF approximation [50], we introduce random noises to the initial state and so that the dynamics of spontaneous topological defects can be studied by the MF theory. As the interaction strength g 22 is quenched according to equation (5), unlike the static bifurcation, in which the bifurcation exactly occurs at the critical point, the dynamical bifurcation takes place after the system crossing the critical point.
The bifurcation delay, b g g d c 22 22 | |ˆ|, is obtained by analyzing spatial fluctuation of the local spin polarization . Before the bifurcation, there is no spatial fluctuation, i.e. Δ J z =0. The dynamical bifurcation occurs at g 22 * where ΔJ z reaches a small nonzero value δ J , which is chosen as 0.05 in our calculations. Actually, based upon our calculations, similar conclusions can be obtained for other small δ J between 0.05 and 0.2.
In figure 3, we show the dependence of the bifurcation delay on the quench time. Clearly, there is a significant delay between the growth of local spin fluctuation and the static phase transition, see the inset of figure 3. For a smaller quench time, the system has a larger bifurcation delay b d . In our numerical results, it is clearly illustrate that the average bifurcation delay b d follows a power-law scaling with respect to the quench time t µwith the scaling exponent d 2 =0.4937. According to the KZ mechanism, we have  in figure 5. One can find that the average domain size increases with the quench time. In each run, we count the number of domains N d by identifying the number of zero crossings of J z at the end of evolution, when the domain shape becomes stable. The number of domains for different τ Q are shown in figure 4. We observe that the average domain number follows a power-law scaling N d Q d 1 t µwith the scaling exponent d 1 =0.2510. According to the KZ mechanism, we have Combining equations (22) and (23), we obtain the static correlation length critical exponent v=0.5084 and the dynamical critical exponent z=2.0172. These exponents well consist with the ones (v=1/2, z=2) obtained from the Landau critical velocity and the correlation length. These critical exponents are in the same universal class as in the previous work [11,13,16,17,34]. However, in a binary BEC with Rabi coupling [5][6][7][8] across the phase separation driven by quenching the Rabi coupling strength, the exponents are given as (v=1/2, z=1). This means that, although the phase separation can be induced by quenching either the Rabi coupling or the atom-atom interaction, they belong to different universal classes.

Summary and discussion
In summary, we have investigated the non-equilibrium dynamics across a phase separation transition in a mixture of two-component BECs. Through analyzing the Bogoliubov excitation spectrum, we find that the  Landau critical velocity vanishes at the critical point, which results in the spontaneous superfluidity breakdown and creation of elementary excitations, and we extract the critical exponents from the Landau critical velocity and the correlation length near the phase transition. When the system is driven across the critical point with a finite quench rate, the domains are spontaneously created due to the vanishing of the spontaneous superfluidity breakdown. On the other hand, by numerically simulating the GPEs, we find that domains appear after the system crossing the critical point g c 22 , and we count the domains number at the end of evolution. Through quenching the system with various τ Q , we find two universal scalings of the average bifurcation delay and average domain number versus the quench rate. We then numerically extract two critical exponents (v=1/2, z=2) from two universal scalings.
Although all this work and the previous ones [5][6][7][8] study miscible-immiscible transitions, there are significant differences between their KZ scalings. Firstly, the miscible-immiscible transitions are induced by different sources. In the previous works [5][6][7][8], the miscible-immiscible transitions are controlled by the Rabi coupling. On the other hand, in our system, the miscible-immiscible transition is controlled by the interaction coefficient. Secondly, the pictures for understanding the KZ mechanism are different. In the previous works [5][6][7][8], the KZ mechanism is understood via the gapless excitations. However, in our manuscript, the KZ mechanism is understood via the divergence of correlation length derived from the Landau critical velocity. So far, we analytically extract the dynamical critical exponent from the Landau critical velocity and determine the static correlation length critical exponent from the correlation length derived from the Landau critical velocity. Thirdly, and most importantly, the KZ scaling exponents are different. In our manuscript, the scaling exponents are given as (v = 1/2, z = 2). However, in the previous works [5][6][7][8], the scaling exponents are given as (v = 1/2, z = 1). This means that the miscible-immiscible transition driven by quenching the atom-atom interaction and the miscible-immiscible transition driven by quenching the Rabi coupling belong to different universal classes. In the previous works [5][6][7][8], the KZ mechanism originates from the gradually vanishing of the Higgs modes when the system approaches to the critical point, that is, the gapless Higgs modes at the critical point. However, in our work, the long-wavelength excitations are always gapless, and the KZ mechanism originates from the spontaneous superfluidity breakdown induced by the dynamical instability.
Based upon current available techniques for atomic BECs, it is possible to probe the above universal scalings. The mixture of two-component BECs can be prepared with different atomic species [38,39], isotopes [40] or spin states [41][42][43]. The phase separation has been observed in several experiments [40,43] and the high tunability of the interaction strength makes the quenching across the phase separation transition possible. Furthermore, the bifurcation delay and the size of domains can be extracted by measuring the local spin polarization via the time-of-flight. Considering an ensemble of 87 Rb atoms in a ring trap with the total number of atoms N=2×10 6 and the transverse trapping frequency 2 1 w p =´kHz, L=96 in our simulation corresponds to 32.6 um, and g 11 =0.5 corresponds to the scattering length a 11 =13.54 nm.

Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grants No. 11874434, No. 11574405). Figure 5. The domain formation in the real-time dynamics for different quench times τ Q . When the system is quenched from miscible to immiscible phases at a finite τ Q , due to the vanishing of Landau critical velocity at the critical point, elementary excitations emerge and domains are spontaneously created. The parameters are chosen as N=2×10 6 , L=96, g 11 =g 12 =0.5, and τ Q ={10 4 , 10 5 , 10 6 }.