Flavour-selective localization in interacting lattice fermions via SU(N) symmetry breaking

A large repulsion between particles in a quantum system can lead to their localization, as it happens for the electrons in Mott insulating materials. This paradigm has recently branched out into a new quantum state, the orbital-selective Mott insulator, where electrons in some orbitals are predicted to localize, while others remain itinerant. We provide a direct experimental realization of this phenomenon, that we extend to a more general flavour-selective localization. By using an atom-based quantum simulator, we engineer SU(3) Fermi-Hubbard models breaking their symmetry via a tunable coupling between flavours, observing an enhancement of localization and the emergence of flavour-dependent correlations. Our realization of flavour-selective Mott physics opens the path to the quantum simulation of multicomponent materials, from superconductors to topological insulators.

A large repulsion between particles in a quantum system can lead to their localization, as it happens for the electrons in Mott insulating materials. This paradigm has recently branched out into a new quantum state, the orbital-selective Mott insulator, where electrons in some orbitals are predicted to localize, while others remain itinerant. We provide a direct experimental realization of this phenomenon, that we extend to a more general flavour-selective localization. By using an atom-based quantum simulator, we engineer SU(3) Fermi-Hubbard models breaking their symmetry via a tunable coupling between flavours, observing an enhancement of localization and the emergence of flavour-dependent correlations. Our realization of flavour-selective Mott physics opens the path to the quantum simulation of multicomponent materials, from superconductors to topological insulators.
Interactions shape the collective behavior of manyparticle quantum systems, leading to rich phase diagrams where conventional and novel phases can be induced by a controlled variation of external stimuli. The most direct example is perhaps the Mott insulator, where the large repulsion amongst particles leads to an insulating state despite the non-interacting system being a metal or a superfluid 1 . In solid-state physics the interest in Mott insulators is reinforced by the observation that the proximity to a Mott transition is a horn of plenty where a variety of novel and spectacular phases can be observed, high-T c superconductivity 2,3 being only the tip of an iceberg 4 .
In recent years it has become clear that the standard SU(2) Fermi-Hubbard model is only a specific example, as many interesting materials require a description in terms of "multicomponent" Hubbard models, e.g. when the conduction electrons have an additional orbital degree of freedom. These systems are not merely more complicated, rather they host new phenomena, challenging the standard paradigm of Mott localization 5 . Indeed, when the symmetry between orbitals is broken, by some field or internal coupling, electrons in specific orbitals (or some combinations of them) can be Mott-localized while others remain itinerant, leading to surprising "orbital-selective Mott insulators" 6 .
Orbital-selective Mott physics has become a central concept for the description of a new class of high-T c superconductors based on iron 7 , as it can describe the anomalies of the metallic state 8,10 and the orbital character of superconductivity 11 in those systems. However, a clean observation of selective Mott transitions is intrinsically hard in solid-state systems, because of the limited The global symmetry is broken by a coherent Rabi driving Ω between two internal states. b) In the absence of coupling, the system experiences a phase transition from an SU(3) metal to a Mott insulator as the hopping is reduced. c) The Rabi driving lifts the degeneracy between the states and, in a dressedbasis picture, causes them to acquire different energies. The competition with the hopping can drive a transition from a metal to an insulator already in the noninteracting case. experimental control over the microscopic parameters and the orbital degree of freedom. The paradigm of selective Mott physics is itself far from being fully explored and has the potential to become a powerful framework to understand a variety of phenomena, from superconductivity to topological properties, in multicomponent quantum materials.
In this work we take a broader perspective and treat orbital-selective Mott physics as an example of the general concept of "flavour-selective Mott localization" 12 , which can be realized in a variety of multi-flavour systems, where the flavour can be the spin, the orbital or any other quantum number. We realize a minimal instance of this phenomenon by means of an atomic quantum simulator based on the optical manipulation of nuclear-spin mixtures of ultracold two-electron 173 Yb atoms. This platform allows the realization of multicomponent systems with global SU(N) interaction symmetry 13,14 , as in recent works reporting the realization of SU(N) quantum wires 15 , SU(N) Mott insulators 16,17 and, more recently, SU(N) quantum magnetism 18 . Here we introduce a novel approach to break the symmetry of SU(3) Fermi-Hubbard systems in a controlled way, which allows us to go beyond the investigations in the solid state and to observe directly the two key signatures of selective-Mott physics 12 : an overall enhancement of Mott localization and the onset of flavour-selective correlations.
The experiment is performed with three-component ultracold 173 Yb Fermi gases with total atom number up to N = 4 × 10 4 and temperature T 0.2T F . The atoms are trapped in a cubic 3D optical lattice (lattice constant d = λ/2 = 380 nm), which realizes the multi-flavour Hubbard HamiltonianĤ (1) where α, β ∈ {1, 2, 3} indicate the fermionic flavours (corresponding to nuclear spin states m = +5/2, +1/2 and −5/2, respectively), t is the tunnelling energy between nearest-neighboring sites i, j , U is the onsite repulsion energy between two atoms of different flavours, andV T = κ i,α R 2 in αi describes the effects of a slowlyvarying harmonic trapping potential (where R i is the distance of site i from the trap center). In the absence of the 4th term of Eq. (1),Ĥ has an intrinsic global SU(3) symmetry, which is ensured by the invariance of atom-atom interactions on the spin state and by the realization of spin-independent optical potentials, i.e. U , t and κ do not depend on α. This symmetry is explicitly broken by the 4th term, which describes a coherent onsite coupling between flavours |1 and |2 . This coupling is provided by a two-photon Raman process with Rabi frequency Ω/h 20 (where h is the Planck constant). At the single-particle level, this coupling lifts the degeneracy between the flavours, creating two dressed combinations |± = (|1 ± |2 ) / √ 2, energy-shifted from |3 by ±Ω/2, as sketched in Fig. 1c.
We use an adiabatic preparation sequence to produce an equilibrium state of the atomic mixture in the optical lattice, with equal state populations N 1 = N 2 = N 3 = N/3 and in the presence of the coherent coupling Ω between states |1 and |2 20 . In order to characterize the degree of localization of the particles in the lattice, we measure the number of atoms N d in doubly-occupied sites (called "doublons" in the following) with photoassociation spectroscopy, as in previous experiments that demonstrated the onset of Mott localization for ultracold fermions 16,19,20 .
In Fig. 2a we report typical measurements of N d as a function of the total atom number N . In a harmonically trapped system, the rate of change of N d vs N gives information on the core compressibility 19,21 . A vanishing value of N d over an extended range of N signals the presence of an incompressible state with one atom per site in the center of the trap (since adding particles does not lead to a proportional increase of doublons), while the critical N above which N d then departs from zero can be connected to the magnitude of the energy gap protecting the localized phase. (1). Fig. 3a shows a zero-temperature phase diagram obtained from Dynamical Mean-Field Theory (DMFT) 22 for the homogeneous system (V T = 0) with a uniform 1/3 filling (one atom per site) and equal number of particles for each flavour. The phase diagram clearly has the same shape of the experimental one, showing that the Hubbard U and the coupling Ω cooperate in driving the system from a metallic phase to a more localized state. Between the standard metal and the Mott phase, we find a region where the degree of correlations (as measured by the quasiparticle weight) is strongly selective, the coupled flavours being much more localized than the uncoupled one. It is evident that both selective and global Mott localization occur at a smaller U if Ω is included.
In order to connect this effect to the experimentally measured signal, in Fig. 3b we show the doublon fraction f d obtained from the homogeneous DMFT results after a local-density approximation (LDA) analysis, to take into account the effect of the harmonic trapping in V T .
The reduction of f d with increasing Ω is in agreement with the experimental observations reported in Fig. 2. The lack of a quantitative matching with Fig. 2d can be attributed to imperfections in the initial state preparation and to the finite temperature of the experiment, resulting in an average entropy per particle S/N ≈ 2.5k B , which is known 21 to produce an effect on the double occupancies in a trapped system already at Ω = 0.
In Fig. 3c we finally show the result of a zerotemperature DMRG calculation of f d for a harmonically trapped 1D system 20 . Although a quantitative agreement with the experimental data should not be seeked (because of the different dimensionality and the finite temperature of the experimental realization), the overall behavior, i.e. a reduction of doublons for increasing U and Ω, is correctly captured and agrees also with DMFT. This also indicates that the phenomena we are exploring are generic and qualitatively independent on the dimensionality.
In addition to the enhancement of Mott localization, the DMFT phase diagram of Fig. 3a shows that Ω is expected to result in a flavour-dependent localization of the many-body system. The flavour-selective behavior can be experimentally detected by resolving the spin character of the doublons, i.e. counting how many atoms form doublons in each of the three pairs |12 , |23 , |31 by means of state-selective photoassociation at a high magnetic field 20 . In Fig. 4 we plot the quantity γ(12) = (12) is the number of atoms forming doublons in the |12 channel, as a function of Ω and fixed U = 2.6D. The measured value at Ω = 0 agrees with the expectation for an SU(3)-symmetric system, for which γ(12) = 1/3. As Ω is increased and the SU(3) symmetry is broken, γ(12) diminishes, eventually approaching zero for Ω ≈ D. The doublons acquire a strongly flavour-selective behavior.
This suppression of |12 doublons is triggered by the polarization effect in the rotated |± basis, which can be understood, at a qualitative level, already from the simplified, noninteracting case sketched in Fig. 1c: while |23 and |31 doublons can be formed by two fermions in the lowest single-particle states |+ and |3 , |12 doublons can be formed only if the two fermions occupy states |+ and |− , therefore with an additional energy cost Ω/2.
Interactions then increase this effect 12 , leading to strong flavour-selective results even for small values of Ω.
The crosses in Fig. 4 are the result of a DMFT calculation in which the flavour populations are kept equal and the harmonic trapping V T has been taken into account in a LDA approach. The results of this numerical calculation are in good agreement with the experimental findings, with a larger degree of selectivity for the theoretical calculation. We argue that a better agreement could be seeked by including finite-temperature effects in the calculation: indeed, we expect the state-selective behavior to be reduced by the thermal occupation of higher-energy states, leading to an effective reduction of the polarization in the |+ , |− basis.
We note that the experimental observation of a finite number of |12 doublons at small, but finite Ω is an indication of the validity of the protocol used for the preparation of the atomic state, which is different for Ω = 0 and Ω > 0 20 . We also note that, despite the polarization in the dressed |+ , |− basis discussed before, we have verified that the populations of the bare states |1 , |2 , |3 remain always equal under all the experimental conditions we have considered. This is indeed an important aspect of our experiment. If the populations of the bare states were not fixed (and in particular N 3 was left free to adapt), the dressed states would be populated according to the scheme in Fig. 1c, leading to N 1 = N 2 > N 3 , which favours flavour-selective physics already in the non-interacting system. The experimental state preparation procedure counteracts the trivial differentiation between flavours by forcing an even occupation. Therefore, we conclude that the flavour selectivity we observed is essentially due to quantum correlations induced by interactions.
In the theoretical DMFT calculations, the equal population constraint is enforced including an external field h 20 that favours the occupation of |3 in such a way to match the experimental condition N 1 + N 2 + N 3 = N/3. Comparing with Ref. 12, where the populations were left free, we observe that the quantum correlations leading to the selective regime survive to the inclusion of the field h. In general terms, single-particle effects trigger flavour selectivity, but the inclusion of interactions strongly enhances the differentiation, turning a minor modulation of kinetic energy into a quantitative phenomenon which can also lead to a selective Mott transition 24 . This is a very general framework, underlying many investigations of multicomponent models.
Exploiting an idealized quantum simulator we have obtained a clear-cut evidence for correlation-induced flavourselective physics, where the SU(N) symmetry-breaking coupling Ω is only the trigger of a flavour-selective phenomenon which is fundamentally driven by correlation effects. Our first realization of multicomponent Hubbard physics with coherent internal couplings opens new paths for the quantum simulation of new classes of materials ranging from high-temperature superconductors to interacting topological insulators as described by the Bernevig-Hughes-Zhang model 23 . As the coupling is realized with a nonzero momentum transfer (i.e. non collinear Raman beams), a full range of possibilities will emerge, including the study of magnetic crystals 25  We measure double occupancies (doublons) in the 3D lattice exploiting a one-color photoassociation (PA) process which transfers the atoms in doubly-occupied lattice sites into highly excited molecular states 19 . The photoassociated atoms are rapidly lost from the lattice due to the fast decay of the molecular states, leaving a sample consisting only of singly-occupied sites. The number of double occupancies is thus inferred by first measuring the total number of atoms in the sample and then subtracting the number of atoms in singly occupied lattice sites remaining after the PA pulse. For the regime of atom number considered in this work we expect the number of triply-occupied sites to be negligible.
Our measurement procedure starts with a rapid freeze of the atomic density distribution that we realize increasing the lattice depth up to V 0 = 25E rec (E rec = h 2 /8md 2 , where h is the Planck constant, m the atomic mass and d the lattice spacing) in a time 1.5 ms. At the final lattice depth tunnelling between lattice sites is completely negligible. Two different schemes are then employed for the PA excitation depending on the specific measurement that we want to perform.
In order to measure the total number of doublons independently on their spin composition (results of Fig.  2 of the main text), a 5 ms long, σ − polarized PA pulse is shone on the atomic sample with an intensity I P A = 90 mW/mm 2 . The PA pulse excites a molecular line that is 796.2 MHz red-detuned with respect to the 1 S 0 (F = 5/2) → 3 P 1 (F = 7/2) atomic transition 16 . During the process, atoms are subjected to a low magnetic field B = 3 G which defines a quantization axis for the spin, but is low enough not to unveil the Zeeman substructure of the molecular line, as it is shown in the spectrum reported in Fig. S1a.
In the presence of a higher external magnetic field, the PA transition frequency is shifted accordingly to the projection M T of the molecular total angular momentum on the quantization axis S1 . M T is a quantum number conserved during the PA process, with a value given by the relation M T = m 1 + m 2 + σ, where m 1,2 are the nuclear spin projections of the the two atoms composing the molecule and σ = −1 is the angular momentum transferred by the PA photon. As long as the sum m 1 + m 2 unambiguously characterizes an atomic pair, the molecular Zeeman substructure can thus be employed as a tool to distinguish between doublons with different spin flavours. Taking this into consideration, in order to detect only doublons with a particular spin composition {m 1 , m 2 }, we apply a magnetic field B = 66 G, which induces a Zeeman shift of the order of several MHz between molecular lines with different M T projection. In this case the spectrum is complicated by the presence of a plethora of features originated by the Zeeman splitting of different transitions, which makes it difficult to label individual lines. The association of the PA lines to a particular spin flavour has been realized preparing different SU(2) mixtures and acquiring individual spectra for each of them, as it is shown in Fig. S1b. In particular, at this magnetic field, we identify three strong PA lines at 778. For the measurements reported in Fig. 4 of the main text, γ(12) is calculated from a spin-resolved measurement of doublons in the |12 and |31 channels as γ(12) = N d (12)/(N d (12) + 2N d (31)), taking into account N d (23) = N d (31). This assumption is justified by the equal population of the two flavours |1 and |2 after the adiabatic preparation sequence described in Sec. S.V, as shown in Fig. S5. In Fig. S2 we compare the points obtained with the analysis presented in the main text and those given by the unsupervised K-Means analysis, taking the same number of clusters. The good agreement between them ensures that the specific choice of data analysis procedure does not introduced biases on the results. We note that the K-Means method implies only a constraint on the number of clusters, allowing for a different width among them, differently from the main analysis procedure, where the bin size was fixed.
As the datasets for the flavour-selective double occupancies are rather large, we can increase the number of groups K to pinpoint the specific range of atom numbers above which doublons form. The flavour-selective averages obtained with the K-Means analysis are fitted with the piecewise function: where N 0 defines the threshold atom number below which atoms are localised. In Fig. S3 we show the K-Means averages for the coupled flavours |12 and the uncoupled ones |23 , together with the results of the fit. The difference between the fitted thresholds for the two doublon flavours ∆N 0 = N 0 (12) − N 0 (23) is shown in the inset of Fig. S3 for different values of Ω/D. This analysis demonstrates that, by increasing the Raman coupling intensity (i.e. lifting the flavour degeneracy more), the threshold for |12 doublons moves to a higher number of atoms, showing an increased localisation in the centre of the trap for the coupled flavours. This result, derived from the unsupervised machine-learning analysis, is consistent with the measurements reported in Fig. 4, and provides another strong indication in support of flavour-selective localization for the symmetry-broken SU(3) Hubbard Hamiltonian.

