The phase-separation mechanism of a binary mixture in a ring trimer

We show that, depending on the ratio between the inter- and the intra-species interactions, a binary mixture trapped in a three-well potential with periodic boundary conditions exhibits three macroscopic ground-state configurations which differ in the degree of mixing. Accordingly, the corresponding quantum states feature either delocalization or a Schrödinger cat-like structure. The two-step phase separation occurring in the system, which is smoothed by the activation of tunnelling processes, is confirmed by the analysis of the energy spectrum that collapses and rearranges at the two critical points. In such points, we show that also Entanglement Entropy, a quantity borrowed from quantum-information theory, features singularities, thus demonstrating its ability to witness the double mixining-demixing phase transition. The developed analysis, which is of interest to both the experimental and theoretical communities, opens the door to the study of the demixing mechanism in complex lattice geometries.

where j = 4 ≡ 1 due to the periodic boundary conditions, T a and T b are the hopping amplitudes, U a and U b the intra-species repulsive interactions, and W represents the inter-species repulsion. We shall consider a mixture with twin condensates where U a = U b = U and T a = T b = T. Slight deviations U a ≠ U b and T a ≠ T b from this ideal case, as one expects to observe in real systems, can be shown to not alter in a significant way our results.

Results
The study of the equilibrium configurations of a bosonic binary mixture in a three-well potential (trimer), reveals how, according to the ratio between the interspecies and the intraspecies repulsion, W/U, three different families of qualitatively different ground states emerge in the zero-temperature regime. The main effect analysed in this paper, the two-step process where different mixing-demixing phase transitions are triggered by a change in W/U, is (a) investigated by means of a semiclassical approach to determine the ground-state configurations, (b) shown to be related to the collapse and rearrangement of the energy spectrum, (c) characterized in terms of quantum-correlation properties (entanglement entropy) between two parts of the systems.

