Magnetism and domain formation in SU(3)-symmetric multi-species Fermi mixtures

We study the phase diagram of an SU(3)-symmetric mixture of three-component ultracold fermions with attractive interactions in an optical lattice, including the additional effect on the mixture of an effective three-body constraint induced by three-body losses. We address the properties of the system in $D \geq 2$ by using dynamical mean-field theory and variational Monte Carlo techniques. The phase diagram of the model shows a strong interplay between magnetism and superfluidity. In the absence of the three-body constraint (no losses), the system undergoes a phase transition from a color superfluid phase to a trionic phase, which shows additional particle density modulations at half-filling. Away from the particle-hole symmetric point the color superfluid phase is always spontaneously magnetized, leading to the formation of different color superfluid domains in systems where the total number of particles of each species is conserved. This can be seen as the SU(3) symmetric realization of a more general tendency to phase-separation in three-component Fermi mixtures. The three-body constraint strongly disfavors the trionic phase, stabilizing a (fully magnetized) color superfluid also at strong coupling. With increasing temperature we observe a transition to a non-magnetized SU(3) Fermi liquid phase.


Introduction
Cold atoms in optical lattices provide us with an excellent tool to investigate notoriously difficult problems in condensed matter physics [1,2]. Recent progress towards this goal is exemplified by the experimental observation of the fermionic Mott insulator [3,4] in a binary mixture of repulsively interacting 40 K atoms loaded into an optical lattice, and of the crossover between Bardeen-Cooper-Schrieffer (BCS) superfluidity and Bose-Einstein condensation (BEC) [5,6,7] in a mixture of 6 Li atoms with attractive interactions.
At the same time, ultracold quantum gases also allow us to investigate systems which have no immediate counterparts in condensed matter. This is the case for fermionic mixtures where three internal states σ = 1, 2, 3 are used, instead of the usual binary mixtures that mimic the electronic spin σ =↑, ↓. These multi-species Fermi mixtures are already available in the laboratory, where three different magnetic sublevels of 6 Li [8,9,10,11] or 173 Yb [12], as well as a mixture of the two internal states of 6 Li with a lowest hyperfine state of 40 K [13] have been successfully trapped. In the case of Alkali atoms, magnetic or optical Fano-Feshbach resonances can be used to tune magnitude and sign of the interactions in the system, and in the case of Ytterbium or group II atoms, it is possible to realise three-component mixtures where the components differ only by nuclear spin, and therefore exhibit SU(3) symmetric interactions [14,15,16]. Moreover, loading these mixtures into an optical lattice would give experimental access to intriguing physical scenarios, since they can realize a three-species Hubbard model with a high degree of control of the Hamiltonian parameters.
Multi-species Hubbard models have attracted considerable interest on the theoretical side in recent years. First studies were focused on the SU (3)-symmetric version of the model with attractive interaction. By using a generalized BCS approach [17,18,19], it was shown that the ground state at weak-coupling spontaneously breaks the SU (3) ⊗ U (1) symmetry down to SU (2) ⊗ U (1), giving rise to a color superfluid (c-SF) phase, where superfluid pairs coexist with unpaired fermions. Within a variational Gutzwiller technique [20,21] the superfluid phase was then found to undergo for increasing attraction a phase transition to a Fermi liquid trionic phase, where bound states (trions) of the three different species are formed and the SU (3)-symmetry is restored. More recently [22,23], the same scenario has been found by using a self-energy functional approach for the half-filled model on a Bethe lattice in dimension D = ∞. It was suggested [24] that this transition bears analogies to the transition between quark superfluid and baryonic phase in the context of Quantum Chromo Dynamics.
Both the attractive and the repulsive version of the model was addressed by numerical and analytical techniques for the peculiar case of spatial dimension D = 1 [25,26,27,28], while Mott physics and instabilities towards (colored) density wave formation have been found in the repulsive case in higher dimensions [17,29,30]. It is important to mention that substantial differences are expected in the attractive case at strong coupling when the lattice is not present [31,32]. Those differences are essentially related to the influence of the lattice in the strong coupling limit in the three-body problem, favoring trion formation [33,34] with respect to pair formation in the continuum, as was shown in Ref. [32,35,36].
Here we consider the SU (3)-symmetric system in a lattice for D ≥ 2 in the presence of attractive two-body interactions by combining dynamical mean-field theory (DMFT) and variational Monte Carlo (VMC). We analyze several cases of interest for commensurate and incommensurate density. Ground state, spectral, and finite temperature properties are addressed. More specifically we focus on the transition between color superfluid and trionic phase and on a better understanding of the coexistence of magnetism and superfluidity in the color superfluid phase already predicted in the SU (3) symmetric case [20,21] but also when the SU (3)-symmetry is explicitly broken [37]. We show that the existence of a spontaneous magnetization leads the system to separate in color superfluid domains with different realizations of color pairing and magnetizations whenever the total number of particles in each hyperfine state is conserved. This would represent a special case, due to the underlying SU (3) symmetry, of a more general tendency towards phase separation in three-component Fermi mixtures. We point out that all this rich and interesting physics arises merely from having three components instead of two. Indeed the analogous SU (2) system would give rise to the more conventional BCS-BEC crossover, where the superfluid ground state evolves continuously for increasing attraction [38]. Moreover in the SU (2) case superfluidity directly competes with magnetism [39].
The case under investigation can be realized with ultracold gases by loading a three-species mixture of 173 Yb [12] or another group II element such as 87 Sr into an optical lattice, or alternatively using 6 Li in a large magnetic field. However, some realizations with ultracold atoms are plagued by three-body losses due to three-body Efimov resonances [8,9,11], which are not any more Pauli suppressed as in the twospecies case. The three-body loss properties and their dependence on the magnetic field have been already measured for 6 Li [8,9,11], while they are still unknown for three-component mixtures of certain group-II elements. Loading a gas into an optical lattice could be used to suppress losses, as a large rate of onsite three-body loss can prevent coherent tunneling processes from populating any site with three particles [40]. As proposed in Ref. [40] for bosonic systems, in the strong loss regime a Hamiltonian formulation is still possible if one includes an effective hard-core three-body interaction, which leads to new interesting physics [41]. The effect of this dynamically generated constraint on the fermionic system in D = 1 with attractive interactions was studied in Ref. [27], where it was shown that the constraint may help to stabilize the superfluid phase in some regions of the phase diagram.
For these reasons we also study the effect of including a three-body constraint in the model, as representative of an SU (3) symmetric mixture in the strong-loss regime. The asymmetric case in the strong loss regime, which is directly relevant for experiments on 6 Li close to a Feshbach resonance, has been already addressed in a separate publication [42].
The paper is organized as follows: in the following sections we first introduce the model (Sec. 2) and then the methods used (Sec. 3). Later on we present our results, focusing first on the unconstrained system (Sec. 4), for commensurate and incommensurate densities and then on the effects of the three-body constraint (Sec. 5).
The emergence of domain formation within globally balanced mixtures is discussed in detail in Sec. 6. Final remarks are drawn in Section 7.

