Direct probing of the Mott crossover in the SU($N$) Fermi-Hubbard model

The Fermi-Hubbard model (FHM) is a cornerstone of modern condensed matter theory. Developed for interacting electrons in solids, which typically exhibit SU($2$) symmetry, it describes a wide range of phenomena, such as metal to insulator transitions and magnetic order. Its generalized SU($N$)-symmetric form, originally applied to multi-orbital materials such as transition-metal oxides, has recently attracted much interest owing to the availability of ultracold SU($N$)-symmetric atomic gases. Here we report on a detailed experimental investigation of the SU($N$)-symmetric FHM using local probing of an atomic gas of ytterbium in an optical lattice to determine the equation of state through different interaction regimes. We prepare a low-temperature SU($N$)-symmetric Mott insulator and characterize the Mott crossover, representing important steps towards probing predicted novel SU($N$)-magnetic phases.

The Fermi-Hubbard model (FHM) is a cornerstone of modern condensed matter theory. Developed for interacting electrons in solids, which typically exhibit SU(2) symmetry, it describes a wide range of phenomena, such as metal to insulator transitions and magnetic order [1]. Its generalized SU(N )-symmetric form, originally applied to multi-orbital materials such as transition-metal oxides [2], has recently attracted much interest owing to the availability of ultracold SU(N )symmetric atomic gases. Here we report on a detailed experimental investigation of the SU(N )symmetric FHM using local probing of an atomic gas of ytterbium in an optical lattice to determine the equation of state through different interaction regimes. We prepare a low-temperature SU(N )symmetric Mott insulator and characterize the Mott crossover, representing important steps towards probing predicted novel SU(N )-magnetic phases [3].
Although the Fermi-Hubbard Hamiltonian has been the object of a large number of studies in the past decades, reaching a complete understanding has remained an elusive task, even for the spin-1/2 case. For repulsive interactions, the SU(2) FHM is known to give rise to a paramagnetic Mott insulator, while antiferromagnetic order emerges below the Néel temperature. Moreover, the FHM is believed to capture the essential physics of d -wave superfluidity in high-temperature superconductors [1,4]. The development of experimental implementations of the three-dimensional (3D) FHM with ultracold atoms has provided a new approach for advancing our understanding of strongly correlated fermions in lattices [5].
The recent realization of degenerate gases of strontium and ytterbium in optical lattices allows to go beyond the conventional spin-1/2 FHM and to access the largely unexplored physics of the more general SU(N ) symmetry [3,6]. In contrast to the well-studied SU(2) FHM, much less is known about the extended SU(N > 2)-symmetric case. Calculations are mostly limited to T = 0, low dimensions or approximated correlations, but point to very rich phase diagrams including novel magnetic and spin liquid phases [6][7][8][9][10][11][12][13][14]. Moreover, accurate predictions of thermodynamic quantities are even harder to obtain than for the SU(2) case, as most numerical algorithms fail for larger N due the unfavourable scaling of the Hilbert space.
Fermionic 173 Yb has nuclear spin I = 5/2, but no electronic angular momentum in the ground state (J = 0). Therefore, nuclear and electronic angular momenta are decoupled, making the atomic interactions independent of the nuclear spin state and the system SU(N )symmetric, with N ≤ 2I + 1 = 6 being the number of populated nuclear spin states [15]. One important manifestation of large N is an enhanced Pomeranchuk effect, leading to a suppression of particle-hole excitations in the lattice [16]. Owing to this, first evidence of an incompressible phase in an SU(6) 173 Yb gas has recently been reported, at an entropy level that does not support a Mott insulator in lower spin gases [17].
The SU(N )-symmetric FHM for arbitrary spin multiplicity N can be written as [3] Here, i, j denote neighbouring lattice site indices, t is the tunnelling matrix element between them, U is the on-site interaction and V i is a position-dependent energy offset which accounts for the confining potential. We denote with t * = 12 t the kinetic energy associated with the bandwidth of the 3D lattice. The operatorĉ iσ annihilates a fermion at site i with spin index σ = 1, ..., N andn iσ =ĉ † iσĉiσ are the respective number operators. In this work, we measure and analyse the equation of state (EoS) of an SU(N ) Fermi gas in a cubic optical lattice, for temperatures above the magnetic ordering temperature. The EoS is a thermodynamic relation that contains all the macroscopic properties of the system and, due to its generality, is particularly well suited to benchmark numerical simulations [18,19]. We take advantage of the confining potential and the local density approximation (LDA) to map the trapped heterogeneous gas with fixed number of particles to a locally homogeneous gas in the grand canonical ensemble with µ i = µ 0 − V i , where µ 0 is the chemical potential at the centre of the trap. The validity of the LDA for fermions trapped in optical lattices above the Néel temperature was previously verified by numerical calculations [20]. We determine the EoS for the density n(µ, T, N, U, t * ) over a wide 2 range of parameters. In particular, we focus on the highest spin multiplicity of our system N = 6 and on the case N = 3, which was the subject of several theoretical studies [10][11][12][13][14]21]. By deriving the local compressibility directly from the measured EoS we are able to detect the emergence of the incompressible Mott phase.
In our experiment, we start by preparing a degenerate Fermi gas of 173 Yb with initially N = 6 equally populated spin components via evaporative cooling in a crossed optical dipole trap (see Methods for details). We then set N by removing individual spin components. The Fermi gas is then loaded into the lowest energy band of a 3D optical lattice with cubic symmetry and lattice spacing d = λ/2 operating at a wavelength of λ = 759 nm. We vary the lattice depth between 3 E r and 15 E r , with E r = h 2 /2mλ 2 being the recoil energy, a range for which the tight-binding approximation is valid. Adjusting the lattice depth allows us to tune the system from U/t * = 0.128 +0.004 −0.008 to U/t * = 11.0 +1.1 −1.0 , spanning a range of two orders of magnitude.
In order to measure the local atomic density in the trap, we probe the cloud by performing in situ spinindependent absorption imaging along the z-axis of the lattice. The spatial resolution is approximately 1.2 µm ≈ 3.2 d. Due to the high optical density of the trapped cloud, saturated absorption imaging at high light intensity is used (see Methods and ref. [22] for details). The resulting integrated two-dimensional density distributioñ n(x, y) is shown in Fig. 1a. Exploiting the carefully characterized geometry of the trap configuration, we determine the local 3D density n(r, y) by performing an inverse Abel transform with y as the symmetry axis. [23]. The transformation is stable for pixels with radial distance r ≥ 7 µm from the symmetry axis. In Fig. 1b, the azimuthally averaged density is shown as a function of the distance to the trap centre for different values of the lattice depth. For the averaging, data within ±4.8 µm along the y-axis was used.
Two regimes can be distinguished: For low lattice depths, where the interaction energy is smaller than the kinetic energy (U < t * ), the density decreases smoothly from the centre to the edge of the trap, indicating that the system is compressible everywhere. For high lattice depths (U > t * ), we observe the formation of a plateau of constant density, as expected for a Mott insulating state [24]. In this region, the atomic density computed using an independent calibration is compatible with the expected value of one fermion per lattice site (see Methods). Due to the higher precision of this measurement, the density of an SU (6) insulator in this region is used as the density reference in this paper.
To determine thermodynamic quantities, the EoS is obtained by relating the 3D atom density distribution and the chemical potential: is done for different values of interaction strength U/t * and number of spin components N = 3, 6. In Fig. 2, examples for small (U/t * = 0.128), intermediate (U/t * = 0.89) and large (U/t * = 3.6) interaction strengths are given.
In the presence of a weak lattice (U t * ), the system is in a normal metallic state, well described by the SU(N ) Fermi-liquid theory [25] with enhanced mass and interactions due to the presence of the lattice. In particular, for low densities nd 3 1, interactions can be neglected and the EoS is well approximated by that of a non-interacting Fermi gas, as shown in Fig. 2a. We use this observation to determine µ 0 of the gases with U/t * ≤ 0.89 (see Methods for details). For high densities nd 3 1 the approximation fails and an SU(N ) theory that considers interactions in a weak lattice would be required. For deep lattice potentials, tunnelling is highly suppressed (t k B T ) and interactions are strong (U t * ). In this regime, the lattice sites can approximately be regarded as independent. We construct a low-tunnelling model using a high-temperature series expansion (HTSE) up to O(t/k B T ) 2 [16,26] (see Methods for details). This model fits the measured EoS data well, with T = 0.13(1) U/k B and T = 0.18(1) U/k B for SU (6) and SU (3) respectively (see Fig. 2c). The lower temperature in the SU(6) case is due to both the lower initial temperature and the stronger Pomeranchuk effect. As the HTSE fit sensitivity to the trap confinement is higher than the precision of the independent trap calibration, we allow a variation of the confinement parameter in the model within the bounds given by the calibration.
The role of the SU(N ) symmetry in the EoS is simple in the limit of negligible interactions: the density scales linearly with N for fixed chemical potential. This dependence is well reproduced by the data as shown in Fig.  2a. For high lattice depths, and therefore larger U/t * , the effect is two-fold: Firstly, the Pomeranchuk effect leads to lower temperatures in systems with higher N [16,17]. Secondly, for a given finite temperature, the EoS in the strongly interacting regime is N -dependent in a non-trivial way due to the different quantum statistics [13]. To illustrate this, we compare the low-tunnelling model for the SU(2) case and the EoS for SU(6) in Fig.  2c. As shown in the inset, a measurable discrepancy is present due to the different N -dependent quantum statistics, which cannot be attributed to different temperature T or central chemical potential µ 0 .
In the regime of intermediate interaction strength, comparable to the kinetic energy (U ≈ t * ), the system is in a strongly correlated many-body state. No exact determination of the EoS is available yet to compare to our data, displayed in Fig. 2b. Nevertheless, having modelfree access to the EoS allows to directly determine the local compressibilityκ = n 2 κ = ∂n/∂µ| T . This can be used to probe the emergence of the incompressible phase and to study the crossover between metal and Mott insulator, as shown in Fig. 3. For the strongly interacting case we distinguish a metallic outer layer 0 < nd 3 < 1 and a metallic core 1 < nd 3 , separated by a Mott shell with nd 3 = 1. In this regime, the low-tunnelling model compares well with our results, with very small compressibilities in the Mott regions.
To characterize the onset of the Mott phase, we determine the minimum of compressibilityκ min around unit filling, in a range nd 3 ∈ [0.85, 1.15], as a function of interaction parameter (Fig. 3). We observe a suppression of κ min by roughly one order of magnitude when increasing the interaction strength. For large U/t * the compressibility saturates at a minimum value and the Mott shell is formed. The behaviour is consistent with numerical calculations [21].  Fig. 2c. c, minimal compressibilityκmin/κ0 as a function of interaction strength, whereκ0 is the compressibility of a non-interacting SU(6) Fermi gas at T = 0 and n = 1/d 3 . The minimum of the compressibilityκ(n) in an interval around unit filling (shaded area in insets) is used asκmin. Error bars denote the confidence intervals of the linear regression used to obtaiñ κ.
To further characterize the state in the lattice, we estimate the entropy per particle s of the SU(6) gas. We determine the entropy before and after a round trip, consisting of loading the atoms into the lattice and back into the dipole trap, by measuring the degeneracy parameter T /T F in the dipole trap, where T F is the Fermi temperature. This gives a lower and an upper estimate for s in the lattice (see Fig. 4). For low lattice depths, we observe a constant entropy rise associated with the transfer from the bulk into the lattice and back. Above U ≈ t * , with the appearance of the incompressible phase, we find an additional entropy increase. This increase also persists for very low ramp speeds, which could indicate diminished adiabaticity due to the insulating phase crossover hindering mass flow. However, we observe that also for these high lattice depths, the entropy is at or below the maximum spin entropy even after a full round trip that includes the reversed second ramp. This indicates that the orientations of the spins in the lattice are likely not fully random, allowing for the presence of partial spin correlations in the system.
In the regime of validity of the low-tunnelling model, lower entropy corresponds to lower temperature and Mott shells with sharper edges. Although fitting the model provides a good match to the measured SU(6)-EoS data as shown in Fig. 2c, the obtained temperature is significantly higher than that expected for the very low entropy measured in the round-trip experiment. In contrast, in the SU(3) case the measured round-trip entropy is higher than the saturated spin entropy and the temperature obtained from the in situ HTSE fit is in the expected range.
An interpretation for the SU(6) case could be that the suppression of mass flow close to the Mott crossover prevents the formation of sharp Mott shell edges, yielding higher fitted temperatures. Nevertheless, the measured low entropy after the round-trip hints that this non-adiabatic process would have to be mostly reversible, except for the observed increase seen in Fig. 4. Understanding a possible edge-softening contribution of entropy stored in that metallic region would require a quantitative solution of the SU(6) 3D theory at low filling and including spin correlations, a challenge for current numerical methods.
In conclusion, we use in situ probing to measure the equation of state of an SU(N )-symmetric Fermi gas in a 3D lattice. We obtain a very low compressibility in the Mott insulating phase and measure an entropy below that of uncorrelated spins. These findings make the system a promising starting point towards novel magnetically ordered many-body states in highly spin-symmetric systems [3,7,8,13,14,25]. An incompressible phase The entropy of the gas in the lattice is constrained to the region between these two measurements. The entropy of a system with our parameters and fully random spin orientations is close to s/kB = ln 6 (red line). Error bars and blue shaded region are the s.e.m. of the temperature measurements.
with filling nd 3 = 1 is expected to faithfully realize the SU(N )-symmetric Heisenberg model, paving the way towards studying unexplored systems with reduced dimensions such as chains [27][28][29] and new magnetic phases in square lattices [7].

