The mixing-demixing phase diagram of ultracold heteronuclear mixtures in a ring trimer

TL;DR

This paper maps the complete phase diagram of mixing and demixing in ultracold heteronuclear mixtures confined in a ring trimer, combining Bose-Hubbard and Gross-Pitaevskii models to guide experimental realizations.

Contribution

It provides a comprehensive analysis of the mixing-demixing transitions in heteronuclear mixtures using both quantum and mean-field models, extending to realistic experimental conditions.

Findings

The phase diagram depends only on two effective parameters.

The mean-field analysis confirms a rich variety of transition scenarios.

Proposes an experimental setup with Na and K bosonic mixtures.

Abstract

We derive the complete mixing-demixing phase-diagram relevant to a bosonic binary mixture confined in a ring trimer and modeled within the Bose-Hubbard picture. The mixing properties of the two quantum fluids, which are shown to be strongly affected by the fragmented character of the confining potential, are evaluated by means of a specific indicator imported from Statistical Thermodynamics and are shown to depend only on two effective parameters incorporating the asymmetry between the heteronuclear species. To closely match realistic experimental conditions, our study is extended also beyond the pointlike approximation of potential wells by describing the systems in terms of two coupled Gross-Pitaevskii equations. The resulting mean-field analysis confirms the rich scenario of mixing-demixing transitions of the mixture and also constitutes an effective springboard towards a viable experimental realization. We additionally propose an experimental realization based on a realistic optical-tweezers system and on mixtures of bosonic $^{23}\mathrm{Na}$ and $^{39}\mathrm{K}$, thanks to the large tunability of their intra- and inter-species scattering lengths.

Equations

$$H=-T_a\sum _{j=1}^3\left(a_{j+1}^{†}{a}_j+a_j^{†}{a}_{j+1}\right)+\frac{U_a}{2}\sum _{j=1}^3{n}_j(n_j-1)-T_b\sum _{j=1}^3\left(b_{j+1}^{†}{b}_j+b_j^{†}{b}_{j+1}\right)+\frac{U_b}{2}\sum _{j=1}^3{m}_j(m_j-1)+W\sum _{j=1}^3{n}_j{m}_j,$$

$$\begin{array}{rl}V=& -2N_a{T}_a\left(\sqrt{x_1{x}_2}+\sqrt{x_2{x}_3}+\sqrt{x_3{x}_1}\right)-2N_b{T}_b\left(\sqrt{y_1{y}_2}+\sqrt{y_2{y}_3}+\sqrt{y_3{y}_1}\right)+\\ & +\frac{U_a{N}_a^2}{2}\left(x_1^2+x_2^2+x_3^2\right)+\frac{U_b{N}_b^2}{2}\left(y_1^2+y_2^2+y_3^2\right)+W{N}_a{N}_b\left(x_1{y}_1+x_2{y}_2+x_3{y}_3\right)\end{array}$$

$$V_{\ast }=\frac{1}{2}\left(x_1^2+x_2^2+x_3^2\right)+\frac{{\beta }^2}{2}\left(y_1^2+y_2^2+y_3^2\right)+\alpha \beta \left(x_1{y}_1+x_2{y}_2+x_3{y}_3\right),$$

$$\alpha =\frac{W}{\sqrt{U_a{U}_b}},\beta =\frac{N_b}{N_a}\sqrt{\frac{U_b}{U_a}}.$$

$$x_1=0,x_2=\frac{{\alpha }^2-\alpha \beta -2}{{\alpha }^2-4},x_3=\frac{\alpha \beta -2}{{\alpha }^2-4},$$

$$y_1=\frac{{\alpha }^2\beta -\alpha -2\beta }{({\alpha }^2-4)\beta },y_2=0,y_3=\frac{\alpha -2\beta }{\left({\alpha }^2-4\right)\beta }.$$

$$x_1=0,x_2=\frac{1}{2},x_3=\frac{1}{2},$$

$$y_1=1,y_2=0,y_3=0.$$

$$x_1=\frac{1}{3}(1-2\alpha \beta ),x_2=\frac{1}{3}(1+\alpha \beta ),x_3=\frac{1}{3}(1+\alpha \beta ),$$

$$y_1=1,y_2=0,y_3=0.$$

$$S_{mix}=-\frac{1}{2}\sum _{j=1}^3\left(x_j\log \frac{x_j}{x_j+y_j}+y_j\log \frac{y_j}{x_j+y_j}\right).$$

$$V_{\ast ,1}=\frac{1}{6}\left(2\alpha \beta +{\beta }^2+1\right),V_{\ast ,2}=\frac{\left({\alpha }^2-2\right){\beta }^2+{\alpha }^2-2\alpha \beta -2}{2\left({\alpha }^2-4\right)},$$

$$V_{\ast ,3}=\frac{1}{4}\left(2{\beta }^2+1\right),V_{\ast ,4}=\frac{1}{6}\left[\beta \left(-2{\alpha }^2\beta +2\alpha +3\beta \right)+1\right].$$

$$\alpha <\sqrt{\left(1+\frac{9}{2}\frac{T_a}{U_a{N}_a}\right)\left(1+\frac{9}{2}\frac{T_b}{U_b{N}_b}\right)}.$$

