Ground-State Properties and Phase Separation of Binary Mixtures in Mesoscopic Ring Lattices

We investigated the spatial phase separation of the two components forming a bosonic mixture distributed in a four-well lattice with a ring geometry. We studied the ground state of this system, described by means of a binary Bose–Hubbard Hamiltonian, by implementing a well-known coherent-state picture which allowed us to find the semi-classical equations determining the distribution of boson components in the ring lattice. Their fully analytic solutions, in the limit of large boson numbers, provide the boson populations at each well as a function of the interspecies interaction and of other significant model parameters, while allowing to reconstruct the non-trivial architecture of the ground-state four-well phase diagram. The comparison with the L-well (L=2,3) phase diagrams highlights how increasing the number of wells considerably modifies the phase diagram structure and the transition mechanism from the full-mixing to the full-demixing phase controlled by the interspecies interaction. Despite the fact that the phase diagrams for L=2,3,4 share various general properties, we show that, unlike attractive binary mixtures, repulsive mixtures do not feature a transition mechanism which can be extended to an arbitrary lattice of size L.


Introduction
The demixing effect in bosonic binary mixtures trapped in optical lattices has attracted a fair amount of attention in the last two decades due to both the rich phenomenology stemming from the spatial separation of mixture components [1,2] and to the considerable complexity that the demixing mechanism features. The latter, basically dominated by the interplay of density-density interspecies and intraspecies repulsions with the particle number of each component, bears memory of the fragmentation of the boson fluid in lattice multiwell structures and of the deeply quantum character of microscopic bosonboson interactions. These features were taken into account by representing boson mixtures and the relevant phenomenology within the well-known second-quantization picture [3] characterizing the Bose-Hubbard (BH) model [4,5]. Experimentally, binary mixtures were realized by either combining different atomic species [6] or by means of components corresponding to different internal states [7,8] of a single atomic species.
Spatial phase separation greatly enriches the phase diagram of mixtures [9,10] characterized by new types of superfluidity and incompressible states [11,12], rules the coexistence of different phases in the presence of trapping potential [13,14], it is involved in the formation of quantum emulsions [15,16] and represents a robust property of mixtures able to survive at non-zero temperature and in spite of the confinement effect [17].
Recently, the miscibility properties of binary mixtures were investigated in small-size optical lattices to obtain a deeper insight into the influence of space fragmentation on the demixing mechanism. The interest for this class of systems was doubly motivated. First, they are within the reach of current experimental techniques. Designed almost fifteen years ago [18], the realization of traps with a ring geometry has been studied in [19,20] to support the development of atomtronics devices, while the technique used for realizing a double well [21,22] indeed supplies a viable scheme to construct linear open arrays of potential wells. Second, the relatively small number of local order parameters necessary to describe the microscopic state of the system often allows one to analytically approach the study of low-energy states and of their properties despite the well-known non-integrable character of BH lattice systems. This, in turn, makes it possible to obtain new information about the miscibility properties and the separation mechanism.
The critical condition for which the interplay of tunneling parameter, boson-boson interactions and boson numbers of the two components links demixing effect and spectral collapse has been analytically investigated in the double-well system [23,24] in the case of both repulsive and attractive components, while the role of inhomogeneity on demixing, due to the trapping effect, has been evidenced by using multiconfigurational time-dependent Hartree method [25]. Purely quantum indicators, such as entanglement and residual entropies [26], and a nontrivial geometric phase [27] have further confirmed the critical behavior distinguishing spatial separation, while the impact of a chaotic dynamical behavior on the robustness of spatial separation has been explored in both lattice structures [28] and in a harmonic trap [29].
More recently, demixing has been investigated for a ring triple well in asymmetric configurations where interspecies and intraspecies repulsions are different, and the imbalance of the two boson populations is non-zero. The higher degree of space fragmentation due to the third well has been clearly shown to imply a more structured phase diagram, and implicitly, a more complex separation mechanism [30,31]. A parallel study on attractive mixtures, whose production has been recently studied in [32] together with the formation of droplet states, has revealed how the transition from mixed fully delocalized components to mixed components confined in a single well (supermixed state) [33] is controlled by the interspecies interaction and the population imbalance, but it is totally independent from the number of potential wells forming the lattice.
This leads to the central problem addressed in this paper. Unlike the attractive mixture, the demixing mechanism of repulsive components seems to feature a strong dependence from boson fragmentation, namely from the number of lattice sites. In particular, any prediction about the structure of the four-well phase diagram can be evinced, apparently, neither from the knowledge of the three-well case nor from some simple or intuitive rule working for an arbitrary site number. This message already emerges from the comparison of the two-well with the three-well phase diagram. Even in the simple case of identical components, the transition from the mixed to demixed state features the occurrence of a third intermediate phase absent in the two-well case. With L = 4 (and presumably, with a larger number L of wells), the understanding of this mechanism becomes even more problematic.
In this paper, we applied the coherent-state (CS) variational picture to the secondquantized BH model describing the four-well system. This allows us to obtain the groundstate equations and thus to reproduce the mechanism governing the formation of the new phases starting from the weak-interaction regime in which the system was fully mixed and delocalized. We shall compare the phase diagram of the three-well BH model (trimer) with that of the two-well (dimer) and the four-well BH model, with the aim to evidence both the common features shared by the L-well systems for L = 2, 3, 4, and the possible diversity they exhibit when demixing, the transition from to full-mixing regime to the complete space separation, is enacted.
Our exploration will be developed by considering the large-population limit, a condition equivalent to assuming vanishingly small effective hopping amplitude in the equations describing the mixture ground state. This approach has been proven to be very effective in investigating both of the three-well phase diagrams [30,31], and the formation of supermixed states in attractive mixtures [33]. The fact that it allows a completely analytic derivation of the phase diagram represents its main feature.
The CS picture of the BH model and the dynamical equations describing the systems in terms of local order parameters is reviewed in Section II. Then, the demixing mechanism of the dimer case, never explored within the current CS-based analytic approach, is discussed in Section III, both to highlight phase separation in its most elementary form, and to exemplify the phase diagram derivation method then applied to the four-well system. Subsection IIA is devoted to revisiting the results concerning the three-well system. Section IV provides a detailed description of the solutions corresponding to five phases occurring in the four-well phase diagram and to the reconstruction of the latter based on the comparison of the ground-state energy of each phase. Finally, Section V is devoted to investigating the onset of critical behaviors and to their characterization by means of the mixing entropy and of the location entropy. The comparison of the properties of the L-well phase diagram for L = 2, 3, 4 completes our discussion.