Model
Three-component Fermi mixtures with attractive two-body interactions loaded into an optical lattice are well described by the following Hamiltonian where σ = 1, 2, 3 denotes the different components, J is the hopping parameter between nearest neighbor sites i, j , µ σ is the chemical potential for the species σ and U σσ < 0. We introduced the onsite density operators n iσ = c † iσ c iσ . The three-body interaction term with V = ∞ is introduced to take the effects of three-body losses in the strong loss regime into account according to Refs. [27,40]. V = 0 corresponds to the case when three-body losses are negligible. While the model and the methods are developed for the general case without SU (3)-symmetry, in this paper we concentrate on the SU (3)symmetric case reflected by species-independent parameters U σσ = U, µ σ = µ. ( In this case the Hamiltonian (1) reduces to an SU (3) attractive Hubbard model if V = 0. Note that the three-body interaction term is a color singlet and thus does not break SU (3) for any choice of V . On the basis of previous works, the ground state of the unconstrained model is expected to be, at least in the weak coupling regime, a color superfluid, i.e. a phase where the full SU (3) ⊗ U (1) symmetry of the Hamiltonian is spontaneously broken to SU (2)⊗U (1) [17,18]. As shown in [17,18], it is always possible to find a suitable gauge transformation such that pairing takes place only between two of the natural species σ, σ and in this paper we choose a gauge in which pairing takes place between the species σ = 1 and σ = 2 (1 − 2 channel), while the third species stays unpaired. Whenever the SU (3)-symmetry is explicitly broken, only the pairing between the natural species is allowed to comply with Ward-Takahashi identities [37]. This reduces the continuum set of equivalent pairing channels of the symmetric model to a discrete set of three (mutually exclusive) options for pairing, i.e. 1 − 2, 1 − 3 or 2 − 3. In this case the natural choice would be that pairing takes place in the channel corresponding to the strongest coupling when the mixture is globally balanced. We can always relabel the species such that strongest attractive channel is the channel 1 − 2.
Other pairing channels can be studied via index permutations of the species. Therefore the formalism developed here is fully general and includes both the symmetric and nonsymmetric case, while only in the SU (3)-symmetric case our approach corresponds to a specific choice of the gauge.

Methods
In order to investigate the model in Eq. (1) in spatial dimensions D ≥ 2 we use a combination of numerical techniques which have proven to give very consistent results for the non-symmetric case [42]. In particular, we use dynamical mean-field theory (DMFT) for D ≥ 3 and variational Monte Carlo (VMC) for D = 2. DMFT provides us with the exact solution in infinite dimension and a powerful (and non-perturbative) approach in D = 3, which has the advantage of being directly implemented in the thermodynamic limit (without finite size effects). VMC allows us to incorporate also the effect of spatial fluctuations which are not included within DMFT, even though the exponential growth of the Hilbert space limits the system sizes that are accessible.

DMFT
Dynamical mean-field theory (DMFT) is a non-perturbative technique based on the original idea of Metzner and Vollhardt who studied the limit of infinite dimension of the Hubbard model [43]. In this limit, the self-energy Σ(k, ω) becomes momentum independent Σ(k, ω) = Σ(ω), while fully retaining its frequency dependence. Therefore the many-body problem simplifies significantly, without becoming trivial, and can be solved exactly. In this sense DMFT is a quantum version of the static mean-field theory for classical systems, since it becomes exact in the same limiting case (D = ∞) and can provide useful information also outside of this limit, fully including local quantum fluctuations. In 3D, assuming a momentum independent self-energy, has proved to be a very accurate approximation for many problems where the momentum dependence is not crucial to describe the physics of the system such as the Mott metal-insulator transition [44] where the frequency dependence is more relevant than the k dependence.
3.1.1. Theoretical setup for SU (3) model with spontaneous symmetry breaking -In this work, we generalize the DMFT approach to multi-species Fermi mixtures in order to describe color superfluid and trionic phases, which are the expected phases occurring in the system. The theory can be formulated in terms of a set of self-consistency equations for the components of the local single-particle Green functionĜ on the lattice. Since we are dealing here with superfluid phases involving also anomalous components of the Green function, we use a compact notation in terms of mixed Nambu spinors ψ = (c 1 , c † 2 , c 3 ), where we already assumed that pairing takes place only between the first two species, as explained in the previous section, and we omit the subscript i (spatially homogeneous solution). We reiterate that this specific choice is valid without loss of generality in the SU (3)-symmetric model, and has the same status as fixing the phase of a complex condensate order parameter in theories with global phase symmetry.
In practice the original lattice model (1) in the DMFT approach can be mapped, by introducing auxiliary fermionic degrees of freedom a † lσ , a lσ , on a Single Impurity Anderson Model (SIAM), whose Hamiltonian reads where the Anderson parameters ε lσ , V lσ , W l have to be determined self-consistently. Selfconsistency ensures that the impurity Green function of the SIAM is identical to the local component of the lattice Green function. The components of the non-interacting Green function for the impurity site, which represent the dynamical analog of the Weiss field in classical statistical mechanics, can be expressed in terms of the Anderson parameters as where ζ l,σ = −iω n + ε lσ . The self-consistency equations for the local Green functions now have the form where M is the number of lattice sites,Ĝ latt (k, iω n ) =Ĝ latt (ε k , iω n ) is the lattice Green function within DMFT and D(ε) is the density of states of the lattice under consideration. The independent components ofĜ latt (k, iω n ) have the form where ζ σ = iω n + µ σ − Σ σ (iω n ) and the self-energy can be obtained by the following local Dyson equationΣ(iω n ) =Ĝ −1 Once a self-consistent solution has been obtained, the impurity site of the SIAM represents a generic site of the lattice model under investigation. Therefore several static thermodynamic quantities can be directly evaluated as quantum averages of the impurity site. As evident from the previous equations, DMFT is explicitly formulated in a grand canonical approach where the chemical potentials µ σ are given as input and the onsite densities n σ = c † σ c σ are calculated.

Calculated observables and numerical implementation -
To characterize the different phases, we evaluated several static observables such as the superfluid (SF) order parameter P = c 1 c 2 , the average double occupancy d σσ = n σ n σ and the average triple occupancy t = n 1 n 2 n 3 . As suggested in Refs. [17,18,37], in order to gain condensation energy in the c-SF phase, it is energetically favorable to induce a finite density imbalance between the paired species (1 − 2 in our gauge) and the unpaired fermions. To quantitatively characterize this feature we introduce the local magnetization m = n 12 − n 3 where n 12 = n 1 = n 2 .
From the normal components of the lattice Green functions in Eqs. (10), (11) and (13) we can extract the DMFT momentum distribution and the average of kinetic energy per lattice site It is evident from the expression ofĜ latt (k, iω n ) given in Eqs. (10), (11) and (13) that n σ (k) only depends on the momentum k through the free-particle dispersion ε k of the lattice at hand. The internal energy per lattice site E can then be obtained as is the average potential energy per lattice site. Solving the DMFT equations is equivalent to solving a SIAM in presence of a bath determined self-consistently. We use Exact Diagonalization (ED) [45], which amounts to truncating the number of auxiliary degrees of freedom a † lσ , a lσ in the Anderson model to a finite (and small) number N s − 1. In this way the size of the Hilbert space of the SIAM is manageable and we can exactly solve the Anderson model numerically. Here we would like to point out that this truncation does not reflect the size of the physical lattice but only the number of independent parameters used in the description of the local dynamics. Therefore we always describe the system in the thermodynamic limit (no finite-size effects). We use the Lanczos algorithm [46] to study the ground state properties (up to N s = 7) and full ED for finite temperature (up to N s = 5). Due to the increasing size of the Hilbert space (σ = 1, 2, 3 instead of σ =↑, ↓) in the multicomponent case the typical values of N s which can be handled sensibly is smaller than the corresponding values for the SU (2) superfluid case. However, in thermodynamic quantities, we found indeed only a very weak dependence on the value of N s and the results within full ED at the lowest temperatures are in close agreement with T = 0 calculations within Lanczos.
A definite advantage of ED is that it allows us to directly calculate dynamical observables for real frequencies without need of analytical continuation from imaginary time. In particular, we can directly extract the local single-particle Green function G σ (ω) and the single-particle spectral function

Variational Monte-Carlo
The variational Monte Carlo (VMC) techniques described in this subsection can be used to calculate the energies and correlation functions of the homogeneous phases at T = 0 in a canonical framework. The basic ingredients of the VMC formalism are the Hamiltonian and trial wavefunctions with an appropriate symmetry. In principle, the formalism presented here can be applied to any dimension, even though here we use it specifically to address the system on a two-dimensional square lattice. The canonical version of Hamiltonian (1) for three-components fermions with generic attractive interactions is given by where the three-body constraint is imposed by using the projector P 3 = i (1 − n i,1 n i,2 n i,3 ) and in the unconstrained case we set P 3 equal to the identity.
Practical limitations do not permit a general trial wave function equally accurate both for the weak-and the strong-coupling limit. Due to this reason we introduce different trial wavefunctions for different coupling regimes.
In the weakly interacting limit, which we operatively define as |U σσ | ≤ 4J = W/2, we use the full Hamiltonian (20) along with the weak-coupling trial wavefunction defined in the next subsection. Here W = 2DJ is the bandwidth. At strong-coupling this wavefunction results in a poor description of the system. In order to gain insight into the strong coupling regime, we derive below a perturbative Hamiltonian to the second order in J/U σσ , which we will combine with a strong-coupling trial wavefunction. Again the strong-coupling wavefunctions are incompatible with the Hamiltonian (20), as will be clarified below. We can therefore address confidently both limits of the model while at intermediate coupling we expect our VMC results to be less accurate.

3.2.1.
Strong coupling Hamiltonian, constrained case -In order to derive a perturbative strong-coupling Hamiltonian for the constrained case we make use of the Wolff-Schrieffer transformation [47] H pert = P D e iS He −iS P D (21) and keep terms up to the second order in J/U σσ . In the expression above, P D is the projection operator to the Hilbert subspace with fixed numbers of double occupancies in each channel (N 12 , and e iS is a unitary transformation defined in Appendix A. So we obtain the perturbative Hamiltonian (see Appendix A), which reads: For the case where the SU (3)-symmetry is restored (U σσ = U ), the perturbative Hamiltonian can be written in a compact notation where the double occupancy operator is now defined as . Now, rather than conserving the number of double occupancies N σσ d in each channel, only the total number N d,0 = N 12 is conserved due to the SU (3)-symmetry. Indeed Eq. (23) contains terms where the tightly bound dimers are allowed to change the composition through second order processes. Thus, the SU (3)-symmetric case, in contrast to the case with strongly anisotropic interactions is qualitatively different from the Bose-Fermi mixture, because the bosons -tightly bound dimers -can change composition as described above, while such a process was not allowed in the case of the strong anisotropic interactions. We also notice that the last of the ∼ J 2 /U terms contributes only when N d,0 < N/2.

Strong coupling
Hamiltonian, unconstrained case -Without the 3-body constraint three fermions with different hyperfine states can occupy the same lattice site and we expect them to form trionic bound states at sufficiently strong coupling. Correspondingly the many-body system should be in a trionic phase with heavy trionic quasiparticles, as mentioned in previous studies [17,18,22,23]. Therefore we expect that our perturbative approach can provide a description of the trions in the strong coupling limit.
First we consider the extreme case J = 0. In this limit formation of local trions takes place, i.e each site is either empty or occupied by three fermions with different hyperfine spins. Their spatial distribution is random, because any distribution of trions will have the same energy. For finite J with J |U σ,σ | the hopping term can break a local trion, but this would result in a large energy penalty.
According to perturbation theory up to third order we could have two different contributions: (i) one of the fermions hops to one of the neighboring sites and returns back to the original site (second order perturbation), (ii) all three fermions hop to the same nearest neighbor site (third order perturbation). As we show below, due to the first process there is an effective interaction between trions on nearest neighbor sites. Also due to this process the onsite energy has to be renormalized. The second process (ii) describes the hopping of a local trions to a neighboring site.
After straightfroward calculations (see Appendix A) we obtain that the effective interaction between two trions on neighboring sites is For the SU (3)-symmetric case this expression is simplified and we obtain Therefore the nearest-neighbour interaction between trions is repulsive in the SU (3)symmetric case.
For the hopping coefficient we obtain where σ, σ and σ are different from each other in the sum. In the SU (3)-symmetric case, the expression again simplifies to So we obtain the following effective Hamiltonian [48] Here t † i is the creation operator of a local trion at lattice site i and n T i = t † i t i is the trionic number operator. Because the effective hopping of trions results from a third order process and the interaction from second order, more precisely J ef f = J/|U | · V ef f , the effective trion theory is interaction dominated. Since the interaction describes nearest-neighbour repulsion, the strong coupling limit clearly favors a checkerboard charge density wave ground state at half-filling ‡, which we will discuss in more detail in Sec. 4.

Trial wavefunctions:
In order to describe a normal Fermi liquid phase without superfluid pairing, we use the following trial wavefunction where |0 is the vacuum state and ε k,σ = −2J(cos(k x ) + cos(k y )) for a 2D square lattice with only nearest-neighbor hopping. The dependence on the densities is included in the value of the non-interacting Fermi energy ε F,σ . The wavefunction above has no variational parameters except for the choice of Jastrow factor which takes into account the effect of the interaction. Here ν 3 and ν c are variational parameters and i,j is summation with nearest neigbours. The weak-coupling version of the wavefunctions presented in this part is obtained by setting P D equal to unity.
We also consider the broken symmetry SU (2) ⊗ U (1) phase with s-wave pairing in the 1 − 2 channel, whose trial wavefunction is given by where In this case, in addition to the Jastrow factor J , we haveμ and ∆ 0 as additional variational parameters. The s-wave gap function ∆ s (k) = ∆ 0 has no k dependence. This parametrization of ∆ s (k) leads upon Fourier transform to a singlet symmetric pairing orbital φ s (r 1 , r 2 ) = φ s (r 2 , r 1 ).
In practice the optimization parameter ∆ 0 depends on the density n as well as on the coupling strength U . Also, even at the same coupling strength U , the ∆ 0 can be qualitatively different for the weak and the strong coupling ansatz (in the intermediate regime U ≈ −5J). On the other hand, the parameterμ depends mostly on n (and only weakly on U ). The general tendency we observe is that ∆ 0 is suppressed beyond the filling density n 1 in the presence of the constraint. Within a BCS mean-field theory approach, the condensation energy E cond is easily related to the order parameter ∆ 0 , being E cond ∝ ∆ 2 0 . We however calculate it explicitly from the definition by comparing ‡ Despite our fermions are charge neutral, we use sometimes the expression charge density wave in analogy with the terminology commonly used in condensed matter physics.
the ground state energies of the normal and the superfluid phases for the same density. Therefore we define where We also calculate the order parameter P that characterizes the superfluid correlation by considering the long range behavior of the pair correlation function where and M is the total number of the lattice sites. Finally, in order to describe the trionic Fermi liquid phase we can use the following trial wavefunction In this case the Jastrow factor Here ν t is a variational parameter and i,j is summation over nearest neigbours.

Results: SU (3) attractive Hubbard Model
We first consider the SU (3) attractive Hubbard model described by the Hamiltonian (1) with V = 0. In a physical realization with ultracold gases in optical lattices, this corresponds to a situation where three-body losses are negligible. In order to address the effects of dimensionality and of particle-hole symmetry, we analyze several cases of interest, namely (i) an infinite-dimensional Bethe lattice in the commensurate case (halffilling), (ii) a three-dimensional cubic lattice and (iii) a two-dimensional square lattice, the latter two in the incommensurate case. In order to simplify comparison of results on different dimensions, we rescaled everywhere the energies by the bandwidth W of the specific lattice under consideration. For a Bethe lattice in D = ∞ the bandwidth is related to the hopping parameter by W = 4J, while for a D dimensional hypercubic lattice it is W = 4DJ.

Bethe lattice at half-filling
We first consider the infinite dimensional case, for which the DMFT approach provides the exact solution of the many-body problem whenever the symmetry breaking pattern of the system can be correctly anticipated. For technical reasons we consider here the Bethe lattice in D = ∞, which has a well defined semicircular density of states, given by the following expression The simple form of the self-consistency relation for DMFT on the Bethe lattice introduces technical advantages, as explained below. Moreover, we can directly compare our results with recent calculations for the same system within a Self-energy Functional Approach (SFA) [22,23].
In the absence of three-body repulsion, the Hamiltonian (1) is particle holesymmetric whenever we choose µ = U . In this case the system is half-filled, i.e. n σ = 1 2 for all of σ and n = σ n σ = 1.5.
We first consider the ground state properties of the system which we characterize via the static and dynamic observables defined in Sec. 3. For small values of the interaction (|U | W ), we found the system to be in a color-SF phase, i.e. a phase where superfluid pairs coexist with unpaired fermions (species 1-2 and 3 respectively in our gauge) and the superfluid order parameter P (plotted in Fig. 1 using green triangles) is finite. This result is in agreement with previous mean field studies [17,18], as expected since DMFT includes the (static) mean-field approach as a special limit, and with more recent SFA results [22,23]. By increasing the interaction |U | in the c-SF phase, P first increases continuously from a BCS-type exponential behavior at weak-coupling to a non-BCS regime at intermediate coupling where it shows a maximum and then starts decreasing for larger values of |U |. This non-monotonic behavior is beyond reach of a static mean-field approach and agrees perfectly with the SFA result [22,23]. As explained in the introduction, the spontaneous symmetry breaking in the c-SF phase is generally expected [20,21,37,49] to induce a population imbalance between the paired channel and the unpaired fermions, i.e. a finite value of the magnetization m in Eq. (15). It is however worth pointing out that, due to particle-hole symmetry, the c-SF phase at half-filling does not show any induced population imbalance, i.e. m = 0 for all values of the interaction strength. As discussed in the next subsection, the population imbalance is indeed triggered by the condensation energy gain in the paired channel. This energy gain, however, cannot be realized at half-filling where the condensation energy is already maximal for a given U .
Further increasing |U |, we found P to suddenly drop to zero at |U | = U c,2 ≈ 0.45W , signaling a first order transition to a non-superfluid phase. This result is in good quantitative agreement with the SFA result in Ref. [22,23], where a first order transition to a trionic phase was found, while a previous variational calculation found a second order transition [20,21].
In this new phase we were not able to stabilize a homogeneous solution of the DMFT equations with the ED algorithm [45]. Such a spatially homogeneus phase would correspond to having identical solutions, within the required tolerance, at iteration n and n + 1 of the DMFT self-consistency loop. In the normal phase (|U | > U c,2 ) instead, we found a staggered pattern in the solutions and convergence is achieved if one applies (Color online) SF order parameter P (green triangles) and CDW order parameter C (blue squares) plotted as a function of interaction strength |U |/W on the Bethe lattice in the limit D → ∞ at half-filling and T = 0. C > (empty squares) and C < (full squares) correspond to calculations starting from a superfluid or a trionic charge density (t-CDW) initial conditions, and similar for P > (note that P < is always vanishing and therefore not shown). In the inset we compare the ground state energies of the c-SF and t-CDW phases. (Unconstrained, i.e. V = 0) a staggered criterion of convergence by comparing the solutions in iteration n and n + 2. This behavior is clearly signalling that the transition to a non-superfluid phase is accompanied by a spontaneous symmetry breaking of the lattice translational symmetry into two inequivalent sublattices A and B. In a generic lattice a proper description of this phase would require solving two coupled impurity problems, i.e. one for each sublattice, and generalizing the DMFT equations introduced in the previous section. In the Bethe lattice instead the two procedures are equivalent §.
In the new phase the full SU (3)-symmetry of the hamiltonian is restored and we identify it with as a trionic Charge Density Wave (t-CDW) phase. In order to characterize this phase, we introduce a new order parameter which measures the density imbalance with respect to the sublattices A (majority) and B (minority), i.e.
where n A ≡ n σ,A and n b ≡ n σ,B for all σ and C = 0 in the c-SF phase because the translational invariance is preserved. The evolution of the CDW order parameter C in the t-CDW is shown in Fig. 1 using blue squares. At the phase transition from c-SF to t-CDW phase, P goes to zero and C jumps from zero to a finite value. Then C increases further with increasing attraction |U | and eventually saturates at C = 1/2 for |U | → ∞. Motivated by these findings, we considered more carefully the region around the transition point. Surprisingly we found that upon decreasing |U | from strong-to weak-coupling the t-CDW phase survives far below U c,2 down to a lower § On the Bethe lattice the sublattices A and B are completely decoupled from each other at a given step n. critical value U c1 0, revealing the existence of a coexistence region in analogy with the hysteretic behavior found at the Mott transition in the single band Hubbard model [44]. In the present case, however, we did not find any simple argument to understand which phase is stable and had to directly compare the ground state energy of the two phases in the coexistence region to find the actual transition point. In the Bethe lattice, the kinetic energy per lattice site K in the c-SF and t-CDW phases can be expressed directly in terms of the components of the local Green functionĜ(iω n ), which is straightforwardly determined by DMFT. The potential energy per lattice site V is given by , where the index indicates the sublattice. By generalizing analogous expressions valid in the SU (2) case [38,50,51], we obtain and Results shown in the inset of Fig. 1 indicate that the t-CDW phase is stable in a large part of the coexistence region and that the actual phase transition takes place at |U | = U c ≈ 0.2W . The good agreement between our findings and the SFA results in Ref. [22,23] concerning the maximum value of the attraction U c2 where a c-SF phase solution is found within DMFT would suggest that this value is indeed a critical threshold for the existence of a c-SF phase. On the other hand we also proved that the c-SF phase close to U c2 is metastable with respect to the t-CDW phase and therefore the existence of the threshold could equally results from an inability of our DMFT solver to further follow the metastable c-SF phase at strong coupling. The disagreement between our findings and Ref. [22,23] for what concerns the existence of CDW modulations in the trionic phase is clearly due to the constraint of homogeneity imposed in the SFA approach of Ref. [22,23] in order to stabilize a (metastable) trionic Fermi liquid instead of the t-CDW solution. In our case, this was not an issue due to the fact that the iterative procedure of solution immediately reflects the spontaneous symmetry breaking of the translational invariance and does not allow for the stabilization of an (unphysical) homogeneous trionic Fermi liquid at half-filling.
On the other hand, the necessary presence of CDW modulation in the trionic phase at half-filling, at least in the strong-coupling limit, can be easily understood based on general perturbative arguments. Indeed, as pointed out in Sec. 3, in the strongcoupling trionic phase where J/|U | 1, the system can be described in terms of an effective trionic Hamiltonian (28). In this Hamiltonian the effective hopping J ef f of the trions is much smaller than the next-neighbor repulsion V ef f between the trions Due to the scaling of the hopping parameter required to obtain a meaningful limit D → ∞, i.e. J → J/ √ z where z is the lattice connectivity, one finds J ef f → 0 in this limit, i.e. the trions become immobile while their nextneighbor interaction term survives. In this limit, the Hamiltonian is equivalent to an antiferromagnetic Ising model (spin up corresponds to a trion and spin down corresponds to a trionic-hole). At half-filling, clearly the most energetically favorable configuration is therefore to arrange the trions in a staggered configuration [52]. Moreover, due to quantum fluctuations, if we decrease the interaction starting from very large |U |, the spread of a single trion (which is proportional to J 2 /U ) increases and it is not a local object any more. In this case the trionic wave-function extends also to the nearest neighboring sites [34], as sketched in Fig. 2. This interpretation is in agreement with the observed behavior of the CDW order parameter C in Fig. 1. Indeed, at large |U |, C asymptotically rises to the value C = 1/2, corresponding to the fully local trions in a staggered CDW configuration. The presence of the CDW also explains the anomalously large value of residual entropy per site s res = k B ln 2 found when imposing a homogeneous trionic phase as in Ref. [22,23]. At strong-coupling in finite dimensions, even though the trions have a finite effective hopping J ef f , one would still expect that the augmented symmetry at half-filling favors CDW modulations with respect to a trionic Fermi liquid phase. In D = 1, 2 it is indeed known [17,27] that the CDW is actually stable with respect to the SF phase at half-filling for any value of the interaction, in contrast to the SU (2) case where they are degenerate [38]. Our results prove that in higher spatial dimensions this is not the case and there is a finite range of attraction at weak-coupling, where the c-SF phase is actually stable.
Further confirmation of the physical scenario depicted above is provided by the analysis of the single-particle spectral function ρ σ in the c-SF and t-CDW phases shown in Fig. 3. In the c-SF phase (Fig. 3(a)), the spectrum shows a gapless branch due to the presence of the third species which is not involved in the pairing, while the spectral function for species 1 (2 is identical) shows a gap. The situation is totally different in the t-CDW phase (Fig. 3(b)), where the spectral functions for the three species are identical but the lattice symmetry is broken into two sublattices. If we plot the spectral functions for the two sublattices (corresponding to two successive iterations in our DMFT loop) a CDW gap is visible. We would like to note that the sharply peaked structure of the spectrum is due to the finite number of orbitals in the ED algorithm. However, the size of the gap should not be affected significantly by the finite number of orbitals. Interestingly for |U | = 0.75W the size of the energy gap ∆ gap ≈ W is in very close agreement with the value obtained within SFA for the same value of the interaction [22,23], indicating that the gap most likely is only weakly affected by CDW ordering.   In order to characterize the system at finite temperature, we studied the evolution of the SF order parameter P as a function of temperature in the c-SF phase for different values of the coupling (Fig. 4(a)) and analogously for the CDW order parameter C in the t-CDW phase (Fig. 4(b)). The superfluid-to-normal phase transition at T SF c (U ) is also mirrored in the behavior of the spectral function for increasing temperature. The results shown in Fig. 5 indicate that the superfluid gap in the spectral function closes for T > T SF c (U ), signaling the transition to a normal homogeneus phase without CDW modulations.
At finite temperatures we also found a coexistence region of the trionic CDW wave phase and the color superfluid or normal homogeneous phases in a finite range of the interaction U (U c1 < |U | < U c2 at T = 0). We however leave a thorough investigation of the stability range of the t-CDW phase at finite temperature to future study, together with its dependence on the distance from the particle-hole symmetric point and on the dimensionality. Due to this coexistence region, we define the two critical temperatures T SF c (U ) and T CDW c (U ) plotted in the phase diagram in Fig. 6, where P (T ) |U and C(T ) |U vanish respectively above the c-SF phase and t-CDW phase. In agreement with the results obtained within SFA [22,23], we also found that the critical temperature T SF c (U ) has a maximum at T SF c /W ≈ 0.025 for |U |/W = 0.4. This is also in qualitative agreement with the SU (2) case [38], where the critical temperature has a maximum at intermediate couplings. Due to the presence of the CDW modulations in the trionic phase which are ignored in Ref. [22,23], we found also a second critical temperature T CDW c where charge density wave modulations in the trionic phase disappear.

Incommensurate density
In this section we consider the system for densities far from the particle-hole symmetric point. Specifically we investigate, using VMC and DMFT respectively, the implementation of the model (1) on a simple-square (cubic) lattice in 2D (3D) with tight-binding dispersion, i.e. k = −2J i=x,y(,z) cos(k i a), where a is the lattice spacing. In particular, we will find that away from the particle-hole symmetric point in the c-SF phase, the superfluidity always triggers a density imbalance, i.e. a magnetization. In order to address this feature quantitatively, we studied the system by adjusting the chemical potential µ in order to fix the total density n = σ n σ , allowing the system to adjust spontaneously the densities in each channel. Due to the spontaneous symmetry breaking of the SU (3) symmetry of the Hamiltonian in the color superfluid phase, it is indeed possible that, for a given chemical potential µ 1 = µ 2 = µ 3 = µ, the particle densities for different species may differ. If such a situation occurs, the systems shows a finite onsite magnetization m. As a more technical remark, we add that the choice of pairing channel, as explained in Sec. 3.1.1, is done without loss of generality: A specific choice will therefore determine in which channel a potential magnetization takes place, but not influence its overall occurrence. Here, since we fix the pairing to occur between species 1 and 2, we found a nonzero value of the magnetization parameter m = n 12 − n 3 , where n 12 = n 1 = n 2 . Therefore the paired channel turns out (spontaneously) to be fully balanced, while there is in general a finite density imbalance between particles in the paired channel with respect to the unpaired fermions.
The implications of the results presented in this subsection and in Sec. 5 for cold atom experiments, where the total number of particles of each species N σ = i n i,σ is fixed, will be discussed in Sec. 6. Combining the grand canonical DMFT results with energetic arguments based on canonical VMC calculations, we show that the system is generally unstable towards domain formation.
We first consider in Fig. 7 how the ground state properties of the 3D system evolve by fixing the coupling at |U |/W = 0.3125, where the system is always found to be in the c-SF phase for any density. We consider only densities ranging from n = 0 to half-filling n = 1.5. The results above half-filling can be easily obtained exploiting a particle-hole transformation. In particular one easily obtains where t and d are the average triple and double occupancies. The superfluid order parameter P increases (decreases) with the density for n < 1.5 (n > 1.5) and is maximal at half-filling. The average triple occupancy is instead a monotonic function of the density. Below half-filling, the magnetization m first grows with increasing density, then reaches a maximum and eventually decreases and vanishes at half-filling in agreement with the findings in the previous subsection. This means that in the c-SF phase for a fixed value of the chemical potential µ the system favors putting more particles into the paired channel than into the unpaired component. For n > 1.5 the effect is the opposite and m < 0. This behavior can be understood by considering that the equilibrium value of the magnetization results from a competition between the condensation energy gain in the paired channel on one side and the potential energy gain on the other side. Indeed the condensation energy found as a function of the density of pairs has a maximum at halffilling. For example in the weak-coupling BCS regime E cond is proportional to P 2 [49]. Therefore the condensation energy gain will increase by choosing the number of particles in the paired channel as close as possible to half-filling. On the other hand, for a fixed total density n, this would reduce or increase the unpaired fermions and consequently the potential energy gain, which is maximal for a non-magnetized system since U is negative. The competition between these opposite trends eventually determines the value of the magnetization in equilibrium, which is finite and rather small at this value of the coupling (see inset in Fig. 7). At half-filling no condensation energy gain can be achieved by creating a density imbalance between the superfluid pairs and the unpaired fermions since the condensation energy is already maximal. Therefore the spontaneous symmetry breaking in the color superfluid phase does not result necessarily in a density imbalance, which is however triggered by a condensation energy gain for every density deviation from the particle-hole symmetric point. We now consider the same system for fixed total density n = 1 and study the ground state properties as a function of the interaction strength |U | (see Fig.8). For weak interactions the system is in a c-SF phase. Upon increasing |U |, the order parameter P first increases and then shows the dome shape at intermediate couplings which we already observed for the half-filled case. Away from the half-filling, the value where P reaches its maximum is shifted to lower values of the interaction strength. The triple occupancy t, on the other hand monotonically increases with |U |. Interestingly the magnetization m(U ) has a non-monotonic behavior. At weak-coupling, magnetization m(U ) grows with increase of the interaction strength. For increasing coupling, m has a maximum and then decreases for larger |U |, indicating a non-trivial evolution due to competition between the condensation energy and the potential energies for increasing attraction. The spontaneous breaking of the SU (3)-symmetry is also well visible in the behavior of the double occupancies. Indeed in the c-SF for n < 1.5 we find d 12 > d 13 = d 23 . The difference d 12 − d 23 is however non-monotonic in the coupling and seems to vanish at |U |/W ≈ 0.35. Our interpretation is that beyond this point the SU (3)-symmetry is restored and the system undergoes a transition to a Fermi liquid trionic phase. Indeed for |U |/W > 0.35 we did not find any converged solution within our DMFT approach, neither for a homogeneous nor for a staggered criterion of convergence. This result is compatible with the presence of a macroscopically large number of degenerate trionic configurations away from the half-filling. A finite kinetic energy for the trions would remove this degeneracy, leading to a trionic Fermi liquid ground state. This contribution is however beyond the DMFT description of the trionic phase where trions are immobile objects. We can address the existence of a Fermi liquid trionic phase at strong-coupling using the VMC approach in 2D, which we will discuss in the following.
As already mentioned in Sec. 3, we use different trial wavefunctions to study the behavior of the system in the weak-(|U | ≤ W/2) and the strong-coupling (|U | > W/2) regimes. At weak-coupling the magnetization is expected to be very small and we can consider the results for the unpolarized system with n 1 = n 2 = n 3 to be a good approximation of the real system which is in general polarized. We found indeed that for |U | ≤ W/2 the system is in the c-SF phase with a finite order parameter P . As shown in Fig. 9, we obtain that P (U ) has a similar dome shape as in the 3D case. Unfortunately, we cannot directly address the trionic transition within this approach since it is expected to take place at intermediate coupling where both ansatz wave functions are inaccurate. We can however consider the system in the strong-coupling limit by using the effective trionic Hamiltonian of Eq. 28. In this way we can study the Fermi liquid trionic phase which we characterize by evaluating the quasiparticle weight, averaged over the Fermi surface Here Z k is extracted from the jump in the momentum distribution at the Fermi surface, which we approximate as where ∆k x (∆k y ) is the translational vector along the x (y) direction in the reciprocal lattice. In Fig. 10 we plot Z as a function of interaction strength |U |/W . By combining DMFT and VMC results we therefore have strong evidence of the system undergoing a phase transition from a magnetized color-superfluid to a trionic Fermi liquid phase at strong-coupling, when the density is far enough from the particlehole symmetric point.

Results: Constrained System (V = ∞)
As referred to in the introduction, actual laboratory implementations of the model under investigation using ultracold gases are often affected with three-body losses, which are not Pauli suppressed as in the SU (2) case. As discussed in Ref. [8], the three-body loss rate γ 3 shows a strong dependence on the applied magnetic field. Therefore the results presented in the previous section essentially apply to the case of cold gases only whenever three-body losses are negligible, i.e. γ 3 J, U . In the general case, in order to model the system in presence of three-body losses, one needs a non-equilibrium formulation where the number of particles is not conserved. However, as shown in Ref. [40], in the regime of strong losses γ 3 J, U , the probability of having triply occupied sites vanishes and the system can still be described using a Hamiltonian formulation with a dynamically-generated three-body constraint. To take it into account in our DMFT formalism, we introduce a three-body repulsion with V = ∞. Within VMC we directly project triply occupied sites out of the Hilbert space. We stress that finite values of V do not correspond to real systems with moderately large γ 3 since then real losses occur and a purely Hamiltonian description does not apply any more; only the limits γ 3 J, U and γ 3 J, U lend themselves to an effective Hamiltonian formulation.

Ground State Properties
In order to address how the system approaches the constrained regime with increasing V , we first used DMFT to study the ground-state properties of the model in 3D as a function of the three-body interaction V for a fixed value of the total density n = 0.48 and the two-body attraction |U |/W = 0.3125. We found that the average number of triply occupied sites t = n 1 n 2 n 3 (not shown) vanishes very fast with increasing V . The SF order parameter P and the densities in the paired and unpaired channels approach their asymptotic values already for V ≈ 3W or V ≈ 10|U |, as shown in Fig. 11 Therefore, we assume that we can safely consider the system to be in the constrained regime whenever V is chosen to be much larger than this value. Both the densities n σ and the superfluid order parameter P are strongly affected by the three-body interaction (see Fig. 11). For this value of the interaction, P and m are strongly suppressed by the three-body repulsion, even though both eventually saturate to a finite value for large enough V . However, as shown below, this suppression of the magnetization and SF properties is specific to the weak-coupling regime and for larger values of |U | both the SF order parameter P and the magnetization m are instead strongly enhanced in the presence of large V .
We now investigate the constrained case (setting V = 1000J ≈ 80W within the DMFT approach) where the total density is fixed as above to n = 0.48. Large values of the density imply an increase of the probability of real losses over a finite interval of time. Therefore we restrict ourselves to a relatively low density which is meant to be representative of a possible experimental setup.
We study the evolution of the ground state of the system in 2D and 3D as a function of the two-body interaction strength U . DMFT results in Fig. 12 show that in the threedimensional system the trionic phase at strong coupling is completely suppressed by the three-body constraint and the ground state is found to be always a color superfluid for any value of the attraction. This remaining c-SF phase shows however a very peculiar behavior of the magnetization m as a function of the attraction U . Indeed the magnetization m = n 12 − n 3 (n 12 = n 1 = n 2 ) steadily increases for increasing interaction and n 3 ≈ 0 (m ≈ n 12 ≈ n/2) already for U ≈ 12J = W .
Our explanation is that the three-body constraint strongly affects the energetic balance within the c-SF phase. Indeed, in the absence of V the magnetization was shown to be non-monotonic and to vanish in the SU (3)-symmetric trionic phase at strong-coupling. Now instead in the same limit the fully polarized c-SF system has a smaller ground state energy for fixed total density n. This result is fully confirmed by the VMC data for the 2D square lattice. As shown in the next section, combining these results essentially implies that a globally homogeneous phase with m = 0 is unstable in the thermodynamic limit with respect to domain formation whenever the global particle number in each species N σ = i n i,σ is conserved. By using the canonical ensemble approach of VMC, we can indeed address also metastable phases and study the effect on the energy of a finite magnetization for fixed total density n = 0.48. In particular we study the energy difference between the magnetized system and the unpolarized one with the same n, i.e. ∆E(m) = E(m) − E(0). Results shown in Fig. 13 indicate that at strong-coupling the energy decreases for increasing magnetization and the minimum in the ground state energy corresponds to the fully polarized system. In the inset of Fig. 13, we show ∆E as a function of the interaction strength for the fully polarized c-SF at strong-coupling, which decreases as ∆E ∼ 1/|U |. We also investigated the system in the weak-coupling regime, where our calculation shows that ∆E(m) has a minimum for very small values of the magnetization (not shown). This indicates that also in 2D the c-SF ground state at weak-coupling is partially magnetized, in complete agreement with the three-dimensional results. For weak coupling we approximate that the system is not magnetized (green lines and circles), while for strong coupling we assume that the system is fully polarized, i.e. contains only pairs (dashed blue line and squares). The dotted line corresponds to the superfluid order parameter in the atomic limit P ∞ (Constrained case) Within DMFT the order parameter P in the c-SF ground state shown in Fig. 12 is also increasing with |U | and saturates at strong coupling to a finite value, which we found to be in agreement with the asymptotic value in the atomic limit for the SU (2) symmetric case [38] The total number of double occupancies d is also an increasing function of |U | and saturates for very large |U | to the value n 12 = n/2 as in the strong coupling limit for the SU (2) symmetric system. This means that in the ground state the strong coupling limit of the SU (3) model is indistinguishable from the SU (2) case for the same total density n and two-body interaction U . As we will show in the next subsection, this is not any more true if we consider instead finite temperatures.
Similar considerations on the superfluid properties in the ground state apply to the two-dimensional case studied within the VMC technique. As the magnetization in the weak-coupling regime is very small, we approximated it to zero and consider an unpolarized system within the weak-coupling ansatz, while at strong-coupling we directly consider the system as fully polarized, i.e. containing only pairs. As visible in Fig. 14, P shows a similar behavior to the 3D case. Indeed at weak-coupling both, DMFT and VMC, show a BCS exponential behavior in the coupling, while at strongcoupling P converges to a constant.
Within VMC we also studied the condensation energy as explained in Sec. 3. Fig.  14b shows that the condensation energy first increases with the interaction strength U as expected in BCS theory, while it decreases as 1/U at strong-coupling as expected in the BEC limit for the SU (2) case [38]. Despite the fact that we cannot reliably address the intermediate region, there are also indications that the condensation energy has a maximum in this region.

Finite temperatures
We also investigated finite-temperatures properties for the three-dimensional case using DMFT. In Fig. 15, we show the evolution at finite temperature T of the SF order parameter P and of the magnetization m at fixed values of the interaction U . At low temperatures, the system is superfluid and the magnetization finite. With increase of the temperature, both P and m decrease and then vanish simultaneously at the critical temperature T = T c (U ). This clearly reflects the close connection between superfluid properties and magnetism in the SU (3)-symmetric case and is markedly different from the strongly asymmetric case which we studied in Ref. [42], where the density imbalance survives well above the critical temperature.
It is however remarkable that for |U | > U m ≈ W , m(T ) and P (T ) clearly show in Fig. 15 the existence of a plateau at finite T , indicating that the system stays in practice fully polarized in a finite range of temperatures. This allows us to define operatively a second temperature T p (U ) below which the system is fully polarized, while for T > T p instead the magnetization decreases and eventually vanishes at T c .
We summarize these results in the phase diagram in Fig. 16. Inside the region marked in orange (|U | > U m and T < T p ) the system is fully polarized and therefore identical to the SU (2) superfluid case. As we will see in the next section, in a canonical ensemble where the total number of particles N σ of each species is fixed, this analogy is not any more correct and we have to invoke the presence of domain formation to reconcile these findings with the global number conservation in each channel. Outside this region and below T c (solid blue line in Fig. 16), the c-SF is partially magnetized and therefore intrinsically different from the case with only two species. This is also visible in the behavior of the critical temperature where the SU (3)-symmetry is restored in the normal phase. We found indeed that the critical temperature first increases with the interaction strength |U |, similarly to the SU (2) case. Then for |U | = U m , the critical temperature T c suddenly changes trend and for larger |U | a power-law decrease T c ∝ 1/|U | occurs as shown in Fig. 16. In the SU (2) symmetric case this power-law behavior only appears for very large |U | (bosonic limit) [38], while in the SU (3) case this regime occurs immediately for |U | > U m . The smooth crossover in T c (U ) and the maximum of in the critical temperature characteristic of the SU (2) case, here are replaced by a cusp at |U | = U m , which marks the abrupt transition from one regime to the other.

Domain Formation
One of the main results of this work is the close connection between superfluidity and magnetization in the c-SF phase. Indeed we found that in the c-SF phase, away from the particle-hole symmetric point, the magnetization is always non-zero. On the other hand ultracold gas experiments are usually performed under conditions where the global number of particles N σ = i n i,σ in each hyperfine state is conserved, provided spin flip processes are suppressed. The aim of this section is to show that domain formation provides a way to reconcile our findings with these circumstances. In particular, combining DMFT and VMC findings, we will show that a globally homogeneous c-SF phase is unstable with respect to formation of domains with different c-SF phases in the thermodynamic limit.
To be more specific, we will consider the case when the global numbers of particles in each species are the same, i.e. N 1 = N 2 = N 3 = N/3, at T = 0, though the discussion can be easily generalized to other cases. The simplest solution compatible with N σ = N/3 is clearly a non-polarized c-SF phase with energy E hom per lattice site. This phase is actually unstable and therefore not accessible in a grand canonical approach like DMFT, where we fix the global chemical potential µ and calculate the particle densities n σ as an output. Since, as shown in Sec. 4 and Sec. 5, the system is spontaneously magnetized in the color superfluid phase out of half-filling, there is no way to reconcile the DMFT result with the global constraint N σ = N/3 assuming the presence of a single homogeneous phase. The VMC approach, on the other hand, operates in the canonical ensemble, and it can be used to estimate the ground state energy per lattice site for specific trial configurations. For the homogeneous configuration, we have E hom = E(m = 0) n , where n = N/M and M is the number of lattice sites.
Let us now contrast this situation with the spatially non-uniform scenario in which we have many color superfluid domains in equilibrium. Each of these domains corresponds to one of the solutions obtained above, and therefore this phenomenon can be seen as a special form of phase separation. For two or more phases to be in thermodynamic equilibrium with each other at T = 0, they need to have the same value of the grand potential per lattice site Ω = E − µn for the same given value of the chemical potential µ, while the onsite density of particles for each species n σ can be ; a specific example of a phase-separated configuration is plotted. Increasing the attraction strength reveals substantial differences between the two cases: (b) In the constrained case, domain formation persists to strong coupling, in parallel to the 3-component asymmetric situation [42]. The unpaired species are expelled from the paired regions, pairing up in other spatial domains. (c) In the unconstrained case instead, a spatially homogeneous trionic phase emerges [20,21]. different in the different phases. Possible candidate phases for the system considered in this paper are suggested by the underlying SU (3)-symmetry. Indeed if we consider c-SF solutions corresponding to different gauge fixing, i.e. with pairing in different channels, they will have the same total onsite density n and therefore the same energy and grand potential, since they correspond to different realizations of the spontaneously broken symmetry. If we consider for simplicity only the three solutions with pairing between the natural species sketched in Fig. 17, then this mixture of phases has globally the same number of particles N σ = N/3 in each hyperfine state whenever we choose the fraction of each phase in the mixture to be α = 1/3 and n = N/M in each domain. In fact in each domain we have the same densities n p in the paired channel and n u for the unpaired fermions, even though they involve different species in different domains. This scenario is therefore compatible with the global number constraint N σ = N/3 and we can compare its energy with the energy E hom of the globally homogeneous c-SF phase. The VMC calculations reported in Fig. 13 clearly indicate that for a fixed onsite density n, the ground state energy per lattice site is lower by having a finite magnetization, i.e. E(m) n < E(0) n and therefore E hom > E phase−separated = α 3 i=1 E i = E(m) and E i = E(m) is the energy per lattice site in the i-th domain. Thus a globally homogeneous c-SF phase has higher energy than a mixture of polarized domains with the same N σ and is therefore unstable with respect to phase separation.
It should be noted however, that the configuration sketched in Fig. 17 only represents the simplest possible scenario compatible with the global boundary conditions N σ = N/3. Indeed in the SU (3)-symmetric case we have continuous set of equivalent solutions, since solutions obtained continuously rotating the pairing state from 1-2 to a generic linear combination of species have the same energy and are therefore equally good candidates for the state with domain formation. Moreover, it is well known that having a continuous symmetry breaking is intrinsically different from the discrete case, because of the presence of Goldstone modes [17]. In large but finite systems, the surface energy at the interface between domains, which is negligible in the thermodynamic limit, will become relevant. On one hand a continuous symmetry breaking allows the system to reduce the surface energy cost through an arbitrarily small change of the order parameter from domain to domain, pointing toward a scenario where a large number of domains is preferable in real systems. On the other hand, when the system is finite, increasing the number of domains decreases their extension, reducing the bulk contribution which eventually defines number and size of the domains at equilibrium. Based on our current approaches, we cannot address the issue of what is the real domain configuration in a finite system, neither the question if different scenarios with microscopical modulations of the SF order parameter take place [53,54]. Similar conclusions concerning the emergence of domain formation in the c-SF phase have been already drawn in [20,21,37] and also in a very recent work [49], which addresses the same system in continuum space.
In real experiments both finite-size effects and inhomogeneities due to the trapping potential could play an important role in the actual realization of the presented scenario. Furthermore, as the SU (3)-symmetry in the cold atomic systems is not fundamental but arises as a consequence of fine-tuning of the interaction parameters, imperfections will also arise from slight asymmetries in these parameters. We have shown before [42] that in the strongly asymmetric limit, phase separation is a very robust phenomenon. We may therefore conjecture that interaction parameter asymmetries favor this scenario.
The combination of the findings in the present paper on the SU (3) case with those on the strongly-asymmetric case in [42] suggests that phase-separation in globally balanced mixtures is a quite general feature of three-species Fermi mixtures. However, the phases involved are in general different in different setups. In the stronglyasymmetric case in presence of a three-body constraint, the color superfluid phase undergo a spatial separation in superfluid dimers and unpaired fermions [42]. In this case, the presence of the constraint is crucial to the phase-separation phenomenon, as testified by its survival well above the critical temperature for the disappearance of the superfluid phase [42]. In the fully SU (3) symmetric case instead, the presence of the constraint only modifies the nature of the underlying color superfluid phase favoring fully polarized domains at strong coupling. The formation of many equivalent color superfluid domains can be seen as a special case of phase separation reflecting the SU (3) symmetry. In this case the phase separation phenomenon is strongly connected to the superfluid and magnetic properties of the color superfluid phase and it is expected to disappear at the critical temperature T c and for the peculiar particle-hole symmetric point at half-filling in the unconstrained case.

Conclusions
We have studied a SU(3) attractively interacting mixture of three-species fermions in a lattice with and without a three-body constraint using dynamical mean-field theory (D ≥ 3) and variational Monte Carlo techniques (D = 2). We have investigated both ground state properties of the system and the effect of finite temperature and find a rich phase diagram. For the unconstrained system, we found a phase transition from a color superfluid state to a trionic phase, which shows additional charge density modulation at half-filling. The superfluid order as well as CDW disappear with increasing temperature.
In the presence of the three-body constraint, the ground state is always superfluid, but for strong interactions |U | > U m the system becomes fully polarized for fixed total density n. It is remarkable that according to our calculations the system stays fully polarized in a range of low temperatures. For high temperatures a transition to the non-superfluid SU (3) Fermi liquid phase is found. The critical temperature has a cusp precisely at U m . This is in contrast to the SU (2)-symmetric case, where a smooth crossover in the critical temperature takes place.
The c-SF phase shows an interesting interplay between superfluid and magnetic properties. Except in the special case of half-filling, the c-SF phase always implies a spontaneous magnetization which leads to domain formation in balanced 3-component mixture.
The kinetic energy operator can be split in several contributions, where the subscripts indicate the change in the total number of double occupancies (N d,0 = N 12 d +N 23 d +N 13 d ), i.e.
(n i,σ h i,σ + h i,σ n i,σ )c † i,σ c j,σ (n j,σ h j,σ + h j,σ n j,σ ), (A.2) Here n iσ = c † i,σ c i,σ , h i,σ = 1 − n iσ and σ =σ =σ = σ. We note that whereas K 0 preserves the total double occupancy N d,0 , it contains two different types of terms: (i) terms that also preserve double occupancy in each channel N σσ d ( K a 0 part) and (ii) terms that change the double occupancy in two different channels such that the total double occupancy stays unchanged (K b 0 part). Thus, we can write We can also decompose the operators that change the total number of double occupancies into K 1 = K 12 1 + K 23 1 + K 13 1 , (A.6) where the superscripts give the type of double occupancies that are being created or destroyed. The canonical transformation can be written as an expansion to the second order Using the relation [V, K σσ ±1 ] = ±U σ,σ K σσ m and applying the projection P D , we arrive at ).
(A.11) 0. One can easily calculate that | iσ|H|t 0 | 2 = J 2 and E t 0 − E iσ = U σσ + U σσ , where σ = σ = σ = σ. So we obtain ∆E = zJ 2 U 12 + U 13 + zJ 2 U 12 + U 23 + zJ 2 U 13 + U 23 , (A. 15) where z is the number of the nearest neighbor lattice sites. The calculation above assumes that neighboring sites of a trion are not occupied. If one of the neighboring sites is occupied by another trion, then the energy gain per trion is given by The effective interaction between two trions on neighboring sites is therefore For the SU (3)-symmetric case this expression is simplified and we obtain Therefore the nearest neighbor interaction between trions is repulsive in the SU (3)symmetric case. The next step is to calculate the effective hopping of the trions. For this purpose one has to use third order perturbation theory Here |t 0 and |t 1 define local trions on lattice site 0 and the neighboring lattice site 1 respectively, |σ defines a state where a fermion with spin σ occupies the lattice site 1, and two other fermions are occupying the lattice site 0. Conversely |σσ defines a state where two fermions with spins σ and σ occupy the lattice site 1. On the lattice site 0 we have only a fermion with spin σ = σ, σ . For any σ and σ the matrix elements are given by t 0 |H|σ = σ|H|σσ = σσ |H|t 1 = −J, E t 0 − E σ = U σσ + U σσ and E t 1 − E σσ = U σσ + U σ σ , where σ, σ and σ are three different hyperfine-spins. So we obtain . (A.20) where σ, σ and σ are different from each other in the sum. In the SU (3)-symmetric case, the expression again simplifies to So we obtain the following effective Hamiltonian [48] H ef f = −J ef f i,j Here t † i is the creation operator of a local trion at lattice site i and n T i = t † i t i is the trionic number operator.