$$-\frac{{\hslash }^2}{2m_a}{\nabla }^2{\phi }_a+g_a|{\phi }_a{|}^2{\phi }_a+g_{ab}|{\phi }_b{|}^2{\phi }_a+V_{ext,a}{\phi }_a={\mu }_a{\phi }_a$$

$$-\frac{{\hslash }^2}{2m_b}{\nabla }^2{\phi }_b+g_b|{\phi }_b{|}^2{\phi }_b+g_{ab}|{\phi }_a{|}^2{\phi }_b+V_{ext,b}{\phi }_b={\mu }_b{\phi }_b$$

$${\int }_{-\infty }^{+\infty }|{\phi }_a{|}^2{\mathrm{d}}^{3}x=N_a,{\int }_{-\infty }^{+\infty }|{\phi }_b{|}^2{\mathrm{d}}^{3}x=N_b$$

$$g_a=\frac{4\pi {\hslash }^2{a}_a}{m_a}=1.01\times 10^{-50}J\text{and}{g}_b=\frac{4\pi {\hslash }^2{a}_b}{m_b}\in [8.67,9.20]\times 10^{-52}J$$

$$g_{ab}=\frac{2\pi {\hslash }^2{a}_{ab}}{{\mu }_{ab}}\in [1.57,7.44]\times 10^{-51}J$$

$$\mathcal{H}[{\phi }_a,{\phi }_b]=\int {\phi }_a^{\ast }\left[-\frac{{\hslash }^2}{2m_a}{\nabla }^2+V_{ext,a}\right]{\phi }_a{\mathrm{d}}^{3}x+\frac{g_a}{2}\int |{\phi }_a{|}^4{\mathrm{d}}^{3}x+$$

$$+\int {\phi }_b^{\ast }\left[-\frac{{\hslash }^2}{2m_b}{\nabla }^2+V_{ext,b}\right]{\phi }_b{\mathrm{d}}^{3}x+\frac{g_b}{2}\int |{\phi }_b{|}^4{\mathrm{d}}^{3}x+g_{ab}\int |{\phi }_a{|}^2|{\phi }_b{|}^2{\mathrm{d}}^{3}x,$$

$${\psi }_a={\psi }_a(x,y,z)={\psi }_{a,x}(x){\psi }_{a,y,z}(y,z)={\psi }_{a,x}(x)\frac{1}{\sqrt{L_y}}\frac{1}{\sqrt{L_z}}$$

$${\psi }_b={\psi }_b(x,y,z)={\psi }_{b,x}(x){\psi }_{b,y,z}(y,z)={\psi }_{b,x}(x)\frac{1}{\sqrt{L_y}}\frac{1}{\sqrt{L_z}}$$

$${\phi }_a=\sqrt{N_a}{\psi }_a,{\phi }_b=\sqrt{N_b}{\psi }_b$$

$$-\frac{{\partial }^2}{\partial x^2}{\psi }_a+{\gamma }_a|{\psi }_a{|}^2{\psi }_a+{\gamma }_{ab}|{\psi }_b{|}^2{\psi }_a+{\tilde{V}}_{ext,a}{\psi }_a={\tilde{\mu }}_a{\psi }_a$$

$$-\frac{{\partial }^2}{\partial x^2}{\psi }_b+{\gamma }_b|{\psi }_b{|}^2{\psi }_b+{\gamma }_{ba}|{\psi }_a{|}^2{\psi }_b+{\tilde{V}}_{ext,b}{\psi }_b={\tilde{\mu }}_b{\psi }_b$$

$${\gamma }_a=\frac{2m_a{g}_a{N}_a}{L_y{L}_z{\hslash }^2},{\gamma }_{ab}=\frac{2m_a{g}_{ab}{N}_b}{L_y{L}_z{\hslash }^2},{\tilde{\mu }}_a=\frac{2m_a}{{\hslash }^2}{\mu }_a,$$

$${\gamma }_b=\frac{2m_b{g}_b{N}_b}{L_y{L}_z{\hslash }^2},{\gamma }_{ba}=\frac{2m_b{g}_{ab}{N}_a}{L_y{L}_z{\hslash }^2},{\tilde{\mu }}_b=\frac{2m_b}{{\hslash }^2}{\mu }_b.$$

$${\tilde{V}}_{ext,a}=-\frac{2m_a}{{\hslash }^2}{P}_a\sum _{i=1}^3{e}^{-\frac{2(x-x_i{)}^2}{{\sigma }^2}},{\tilde{V}}_{ext,b}=-\frac{2m_b}{{\hslash }^2}{P}_b\sum _{i=1}^3{e}^{-\frac{2(x-x_i{)}^2}{{\sigma }^2}}$$

$$\frac{{\partial }^2{\psi }_a}{\partial {\tau }^2}+\eta \frac{\partial {\psi }_a}{\partial \tau }=\frac{{\partial }^2}{\partial x^2}{\psi }_a-{\gamma }_a|{\psi }_a{|}^2{\psi }_a-{\gamma }_{ab}|{\psi }_b{|}^2{\psi }_a-{\tilde{V}}_{ext,a}{\psi }_a+{\tilde{\mu }}_a{\psi }_a$$

Full text

Introduction

Ultracold mixtures of different atomic species have a long experimental history both for studies in bulk [1] and in optical lattices [2, 3, 4, 5] and allowed the investigation of a plethora of intriguing phenomena such as entropy exchange[6], miscibility in degenerate gases[7], beyond-mean-field effects[8, 9, 10], universality of three-body losses[11, 12, 13] and ultracold chemistry [14].

