Equation of State and Thermometry of the 2D SU($N$) Fermi-Hubbard Model

We characterize the equation of state (EoS) of the SU($N>2$) Fermi-Hubbard Model (FHM) in a two-dimensional single-layer square optical lattice. We probe the density and the site occupation probabilities as functions of interaction strength and temperature for $N = 3, 4$ and 6. Our measurements are used as a benchmark for state-of-the-art numerical methods including determinantal quantum Monte Carlo (DQMC) and numerical linked cluster expansion (NLCE). By probing the density fluctuations, we compare temperatures determined in a model-independent way by fitting measurements to numerically calculated EoS results, making this a particularly interesting new step in the exploration and characterization of the SU($N$) FHM.

The interest in the square lattice SU(2) Fermi-Hubbard Model (FHM) has been historically driven by its suitability to describing cuprate superconductors, owing to their layered character and exceptionally simple band structure near the Fermi surface.For other more complex and multi-orbital materials, however, descriptions with N > 2 spin components have long been used, which, in addition to being of fundamental interest, provide an elegant approximation of degenerate orbitals using a higher symmetry group.Larger N systems, in particular in 2D geometries, are relevant for describing the physics of transition-metal oxides [1][2][3], orbitally-selective Mott transitions [4][5][6][7], graphene's SU(4) spin valley symmetry [8], twisted-bilayer graphene [9][10][11][12], the Kondo effect [13,14], heavy fermion behavior [15], and achieving robust itinerant ferromagnetism [16,17].The SU(N ) FHM is a special case of the N > 2 models that enjoys a higher symmetry group that stabilizes quantum fluctuations [18], making it a fertile ground for theory, and constituting a baseline to more complex multi-orbital models.The determination of the N > 2 equation of state (EoS) of the SU(N ) FHM is an important milestone in the attempt of understanding its properties.However, the exponential scaling of the Hilbert space with N and the increased severity of the fermion sign problem [19] make its numerical simulation more challenging than the N = 2 case [20][21][22][23].
In this Letter, we probe the equation of state of the two-dimensional SU(N ) FHM in a square lattice at intermediate temperatures in both the metallic and the Mott regime and compare our results with numerical calculations.In particular, we determine the in-lattice temperature and entropy by fitting experimental data using numerical methods such as determinantal quantum Monte Carlo (DQMC) [61,62] and numerical linked cluster expansion (NLCE) [63,64].We additionally determine the entropies in the 2D bulk before loading and after unloading from the lattice potential, and separately characterize the system inside the lattice with a thermometry relying on the fluctuation-dissipation theorem (FDT) based on the measurement of density fluctuations, without requiring modeling by theory.
The SU(N ) FHM Hamiltonian is given by: where ĉ † iσ and ĉiσ represent the fermionic creation and annihilation operators at site i with spin σ ∈ {1 . . .N }, niσ = ĉ † iσ ĉiσ is the number operator, ⟨i, j⟩ denotes nextneighbor lattice sites, t is the hopping amplitude, U is the on-site interaction strength and µ denotes the chemical potential, which absorbs the contribution of the trap confinement in the local density approximation (LDA) [66].
In this work, we directly probe the local density, components of the site-occupation distribution, and the density fluctuations within the detection resolution of a few lattice sites.By differentiating the density with respect to the local chemical potential, we evaluate the isothermal compressibility κ = ∂n/∂µ| T .Crucially, we implement a 2D single-layer SU(N ) ensemble that we probe with perpendicular high-resolution absorption imaging.This avoids integrating over inhomogeneous stacks of 2D systems [67,68] and allows us to directly access the density profile without complex reconstruction techniques required in 3D [37] and access density fluctuations as an additional thermodynamic in-situ observable.
In our experiment, we start by loading a degenerate Fermi gas of 173 Yb with tunable N ≤ 6 equally populated components [see Fig. 1(a)] and an entropy per particle s/k B > ∼ 1.0 into the single, horizontal layer of a vertical lattice.In this layer, we adiabatically ramp up a 2D square lattice potential with a wavelength of λ = 759 nm and a spacing of d = λ/2 [see Fig. 1(b)].By modifying the lattice depth, we can tune the strength of the interactions.We measure the density distribution using in-situ, saturated absorption imaging with a spatial resolution of approximately 2 µm ≈ 5d [65].
The measured 2D density n(x, y) of an SU(6) ensemble is shown in Fig. 1(c) for different interaction strengths and the same initial state preparation in the 2D bulk (the potential without in-plane lattices).Because of the harmonic confinement generated by the Gaussian profile of the lattice beams, the chemical potential varies across the trap, sampling different regions of the EoS.For increasing interactions, and in particular when the on-site interaction is larger than the square lattice bandwidth (U/t > ∼ 8), we observe the emergence of plateaus at integer density which we associate with an incompressible regime, a signature of a Mott insulating state.
As a distinctive probe of number squeezing effects in and close to the Mott regime, we determine the occupation number distribution by measuring the parityprojected density.After tuning U/t, we freeze the motion of the atoms by rapidly increasing the lattice depth and applying a photoassociation beam [36], which converts on-site pairs into excited-state molecules that are subsequently lost.The process removes > 99% of the on-site pairs and ≈ 5% of the remaining atoms [65].Fig. 1(d) shows the distribution of the singly-occupied sites corresponding to the same states of Fig. 1(c).The increase in depletion in the center with increasing interaction strength is a consequence of number squeezing to a high atom pair fraction.
To access different spin degeneracies, we prepare N < 6 ensembles by removing spin components using optical pumping [65].In Fig. 2(a-d), we show the EoS as a function of the local chemical potential µ/U for N ∈ {6, 4, 3}.The chemical potential at a given location is calculated from the potential of the trap µ(x, y) = µ 0 − 1 2 κ x x 2 + κ y y 2 .The exact shape of the potential is determined by fitting the density, where the trap frequencies are left as free parameters.We use a combined fit of the densities for N = 3, 4 and 6 for each separate U/t, but verify that separate fits for each N return values compatible with those of the combined fit.The fit of the EoS is performed in two dimensions, leaving as free parameters the temperatures T (U/t, N ) and the chemical potential µ 0 (U/t, N ) at the center of the trap.The theoretical density is convolved with the reconstructed point spread function (PSF) [65] to take into account the imaging imperfections.
For the EoS of Fig. 2, each spin mixture has been prepared with the same initial entropy per particle s/k B = FIG.2. Equation of state for the SU(N ) Fermi-Hubbard Model with N = 6 (blue), N = 4 (purple) and N = 3 (red).(a)-(d) Density (circles) and singly-occupied sites (diamonds) as a function of the chemical potential.Data for N = 4 and N = 6 have been offset by 0.5 nd 2 and 1.0 nd 2 along the vertical axis, respectively.Continuous lines associated to the density curves correspond to the fit of the EoS calculations to the total density as described in the text.The theory used for the fit is DQMC for U/t = 2.3 (1) and NLCE for the other values of U/t.The results from this fit model are also used to calculate the expected pair / single site distribution measurement.The chemical potential is defined with respect to the reference half filling [nd 2 (µ = 0) = N/2].For each U/t, we fit the average of 15 frames with similar atom number after center of mass alignment [65].Error bars are the standard error of the mean (s.e.m.).(e) Temperature according to the fit of the EoS shown in Fig. 2(()a-d).(f)-(h) Entropy per particle.Horizontal line: entropy in the 2D bulk before loading into the lattice; triangles, squares, hexagons: entropy in the lattice according to the fit of the EoS; small circles: entropy in the 2D bulk after a round-trip experiment.The entropy in the bulk takes into account the effect of the interactions and the 3D anisotropic density of states [65].Error bars correspond to the s.e.m. of the fit results.
1.2(1) in the 2D bulk before ramping up the lattice.We fit and benchmark NLCE and DQMC [65] which are commonly used state-of-the-art methods for finitetemperature SU(2) Hubbard models in the regime we are considering but have only recently been extended and applied to the SU(N ) experimental regime, which requires calculations away from nd 2 = 1 [23,41].This is, to our knowledge, the first application of SU(N ) NLCE to noninteger filling, and to the calculation of the occupancy distributions.Moreover, compared to previous works, the calculation has been extended to higher orders [65] to ensure a better convergence at low temperatures.For U/t = 7.5(4) and 10.4(6) we fit both DQMC and NLCE, observing an excellent agreement between the theory [65] and the experiment.For U/t = 33(2), we use NLCE and a high-temperature series expansion (HTSE-2), observing also in this case an excellent agreement [65].For U/t = 2.3(1), the temperature lies below the range of convergence of NLCE and we resort to DQMC alone.In Fig. 2(a-d), for the cases in which we fit more than one model, we only plot the NLCE results, because the lines would overlap.
In addition to the total density, in Fig. 2(a-d) we also characterize the distribution of on-site occupation numbers by removing doublons using the pair removal process described above.Experimental measurements (diamonds) are compared with the NLCE prediction (lines) based on the fit of the density, without additional free fit parameters, and agree well with the experimental data whenever available.As opposed to the N = 2 case, where only double occupancies are allowed, higher occupancies occur for N > 2. Although the numbers of these occupancies are small for the results considered at the temperatures and chemical potentials presented here, the photoassociation technique can be used to probe triple occupancies and their dynamics [69].
The harmonic confinement of the trap returned by the density fit can be compared to the confinement obtained from an independent measurement of the oscillatory motion of the atoms in the combined dipole potentials [65].We find a discrepancy between about 13 % for U/t = 7.5(4) and 40 % for U/t = 33(2), which is not fully explained by tolerances or the trap loading model.A possible contribution could be a lack of adiabaticity during the loading into the lattice [70][71][72].However, neither varying the speed of the lattice ramps up to a fac-tor of four (up to 1 s length) nor variations of the atom number lead to significant changes in the fit results.This would require the non-adiabatic effects to produce minimal changes in density and parity profiles [65].
In Fig. 2(e) we plot the temperatures obtained by the fits of the EoS.We observe a smaller temperature for larger N , a behavior expected due to the Pomeranchuk effect [30,73], but somewhat weaker than the ideal theoretical predictions [21,73] with the temperatures for N = 4 and 6 differing from each other by up to 20 %.We interpret this as a consequence of the heating not depending on N during the loading process, resulting in different entropies in the lattice.This is supported by the results presented in Fig. 2(f-h).We find that, despite the initial entropy in the 2D bulk before loading into the lattice being independent of N , the entropy returned by the fit of the EoS is larger for larger N , which explains the weakening of the Pomeranchuk effect.We also determine the entropy per particle in the 2D bulk after a round-trip experiment, which adds an inverted ramp back to the 2D bulk system.In this case, we obtain entropies comparable to those reported by the fit in the lattice for N = 3 and N = 4 but smaller for N = 6, similar to previous observations [37] and potentially indicating nonadiabatic effects in the preparation or return ramp.
Complementary to the measurement of the EoS, the new possibility to directly access the density in the 2D SU(N )-system allows us to probe the density fluctuations.For an integration area of size A ≫ d 2 , the variance of the detected atom number is related to the isothermal compressibility κ and the temperature T through the fluctuation-dissipation theorem (FDT) [74]: By measuring the density fluctuations, we can access the temperature with spatial resolution, largely independently from the EoS models and without relying on fits of the potential parameters and to the theory [60].In Fig. 3(a) we show such density fluctuations as a function of the chemical potential for different U/t values and N = 6.For strong interactions, we observe a reduction of the fluctuations in the proximity of nd 2 = 1, where we expect an incompressible Mott-insulating regime.Notably, the fluctuation amplitude is determined by area integration as described in Eq. ( 2), and therefore agrees with the thermodynamic fluctuations from the FDT as opposed to the expected on-site fluctuations δn 2 0 = ⟨n 2 ⟩ − ⟨n⟩ 2 (grey dashed line).This discrepancy illustrates the role of nonvanishing short-range density correlations.
The FDT holds locally for each density.In thermal equilibrium, the ratio between the fluctuations and the compressibility is constant.We use the FDT to check this assumption and extract the temperature of the system.For this purpose, we determine the isothermal compressibility κ(µ) directly from the density profile data with three-point differentiation and fit the temperature T FDT as the proportionality factor between the fluctuations and the compressibility.In Fig. 3(b) we compare T FDT (diamonds) with the temperature T EoS (hexagons) returned by the fit of the EoS.We observe a good agreement for all interactions.Moreover, we see that the residuals of the FDT analysis typically show similar temperatures at the center and at the edge of the cloud, indicating that there are no strong deviations from thermal equilibrium [65].
In conclusion, we report the measurement of the equation of state of the 2D SU(N ) FHM across the Mott crossover for temperatures comparable with or below the hopping energy and we compare the experiment with state-of-the-art numerical models.Moreover, with direct access to a single 2D plane system, we can independently determine temperatures in the experiment with spatial resolution using density fluctuation analysis, which al-lows one to e.g.cross-check thermal equilibrium.This measurement characterizes the EoS also in regimes hard to reach by current numerical methods.When compared to the experimental data, we find the theoretical calculations describe well the properties of the SU(N ) gas for the applicable range of temperatures.The temperature measurements also indicate that thermal equilibration is not inhibited even in the case of deep lattices in a temperature range where the onset of spin correlations between sites is expected.
The implementation of the directly-accessible 2D ensemble, together with the accompanying theoretical description, paves the way towards more direct quantum simulation of the typically-2D models of interest in naturally occurring systems with SU(N > 2) representations such as transition metal oxides and orbitally-selective Mott transitions.An intriguing example is the case of cerium volume collapse, where there is a long-standing debate whether the single orbital Hubbard model (N = 2) or the double-orbital Hubbard model (N = 4) [75][76][77][78] is the correct description.While in the condensed matter examples the SU(N ) symmetry is typically only approximately realized, cold atom representations provide an essentially exact realization of SU(N ), allowing to implement fully SU(N )-symmetric and previously purely theoretical models.It should even be possible to smoothly connect both regimes in a continuous way by controlled symmetry breaking using e.g.optical state manipulation or state-dependent potentials [32,39,79], but more generally alkaline-earth-like quantum simulations of SU(N ) FHM can provide insight into the validity of the SU(N ) approximation in more realistic models.Our experiment starts by loading a spin-balanced unpolarized mixture (N = 6) of approximately 1.6 × 10 6 173 Yb atoms from a magneto-optical trap into a crossed optical dipole trap (XODT), where we perform forced evaporation to T /T 3D F < 0.2.We then perform a second stage of evaporation with an optical gradient, leading to an ensemble of N p ≈ 2 × 10 3 atoms in the central plane of a vertical lattice with wavelength λ = 759.3nm and lattice spacing d vert = 3.9(3) µm.The vertical band gap is 3.95(1) kHz, determined from a measurement of the 0 → 2 band resonance with parametric modulation.In this configuration, the coefficients of the harmonic confinement are (κ x , κ y )d2 = [1.35(3),2.23( 6)] × h.Numerical simulations [1] predict that for N p < ∼ 5 × 10 3 , nearly all the atoms should be loaded in the single central plane of the lattice.We verify this assumption with a momentum refocussing technique similar to the one described in Ref. [2].
Mixtures with N = 4 are prepared in the XODT before the evaporation by optically pumping the population of the two nuclear spin states m F ∈ {±1/2} to the other spin states with four circularly-polarized pulses on the1 S 0 → 3 P 1 intercombination line at 50 G.Mixtures with N = 3 are prepared in a similar fashion by pumping the spin components m F ∈ {±1/2, −3/2} to the other states with five pulses.In both cases, we check the balancing of the spin components with an Optical Stern-Gerlach (OSG) technique in time of flight [3].The standard deviation of the population of the selected spin components is below 5 % per component for SU (6) and SU(3) and below 8 % for SU (4).The residual fraction of unwanted spin components is below 5 % per component.

