Residual Entropy and Critical Behavior of Two Interacting Boson Species in a Double Well

Motivated by the importance of entanglement and correlation indicators in the analysis of quantum systems, we study the equilibrium and the bipartite residual entropy in a two-species Bose–Hubbard dimer when the spatial phase separation of the two species takes place. We consider both the zero and non-zero-temperature regime. We present different kinds of residual entropies (each one associated with a different way of partitioning the system), and we show that they strictly depend on the specific quantum phase characterizing the two species (supermixed, mixed or demixed) even at finite temperature. To provide a deeper physical insight into the zero-temperature scenario, we apply the fully-analytical variational approach based on su(2) coherent states and provide a considerably good approximation of the entanglement entropy. Finally, we show that the effectiveness of bipartite residual entropy as a critical indicator at non-zero temperature is unchanged when considering a restricted combination of energy eigenstates.


Introduction
Systems formed by gases of ultracold bosons trapped in homogenous arrays of potential wells (optical lattices) [1] have attracted, in the last two decades, an enormous amount of attention due to the rich variety of phenomena they feature at zero temperature [2,3]. The physical properties of such quantum fluids have been shown to be mainly determined by the competition of boson-boson repulsive interactions with the tunneling effect between adjacent wells, causing the boson mobility through the lattice. Among many effects observed in such systems, one of the most significant is the famous superfluid-insulator transition in which, for a boson-boson interaction strong enough, the boson mobility is quantum-mechanically inhibited when the boson density takes integer values [4].
In this framework, introducing in the lattice a second bosonic species interacting with the primary one has allowed the realization [5,6] of the so-called binary quantum fluids whose complex phenomenology has revealed unexpected effects and behaviors. These are, for example, the formation of new types of insulating phases and superfluidity [7,8], quantum emulsions exhibiting a glassy character [9,10], the deformation of the insulating (Mott) domains accompanying the formation of polaron excitations [11,12], the presence of interspecies entanglement [13], and the spatial separation of the two species (demixing effect) [14,15].
The simplest possible lattice system in which the interplay of two species can be studied is represented by the two-species dimer (TSD), namely, a mixture trapped in a lattice with two wells (dimer). This system, sufficiently simple to allow the use of standard analytic approaches, is, however, complex enough to exhibit the space-localization effects distinguishing two-species mixtures in larger lattices at zero temperature. The TSD has allowed both a thorough study [16] of such behaviors when the interspecies (repulsive) interaction W > 0 is varied and the analytic derivation of the critical value of W for which the mixed species, equally distributed in the two wells, localize in two separated wells. Such an effect, called delocalization-localization (DL) transition, also takes place in the attractive case (W < 0) but in the final state both species occupy the same well for W strong enough. The DL transition, characterizing the ground state of the TSD, and the spectral collapse related to this phenomenon have been studied numerically in Ref. [16] and by means of the continuous-variable approach in Ref. [17]. The critical behavior of the TSD has been confirmed by resorting to quantum-correlation indicators such as the Fisher information, the coherence visibility and the entanglement entropy (EE). The latter, in particular, has proved particularly sensitive in detecting the macroscopic changes in the ground-state structure both for repulsive and for attractive interspecies interaction.
This motivates our interest for the entanglement entropy and, more in general, for the bipartite residual entropy at non-zero temperatures in the TSD. In this paper, we perform a systematic study of this correlation property effecting numerical simulations which include non-zero temperatures. We begin with studying the zero-temperature regime. To check the robustness of the EE, we determine its dependence from the interspecies interaction both by considering the exact ground state (calculated numerically) and by representing the ground state in terms of atomic coherent states (CS). The CS picture is interesting since, in addition to allow for fully-analytic calculations, its semiclassical character approximates the system ground state in a form closer to the preparation of the system in real experiments. As is well known, the EE describes the entanglement property of a physical system through the Von Neumann entropy of a suitably defined sub-system. In the sequel, we calculate the EE by partitioning the system in sub-systems such as (i) the left-well and the right-well bosons; (ii) bosons with zero and non-zero momentum; and (iii) the species-A and species-B bosons.
As noted above, a number of new quantum phases has been predicted in the last fifteen years whose distinctive feature is to manifest at zero temperature. On the other hand, after the realization of optical lattices trapping ultracold atoms, it has become more and more evident that reducing (and measuring) the temperature on the nanoscale is an outstanding problem [18]. For this reason, the detection of zero-temperature phase transitions such as the space separation in bosonic mixtures (or its simpler dimer version, the DL transition) must more realistically rely on indicators which are reminescent of the critical behavior of the system even when temperature is non-zero [19].
In this perspective, achieving a good control of the correlation properties for a system undergoing the DL transition at non-zero temperature certainly represents a useful tool for its observation in future experiments. For this reason, we have thoroughly explored the residual-entropy behavior at non-zero temperature and have tried to understand the effect of thermal fluctuations in regimes where they compete with quantum fluctuations. The bipartite residual entropy has been calculated numerically by exploiting the knowledge of the TSD exact spectrum. As in the zero-temperature case, we consider the reduced thermal density matrix for three different partition schemes of the system based on separating space modes, momentum modes and atomic species. Finally, to further test the bipartite residual entropy as a critical indicator, we have compared the exact bipartite residual entropy with that calculated using a restricted range of energy levels around the expected average energy.
The paper is organized as follows. In Section 2, we introduce the TSD model and review the DL transition discussing the change of structure it induces in the ground state and the spectral collapse, a significant property that marks the transition. Section 3 is devoted to defining the equilibrium and the bipartite residual entropy and the relation thereof with the EE. Section 4 contains the results of our numerical calculations of bipartite residual entropy within the previously discussed partition schemes at zero and non-zero temperature. In Sections 5 and 6, we compute the bipartite residual entropy in the coherent-state and in the restricted-basis approach, respectively. Section 7 is devoted to concluding remarks.