Coherent-State Picture of BH Model
Binary mixtures loaded in a ring lattice formed by L potential wells are well described by the Bose-Hubbard Hamiltonian [3]: wheren ai =â + iâ i andn bi =b + ib i are boson number operators, W is the interspecies repulsion strength describing the coupling between the bosons of two components of the mixture, and:Ĥ is the single-component BH model [5]. Parameters U c and J c represent the boson-boson intraspecies repulsion and the hopping amplitude, respectively, for the component c = a, b. Index i ranges within [1, L], where L is the site number of the ring lattice, and the bosonmode operatorsĉ + i andĉ i satisfy standard commutators [ĉ ,ĉ + i ] = δ i . A distinctive feature of model (1) is that the number operators of each boson component: are conserved quantities, where [Ĥ,N a ] = [Ĥ,N b ] = 0, whose eigenvalues N a and N b describe the boson numbers of the two components.
If the system is in the superfluid regime, an assumption guaranteed by the condition that J c /U c (with c = a, b) is sufficiently large, then the dynamical evolution of the bosonic fluid can be represented by means of the coherent-state variational picture [34,35]. In this picture the state of the system is expressed by means of the factorized trial state: and the Glauber coherent states |a i and |b i describe the physical state of the two components at the i-th well. Within the Glauber-state representation, in fact, the coherent-state labels a i = a i |â |a i , b i = b i |b |b i , with a i , b i ∈ C, can be interpreted as local order parameters while the expectation value n ci = c i |ĉ + iĉ i |c i = |c i | 2 (c = a, b) represents the average local boson population. The time evolution can be shown to be governed by the discrete nonlinear Schrödinger equations: The latter are easily derived from the Hamilton equations: associated to the effective Hamiltonian: one obtains by applying the time-dependent variational procedure to model (1) (see, for example, [34]). If the two components are decoupled, namely W = 0, one obtains two independent systems in which the uniform configuration featuring: with a i = a and b i = b, can be easily shown to correspond to the ground state, as witnessed by the fact that the effective total energy Φ|H c |Φ − µ c N c reaches its minimum value. The motion equations are reduced to the equation pair: determining µ a , µ b thanks to the conservation of the boson numbers N a = L|a| 2 , N b = L|b| 2 . For coupled components (W > 0), the study of the three-well system has clearly shown how the ground state does not necessarily coincide with the uniform configuration so that, in general, the condition a i = a and b i = b cannot be applied any longer. For W > 0, the system: µ a a = U a |a | 2 + W|b | 2 a − J a (a +1 + a −1 ), obtained by substituting the collective-frequency modes (4) in Equations (2) and (3), provides the ground state of the system. Its highly nonlinear character, in general, requires the search of numerical solutions. Interestingly, however, the analysis developed in [30,31] reveals that analytical solutions can be found if one considers the case of population N a and N b to be sufficiently large. This becomes evident by effecting the substitutions The previous equations take the form: where the local densities ρ = |a | 2 /N a ≤ 1, and η = |b | 2 /N b ≤ 1 obey the boson-number constraints.
and the effective chemical potentials and hopping amplitudes: together with the new parameters: are introduced. Parameters α and β are particularly useful because they embody the significant information about the mixture: α > 1 (α < 1) distinguishes the strong (weak) interspecies-interaction regime, while β measures the asymmetry of the mixture with respect to the twin-component case Note that, in principle, both α and β range in the interval [0, ∞]. For β > 1, however, one can swap the labels a and b of the two components, thus, recovering a parameter choice satisfying β ≤ 1. Then, hereafter, the latter inequality is assumed to define the range of β.
For sufficiently large N a and N b , one easily sees how the effective amplitudes τ a and τ b become negligible quantities so that Equations (5) and (6) reduce to: This strategy enables one to reconstruct the essential structure of the ground-state phase diagram in which the large-population condition, entailing τ a , τ b → 0, can be interpreted as a sort of thermodynamic limit according to the statistical-mechanical approach discussed in [36,37]. Its effectiveness is attested by the simple derivation of the phases' landscape for the three-well system for both repulsive [31] and attractive [33] interspecies interaction. It is worth mentioning a similar version of this scheme, recently applied to study the critical properties of few-fermion systems [38], where the limit N a , N b >> 1 is replaced by the equivalent assumption that intraspecies interactions are sufficiently strong.