B. Lattice loading and Hubbard parameters
The two orthogonal in-plane lattices have the same depth and they are tuned to the desired value following the ramp profile shown in Fig. S1.The ramp speed has been chosen to minimize the entropy in the roundtrip experiment.After reaching the desired lattice depth, we quickly switch off the vertical lattice to avoid lightassisted collisions and to perform in situ imaging.
We calibrate the lattice depths along the two directions by measuring the band gap between the lowest and the second excited band with parametric modulation.U and t are obtained by numerically solving the band structure for the measured band gap and U is calculated as the Wannier overlap in the lowest band [4] The lattice depth and the hopping along the spatial x and y directions are the same.
ported in Tab.S1.We cross-check the calculation by directly measuring U with modulation spectroscopy for V = 13 E rec (U/t ≈ 43) and 19 E rec (U/t ≈ 178).For this measurement, we first tune the lattice depth to the desired value.We then modulate the depth of both lattices with an amplitude of 2 % to 6 % and at the same time apply the pair removal pulse.For modulation frequencies resonant with U , atom losses are enhanced.For 13 E rec , we measure U = h • 886(8) Hz for an expected value of h • 856(8) Hz.For 19 E rec , we measure U = h • 1065(8) Hz for an expected value of h • 1069(9) Hz (Fig. S2).

