Immiscible and miscible states in binary condensates in the ring geometry

TL;DR

This paper investigates the existence and stability of mixed and demixed binary Bose-Einstein condensate states in ring geometries, analyzing their stability through eigenvalue spectra and simulations across 1D and 2D configurations.

Contribution

It provides a comprehensive analysis of stability for miscible and immiscible states in ring-shaped condensates, including analytical and numerical results for various vorticity configurations.

Findings

Stable radially-demixed states can carry arbitrary vorticity.

Increasing vorticity can enhance the stability of certain states.

Demixed states include single, double, and triple-peak structures.

Abstract

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.

Equations

$$i{\varphi }_t+\frac{1}{2m_1}{\nabla }^2\varphi -\left(g_1{\left|\varphi \right|}^2+g_{12}|\psi {|}^2\right)\varphi$$

$$i{\psi }_t+\frac{1}{2m_2}{\nabla }^2\psi -\left(g_2{\left|\psi \right|}^2+g_{12}|\varphi {|}^2\right)\psi$$

$$\varphi \left(r=r_{\mathrm{outer},\mathrm{inner}}\right)=\psi \left(r=r_{\mathrm{outer},\mathrm{inner}}\right)=0.$$

$$r_{\mathrm{inner}}\equiv 1,$$

$$w\equiv r_{\mathrm{outer}}-1.$$

$$N=\int \int \left(|\varphi {|}^2+|\psi {|}^2\right)dxdy\equiv N_{\varphi }+N_{\psi },$$

$$E=\int \int \left[\frac{1}{2m_1}|\nabla \varphi {|}^2+\frac{1}{2m_2}|\nabla \psi {|}^2+\frac{1}{2}(g_1|\varphi {|}^4+g_2|\psi {|}^4)+g_{12}|\varphi {|}^2|\psi {|}^2\right]dxdy,$$

$$\varphi \left(x,t\right)=a\left(x,t\right)\exp \left(i\chi \left(x,t\right)\right),\psi \left(x,t\right)=b\left(x,t\right)\exp \left(i\eta \left(x,t\right)\right),$$

$$a_t+\frac{1}{2m_1}a{\chi }_{xx}+\frac{1}{m_1}{a}_x{\chi }_x=0,$$

$$b_t+\frac{1}{2m_2}b{\eta }_{xx}+\frac{1}{m_2}{b}_x{\eta }_x=0,$$

$$-a{\chi }_t+\frac{1}{2m_1}{a}_{xx}-\frac{1}{2m_1}a{\chi }_x^2-g_1{a}^3-g_{12}{b}^{2}a=0,$$

$$-b{\eta }_t+\frac{1}{2m_2}{b}_{xx}-\frac{1}{2m_2}b{\eta }_x^2-g_2{b}^3-g_{12}{a}^{2}b=0.$$

$$\chi =-{\mu }_{1}t+sx,\eta =-{\mu }_{2}t-sx,$$

$${\mu }_1=g_1{a}_0^2+g_{12}{b}_0^2+\left(s^2/2m_1\right),{\mu }_2=g_2{b}_0^2+g_{12}{a}_0^2+\left(s^2/2m_2\right).$$

$$a\left(x,t\right)$$

$$b\left(x,t\right)$$

$$\chi \left(x,t\right)$$

$$\eta \left(x,t\right)$$

$$\left|\begin{array}[]{cccc}\sigma+i\frac{s}{m_{1}}p&0&-\frac{p^{2}}{2m_{1}}&0\\ 0&\sigma-i\frac{s}{m_{1}}p&0&-\frac{p^{2}}{2m_{2}}\\ -\frac{p^{2}}{2m_{1}}-2g_{1}a_{0}^{2}&-2g_{12}a_{0}b_{0}&-\sigma-i\frac{s}{m_{1}}p&0\\ -2g_{12}a_{0}b_{0}&-\frac{p^{2}}{2m_{2}}-2g_{2}b_{0}^{2}&0&-\sigma+i\frac{s}{m_{2}}p\end{array}\right|=0.$$

$$p^8+4\left(g_1{a}_0^2{m}_1+g_2{b}_0^2{m}_2\right)p^6+4\left[{\sigma }^2{m}_1^2+{\sigma }^2{m}_2^2-4a_0^2{b}_0^2{m}_1{m}_2\left(g_{12}^2-g_1{g}_2\right)\right]p^4$$