S.IV. IMPLEMENTATION OF RAMAN COUPLING
The coherent coupling between states |1 (m = 5/2) and |2 (m = 1/2) is realized by exploiting a two-photon σ + /σ − Raman transition. The Raman coupling is implemented with two co-propagating λ = 556 nm laser beams characterized respectively by angular frequencies ω and ω + δω. In order to reduce the inelastic photon scattering rate, the two beams are 1.754 GHz blue-detuned with respect to the 1 S 0 → 3 P 1 (F = 7/2) intercombination transition. A 150 G magnetic field is used to define a quantization axis and to remove the degeneracy between the six states of the 173 Yb ground-state manifold, which are split according to their nuclear spin m by 207 × m Hz/G.
The σ + /σ − coupling between m = +5/2 and m = +1/2 is obtained by setting the polarization of the two beams to be orthogonal with respect to the quantization axis and adjusting the frequency difference δω/2π in order to compensate the Zeeman splitting and the residual Raman light shift between the two states. As explained in detail in the supplementary material of Ref. S3, this choice for the polarization makes the Raman light shifts largely spin-dependent, thus making the m = +1/2 ↔ m = −3/2 coupling nonresonant and effectively excluding state m = −3/2 from the dynamics. For a similar reason, also the m = −5/2 ↔ m = −1/2 coupling is nonresonant and state |3 (m = −5/2) does not participate to the Raman dynamics.
On the basis of the discussion above, the Raman coupling on the basis formed by the three states {|1 , |2 , |3 } can be described by the following 3x3 rotating-waveapproximated Hamiltonian where Ω is the angular Rabi frequency associated to the Raman coupling and δ is the two-photon detuning with respect to the resonance frequency for the |1 ↔ |2 transition. No detuning is indicated for state |3 as it is decoupled from the Raman process. In the resonant case δ = 0,Ĥ R corresponds to the third term in the many-body Hubbard Hamiltonian in Eq. (1) of the main text. The Rabi frequency Ω can be experimentally changed by adjusting the intensities I 1 and I 2 of the two Raman beams, according to the relation Ω ∝ √ I 1 I 2 . In order to assess the value of Ω for given values of the Raman beam intensities I 1 and I 2 , we induce resonant Rabi oscillations between states |1 and |2 and determine Ω from a sinusoidal fit of the experimental data, as shown in Fig.  S4.