Phase Separation in the BH Dimer
We derive the characteristic phases of the two-well system for τ a = τ b → 0 and compare them with those occurring in the three-well phase diagram. The two-well case allows us both to highlight the spatial separation mechanism in its most elementary form, and to exemplify the method we apply to determine the ground-state configurations when α and β are varied. In the two-well case, Equation (9) feature three types of solutions: which describe fully mixed, partially demixed and fully demixed components, respectively.
Full-mixing solution. In the first case, after suppressing factors a , and b in system (9), the latter reduces to: entailing twin densities for each component: In view of this, constraints (7) then impose that: thus providing the conditions whereby unknown parameters λ a , λ b can be fixed. One obtains λ a = (1 + αβ)/2 and λ b = (1 + α/β)/2. Equation (11) describe the well-known uniform solution which exists for any value of α and β but loses its status as the ground state for α > 1. In the large-population limit, this inequality in fact coincides with the condition α > 1 + 2τ → 1 (twin components with τ a = τ b = τ were assumed) for which the perturbation of the uniform solution within the Bogoliubov scheme develops unstable oscillations [23,39] and then can no longer represent the ground state. The second-type and third-type solutions, however, offer alternative candidates for the ground-state configuration.
Single-well-mixing solution. In the second case, the assumption that component b is absent in the second well (y 2 = 0 → η 2 = 0), so that only the first well features mixing, implies that: with the constraints 1 = ρ 1 + ρ 2 and 1 = η 1 . The solution is then given by showing how the increasing demixing corresponds to a residual single-well-mixing (SWM) in one of the two wells. Of course, a perfectly equivalent solution is obtained by assuming b 1 = 0 in place of b 2 = 0 which simply entails the exchange of the analytic expressions describing ρ 2 and ρ 1 . The crucial information emerging from solution (12) is that the relevant definition domain must satisfy the inequality β < 1/α to prevent a negative density ρ 1 .
Full-demixing solution. Finally, we consider the third case a 2 = 0 = b 1 , a 1 = 1 = b 2 , which features complete space separation. Imposing the constraints (7) gives: As a consequence, system (9) is easily satisfied in that the equations for a 2 and b 1 are simply eliminated, while the remaining two equations λ a = ρ 1 ≡ 1, λ b = η 2 ≡ 1 serve to fix the effective chemical potentials. The resulting scenario is summarized in Figure 1, where regions 1, 2 and 3 refer to ground-state phases associated to solutions (11)-(13), respectively. The comparison of the three solutions shows that, while SWM solution (12) coalesces with solution (13) along the hyperbole β = 1/α, solution (12) does not match solution (11) for α = 1, thus revealing a discontinuous behavior. Inequality α < 1, rewritten as W < √ U a U b , shows how region 1 represents the weak-interaction regime in which the inter-species repulsion is too small to cause spatial separation. However, the possibility to directly reach region 3 with fully demixed components is given only for β In any other case, β < 1, the crossing of region 2, the intermediate partial-mixing phase, represents the unavoidable process whereby the residual SWM is suppressed by increasing α and region 3 is then reached. This circumstance highlights the important role of the asymmetry degree of the system which, through the inequality β > 1/α, determines the critical value W = √ U a U b /β for which the interaction is strong enough to trigger full demixing.
The ground-state map, namely, the phase diagram in Figure 1, is derived by combining the knowledge of the energies corresponding to each solution with the information about their definition domains. These energies read (Ē σ = E σ /U a N 2 a , σ = 1, 2, 3): One easily checks that:Ē More interestingly, one findsĒ 1 <Ē 2 for α < 1, while for α > 1 one hasĒ 1 >Ē 2 . Then, the ground-state status is transferred to solution (12) for α > 1. This property is confirmed, in parallel, by the fact that the inequality: (12) and (13) are well defined), then excluding the possibility that solution (13) is the ground state. The latter inherits this role in D 3 = {(α, β) : β < 1, β > 1/α} (region 3) whereĒ 3 is the minimum energy.