$$+16\left(g_1{a}_0^2{m}_2+g_2{b}_0^2{m}_1\right)m_1{m}_2{\sigma }^2{p}^2+16{\sigma }^4{m}_1^2{m}_2^2=0.$$

$$p=n/r,n=0,\pm 1,\pm 2,...$$

$$g_{12}^2-g_1{g}_2>{\left(g_{12}^2-g_1{g}_2\right)}_{\mathrm{cr}}\equiv \frac{r^{-2}\left[4\left(m_1{g}_1{a}_0^2+m_2{g}_2{b}_0^2\right)+r^{-2}\right]}{16m_1{m}_2{a}_0^2{b}_0^2}.$$

$${\left(g_{12}^2-g_1{g}_2\right)}_{\mathrm{cr}}{|}_{g_1=g_2=0}=\frac{r^{-4}}{16m_1{m}_2{a}_0^2{b}_0^2}.$$

$$m_1=m_2\equiv m,g_1=g_2\equiv g,g_{12}\equiv 1,a_0=b_0,$$

$$16m^4{\sigma }^4+\left(32s^2+32gm{a}_0^2+8p^2\right)m^2{p}^2{\sigma }^2+16\left[{\left(gm{a}_0^2-s^2\right)}^2-m^2{a}_0^4\right]p^4+8\left(gm{a}_0^2-s^2\right)p^6+p^8=0.$$

$$16\left[{\left(gm{a}_0^2-s^2\right)}^2-m^2{a}_0^4\right]p^4+8\left(gm{a}_0^2-s^2\right)p^6+p^8\equiv {\left[4\left(gm{a}_0^2-s^2\right)+p^4\right]}^2-16m^2{a}_0^4{p}^4.$$

$$16{\Sigma }^2+8\left(4\gamma +4g+P\right)P\Sigma +16\left[{\left(g-\gamma \right)}^2-1\right]P^2+8\left(g-\gamma \right)P^3+P^4=0.$$

$$\Sigma \equiv {\sigma }^2/a_0^4,P\equiv p^2/\left(m{a}_0^2\right),$$

$$\gamma \equiv s^2/\left(m{a}_0^2\right)\mathrm{and}g.$$

$$P_n=\frac{n^2}{m{a}_0^2{r}^2},b=0,1,2,....$$

$$16{\Sigma }^2+8\left(4\gamma +4g\right)P\Sigma +16\left[{\left(g-\gamma \right)}^2-1\right]P^2=0.$$

$$g\ge 1+\gamma ,\mathrm{or}\gamma \ge 1+g.$$

$$P=0\mathrm{and}P=4\left(\gamma -g\pm 1\right).$$

$$\varphi \left(r,\theta ,t\right)={\Phi }_S\left(r\right)\exp \left(-i{\mu }_{1}t+i{S}_1\theta \right),\psi \left(r,\theta ,t\right)={\Psi }_S\left(r\right)\exp \left(-i{\mu }_{2}t+i{S}_2\theta \right),$$

Full text

I Introduction

Superfluid mixtures are currently routinely probed in experiments with ultracold atomic gases. In addition to Bose-Bose mixtures of different isotopes and atomic species myatt_burt_1997 ; hall_matthews_1998 ; maddaloni_modugno_2000 ; mertes_merrill_2007 ; papp_pino_2008 ; sugawa_yamazaki_2011 ; modugno_modugno_2002 ; thalhammer_barontini_2008 ; mccarron_cho_2011 ; lercher_takekoshi_2011 ; pasquiou_bayerle_2013 ; wacker_jorgensen_2015 ; Petrov_2015 ; wang_li_2016 ; Cabrera_2018 ; Cabrerera2 ; Semeghini_2018 ; Semeghini2 , experimentalists have in the past few years created condensates with a spin degree of freedom stamper-kurn_ueda_2013 , also implementing spin-orbit coupling which gives rise to exciting new states Galitski_2013 ; Zhou_2013 ; Goldman_2014 ; Zhai_2015 ; Malomed_2018 ; moreover, recent achievements have led to the generation of doubly-superfluid Bose-Fermi mixtures ferrier-barbut_delehaye_2014 , in which both components are condensed, a state so far inaccessible in other settings (such as superfluid helium). Although the stability and phase diagrams of such systems have been extensively studied in the course of more than 20 years  esry_greene_97 ; pu_bigelow_1998 ; ao_chui_98 ; timmermans_98 ; ohberg_1999 ; trippenbach_goral_2000 ; delannoy_murdoch_2001 ; riboli_modugno_02 ; jezek_capuzzi_02 ; svidzinsky_chui_2003 ; kasamatsu_tsubota_2004 ; ronen_bohn_2008 ; takeuchi_ishino_2010 ; suzuki_takeuchi_2010 ; mason_aftalion_2011 ; aftalion_mason_2012 ; wen_liu_12 ; pattinson_billam_2013 ; hofmann_natu_2014 ; pattinson_parker_2014 ; edmonds_lee_2015a ; edmonds_lee_2015b ; lee_proukakis_2016 ; liu_pattinson_2016 ; lee_jorgensen_2016 ; lee_jorgensen_2018 , even simple hetero-species binary mixtures still reveal unexpected features, such as the role of the trap sag, atom number and kinetic energy contribution to the extent of miscibility in trapped configurations lee_jorgensen_2016 ; wen_liu_12 ; Richaud_2019 , and nontrivial effects of the expansion on the mixtures’ dynamics lee_jorgensen_2018 ; Yongyao_2018 .