Demixing transitions. It is advantageous to investigate spatial phase separation by reformulating
Hamiltonian (1) within the continuous variable picture (CVP) described in the section Methods. This allows one to replace Ĥ with an effective Hamiltonian operator typically providing a satisfactory representation of low energy states. In particular, the CVP approach allows one to reduce the search for the ground state of Hamiltonian (1) to the one for the minimum of the following effective potential: where variables x i and y i (i = 1, 2, 3) are normalized boson numbers, entailing that x i , y i ∈ [0, 1] and that z 1 + z 2 + z 3 = 1, (z = x, y) due to particle number conservation and T a = T b =:T, U a = U b =:U. Notice that only many-body states with the same fixed number N = M of bosons for the two species are considered. The configuration x y ( , ) → → which minimizes the effective potential V * is determined, at first, when the tunnelling T is suppressed. In this situation, in fact, it is possible to carry out a fully analytic study (see Methods section), capable of illustrating the physics of the problem in a particularly simple and effective way. The results, sketched in Fig. 1, can be summarized as follows: (i). For small W/U values (more specifically, if W/U ∈ [0, 1)), the two species are delocalized and uniformly distributed in the trimer, thus entailing a perfect mixing. (ii). (For moderate W/U values (namely, if W/U ∈ (1, 2)), the demixing occurs in two out of the three wells, while, in the third well, the species are still mixed. Of course, due to the symmetry of the trimer system, no well is favoured compared to the others, and, as a consequence, there are three configurations → → x y ( , ) which minimize the effective potential V * . Among such configurations, which are equal up to a cyclic permutation of the site labels, the one that we have illustrated in Fig. 1 is where each |ψ j 〉 corresponds to a different macroscopic configuration. Our analysis, which is semiclassical, breaks the symmetry of the problem and considers just one side of the cat, the one defined above and sketched in Fig. 1. It is interesting to observe that, crossing the critical point W/U = 1, the populations of site 1 and site 2 exhibit discontinuities, i.e. jumps from the uniform configuration x 1 = y 1 = x 2 = y 2 = 1/3 to the characteristic values 2/3 and 0. Conversely, functions x 3 (W/U) and y 3 (W/U) are continuous in W/U = 1. Increasing the ratio W/U from 1 to 2, the third site (where the two bosonic species coexist) gradually looses bosons, while the first and the second site (each one hosting just one species) become more and more populated. For W/U → 2 − , the first well hosts 3/4 of species-b bosons, while the second well hosts 3/4 of species-a bosons. The remaining parts coexist in the third well in equal measure. (iii). For big W/U values (i.e., for W U / 2 > ), the two species completely demix. One of them conglomerates in one site, while the other one spreads on the remaining two sites. This situation corresponds to 6 possible scenarios, i.e. to 6 configurations x y ( , ) → → minimizing the effective potential V * . In fact, if configuration ) is a minimum for V * , also the two configurations obtained by cyclic permutations of the site labels are minima for V * . Moreover, as the species we are considering feature the same hopping amplitude T and the same on-site repulsion U, the configurations where species-a conglomerates must not be favoured to the configurations where it is species-b the one that conglomerates. In other words, the energy of the system remains equal after the exchanges x i ↔ y i (i = 1, 2, 3). The presence of 6 minimum-energy configurations in the classical analysis corresponds to a quantum ground state Ψ 0 which is a six-faced Schrödinger cat, i.e. a state of the type where each |ψ j 〉 is a state corresponding to a certain macroscopic configuration. In Fig. 1 we have broken such symmetry and plotted just one of the 6 possible configurations. As illustrated in Fig. 2, for non-zero hopping amplitudes T, the critical points move rightward (i.e., the transitions occur at bigger values of W/U). It is possible to show with a fully analytic computation based on the request that Hessian matrix (11) must be positive definite (see Methods section), that the first mixing-demixing phase transition onsets at a value which is greater than 1, mirroring the fact that tunnelling T favors the mixed phase and delays the occurrence of the first demixing. On the other hand, one can notice that the effect of T is smaller and smaller as the number of bosons, N, increases and, in the thermodynamic limit, it becomes negligible. In such limit, the scenario is again the one depicted in Fig. 1.
Collapse and rearrangement of the energy spectrum. The dramatic change in the ground state qualitative structure occurring at the two critical values of W/U is mirrored by a deep rearrangement of the energy spectrum across these two transition points. As shown in Fig. 3, the energy levels (obtained by means of an exact diagonalization of Hamiltonian (1)) tend to collapse in the neighborhood of the critical values, meaning that the (Notice that, wherever two lines overlapped, they have been slightly displaced for purely graphical reasons). One can recognize the presence of three different quantum phases: (i) for W/U ∈ [0, 1) the two condensates are fully mixed and completely delocalized; (ii) for W/U ∈ (1,2) they are separated in two out of the three wells while they coexist in the third one; (iii) for W/U > 2 the demixing is complete, meaning that one species conglomerates in a well, while the other species distributes in the remaining two wells. It is clear that W/U = 1, 2 represent critical values where two different mixing-demixing phase transitions occur. Data have been obtained by means of a fully analytic study of the global minimum of effective potential (7) (see Methods section).   Entanglement entropy as a critical indicator. As is well known, the entanglement entropy (EE) quantifies the quantum correlation between two parts of a physical system through the Von Neumann entropy of a suitably defined sub-system 36,37 . In this work, the physical system is of course the binary mixture, the two parts correspond to the two different atomic species and, therefore, the EE describes the quantum correlation between them (see Methods section for a detailed mathematical definition). As shown in Fig. 4, this indicator strongly depends on the specific phase of the system which, in turn, is determined by the ratio W/U. In the completely mixed phase, i.e. for W/U smaller than the critical value (5), the EE steadily increases as W/U increases, despite the fact that semiclassical boson populations remain constant (uniformly distributed on the three wells). Notice that EE → 0 for W/U → 0 because, in this limit, the species are decoupled. At the two critical values of W/U, the EE features peaks, which not only indicate that the two mixtures are strongly entangled, but also allow to clearly distinguish the three different phases. Such peaks are sharper if T/U is smaller, that means that the effect of the hopping amplitude is to smooth the mixing-demixing transition. For W/U greater than the second critical value, the demixing gets more and more complete and the EE approaches the limiting value log (6) 2 . According to standard notions of Quantum Information Theory and to applications thereof in the field of ultracold bosons 35,36 , the emergence of this value can be appreciated by recalling that there are 6 possible ways of realizing a completely demixed semiclassical configuration which, in turn, equally contribute to the formation of a unique quantum ground state of the type (4).