Comparison with the BH Trimer
The ground-state phase diagram of the BH trimer is illustrated in Figure 2 and includes four phases. These correspond to the definition domains of four different solutions to Equation (9) which represent the lowest energy configurations. Their derivation, discussed in Ref. [31], involves six local densities since = 1, 2, 3.  (15) and (16)) separate phase 1 (uniform solution) and phase 3 (fulldemixing solutions (14)). The comparison with the 2-well phase diagram in Figure 1 shows how the presence of the third well causes the onset of phase 2, representing a second independent way to realize partial demixing.
Phase 1 corresponds to the uniform solution characterized by full mixing and the local densities ρ = η = 1/3, with = 1, 2, 3, of the two components among the three wells. As for the two-well system, the border α = 1 of phase 1 coincides with the upper limit of the weak interaction regime: for α < 1, the inter-species repulsion is too small to cause spatial separation.
Phase 3 in Figure 2 represents the opposite regime characterized by complete demixing: one species clots in one site, while the other spreads the remaining two. One has, for example: This phase is delimited by 1/α < β < α/2. It occurs for values of α sufficiently large and includes both symmetric (β = 1) and asymmetric components (0 < β < 1). Phase 3 shows how, except for β = 1/2, the transition from the fully mixed phase 1 to the fully demixed phase 3 cannot occur without the crossing of intermediate phases (either 2 or 4) in which mixing is progressively suppressed. Phases 2 and 4 intercalate phases 1 and 3. Phase 4, delimited by α > 1 and β < 1/(2α), is quite naturally associated to the phase 2 of Figure 1. Partial demixing of phase 4 exhibits the same space structure of the SWM solution (12): The minority component clots in a single site together with a small fraction of the majority component, however, the remaining fraction of the majority component occupies the other two wells. This is exemplified by the solution: The symmetry of Equation (9) under simple permutations of the local densities shows how, altogether, one obtains three, six and three equivalent realizations of a given solution to (14)-(16), respectively. Phase 2 gives the possibility to follow a second independent path to reach phase 3 from phase 1. It exhibits partial demixing, namely the ground-state solution features two wells with complete demixing (the two components occupy different wells) while the remaining fractions of the two components are mixed in the third well: This phase is delimited by 1 < α < 2β, thus corresponding to intermediate values of the effective interaction α and not overly asymmetric components (1/2 < β < 1). As observed above, phase 4 is essentially equivalent to phase 2 in Figure 1 due to the similarity of the relevant SWM solutions (12) and (15), and of the hyperbolic curve constituting their upper phase boundary. The significant change exhibited by the trimer phase diagram is thus the appearance of phase 2 for 1/2 < β < 1, a ground-state solution absent in the dimer phase diagram. The presence of the third well triggers, for 1/2 < β < 1, a parallel mechanism to reach the full-demixing regime in which the reduction in mixing coexists with the fact that each component is distributed on at least two of three wells. This mechanism significantly modifies for β < 1/2 when, similar to the dimer behavior, the minority component is confined in a single well.