S.V. ADIABATIC STATE PREPARATION
A degenerate Fermi gas of 173 Yb atoms, initially confined in a crossed optical dipole trap with harmonic trapping frequencies ω x,y,z = 2π × {55, 96, 73} Hz, is transferred in a simple-cubic 3D optical lattice using a two-step 3s-long adiabatic loading procedure S4 . The optical lattice is described by a potential energy V (x, y, z) = V 0 sin 2 (πx/d) + sin 2 (πy/d) + sin 2 (πz/d) , where d = λ/2 is the lattice spacing and V 0 is the lattice energy depth along each of the three principal axes.
In the first 2 seconds, the lattice intensity is ramped up with a first spline ramp from V 0 = 0 to V 0 = 4E rec (E rec = h 2 /8md 2 is the recoil energy). After this first step, the lattice intensity is further increased with a 1slong spline ramp to the final value V 0 ranging from 4E rec to 15E rec . During the lattice loading, the depth of the crossed dipole trap is progressively reduced in such a way to keep the harmonic trapping frequencies constant along the three principal axes, independently from the value of V 0 .
For the measurements at Ω = 0 the loading procedure described above is applied to a balanced 3-component mixture of 173 Yb atoms in the three spin states |1 , |2 and |3 . The mixture is prepared before the lattice rampup procedure with a sequence of optical pumping pulses on the 1 S 0 → 3 P 1 transition, following the techniques discussed in Ref. 15. The populations of the three spin states are all equal, N 1 = N 2 = N 3 = N/3, with experimental imperfections on the order of a few % at most.
For the measurements at Ω = 0 the loading procedure starts with a 2-component unbalanced mixture of atoms in states |1 and |3 . With a proper choice of optical pumping pulses we adjust the initial populations to be N 1 = 2N/3 and N 3 = N/3. After the lattice loading, we switch on the Raman beams, initially far detuned from any two-photon transition, and perform an adiatic frequency sweep to bring them resonant with the |1 ↔ |2 transition. The resonant condition is reached by means of an exponential frequency sweep of the form that reduces the two-photon detuning from δ 0 Ω to δ = 0 in a time T and time constant τ (see an example in Fig. S6a). This procedure is an adiabatic passage that brings an atom in state |1 to the lowest-energy dressed state of the Raman-coupled system |+ = (|1 + |2 ) / √ 2 (see also the sketch in Fig. 1 of the main text). We note that, due to the initial spin unbalance and because of the equal-weighted admixture of states |1 and |2 in the Raman-dressed states, at the end of the detuning ramp the population is equally distributed among the spins, N 1 = N 2 = N 3 = N/3, as it was natively in the loading protocol employed for Ω = 0.
The ramp parameters δ 0 , T and τ are chosen according to a numerical analysis in which we solve the timedependent Schrödinger equation associated to the Raman Hamiltonian in Eq. (S.2), verifying that at the end of the ramp the initial state is effectively projected onto the Hamiltonian lowest energy eigenstate. Experimentally, we check the adiabaticity of this procedure by verifying that the initial unbalanced |1 -|3 mixture can be recovered with a reversed detuning ramp following that in Eq. (S.3). To further assess the loading fidelity, we verify the time stability of the spin populations at the end of the loading ramp, as shown in Fig. S5. In both the cases, we measure population differences of a few % at most, validating the effectiveness of this loading procedure.

