Immiscible and miscible states in binary condensates in the ring geometry

We report detailed investigation of the existence and stability of mixed and demixed modes in binary atomic Bose-Einstein condensates with repulsive interactions in a ring-trap geometry. The stability of such states is examined through eigenvalue spectra for small perturbations, produced by the Bogoliubov-de Gennes equations, and directly verified by simulations based on the coupled Gross-Pitaevskii equations, varying inter- and intra-species scattering lengths so as to probe the entire range of miscibility-immiscibility transitions. In the limit of the one-dimensional (1D) ring, i.e., a very narrow one, stability of mixed states is studied analytically, including hidden-vorticity (HV) modes, i.e., those with opposite vorticities of the two components and zero total angular momentum. The consideration of demixed 1D states reveals, in addition to stable composite single-peak structures, double- and triple-peak ones, above a certain particle-number threshold. In the 2D annular geometry, stable demixed states exist both in radial and azimuthal configurations. We find that stable radially-demixed states can carry arbitrary vorticity and, counter-intuitively, the increase of the vorticity enhances stability of such states, while unstable ones evolve into randomly oscillating angular demixed modes. The consideration of HV states in the 2D geometry expands the stability range of radially-demixed states.

and BEC [52,[77][78][79][80][81][82][83]. Most essential are analytical results for stability of these states. Sec. III presents the key results, showing various types of mixed and demixed states in 1D, characterized by different numbers of peaks in them, and both mixed and demixed 2D states. The latter ones include both radially-demixed modes, with different vorticities, and their azimuthally-demixed counterparts. Such states are obtained by means of the imaginarytime-propagation method, applied to the coupled Gross-Pitaevskii equations (GPEs). We also address effects of the strength of the repulsive intra-component interaction, annular width, and embedded vorticity on the existence and stability of different states. A noteworthy finding is that the stable radially-demixed states can exist with arbitrary vorticity. Our findings are summarized in Sec. IV.

II. THE MEAN-FIELD MODELS
At low temperatures, a binary condensate mixture is well described by the mean-field theory for the set of wave functions φ and ψ of the two components. Here we address the system (e.g., a heteronuclear one) which does not admit linear interconversion (Rabi and/or spin-orbit coupling) between the components. The wave functions obey the system of GPEs with nonlinear terms accounting for self-(intra-species) and cross-(inter-species) interactions. In the normalized form, the GPE system is written as iψ t + 1 2m 2 ∇ 2 ψ − g 2 |ψ| 2 + g 12 |φ| 2 ψ = 0, where m 1,2 are scaled atomic masses, g 1,2 are coefficients of self-interaction in species φ and ψ, and g 12 > 0 is the cross-interaction coefficient. In this work, the analysis is restricted to repulsive interactions, with g 1,2,12 > 0. Then, condition g 12 = √ g 1 g 2 separates the mixing ( √ g 1 g 2 > g 12 ) and phase-separation (demixing, √ g 1 g 2 < g 12 ) regimes in free space [84]. This criterion is modified by the presence of a trapping potential, which tends to enhance the miscibility [42,50,85].
Equations (1) are supplemented by b.c. set at rigid edges, r = r outer and r = r inner of the annular area filled by the condensate (r is the radial coordinate): φ (r = r outer,inner ) = ψ (r = r outer,inner ) = 0. ( By means of scaling, we fix and define the annulus' width, w ≡ r outer − 1. These b.c. imply that the annular area is confined by rigid circular potential walls, as in a recent experiment [86] (performed for a gas of fermions). The total norm of the 2D system is where the integration is performed over the annular region, or over the circumference, in the 1D limit, which corresponds to very tight confinement in the radial direction [see Eq. (33)] below. The energy (Hamiltonian) of the coupled system is which is accordingly reduced in the 1D limit.
A. The one-dimensional setting To define the 1D limit, we assume that the single coordinate, x, running along the ring of radius r = 1 [which is fixed by scaling in agreement with Eq. (3)], takes values 0 ≤ x ≤ 2π. Then, the substitution of solutions in the Madelung form, leads to the system of four real equations for the amplitudes and phases:

The analytical approach in the 1D case
Choosing the constant amplitudes of the two states as a 0 and b 0 respectively, we obtain CW (continuous-wave) solutions of the HV type of Eqs. (8)- (11), Here integer s determines the opposite vorticities in the two components, without introducing net phase circulation.
To address the important issue of the stability of the HV-CW state, or the zero-vorticity one in the case of s = 0, perturbed solutions to Eqs. (8)- (11) are looked for as where σ is the instability growth rate (which may be complex), p is a real wavenumber of the perturbations, while δa, δb and δχ, δη are their infinitely small amplitudes. The substitution of these expressions in Eqs. (8) - (11) and linearization [i.e., the derivation of the respective Bogoliubov -de Gennes (BdG) equations] yields the following dispersion equation for σ(p): Next, we consider two separate cases, depending on the value of s.
(recall we here fix r = 1 by means of scaling). The onset of the transition to the immiscibility (i.e., instability against the phase separation) is signalled by condition σ (p = ±1/r) = 0. It follows from Eq. (16) that this instability takes place at Note that even in the case of g 1 = g 2 = 0 (no self-repulsion), Eq. (18) yields a finite threshold for the onset of the phase-separation instability: This result explicitly demonstrates that periodic b.c. provide for partial stabilization of the mixed state, in comparison with the infinite free space, cf. Ref. [84].
Alternatively, the free term in Eq. Further, Eq. (21) can be cast in the rescaled form, with which implies measuring σ 2 and p 2 in their natural units, and demonstrates that the equation depends on two parameters only, Typical examples of the Σ(P ) dependence for unstable and stable HV states, produced by Eq. (22), are presented in Figs.1(a) and (b), respectively. Note that the underlying quantization condition (17) implies that P takes, in fact, only discrete values: The stability condition means that, for all discrete values of P , given by Eq. (25), Eq. (22) must produce negative real solutions for Σ. Full consideration of the stability conditions following from Eq. (22) is too cumbersome for the analytical investigation. Nevertheless, for the infinite system [r → ∞, i.e., considering P as a continuous variable, Σ. Panel (a) shows that Σ vanishes at P = 0 and P = 6.12, as predicted by Eq. (28).
rather than the discrete one, defined by Eq. (25)], it is easy to obtain the stability condition in the limit of P → 0, for which Eq. (22) amounts to It is easy to see that Eq. (26) produces stable solutions, i.e., real Σ < 0 [see Eqs. (14) and (23)], under condition |g − γ| ≥ 1, i.e., in either of the two cases: According to Eq. (24), conditions (27) hold in the case of a relatively high nonlinearity (large g, or the atom number), or large hidden vorticity, s 2 , which appears in Eq. (24), see further details below. Further, it is possible to find values of P at which Σ vanishes: substituting Σ = 0 in Eq. (22), one obtains P = 0 and P = 4 (γ − g ± 1) .
If Eq. (28) yields P ≤ 0, i.e., g ≥ 1 + γ, cf. Eq. (27), this implies that the HV states are completely stable both for the infinite system and the ring (since, by definition, P may only be positive). Note the modulational stability of uniform HV states with periodic boundary conditions was studied in Ref. [78] for the case of the attractive nonlinearity (on the contrary to the repulsive nonlinearity considered here), for which it was found that the HV-CW states can never be stable.

B. The two-dimensional setting
Stationary solutions to Eq. (1) are looked for in the general form: where (r, θ) are the polar coordinates, µ 1,2 chemical potentials of the two components, S 1,2 = 0, 1, 2, ... their vorticities [87], and real wave functions Φ and Ψ obey the radial equations: To address its stability, we replace the stationary solutions with perturbed ones: where l is an integer azimuthal index of the perturbation with components u 1,2 , v 1,2 , and σ is the instability growth rate. Linearization around the stationary solutions leads to the BdG equations for the two-component condensate: where the prime stands for d/dr. Instabilities are predicted when numerical solution of Eq. (32) produces eigenvalues with Re(σ) = 0. In the 1D version of Eq. (32), d 2 /dr 2 is replaced by d 2 /dx 2 , and terms ∼ 1/r and 1/r 2 are absent. Previously, BdG equations were addressed in the annular geometry defined not by the rigid boundaries, as per Eq. (2), but by weak confinement constructed as the sum of Gaussian and harmonic oscillator potentials [67]. BdG equations for two-component condensates were also studied in other settings, including free space [88,89], 1D configurations [62], and a full analytical solution [64].

III. RESULTS AND DISCUSSION
A. The one-dimensional regime Stationary solutions to Eqs.(1) were produced numerically by means of the imaginary-time-evolution method, using different inputs. Then, stability of these solutions was identified through the calculation of their eigenvalue spectra, using the 1D version of Eq. (32), and further verified by direct numerical simulations of the perturbed evolution. The system conserves the total norm, i.e., scaled number of atoms. N total = N φ + N ψ ≡ 2π 0 |φ| 2 + |ψ| 2 dx. Below, we report numerical results obtained for the basic symmetric states, with m 1 = m 2 = 1 (fixed by scaling), g 1 = g 2 ≡ g, g 12 = 1 [also fixed by scaling, cf. Eq. (20)] and equal 1D norms in the two components, In the miscible phase, the two components of the condensates overlap with each other, whereas they spatially separate in the immiscible phase. A measure to characterize these phases is the overlap integral Λ : In this work, we identified demixed and mixed states as those with Λ = 1 and Λ = 1, respectively. As expected [84,90], demixed states exist only when the cross-repulsion is stronger than the self-repulsion, i.e., g g 12 = 1. They are characterized by local density peaks in each component, located so that a peak in one component coincides with a density minimum in the other, see insets to Figs. 2(b)-(d). Overlap integral (34) for 1D single-peak demixed modes is displayed in Fig. 2(a), as a function of self-repulsive coefficient g, for a fixed norm, N = 10. In this figure, the demixed single-peak mode terminates at g cr = 0.845, only the uniformly mixed state existing at g > g cr . This numerically identified critical value exactly coincides with the analytical prediction given by Eq. (18). Further, the existence area for demixed single-peak and mixed modes is presented in Fig. 2(b). The boundary between them, analytically predicted by Eq. (18), also exactly coincides with the numerically found counterpart, shown by the black curve in Fig. 2(b). The single-peak demixed modes are completely stable in their existence domain, which is consistent with earlier findings [68].
Stability and existence areas of demixed double-and triple-peak modes are displayed in parameter plane (g, N ) in Figs. 2(c) and (d), respectively, and typical examples of such modes are displayed in their respective insets. An essential finding is that, unlike the single-peak modes which are always stable, states with two and three peaks feature instability areas in Figs. 2(c) and (d), being stable only for a sufficiently large norm.
The numerical analysis reveals two instability scenarios for the double-peak mode. If it is taken in the area far from the stability boundary in Fig. 2(c), the real parts of the corresponding eigenvalues σ are relatively large [see Eq. (31)], and the mode spontaneously transforms into an oscillating single-peak state, see Fig. 3(a1-a3). If the unstable mode is selected close to the instability boundary, with smaller real parts of the eigenvalues, it oscillates around itself, rather than transforming into a single-peak state. Similar to the double-peak states, unstable triple-peak ones transform into oscillating single-peaks modes far from the corresponding stability boundary, and persistently oscillate around themselves, if taken close to boundary, see Fig. 3(b1-b3).
We also simulated collision between single-peak demixed components, set in motion by applying opposite kicks to them: Figure 4 shows that, for the kick small enough (k = 0.5), the peaks in the two components periodically bounce back from each other, which is accompanied by some randomization of the patterns. On the other hand, under the action of a strong kick (e.g., k = 5), the moving components pass through each other for about five times, but eventually suffer randomization too, as shown in Fig. 4(b1,b2). Under the action of a still stronger kick, k = 10, the components kept passing through each as long as the simulations were run, see Fig. 4(c1,c2). It is also relevant to simulate evolution of unstable mixed (uniform) states, which is displayed in Fig.5. The instability triggers periodic transformations between the mixed state and a single-peak demixed one, with the period ≈ 10 in this case.

B. Two-dimensional regime
Focusing on the phase-separation scenarios, we identify two different types of 2D demixed modes, namely, those which may be defined as demixed in the radial direction (cf. Ref. [66]), and azimuthally demixed ones, cf. Refs. [62,66,67]. In previous works, similar scenarios of the phase separation were also reported for other binary systems, which include rotation [62,68] and spin-orbit coupling [70].
First, we address 2D zero-vorticity states, including mixed and radially-demixed ones, which may be both stable and unstable (at larger and smaller values of the norm, respectively), as shown in Fig. 6. A typical example of the evolution of unstable 2D radially-demixed states with vorticities S 1,2 = 0 is shown in Fig. 7. It is observed that the unstable state spontaneously evolves into an azimuthally-demixed one.
A noteworthy finding is that the system supports stable 2D radially-demixed states with arbitrarily high vorticities S 1 = S 2 ≡ S. We first analyze 2D demixed states with S = 0 and S = 5 in the parameter space of (g, N, Λ), see  8. It is seen that the solutions are stable (similar to what was found above for other configurations) above a threshold value of the norm, N > N th . We stress that the stability threshold is much lower for S = 5 than for S = 0 [note different scales of vertical axes in Figs. 8(a) and (b)].
To further explore how the vorticity affects the stability of the 2D demixed states, we define the atomic density, [recall that the inner radius of the annulus is scaled to be 1, as per Eq. (3)], and display the stability-threshold value of n as a function of S in Fig. 9(a). A salient feature is the steep drop of n th while S increases from 1 to 2, which is followed by gradual decrease of the threshold with further increase of N . Thus, the vorticity helps to strongly stabilize the axially symmetric states in the annular domain. It is also relevant to investigate an effect of the annulus' width w, defined as per Eq. (4), on the stability. For the zero-vorticity radially-demixed states, the respective stability diagram in parameter panel (n, w) is presented in Fig.  9(b)). It is seen that the stability area strongly broadens with the increase of w, i.e., as it might be expected, stable radially-demixed modes prefer broad annular domains.

Azimuthally-demixed modes (with S = 0)
The 2D setting supports, as well, stable modes which are phase-separated in the azimuthal direction (with zero vorticity) [68,70], an example of such modes can be seen in Fig. 10. These modes are related to their 1D counterparts displayed above in the insets of Figs. 2(b,c,d), and unstable 2D radially-demixed modes transform into them (in an excited oscillating state), see Fig. 7.
To illustrate the evolution of those 2D azimuthally-demixed states which are unstable, we display the evolution of an unstable double-peak state with a small total norm, N = 20 in Fig.11 (similar to other states considered here, they tend to be unstable for relatively small values of N ). It first evolves into a pattern with unequal heights of two peaks, and then restored the original configuration with equal peaks. After several cycles of such shape oscillations, it finally settles into an oscillating single-peak state. The same happens with unstable triple-peak 2D states. This kind of dynamics resembles what was observed above for unstable double-and triple-peak states in the 1D geometry, see Fig. 3.
On the other hand, we have not found any azimuthally-demixed states with nonzero vorticity. Finally, it makes sense to address demixed modes in the full circle, with r inner = 0, instead of the annulus, cf. Eq. (3). It has been found that radially-demixed modes may (quite naturally) exist in the latter case, while no azimuthally-demixed states were found. C. Hidden-vorticity states

The 1D setting
A typical example of the HV state, predicted by analytical solution (12), is presented in Fig. 12. This particular HV state is an unstable one, as illustrated in Fig. 1(a) by the dependence of its instability growth rate on the perturbation wavenumber, which is predicted by Eq. (22); for comparison, Fig. 1 (b) exhibits the same analytical result for a stable HV state. Simulations demonstrate that the evolution transforms unstable 1D HVs into stable single-peak demixed states, with some intrinsic oscillations (not shown here in detail).

The 2D setting
We have also numerically produced 2D radially-demixed HV states, example of which, with S 1,2 = ±1 and ±5, are displayed in Figs. 13(a) and (b), respectively. The same setting may also support 2D mixed HV states, which we do not consider here in detail, as the demixed states seem more interesting. Results for the stability of the 2D radiallydemixed HV modes with the same values of S 1,2 are summarized in Figs. 13(c) and (d). An obviously interesting  conclusion following from the latter plots is that the increase of the hidden vorticity, |S 1,2 |, leads to stabilization of the the HV states [note that difference in the scales of vertical axes in panels (c) and (d)].
Finally, comparing the total energy of different 2D mixed and demixed states [see Eq. (6)], which share equal values of the total norm and angular momentum, we have concluded that the single-peak azimuthally-demixed states realize the lowest energy, i.e., the system's ground state, while the totally mixed configuration has the highest energy.

D. Physical Estimates
To translate the scaled units into the physical ones, we consider the binary condensate of 87 Rb atoms in two different spin states, such as ones with F = 1, m F = 1 and F = 1, m F = 0, and use the same parameters as experiments performed with the two-components condensate in a ring [71], with the radius 12 µm, and the scattering length a s ∼ 10 nm [91]. We conclude that the stable effectively 1D modes predicted by the present analysis may have the actual transverse thickness ∼ 3 µm, containing up to ∼ 10 4 atoms, while the stable 2D modes, predicted for the same outer radius, 12 µm, and the inner one 4 µm, contain 10 4 ∼ 10 5 atoms.

IV. CONCLUSION
We have studied the stability and phase diagram of the two-component BEC loaded in the 2D annular potential box, as well as its 1D limit form corresponding to a ring. The system was analyzed in the framework of the mean-field approximation, based on coupled Gross-Pitaevskii equations with repulsive intra-species and inter-species interactions.
In the 1D setting, the demixed (phase-separated) states are identified as single-, double-and triple-peak modes, with density peaks in one component coinciding with density minima in the other one. The 1D single-peak demixed states are all stable, while the double-and triple-peak ones are stable only above critical values of the total norm, N . The unstable double-and triple-peak modes oscillate around themselves when they are located close to the instability boundary, or spontaneously transform into stable single-peak states deeper in the unstable domain of the parameter space. Collisions between two components of stable demixed single-peak states were studied too, by applying opposite kicks to the components. The simulations demonstrate that the weakly kicked components repeatedly bounce from each other, suffering gradual chaotization, while fast ones pass through each other. If the kicks are moderately strong, the components originally pass through each other, and then evolve into the bouncing regime. The evolution of unstable 1D mixed (spatially uniform) modes shows periodic transitions between the mixed state and single-peak demixed ones.
In the 2D setting, we have found both radially-and azimuthally-demixed states, with unstable radially-demixed ones found to evolve into their azimuthally-demixed counterparts. An essential finding is that the system supports radially-demixed modes with arbitrarily large overall vorticity S, which are stable above the threshold value of the norm, N th . The increase of S leads to stabilization of the modes (decrease of N th ), with a dramatic drop, following the transition from S = 1 to S = 2, in Fig. 9(a). The stability area gradually broadens with the increasing of the annulus' width, w, in Fig. 9(b). Similar to the 1D demixed states, 2D azimuthally-demixed ones are also identified as single-, double-and triple-peak modes. Unstable 2D double-and triple-peak azimuthally-demixed states (those with relatively small norms) evolve into oscillating single-peak modes. In the solid circle, taken instead of the annulus, only radially-demixed modes are found.
Lastly, both 1D and 2D HV (hidden-vorticity) states, with opposite vorticities in the two components, have been addressed too. The stability region for 1D HV modes was found analytically, and fully confirmed by the numerical analysis. Unstable 1D HV modes with components vorticities S 1,2 = ±1 showed evolve into oscillating single-peak demixed modes. The stability domain for 2D radially-demixed HV modes expands with the increase of the hidden vorticity, |S 1,2 |.