The Model and the Ground-State Properties
An effective description of ultracold bosons trapped in homogenous arrays of potential wells is provided by the Bose-Hubbard (BH) model [4] in which local boson operators A i and A + i represent the microscopic annihilation and creation processes, respectively, at the ith well. The experimental realization of this model is currently achieved by means of well-known optical-trapping techniques [3,20]. These, by combining counter-propagating laser beams, cause the formation of (optical) lattices, the sites of which correspond to effective local potentials attracting bosons. In the simplest possible case of a two-site lattice (a double potential well), the BH Hamiltonian is given by where L and R refers to the left and right well, respectively, U a is the boson-boson interaction and J a is the hopping amplitude controlling interwell boson exchange. The boson operators A L , A + L , A R , and A + R satisfy the standard commutator [A σ , A + σ ] = 1 with σ = L, R. If, in addition to species A, a second species B is introduced, the spatial modes become four, A L , A R , and B L , B R , for the species A and B, respectively. The resulting mixture is thus described by the two-species dimer Hamiltonian [16] in which, apart from the single-species BH Hamiltonians H a and H b , the significant term is that depending on interspecies interaction W. This couples the two species through the boson local populations described by the number operators N σ = A + σ A σ and M σ = B + σ B σ with σ = L, R. When the interspecies interaction W becomes sufficiently strong, the two interacting species trapped in a double-well potential feature macroscopic localization effects. In particular, a repulsive interaction tends to spatially separate the species into different wells while an attractive interaction tends to confine both species in the same well. This represents the DL transition. In the first case, this is characterized by an almost complete localization of the two species in different wells, and thus by a demixing effect, whereas, in the second case, the attractive interaction leads to a "supermixed" state with a localization of both species in a single well.
Such effects are confirmed by the numerical calculation of the ground state for different values of W. To see this, we note that the energy eigenstates can be suitably represented in the basis of space-mode Fock states where labels n σ and m σ , describing the local boson populations, are the eigenvalues of number operators N σ and M σ , respectively. The parametrization n L = i, m L = j, n R = N − i and m R = M − j has been assumed to include the property that both operator N = N L + N R and operator M = M L + M R (representing the total boson numbers of the two species) commute with Hamiltonian H and thus are conserved quantities. The factorized form of Eqaution (2) aims to better distinguish left-well from right-well populations. A generic quantum state is then represented as Determining the energy eigenstates thus amounts to calculating coefficients w i,j for which the eigenvalue equation H E⟩ = E E⟩ is fulfilled. For values of W small enough, the ground state E 0 ⟩ is approximated in terms of su(2) coherent states [21] in the repulsive and attractive case, respectively, well illustrating the space-localized distributions emerging from the delocalization-localization transition [16] and leading to Schrödinger cats with strongly localized component states. Figure 1, obtained by numerically calculating the ground state in the repulsive case for different W, supplies us with an exact description of the DL transition and of the macroscopic changes in the ground-state structure. A similar behavior characterize the DL transition in attractive case, but the two emerging peaks finally localize around i = j = 0 an i = 30, j = 40. The critical behavior of the DL transition has been studied analytically by resorting to the semiquantum approach where boson number operators are approximated in terms of continuous variables [17]. This method has provided the critical value of W at which the transition takes place in the case of twin species (J a = J b = J, U a = U b = U). In this approach, the Fock states essentially become wave functions depending on the new continuous variables while, for energies low enough, the energy-eigenvalue equation takes the form of the Schrödinger problem for a multidimensional harmonic-oscillator Hamiltonian. The extremal points of the corresponding potential allow one to determine the ground-state configuration, and, in particular, to find the formula defining, for large boson numbers (N = M ≫ 1), the transition critical point in the parameter space. Interestingly, when W approaches this critical value, the energy spectrum has been shown to undergo a collapse in which the inter-level separation tends to zero. This spectral collapse can be seen as the hallmark of the dynamical transition which features the macroscopic change in the structure both of the ground state (see the previous discussion) and, more in general, of the low-energy excited states described in Ref. [17]. The generalized version of the previous formula for a mixture in a L-well ring lattice has been derived in Ref. [22].