Fragmentation of bosonic mixtures in arrays of potential wells constitutes an important factor in determining their mixing properties and plays a crucial role, together with the repulsion between the mixture’s components, in triggering spatial phase separation. A large amount of recent theoretical studies have revealed a colorful phenomenology ranging from phase separation mechanism [15, 16], demixing of dipolar mixtures[17], and quantum emulsions [18, 19] to the critical properties of quasiparticles spectrum across the transition[20], thermal suppression[21] and the role of entanglement[22] in phase separation, and cascades of quantum phase transitions of vector bosons[23].

Recently, phase separation has been investigated for a binary mixture in a two-well trap (dimer). This constitutes the simplest geometry whose fragmented character exacerbates the competition between intra- and inter-species repulsions[24] therefore affecting the occurrence of spatial separation. In this system, in fact, the transition from a mixed to a fully demixed ground-state, which transparently emerges within the Continuous Variable Picture (CVP)[25], is confirmed by the analysis of genuinely quantum indicators, such as the energy spectrum, the entanglement between the species, and the quantum tomography of the ground state [26, 22]. The study of phase-separation mechanism has been extended to the ring-trimer geometry in the simple case of twin species [27], i.e. species featuring the same hopping amplitude $T$, the same on-site interaction $U$, and the same number of atoms $N$. Despite its simplicity, this toy model has already shown a non-trivial phenomenology, marked by the presence of an unexpected intermediate phase which is in between a perfectly mixed and a fully demixed phase.

In this work, we aim at ferrying this intriguing but maybe too simplistic model towards an actual experimental realization. In doing so, we take into account the presence of asymmetries $T_{a}\neq T_{b}$, $U_{a}\neq U_{b}$ and $N_{a}\neq N_{b}$ between atomic species a and b, a circumstance which makes the model we are dealing with more realistic. The increased number of variables that come into play therefore makes the developed analysis remarkably more articulated and comprehensive than the original one, but significantly broadens its scope. Such study results in a rather rich two-dimensional phase diagram, which includes four different phases, each one marked by a specific boson redistribution among the wells and featuring a different degree of mixing. In order to corroborate our previsions, we extend our analysis beyond the pointlike-site approximation (reducing the local potential traps to dimensionless sites) which intrinsically characterizes both the Bose-Hubbard model and the CVP. To this end, we present a mean-field (MF) approach based on the solution of two coupled stationary Gross-Pitaevskii equations (GPEs). This technique, which allows one to take into consideration the boson distributions not only among but also within the wells, lends itself to direct interlink with a real-world experiment and therefore represents an effective bridge between the quite abstract (but rather clarifying) CVP model and a real physical system.

The model

In this paper, we consider a bosonic binary mixture confined in a three-well potential with periodic boundary conditions (i.e. a ring trimer). This is well described by the Bose-Hubbard Hamiltonian

[TABLE]

a second-quantized model where operator $a_{i}$ ($a_{i}^{\dagger}$) destroys (creates) a species-a boson in site $i\in\{1,2,3\}$ (with $i=4\equiv 1$ due to the ring geometry). The same holds for operators $b_{i}$ and $b_{i}^{\dagger}$, which destroy and create species-b bosons. The bosonic character of these operators is enforced by standard commutation relations, i.e $[a_{i},a_{\ell}^{\dagger}]=\delta_{i,\ell}=[b_{i},b_{\ell}^{\dagger}]$, $[a_{i},b_{\ell}^{\dagger}]=0$. After introducing number operators $n_{i}=a_{i}^{\dagger}a_{i}$ and $m_{i}=b_{i}^{\dagger}b_{i}$, one can observe that the total numbers of bosons of each species ($N_{a}=\sum_{i}n_{i}$ and $N_{b}=\sum_{i}m_{i}$ respectively) constitute two independent conserved quantities, namely $[N_{a},H]=[N_{b},H]=0$. As regards model parameters, $T_{a}$ and $T_{b}$ are hopping amplitudes, $U_{a}$ and $U_{b}$ represent intra-species repulsive interactions, and $W$ represents the inter-species repulsion.

Semiclassical approach towards demixing.

The CVP, originally introduced to investigate the spatial fragmentation of a condensate in a two-well potential [28], is a semi-classical approximation scheme based on the replacement of the inherently discrete quantum numbers associated to the Fock-state basis with continuous variables. This technique has proved to be particularly effective in capturing the essential critical behaviour of complex many-body systems [29, 30, 31]. It also allows one to derive an effective semi-classical Hamiltonian $H_{eff}$ which well reproduces the low-energy physics of the original quantum model [32, 33, 34] and to turn the search for the ground state thereof into that of the minimum of the potential $V$ provided by $H_{eff}$. In the current case, the CVP associates to Hamiltonian (1) the potential

[TABLE]

which comes with the two constraints $\sum_{i}x_{i}=\sum_{i}y_{i}=1$, enforcing particle number conservation in both condensed species. Notice that variables $x_{i}:=n_{i}/N_{a}$ and $y_{i}:=m_{i}/N_{b}$ represent normalized boson populations and are regarded as continuous in view of the fact that the total numbers of bosons $N_{a}$ and $N_{b}$ are assumed to be large.