C. Pair removal
The pair-removal operation is performed with a photoassociation beam working on the 1 S 0 → 3 P 1 intercombination line, resonant with a molecular transition detuned by −599.28 (8) MHz with respect to the singleparticle transition with a bias magnetic field of 1 G.The beam is coplanar with the lattices and has a power of 30 mW.For pair removal, we first ramp the lattice to the desired U/t value according to the ramps described in Sec.S.I B, quench the depth to 30 E rec in 0.5 ms and then apply the photoassociation pulse for 10 ms.S1.Hubbard parameters for our system.V is the lattice depth along one direction (both lattices have the same depth).For U/t = 2.3(1), the next-nearest-neighbor hopping is h • 12 Hz.
The on-site atom number after pair removal is: where p α is the probability of having α particles on a given site (αp α is the average number of particles on a site with occupation α and n = N α=0 αp α ) and we neglect the tunneling during the quench.We correct Eq.S.1 to take into account the experimental imperfections of the photoassociation beam: where ⌊•⌋ represents the floor function, γ s and γ s + γ d represent the decay rate of the singlons and the doublons and we neglect the fast decay of the states with more than two particles per site [5].For γ s , γ d → 0 we recover ñPR → n PR .

Efficiency calibration
We calibrate the efficiency of the pair removal by looking at the atom losses as a function of the pulse duration in the deep Mott insulating regime [U/t = 33 (2)].Neglecting the sites occupied by more than two particles, we fit a double-exponential model: where N p = N s + N d is the total atom number, N s the number of singlons, N d the number of doublons.We obtain 1/γ d = 1.2(2) ms and 1/γ s = 200 (11) ms, which we find to be independent of N inside the uncertainties (see Fig. S3).We cross-checked the value of γ s by repeating the same experiment in a small Mott insulator with ≈ 1 × 10 3 atoms and no doublons (inset of Fig. S3).For a pulse duration of 10 ms, we remove all the doublons within our detection sensitivity and about 5 % of the singlons.We take this into account when calculating the theory curves in Fig. 2. D. Characterization of the potential

