Spatial Phase Separation of a Binary Mixture in a Ring Trimer

TL;DR

This paper explores the phase separation phenomena in a binary bosonic mixture confined in a ring trimer, combining semiclassical and quantum analyses to reveal critical behavior, entanglement, and persistent demixing in dynamical regimes.

Contribution

It introduces a comprehensive study of phase separation in a ring trimer, including quantum indicators and dynamical evolution, extending beyond ground state analysis.

Findings

Critical phenomena in phase separation are identified.

Quantum indicators corroborate semiclassical predictions.

Persistent demixing occurs in chaotic dynamical regimes.

Abstract

We investigate the phase separation mechanism of bosonic binary mixtures in spatially-fragmented traps, evidencing the emergence of phases featuring a different degree of mixing. The analysis is initially carried out by means of a semiclassical approach which transparently shows the occurrence of critical phenomena. These predictions are actually corroborated by the study of genuinely quantum indicators, including, but not limited to, the energy levels' structure and the entanglement between the species. The scope of our work goes also beyond the ground state's properties, as it comprises excited states and the dynamical evolution thereof. In particular, after introducing an indicator to monitor the degree of mixing, we show that several dynamical regimes feature persistent demixing in spite of their remarkably chaotic character.

Tables

Table 1. Table 1: Boson populations in the partially (de)mixed phase, in the limit N → + ∞ → 𝑁 N\to+\infty .

$x_1=0$				$x_2=\frac{1+\frac{W}{U}}{2+\frac{W}{U}}$				$x_3=\frac{1}{2+\frac{W}{U}}$
$y_1=\frac{1+\frac{W}{U}}{2+\frac{W}{U}}$				$y_2=0$				$y_3=\frac{1}{2+\frac{W}{U}}$

Equations

$$\hat{H}=-T_a\left(a_L{a}_R^{†}+a_R{a}_L^{†}\right)+\frac{U_a}{2}\left[n_L(n_L-1)+n_R(n_R-1)\right]$$

$$-T_b\left(b_L{b}_R^{†}+b_R{b}_L^{†}\right)+\frac{U_b}{2}\left[m_L(m_L-1)+m_R(m_R-1)\right]+W\left(n_L{m}_L+n_R{m}_R\right),$$

$$V_{dimer}=-2N{T}_a\sqrt{x_L{x}_R}-2M{T}_b\sqrt{y_L{y}_R}+\frac{U_a{N}^2}{2}\left(x_L^2+x_R^2\right)+\frac{U_b{M}^2}{2}\left(y_L^2+y_R^2\right)$$

$$+WNM\left(x_L{y}_L+x_R{y}_R\right)$$

$$\frac{W}{U}<1+\frac{2T}{UN}.$$

$$x_L=\frac{1}{2}+\frac{1}{2}\sqrt{1-{\left[\frac{2T}{N(W-U)}\right]}^2},x_R=1-x_L$$

$$y_L=\frac{1}{2}-\frac{1}{2}\sqrt{1-{\left[\frac{2T}{N(W-U)}\right]}^2},y_R=1-y_L,$$

$${\omega }_{1,2}=\sqrt{T\left(T+N\frac{U\pm W}{2}\right)}$$

$$E{E}_{V_1-V_2}=-{\mathrm{Tr}}_{V_1}({\hat{\rho }}_{V_1}{\log }_2{\hat{\rho }}_{V_1})$$

$${\hat{\rho }}_{V_1}={\mathrm{Tr}}_{V_2}\left({\hat{\rho }}_0\right)$$

$$|{\psi }_0⟩=\frac{1}{\sqrt{2}}\left(|N,0{⟩}_L|0,N{⟩}_R+|0,N{⟩}_L|N,0{⟩}_R\right)$$

$${\hat{\rho }}_T=\frac{1}{\mathcal{Z}}\sum _n{e}^{-\beta E_n}|{\psi }_n⟩⟨{\psi }_n|,\text{where}\mathcal{Z}=\sum _n{e}^{-\beta E_n}\text{and}\hat{H}|{\psi }_n⟩=E_n|{\psi }_n⟩$$

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

$$+W\sum _{j=1}^3{n}_j{m}_j$$

$$V_{trimer}=-2N{T}_a\left(\sqrt{x_1{x}_2}+\sqrt{x_2{x}_3}+\sqrt{x_3{x}_1}\right)+\frac{U_a{N}^2}{2}\left(x_1^2+x_2^2+x_3^2\right)$$