Phase Diagram of the Four-Well System
When the ring lattice is formed by four wells, the equation system (9) depends on eight local densities ρ and η with ∈ [1,4]. The number of solutions considerably increases and before reconstructing the ground-state map, one must identify the set of solutions which cannot represent the ground state whilst representing excited states. For this reason, below we focus our attention on the four solutions which, in addition to the expected uniform solution, result in the formation of the four-well phase diagram. Details about their derivation and the presence of other solutions not fitting the minimum-energy condition are discussed in Appendix A. To better highlight the structure of the phase diagram described in Figure 3 we explicitly mention, at each step, the link between phases and solutions.  (17)) and phases 6 and 7 (full-demixing solutions (20) and (21)) are intercalated by phases 3 and 4 (partial-demixing solutions (18) and (19)) which highlight two different ways to enact spatial phase separation. The comparison with Figure 2 shows how the presence of the fourth well triggers the onset of phases 6 and 7.
Fully mixed solution. This is the uniform solution, which corresponds to phase 1, and is obtained from configuration (A1). It exhibits the usual form: SWM solution. This solution describes phase 3. The space distribution obtained from configuration (A3) predicts that the minority component clots in a single well where a small fraction of the majority component survives causing residual mixing: This configuration, defined for: also characterizes both the two-well (see solution (12)) and the three-well cases (see solution (15)). As discussed in the three-well case, the complete confinement of one component in a single well allows one to regard (18) as the configuration that, before reaching full demixing, features the most pronounced demixing degree.
Weak-demixing solution. Phase 4 is described by the space distribution ensuing from the trial solution (A4). The distribution has the form: and is defined for β < 1/α. It describes a double-dimer structure in which the distribution of the two components characterizing wells 1 and 2 is equal to the distribution in the wells 3 and 4 and corresponds to the dimer solution (12). Such a state, by definition, cannot occur in a system with three wells. Note that going from the two-well to the four-well case highlights how the SWM solution (12) (phase 2 in Figure 1) is now replaced by two different intermediate configurations (solutions (18) and (19)) instead of one due to the larger number of wells. They exhibit different realization of partial demixing.
Twin-dimer full-demixing solution. This is associated to phase 6 and is obtained from configuration (A6). The relevant components' distribution reads: It essentially represents the fully demixed solutions (see solution (13)) for a pair of twin two-well systems.
Other solutions. The equation system (9) provides two further solutions that are physically meaningful and are characterized by partial demixing. These are described by formulas (A9) and (A12) in Appendix A. However, the comparison of their energies with those corresponding to the previous solutions reveals how they represent in any case excited states, and thus cannot be associated to any phase in the four-well phase diagram.

Phase Diagram Derivation and Characteristic Ground-State Energies
The ground-state energiesĒ k = E k /(U a N 2 a ) (index k corresponds to the phase index) of the five phases in Figure 3 are easily obtained by substituting the corresponding expressions of the local densities in the four-well Hamiltonian with τ a = τ b = 0. The boundaries of each phase are recognized by determining the domain in which, for a given k, the conditionĒ k ≤Ē i is verified for any other i = k. As for the dimer case, this information combined with the definition domains of each solution allows one to reconstruct the phase boundaries visible in Figure 3.
The results obtained by applying this procedure are summarized in the following list which describes the ground-state energies and the relevant definition domains: with D 6 = {α ∈ [1, √ 3] : 1/α < β < 1} and D 6 = {α > √ 3 : 1/ √ 3 < β < 1}, and: with: The phase diagram shown in Figure 4 displays the ground-state domains corresponding to this energy scenario. For the sake of clarity, we briefly sketch the essential steps that allows one to reconstruct the four-well phase diagram. One begins by considering the uniform solution (17) which identifies with the ground state for α < 1 since one easily shows that E 1 <Ē k for any k = 1. Then, the fact thatĒ 1 =Ē 3 =Ē 4 , for α = 1, suggests that phases 3 and 4 associated to solutions (18) and (19), respectively, could be the ground state for α > 1. Concerning phase 3, solutions (18) is confirmed to be the ground state because inequalities E 3 <Ē 4 ,Ē 6 ,Ē 7 are verified in the domain α > 1, and β < 1/(3α). In parallel, the validity domains of inequalitiesĒ 4 <Ē 6 andĒ 4 <Ē 7 are found to be: respectively, which highlights, together with α > 1, the definition domain of phase 4 in which solution (19) is the ground state. Note that the two boundaries β = 1/α and β = f (α) (implicitly defined by the previous inequalities) are such that 1/α = f (α) = 1/ √ 3 for α = √ 3. The value β = 1/ √ 3 has a second important role: it represents the straight line separating phases 7 and 6. This clearly emerges from inequalityĒ 7 <Ē 6 whose validity domain is described by β < 1/ √ 3. The recognition of the definition domains in which E 6 and E 7 are the ground-state energies completes the phase diagram exploration.