Results

We investigate the ground-state properties of a bosonic binary mixture in a ring trimer, showing that they can be conveniently represented by means of a two-dimensional phase diagram which is, in turn, spanned by two specific effective variables representing functions of the original-model parameters. Each phase features a characteristic boson-population distribution, which results in a different degree of mixing and in a different expression of the ground-state energy.

The phase diagram.

The presence of four different phases emerges in a crystal-clear way when one considers the large-populations limit (this can be seen as sort of thermodynamic limit according to the scheme highlighted in [33]), i.e. when ratios $U_{a}N_{a}/T_{a}$ and $U_{b}N_{b}/T_{b}$ are big enough to make hopping terms negligible. In this limit, in fact, the (rescaled version of) effective potential (2), namely $V/(U_{a}N_{a}^{2})$, reduces to

[TABLE]

an expression including only two effective parameters,

[TABLE]

The former corresponds to the ratio between the inter-species and the (geometric average of) the intra-species repulsive interactions while the latter accounts for the degree of asymmetry between the two species. This reduced set of parameters proves to be the most natural one to investigate the occurrence of mixing-demixing transitions in the system and is therefore employed to build the phase diagram.

In this regard, the region $\mathcal{D}=\{(\alpha,\beta)\,:\,0\leq\alpha\leq 3,\,0\leq\beta\leq 1\}$, which is depicted in Fig. 1 already includes all four phases that the system’s ground state can exhibit (notice, in fact, that if $\beta>1$, one can swap species labels $a$ and $b$ and therefore come back to region $\mathcal{D}$). In each phase, the configuration $(\vec{x},\vec{y})$ which minimizes effective potential (3) is marked by a specific dependence on the effective model parameters (4). Analytic expressions $x_{j}(\alpha,\beta)$ and $y_{j}(\alpha,\beta)$ have been derived by means of an exhaustive exploration of the polytope-like domain of function $(\ref{eq:V_*})$ (see [27] for a detailed description of the method), and are presented here, together with the essential features of the associated phase:

In phase 1, being $\alpha<1$, the inter-species repulsion is too small to trigger spatial phase separation and the two species are uniformly distributed among the three wells, i.e. $x_{j}=y_{j}=1/3,\,\forall j=1,2,3$. 2. 2.

Phase 2 occurs for $1<\alpha<2\beta$, that means for intermediate values of $\alpha$ and not too asymmetric species. In this region, one has complete demixing in two wells and mixing in the remaining one. The explicit expression of the minimum-energy configuration as functions of model parameters (4) reads

[TABLE]

[TABLE] 3. 3.

Phase 3, occurring for sufficiently high values of $\alpha$ and not too asymmetric species (more specifically, for $1/(2\alpha)<\beta<\alpha/2$), features a completely demixed scenario, such that species $b$ clots in one well, while species $a$ equally spreads in the other two sites:

[TABLE]

[TABLE] 4. 4.

Phase 4 occurs for sufficiently high values of $\alpha$ and sufficiently asymmetric species (i.e. for $\alpha>1$ and $\beta<1/(2\alpha)$). A good degree of asymmetry can be achieved, for example, if species $b$ constitutes an impurity with respect to the majoritarian species $a$, meaning that $N_{b}\ll N_{a}$. The hallmark of this phase is that species $b$ clots in one site, while species $a$ spreads in all three sites, but in different proportions:

[TABLE]

[TABLE]

An illustrative summary of all presented minimum-energy configurations is provided in Fig. 2. Concerning the ones discussed at points 2, 3 and 4, we remark that, due to the $Z_{3}$ symmetry of the system, they are not unique and that other isoenergetic configurations can be obtained by cyclic permutations of site indexes.

Quantum-mechanically, the ground state is degenerate only in the infinite-population limit because, as soon as they come into play, hopping terms lift the degeneracy and the ground state’s structure gets that of a Schrödinger cat [30, 35, 27, 22]. For example, in phase 3, the ground-state is of the type $|E_{0}\rangle\approx\frac{1}{\sqrt{3}}\left(|0,\frac{1}{2},\frac{1}{2}\rangle_{a}|1,0,0\rangle_{b}+|\frac{1}{2},0,\frac{1}{2}\rangle_{a}|0,1,0\rangle_{b}+|\frac{1}{2},\frac{1}{2},0\rangle_{a}|0,0,1\rangle_{b}\right)$, an expression where symbol “$\approx$" reminds that minor contributions coming from Fock states with different boson distributions and activated by the non-zero hopping amplitudes have been understood.

As clearly illustrated in Fig. 2, for $\beta<1$, the minimum-energy configuration $(\vec{x},\vec{y})$, as a function of model parameters $\alpha$ and $\beta$, is discontinuous at transitions 1-2 and 1-4 while it is continuous at transitions 2-3 and 3-4. Nevertheless, in the case of perfectly symmetric species, $\beta=1$, the system gains an additional symmetry (consisting in the interchangeability of species’ labels in formulas (1) - (3)) and a further discontinuity at transition 2-3 appears. The twin-species limit, widely discussed in [27], therefore features a qualitatively different critical behavior.

Phases and degree of mixing.