$$-2M{T}_b\left(\sqrt{y_1{y}_2}+\sqrt{y_2{y}_3}+\sqrt{y_3{y}_1}\right)+\frac{U_b{M}^2}{2}\left(y_1^2+y_2^2+y_3^2\right)$$

$$+W(x_1{y}_1+x_2{y}_2+x_3{y}_3)$$

$${\omega }_{1,2}=\sqrt{3T\left(3T+\frac{2UN}{3}\pm \frac{2NW}{3}\right)}$$

$$\frac{W}{U}<1+\frac{9T}{2UN}$$

$$E{E}_{a-b}=-{\mathrm{Tr}}_a\left({\hat{\rho }}_a{\log }_2{\hat{\rho }}_a\right).$$

$$\sqrt{n_j}{e}^{i{\varphi }_j}\approx a_j,\sqrt{m_j}{e}^{i{\psi }_j}\approx b_j,$$

$$S_{mix}=-\frac{1}{3}\sum _{j=1}^3\left[\frac{n_j}{n_j+m_j}\log \left(\frac{n_j}{n_j+m_j}\right)+\frac{m_j}{n_j+m_j}\log \left(\frac{m_j}{n_j+m_j}\right)\right].$$

Full text

1 Introduction

The localization of two interacting condensed species in different spatial regions is a phenomenon that occurs in presence of strong interspecies repulsive couplings. The mechanism underpinning this spatial phase separation has been deeply investigated in systems of ultracold bosons within the mean-field approach to condensate’s dynamics [1, 2, 3]. It has been evidenced that a number of elements, including the depth of the trapping potential, the number of particles and the intraspecies interactions, influence the emergence of the mixing-demixing transition and the way it develops [4]. From the experimental side, the unprecedented control of quantum matter has allowed to realize multicomponent [5, 6, 7, 8] and spatially-fragmented [9, 10, 11] condensed systems whose rich phenomenology, in turn, rekindled the interest of the theoretical community for a deep comprehension of the demixing transition in the presence of optical lattices [12, 13].

Systems made up of two condensed species (be they heteronuclear [14] or homonuclear [15]) indeed offer a plethora of observable physical phenomena ranging from the dynamics of solitons [16, 17, 18] to the formation of entanglement between the species [19, 20], not to mention the mechanisms that directly account for phase separation [21, 22, 23].

In this context, the present manuscript aims to discuss the phenomenon of spatial phase separation in the two simplest possible geometries, the dimer and the trimer, whose fragmented character strongly affects the mixing properties of the bosonic binary mixture. Such systems are thoroughly analyzed in section 2 and 3 respectively. In both sections we start presenting the Bose-Hubbard (BH) Hamiltonian, the second-quantized model which well captures the essence of the physical problem. We then switch to a semiclassical perspective, as we think that it allows for the emergence of the mixing-demixing transition at its most transparent. We support the predictions provided by the aforementioned approximation scheme by introducing genuinely quantum indicators. The latter, which include the energy spectrum and the entanglement entropy (EE), prove to be singular where anticipated by the semiclassical approach. Eventually, in section 3.4, we extend our analysis beyond the ground-state physics by providing an overview on the interplay between chaotic dynamics and persistent demixing.

2 The binary mixture in a double well

In order to effectively model a system of two condensed species in a two-well potential, we introduce operators $a_{L}$, $a_{R}$, $b_{L}$ and $b_{R}$. Operator $a_{L}$ $(a_{R})$ destroys a species-a boson in the left (right) well. Similarly, operator $b_{L}$ $(b_{R})$ destroys a species-b boson in the left (right) well. Standard bosonic commutators hold, namely $[a_{X},a_{Y}^{\dagger}]=\delta_{X,Y}=[b_{X},b_{Y}^{\dagger}]$ and $[a_{X},b_{Y}]=$ $[a_{X},b_{Y}^{\dagger}]=0$, with $X,Y\,\in\{R,L\}$. In this framework, the BH Hamiltonian associated to the system reads

[TABLE]

[TABLE]