Discussion
We have investigated the two-step phase separation of a binary mixture in a trimer, highlighting the presence of an intermediate phase, which is neither completely mixed nor completely demixed. This element originates from the fact that the number of sites is odd and it is absent in the two-well system 28,29 , where there is a one-step transition from a fully mixed to a fully demixed configuration.
The problem of determining the ground state has been reduced, by means of the CVP approach (see Methods section), to the search for the minimum-energy configuration of effective potential (2). The latter features an effective energy landscape whose minimum points are formed or disappear depending on the ratio W/U. We start analyzing the situation where the tunnelling T is suppressed. This case is of particular importance not only because it is associated with fully analytic results, but also because it captures the system behavior in the large N limit. As explained in the Methods section, the search for the global minimum of function (2) is not trivial, as one has to take into account that the global minimum is not necessarily a stationary point. A careful analysis, based on a systematic exploration of the domain boundary leads to the identification of three different phases, which qualitatively differ in the degree of mixing and spatial localization (see Fig. 1).
The activation of the hopping amplitude T makes the model more realistic but does not distort the scenario, which is still marked by two critical points. From a technical point of view, the search for the minimum energy configuration is simpler because T constitutes a regularizing term which prevents the points of minimum from falling on the boundary. The exploration of the domain boundary is therefore no longer necessary and it is enough to look just for local points of minimum in the interior of the domain. This conceptual simplification comes with the impossibility of analytically solving the resulting system of (10). With reference to Fig. 2, where numerical results are plotted, one can observe that bigger W/U values are needed to trigger the demixing transitions which, in turn, are less abrupt. To test their reliability, we have compared CVP's predictions to the results obtained by means of exact numerical diagonalization. Figure 5 shows that the two mixing-demixing phase transitions occur at the same critical values of W/U and, more generally, that the semiclassical approach and the purely quantum scenario are in good agreement. Notice that we chose not to plot the expectation values of boson populations in the ground state (e.g. 〈Ψ 0 |N j |Ψ 0 〉) because one would obtain the average over the different macroscopic configurations constituting the cat-states and, therefore, the constant value N/3, irrespective of the ratio W/U. Instead, we have determined, sweeping the parameter W/U, the most probable configuration i.e. the Fock state associated to the biggest (square modulus of the) coefficient in the decomposition of the ground state. It is interesting to notice that function (2) features a unique minimum in the uniform and fully mixed phase, three isoenergetic points of minimum in the intermediate phase and six isoenergetic points of minimum in the fully demixed phase. These degeneracies in the semi-classical energy landscape are resolved, quantum mechanically, with the formation of cat-states of the type (3) and (4) respectively.
The phase diagram of the binary mixture on a three-well potential, sketched in Fig. 6, has been determined numerically solving systems (10) and (12), and checking that the obtained stationary point actually is a minimum point for effective potential (2). Purple solid lines represent critical condition (5) derived analytically while no analytic expression can be found to define the orange dashed lines due to the complex nonlinear character of (10). These lines have been plotted on the basis of numerical data (see Methods section for details).
The double critical behaviour exhibited by the system upon a variation of W/U can be evidenced also from the energy-spectrum standpoint 29 . As shown in Fig. 3, whose data have been obtained by means of an exact (numerical) diagonalization of Hamiltonian (1), the energy levels tend to collapse in the neighborhood of the two critical points and have different arrangements in the three phases. One can recognize, for example, families of energy levels which are clearly distinct in the first two phases while tend to overlap in the third one. Approaching the limiting case T → 0 (note that this can be shown to be equivalent to taking N → +∞), the collapses are more and more evident and the points where they take place tend to W/U = 1 and W/U = 2. The low-energy excitations spectrum of a binary mixture in the uniform and completely mixed phase was determined for a generic L−site ring by means of a group-theoretic approach 27 combined with the Bogoliubov scheme. Choosing L = 3 sites, the condition corresponding to the collapse of a Bogoliubov frequency exactly reproduces (5), which gives the onset of the first mixing-demixing phase transition. One can thus appreciate the fact that the range of validity of the quasi-particles Bogoliubov spectrum exactly corresponds to the region of parameters space where the system is in the uniform and fully mixed phase. This analysis is confirmed by the behaviour of the ground state energy  Fig. 7, show, for T/U → 0 (or, equivalently, for N → +∞), the emergence of two critical points at W/U = 1 and W/U = 2. In such points, the ground state energy features discontinuities in the first derivative, a circumstance which suggests the onset of a quantum phase transition. The  systematic study of this double phase transition will be developed elsewhere in the same spirit of the localization-delocalization transition in the Bose-Hubbard trimer 38 .
Quantities traditionally belonging to quantum information theory have been used to highlight the occurrence of quantum phase transitions [39][40][41][42][43][44] . Among the others, Entanglement Entropy has proved to be particularly effective in detecting the quantum criticality of the localization transition in a bosonic Josephson junction 45 . This indicator, which measures the degree of entanglement between two parts of a given system through the Von Neumann entropy of a reduced subsystem, has been shown to be sensible also to the mixing-demixing phase transition occurring in a two-species bosonic mixture loaded in a dimer 28 and to be robust with respect to the change of partition scheme 36 . As explained in the Methods section, one starts from the density matrix associated to the ground state, i.e. from ρ ψ ψ = 0 0 0 . Then, a partition scheme aimed at identifying two sub-systems is chosen. In our case one sub-system corresponds to species-a bosons while the other one corresponds to species-b bosons. The density matrix of a reduced subsystem, Tr ( ) , is obtained by tracing out the degrees of freedom of the other one (see formula (13)). One thus remains with the reduced density matrix of a mixed state and the Von Neumann entropy there of (14) corresponds to the entanglement between the two sub-systems, namely between the two bosonic species.
By means of the exact (numerical) diagonalization of Hamiltonian (1), we have computed the EE as a function of W/U for different values of T/U. The results, illustrated in Fig. 4 and discussed in the previous section, not only demonstrate the ability of EE to capture the double critical behaviour of the system upon a variation of W/U, but also witness the formation of 6-faced cat-like states of the type (4) for large W/U values. Notice, in fact, that, for W/U → +∞, the EE tends to the limiting value log (6) 2 , where the number 6 is the number of different semiclassical configurations x y ( , ) → → which minimize effective potential (2) when  W U / 3. In this work, we have evidenced a phase separation mechanism considerably more complex than the one occurring in a two-well system 28,29,36 . The presence of an intermediate phase, which stands between the uniformly mixed configuration and the completely demixed phase, opens the door to the study of the phase separation mechanism in rings having an even or an odd number of sites and, more generally, in more complex lattice topologies. An illustrative example is shown in Fig. 8, where we consider a 4-well ring lattice. The CVP approach leads to a generalization of effective potential (2) which, in the case of vanishing tunnelling and at the transition point W/U = 1, can be shown to be minimized by any macroscopic configuration of the type y j = 1/2 − x j , with j = 1, 2, 3, 4 (four of them are shown in Fig. 8). As soon as > W U / 1, this degeneracy is broken and one macroscopic configuration gets energetically more favorable. The determination (a) of the ground state configuration when ). This scheme leads to a new effective Hamiltonian written in terms of coordinates x j , y j and of their generalized conjugate momenta 38 . Accordingly, in the CVP framework, the relevant eigenvalue problem ˆ= H E EE (after expanding quanity H Ê up to the second order) can be recast as