An effective indicator to quantify the degree of mixing of two different species in discretized domains is the entropy of mixing ($S_{mix}$). Originally introduced in the context of macromolecular simulations[36], this measure has been recently introduced in the realm of ultracold atoms in order to investigate the link between chaotic dynamical regimes and mixing properties of a bosonic binary mixture in a ring trimer [37]. According to the definition given in [36], the entropy of mixing associated to a certain minimum-energy configuration $(\vec{x},\vec{y})$ reads

[TABLE]

In comparison to the one used in [38, 37], this formula is based on normalized populations $x_{j}$ and $y_{j}$ (rather than $n_{j}$ and $m_{j}$) and is therefore more suited to describe binary mixtures where the species feature a particle-number imbalance (moreover the contribution of each site $j$ to the total entropy of mixing is not fixed, but it is weighted by the fractions of particles present therein). As shown in Fig. 3, $S_{mix}$ is zero in phase 3 (perfect demixing) while achieves the maximum possible value, $\log 2\approx 0.69$ in phase 1 (perfect mixing). More generally, $S_{mix}$, as a function of model parameters $\alpha$ and $\beta$, mirrors the criticalities exhibited by the minimum-energy configuration, i.e. it is discontinuous at transition 1-2 and 1-3 while it is continuous at transitions 2-3 (apart from the special case $\beta=1$) and 3-4. For this reason, it constitutes a valid indicator to capture the occurrence of mixing-demixing transitions. Note that $S_{mix}$ is particularly advantageous in that it allows one to represent the system critical behavior by avoiding its description in terms of six boson populations $x_{j}$ and $y_{j}$.

Phases and free energy.

Being the temperature zero, the free energy $F=E-TS$ coincides with the internal energy $E$, i.e. the ground-state energy. The latter, within the CVP approach, can be computed by means of the effective potential (3) which, for the four boson distributions corresponding to the four phases, exhibits a specific dependence law on model parameters $\alpha$ and $\beta$ given by

[TABLE]

[TABLE]

The graphic representation of these expressions (see first panel of Fig. 4), shows that $V_{*}$ is indeed continuous everywhere and, in particular, across the transitions. Nevertheless, in agreement with the previous mixing-entropy analysis, the non-analytic character of the ground-state energy does emerge if one computes the first and the second derivative of $V_{*}$ with respect to $\alpha$, which is regarded as a control parameter. As depicted in the second and in the third panel of the aforementioned figure, in fact, at transitions 1-2 and 1-4, the first derivative is discontinuous, while at transitions 2-3 and 3-4, the first derivative is continuous but the second one is discontinuous (except for the special case $\beta=1$, where the first derivative is discontinuous at transition 2-3).

Finite-size effects on the mixing-demixing transitions

As already mentioned, the four different phases emerge at their clearest only in the large-populations limit, i.e. when $U_{a}N_{a}/T_{a}\to\infty$ and $U_{b}N_{b}/T_{b}\to\infty$. Such phases are still recognizable also in the more realistic case of finite-size systems (i.e. systems with limited numbers of atoms and featuring non-vanishing hopping terms), although the phase diagram presented in Fig. 1 gets blurred and deformed. The effects of finite ratios $U_{a}N_{a}/T_{a}$ and $U_{b}N_{b}/T_{b}$ can be summarized as follows:

• •

The minimum-energy configuration $(\vec{x},\vec{y})$, regarded as a function of model parameters $\alpha\in(0,3)$ and $\beta\in(0,1)$, is continuous and, as a consequence, so is the entropy $S_{mix}$ (comparing Fig. 3 and Fig. 5, one can notice that the jump discontinuities give way to smooth junctions).

• •

The fully-mixed phase is favored by the presence of non-negligible hopping amplitudes, its border being given by the inequality

[TABLE]

The latter represents the condition under which the Hessian matrix associated to effective potential (2) (and evaluated at point $x_{j}=y_{j}=1/3$, with $j=1,\,2,\,3$) is positive definite. Walking away from the large-populations limit, the right-hand term rises above the value $1$, thus determining an enlargement of phase 1 at the expenses of the neighboring phases (see the enlargement of the blue region in Fig. 5 moving from the left panel to the right panel). It is worth mentioning that (the twin-species limit of) condition (6) was found to define the region of parameters’ space where the spectrum of Bogoliubov quasiparticles is well-defined and not degenerate[39].

• •

As intuition suggests, the fully demixed phase is mined by the presence of hopping processes between the wells and therefore occurs for higher values of $\alpha$ (see the shrinking of the purple region in Fig. 5 moving from the left panel to the right panel). This clearly reminds the physics of the Superfluid-to-Mott-insulator transition [40], where a stronger on-site interaction is required to observe the transition when a larger hopping rate is present.

• •

Interestingly, increasing the hopping amplitudes, phase 4 not only enlarges, but invades the $\beta>1/2$ region, crushing and shrinking phases 2 and 3 (moving from the left to the right panel of Fig. 5, one can see that the green region enlarges at the expenses of the orange and of the purple regions).

Toward an experimental realization: a "real" mixture and beyond point-like potential wells

In the previous section we have discussed the mixing properties of an asymmetric binary mixture within the CVP of the BH model. In order to get closer to an experimental realization, one has to consider also the non-pointlike character of potential wells and the actual tunability of the scattering lengths in atomic systems.

