Phase separation can be stronger than chaos

We investigate several dynamical regimes characterizing a bosonic binary mixture loaded in a ring trimer, with particular reference to the persistence of demixing. The degree of phase separation is evaluated by means of the"Entropy of mixing", an indicator borrowed from Statistical Thermodynamics. Three classes of demixed stationary configurations are identified and their energetic and linear stability carefully analyzed. An extended set of trajectories originating in the vicinity of fixed points are explicitly simulated and chaos is shown to arise according to three different mechanisms. In many dynamical regimes, we show that chaos is not able to disrupt the order imposed by phase separation, i.e. boson populations, despite evolving in a chaotic fashion, do not mix. This circumstance can be explained either with energetic considerations or in terms of dynamical restrictions.


Introduction
The demixing of the condensed species constituting a binary bosonic mixture, i.e. their localization in different spatial regions, is a process that can be triggered by the presence of strong inter-species repulsive interactions. The underlying mechanism of phase separation (also called species separation) has been thoroughly investigated within the field of ultracold bosons by means of the mean-field representation of condensate dynamics [1,2,3]. The onset of the demixing transition has been shown to depend on a number of factors, including the shape of the trapping potentials, the number of bosons in each atomic species and the interaction parameters. At the same time, considerable attention has been devoted to the dynamical-stability analysis and to the appearance of excited states in regimes close to the transition point. [4,5,6,7].
Engineering of experimental setups capable of trapping boson mixtures [8,9,10] in optical lattices [11,12,13] has stimulated the interest of the theoretical community for the understanding of demixing effect in the presence of spatial fragmentation. In this context, beyond phase separation [14,15,16], a rich variety of phenomena have been highlighted including (but not limited to) the formation of boson currents in ring lattices, [17], the appearance of novel phases exhibiting magnetic-like properties [18,19], quantum emulsions [20], the emergence of polaron excitations [21,22], the entanglement between the species [23,24], the collision of a condensate with impurities [25,26,27], the formation of immiscible solitons [28] and the modulation instability in the phase separation [29]. Systems consisting of two bosonic species and optical lattices with ring geometry are within the reach of current experimental setups. In particular, heteronuclear and homonuclear mixtures have been realized in [8,30] and [9,31] respectively, while the ring-lattice geometry has been designed and employed by [32,33]. Moreover, population oscillations can be monitored with the same techniques used to reveal the self-trapping effect of a single condensate in a two-well system [34,35].
Recently, the phase separation mechanism of a binary bosonic mixture has been investigated in the two simplest but non-trivial lattice geometries: the double well [24,36,37] and the ring trimer [38]. It has been evidenced that in the former case only one kind of transition occurs, while in the latter case the demixing develops in two steps, i.e. an intermediate neither fully mixed nor completely demixed phase exists. Moreover, it has been shown that in the parameters' space of the model, the critical point where phase separation occurs is characterized by the collapse and rearrangement of the energy levels and by singularities in the entanglement entropy between the two species [24,38].
Provided that the number of bosons is large, one can consider the semiclassical counterpart of the quantum system and investigate the dynamics generated by the resulting set of discrete nonlinear Schrödinger equations [39,40,41]. The latter, which constitute the discrete analogue of the Gross-Pitaevskii equation, feature an extraordinary rich scenario of dynamical regimes. Already in the relatively simple case of a single condensed species loaded in a ring trimer, one can identify both stable and unstable regimes including vortices, dimerlike states and chaotic oscillations [42,43,44]. In [45], the non-integrable character of low dimensional circuits has been evidenced and chaos has been shown to support the persistence of superfluidity.
In this work we investigate, by means of a semiclassical approach, the dynamics of a bosonic binary mixture loaded in a ring trimer, emphasizing its relation with the entropy of mixing and the persistence of spatial phase separation. After identifying three classes of stationary configurations featuring an high degree of demixing and after developing the energetic-and the linear-stability analysis, we simulate the dynamics of thousands of trajectories starting in the vicinity of fixed points (FPs). These simulations (i.e. the numerical solutions of motion equations (6)) not only allow one to compute the first Lyapunov exponent, an indicator which allows to distinguish between regular and chaotic trajectories, but also give the possibility to monitor the degree of mixing of the two condensed species. Contrary to expectations, we show that there are several dynamical regimes where chaos, despite present, is not able to disrupt the order imposed by spatial phase separation. This circumstance will be clarified both with dynamical considerations and in terms of energy conservation.
The outline of the manuscript is as follows: in section 2 we present the model and its semiclassical counterpart. In section 3, we identify three notable classes of stationary configurations featuring spatial phase separation. Section 4 is devoted to the analysis of their energetic and linear stability. In section 5, we perform our numerical simulations and we compute the first Lyapunov exponent. An indicator to quantify the degree of mixing of the atomic species is presented in section 6. In section 7, we discuss some meaningful dynamical regimes, putting particular emphasis on those ones where chaos and persistent demixing coexist. Eventually, section 8 is devoted to concluding remarks.