where $n_{X}=a_{X}^{\dagger}a_{X}$ ($m_{X}=b_{X}^{\dagger}b_{X}$) is the numbers of species-a (species-b) bosons in well $X\in\{R,L\}$. This model includes tunnelling terms $T_{a}$ and $T_{b}$, intraspecies interactions $U_{a}$ and $U_{b}$, and a density-density interspecies coupling $W$. The boson number of each atomic species, namely $N:=n_{L}+n_{R}$ and $M=m_{L}+m_{R}$ represent two independent conserved quantities, being $[N,\hat{H}]=[M,\hat{H}]=0$.

2.1 CVP approach to the mixing-demixing transition

In order to investigate the miscibility properties of the ground state of Hamiltonian (1), we switch to a description of boson populations in terms of continuous variables $x_{X}:=n_{X}/N$ and $y_{X}:=m_{X}/M$, with $X\in\{L,R\}$. This technique, successfully used in the study of spatially-fragmented BECs [24, 25, 26, 27], allows one to recast the computation of the ground state of quantum Hamiltonian (1) into the search for the global minimum of the classical effective potential

[TABLE]

[TABLE]

In the following, we will assume, for simplicity, $T_{a}=T_{b}=:T$, $U_{a}=U_{b}=:U$ and $N=M$ but we note in advance that small deviations from this ideal scenario do not substantially affect our analysis. As the intuition suggests, if the interspecies repulsion $W$ is small compared to the intraspecies repulsion $U$, than the two species are mixed and uniformly distributed in the two wells. More specifically, one can verify that the uniform configuration $x_{X}=y_{X}=1/2$ (with $X\in\{L,R\}$) is indeed a global minimum of effective potential (2) provided that

[TABLE]

Above this critical value, the two species start to separate in different spatial regions, i.e. to localize in different wells. The analytic expressions of boson populations in the two wells,

[TABLE]

[TABLE]

show that, for $W\to+\infty$, the two species completely demix, meaning that species-a bosons occupy just the left well while species-b bosons localize just in the right well (namely, $x_{L}=1$, $x_{R}=0$, $y_{L}=0$, $y_{R}=1$). We remark that, due to the system symmetry, formulas (4) and (5) represent just one of the two minimum-energy configurations, the other one being obtained swapping labels “$L$” and “$R$”. It is interesting to notice (see Fig. 1) that, all other things being equal, the bigger the tunnelling amplitude $T$, the bigger the interspecies repulsion $W$ needed to trigger species demixing, the smoother the transition from the fully mixed to the completely demixed configuration.

2.2 Spectrum collapse at the critical point

The energy level distribution considerably varies according to the ratio $W/U$. As depicted in Fig. 2, in fact, the energy levels’ structure is very different at the two sides of the critical point and exhibits a collapse where the mixing-demixing transition occurs.

This circumstance, evidenced by means of an exact numerical diagonalization of Hamiltonian (1), can be better highlighted if one resorts to the Bogoliubov technique for the computation of quasi-particles frequencies. Within this approximation scheme [25], and in agreement with the results of the CVP approach, characteristic frequencies

[TABLE]

are well defined, in fact, only in that region of parameters’ space identified by inequality (3).

2.3 Entanglement entropy as a critical indicator

The study of the correlation properties between two partitions of the system’s Hilbert space has proved to be an effective indicator for detecting the occurrence of mixing-demixing transitions [20, 27]. The Entanglement Entropy (EE) is the indicator commonly used to quantify these quantum correlations and the standard way to compute it involves the subdivision of the space of states into two partitions, $V_{1}$ and $V_{2}$. More specifically, the entanglement between $V_{1}$ and $V_{2}$ is given by

[TABLE]

and corresponds to the Von Neumann entropy of the reduced density matrix

[TABLE]

which is obtained, in turn, tracing out the degrees of freedom of $V_{2}$ from the ground state’s density matrix $\hat{\rho}_{0}=\lvert\psi_{0}\rangle\langle\psi_{0}\rvert$. We remark that the partitioning-step is a crucial one and that a strong entanglement in a certain partition may correspond to a weak (or null) entanglement in another one. In [20], we have investigated three different kinds of EE according to three different ways of partitioning the two-species-dimer system: i) $EE_{L-R}$, between the site modes; ii) $EE_{s-d}$ between momentum modes and iii) $EE_{a-b}$ between species modes. Interestingly, we have evidenced that they build up, respectively, due to the presence of a hopping amplitude $T$ between the sites, because of an on-site repulsion $U$ and owing to the interspecies coupling $W$. On top of that, as clearly illustrated in Fig. 3, all three indicators feature singularities (either peaks or jumps) where the mxing-demixing transition takes place.