Preparation of SU(N ) Fermi gases
Approximately 3.5 million ground-state atoms of 173 Yb, with equal populations in the nuclear spin states, are loaded into a crossed optical dipole trap operating at a wavelength of λ = 1064 nm. Forced evaporative cooling produces a degenerate Fermi gas of 5000 atoms per spin state at temperature T = 0.07(1) T F , where T F is the Fermi temperature. The gas is very weakly interacting since k F a s 0.07, with a s = 199.4 a 0 being the s-wave scattering length and a 0 the Bohr radius. We perform state preparation by driving the 1 S 0 → 3 P 1 optical transition to remove unwanted nuclear spin states from the trap, in the presence of a homogeneous magnetic field that lifts the spin state degeneracy. Using this technique we generate an SU(3) Fermi gas with a temperature of 0.15(1) T F , with a residual fraction of unwanted spin components below 5% per component.

Loading into the optical lattice
At the end of evaporation, the atoms are transferred into a cubic optical lattice in two steps. The sample is first loaded into a shallow lattice with depth V = 3 E r in 150 ms, avoiding band excitations. We then ramp the lattice depth to the final value, between V = 3 E r and V = 15 E r , in 150 ms. The atoms are trapped in the combined harmonic confinement produced by the lattice beams and the crossed optical dipole trap. The confinement frequencies vary between ω x,y,z = 2π × (31, 42, 183) Hz at 3 E r , ω x,y,z = 2π × (21, 33, 183) Hz at 7 E r and are approximately constant for V > 7 E r . We verified the validity of the harmonic approximation for our experimental trap configuration, by taking into account the combined Gaussian beam profiles. The variation of U/t * in the region occupied by the atoms is estimated to be below 8%.
Dilute samples are loaded into the lattice in order to minimize losses and heating from three-body recombination, which is not suppressed by Pauli blocking for N > 2 (see Supplementary Information for more details).