Trap frequencies
Across the region sampled by the atomic cloud, the confinement due to the Gaussian lattice beams is approximately harmonic, such that µ(x, y) = µ 0 − 1 2 κ x x 2 + κ y y 2 .The two coefficients are determined by a fit of the EoS for each lattice depth with free parameters T /t, µ 0 /t, κ x d 2 /t, κ y d 2 /t.We verify that the fourth order anharmonic correction coming from the Gaussian profile, calculated for the measured beam waist and power of the beams, is negligible for the size of the cloud.However, the values obtained for κ x and κ y do not fully agree with the independent measurement of the frequency of the oscillatory motion of the atoms in the combined potential.More specifically, we load a spin-polarized cloud in the combined potential of the vertical lattice and one in-plane lattice at a time.We measure the center of mass (COM) oscillation after an initial displacement of about 5 d along the lattice direction and we verify that the Rayleigh range contribution to the combined potential is negligible.A comparison between the two sets of values can be found in Fig. S4.

Fit comparisons
A possible explanation for the discrepancy in the shape of the inferred potential described in Sec.S.I D 1 might be a lack of adiabaticity and equilibration during the tuning of U/t.However, the temperature measured with the fluctuation-dissipation theorem matches the one obtained from the fit of the EoS.This is not the case if we use the chemical potential coming from the oscillatorymotion calibration.In particular, we verify that for this potential, N = 6 and U/t = 33(2) we obtain a temperature inside the regime of convergence of our theory models when we use the FDT, but which does not quan- titatively match the temperature returned by the fit of the EoS with the same potential.
Moreover, we verify that the values of κ x and κ y determined by the EoS fit are robust against atom number variation (see Fig. S5).
Finally, we observe that κ x and κ y are also insensitive to the lattice ramp time duration.In Fig. S6 we show a test which consists of tuning the lattice depth of an SU (6) sample from zero to 13 E rec (U/t ≈ 44) with a linear ramp of different durations ∆t between 300 ms and 1 s.For each ∆t, we fit the EoS with secondorder high-temperature series expansion (HTSE-2) and leave (T /t, κ x d 2 /t, κ y d 2 /t, µ 0 /t) as free fit parameters.We observe that the values of κ x and κ y returned by the fit for different ∆t are compatible among each other [κ x d 2 /t = 0.151(2), κ y d 2 /t = 0.204(3)].However, they are incompatible with the ones predicted by the independent "oscillatory-motion" calibration described in Sec.S.I D 1 (κ osc x d 2 /t ≃ 0.109, κ osc y d 2 /t ≃ 0.131).If we use these values for the fit, it fails for ∆t = 0.3 s and 0.5 s and it returns high residuals for ∆t = 1 s, failing to reproduce the cloud shape especially in the center [Fig.S6(c)].We conclude that if the mismatch between the two different harmonic potential parameters comes from a lack of equilibration during the tuning of the lattice depth, the time scale to achieve this equilibration significantly exceeds the experimentally accessible timescales.