Configurations which keep drawing growing interest in studies of ultracold atomic gases are based on the annular, alias ring-trap, geometry gupta_murch_2005 ; arnold_garvie_2006 ; ryu_andersen_2007 ; ramanathan_wright_2011 ; eckel_lee_2014 ; corman_chomaz_2014 ; kumar_eckel_2017 . These configurations are interesting as they lead to closed geometries with controlled flows, that are also of potential use to the emerging field of atomtronics eckel_lee_2014 ; Amico_2017 . In this context, mixtures of atomic condensates in toroidal traps and the possibility of sustaining stable persistent currents in them have been previously considered in smyrnakis_bargi_2009 ; shimodaira_kishimoto_2010 ; bargi_malet_2010 ; anoshkin_wu_2013 ; yakimenko_isaieva_2013 ; mason_2013 ; abad_sartori_2014 ; white_hennessy_2016 ; Yakimenko_2015 ; white_zhang_2017 , and the corresponding experimental observation beattie_moulder_2013 has helped to clarify some issues, also raising new questions, such as expansion of the stability area for such states. The aim of the present work is to perform a full classification of accessible stable mixture states in such a geometry, including examination of their stability and decay channels of their unstable counterparts, both in the absence and presence of overall rotation. Given the potential significance of multi-component states in ring-shaped traps for applications such as rotational sensors, such a classification is relevant. It can also assist in developing methods for control of such mixtures in the experimental work which is currently going on in many laboratories.

Specifically, in this work, we construct a binary Bose-Einstein-condensate (BEC) system trapped in an annular geometry, through the mean-field analysis in the presence of inter- and intra-species interactions, whose parameters are varied in broad limits. After analyzing the corresponding one-dimensional (1D) problem, we focus on the more experimentally-relevant 2D annular structure. We implement periodic boundary conditions (b.c.) in the azimuthal direction in the 1D case, and zero b.c. at inner and outer boundaries of the 2D annular structure. The latter b.c. set enables one to study how the annular structure affects density patterns in the repulsive bosonic mixtures, as a result of the existence of different demixed and mixed states, and their stability.

This paper is structured as follows. First, Sec. II introduces and analytically considers our basic model for the mixtures in both 1D and 2D geometries, and presents analytical results for spatially uniform 1D mixed solutions, with zero and hidden vorticities (HV), the latter implying opposite topological charges in the two components, which makes it possible to construct stable binary vortex states with zero total angular momentum in nonlinear optics desyatnikov_kivshar_2001 ; desyatnikov_2005 ; bigelow_park_2002 ; Leykam_2013 ; Salgueiro_2016 and BEC mihalache_mazilu_2006 ; brtka_gammal_2010 ; Yakimenko_2012 ; He_2012 ; Linghua_2013 ; Ishino_2013 ; Hoashi_2016 ; Yongyao_2018 . 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 imaginary-time-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 $\phi$ and $\psi$ 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

[TABLE]

where $m_{1,2}$ are scaled atomic masses, $g_{1,2}$ are coefficients of self-interaction in species $\phi$ and $\psi$, 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}=\sqrt{g_{1}g_{2}}$ separates the mixing ($\sqrt{g_{1}g_{2}}>g_{12}$) and phase-separation (demixing,$\sqrt{g_{1}g_{2}}<g_{12}$) regimes in free space Mineev_1974 . This criterion is modified by the presence of a trapping potential, which tends to enhance the miscibility Merhasin_2005 ; lee_jorgensen_2016 ; wen_liu_12 .