High-intensity imaging
The atom cloud has a typical optical density around 2 in the trap centre. In order to perform in situ absorption imaging we saturate the imaging transition 1 S 0 → 1 P 1 with light intensity I = 15 I sat , where I sat = 60 mW/cm 2 . Details on the calibration procedure are in the Supplementary Information. The optical resolution is about 1.2 µm, determined as described in ref. [30]. The imaging pulse has a duration of 5 µs, sufficiently short to avoid atoms to escape from the focal plane and to minimize the Doppler shift due to photon scattering. Imaging is performed in the absence of magnetic fields in order to have spin-independent detection. Nevertheless, a 6% difference between the densities of the Mott plateaus is found comparing the SU(6) and SU(3) gases, likely caused by the different line strengths within the hyperfine substructure of the imaging transition.

Analytical models
In the absence of interactions, the EoS of a Ncomponent harmonically confined Fermi gas is where Li 3/2 denotes the poly-logarithm of order 3/2. This expression is also valid in the presence of a weak lattice potential, provided that the effective mass associated with the dispersion of the lowest band is used.
In the strongly interacting regime, the low-tunnelling model is obtained using a high-temperature series expansion (HTSE) of the grand canonical potential up to second order in t/k B T [16,17,26]: with z being the number of next neighbours in the lattice, x = e βµ , y = e −βU and Ω 0 (µ, T ) the grand potential and Z 0 the partition function in the atomic limit.