Anharmonicities
We characterize the anharmonic corrections of the potential by an ensemble with a high temperature T ≫ t into the lattice and applying the atomic limit model with t = 0.In this way, we can generate a spatial map of the chemical potential and observe that the functional modeling as a harmonic trap is a good approximation.In particular, we estimate the anharmonic corrections being less than 0.03 µ 0 over the whole region of interest in the deep Mott insulating regime.

E. Imaging techniques
We use saturated in situ absorption imaging on the stretched 1 S 0 → 1 P 1 transition with circularly polarized light in the presence of a small magnetic field offset ≃ 1 G and a pulse duration of 5 µs.We determine the density with the modified Lambert-Beer law [6], which accounts for the high-intensity effects: where n(x, y) is the density at pixel position (x, y) and I in = I in (x, y) and I out = I out (x, y) are respectively the incident light and the light after absorption.We calibrate the effective saturation intensity I eff sat by varying I in /I sat between 2 and 8 (where I sat = πhcΓ/(3λ 3 ) ≃ 60 mW/cm 2 is the saturation intensity of the transition with wavelength λ and linewidth Γ) and minimizing the variation of the density profile as described in Ref. [6].We find I eff sat /I sat = 3.0(2) with variations of less than 5% between spin mixtures.We obtain compatible values in the 3D dipole trap, in the 2D bulk and in the deep Mott insulating regime with variations smaller than 5%.The effective cross section σ is calibrated for each spin mixture by determining the minimum of the compressibility κ = ∂n/∂µ near the insulating regime at nd 2 ≃ 1 (see also Fig. S10).We obtain σ/σ 0 = [0.310(3),0.320(3), 0.321(3)] respectively for N = [3,4,6], where σ 0 = 3λ 2 /(2π) is the photon resonant cross section.As an independent measurement, we extract the effective cross section from the density shot noise of a thermal sample in the 2D bulk according to the method described in Ref. [7].In particular we make use of the fact that ⟨|δn 0 (k)| 2 ⟩ is proportional to σN p M 2 (k), where δn 0 are the measured local density fluctuations per pixel, k is a vector in Fourier space, N p is the total atom number and M(k) is the modulation transfer function of the imaging system.By modeling M with aperture, intensity attenuation and aberrations as free fit parameters, we determine both σ and M [7].We obtain σ/σ 0 = [0.35(1),0.38(1), 0.384 (7)] respectively for N = [3,4,6].We attribute the differences between the values to cooperative optical response effects at high densities [8,9].

PSF modeling
The point spread function (PSF) as reconstructed from M [7] is shown in Fig. S7.The imaging imperfections are taken into account by convolving the 2D profile of the density predicted by the theory with the PSF.We estimate the systematic error in the determination of the EoS by varying the HWHM of the reconstructed PSF.
Testing the sensitivity of the fit parameter results, we find that, in the deep Mott insulating regime, where the effects of the PSF are most relevant, a variation of 50 % in the size of the HWHM causes only a change of 4 % in the entropy [for U/t = 33(2) and N = 6].

F. Fit method and parameters for the EoS
We fit a 2D model of the density with fixed κ x , κ y and cross section determined with a separate fit as described in Sec.S.I D 1 and S.I E. Due to small position fluctuations during imaging, we perform an initial fit of each image to determine the center of the trap.We use this to align the centers of the distributions and perform a fit of the average of several realizations with similar atom number (15 realizations for the dataset of Fig. 2).
In Tab.S2 we report the values returned by the EoS fit shown in Fig. 2 and a comparison between different numerical methods.
The fit results do not take into account the uncertainties on U/t.By fitting different values U/t to the same data, we estimate the systematic error on the entropy per particle and the temperature to be about 1% and 0.015 U respectively.

G. Equation of state for higher temperatures
In Fig. S8 we probe the equation of state as a function the initial entropy for N = 6.We increase the entropy by holding the atoms in the 2D bulk before the lattice ramps.In this way, we can vary the entropy per particle s between approximately 1 k B and 2 k B (grey squares).After loading into the lattices, we fit the density and compare different numerical methods, including DQMC, NLCE and HTSE.For high initial entropy and large interactions, the methods agree very well with each other.For lower initial entropy, first HTSE and then NLCE begin to deviate.We also perform a round-trip experiment as described in the main text (grey crosses).For large initial entropies, we observe a good agreement between the entropies in the lattice and after the round trip experiment, indicating that the increase is mainly occurring during the ramp up.

H. Equation of state in the bulk
The entropy per particle in 2D prior to turning on the in-plane lattices is calculated from a fit of the equation of state for a weakly-interacting Fermi gas.In particular, we obtain the temperature from the fit of the in situ density according to: where Ω is the grand potential, Li 1 is the polylogarithm of order 1, β = 1/(k B T ), and m is the mass of one 173 Yb atom.The harmonic potential is taken into account in local density approximation.The effect of interactions is taken into account by introducing a correction to the chemical potential [10]: where E F is the Fermi energy and g is the interaction parameter: Here, k F = 4πn/N is the local Fermi vector and a 2D ≃ 2.5 × 10 −3 a 0 is the 2D scattering length.The density-dependency of the interaction parameter is evaluated self-consistently in the fit routine similarly to what has been done in 3D in Ref. [11].We convolve the predicted density profile with the PSF to take into account imaging imperfections.The entropy is calculated from the temperature returned by the fit and the spectrum of the 3D non-isotropic harmonic oscillator [12].We determine a correction to the entropy due to the interactions between 5 % and 20 % compared to the one returned by the fit of a non-interacting Fermi profile.Fit parameters for the EoS shown in Fig. 2 and comparison between different methods.For each sample, we consider 15 shots postselected according to the total atom number.For HTSE and NLCE methods, the number indicates the order.The agreement between two consecutive orders indicates that NLCE has converged.In Fig. 2 we plot the highest NLCE order if DQMC otherwise.The uncertainties correspond to the values returned by the fit and do not take into account additional systematic errors.

