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

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 the bosonic mixture 23Na + 39K, thanks to the large tunability of their intra- and inter-species scattering lengths.

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 which comes with the two constraints ∑ = ∑ = x y 1 i i i i , 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 UN /( ) a a 2 , reduces to an expression including only two effective parameters, 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 Fig. 1 already includes all four phases that the sys-www.nature.com/scientificreports www.nature.com/scientificreports/ tem's ground state can exhibit (notice, in fact, that if β > 1, one can swap species labels a and b and therefore come back to region ). In each phase, the configuration → → x y ( , ) which minimizes effective potential (3) is marked by a specific dependence on the effective model parameters (4). Analytic expressions x j (α, β) and y j (α, β) have been derived by means of an exhaustive exploration of the polytope-like domain of function (3) (see 27 for a detailed description of the method), and are presented here, together with the essential features of the associated phase: 1. In phase 1, being α < 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, ∀j = 1, 2, 3. 2. Phase 2 occurs for 1 < α < 2β, that means for intermediate values of α 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 3. Phase 3, occurring for sufficiently high values of α and not too asymmetric species (more specifically, for 1/ (2α) < β < α/2), features a completely demixed scenario, such that species b clots in one well, while species a equally spreads in the other two sites: 4. Phase 4 occurs for sufficiently high values of α and sufficiently asymmetric species (i.e. for α > 1 and β < 1/ (2α)). A good degree of asymmetry can be achieved, for example, if species b constitutes an impurity with respect to the majority species a, meaning that  N N b 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:  Figure 1. Phase diagram of a binary mixture in a ring trimer. Each of the four phases corresponds to a different functional relationship between the minimum-energy configuration and effective model parameters (4). This circumstance has a direct impact on the mixing properties of the system and on its ground-state energy. Red dashed (solid) line represents a mixing-demixing transition across which the components x j and y j of the minimum-energy configuration feature a jump discontinuity (are continuous). The plot has been obtained by means of a fully-analytic minimization of potential (3) under the constraints www.nature.com/scientificreports www.nature.com/scientificreports/ αβ αβ αβ 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 22,27,30,35 . For example, in phase 3, the ground-state is of the type | 〉 ≈ + + ( ) , an expression where symbol "≈" 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 β < 1, the minimum-energy configuration → → x y ( , ), as a function of model parameters α and β, 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, β = 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 → → In comparison to the one used in 37,38 , 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, log2 ≈ 0.69 in phase 1 (perfect mixing). More generally, S mix , as a function of model parameters α and β, 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 www.nature.com/scientificreports www.nature.com/scientificreports/ β = 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 α and β given by 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 α, 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 β = 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 → ∞ and U b N b /T b → ∞. 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:  www.nature.com/scientificreports www.nature.com/scientificreports/ • The minimum-energy configuration → → x y ( , ), regarded as a function of model parameters α ∈ (0, 3) and β ∈ (0, 1), is continuous and, as a consequence, so is the entropy S mix (comparing Figs 3 and 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 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 α (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 β > 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).
Quantum Results. The CVP approach can be safely applied when the number of particles is large enough, i.e. when it is possible to replace quantum mechanical operators with their -number counterpart. Conversely, for small boson populations, the adoption of this approximation scheme is less justified, and a fully quantum approach is needed. Such approach, involving the exact numerical diagonalization of Hamiltonian (1), can be computationally demanding since the dimension of the fixed-boson-number Hilbert space of states (the latter being of the type | → → 〉 = | n m n n n , : , 3 ), rapidly increases with increasing number of particles. Devoting the necessary attention to the numerical protocol, one can still handle the matrix representation of Hamiltonian (1) when N a = N b = 15 and therefore accomplish its exact diagonalization, thus obtaining the energy levels i and the energy eigenstates |ψ i 〉. From these outputs, one can build several indicators which, in spite of the still small numbers of employed particles, already highlight the occurrence of the same critical phenomena and mixing-demixing transitions predicted within the CVP scheme. Among the various possibilities, and in order to offer a direct comparison with its CVP-counterpart, we introduce the quantum version of the entropy of mixing,  www.nature.com/scientificreports www.nature.com/scientificreports/ upon substitutions n i /N a → x i and m i /N b → y i . Sweeping 41 model parameters W and U b , one can effectively reconstruct the behaviour of  S mix in the (α, β)-plane. The result, illustrated in Fig. 6, witnesses the presence of the discussed four quantum phases also in this purely quantum scenario and, being in great agreement with its CVP-counterpart (compare with Fig. 5), validates the predictions given by the semiclassical approximation scheme. In particular, the latter proves to be rather accurate already with quite limited numbers of particles (N a = N b = 15).

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 42 , 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 43 , very narrow resonances far beyond experimental magnetic field stability 43 , overlapping resonances 44 , 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 Na + 39 K constitutes an ideal mixture. The spin mixture | f = 1, m f = −1〉 Na + | f = 1, m f = −1〉 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 44 (see left panel of Fig. 7). Additionally, in this range of magnetic field, three-body losses 45 are limited thanks to the small value of a Na,K . In our case, a change of the magnetic field corresponds to a shift along a line in the phase diagram of Fig. 9, where α and β (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 β 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 46 . The trimer potential consists (see right panel of Fig. 7) 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 47 . Being the wavelength 1064 nm, we consider realistic values for the width of the three beams (σ = 1,00 μm) and for the distance of their centers, d = 2 μ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 Figure 6. Entropy of mixing associated to the ground state |ψ 0 〉, obtained by means of formula (7) sweeping model parameters α and β. The following values/intervals have been chosen:   Na K mixture given in the Bohr radii a 0 . Yellow dashed, red dotted and blue continuous line refer to sodium a Na , potassium a K and sodium-potassium a Na,K scattering lengths. Labels refer to points considered in the next section.
Note the approximately constant intra-species scattering lengths and the inter-species a Na,K varying from <a K to >a Na . Right panel: Schematic representation of the three gaussian traps constituting the ring-trimer potential. Yellow dashed line represents the equilateral triangle formed by the trap's center. Such triangle can be associated to the curvilinear abscissa x present in equations (12) and (13) (and therefore automatically endowed with periodic boundary conditions). 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 (8) and (9), which minimize the Hamiltonian functional  b  ab  a  b   2  2  ,  3  43   2  2  ,  3  43  2  2 3 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. 8. 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  (8) and (9) into functional (11)). Rectangular (circular) markers indicate the points where the first (second) derivative of E 0 with respect to B is discontinuous. Numbers 1-4 not only refer to the phase crossed by the pathway but specifically correspond to the minimum-energy solutions 1-4 illustrated in Fig. 8. Right column's panels: entropy of mixing computed along each pathway. Normalized boson populations appearing in formula (5) are computed by integrating the (square modulus of) the solutions of equations (8) and (9) in the three spatial domains associated to the potential wells. (2019) 9:6908 | https://doi.org/10.1038/s41598-019-43365-6 www.nature.com/scientificreports www.nature.com/scientificreports/ 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 48 , 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 (11) 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 α and β (see formulas (4)) correspond to, respectively, quantities g g g / ab a b and N N g g / / b a b 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. 9 we also plot the four phases by mapping β into N b /N a and α to the magnetic field B thanks to the known values of scattering lengths shown in Fig. 7. Accordingly, gray lines in the both left panels of Fig. 9 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 (8) and (9) which minimize energy functional (11). Two examples of this analysis are shown in Fig. 9, where we plot the ground state energy  ϕ ϕ = E : [ , ] , 0 , 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 dE 0 /dB is discontinuous across transitions 1-2 and 1-4 (see rectangular markers in Fig. 9). Moreover, although not visually obvious, one can verify that the second derivative d 2 E 0 /dB 2 is discontinuous at transitions 2-3 and 3-4 (see circular markers in Fig. 9). 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. 9 we plot the the mixing entropy S mix relative to the two paths previously considered and the presence of different mixed/demixed 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 α ↔ g g g / ab a b and β ↔ N N g g / / b a b 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 44 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 β. 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 m 6 /s, taken from experimental measurements performed the Hannover setup and compatible with literature values 49 . www.nature.com/scientificreports www.nature.com/scientificreports/ Methods Numerical solution of stationary Gpes. Equations (8) and (9) constitute a non-linear, coupled, eigenvalue problem. Its solution is obtained by means of a numerical technique based on the damped oscillating particle method 50 . Preliminary, in order to lighten the computational burden, we enact a dimensional reduction (in the spirit of the scheme discussed in 51,52 ) to turn the original genuinely 3D problem into a 1D one. This is obtained by means of factorizations where ψ a and ψ b are such that and are therefore normalized to unity, as a consequence of normalization conditions (10). These substitutions lead to equations (notice that wavefunctions labels "a, x" and "b, x" are understood) whose effective parameters are where L y and L z are chosen in a such a way that L y ≈ L z ≈ 2σ ≈ 2 μm. Moreover, the optical-potential terms read: 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. 7).
Resorting to the damped oscillating particle method 50 , one can turn the search for the ground state of stationary non linear Schrödinger equations (12) and (13) into the steady-state solution of the following damped oscillators: where the possible differences between left and right members of equations (12) and (13) act as forcing terms, and where the presence of a non-zero damping η guarantees the convergence to a steady-state solution. Upon recasting equations (14) and (15)  one can solve the resulting dynamical system with standard finite difference methods 53 . More specifically, we made use of the leapfrog integrator 54 , iterating it until equations (12) and (13) 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 55 . We consider also the possibility to obtain trimer potentials based on triangular www.nature.com/scientificreports www.nature.com/scientificreports/ 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 ω a ≈ 2π × 5600 Hz, ω b ≈ 2π × 6900 Hz. The depth of the wells is chosen to allow the hopping times {1/T a , 1/T b } to be smaller than 20 ms for both species; estimated hopping rate are evaluated by calculating overlap integrals and one has {T a , T b } ≈ 2π × {1000, 50} Hz. Excessively large hopping rate would smooth out the phases previously discussed, while too small will not allow mixing/demixing 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 56 and allows a smooth tuning of the scattering length well within the resolution requested to observe the phase transition shown in Fig. 9. The optical-potential stability can be assured by advanced feedback techniques 57 while the stability of the atom number by real-time analysis 58 .