Moreover, one can notice that, for large interspecies interactions, the entanglement between the species $(EE_{a-b})$ tends to the limiting value $1=\log_{2}2$, a number which is reminiscent of the two-sided structure of the system’s ground-state. In that circumstance, in fact, the latter, can be well approximated by the Schrödinger-cat-like state

[TABLE]

where state $\lvert n,m\rangle_{L}\lvert N-n,N-m\rangle_{R}$ is such that there are $n$ species-a and $m$ species-b bosons in the left well.

The concept of EE, which is a quantity strictly defined for pure states, can be readily generalized to non-zero temperatures by substituting $\hat{\rho}_{0}$ with the thermal density matrix

[TABLE]

in formulas (6) and (7). This indicator, which goes under the name of bipartite residual entropy, despite being affected by classical correlations between the subsystems and by the classical entropy, is still able to detect the occurrence of mixing-demixing transitions when these occur at finite temperature [20].

3 The binary mixture in a ring trimer

We turn to the analysis of the mixing properties of a bosonic binary mixture in a ring trimer, a system where the presence of an additional well discloses a phenomenology much richer than the one discussed in Section 2. Similarly to what discussed in relation to the two-species dimer, we start with introducing the second-quantized operator $a_{i}$ ($b_{i}$) which annihilates a species-a (species-b) particle in the $i-$th well (where $i=1,2,3$) and whose bosonic character is expressed by the following commutation relations: $[a_{i},a_{j}^{\dagger}]=\delta_{i,j}$, $[b_{i},b_{j}^{\dagger}]=\delta_{i,j}$. Within this formalism, the model of a binary mixture in a three-well potential with periodic boundary conditions is effectively portrayed by two BH Hamiltonians (each one being associated to a bosonic species and involving three spatial modes) coupled by an interspecies repulsion term. This reads:

[TABLE]

[TABLE]

where the periodic boundary conditions implied by the ring geometry are such that $j=4\equiv 1$. After recalling that $n_{j}=a_{j}^{\dagger}a_{j}$ and $m_{j}=b_{j}^{\dagger}b_{j}$, we remark that the total number of bosons in the two species, $N=\sum_{j}n_{j}$ and $M=\sum_{j}m_{j}$, constitute, as in the dimer case, two independent conserved quantities. Terms proportional to the hopping amplitudes $T_{a}$ and $T_{b}$ favour particles’ delocalization in the ring geometry; conversely, terms proportional to intraspecies repulsive interactions $U_{a}$ and $U_{b}$ favor particles’ localization. Eventually, the interspecies-repulsion term ($\propto W$) is responsible for the spatial phase separation of the two bosonic species.

3.1 CVP approach to the two-step demixing transition

As explained in Section 2.1, an effective way to investigate the phase separation mechanism is to switch to the CVP, an approximation scheme where the inherently discrete eigenvalues of number operators $n_{j}$ and $m_{j}$ are replaced by continuous variables $x_{j}=n_{j}/N$ and $y_{j}=m_{j}/M$. In this framework, the ground state of quantum Hamiltonian (9) corresponds to the global minimum of classical potential

[TABLE]

[TABLE]

[TABLE]

subject to the constraints $0<x_{j}<1$, $0<y_{j}<1$ $\forall j$, $\sum_{j}x_{j}=\sum_{j}y_{j}=1$. For the sake of simplicity, in the following we will assume that $U_{a}=U_{b}=:U$, $T_{a}=T_{b}=:T$ and $N=M$, although small offsets from this ideal condition do not substantially affect the results.

The basic difference with respect to the two-species dimer is that, in this more complex geometry, tuning the ratio $W/U$ one triggers two mixing-demixing transitions, rather than only one. In other words, one can recognize three different phases, which differ in the degree of mixing:

Fully mixed phase. If $W/U\in\left(0,\,1+\frac{9T}{2UN}\right)$ then the two species are uniformly distributed in the ring trimer, meaning that $x_{j}=y_{j}=1/3\,\forall j$. 2. 2.