Despite the large tunability of the scattering length via magnetic Feshbach resonances [41], this knob does not allow for an independent tuning of the intra- and the inter-species scattering lengths. Some mixtures have extremely large background scattering lengths [42], very narrow resonances far beyond experimental magnetic field stability [42], overlapping resonances [43], and promising broad inter-species resonances where, however, the inter-species scattering length is excessively large or has the wrong sign [7].

In this prospective, ${}^{23}\mathrm{Na}+^{39}\mathrm{K}$ constitutes an ideal mixture. The spin mixture $|f=1,m_{f}=-1\rangle_{Na}+|f=1,m_{f}=-1\rangle_{K}$, where $f$ is the atomic total angular momentum and $m_{f}$ its projection on the quantization axis, has been recently doubly Bose-condensed and presents a favorable window of magnetic field (from 90 to 150 G) where the combination of two intra- and two inter-species resonances allow for a smooth tuning of the scattering lengths [43] (see left panel of Fig. 6). Additionally, in this range of magnetic field, three-body losses [44] are limited thanks to the small value of $a_{\mathrm{Na,K}}$. In our case, a change of the magnetic field corresponds to a shift along a line in the phase diagram of Fig. 8, where $\alpha$ and $\beta$ (horizontal and vertical axis respectively) have been rewritten in terms of MF parameters. In the region of interest, the scattering lengths of sodium and potassium are almost constant and the line will be almost horizontal, crossing the different phases at a value $\beta$ mainly determined by the atom-number ratio.

From now on, the formulas will be strictly linked with the real experimental setup and labels $a$ and $b$ will stand for sodium and potassium respectively. Parameters $m_{a}=22.99$ u and $m_{b}=39.10$ u are therefore the atomic masses of such elements, where u stands for unified atomic mass unit.

Recent developments in single-site resolution detection of single atoms and the realization of diffraction limited potentials have pushed forward the capabilities of handling few atoms in tweezers systems [45]. The trimer potential consists (see right panel of Fig. 6) of three gaussian traps whose centers correspond to the vertexes of an equilateral triangle with edge $d$. The realization of such a potential is based on the projection of three dipole-trap beams from the vertical direction of our experimental apparatus which features a large numerical aperture objective with demonstrated high resolution capability [46]. Being the wavelength 1064 nm, we consider realistic values for the width of the three beams ($\sigma=1,00\,\mu m$ ) and for the distance of their centers, $d=2\mu m$. For simplicity, we consider monochromatic tweezers, despite the large difference between sodium and potassium trap frequencies resulting from the different atomic polarizability. A bichromatic sheet of light on the horizontal direction or a time-averaged light potential provide confinement in the vertical direction and compensate the presence of a gravitational sag.

The GPE solution in a real system

In pursuing an actual experimental realization, one must interface the rather essential BH-like model (1) described within the CVP (see effective potential (2)) with experimentally-accessible parameters and measurable quantities. The bridge between theory and experiment is provided by a MF treatment of the problem. This approach also comes with a number of improvements to the CVP and, while confirming the fundamental results obtained therein, shines light on some possibly tricky aspects which were set aside by the pointlike approximation of potential wells and inherently present both in the BH model and in the CVP.

In the MF picture, the system under analysis can be effectively modeled by two stationary GPEs

[TABLE]

[TABLE]

with normalization conditions

[TABLE]

enforcing particle-number conservation in both atomic species. Formulas

[TABLE]

give the intra-species interactions and

[TABLE]

the inter-species coupling being $\mu_{ab}=(m_{a}^{-1}+m_{b}^{-1})^{-1}$ the reduced mass. These last two quantities can be tuned by means of the applied magnetic $B$ field, which, acting on scattering lengths $a_{b}$ and $a_{ab}$, therefore plays the role of the effective control parameter, capable of triggering the mixing-demixing transitions. Here we use $a_{a}=52\,a_{0}$, $a_{b}\in[7.6,8.1]\,a_{0}$ and $a_{ab}\in[10.2,48.5]\,a_{0}$. Term $V_{ext}(\vec{x})$ represents the optical potential corresponding to the ring-trimer geometry. One can generate (experimentally achievable) parameters’ vectors ($N_{a}$, $N_{b}$, $g_{a}$, $g_{b}$, $g_{ab}$) both by tuning the magnetic field $B$ (the latter has an impact on ratios $g_{ab}/g_{b}$ and $g_{b}/g_{a}$) and by setting the boson population asymmetry $N_{b}/N_{a}$. On the basis of such sets of parameters, we determined (see Methods for details concerning the numerical simulation) the eigenfunctions of coupled, non-linear, eigenvalue problems (7) and (8), which minimize the Hamiltonian functional

[TABLE]

[TABLE]

representing, in turn, the total energy of the system. Being qualitatively different, such minimum-energy solutions can be classified into four categories which indeed correspond to the already discussed four phases. In order to better appreciate the great agreement between the results obtained by means of the CVP and those obtained within the MF treatment, a representative of each of these categories is shown in Fig. 7.