Equilibrium Entropy and Bipartite Residual Entropy
The third law of thermodynamics states that a perfect crystal at temperature T = 0 exhibits entropy S = 0. This entropy is defined as the Equilibrium Entropy S eq . However, several physical systems ranging from, e.g., water ice [23,24], carbon monoxide [25], highly pressurized liquid-helium [26], glass systems [27,28], proteins [29], and even black-holes [30,31], seem to manifest a residual content of information (corresponding to a residual entropy) for T → 0. The presence of such residual entropy S R has been generally associated with residual degrees of freedom at T = 0 such as, among others, ground-state degeneracy, residual structural disorder, geometrical frustration and entanglement. These physical phenomena act as sources of uncertainty preventing the possibility to acquire knowledge on the exact state of the system, thus resulting as possible sources of information (i.e., a finite, residual value of the entropy).
For quantum systems, the residual entropy can be related to the presence of entanglement in the ground-state through the entanglement entropy (As an example, consider a two-spin system such that its ground state is If, on one hand, Ψ 0 ⟩ is a pure state with no information at all if measured in the energy basis, on the other hand, the outcome of a measurement is not certain anymore if the spins of the two particles are measured. In that case a 50% chance to measure either ↑ ↓⟩ or ↓ ↑⟩ there exists, hence suggesting a residual information content of one bit at T = 0. For this reason in quantum systems, the entanglement measure in the ground state by means of the entanglement entropy is associated to a bipartite residual entropy.). Entanglement entropy is a measure of the "amount of entanglement" in the system. A standard and accepted way to quantify entanglement is through the Von Neumann entropy. What is measured by the bipartite entanglement entropy is the mutual information shared between two partitions of the physical system (e.g., Alice and Bob). Givenρ the density matrix of the system, and defining two partitions A, B of the Hilbert space H such that H = H A ⊗ H B , the bipartite Von Neumann entropy is defined as [32] is the reduced density matrix of partition A (B) obtained tracing out the degrees of freedom of B (A). From now on, we will refer to this indicator as bipartite residual entropy. Notice that, for the same system, in principle, there exists infinitely many possible ways to partitions the Hilbert space H in two parts. This leads to the consideration that, since the choice of the partition is arbitrary, the measure of entanglement, i.e., the bipartite residual entropy, cannot have a global character by definition. We shall see how this is indeed the case in Section 4 (and, more specifically, in Section 4.4) when we will compute the bipartite residual entropy for the TSD for different choices of the partition A-B.

Equilibrium Entropy in the TSD
According to quantum statistical mechanics [33,34], the expression of the equilibrium entropy S eq (T) can be derived from the expression of the density operator as where the (canonical) density operator [35] at finite temperature is defined aŝ with E n representing the energy eigenvalue associated to the energy eigenstate Ψ n ⟩. Combining Equations (6) and (7), one finds the explicit expression for the equilibrium entropy with ρ n = e −βE n Z.
Since the ground-state of Hamiltonian Equation (1) cannot be degenerate [34,36], expression Equation (8) is a good definition of equilibrium entropy for the TSD as, for T = 0, it exactly satisfies S eq (0) = 0. In Figure 2, we show the equilibrium entropy computed for the TSD as a function of the interspecies interaction W J and effective temperature Tk B J. At T = 0, one clearly sees that S eq = 0 for all values of W J (black-dashed line). By sufficiently increasing the temperature, two peaks appear at the boundary of the central region where the phase transitions between the mixed and demixed phases occur ( W J ≈ 0.16). Their presence is related to the energy-level spectral collapse: this, in fact, entails that a bigger and bigger number of energy levels significantly contributes to the thermal density matrix used to calculate this quantity. Such peaks progressively vanish due to fluctuations when the temperature increases. In Section 4, we will show that, in general, bipartite residual entropy exhibits similar features. According to Figure 2, the equilibrium entropy tends to S eq = 1 for T ≠ 0 and if W J is large enough (plot tails). This reflects the fact that two dominating states (those corresponding to the lowest energies E 1 and E 0 ) provide contributions of about 1 2 log 2 2 to the limiting value S eq = 1. It is important to notice that, in both tails, the first excited level E 1 can be shown to tend to the ground state energy E 0 as a consequence of the spectral collapse characterizing the TSD. Accordingly, the smallest non-zero temperature that has been considered (Tk B J = 10 −4 ) is large enough to populate in a nearly equal way both the ground state and the first excited level because their separation E 1 − E 0 becomes smaller and smaller for large W J.
Notice also that, in this regime, the splitting between E 0 and E 1 decreases exponentially (with the number of particles) to a point that may lie below the actual experimental limit (see Appendix B for details). We also note that, at high temperatures, the equilibrium entropy tends to the value S eq ≃ 10.31 = log 2 D, where D is the number of energy levels (i.e., the dimension of the Hilbert space), showing the fact that all the energy eigenstates are equiprobable with probability 1 D.