A binary mixture in a ring trimer
The second-quantized Hamiltonian describing a bosonic mixture of two atomic species in a three-well potential (with periodic boundary conditions) iŝ where j = 4 ≡ 1 due to the ring geometry. This is a typical Bose-Hubbard Hamiltonian, where T a and T b are the tunnelling amplitudes, U a and U b represent intra-species repulsive interactions and W corresponds to the inter-species repulsion. Creation and destruction operators satisfy usual bosonic commutators, namely j ] = 0.n j =â † jâ j andm j =b † jb j are number operators and their sums,N = 3 j=1n j andM = 3 j=1m j respectively, constitute two independent conserved quantities, being [N ,Ĥ] = [M ,Ĥ] = 0. Provided that the number of bosons is sufficiently high [46], it is possible to replace field operators in Hamiltonian (1) with local order parameters, [39,47]. Such substitutions, which explicitly read allow one to cast the quantum dynamics generated by Hamiltonian (1) in a classical form, that is It is convenient to express local order parameters in terms of number of bosons and local phase [45,46], i.e. a j = √ n j e iφ j and b j = √ m j e iψ j . One thus obtain the following classical Hamiltonian which, in turn, after setting = 1, entails the following motion equationṡ

Notable demixed stationary configurations
The exhaustive study of all possible stationary configurations (i.e. configurations having a trivial time dependence) which the system admits goes beyond the scope of this work, as the scenario is extraordinarily branched and rich. Already with a single condensed species confined in a ring trimer, several classes of stationary states (e.g. vortex, π and dimerlike states) have been evidenced [42].
In this work, we put the focus onto the miscibility properties of the two condensed species, and we start our analysis from some notable stationary configurations which feature phase separation. According to the theory of discrete nonlinear Schrödinger equations [48], substitutions φ j → Φ j + λ a t and ψ j → Ψ j + λ b t, where λ a and λ b represent collective angular frequencies of condensates' phases, constitute a preliminary step in the search for stationary configurations (the presence of two independent collective frequencies λ a and λ b , follows from the density-density form of the interspecies interaction of Hamiltonian (1)). These substitutions, in fact, allow one to recast Hamilton equations (3) and (4) into the following dynamical system with j = 1, 2, 3. Looking for FPs of the latter (which therefore correspond to stationary solutions of equations (3) and (4)), together with the two constraints 3 j=1 n j = N and 3 j=1 m j = M , one finds three classes of configurations which, in the limit T a , T b → 0, feature perfect demixing. They are schematically illustrated in figure 1 (upper row) and described below: (i) Dimer -Soliton: Condensate A is equally subdivided in two wells, the phases therein being the same, while the third well contains all the condensate B. (ii) Single Depleted Well (SDW) -Soliton: Condensate A is equally subdivided in two wells but, contrary to the previous case, the relative phase Φ 3 − Φ 1 between such wells is π. The third well contains all the condensate B.
(iii) Soliton -Soliton: One well contains all the condensate A while an other well contains all the condensate B.
Upon activation of hopping amplitudes T a and T b , FPs slightly deviate from the aforementioned ones, as some bosons move from the macroscopically occupied wells to neighbouring ones. The new scenario of FPs, thus moderately blurred by the presence of non-zero tunnelling processes, is pictorially sketched in the lower row of figure 1 and fully discussed in Appendix A. x j , y j x j , y j x j , y j 1 1 1 2 2 2 3 3 3 x j , y j x j , y j x j , y j 1 1 1 2 2 2 3 3 3 Dimer -Soliton SDW -Soliton Soliton -Soliton Figure 1. The three families of stationary, demixed, configurations for zero (upper row) and non-zero (lower row) hopping amplitudes T a and T b . Lower row displays in an exaggerate but illustrative manner the deviations from the zero-tunnelling scenario. Numbers on the horizontal axis correspond to wells' labels, while the height of the histograms represents normalized populations x j = n j /N and y j = m j /M . Parallel (antiparallel) arrows stand for "in-phase" (antiphase) condensates. Wherever not explicitly defined, condensate phases assume different values according to different choices of model parameters (see Appendix A).
In the following, we shall assume that the two condensates feature the same dynamical parameters, i.e. U a = U b =: U , T a = T b =: T and N = M . Nevertheless, we note in advance that the presence of small deviation T a = T b , U a = U b and N = M , from the previous ideal conditions (deviations which could be present in a real experimental setup), have been proved (by means of numerical simulations) not to significantly affect the developed analysis and the dynamical scenario discussed in the following. Both the hopping amplitudes and the interaction strengths depend on the lattice constants and on the scattering length [49], which are tunable parameters of the experimental setups mentioned in the introduction. In particular, interaction strengths can be controlled by means of Feshbach resonances [50,51].