With reference to such figure, the first panel represents the category of ground-states belonging to phase 1, the one exhibiting perfect mixing. The two condensates, in fact, clearly overlap in each of the three wells and both boson distributions share the same $Z_{3}$ symmetry of the system. Notice also that, correctly, sodium (depicted in blue), being lighter, is more delocalized with respect to potassium (depicted in red) which, being heavier, is more strongly confined. The second panel depicts a ground-state belonging to phase 2, where one has demixing in two wells and mixing in the remaining one. Notice, in fact, that the first (second) well includes just potassium (sodium), while the third well includes both. The third panel shows a fully-demixed ground state, i.e. a representative of phase 3: potassium occupies the first well while sodium spread in the second and in the third one. Eventually, the fourth panel shows a representative of phase 4, the one characterized by the clotting of potassium in a single well and by the presence of a non-zero fraction of sodium therein.

Importantly, for all parameters’ sets used in our numerical simulations, we observed that the phase separation mechanism is such that the two condensates redistribute among the wells and not within the wells. Nevertheless, since in a single harmonic trap different phase separation mechanisms (e.g. hemispheric-like or spheric-shell-like) can be triggered by different potential and interaction strengths[47], we expect the interplay of these parameters to play an even more crucial role in multiple-trap systems. Due to its remarkable complexity, the analysis of this phenomenology will be discussed in a separate work.

It is worth noticing that, when trap parameters and interaction strengths are such that phase separation does occur via an among-well boson redistribution, the terms proportional to $g_{a}$, $g_{b}$ and $g_{ab}$ in energy functional (10) quite naturally exhibit the same structure of the terms proportional to $U_{a}$, $U_{b}$ and $W$ in effective potential (2). In this sense, parameters $\alpha$ and $\beta$ (see formulas (4)) correspond to, respectively, quantities $g_{ab}/\sqrt{g_{a}g_{b}}$ and $N_{b}/N_{a}\sqrt{g_{b}/g_{a}}$, and therefore span the MF counterpart of the phase diagram of Fig. 1. To link the phase diagram to experimental accessible parameters, in Fig. 8 we also plot the four phases by mapping $\beta$ into $N_{b}/N_{a}$ and $\alpha$ to the magnetic field $B$ thanks to the known values of scattering lengths shown in Fig. 6. Accordingly, gray lines in the both left panels of Fig. 8 play the role of the red boundaries separating the four phases in Fig. 1.

As $g_{b}$ and $g_{ab}$ depend on the applied magnetic field $B$, by sweeping the latter, one can generate pathways through the phase diagram. Moreover, one can vertically translate them simply by tuning the ratio $N_{b}/N_{a}$. For many such pathways, practically constituting extended sets of n-tuples ($N_{a}$, $N_{b}$, $g_{a}$, $g_{b}$, $g_{ab}$), we found the solutions of equations (7) and (8) which minimize energy functional (10). Two examples of this analysis are shown in Fig. 8, where we plot the ground state energy $E_{0}:=\mathcal{H}[\varphi_{a,0},\varphi_{b,0}]$ as a function of the actual control parameter $B$ along two specific pathways through the phase diagram. In remarkable analogy with the discussion relevant to effective potential (3) within the CVP, it is clear that $E_{0}$ is a continuous function of $B$ and that its first derivative $\mathrm{d}E_{0}/\mathrm{d}B$ is discontinuous across transitions 1-2 and 1-4 (see rectangular markers in Fig. 8). Moreover, although not visually obvious, one can verify that the second derivative $\mathrm{d}^{2}E_{0}/\mathrm{d}B^{2}$ is discontinuous at transitions 2-3 and 3-4 (see circular markers in Fig. 8). Correctly, one does not observe phase transitions exactly at grey boundaries shown in the phase diagrams. This is due to the presence of non-zero hopping terms which smooths and deforms the zero-hopping phase diagram. On the right panels of Fig. 8 we plot the the mixing entropy $S_{mix}$ relative to the two paths previously considered and the presence of different mix/demix phases is confirmed. This is very promising in the outlook of an experimental realization as $S_{mix}$ represents the most obvious quantity obtainable from the measured atom number in each wells, irrespective of the observed permutation. We associate the presence of small fluctuations at small values of $S_{mix}$ in both panels to numerical resolution resulting from the strong demixing and the consequent very small atom number in the wells. One could consider these fluctuations as a statistical error on the simulations in that range of magnetic field strength.

Upon the associations $\alpha\leftrightarrow g_{ab}/\sqrt{g_{a}g_{b}}$ and $\beta\leftrightarrow N_{b}/N_{a}\sqrt{g_{b}/g_{a}}$, our MF analysis, based on realistic model parameters involving sodium and potassium atoms, has therefore led to the same phases and the same mixing-demixing transitions evidenced within the CVP and thus corroborates, both qualitatively and quantitatively, the predictions about the phase-separation mechanism obtained therein, offering, also, a viable path towards an actual experimental realization.

Other experimental aspects

To strengthen the feasibility of our proposal, we consider also a possible experimental sequence and problems related to the limited lifetime of the sample because of three-body losses.

The starting point of the experimental sequence will be the creation of the degenerate samples in a cross dipole trap [43] at a magnetic field of about 154 G and the following ramp of the magnetic field to the zero of the inter-species scattering length at 117 G. This allows one to maximize the lifetime of the sample against three-body losses. The mixture can be adiabatically loaded into the trimer potential by rising the intensity of the tweezers and decreasing the intensity of the crossed dipole trap. A ramp of the magnetic-field strength to the target value follows, and the system is let mix/demix. After a fast ramp up of the tweezers intensity to stop hopping between the wells, the tweezers will be separated and the atom number of each species in the three wells will be detected with absorption imaging after a short time of flight. We expect that on each single realization, the system will lie in one of the possible cyclic permutations of the typical boson distributions. Post analysis will reveal the mixing/demixing depending on the applied magnetic field.