Bipartite Residual Entropy in the TSD
In Ref. [16], we showed that the TSD manifests non trivial entanglement properties (relevant to the boson distribution in the two wells) in the ground-state suggesting the presence of a bipartite residual entropy at T = 0. This residual information at T = 0 is not grasped by Equation (8) as it exhibits S eq (0) = 0.
A consistent and different definition of the entropy is therefore required in order to be able to correctly describe the residual quantum information hidden in the ground-state structure. This can be naturally done by extending the definition of entanglement entropy Equation (5) at finite T in the way suggested by the expression for the equilibrium entropy Equation (8). We will call this quantity bipartite residual entropy at finite temperature S R (T) in order to distinguish it from the equilibrium entropy S eq (T) of expression Equation (8).
The key difference between definitions Equations (5) and (6) lies in the fact that, in the entanglement entropy, a reduced density operatorρ A is used. Given a partition of the Hilbert space, the reduced density operator of H A is obtained by tracing out the degrees of freedom of H B . The idea is then to compute the reduced density matrix of the thermal density operatorρ defined in Equation (7), and then to use the new density operator for computing the bipartite residual entropy at finite T. Although the partition of the Hilbert space is obviously independent from the choice of the basis in which the density operator is represented, to perform the calculation described above is convenient express the density operator Equation (7) in an alternative suitable basis for the partition A-B one has chosen. From the practical point of view, a suitable choice of the basis can give easy access to a partition that in another basis would be really hard to handle computationally. An example of this is shown in Section 4 when we consider the partition between the momentum modes.
Let's expand density operator Equation (7) in a convenient basis { φ i ⟩} for the choice of the partition. To do so, we expand the energy eigenstate substitute it in expression Equation (7), and obtain the new expression for the density operator [37] ρ ≡ρ(T) = where Notice that coefficients C i,j,i ′ ,j ′ (T) contain both thermal and quantum information as they are obtained by thermal-averaging the quantum amplitudes w i,j,n w * i ′ ,j ′ ,n of each energy eigenstate Ψ n ⟩. By tracing over the degrees of freedom of H B is possible to derive the expression of the reduced density operatorρ A (T)ρ The bipartite residual entropy at finite temperature S R (T) is then defined as The details of this calculation, together with the results of the computation of Equations (12) and (13) for different choices of the partition, are discussed in Section 4.

Bipartite Residual Entropy at Zero and Finite Temperature
As already mentioned, "bipartite entanglement" is well defined when the way to partition the system with respect to a certain physical property is specified. Investigating specific properties of a given system leads to consider specific kinds of entanglement. An effective and standard way to quantify the bipartite residual entropy is to compute the Von Neumann entropy according to the scheme discussed in the previous Section. Of course, once the partition is fixed, the computation of the Von Neumann entropy relevant to the reduced density matrix (bipartite residual entropy) is independent on the basis chosen to represent physical states, namely, S(ρ) = S(UρU † ) for any unitary transformation U that enacts the change of basis.
In the sequel, we consider three different kinds of bipartite residual entropy, each one associated with a different way of partitioning the system. First, we consider the quite natural partition in terms of left-well bosons and right-well bosons suggested by the representation of physical states in the space-mode Fock basis Equation (2). Then, by representing physical states in the momentum-mode Fock basis, we partition the system in terms of zero-momentum and non-zero-momentum bosons. Finally, we consider the partition of the system distinguishing species-A from species-B bosons, which is again suggested by definition Equation (2) where populations n L , n R and m L , m R refer to species A and B, respectively. In all three cases, we present the results, obtained numerically, both for the zero-temperature scenario, when only the ground state ψ 0 ⟩ is involved and for the finite-temperature configuration, when the system is naturally described by means of a thermal density matrix. It is worth remarking that at T = 0, the bipartite residual entropy reduces to entanglement entropy because classical correlations are suppressed.