Stability of stationary demixed states
Energetic-stability and linear-stability analysis are standard but powerful tools to investigate the qualitative dynamical behaviour of the system in the vicinity of a FP [52,53]. In view of an experimental realization, the developed analysis plays an important role, since it is impossible to prepare the system in a state which exactly coincides to one of the aforementioned stationary configurations. Preliminary, it is convenient to introduce vector and to write dynamical system (5) in the compact forṁ where is the standard symplectic matrix, is the effective Hamiltonian and ∇H = (∂ Φ 1H . . . , ∂ m 3H ).

Energetic stability
An effective way to determine whether a FP z * is energetically stable or not is to study the signature of the relevant Hessian matrix According to Lagrange-Dirichlet Theorem, a FP z * is energetically stable if H( z * ) is positive or negative definite [53] (to be more precise, in the same spirit of [45], one has to exclude the pair of vanishing eigenvalues corresponding to the two conserved quantities or, equivalently, consider a 8 × 8 Hessian matrix obtained after explicitly introducing constraints Φ 1 = Ψ 1 = 0, n 1 = N − n 2 − n 3 and m 1 = M − m 2 − m 3 ). A point exhibiting energetic instability is therefore neither a local minimum nor a local maximum of the energy functionH [45]. With reference to the first row of figure 2, one can observe that no FP of the class "SDW -Soliton" is energetically stable, each of them being a multidimensional-saddle point for Hamiltonian functionH (see second panel). Interestingly, the energetically-stable region relevant to FPs of the class "Dimer -Soliton" (see first panel) exactly corresponds to one of the three kinds of ground states that were found and discussed in [38]. In fact, all FPs z * in such region (depicted in green) are indeed global minima of functionH. Eventually, observing the third panel, one can recognize the presence of an energetically-stable region for moderately low values of W/U (depicted in blue). In this region, FPs z * are local maxima of functionH. White regions in figure 2 correspond to those values of W/U and T /(U N ) for which FPs belonging to a given class do not exist. Such regions stand in between different sub-classes which differ in the relative phases Φ j − Φ j−1 between the wells. Each sub-class is a portion of parameters' space where stationary configurations share common features (e.g. the relative phases) and which are delimited by white regions. Fox example, the first class includes three sub-classes which, in turn, include points 1A, 1B and 1C respectively.

Linear stability
The linear stability (also called dynamical stability [45]) of a FP z * of motion equations (6), (i.e. a configuration such that˙ z * = 0) is determined by the eigenvalues of Jacobian matrix [52] More precisely, as dynamical system (6) is a Hamiltonian one, a FP z * is said to be linearly stable (or elliptic) if all eigenvalues of J i,j ( z * ) are purely imaginary; conversely, it is said to be linearly unstable if at least one (pair of) eigenvalues of matrix (9) has non-zero real part. In the second row of figure 2, obtained by sweeping model parameters W/U and T /(U N ), we have represented, for each of the three classes of stationary points characterized by demixing, the largest real part among the eigenvalues of matrix (9). One can notice wide regions of the parameters' space where FPs of the class "Dimer -Soliton" are linearly stable (represented in dark blue). Interestingly, while in the first sub-class (the one including point 1A), all FPs are linearly stable, in the remaining two sub-classes (respectively including points 1B and 1C) there are regions featuring linear stability and regions featuring linear instability. FPs of the class "SDW -Soliton" (see second panel), are mostly linearly unstable, excepts for a tiny triangular-like region existing only for W/U > 2 (notice that FPs obtained in the unphysical situation T = 0 are linearly stable too). As shown in the third panel, the vast majority of FPs belonging to the class "Soliton -Soliton" are linearly stable, except for those ones in a narrow band confining with the white region and those ones in a needle-like region present for W/U ≈ 2 and moderately high values of T /(U N ).