Binning, PSF and shot noise
We spatially bin each frame in square "superpixels" and compute the chemical potential and the local density variance in each of them.
We take into account the PSF originating from the binning by measuring the density fluctuations for a SU (6) thermal cloud (T ≫ T F ) in the bulk with the same binning.For a non-interacting cloud with point-like PSF and T /T F ≫ 1, var(n)/⟨n⟩ ≃ 1.The binning lowers this ratio, which can be used as a calibration scaling factor for the fluctuations in the lattice.However, we also take into account additional corrections due to FIG.S8.Equation of state for high temperature for N = 6.By varying the hold time in the 2D bulk before ramping up the lattices, we increase the initial entropy per particle (grey squares).After loading into the lattice, we perform a fit of the density obtained by the average of 15 frames with Np ≈ 2.2 × 10 3 and a s.e.m. of 15 to 35 for each combination of U/t and hold time, after atom number postselection and COM alignment.We fit the same datapoint with different methods and compare the results.We fit DQMC (blue circles), NLCE of order 4 (orange diamonds) and HTSE of order 2 (green triangles).After loading into the lattice we perform a roundtrip experiment by symmetrically ramping down the lattices and measuring the entropy in the bulk (grey crosses).For this measurement, the vertical bandgap was h × 3.63(1) kHz and the confinement's parameters κx and κy have been calibrated according to the same method presented in the main text (for N = 6 only).
the temperature and the interactions in the bulk.We do so by fitting the cloud with a weakly-interacting model (see Sec. S.I H) and determine T /T F = 1.02 (1).For this temperature and interaction strength, we expect var(n)/⟨n⟩ ≡ ξ ∞ = 0.69 (1) for small densities (see Fig. S9).The binning to finite-size superpixels introduces an additional correction ξ ∞ → ξ i , where i is the linear size of the superpixel.For 4 px × 4 px superpixels, the size that we use in Fig. 3, we measure the correction ξ 4 = 0.39 (1).
In the lattice, we first do the binning, then subtract the photon shot noise as the average density at the edge of the region of interest and finally rescale the variance by the factor ξ 4 σ bulk /σ lat −1 , which also takes into account the correction to the cross section described in Sec.S.I E.

Compressibility
To calculate the isothermal compressibility κ = ∂n/∂µ| T , we first bin the experimental data in 1D as n(µ) and then numerically perform a 3-point differentiation.Fig. S10 shows the compressibility for the N = 6 dataset of Fig. 3 and a comparison with the datasets for N = 3 and 4 shown in Fig. 2.

FDT Thermometry
In order to extract the temperature from the FDT, we linearly interpolate the compressibility (Sec.S.I I 2) to match the binning grid of the density fluctuations (Sec.S.I I 1) and calculate the local temperature as the ratio between the local variance and compressibility.The global temperature T FDT is calculated as weighted average of the local temperature as a function of the chemical potential for nd 2 > 0.05.In Fig. S11, we show the local temperature and compare it to T FDT and T EoS for N = 6.Isothermal compressibility κ = ∂n/∂µ| T for N = 3 (red), N = 4 (purple), and N = 6 (blue).κ0 is the compressibility of a SU(6) non-interacting gas at zero temperature.For each combination of N and U/t, we numerically differentiate the density with respect to the chemical potential.The dataset is the same as Figs. 2 and 3. Points correspond to the differentiation of the experimental data and lines correspond to the differentiation of the theoretical curves.

S.II. NUMERICAL METHODS
In this section we present details of the DQMC and NLCE calculations used for results in the main text, including estimates of systematic errors, and a derivation of the HTSE.