Once the interaction parameters have been fixed, mixing and demixing properties of the mixture do not depend on the absolute number of atoms but on their ratio through $\beta$. Envisioning an experimental realization, the atom number becomes critical as appearing with the third power into the three-body loss rate coefficient, which rules the lifetime of the mixture in the tweezers’ wells. In the calculations we chose the order of magnitude of atom numbers in such a way that the lifetime exceeds the smallest of the two hopping time ($<20$ms, see Methods) of an order of magnitude. To this goal we calculate three-body densities from GPE simulations and consider three-body loss rate coefficients on the order of $10^{-41}$ m6/s, taken from experimental measurements performed the Hannover setup and compatible with literature values [48].

Methods

Numerical solution of stationary GPEs.

Equations (7) and (8) constitute a non-linear, coupled, eigenvalue problem. Its solution is obtained by means of a numerical technique based on the damped oscillating particle method[49]. Preliminary, in order to lighten the computational burden, we enact a dimensional reduction (in the spirit of the scheme discussed in [50, 51]) to turn the original genuinely 3D problem into a 1D one. This is obtained by means of factorizations

[TABLE]

[TABLE]

where $\psi_{a}$ and $\psi_{b}$ are such that

[TABLE]

and are therefore normalized to unity, as a consequence of normalization conditions (9). These substitutions lead to equations

[TABLE]

[TABLE]

(notice that wavefunctions labels "$a,x$" and "$b,x$" are understood) whose effective parameters are

[TABLE]

[TABLE]

where $L_{y}$ and $L_{z}$ are chosen in a such a way that $L_{y}\approx L_{z}\approx 2\sigma\approx 2\mu m$. Moreover, the optical-potential terms read:

[TABLE]

where $x_{i}$’s are the positions of the centers of each gaussian trap, expressed in terms of the curvilinear abscissa describing the perimeter of the equilateral triangle formed by the $x_{i}$’s themselves (see yellow line in the right panel of Fig. 6).

Resorting to the damped oscillating particle method[49], one can turn the search for the ground state of stationary non linear Schrödinger equations (11) and (12) into the steady-state solution of the following damped oscillators:

[TABLE]

[TABLE]

where the possible differences between left and right members of equations (11) and (12) act as forcing terms, and where the presence of a non-zero damping $\eta$ guarantees the convergence to a steady-state solution. Upon recasting equations (13) and (14) into pairs of first-order equations by means of substitutions

[TABLE]

one can solve the resulting dynamical system with standard finite difference methods[52]. More specifically, we made use of the leapfrog integrator [53], iterating it until equations (11) and (12) are satisfied up to relative errors $<10^{-3}$ (in order to quantify the mismatch between left and right members, the 2-norm has been employed).

The optical potential

Acousto-optical modulators will generate the three beams and particular attention will be dedicated to calibrate the three wells in order to verify their being identical. To this goal, the use of digital mirror device (DMD) can be fruitful to compensate aberrations and defects in the potentials as already shown in other experiments[54]. We consider also the possibility to obtain trimer potentials based on triangular optical lattice. In this case the tighter confinement suggests that a much smaller atom number per site must be considered.

Species-specific trap frequencies are considered in the calculations as well as the different hopping rate coming from the different height of the inter-well barriers. These correspond to $\omega_{a}\approx 2\pi\times 5600$ Hz, $\omega_{b}\approx 2\pi\times 6900\,Hz$. The depth of the wells is chosen to allow the hopping times {1/$T_{a}$,1/ $T_{b}$} to be smaller than 20ms for both species; estimated hopping rate are evaluated by calculating overlap integrals and one has {$T_{a},T_{b}$} $\approx 2\pi\times\{1000,50\}\,Hz$. Excessively large hopping rate would smooth out the phases previously discussed, while too small will not allow mix/demix in a time shorter than the mixture lifetime. Real hopping rate will be experimentally determined to overpass limitations coming from simple overlapping integral approximation.

Experimental stability

One requires a large experimental stability of the three quantities involved in the mix/demixing phases. These are the atom number (preparation and detection), the trimer potential depth and the magnetic field strength. This last is below 30 mG in the Hannover setup [55] and allows a smooth tuning of the scattering length well within the resolution requested to observe the phase transition shown in Fig. 8. The optical-potential stability can be assured by advanced feedback techniques [56] while the stability of the atom number by real-time analysis [57].

Acknowledgements

A.R. and V.P. wish to thank Magnus Ögren and Edoardo Paganoni for their valuable help in the numerical methods. A.R. and A.Z wish to thank Luis Santos and Silke Ospelkaus for fruitful discussions. A.Z. gratefully acknowledges financial support from the European Research Council through ERC Starting Grant POLAR and from the Deutsche Forschungsgemeinschaft (DFG) through CRC 1227 (DQ-mat), project A03 and FOR2247, project E5.

Author contributions statement

A.R. performed analytic and numerical calculations. A.Z. took care of the experimental aspects. A.R., A.Z. and V.P. have analyzed the results and equally contributed in writing and reviewing the manuscript. V.P. supervised the work.

Additional information

Competing interests: The authors declare no competing interests.