Scope of the energetic-and the linear-stability analysis
If a trajectory moves away from a FP, the energetic-and the linear-stability analysis thereof are of little use. For this reason, one should employ other indicators such as the first Lyapunov exponent, which is the gold standard to distinguish regular and chaotic trajectories. A further limitation affecting the energetic-and the linear-stability analysis comes from their local character. More specifically, also in view of an experimental realization, one should pay particular attention to the size of the FP's neighborhood where they are valid. Both aspects are discussed in section 5.
An important remark is in order concerning the traditional criterion to evaluate the linear stability of a FP (i.e. all eigenvalues λ j 's of matrix (9) must be purely imaginary). Interestingly, the latter fails when the characteristic frequencies ω j = I{λ j } satisfy a certain commensurability condition, essentially represented by a diophantine equation (see [53] for details, in particular for the procedure which is used to determine their sign). In this case, in fact, a so-called "elliptic" FP ceases to be the center of an elliptic island and turns unstable. More specifically, it has been proven by Moser [53] that, if the initial configuration z(t = 0) =: z 0 is sufficiently close to a linearly stable FP z * , solutions of the actual non linear system (6) almost always depart from those of the linearized one only extremely slowly, if at all. Nevertheless, this is true only if characteristic frequencies ω j 's, properly taken with a certain sign, do not satisfy the aforementioned commensurability condition. The frequency vectors ω which satisfy such condition constitute a dense set, although of measure zero, excepts in the positive (ω j > 0 ∀j) and negative (ω j < 0 ∀j) quadrants of space ω. These two quadrants exactly corresponds to the regions where energetic stability holds.

Regular and chaotic oscillations of boson populations
For each of the 102729 pairs of model parameters (W/U, T /(U N )), a starting point z 0 very close to the relevant FP z * is chosen in such a way that the relative difference between the vector components of z 0 and the corresponding ones of z * is from 2% to 5% thus emulating what could be achieved in a real experimental set up [34,35,54,55]. Then motion equations (6) are numerically solved ‡ for a series of consecutive time intervals and the first Lyapunov exponent is iteratively computed according to the standard scheme described in [56]. The comparison between the results (see third row of figure 2) and the previously discussed linear-stability analysis (see second row of figure 2) indeed shows that if a FP z * is linearly unstable, than a trajectory starting from a point z 0 close to it is chaotic, i.e. it is associated to a non-zero Lyapunov exponent. Note to the reader: in the following, labels 1A,...,3C are indistinctly used to indicate both a FP or a trajectory starting in a neighborhood thereof. Likewise, the vast majority of FPs featuring linear stability is such that the trajectory originating from a point z 0 close to it is regular. Actually, for a limited number of FPs this is not ‡ Computational resources provided by HPC@POLITO (http://www.hpc.polito.it) true. For example, with reference to the central column of 2, the tiny linearly stable region present in the second row has no counterpart in the third row, all the trajectories therein represented being chaotic. This circumstance can be interpreted in terms of size of the elliptic islands centered around an elliptic FP. As already evidenced in figure  7 of reference [42], such elliptic islands are very small when their center is a FP of the class "SDW -Soliton" and so the distance | z 0 − z * |, despite chosen to be small, is already greater than the islands' characteristic radius. Moreover, it is worth noticing the presence of curved lines featuring a large Lyapunov exponent (e.g., in the first panel of the third row of figure 2, the curve whose bounds are points (0.5, 0.075) and (2, 0) in parameters' space (W/U, T /(U N ))) which are expected to correspond to linearly stable FPs (see first panel of the second row). The seed of chaotic behavior is, in this case, the commensurability of characteristic frequencies ω j 's, which are the imaginary parts of the eigenvalues of Jacobian (9) (see [53] for details, in particular for the procedure which is used to determine their sign). In fact, one can verify that all FPs constituting the aforementioned curve are such that ω 1 = −2ω 2 .

How to quantify mixing and demixing of boson populations
Looking at the first row of figure 1, one can recognize that the three presented configurations feature perfect demixing, as the presence of a condensed species in a certain well always implies the complete absence of the other species. In the second row of the same figure, where the aforementioned ideal configurations are blurred by the activation of tunnelling processes, the two species, despite being still overall separated, feature a small degree of mixing. In fact, in whichever well a certain species is macroscopically present, the other one is nearly absent, yet non zero. In the following we present an indicator to quantify the degree of separation or, to be more precise, the degree of mixing. Entropy of mixing, generally denoted with S mix , is a standard indicator commonly used in Statistical Thermodynamics when investigating miscibility properties of chemical compounds [57,58,59] whose role, in the present work, is played by quantum gases. As the geometry behind the extended Bose-Hubbard model we are investigating is inherently discretized, one has to compute the entropy of mixing in each well and then evaluate the average over the three wells. One therefore obtains that In passing, we observe that the relative phases between the wells play no role in the computation of S mix , but are crucial for its time evolution. The fourth row of figure 2 shows the value of the entropy of mixing for all FPs belonging to the three classes discussed in section 3. As expected, we observe that S mix steadily increases with T /(U N ), since a bigger hopping amplitude blurs the fully demixed configurations depicted in the upper row of figure 1. Contrary to expectations, there are regions where S mix increases for increasing inter-species repulsion W/U (to be more specific, the region featuring 1 < W/U < 2 in the first panel, the region for 0 < W/U < 2 in the second panel and the region featuring 0 < W/U < 1 in the third panel). This circumstance looks counter intuitive, but it is easily explained by recalling that FPs are not necessarily minimum-energy configurations. As illustrated in the second row of figure 2 and discussed in section 4.1, in fact, a FP can also be a local maximum or a saddle point for effective Hamiltonian (7). In the third panel, moreover, we observe that there are values of W/U and T /(U N ) for which S mix tends to the maximum value log 2. This happens because FPs of the class "Soliton -Soliton" have been so blurred by the activation of the hopping amplitude T , that they have almost lost their identifying aspect and turned into a uniform configurations of the type z un (as shown in figure  A3, notice that the stationary configuration continuously varies with respect to model parameters W/U and T /(U N )).

Competition between phase separation and chaotic behaviour
In the same spirit of section 5, we have numerically solved motion equations (6) choosing, for each parameters' pair (W/U, T /(U N )), a starting point z 0 close to the corresponding FP z * . The choice has been made in such a way that the difference between the vector components of z * and those of z 0 is from 2% to 5%. The knowledge of the time evolution of boson populations n j (t) and m j (t) allows one to readily compute S mix (t) and so to monitor the mixing properties of the atomic species. If the initial configuration z 0 lies within a regular island centered around a linearly stable FP z * (and the characteristic frequencies thereof do not match Moser's commensurability condition [53]), the motion consists in small oscillations around z * . Therefore, the entropy of mixing S mix (t) features small oscillations around the constant values S mix ( z * ) which, in turn, is very low for the vast majority of FPs z * belonging to the three classes of notable stationary demixed configurations under considerations (recall the fourth row of figure 2). On the other hand, if z 0 lies outside regular islands, the motion is chaotic, as discussed in section 5. In these circumstances, one would expect the two condensed species to fully mix and thus to quickly loose memory of their initial, demixed character. We show that this is not always the case, as demixing and chaos can coexist indeed. To this purpose, for each simulated trajectory, we have recorded max t {S mix (t)}, where t ranges from t = 0s to t = 50s, a time interval whose width is three orders of magnitudes larger than the smallest characteristic period of populations' oscillations. Overall, 102729 trajectories have been simulated, each one starting from Energetic-stability analysis of FPs z * . Green and blue correspond to energetically stable regions, as they represent, respectively, minima and maxima of Hamiltonian (7). Yellow refers to energetically unstable regions, i.e. to saddle points of Hamiltonian (7). Second row: Linear-stability analysis of FPs z * . The color corresponds to max j {R{λ j }} where λ j 's are the eigenvalues of matrix (9). Dark blue represents linearly stable FPs. Third row: First Lyapunov exponent associated to trajectories starting close to FPs z * . Dark blue (yellow) corresponds to regular motions (highly chaotic trajectories). Fourth row: FPs' entropy of mixing, i.e. S mix ( z * ). Blue (red) represents a remarkable phase separation (an high degree of mixing). Fifth row: Maximum entropy of mixing, i.e. max 0<t<50 {S mix (t)} relevant to trajectories starting close to FPs z * . Blue corresponds to trajectories which feature a very small degree of mixing S mix throughout all the simulated dynamics. Conversely, red is associated to trajectories which, one or more times during the time evolution, feature complete mixing. All panels: White regions correspond to model parameters W/U and T /(U N ) for which no stationary solutions of the type defined in the column's title exist. In each panel, three points have been highlighted in order to facilitate the discussion. Model parameters N = 50 and U = 1 have been chosen. an initial condition z 0 which, in turn, is close to a FP z * . The result is shown in the fifth row of figure 2.
In the attempt to facilitate the comparison among the info provided by the linearstability analysis (second row of figure 2), by the computation of the first Lyapunov exponent (third row), by the evaluation of S mix ( z * ) (fourth row) and of max t {S mix (t)} (fifth row), we highlight and analyze some notable dynamical regimes and we explicitly illustrate them in figures 3-5. We group them according to the regularity of the motion and to the persistence of demixing during the dynamics. We remind that labels 1A,...,3C are indistinctly used to indicate either a FP or a trajectory starting in a neighborhood thereof.
• Regular oscillations of demixed species. The regimes represented in the first column of figure 3 and in the first column of figure 5 are regular, i.e. they feature a vanishing Lyapunov exponent. The initial states (t = 0) lie in elliptic islands centered around linearly stable FPs 1A and 3A, respectively (and their characteristic frequencies do not match Moser's commensurability condition [53]). The time evolution of boson populations consists in small oscillations around FPs. As a consequence, S mix (t) slightly oscillates around constant values S mix (1A) ≈ 0.08 and S mix (3A) ≈ 0.25 respectively, thus witnessing the persistence of a remarkable demixing.
• Fully developed mixing. The regimes illustrated in the second column of figure  3, second column of figure 4 and third column of figure 5 are chaotic, i.e. they are associated to a non-zero Lyapunov exponent. The initial configurations lies in the vicinity of linearly unstable FPs. The onset of chaos completely destroys the original, demixed configurations whose entropies of mixing are approximately equal to the ones relevant to the corresponding FPs, i.e. S mix (1B) ≈ 0.15, S mix (2B) ≈ 0.13 and S mix (3C) ≈ 0.23 respectively. As a consequence, in all three cases, S mix (t) repeatedly reaches the maximum possible value of ≈ 0.69 which, in turn, witnesses the full mixing of the bosonic species.
• Persistent demixing despite chaos. The dynamical regimes depicted in the first column of figure 4 and in the second column of figure 5 consist in small chaotic oscillations around the FPs 2A and 3B respectively. Chaos develops because the initial conditions already lie in the chaotic sea which, in turn, has been shown to surround linearly unstable FPs. However, despite the occurrence of chaos, not only the demixing of bosonic species persists, but also the macroscopic structure of the initial configurations remains unchanged during the dynamics. As a consequence, the oscillations of S mix (t) are chaotic but their amplitude is small, namely it never exceeds critical values ≈ 0.10 and ≈ 0.44 respectively. The dynamical regimes illustrated in the third column of figure 3 and in the third column of figure 4 are chaotic as well, but much less shrunken. It is a fact that, also in these cases, despite the presence of chaos, the two atomic species feature a low degree of mixing for all the simulated dynamics (in fact S mix (t) remains smaller than ≈ 0.29 and ≈ 0.27 respectively). Nevertheless, contrary to trajectories 2A and 3B (which consist in small chaotic oscillations around the equilibrium points), in these cases, chaos disrupts the structure of the original configurations, repeatedly triggering populations inversions. In other words, for what concerns trajectories 1C and 2C, as time goes by, boson populations in each well severely change but always in such a way to preserve the (low) degree of mixing.  Table 7 is intended to summarize the aforementioned results. Interestingly, we observe that the persistence of spatial phase separation marks extended bundles of chaotic trajectories, i.e. chaotic trajectories involving wide patches of starting points and rather extended ranges of model parameters. Such a persistence can be explained either in terms of energy conservation or by recalling the presence of regular islands where chaotic trajectories cannot enter.
Concerning the energy-conservation argument, we remind that the choice of a certain z 0 automatically fixes the trajectory and the constant-energy hypersurface Γ where the trajectory will be confined. Then, we analytically determine the maximum value of the entropy of mixing over the entire Γ,S mix , and note that the entropy of mixing along the trajectory, S mix (t), will never exceedS mix at any time. In other words, (see also Appendix C, where we comment on the structure of phase space), the trajectory will wander wide regions of Γ but, ifS mix is small enough, it will never visit highly-mixed configurations. The valuesS mix are determined over surfaces Γ by maximizing the objective function S mix under the constraints j n j = N , j m j = N and H = H( z 0 ), for each initial state z 0 . Such computation ofS mix as a function of model parameters is based on the well-known method of Lagrange multipliers. The result is shown figure 6. As an example, compare the green domain around point 1C in figure 6 (witnessing persistent spatial phase phase separation) with the light-blue domain around point 1C in the third row of figure 2 (displaying a manifestly chaotic behaviour). The same reasoning holds also for point 2C. Concerning the regular-islands argument, the latter is suggested by several examples of trajectories which, despite exhibiting a chaotic behaviour together with persistent spatial phase separation, are associated to the maximum possible (i.e. energetically accessible) value ofS mix := max Γ {S mix }. For example, the trajectory 3B is of this type since, as shown in figure 6, the constant-energy hypersurface Γ where it is embedded, features the biggest possible entropy of mixing, log 2 (depicted in dark red). A reasonable interpretation of this apparent mismatch is linked to the possible presence of regular islands on Γ [60] (see also Appendix C). IfS mix lies inside such islands, in fact, chaotic trajectories will never have the chance to visit highly-mixed configurations. This interpretation can be applied also to point 2A. A detailed analysis to ascertain the presence of mixed configurations inside regular islands requires an extended work that will be developed elsewhere.

Concluding remarks
We have investigated the dynamics of a bosonic binary mixture loaded in a three-well potential with periodic boundary conditions, its relation with the entropy of mixing and the robustness of spatial phase separation. In general, the developed analysis, even if focused on some particular classes of configurations, has provided a considerable amount of information about dynamical regimes characterized by regular and chaotic behaviours.
In section 2, we have introduced the model describing the mixture in the ring trimer and derived the corresponding semiclassical motion equations. Section 3 has been devoted to the presentation of the three notable classes of stationary configurations featuring demixing, "Dimer -Soliton", "SDW -Soliton" and "Soliton -Soliton". In section 4, we have developed the energetic-and the linear-stability analysis of the previously identified stationary configurations, highlighting their scope and limitations. In section 5, we have explicitly computed the first Lyapunov exponent along 102729 trajectories starting in the vicinity of as many FPs thus clearly identifying regular and chaotic regimes. We have observed that chaos can originate in three different ways: 1) When the trajectory starts in the vicinity of a linearly unstable FP; 2) When the trajectory starts in the neighbourhood of an elliptic FP such that its characteristic frequencies match Moser's commensurability condition [53]; 3) When the initial configuration lies outside the regular island centered around a linearly stable FP. In section 6, we have introduced the entropy of mixing S mix , borrowed from Statistical Thermodynamics [59], to quantify the degree of mixing between the two condensed species.
Eventually, in section 7 we have shown that the chaotic motion of boson populations and demixing can coexist or, in other words, that chaos, despite present, may not be able to completely disrupt the order imposed by phase separation. Such coexistence can occur either because highly mixed states lie in regular islands where the chaotic trajectory cannot penetrate or because the constant-energy hypersurface does not contain mixed states at all. In conclusion, we notice that our study could be of interest both for further theoretical investigations and for future experiments with bosonic mixtures trapped in ring lattices which, as discussed in the introduction, are within the reach of current experimental technology.