Bipartite Residual Entropy for a Partition Characterized by Spatial Modes
Let us start by computing the bipartite residual entropy S R by considering the partition of the TDS in terms of left-well bosons and the right-well bosons. Following Formula Equation (3) The reduced density matrix relevant to the right-well bosons, obtained tracing out the degrees of freedom of the left-well bosons, iŝ In the presence of a non-zero temperature, the density matrix modifies taking into account the contributions of the whole energy spectrum. By following the scheme discussed in Section 3.2, as the T ≠ 0 density matrix is diagonal, one can easily compute the bipartite residual entropy Equation (13) finding with C i,j (T) given by Formula Equation (11). Figure 3 shows how the bipartite residual entropy relevant to right-well bosons varies with respect to W J, for different temperatures. At T = 0, the plot of S R (black dashed line), which represents the entanglement entropy, exhibits two sharp peaks where the DL transitions occur. In the region between the two peaks bosons are delocalized and the quantum fluids fully mixed, the left tail corresponds to supermixed states (states where both species are localized in a single well) and, eventually, the right tail is the region where the two species localize in different wells. Both tails feature a genuinely quantum behavior because the relevant ground states correspond to Schrödinger cats, in which the spatial separation gets more and more pronounced as W J increases (see Formula Equation (4)). In fact, the entanglement entropy asymptotically tends to 1, a value which is reminiscent of the double-edged structure of cat states Equation (4) because both their components contribute to Formula Equation (15) with 1 2 log 2 2. It is worth noticing that, at T = 0, S R is always different from zero in that, even for noninteracting species (W = 0), the presence of a non-zero J couples the left and right modes of either species. As expected, one can show numerically that the height of the central minimum of S R decreases more and more (tending to zero) as the interwell hopping J becomes smaller and smaller.
At temperature T > 0, Figure 3 shows that the bipartite residual entropy is still able to highlight the difference among mixed, demixed and supermixed phases thanks to the presence of the two peaks. In addition to the spectral collapse mechanism, the emergence of the two peaks in the proximity of the W J critical values is caused by the bigger and bigger number of Fock states significantly contributing to increase S R (the same interpretation holds for the two peaks characterizing S R with A-B partition). This effect is suggested by the fact that the two Gaussian-like distributions (see Figure 1b) involve the maximum number of Fock states when they are about to merge and the transitions take place. The effect of a finite temperature is to smooth the DL phase transitions, an effect that can be clearly appreciated observing the decreasing sharpness of the peaks as T is increased. Interestingly, all the tails of the plotted curves tend to the limiting value 1. For example, in the left tail (W J < 0), this means that, in S R one has C 0,0 2 = C N,M 2 ≃ 1 2 while all the other C i,j 2 are vanishingly small. The resulting S R = 1 follows from the fact that there exist two dominating macroscopic configurations N, M⟩ L 0, 0⟩ R and 0, 0⟩ L N, M⟩ R whose correlation is mainly due to quantum entanglement for T → 0 but assumes a more and more classical character for higher temperatures. In the case of the right tail (W J > 0), the same effect is observed, but the dominating components are C N,0 2 and C 0,M 2 . Note that, at fixed temperature, such configurations emerge provided that the interspecies interaction W is strong enough to contrast the temperature-induced disorder. Of course, for a given value of W J, one has larger residual entropies at higher temperatures in that increasing T makes more and more energy eigenstates accessible in the thermal superposition ensuing from Formula Equation (7). We conclude by observing that, at high temperatures, in the central region around W J = 0, S R approaches the limiting value log 2 D ≈ 10.31 because C i,j 2 ≃ 1 D for all (i, j), where D = 1271 = (N + 1)(M + 1) (with N = 30, M = 40) is the dimension of the Hilbert space. This limiting situation reflects the fact that, at high temperatures, S R → S eq (see Figure 2).