S.VI. NUMERICAL VALIDATION OF THE LOADING PROCEDURE
In order to assess the validity of the adiabatic state preparation at the many-body level, we have developed a numerical simulation based on the exact diagonalization of Eq. (1) on a reduced-scale version of our system. We consider N three-component fermions in a 1D lattice with N s sites, retaining all the relevant terms of Eq. (1): hopping, repulsive interactions, Raman coupling between states |1 and |2 . We work in the occupation number basis in which the Hilbert space is constituted by the where we have added the last term to take into account the detuning of the Raman coupling (already introduced in the rotated-wave-approximated Hamiltonian of Eq. (S.2)), that is crucial in the state preparation protocol.
We simulate the effect of our loading procedure for a system composed by N = 3 particles in a lattice with N P = 3 sites. In this case, the dimension of the Hilbert space associated to the system is In order to validate the loading procedure we calculate the fidelity between the time-evolved state |Ψ(t) = e −iĤt/ |Ψ 0 and the target wavefunction, that is defined as the lowest-energy eigenstate |Ψ f of the final Hamiltonian with δ = 0 that is compatible with the constraint on the populations N 1 + N 2 = 2 and N 3 =1. An example is shown in Fig. S6b for the actual detuning ramp used in the experiment to prepare the state at U = 2.6D and Ω = 0.35D. Fig. S6c shows the full spectrum of the system as a function of the Raman detuning. There it is evident that the ground state of the system at large detuning |Ψ 0 is adiabatically connected with the ground state |Ψ f on resonance. We note that, even if the system initially starts with an unbalanced mixture of two species only, after this procedure it behaves as a mixtures of three species N 1 = N 2 = N 3 = 1 all interacting among themselves: in particular, in the presence of tunnelling and interactions between the atoms, the state rotation induced by the detuning sweep does not happen uniformly in all the sites of the lattice, making the atoms initially in state |1 distinguishable and therefore, interacting, at the end of the state preparation sequence.