A. Determinantal Quantum Monte Carlo
Averages of the thermal equilibrium observables are evaluated with Determinantal Quantum Monte Carlo (DQMC) [13,14] on 6 × 6 lattices by introducing N (N − 1)/2 auxiliary Hubbard-Stratonovich fields, one for each interaction term [15,16] 1 .In this method, the inverse temperature β is discretized in steps of ∆τ (Trotter steps).Results in the main text use ∆τ = 0.05/t.We obtain DQMC data for 5 different random initial seeds for U/t = 7.43, 10.38 and for 20 different random seeds  for U/t = 2.34.For each Monte Carlo trajectory we perform 2000 warm up sweeps through the whole space-time lattice and 8000 sweeps for measurements.Furthermore, we use 4 global moves per sweep to ensure ergodicity [21].
In addition to statistical errors, (controllable) systematic errors can arise from the finite Trotter step, finitesize effects, and imperfect equilibration of the Monte Carlo algorithm before measurement.The Trotter error obtained by comparing ∆τ = 0.04/t and ∆τ = 0.05/t results is < ∼ 0.03 t for the temperature and < ∼ 0.06 k B for the entropy across all N and U/t considered.These correspond to relative errors that are < ∼ 8%, and in all cases are smaller than the statistical error bars 2 .The finitesize error assessed by taking the difference of temperature and entropy obtained from fits using results from 4 × 4 and 6 × 6 lattices is < ∼ 0.04 t for the temperature and < ∼ 0.03 k B for the entropy for all N and U/t considered, except for N = 6 at U/t = 2.34 where the errors are 0.07 t and 0.13 k B .The equilibration error is negligible in the homogeneous case and was estimated by comparing results using 2000 and 8000 warm up sweeps.This error is < ∼ 0.002 for the density and < ∼ 0.04 t for the energy across all chemical potentials considered in the homogeneous case.
DQMC results are calculated on a grid of µ and T , which can also introduce systematic errors into the fits of the experimental equation of state measurements to theory.We use a grid with dµ = 0.25 t and dT /t ∼ 0.05 for T /t ≤ 0.6, and dT /t ∼ 0.1 for 0.6 ≤ T /t ≤ 1.5.Effects due to the coarseness of the µ grid were estimated by considering differences of temperature and entropy fits using dµ = 0.25 t and dµ = 0.5 t and are < ∼ 0.007 t and < ∼ 0.05 k B .
2 For results at U/t = 7.43 and 10.38, relative errors are ∼ 1%.Convergence criteria correspond to the lowest T /t for which fluctuations in entropy for two consecutive µ data points still fall within statistical error bars.

B. Numerical Linked Cluster Expansion
Thermodynamic observables are computed using site expansion Numerical Linked Cluster Expansion (NLCE) up to seventh order (seven sites) for SU(3), fifth order for SU(4) and fourth order for SU(6) Fermi-Hubbard models.A brief discussion of the algorithm is presented below.Extensive properties of a lattice L are computed by performing a weighted sum of the property values over all embeddable clusters c [22,23] where P (c) is computed by performing exact diagonalization of the Hamiltonian in Eq. 1 defined over the cluster c.The Hilbert space dimension grows exponentially with both system size N s and number of spin flavors N .We use spin flavor conservation symmetry to reduce the Hilbert space dimension as mentioned in Ref. [15].Note that Eq.S.8 follows from Eq. S.9 applied to c = L.The NLCE works by truncating Eq.S.8 to a finite order (finite number of sites in the clusters included), which yields accurate results when correlation lengths are sufficiently short, as happens when temperature is not too low.The NLCE is typically much more accurate than exact diagonalization employing the same number of sites.See the Appendix of Ref. [15] for a comparison of bare ED and NLCE for the SU(N ) Hubbard model.The NLCE converges to lower temperatures when the density is an integer, and when U/t is larger.For both cases, the sys-tem has shorter-ranged correlations at a fixed temperature and thus the system properties are captured by small clusters.
The only error in the NLCE arises from truncating Eq.S.8 to finite cluster size.We evaluate this by comparing properties calculated with successive orders of the NLCE approximation.When two orders of NLCE consistently agree, generally the result is converged and the value equals that in the thermodynamic limit [23].
The NLCE was computed over T /t grid of 500 points linearly spaced in the range [0.5, 10], i.e. dT /t ∼ 0.02.The µ/t grid was varied for different U/t and N such that the density, ⟨n⟩ ∈ [0, N/2].For each U/t and N , there are 500 µ/t-points.In the main text (Fig. 2), the chemical potential is re-scaled with U , thus the µ-spacing dµ/U varies between 0.01 and 0.02.The choice of the temperature and chemical potential grids was made to get a dense enough grid to get small fit errors to the experimental data (Fig. 2).

C. High Temperature Series Expansion
We also present results from a second-order high temperature series expansion (HTSE) in t/T which is accurate for T > ∼ t.In the following, subscripts are used to indicate the expansion's order in t/T for the different physical quantities ⟨O⟩ ℓ .
1. Zeroth-order HTSE (t = 0) The energy as a function of the occupation of a single site in the grand canonical ensemble is given by,

FIG. 1 .
FIG. 1. Probing the 2D SU(N ) Fermi-Hubbard model with ultracold atoms.(a) The 1 S0 ground state of 173 Yb naturally features an SU(6) symmetry, which can be freely tuned to N ≤ 6 by preparing a suitable combination of the nuclear spin states mF = −5/2, −3/2, . . ., +5/2 (colored circles and arrows).(b) Schematic of the experimental setup showing a gas in a single layer 2D square lattice with harmonic confinement detected with absorption imaging along gravity.(c) Spatial distribution of the density n(x, y) for N = 6.Each cloud shown in the horizontal frame has been prepared with the same initial entropy in the bulk and loaded into the lattice to a different U/t value.(d) Singly-occupied sites after parity projection.Each horizontal frame corresponds to the same state shown in the same column of (c).(e) Density profiles for the data shown in (c) and (d) along the corresponding dashed lines in the left frames.Each image was produced using the averaging of 8 shots after center of mass alignment [65].