Bipartite Residual Entropy for a Partition Characterized by Momentum Modes
Let us introduce the following momentum-mode operators obtained summing and subtracting usual site-mode operators together with the corresponding number operators The reduced density matrix relevant to the sub-system of bosons having non-vanishing momentum (modes D's) is obtained by tracing out the degrees of freedom relevant to the sub-system of bosons having zero momentum (modes S's) For non-zero temperatures, one must consider the contributions of all the energy levels. Making use of the same scheme discussed in Section 3.2, as the reduced density matrix relevant to the thermal superposition is diagonal, the bipartite residual entropy Equation (13) is found to be C n S ,m S (T) 2 log 2 C n S ,m S (T) 2 . Figure 4 shows the bipartite residual entropy characterizing the separation between still and circulating bosons in respect of the ratio W J, for different temperatures. At T = 0, bipartite residual entropy corresponds to entanglement entropy and its plot (black dashed line) exhibits two sharp discontinuities at the two values of W J for which the DL phase transitions occur. Such discontinuities separate three quasi-plateaus corresponding to supermixed, mixed and demixed phases.
The central region (mixed species) features a quite small entanglement between circulating and still bosons. In fact, if the interspecies coupling W is small compared to the tunneling J and if the ratio U J is small enough to guarantee superfluid and delocalized bosons, momentum modes S a and S b are macroscopically occupied, while D a and D b are poorly populated. If the intraspecies repulsion U tends to zero, one can show that the latter momentum modes are not populated at all, and, at T = 0, the EE vanishes for W J = 0.
At finite temperatures, the behavior of the bipartite residual entropy still mirrors the presence of the three quantum phases. Unlike the behaviors of S R discussed in Section 4.1, where S R = 1 associated to outer plateaus showed that system features two dominating space configurations, here, the value S R = 7.2 implies that, for sufficiently large W J, a much larger number of momentum configurations is involved in determining the system correlations. Figure 4 displays a gap between the plateau S R ≈ 6.2 obtained at T = 0 (black dashed lines) and the limiting value S R ≈ 7.2 of the plateaus obtained at T ≠ 0 (colored lines). This is due to the fact that, in the tails, the energy gap between the ground state and the first excited level becomes vanishingly small but remains non-zero and so the lowest non-zero temperature T = 0.1J k B considered in Figure 4 is already enough to populate both the ground state and the first excited level. The activation of the excited level (absent at T = 0) is sufficient to redistribute the boson population, thus causing the jump of S R from 6.2 to 7.2. As noticed for the partition in terms of spatial modes discussed in Section 4.1, (i) the maximum value of S R tends to the extreme value log 2 1271 ≈ 10.31 at high temperature; and (ii) given a certain value of W J, the bipartite residual entropy steadily increases with temperature T because more and more energy eigenstates become statistically accessible.

Bipartite Residual Entropy for a Partition Characterized by Boson Species
A third way to compute the bipartite residual entropy consists in partitioning the system in terms of species-A and species-B bosons. We use the representation in terms of space-mode Fock states, although the momentum-mode Fock basis is equally convenient to the job. Starting from density matrix Equation (14), the reduced density matrix relevant to species-B sub-system is obtained by tracing out the degrees of freedom relevant to species-A sub-system where we have defined The diagonalization ofρ B provides the eigenvalues {λ j } necessary to compute the relevant Von Neumann entropy Figure 5 shows the bipartite residual entropy relevant to species-mode partition scheme as a function of W J, for different temperatures.
As in Figure 3, at zero temperature (black dashed line), two sharp peaks, at which the DL transitions occur, separate the three regions corresponding to the supermixed, mixed and demixed phase. Also in the present case, the outer regions consist of two quasi-plateaus whose height quickly converges to 1, a limiting value which is, once again, reminiscent of the two-component character of cat states Equation (4) (recall that 1 = 2 × − 1 2 log 2 1 2 ). As noted in the previous subsections, one can show that the zero-temperature EE relevant to the space-mode and the momentum-mode separation schemes features a central minimum tending to zero for J → 0 and U → 0, respectively. In the current case, where the species-mode separation is adopted, the vanishing of the minimum of S R is obtained when the two species are non interacting, namely, for W = 0.
When the temperature is switched on, the DL phase transitions become less abrupt and the corresponding peaks in the plots are less sharp. However, as shown in Figure 5, S R still represents an effective indicator of the critical behavior in a non-small temperature range. As for the S R analyzed in Section 4.1, the residual-entropy plot at non-zero temperatures shows that S R → 1 for W J large enough. Once more, the limiting value S R = 1 = log 2 2 (which all colored lines of Figure 5 converge to) highlights how the system features two equiprobable dominating configurations for large interactions. A non-vanishing T disturbs the formation of such configurations since, in the tails, for a given value of W J, the higher the temperature, the more S R differs from S R = 1. As for the other partition schemes, for large T, S R tends to a maximum value, log 2 D B ≃ 5.36, where D B = (M + 1) is the dimension of sub-system-B Hilbert space.

Bipartite Residual Entropy at Zero Temperature
As repeatedly stressed in the previous discussion, in principle, the choice of the partition, is completely arbitrary and independent on the system under examination. It has more to do with the concepts of "observer" and "measure" than with the physical system itself, opening interesting questions on the relation between entropy and quantum information. To emphsize this fact, in Figure 6, we compare the bipartite residual entropy at T = 0 for the three partition schemes considered above and shows how the presence of a non-zero bipartite residual entropy (i.e., of the EE in the ground-state) strongly depends on the choice of the partition. In particular, we notice how a strong entanglement in a partition can result in a weak (or zero) entanglement in another one. This is the case, e.g., of W J = 0 in which the ground-state is strongly entangled if measured through the partition L-R (finite bipartite residual entropy S R ), or completely disentangled if measured through the partition A-B (bipartite residual entropy S R = 0). In other words, in the same physical system, while the knowledge of the state of the system in the left (right) well is strongly correlated with the information on the state in the right (left) well, on the opposite, the knowledge of the species-A state does not produce information on the species-B state.