Discussion
The distinctive property of the four-well phase diagram is the presence of two independent realizations of the full-demixing state associated to solution (20) (phase 6) and solution (21) (phase 7). As a consequence, unlike the three-well phase diagram (Figure 2) where the unique full-demixing regime (phase 3) is reached through phases 2 and 4 for a sufficiently large α, the suppression of mixing in the four-well case can be realized in three independent ways (intuitively described by the phase sequences 1-4-6, 1-4-7 and 1-3-7) leading to the full-demixing phases 6 and 7. More specifically, the emergence of phase 6 significantly modifies the structure of the four-well phase diagram which otherwise would display a structure, with four domains, qualitatively similar to that of the three-well phase diagram for β < 1/ √ 3. A common trait with the two-well and three-well phase diagrams emerges for a sufficiently small β. Similarly to the three-well case for β < 1/2, the four-well phase diagram reproduces, for β < 1/3, the structure of the domain in the two-well phase diagram corresponding to the SWM solution (phase 2). This circumstance suggests that this behavior (and the relevant SWM solution) should also persist for a higher number of wells and sufficiently small β. Note that, for U a = U b , the inequalities defining the upper boundaries β < 1/α, β < 1/2α, and β < 1/3α, of strong-demixing solutions (12), (15) and (18), respectively, become: entailing that, for a given value of α, the access to the lower SWM phase in the three phase diagrams requires that N b is proportional to a decreasing fraction 1/(L − 1) of N a when the well-number L is increased. Another aspect characterizing the four-well phase diagram is the behavior of the local densities at the phase boundaries separating phases and k that we denote with γ k . Considering solutions (17)- (21), one easily discovers that the phase-solution and the phase-k solution match each other only along γ 37 and γ 46 , whereas the ground state exhibits a discontinuous behavior of the local densities when crossing γ 14 , γ 13 , γ 47 , and γ 67 .
The discontinuities observed when going from the full-mixing phase 1 to partialdemixing phases (boundaries γ 14 , γ 13 at α = 1) represents a common feature with the twowell and three-well phase diagrams as well as the continuous character of local densities when crossing the upper boundaries β = 1/[(L − 1)α] of the SWM-solution phases for L = 2, 3, 4. Unlike boundaries γ 47 , and γ 67 of the four-well case, one easily checks that the two-well and three-well phase diagrams do not feature further phase boundaries affected by local-density discontinuity for α > 1.
Interestingly, this feature can be put in relation to the well-known property of thermodynamic potentials that links the presence of a transition to the discontinuities of the partial derivatives of potentials. In the current case, the free energy F = E − TS = E (at temperature T = 0) and its partial derivatives allow one to detect possible critical behaviors at the phase boundaries γ k . Despite the continuous character of energyĒ = E/(U a N 2 a ), well illustrated in Figure 4, one easily discovers that at least one of the following discontinuity conditions: involving first-order derivatives takes place for any pair (k, ) except for (k, ) = (3, 7), (4, 6) (see Figure 5). Concerning the latter cases, however, their continuity is lost when considering second-order derivatives. Instead of representing the behavior of the F derivatives we confirm the presence of the critical behavior by exploiting two indicators of spatial order which, in the case of the three-well systems, were proven to effectively detect the different domains of the phase diagram [31,33]. These are [40,41], the mixing entropy and the location entropy: respectively. The states of the two-well phase diagram well exemplifies the information provided by S mix and S loc . Considering solutions (11)-(13), relevant to the full-mixing, partial-demixing and full-demixing regimes, one finds: respectively, where ln2 is the maximum value of the mixing entropy. As expected, the maximum and minimum mixing degrees are thus associated to the full-mixing and fulldemixing regimes, respectively, while partial demixing corresponds to the intermediate value S mix = 0. Location entropy S loc measures the dispersion of the two components in the two wells. The dispersion reaches its minimum value if both components are localized in the same well. For example, by setting ρ 1 = η 1 = 1, ρ 2 = η 2 = 0m one obtains S loc = 0.
Conversely, the full-mixing state (11) entails maximum dispersion S loc = ln2 while state (13) gives the intermediate value S loc = ln2/2. As applied to the four-well phase diagram, the behavior of S mix , described in Figure 6 clearly shows the discontinuity characterizing the transition from the fully mixed state to the partial-demixing states (phases 3 and 4), when crossing boundaries γ 13 and γ 14 . This indicator also captures the discontinuity relevant to the transition from phase 4 to phase 7 (boundary γ 47 ), but fails to detect the discontinuity between phases 6 and 7. The latter and the relevant boundaries γ 67 are successfully captured by S loc , as shown in Figure 7, which is able to distinguish the different dispersion degree for the two full-demixing states. Quite consistently, neither S mix nor S loc exhibit discontinuous behaviors at γ 37 and γ 46 where, despite the crossing, the local densities of both components remain continuous.  Finite-size effects. The study of mixtures in the ring trimer [31,33] has clearly shown how the information obtained in the limit τ → 0 about the three-well phase diagram ( Figure  2) and its characteristic phases is still qualitatively valid when considering finite-size effects, namely for small but non-zero τ (finite boson populations).
If τ is small, the local-density solutions satisfy Equations (5) and (6) and the jump discontinuities affecting solutions (17)-(21) at certain boundaries γ k are removed, substituted by smooth junctions embodying the continuous dependence of local densities from α and β. This reflects, of course, on the jump discontinuities exhibited by S mix and S loc (function of local densities) which are replaced by rapid but continuous changes inside restricted thin domains. For τ → 0, these domains collapse to the original boundaries γ k , while jump discontinuities of local densities crop up at phase boundaries. In this scenario, the phases of the large-population limit remain well recognizable, confirming the fact that the phase diagrams of Figures 1-3 incorporate the essential information on the ground-state properties of L-well systems.
Conjectures for large-size rings. One can show that certain ground-state configurations suggested by the small-size ring analysis can be extended to rings with a large well number L. To this end, it is useful to remember that ring lattices are distinguished by the fact that L is either even or odd.
Let us first consider the case with odd L = 2s + 1, s ∈ N. We conjecture that, when β is sufficiently close to β = 1, there exists an intermediate partial-demixing phase with a solution exhibiting the same structure of solution (16) (associated to phase 2 in Figure 2) exhibiting almost complete demixing except for a single well. For the sake of simplicity, we assume that the residual mixing occurs at a single well, j = s + 1, so that only ρ s+1 and η s+1 are simultaneously non-zero. This amounts to setting: With this assumption, solving Equation (9) gives ρ 1 = ρ 2 = · · · = ρ s , and η s+2 = η s+3 = · · · = η L with: (remember that L = 2s + 1). One easily checks that the two constraints sρ 1 + ρ s+1 = 1 and η s+1 + sη L = 1 are satisfied. The definition domain is given by which extends the definition domain of solution (16) to L = 2s + 1 > 3. As for the trimer case, the previous solution exhibits a discontinuous behavior at the border α = 1 where it does not match the uniform solution ρ j = 1/L, η j = 1/L, ∀j. Conversely, the energy is expected to be continuous, a behavior confirmed by the fact that the uniform-solution energy:Ē 0 = 1 2 (2s + 1)ρ 2 1 + β 2 (2s + 1)η 2 1 + 2αβ(2s + 1)ρ 1 η 1 , and the solution with the single-well residual mixing (Equation (23)): give, for α = 1:Ē This proves that the two phases merge at the phase-separation line α = 1. On the other hand, for β = αs/(1 + s) (this curve is the (2s + 1)-well version of the boundary separating phases 2 and 3 in Figure 2 for L = 3 ↔ s = 1) one finds: showing how the current solution characterized by L = 2s + 1 at the lower border β = αs/(1 + s) reproduces the full-demixing solution one expects to find in the phase diagram for a sufficiently large α (this corresponds to the transition from phase 2 to phase 3 in Figure 2). With even L, the weak-demixing solution (19) (associated to phase 4 in Figure 3) of the case L = 4 can be easily generalized to case L = 2s. One determines: (defined in 1 < α ≤ 1/β, β ≤ 1) which can be also seen as the s-fold replica of solution (12) for the elementary case L = 2 (well occupied by the majority component and the remaining well with residual mixing). For β → 1/α, one recovers the full-demixing state with ρ 2p = 0 = η 2p−1 while for α → 1, the energy of this solution can be shown to match the energyĒ 0 of the full-mixing uniform state. The problem is that this solution (as well as that described by Equation (23) for odd L) now represents one among many other (unknown) ground-state solutions that must be identified and organized in order to reconstruct (as shown for the case L = 4) the map of adjacent phases forming the phase diagram. This task results to be extremely difficult.
The degree of complexity of the phase diagram for arbitrary L can be understood by considering the full-demixing phases in the large-α regime. As suggested by the case L = 4 (see phases 6 and 7 in Figure 3), these phases are expected to correspond to a number, increasing with L, of different realizations of the full-demixing state. In this case, the ring is essentially formed by p wells occupied by one component and L − p wells occupied by the remaining one (cases p = L and p = 0 are of course excluded). Then, one has L − 2 realizations of demixing: with ρ i = 0 for i ∈ [p + 1, L] and η k = 0 for k ∈ [1, p]. This configuration solves Equation (9). Interestingly, the value of p can be shown to depend on the value of β by imposing the minimum-energy condition. The energy of the previous configuration: entails that dẼ/dp = 0 for p = L/(1 + β). The number of wells associated to the two components is thus linked to the asymmetry parameter β. For U a = U b and interpreting p as an almost continuous parameter, one has β = N b /N a with: determining the numbers of wells associated to the two components in the ring lattice. The solutions emerging from the previous discussion for the two cases with even or odd L (and β close enough to β = 1) and the large-α full-demixing states (25) represent some of the few cases in which one can make a simple conjecture about the structure of the phase diagram in the limit of large L, in addition to 1) the uniform solution (and the relevant phase) defined for α ≤ 1 and for any L, and 2) the phase represented by solutions (12), (15) and (18) for L = 2, 3, 4, respectively, in the opposite region with sufficiently small β, whose general form for an arbitrary L, can be easily shown to be: For an arbitrary L, the solutions and the complex architecture of the corresponding phases occurring at intermediate values of β and α cannot be predicted in a simple way. Intuitively, one expects to find an increasing number of partial-demixing phases occupying the intermediate sector of the phase diagram which connects the full-mixing uniform phase (α < 1) and the full-demixing phases characterized by Equation (26) for a sufficiently large α.