FIG. 3 .
FIG. 3. (a) Measured density fluctuations (blue) for N = 6 as a function of the chemical potential for different interaction strengths.The data points have been obtained from the variance of 15 frames (same as Fig. 2) computed on spatiallybinned probe areas of size ≈ 5.1 × 5.1 d 2 (4 × 4 square camera pixels).The photon shot noise has been subtracted and a PSF correction has been taken into account [65].The green line corresponds to the numerically-differentiated compressibility κ times the temperature TEoS obtained from the EoS-fit of the averaged data, while the black dashed line corresponds to the theory-derived compressibility times TEoS.The vertical line indicates µ(nd 2 = 1).The grey dashed line corresponds to the on-site density fluctuations δn 2 0 = ⟨n 2 ⟩ − ⟨n⟩ 2 calculated with NLCE for TEoS.(b) Comparison of the temperatures TFDT (dark blue diamonds) and TEoS (light blue hexagons).Error bars are the s.e.m.
FIG. S2.Direct measurement of U with modulation spectroscopy for U/t ≈ 178.Red line: expected value of U according to the Wannier overlap.Black line: fit of a Lorentzian function.To enhance the signal-to-noise ratio, data have been evaluated in an elliptical shell where the cloud is mainly in the insulating phase (inset).
FIG.S3.Calibration of the photoassociation beam.Atom number as a function of the pulse duration for N = 3 (red), N = 4 (green) and N = 6 (blue).We first tune the lattice depth to 12 Erec [U/t = 33(2)] and then quench it to 30 Erec before applying the photoassociation pulse.Inset: Same experiment repeated for a smaller atom number, such that nd 2 ≃ 1 in the center.
FIG.S4.Coefficients of the harmonic confinement along the lattice axes (x, y).Circles correspond to the harmonic confinement according to the EoS fit described in the text.Diamonds correspond to an independent calibration measuring the frequencies of the COM oscillations around the center of the trap.Inset: example of COM oscillations for U/t = 33(2).
FIG. S6.Robustness of the EoS fit as a function of the ramp duration.(a) We tune a SU(6) sample to U/t ≈ 44 with a linear ramp of duration ∆t = 0.3 s (green), 0.5 s (red), 1 s (purple).(b) For each ramp duration ∆t, we fit the density profile assuming an harmonic potential with coefficients determined by the "oscillatory-motion" calibration (orange) and by leaving the coefficients free in the EoS fit (blue).Here, we plot results of the fit with respect to the spacial radial profile.For the first two frames, the "oscillatory-motion" fit fails.(c) Comparison of the harmonic potential coefficients.
FIG. S7.(a) Reconstructed point spread function (PSF).x and y have been labeled according to the convention defined in Fig. 1.The color scale is in arbitrary units and varies between -1 (blue) and 1 (red).(b) Cuts along the main axes of the PSF [dashed lines in (a)].
FIG. S9.Fluctuations in the bulk for N = 6.Left: density variance calculated with superpixels of size 8 px × 8 px.Blue points: measured density variance.Grey line: for a thermal, non-interacting cloud var(n)/n = 1.Purple line: theoretical prediction estimated from the temperature returned by a non-interacting Fermi fit [T /TF = 1.37(1)].Green line: theoretical prediction according to a weakly-interacting Fermi fit [T /TF = 1.02(1)].Orange line: linear fit of the measured density fluctuations.Right: slope nvar/n = ξi returned by the linear fit as a function of the superpixel size i [(i × i) px 2 is the superpixel area].Star = 4): superpixel size used in the lattice.

4 .
Fluctuations for N = 3 and 4In Fig.S12and S13 we show the density fluctuations and the comparison between T EoS and T FDT for N = 4 and 3 respectively.In the case of N = 4, we show a dataset with a larger number of frames compared to what has been used in Fig.3(35 frames instead of 15) with similar s.e.m. in the total atom number to narrow down the error bar on T FDT and allow for a more precise comparison.
FIG. S10.Isothermal compressibility κ = ∂n/∂µ| T for N = 3 (red), N = 4 (purple), and N = 6 (blue).κ0 is the compressibility of a SU(6) non-interacting gas at zero temperature.For each combination of N and U/t, we numerically differentiate the density with respect to the chemical potential.The dataset is the same as Figs.2 and 3. Points correspond to the differentiation of the experimental data and lines correspond to the differentiation of the theoretical curves.
FIG. S12.(a) Measured density fluctuations (purple) for N = 4 as a function of the chemical potential for different interaction strengths.The data points have been obtained from the variance of 35 frames according to the same binning and evaluation method of Fig. 3.The green line corresponds to the numerically-differentiated compressibility κ times the temperature TEoS obtained from the EoS-fit of the averaged data.The vertical line indicates µ(nd 2 = 1).The dashed line corresponds to the on-site density fluctuations δn 2 0 = ⟨n 2 ⟩ − ⟨n⟩ 2 calculated with NLCE for TEoS.(b) Comparison of the temperatures TFDT (dark purple diamonds) and TEoS (light purple squares).Error bars are the s.e.m.
FIG. S13.(a) Measured density fluctuations (red) for N = 3 as a function of the chemical potential for different interaction strengths.The data points have been obtained from the variance of 15 frames according to the same binning and evaluation method of Fig. 3.The green line corresponds to the numerically-differentiated compressibility κ times the temperature TEoS obtained from the EoS-fit of the averaged data.The vertical line indicates µ(nd 2 = 1).The dashed line corresponds to the on-site density fluctuations δn 2 = 2 ⟩ − ⟨n⟩ 2 calculated with NLCE for TEoS.(b) Comparison of the temperatures TFDT (dark red diamonds) and TEoS (light red triangles).Error bars are the s.e.m.
. The values of the relevant Hubbard parameters for this work are re- FIG. S1.Time dependence of the lattice depth V and the Hubbard parameters U , κx, κy and t during the lattice ramps.

TABLE S2 .
FIG. S11.Local temperature for the N = 6 dataset shown in Fig.3.Blue line: temperature according to the fit of the EoS.Red line: temperature according to the fit of the FDT.Vertical line: µ(nd 2 = 1).The local temperature for small and large densities is in general comparable, hinting at global thermal equilibrium.