E E
The application of this rather versatile approximation scheme to Hamiltonian (1) (which describes a binary mixture confined in a three-well optical lattice with periodic boundary conditions) leads to an effective eigenvalue problem of the type (6) where Minimum-energy configuration for T = 0: geometric approach on a polytope-like domain. If the tunnelling is suppressed, effective potential (2) simplifies as follows . The domain  is a 4-polytope, consisting in the direct product of two triangular regions in 2  , as shown in Fig. 9 (left panel). Computing the gradient with respect to the four independent variables, one finds that the only stationary point is the one associated to the uniform configuration, i.e x 1 = x 2 = y 1 = y 2 = 1/3, corresponding to The eigenvalues of the associated Hessian matrix show that this stationary point is a minimum provided that W < U. Notice that, in general, computing the stationary points does not necessarily give the global minimum. The latter, in fact, can live on the boundary of the domain , in a point where the four-dimensional gradient is not well defined. So, one has to compare the value of ν at the local minimum, i.e. N 2 (U + W)/3, with the values of ν at the boundary of .
The exhaustive exploration of the domain boundary is a fortiori necessary when W > U i.e. when the stationary point x 1 = x 2 = y 1 = y 2 = 1/3 is no longer a minimum, a circumstance which entails that the global minimum lives on the domain boundary. The complexity of our problem is due to the fact that the domain of ν is four-dimensional and so its boundary is the union of six three-dimensional objects of the type sketched in Fig. 9 (right panel). Each of them corresponds to a Cartesian product × ′  I j x and ×  I j y (where I j and ′ I j , with j = 1, 2, 3 are the edges of triangles x  and y  , respectively) and constitutes the domain of a 3-variable function obtained by introducing an additional constraint in formula (7). In the same spirit, the global minimum of these constrained functions must be searched not only in the interior of their domain (imposing the stationarity of a three-dimensional gradient), but also on their boundaries, by means of an exhaustive exploration thereof. One therefore iterates this process looking for stationary points, at first inside the volume, then on the faces, on the edges (employing, at each step, a lower-dimensional gradient) and, eventually, evaluates the function at the vertices of this 3D region. The resulting set of candidates for the global minimum is such that each local minimum is linked to an existence condition (typically, an inequality in the one-dimensional space W/U) on the sub-domain where it was found. The global minimum, in each interval of the space W/U, is determined by comparing all the possible candidates. In the following we sketch the application of this rather general scheme to our problem.
We start fixing x 1 = 0. Notice that the other possible ways of fixing the first variable, i.e. x 2 = 0, x 2 = −x 1 + 1, y 1 = 0, y 2 = 0 and y 2 = −y 1 + 1 would lead to minimum-energy configurations which are equivalent up to cyclic permutations of the indexes and/or species labels swapping. The resulting constrained function (1 ) is defined on the prism represented in the right panel of Fig. 9. The stationary point of this 3-variables function, computed imposing ν ∇ =(0, 0, 0) x y y , , is a local minimum of ν 3 if W < U but it must be discarded because the corresponding value of ν is greater than (8). Let us focus on the 5 faces of this prism. Each face is fixed by introducing an additional constraint. As an example, face ABC is the domain of function Computing the gradient, y y , 1 2 ∇ one finds that there is one stationary point, which is ( ) . This is a minimum (both eigenvalues of the Hessian matrix bigger than 0) in the domain of interest iff W < U/2. Nevertheless, the corresponding value of ν is smaller than the uniform-configuration value (8) just for W > U. But, in this range, this stationary point is no longer a minimum, so it must be discarded. The analysis of the 5 faces, carried out according to this scheme, leads to the conclusion that, for 1 < W/U < 2, there are 2 minimum-energy configurations Notice that these two configurations are equivalent up to a permutation of site indexes. The first one has been plotted in Fig. 1.
The 9 edges of the solid of Fig. 9 (right panel) are identified by fixing an additional constraint. Their analysis shows that local iso-energetic points of minimum are present in segments CB DF and CF. Such minimum energy configurations are equal up to permutation of site indexes and/or species labels swapping. For example, segment CB is the domain of function 2 of ν, is less than the one found in faces ACFD and CBEF provided that W > 2U. Eventually, the analysis of the 6 vertices does not give any further minimum-energy configuration.
The search for the minimum-energy configuration for T ≠ 0. The minimum-energy configuration is the configuration → → = x y x x x y y y ( , ) ( , , , , , ) 1 2 3 1 2 3 which minimizes effective potential (2). The presence of a non-vanishing T constitutes a regularizing term which makes the exploration of the domain boundary no longer necessary. This fact, verified numerically, can be understood in physical terms because the presence of T ≠ 0 makes configurations where one or more populations are exactly zero impossible. To compute the SCientifiC RePoRTs | (2018) 8:10242 | DOI:10.1038/s41598-018-28573-w minimum-energy configuration, one therefore imposes stationary conditions for effective potential V * in which x 3 = 1 − x 1 − x 2 and y 3 = 1 − y 1 − y 2 take into account the boson number conservation. The resulting system of equations is , , (0, 0, 0, 0) (9) 1 2 1 2 and the explicit form thereof is NT  In general, it is not possible to find all the possible solutions of this system in a closed form and one needs to resort to numerical methods. Boson populations shown in Fig. 2 not only fulfill system (10), but have been checked (by evaluating the eigenvalues of each associated Hessian matrix), to be minimum points of the effective potential (2). As already mentioned, Fig. 2 illustrates a particular family of solutions, meaning that those configurations which are obtained by means of cyclic permutations of site indexes and/or by species-labels swapping, are still solutions of system (10) and therefore points of minimum of (2).
It is worth noticing that a particularly simple and significant solution of system (10) is the one which corresponds to the uniform configuration x 1 = x 2 = y 1 = y 2 = 1/3. This stationary point is a minimum of (2) provided that the associated Hessian matrix is definite positive, a condition which is verified iff W/U < 1 + (9T)/(2UN). This argument proves the correctness of critical value (5), presented in the main text and corresponding to purple solid lines in Fig. 6. An effective way to reduce the complexity of system (9) (and possibly simplify its numerical solution) comes from the symmetry which marks the boson populations when the system is in phase (ii) (see central regions of Fig. 2 panels). In such phase, i.e. for moderate W/U values, one can notice, in fact, that y 2 = x 1 and y 1 = x 2 , two constraints that allow to rewrite system (9) as    The orange dashed lines in Fig. 6, which constitute the border between the intermediate and the fully demixed phase, have no analytical expression and correspond to those solutions of system (12) which are no longer minimum points for effective potential (2), i.e. to stationary points for which at least one Hessian-matrix eigenvalue