Equations (LABEL:phipsi) are supplemented by b.c. set at rigid edges, $r=r_{\mathrm{outer}}$ and $r=r_{\mathrm{inner}}$ of the annular area filled by the condensate ($r$ is the radial coordinate):

[TABLE]

By means of scaling, we fix

[TABLE]

and define the annulus’ width,

[TABLE]

These b.c. imply that the annular area is confined by rigid circular potential walls, as in a recent experiment Hueck_2018 (performed for a gas of fermions).

The total norm of the 2D system is

[TABLE]

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

[TABLE]

which is accordingly reduced in the 1D limit.

II.1 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\leq x\leq 2\pi$. Then, the substitution of solutions in the Madelung form,

[TABLE]

leads to the system of four real equations for the amplitudes and phases:

[TABLE]

II.1.1 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),

[TABLE]

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

[TABLE]

where $\sigma$ is the instability growth rate (which may be complex), $p$ is a real wavenumber of the perturbations, while $\delta a,\delta b$ and $\delta\chi,\delta\eta$ 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 $\sigma(p)$:

[TABLE]

Next, we consider two separate cases, depending on the value of $s$.

II.1.2 Zero-vorticity states, $s=0$

For $s=0$, determinant (15) defining the stability takes the explicit form

[TABLE]

Due to the periodic b.c. set by the ring geometry, $p$ is quantized,

[TABLE]