Calculation of the EE in the Coherent-State Picture
The coherent-state variational approach has found large application in the study of many-body quantum systems [38] since, due to their semi-classical character, they provide an effective description of physical systems and allow one to gain insights into their properties. In addition, from the experimental point of view, coherent states have an important role since their semi-classical character enables one to achieve a realistic approximation of the quantum state describing the real system.
An su(2) coherent state describing single condensate trapped in a dimer is given by [21] ξ L , where 0⟩ = 0, 0⟩ is the boson vacuum state and the normalization condition ξ L 2 + ξ R 2 = 1 must be assumed. Since ⟨ψ a A † σ A σ ψ a ⟩ = N ξ σ 2 , with σ = R, L, is the expectation value of number operator N σ = A † σ A σ then ξ σ 2 represents the fraction of bosons in the well σ. In the following, we employ combinations of coherent states Equation (16) (for a single species in a double well) to approximate the cat structure of the ground state relevant to the TSD system in the strong-interaction regime, both for W J > 0 and for W J < 0.

1.
Supermixing (attractive cat). If the interspecies attraction (W J < 0) is large enough, the two species aggregate together in the same well. Since none of the two wells is privileged with respect to the other, quantum mechanically both configurations are equally probable, and the system lives in both states at the same time. By using the notation of Formula Equation (16), the resulting cat state can be written as where "Loc" stands for "localized" and entails the fact that η c 2 ≪ λ c 2 . Following the scheme discussed in Ref. [39], one can show that the expectation value of the model Hamiltonian reduces to where the local order parameters λ a , λ b , η a , and η b are complex quantities defined as The minimum-energy configuration is reached for φ a = θ a , θ b = φ b and These formulas give the fraction of bosons characterizing the minority component and, correctly, give zero in the limit W → −∞.

2.
Demixing (repulsive cat). If the interspecies repulsion (W J > 0) is large enough, the two condensed species separate in different wells. Similarly to what explained in the previous paragraph, the ground state features a two-sided cat-like structure because the left (right) well can indistinctly host species A (B). Hence, the quantum state consists of an equally-weighted superposition of the two possible arrangements where λ c , η c are such that η c 2 ≪ λ c 2 and (obviously) λ c 2 + η c 2 = 1, with c = a, b. Following the variational approach described in the previous paragraph, and adopting the same conventions, we obtain that the variational energy is minimized for θ a = φ a , θ b = φ b and Parameters x a and x b represent the fractions of bosons, which do not aggregate with the others and thus make the "demixed phase" not ideal. Notice that, correctly, if W → +∞, then x a,b → 0, i.e., the demixing gets more and more complete.
Both for the supermixing and for the demixing scenario, after computing the fraction of bosons in each well, it is possible to reconstruct the cat state by superimposing two coherent states. This procedure, described in Appendix A, allows one to analytically compute the EE between left-well and right-well bosons, at zero temperature [39]. As shown in Figure 7, the result perfectly matches the numerical EE, of course in the validity range of this approximation, i.e., in the whole range of W J except the central region (mixed phase) between the two critical values.

Calculation of the Bipartite Residual Entropy in a Restricted Energy Basis
As already explained, the density operator associated to a thermal mixture of eigenstates iŝ where E n is the energy eigenvalue associated to the energy eigenstate ψ n ⟩, β is (proportional to) the inverse temperature and D is the dimension of the Hilbert space of physical states. From a computational, but also from a conceptual point of view,ρ is the superposition of D different contributions, each one weighted by a different Boltzmann factor. The dimension D rapidly increases with the number of particles hosted in the system its exact value being D = (N + 1)(M + 1). As a consequence, the computation of the thermal density matrix becomes unfeasible even for a relative small number of bosons. By taking advantage of the well-known equivalence between microcanonical and canonical ensemble (see, e.g., [37]), for large numbers of particles, we provide an effective way to approximate a thermal state. To this end, we consider just a restricted set of energy eigenstates, namely those ψ n ⟩ whose energy E n lies in the range are the expectation value of the energy and its standard deviation, respectively. The density matrix relevant to this restricted thermal state is thus constructed by equally-weighting the contributions coming from such ψ n ⟩, i.e.,ρ where N * is the number of energy eigenstates whose energies E n lie in the aforementioned interval.
To test the effectiveness of the bipartite residual entropy as a critical indicator, we consider the partition in terms of left-well bosons and right-well bosons, we set a non-zero value of the temperature and we compare the results obtained from a complete and from a restricted thermal state. The left panel of Figure 8 shows an overall good agreement between such results, especially in the central region (small W J values), while the outermost regions feature step-like discontinuities. The presence of such discontinuities can be understood observing ethe right panel of Figure 8, which shows the fraction of energy states involved in the restricted thermal state, N * D, as a function of W J. As W J increases, in fact, fewer and fewer energy states join the restricted thermal state and their inherently discrete character is reflected by the presence of step-like regions, each one corresponding to the activation of a single energy state.