Partially (de)mixed phase. If $W/U\in\left(1+\frac{9T}{2UN},\,R_{c}\right)$, where $R_{c}=R_{c}(T,U,N)>2$ then one can recognize a rather low degree of mixing in two out of three wells and an almost perfect mixing in the remaining one. It is worth noticing that, in analogy with the two-species dimer (see section 2.1), increasing the ratio $T/(UN)$ has the twofold effect of shifting critical points $1+\frac{2T}{UN}$ and $R_{c}$ and of smoothing the transitions. 3. 3.

Fully demixed phase. If $W/U>R_{c}$ then the two species completely separate, meaning that one species occupies two wells while the second species conglomerates in the remaining site.

This scenario, illustrated in the second row of Fig. 4, turns even clearer and sharper if one considers the semiclassical limit $N\to+\infty$ (or, equivalently, the no-tunnelling limit $T\to 0$). In this circumstance, in fact, a fully-analytic investigation of the two-step mixing-demixing transition can be carried out [27]. Such analysis is rather complex because the non-convex character of effective potential (10) makes it necessary an exhaustive exploration of its domain boundaries. The result is that the partially (de)mixed phase occurs for $W/U\in(1,2)$, a range where boson populations behave as shown table 1.

We observe that the aforementioned configuration, depicted in the first row of Fig. 4, is not the only one which minimizes effective potential (10). The latter, in fact, owing to the system’s symmetry, features three isoenergetic points of minimum which are equal up to cyclic permutations of the site’s indexes. The link between symmetry and number of minimum-energy configurations will be made clearer in section 3.3.

3.2 Spectrum collapse at the mixing-demixing transitions

We have performed an exact numerical calculation of the energy levels of Hamiltonian (9) and the result, shown in Fig. 5, clearly illustrates the spectrum collapse and reorganization across the two critical points.

In this regard, it is worth noticing that, the smaller the ratio $T/(UN)$, the sharper the collapse, the more abrupt the transition.

A deeper physical insight can be gained from the study of quasi-particle frequencies which feature the fully mixed phase. Within the Bogoliubov approximation scheme and making use of the dynamical algebra method, we were able to determine such frequencies for a ring made up of a generic number of sites [28]. In the specific case of a ring trimer, their explicit expression

[TABLE]

shows that they are well defined only if

[TABLE]

an equality constituting the same condition already emerged in the context of the CVP approach (see section 3.1).

3.3 Entanglement between the species

Among the different partition schemes which can be applied to the Hilbert space of states associated to Hamiltonian (9), we focus on the one involving species-a and species-b modes. To be more specific, after numerically computing the ground state $\lvert\psi_{0}\rangle$ of (9) and the relevant density matrix $\hat{\rho}_{0}=\lvert\psi_{0}\rangle\langle\psi_{0}\rvert$, one obtains the reduced density matrix $\hat{\rho}_{a}=\mathrm{Tr}_{b}\left(\hat{\rho}_{0}\right)$ associated to species-a bosons by tracing away the degrees of freedom relevant to species-b bosons. The entanglement between the species can be thus computed as

[TABLE]

With reference to the left panel of Fig. 6, where we illustrate how indicator (11) changes as a function of the ratio $W/U$, we can notice that the entanglement between the species: i) is correctly null when they do not interact, i.e. when $W=0$; ii) features peaks where the mixing-demixing transitions occur; iii) tends to the limiting value $\log_{2}(6)$ for large interspecies repulsions.

Concerning the peaks, in agreement with what already discussed in section 3.1, one can observe that they broaden and move rightward upon an increase of ratio $T/(UN)$, a circumstance which mirrors the fact that the hopping amplitude $T$ has the twofold effect of smoothing the transitions and of delaying their occurrence. As regards the argument of the logarithm, six is the number of classical configurations which minimize effective potential (10). Such configurations, depicted in the right panel of Fig. 6 and which are equal up to species labels exchange and up to cyclic permutations of the site indices, are quantum-mechanically reabsorbed in the formation of a unique ground state. The latter, in fact, exhibits a six-faced Schrödinger-cat-like structure [27] whose explicit expression can be readily obtained properly generalizing equation (8).

3.4 Chaos and persistent demixing

So far, our discussion has been focused on the miscibility properties of the two species when the system is in its lowest energy state. A much broader panorama of physical phenomena emerges if one extends the analysis to excited states. In view of this complexity, which makes the scenario remarkably rich and branched, an exhaustive exploration of all possible dynamical regimes would go beyond the scope of this work. For this reason, we restrict our investigation to the dynamics of highly demixed configurations [29]. More specifically, we numerically simulate extended bundles of trajectories starting in the close vicinity of stationary states featuring a low degree of mixing. This is done in the semiclassical limit, where it is possible to replace second quantized operators $a_{j}$ and $b_{j}$ with local order parameters