(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 $\sigma\left(p=\pm 1/r\right)=0$. It follows from Eq. (16) that this instability takes place at

[TABLE]

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:

[TABLE]

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. Mineev_1974 .

II.1.3 Hidden-vorticity (HV) states, $s\geq 1$

As defined above, HV states carry opposite angular momenta in the two components of the mixture, while the total momentum is zero. In an explicit form, the corresponding equation (15), which determines their stability, takes a very cumbersome form. It becomes relatively simple in the case of full symmetry in Eqs. (8)-(11) and (12)-(13), viz.,

[TABLE]

for which we obtain $\mu_{1}=\mu_{2}=-\left(g+1\right)a_{0}^{2}+s^{2}/\left(2m\right)$. Then, Eq. (15) can be explicitly written as

[TABLE]

Alternatively, the free term in Eq. (21) (the part which does not contain $\sigma^{2}$) can be written as

[TABLE]

Further, Eq. (21) can be cast in the rescaled form,

[TABLE]

with

[TABLE]

which implies measuring $\sigma^{2}$ and $p^{2}$ in their natural units, and demonstrates that the equation depends on two parameters only,

[TABLE]

Typical examples of the $\Sigma(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:

[TABLE]

The stability condition means that, for all discrete values of $P$, given by Eq. (25), Eq. (22) must produce negative real solutions for $\Sigma$. Full consideration of the stability conditions following from Eq. (22) is too cumbersome for the analytical investigation. Nevertheless, for the infinite system [$r\rightarrow\infty$, i.e., considering $P$ as a continuous variable, rather than the discrete one, defined by Eq. (25)], it is easy to obtain the stability condition in the limit of $P\rightarrow 0$, for which Eq. (22) amounts to

[TABLE]

It is easy to see that Eq. (26) produces stable solutions, i.e., real $\Sigma<0$ [see Eqs. (LABEL:pert) and (23)], under condition $\left|g-\gamma\right|\geq 1$, i.e., in either of the two cases:

[TABLE]

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 $\Sigma$ vanishes: substituting $\Sigma=0$ in Eq. (22), one obtains

[TABLE]

If Eq. (28) yields $P\leq 0$, i.e., $g\geq 1+\gamma$, 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. brtka_gammal_2010 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.

II.2 The two-dimensional setting

Stationary solutions to Eq. (LABEL:phipsi) are looked for in the general form:

[TABLE]

where $(r,\theta)$ are the polar coordinates, $\mu_{1,2}$ chemical potentials of the two components, $S_{1,2}=0,1,2,...$ their vorticities old , and real wave functions $\Phi$ and $\Psi$ obey the radial equations:

[TABLE]

To address its stability, we replace the stationary solutions with perturbed ones:

[TABLE]

where $l$ is an integer azimuthal index of the perturbation with components $u_{1,2}$, $v_{1,2}$, and $\sigma$ is the instability growth rate.

Linearization around the stationary solutions leads to the BdG equations for the two-component condensate:

[TABLE]

where the prime stands for $d/dr$. Instabilities are predicted when numerical solution of Eq. (32) produces eigenvalues with Re$(\sigma)\neq 0$. In the 1D version of Eq. (32), $d^{2}/dr^{2}$ is replaced by $d^{2}/dx^{2}$, and terms $\sim 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 abad_sartori_2014 . BdG equations for two-component condensates were also studied in other settings, including free space Shungo_2011 ; CKLaw_2001 , 1D configurations shimodaira_kishimoto_2010 , and a full analytical solution anoshkin_wu_2013 .

III Results and discussion

III.1 The one-dimensional regime

Stationary solutions to Eqs.(LABEL:phipsi) 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_{\mathrm{total}}=N_{\phi}+N_{\psi}\equiv\int_{0}^{2\pi}\left(|\phi|^{2}+|\psi|^{2}\right)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}\equiv g$, $g_{12}=1$ [also fixed by scaling, cf. Eq. (20)] and equal 1D norms in the two components,

[TABLE]

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, which we define here in the 2D form, as it will be used below in the analysis of the 2D setting:

[TABLE]

the reduction of the definition to the 1D limit being obvious. In this work, we identify demixed and mixed states as those with $\Lambda\neq 1$ and $\Lambda=1$, respectively.

As expected Mineev_1974 ; Lingua , demixed states exist only when the cross-repulsion is stronger than the self-repulsion, i.e., $g\leqslant 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_{\mathrm{cr}}=0.845$, only the uniformly mixed state existing at $g>g_{\mathrm{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 white_hennessy_2016 .

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 $\sigma$ are relatively large [see Eq. (LABEL:PtbSol)], 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-peak 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:

[TABLE]

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 $\approx 10$ in this case.

III.2 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. mason_2013 ), and azimuthally demixed ones, cf. Refs. mason_2013 ; abad_sartori_2014 ; shimodaira_kishimoto_2010 . In previous works, similar scenarios of the phase separation were also reported for other binary systems, which include rotation shimodaira_kishimoto_2010 ; white_hennessy_2016 and spin-orbit coupling white_zhang_2017 .

III.2.1 Radially-demixed modes

In the consideration of the 2D setting subject to b.c. (2), we focus, as above, on the scaled symmetric system, with $m_{1}=m_{2}\equiv 1$, $g_{1}=g_{2}\equiv 1$, $g_{12}=1$, and $r_{\mathrm{inner}}=1$, cf. Eqs. (20) and (3). Stability of 2D modes was identified by the computation of the eigenvalue spectra in the corresponding Bogoliubov-de Gennes (BdG) equations (32), and further verified by direct simulations.

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}\equiv S$. We first analyze 2D demixed states with $S=0$ and $S=5$ in the parameter space of $(g,N,\Lambda)$, see Fig. 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_{\mathrm{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,

[TABLE]

[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_{\mathrm{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.

III.2.2 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) white_hennessy_2016 ; white_zhang_2017 , 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_{\mathrm{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.

III.3 Hidden-vorticity states

III.3.1 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).

III.3.2 The 2D setting

We have also numerically produced 2D radially-demixed HV states, example of which, with $S_{1,2}=\pm 1$ and $\pm 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 radially-demixed 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, $\left|S_{1,2}\right|$, 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.

III.4 Physical Estimates

To translate the scaled units into the physical ones, we consider the binary condensate of 87Rb 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 beattie_moulder_2013 , with the radius $\lesssim 12$ $\mu$m, and the scattering length $a_{s}\sim 10$ nm Egorov_2013 . We conclude that the stable effectively 1D modes predicted by the present analysis may have the actual transverse thickness $\sim 3$ $\mu$m, containing up to $\sim 10^{4}$ atoms, while the stable 2D modes, predicted for the same outer radius, $\lesssim 12$ $\mu$m, and the inner one $\lesssim 4$ $\mu$m, contain $10^{4}\sim 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_{\mathrm{th}}$. The increase of $S$ leads to stabilization of the modes (decrease of $N_{\mathrm{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}=\pm 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, $\left|S_{1,2}\right|$.

Acknowledgments

This work was supported, in part, by grants No.11874112 and No.11575063 from NNSFC (China), by grant No. 2015616 from the joint program in physics between NSF and Binational (US-Israel) Science Foundation, and by grant No. 1286/17 from the Israel Science Foundation, Z.C. acknowledges an excellence scholarship provided by the Tel Aviv University. B.A.M. appreciates hospitality of the Joint Quantum Centre (JQC) Durham-Newcastle, during his Visiting Professorship at Newcastle University, and NPP funding from EPSRC, UK (Grant No. EP/K03250X/1).