Conclusions
In this work, we have investigated the equilibrium and the bipartite residual entropy in a two-species Bose-Hubbard dimer at zero and non-zero temperature. In Section 2, we have introduced the model and highlighted the importance of W (the interspecies repulsion) in determining the quantum phase of the system (supermixed for W J ≪ 0, mixed for small W J and demixed for W J ≫ 0). In Section 3, we have introduced the concepts of equilibrium and bipartite residual entropy commenting on the fact that, at zero temperature, the latter corresponds to the entanglement entropy. Section 4 has been devoted to the analysis of the bipartite residual entropy for three different partitions of the total system. In this regard, we have stressed the fact that different ways of partitioning the system into two sub-systems, correspond to different kinds of residual entropies S R . In all three cases, S R features discontinuities where the localization-delocalization phase transitions occur and quasi-plateaus where two dominating macroscopic configurations emerge. Bipartite residual entropy is therefore a valid critical indicator not only at zero temperature (where it corresponds to the entanglement entropy, a purely quantum correlation), but also at higher temperatures, where it is influenced by the classical correlation between the sub-systems and by classical entropy. Interestingly, we have made evidence that, at zero temperature, (i) a non-zero hopping J causes a non-zero entanglement between spatial modes; (ii) the intraspecies interaction U contributes to the entanglement between momentum modes; and (iii) the interspecies interaction W is responsible for the entanglement between species modes.
In Section 5, we have introduced su(2) coherent states and developed a fully-analytic variational approach apt to describe the supermixed and the demixed phases at zero temperature. The superposition of two such coherent states has provided a good approximation of the ground state of the system in a non-small range of W J, as demonstrated by the comparison with the numerical results. In Section 6, we have approximated the complete thermal superposition Equation (10) with an incoherent combination of a reduced number of equally-weighted energy eigenstates and showed that the bipartite residual entropy is still a good critical indicator, well reproducing the exact results obtained numerically.
On the basis of the coherent-state approach derived in Section 5 and in the same spirit of Ref. [39], we compute the entanglement entropy between left-well bosons and right-well bosons. To begin, let us define ρ n,m (i) as the probability of having n bosons of species A and m bosons of species B at site i. The normalization of probability requires that Neglecting the possible presence of cat states (a situation that will be re-inserted a posteriori), a generic coherent state can be written in the factorized form Of course, the normalization conditions ξ L 2 + ξ R 2 = 1, ν L 2 + ν R 2 = 1 must hold. State Ψ⟩ can be recast into the form We calculate the reduced density matrix ρ partitioning the system into two sub-systems (left-well bosons and right-well bosons) and tracing out the degrees of freedom relevant to one of them. Taking into account the orthogonality of the states, the reduced density matrix that originates from a coherent state can be written as where the expressions of coefficients ξ L = ξ L (J, U, W), ξ R = ξ R (J, U, W), ν L = ν L (T, U, W) and ν R = ν R (T, U, W) can be computed within the variational approach. In passing, notice that the probability distribution is correctly normalized, i.e., ∑ This quite general procedure needs to be slightly modified in case one is considering cat states. In fact, the reduced density matrix must take into account the two-sided nature of a cat state and so it must be written as the average of the densities matrices relevant to simple coherent states, namely implying ρ cat n,m = 1 2 ρ side L n,m + ρ side R n,m .

Appendix B. Quasi-Degeneracy of the Ground-State
Due to the spectral-collapse [16,17], for sufficiently strong values of W J, the TSD ground-state may appear quasi-degenerate. However, the degeneracy of the ground-state is only apparent as Hamiltonian Equation (1) is non-degenerate [34,36]. As shown in Figure A1, the energy splitting between the ground-state energy E 0 and the first energy level E 1 decays exponentially as a function of the number of particles per species (N and M). Energy levels E 0 and E 1 differs always by a small, finite, quantity function of the interactions and the number of particles. This has been verified in Figure A1 down to computational limit fixed by the machine precision.