Conclusions
We reformulated the BH-like model describing the binary mixture distributed in a four-well ring within the coherent-state variational picture in which a set of local order parameters describes the microscopic state of the system, and more specifically, the local populations of the two components. This well-known semiclassical approach allows one to derive the ground-state equations and thus reproduce the mechanism governing the formation of the new phases starting from the weak-interaction regime in which the system is fully mixed and delocalized. As discussed in Section 1, the realization of L-well rings, and then the study of their ground-state phases, is within the reach of current trapping techniques [18][19][20]. Various experimental aspects concerning a mixture in a three-well potential (realized by means of tweezer-based trapping techniques [42]) have been discussed in [31] revealing how the mixture 23 Na + 39 K constitutes an almost ideal case which permits a smooth tuning of the interspecies and intraspecies scattering lengths and a considerable reduction in three-body losses. Furthermore, one can identify the ground-state phases of the system through the detection of the boson numbers at each well, for each component, effected by means of absorption-imaging techniques [43,44]. The experimental realization of these systems, in addition to the possibility of validating theoretical predictions on the separation mechanism, will enable the testing of the effectiveness of the BH-model picture in which potential wells, despite their small but finite extension, are approximated to point-like structures. The numerical study of the mixture in the three-well potential within the mean-field picture [31] has confirmed the validity of the BH-picture predictions.
In this paper, we derived the phase diagram of the four-well system and discussed its critical properties. We compared it with those of the two-well and three-well systems, and tried to highlight both possible common features that do not depend on the lattice size L and properties that instead reflect the influence of the degree of fragmentation in the transition from the full-mixing to the full-demixing regime (complete space separation). Our analysis suggests that, in addition to the full-mixing uniform solution, various other phases are expected to survive for a large L. This is the case, for example, of the SWM solutions (12), (15) and (18) (phases 2, 3 and 4 in Figures 1-3, respectively) which correspond to state (27) for a generic L. Likewise, solutions (16) (phase 2 in Figure 2) can be associated, for arbitrary L = 2s + 1 and β 1, to solution (23) with fully demixed components, except for a single well with residual mixing, while solution (24) generalizes, for arbitrary L = 2s, solution (19) (phase 4 in Figure 3), showing s wells occupied by the majority component and s wells with residual mixing.
Solution (25), which describes full-demixing states for arbitrary L and sufficiently large α, represents the extension of solutions (20) and (21) (phases 6 and 7, respectively, in Figure 3), and has revealed that full demixing can be realized in L − 2 different ways related to the asymmetry parameter β. This number, in turn, highlights how the partial-demixing region (featuring intermediate values of β ∈ [0, 1] and α weak enough but such that α > 1), must be equipped with a number of ground-state solutions sufficiently large in order to reproduce, for large α, the transitions to the extended set of full-demixing solutions. The architecture of this region is then expected to be really complex and to feature manifold ways to reach the full-demixing regime whose recognition is far from being trivial (see, e.g., the four-well phase diagrams for a preliminary example). This richness strongly resembles the scenario of ground-state configurations of quantum emulsions [15,45] and suggests the observation of quantum-granularity effects triggered by the interplay of strong effective interactions α, weak hopping amplitudes and fractional values of β [46]. These aspects will be further explored in a separate paper.
Author Contributions: A.C. and A.R. performed analytical and numerical computations. All the authors equally contributed in reviewing the results and in the final discussion. V.P. wrote the article and supervised the work. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

BH
Bose-Hubbard CS Coherent state SWM Single well mixing