Appendix A. Stationary configurations featuring demixing
Figures A1-A3 show in detail how the stationary configuration changes in a non-narrow range of model parameters. They account both for the variations of boson populations n j , m j (different color shades) and for the relative phase between the wells (numbers within the various regions). Notice that white areas, representing the absence of a certain class of FPs, stand in between different sub-classes. The latter differ in the relative phases between the wells and are three in the "Dimer -Soliton" case, two in the "SDW -Soliton" case and in the "Soliton -Soliton" case as well (see figures 3, 4 and 5 respectively). We conclude by observing that the collective phase difference between the two condensed species play no role in the dynamics of the system (see Hamiltonian (2)), so figures 3-5 have been drawn arbitrarily setting Φ 1 = Ψ 1 = 0.
Appendix B. Remarks on the dynamical behaviour of trajectories starting close to a fixed point In the following, we describe the qualitative behaviour of a trajectory which, at t = 0, starts from a point z 0 very close to a FP z * . To this end we consider 4 different kinds of FPs z * corresponding to different kinds of stability/instability.  Linearly-stable FP. The solutions of the linearized system of differential equationṡ y = J( z * ) y induced by matrix (9) correspond to small oscillations around a given equilibrium point z * [52], the characteristic frequencies thereof being ω j = I{λ j }, where λ j 's are the eigenvalues of matrix (9). The effectiveness of the solutions of the linearized equations in representing those of the actual non-linear equations has been discussed in Figure A3. Class of FPs of the type Soliton -Soliton. Each column corresponds to a well while each row to a different condensed species. The color corresponds to the fraction of bosons hosted by the well (see color bars) while numbers 0 and π indicate the phase shift with respect to the first well. In the white region, this kind of stationary configuration does not exist. Energetically-stable FP. Energetic stability is stronger than linear stability, as one can prove that if the initial configuration z 0 is sufficiently close to FP z * , solutions of the actual non linear system (6) remain arbitrarily close to those of the linearized one for all times and, moreover, there are no issues associated to the commensurability of characteristic frequencies ω j 's.
Linearly-unstable FP. A point of linear instability (also called dynamic instability [45]) is such that almost every trajectory will depart from it. For a generic (i.e. non necessarily Hamiltonian) dynamical system, a trajectory starting close to an unstable FP can have any sort of behaviour (e.g. fall towards a FP, towards a periodic orbit, become chaotic, etc.). Since we are dealing with an Hamiltonian system, one knows a priori that the relevant flow in the phase space is incompressible, so the number of possible alternatives decreases. Despite not rigorously proven §, as the number of effective degrees freedom is relatively high (D = 4), one can expect that linearly unstable FPs are surrounded by chaos (see Appendix C). The validity of this reasonable ansatz is confirmed by the explicit calculation of the first Lyapunov exponent along trajectories starting close to FPs z * (compare second and third row of figure 2).
Energetically-unstable FP. An energetically-unstable FP can be linearly stable (but not viceversa), so the qualitative behaviour of a trajectory starting from its neighborhood strongly depends on the linear stability of z * and, as already discussed, on the possible commensurability of characteristic frequencies ω j 's. Although we have restricted our analysis to isolated systems at zero temperature, it is worth mentioning that, if dissipation is introduced or in presence of thermalization processes, an energeticallyunstable system will tend to decay to an energetically-stable state [61].
Appendix C. On the structure of phase space The phase space associated to Hamiltonian dynamical system 6 seems to be 12dimensional as, for each of the two condensed species and for each of the three wells, there are two canonically conjugate dynamic variables: local boson number n j (m j ) and local phase Φ j (Ψ j ). As discussed in section 3, the relative phase between the two condensed species play no role in the dynamics (see Hamiltonian (2)) so one can arbitrarily fix Φ 1 = Ψ 1 = 0. Moreover, as the total number of bosons N = 3 j=1 n j and M = 3 j=1 m j constitute two independent conserved quantities, one can substitute n 1 = N − n 2 − n 3 and m 1 = M − m 2 − m 3 . Therefore, the number of effective dynamical variables is 8, which correspond to D = 4 degrees of freedom. It is well known, also in the field of ultracold atoms [45], that there is a profound difference between systems featuring D = 2 and D > 2 degrees of freedom [62]. Only in the first case, the threedimensional space corresponding to the constant-energy (H = E) hypersurface can be divided by the relevant two-dimensional KAM tori into separated regions. Chaotic trajectories, if present, are therefore always topologically confined by KAM tori. For D > 2, instead, the latter cannot divide the phase space into separated regions (in the same way as a circumference cannot divide the euclidean space R 3 into two parts). All chaotic regions are therefore interconnected by a very slow percolation-like phenomenon which goes under the name of Arnold diffusion [62,45]. Nevertheless, they do not occupy the whole constant-energy space, since regular islands are still present (e.g. the neighbourhoods of energetically-stable FPs), the relative measure thereof being still an open problem [60]. In passing, we mention that, when the number of degrees of freedom of a non linear dynamical system tends to infinite, the measure of regular islands tends to zero and so the constant-energy hypersurface gets completely chaotic, thus justifying the ergodic hypothesis and a microcanonical approach to the problem [60].