Fitting and calibration
The effective optical scattering cross section σ of the atoms is calibrated using the known EoS of the polarized Fermi gas in the dipole trap with known trap frequencies.
We obtain σ/σ 0 = 0.222 ± 0.034 with σ 0 = 3 2π λ 2 . The HTSE model fit to the measured EoS, effectively fitting the amplitude of the plateau, is more sensitive to the optical cross section and yields σ/σ 0 = 0.233 ± 0.003, which is used in this work. The plateau density for the N = 6 data shown is therefore fixed to 1/d 3 , whereas an independent fit for the N = 3 case yields nd 3 = 0.94 ± 0.02, possibly due to residual spin dependence of the imaging technique.

ADDITIONAL INFORMATION
Correspondence and requests for materials should be addressed to S.F.

HIGH-INTENSITY IMAGING CALIBRATION
In order to correctly determine the density of the gas using high-intensity absorption imaging, we use the modified Lambert-Beer law [S1] which accounts for saturation of the imaging transition: where OD is the optical density, I i the incident light intensity, I f the final light intensity after absorption and I sat the saturation intensity. The saturation intensity I sat 60mW/cm 2 is calculated from the linewidth and wavelength of the transition, and additionally verified by an intensity-dependent linewidth measurement. The α-parameter is extracted by varying the saturation parameter I/I sat and adjusting its value in a way that the measured optical density from equation (S.1) is independent of the light intensity [S1]. It accounts for a transition strength lower than 1 and, potentially, a systematic error when determining the value of the light intensity impinging on the atoms.

CHARACTERIZATION OF THREE-BODY LOSSES
In order to avoid three-body losses in the lattice we keep the central average filling below nd 3 1.7. For higher filling factors we observe a fast density decay in the central region of the cloud in addition to a slow density decay due to vacuum losses and technical heating of the lattice beams. We fit a double exponential decay to the core density of the cloud n(t) = n d * e −t/τv + n 3 * e −t(3/τv+1/τ3) (S.2) where n d is the density of singly and doubly occupied sites, n 3 is the density of triply occupied sites and τ v ,τ 3 denote the vacuum lifetime and three-body loss timescale respectively. We extract a loss rate γ = 1/τ 3 = 2.4(3) Hz in a 15E r deep lattice. The scaling of the loss rate as a function of the lattice depth is found to be compatible with the expected scaling of a three-body decay rate