[TABLE]

complex numbers which, in turn, depend on local populations ($n_{j}$ and $m_{j}$) and local phases ($\phi_{j}$ and $\psi_{j}$). After finding the semiclassical version of Hamiltonian (9) and the motion equations that the latter induces, one can therefore indentify different classes of stationary configurations featuring a trivial time evolution and a low value of the indicator

[TABLE]

The latter, which goes under the name of Entropy of Mixing, is commonly employed in Physical Chemistry to quantify the miscibility of chemical compounds [30], but its scope can be effectively extended to the realm of ultracold atoms, where it constitutes a reliable tool for monitoring the mixing of different condensed species.

Interestingly, by exploiting this indicator, we have shown [29] that many dynamical regimes feature persistent demixing in spite of a strong chaotic behaviour. For example, Fig. 7 clearly shows that the two species, despite their chaotic evolution, tend to exclude each other throughout the simulated dynamics. This kind of phenomenology, which has been evidenced by means of an explicit calculation of the first Lyapunov exponent and of indicator (12) along extended sets of trajectories, can be interpreted either in terms of energy conservation or by means of dynamical considerations. Both these motivations, pictorially sketched in the two panels of Fig. 8, are based on the analysis of the constant energy hypersurface and of the orbits embedded therein.

Concerning the “energetic argument”, it is possible to show that the constant-energy manifold where the dynamics develops may not contain highly mixed states at all. In this circumstance, therefore, no matter how large the first Lyapunov exponent is, species mixing will never occur. As regards the “dynamical argument”, after recalling that the constant-energy hypersurface of low-dimensional systems features a coexistence of both regular and chaotic orbits [31], one can notice that, in the system under analysis, the former constitute regular islands where highly mixed states are encapsulated and where chaotic trajectories cannot enter. In this scenario, therefore, mixing is energetically possible, but actually forbidden by dynamical restrictions.

4 Conclusions and future perspectives

In this work we have provided an overview on the rich phenomenology which distinguishes bosonic binary mixtures in low-dimensional systems. In section 2, we started our discussion with the two-well (dimer) system, the simplest possible geometry whose spatially-fragmented character has a remarkable impact on the miscibility properties of the species. By means of a semiclassical approach, we have evidenced the occurrence of a single mixing-demixing transition, from a fully mixed to a completely demixed phase, providing the analytic expression of boson populations for all the possible values of the control parameter $W/U$. To support the semiclassical analysis, we have also investigated genuinely quantum indicators, such as the energy spectrum and the entanglement entropy, highlighting their singular behaviour at the critical point. In section 3, we turned to the trimer. The higher complexity of its geometry (given by the presence of an additional well) is mirrored by a more sophisticated phase-separation mechanism. In this context, in fact, an additional intermediate phase, neither fully mixed nor completely demixed, can be identified. The presence of two critical points is confirmed by the doubly-critical character of the energetic fingerprint and of the entanglement between the species. Interestingly, as in the dimer system, all the indicators that we have considered show that the hopping amplitude $T$ has the twofold effect of smoothing the transitions and of delaying their occurrence. In the last part of the manuscript, in order to go beyond the ground state’s physics, we have extended our analysis to excited states, investigating the interplay between persistent demixing and chaotic behaviour.

We hope to have made it clear that phase separation phenomena in low-dimensional systems constitute an extraordinary rich panorama and their analysis can touch on several different aspects, ranging from entanglement properties to quantum dynamics. The work that has been done so far has shined light just on those systems featuring the most elementary geometries (dimer and trimer). We therefore aim to extend it to more complex geometries where, presumably, multiple intermediate phases will come into play. Moreover, we would like to generalize our analytic and numerical computations to the case of non-symmetric species, i.e. $U_{a}\neq U_{b}$, $T_{a}\neq T_{b}$ and $N\neq M$, a situation where the presence of additional control parameters should offer an even more intriguing physics. Current work is focused on the study of attractive interspecies interactions, which constitutes a realm on its own, as they lead to the emergence of phases featuring supermixed states rather than phase separation. We will report on this topic in a separate, dedicated work.