S.VII. DYNAMICAL MEAN-FIELD THEORY CALCULATIONS
We solve the model in Eq. (1) using Dynamical Mean-Field theory (DMFT) 22 , a non-perturbative theoretical method which maps a lattice model onto an effective impurity model. The interaction of the site with the rest of the lattice is approximated by an effective dynamical bath which is determined self-consistently. This allows to fully capture the quantum dynamics while freezing spatial fluctuations beyond mean-field. We solve the model in a Bethe lattice of bandwidth W , which is known to account correctly for the physics of three-dimensional lattices. The impurity model is solved using an exact-diagonalization solver developed at SISSA S5 which requires to describe the bath in terms of a finite number of levels N b that we fix to 9 in the calculations reported in this work.
DMFT is naturally formulated in a grandcanonical ensemble, where the total density is fixed by a chemical potential µ and the occupations of the various flavours are not fixed. In order to enforce the experimental population constraint N 1 + N 2 = 2N 3 in the presence of the coupling Ω, we have to include an external field which counterbalances the tendency to favour (for total occupation of 1) the occupation of the coupled flavours |1 and |2 . Therefore the Hamiltonian is supplemented by the termŝ H f ields = i −µ αn αi + h (n 1i +n 2i − 2n 3i ) , (S. 6) where µ and h have to be determined self-consistently to reach a total occupation of 1 fermion per site. For a uniform infinite lattice, we have drawn a phase diagram based on the flavour-resolved quasiparticle weight Z α , which is 1 for a non-interacting fermion and it is reduced by the interactions. The global Mott transition is reached when all the Z α vanish, while a sharp crossover is observed between a standard metal and a region of selective correlation where the quasiparticle weight of the coupled species rapidly drops to a value smaller than 0.05.
Besides the calculations for a uniform lattice, we also have used DMFT to estimate the effect of the harmonic trapping potentialV where R i is the distance of the site i from the center of the trap. In order to reduce the computational cost, we have employed a local-density approximation (LDA), which assumes that the local properties of the system can be computed from a uniform model with the local values of the physical parameters, in the present case of the potential given by the harmonic trap. This amounts to have a local value of the chemical potential µ = µ − κ i,α R 2 i .
In order to build the global observables for the trapped system, we have solved the model given by Eqs. (1) and (S.7) for a wide range of values of µ spanning the whole range of densities from 0 to 3. For every value of µ we have found the value of h which gives N 1 = N 2 = N 3 .
The global observables (total number of particles, doublons, flavour-selective doublons) are then obtained integrating the local counterparts over the whole system. In this way we obtain the plots of doublons as a function of the total number of particles shown in the manuscript.

S.VIII. DENSITY MATRIX RENORMALIZATION GROUP CALCULATIONS
We also solve the model in Eq. (1) by means of a Density Matrix Renormalization Group (DMRG) calculation powered by the ITensor library S6 . We fix the size of the system as L = 30 and the total number of particles as N = 21. We simulate the internal states as effective physical sites with the redefinition of proper tunnelling and interaction couplings, which results in a total of N L = 90 lattice sites. We find the ground state of the Hamiltonian in Eq. (1) by performing 35 − 50 DMRG sweeps, with a maximum bond dimension in the interval (20,200). To enforce particle number conservation in the Raman-coupled case, we introduce an additional chemical potential defined as the h term in Eq. (S.7). The number of doubly occupied sites shown in Fig. 3 also includes the number of triple occupancies (giving a maximum contribution of ∼ 10% at small values of U , Ω).