Bad metallic transport in a cold atom Fermi-Hubbard system
Peter T. Brown Debayan Mitra Elmer Guardado-Sanchez Reza Nourafkan Alexis Reymbaut Charles-David Hébert Simon Bergeron A. -M. S. Tremblay Jure Kokalj David A. Huse Peter Schauss Waseem S. Bakr





Department of Physics, Princeton University, Princeton, New Jersey 08544 USA

Department of Physics, Princeton University, Princeton, New Jersey 08544 USA

Department of Physics, Princeton University, Princeton, New Jersey 08544 USA

Département de Physique, Institut Quantique, and Regroupement Québécois sur les Matériaux de Pointe, Université de Sherbrooke, Sherbrooke, Québec, Canada J1K 2R1

Département de Physique, Institut Quantique, and Regroupement Québécois sur les Matériaux de Pointe, Université de Sherbrooke, Sherbrooke, Québec, Canada J1K 2R1

Département de Physique, Institut Quantique, and Regroupement Québécois sur les Matériaux de Pointe, Université de Sherbrooke, Sherbrooke, Québec, Canada J1K 2R1

Département de Physique, Institut Quantique, and Regroupement Québécois sur les Matériaux de Pointe, Université de Sherbrooke, Sherbrooke, Québec, Canada J1K 2R1

Département de Physique, Institut Quantique, and Regroupement Québécois sur les Matériaux de Pointe, Université de Sherbrooke, Sherbrooke, Québec, Canada J1K 2R1
Canadian Institute for Advanced Research, Toronto, Ontario, Canada M5G 1Z8

Faculty of Civil and Geodetic Engineering, University of Ljubljana, SI-1000 Ljubljana, Slovenia
Jožef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia

Department of Physics, Princeton University, Princeton, New Jersey 08544 USA

Department of Physics, Princeton University, Princeton, New Jersey 08544 USA
Present address: Department of Physics, University of Virginia, Charlottesville, Virginia 22904 USA

wbakr@princeton.edu
Department of Physics, Princeton University, Princeton, New Jersey 08544 USA




Strong interactions in many-body quantum systems complicate the interpretation of charge transport in such materials. To shed light on this problem, we study transport in a clean quantum system: ultracold ^6Li in a 2D optical lattice, a testing ground for strong interaction physics in the Fermi-Hubbard model. We determine the diffusion constant by measuring the relaxation of an imposed density modulation and modeling its decay hydrodynamically. The diffusion constant is converted to a resistivity using the Nernst-Einstein relation. That resistivity exhibits a linear temperature dependence and shows no evidence of saturation, two characteristic signatures of a bad metal. The techniques we develop here may be applied to measurements of other transport quantities, including the optical conductivity and thermopower.


Bad metallic transport in a cold atom Fermi-Hubbard system
    Waseem S. Bakr
    December 30, 2023
==========================================================







In conventional materials, charge is carried by quasiparticles and conductivity is understood as a current of these charge carriers developed in response to an external field. For the conductivity to be finite, the charge carriers must be able to relax their momentum through scattering. The Boltzmann kinetic equation in conjunction with Fermi liquid theory provides a detailed description of transport in conventional materials, including two trademarks of resistivity. The first is the Fermi liquid prediction that the temperature-dependent resistivity ρ(T) should scale like T^2 at low temperature <cit.>. The second is that the resistivity should not exceed a maximum value ρ_max, obtained from the Drude relation assuming the Mott-Ioffe-Regel (MIR) limit which states that the mean free path of a quasiparticle cannot be less than the lattice spacing <cit.>. This resistivity bound itself is sometimes referred to as the MIR limit.

Strong interactions can however lead to a breakdown of Fermi liquid theory. One signal of this breakdown is anomalous scaling of ρ with temperature, including the linear scaling observed in the strange metal state of the cuprates <cit.> and other anomalous scalings in d- and f- electron materials <cit.>. Another is the violation of the resistivity bound ρ < ρ_max, which is observed in a wide variety of materials <cit.>. Additionally, interactions may lead to a situation where the momentum relaxation rate alone does not determine the conductivity, in contrast to the semiclassical Drude formula, generalizations of which hold for a large class of systems called coherent metals <cit.>. Approaches introduced to understand these anomalous behaviors include hidden Fermi liquids <cit.>, marginal Fermi liquids <cit.>, proximity to quantum critical points <cit.> and associated holographic approaches <cit.>, and many numerical studies of model systems, most notably the Hubbard <cit.> and t-J <cit.> models.

Disentangling strong interaction physics from other effects, such as impurities and electron-phonon coupling, is difficult in real materials. Cold atom systems are free of these complications, but transport experiments are challenging due to the finite and isolated nature of these systems. Most fermionic charge transport experiments have focused on either studying mass flow through optically structured mesoscopic devices <cit.> or bulk transport in lattice systems <cit.>. Here, we explore bulk transport in a Fermi-Hubbard system by studying charge diffusion, which is a microscopic process related to conductivity through the Nernst-Einstein equation σ = χ_c D, where D is the diffusion constant and χ_c = . (∂ n/∂μ) |_T is the compressibility. This relation requires only the assumption of linear response and absence of thermoelectric coupling <cit.> and does not rest on assumptions concerning quasiparticles.











We realize the 2D Fermi Hubbard model using a degenerate spin-balanced mixture of two hyperfine ground states of ^6Li in an optical lattice <cit.>. Our lattice beams produce a harmonic trapping potential, which leads to a varying atomic density in the trap. To obtain a system with uniform density, we flatten our trapping potential over an elliptical region of mean diameter 30 sites using a repulsive potential created with a spatial light modulator. We superimpose an additional sinusoidal potential that varies slowly along one direction of the lattice with a controllable wavelength (Fig. <ref>, A and B). By adiabatically loading the gas into these potentials, we prepare a Hubbard system in thermal equilibrium with a small amplitude (typically 10%)  sinusoidal density modulation. The average density in the region with the flattened potential is the same with and without the sinusoidal potential. Next, we suddenly turn off the added sinusoidal potential and observe the decay of the density pattern versus time (Fig. <ref>, C and D), always keeping the optical lattice at fixed intensity. We measure the density of a single spin component, n_↑ using techniques described in <cit.>, giving us access to the total density through n=2n_↑. 

We work at average total density n = 0.82(2). This value is close to a conjectured quantum critical point in the Hubbard model <cit.>. Our lattice depth is 6.9(2) E_R, where E_R/h = 14.66 kHz is the lattice recoil, leading to a tunneling rate of t/h = 925(10) Hz. Here h is Planck's constant. We adjust the scattering length, a_s = 1070(10) a_o, by working at a magnetic bias field of 616.0(2)G, in the vicinity of the Feshbach resonance centered near 690G. These parameters lead to an on-site interaction to tunneling ratio U/t = 7.4(8), which is in the strong-interaction regime and close to the value that maximizes antiferromagnetic correlations at half-filling <cit.>.











We observe the decay of the initial sinusoidal density pattern over a period of a few tunneling times. The short timescale ensures that the observed dynamics are not affected by the inhomogeneous density outside of the central flattened region of the trap. To obtain better statistics, we apply the sinusoidal modulation along one dimension and average along the other direction (Fig. <ref>, A and C). We fit the average modulation profile to a sinusoid, where the phase and frequency are fixed by the initial pattern (Fig. <ref>A). The time dependence of the amplitude of the sinusoid quantifies the decay of the density modulation, Fig. <ref>B. Our experimental technique is analogous to that of <cit.>, which studied the decay of a sinusoidally modulated spin pattern in a bosonic system.

The decay of the sinusoidal density pattern versus the wavelength of the modulation becomes consistent with diffusive transport at long wavelengths. In diffusive transport, the amplitude of a density pattern at wave vector k = 2π/λ will decay exponentially with time constant τ = 1/Dk^2, where D is the diffusion constant. We observe exponentially decaying amplitudes with diffusive scaling for wavelengths longer than 15 sites. However, the decay curves are flat at early times, showing clear deviation from exponential decay. For short wavelengths, we observe deviations from diffusive behavior in the form of underdamped oscillations, which can be understood as the damped limit of sound waves. Both of these effects are related to the fact that a density modulation does not instantaneously create a current, as implied by the diffusion equation. Rather, a current requires a finite amount of time to reach an equilibrium value after the creation of a density modulation. 

To unify the description of modulation decay at all wavelengths, we developed a hydrodynamic description that conserves density and has a finite momentum (or current) relaxation rate <cit.>. This approach leads to a differential equation for the density decay,

    ∂_t^2 n + Γ∂_t n + Γ D k^2 n = 0,

where Γ is the momentum-relaxation rate and D is the diffusion constant. This oscillator model crosses over from an underdamped to an overdamped (approximately diffusive) regime at a modulation wavelength 4π√(D/Γ). Instead of assuming that D and Γ are dependent parameters linked through a Drude formula, as would be the case in a system that can be described using quasiparticles, we determine D and Γ from our data at a fixed temperature, simultaneously fitting the amplitude as a function of time for all wavelengths as shown in Fig. <ref>B. 

Our model neglects thermoelectric effects, which affect the measured density response by coupling local energy density modulations, and the resulting temperature gradients, to the particle current. We justify this approximation based on the empirical fact that our simple model fits the data and that we have not been able to detect any measurable temperature modulation in the gas <cit.>. In addition, theoretical work suggests that the thermopower (Seebeck coefficient) is negligible near our doping <cit.>. 











The temperature dependence of D and Γ are the focus of the rest of this paper. The temperature is controlled as follows. After the initial preparation of the cloud, we hold the atoms in the trap or modulate the lattice amplitude for a controlled time to heat the system. To determine the temperature of the cloud after the system has equilibrated, we measure the singles density or local moment, n^s = n_↑ + n_↓ - 2n_↑ n_↓, and the nearest-neighbor correlations between spin-up atoms C^↑(𝐝) = 4 (n_𝐢 + 𝐝,↑ n_𝐢↑ - n_𝐢 + 𝐝,↑n_𝐢↑), where 𝐢 = (i_x, i_y). We compare these quantities to determinantal quantum Monte Carlo (DQMC) simulations to extract the temperature <cit.>. For temperatures at the low end of the range we can access, between 0.3 < T/t < 1, the density correlations are a sensitive thermometer. At higher temperature the singles density becomes a better thermometer. We have compared the temperature of the gas before switching off the potential modulation and after the density modulation has decayed, and find no measurable increase.

As the temperature is lowered, Pauli blocking closes scattering channels, leading to an increased range of diffusion, in agreement with our observations in Fig. <ref>A. At high temperatures, D is expected to saturate, eventually approaching an infinite temperature limiting value <cit.>. The diffusion constant is closely related to the mean-free-path, l, and is often estimated as D = l v /2, where v is the mean quasiparticle velocity <cit.>. Therefore, the MIR limit implies a lower bound on the diffusion constant, D ≳ t a^2 / ħ, where a is the lattice constant. Our measured diffusion constants approach this derived bound at high temperatures, but do not violate it. Because of the difficulty of measuring diffusion constants in materials, this limit has not been tested in real bad metals. We do not compare the measured diffusion constants with theory because determining D requires working in the limit λ→∞ <cit.>, and exact techniques such as diagonalization of finite systems and DQMC are limited to small system sizes. Even determining the infinite temperature limiting value is a non-trivial quantum dynamics problem <cit.>. 

In a clean system like ours momentum relaxation can only occur thanks to umklapp scattering, where a portion of the net momentum in a collision is transferred to the rigid lattice. Nevertheless, the momentum relaxation is strong at our interaction strength which makes determining the temperature dependence of Γ challenging because Γ drops out of the model entirely in the overdamped limit. We find that Γ decreases weakly with decreasing temperature (Fig. <ref>, Inset). This trend may again be understood as Pauli-blocking suppressing momentum relaxation at low temperatures.

We compare the experimental Γ to results from state-of-the-art finite-temperature Lanczos method (FTLM) and dynamical mean-field theory (DMFT) simulations by estimating the momentum relaxation rate as the half-width at half-maximum of the Drude peak in the optical conductivity. The optical conductivity has an additional peak at ω∼ U, but this does not affect Γ significantly <cit.>. Our experimental Γ agrees reasonably with the DMFT results, but exceeds the FTLM results by up to a factor of two. FTLM is an exact technique expected to give correct results at high temperature. One possible explanation for the discrepancy is that Γ is sensitive to the amplitude of the density modulation. To test this, we measured Γ and D versus modulation amplitude <cit.>. We found D is insensitive to the amplitude in the range explored. Γ shows some amplitude dependence but, because of the large error bars we can not conclusively say if this is the source of the discrepancy between experiment and FTLM (Fig. S1).

 To extract a resistivity using the Nernst-Einstein relation, we need the compressibility. It is determined in a separate experiment by measuring the variation of total density versus position in a harmonic trap and converting the position to chemical potential in the local density approximation <cit.>. The measured compressibility increases with decreasing temperature (Fig. <ref>B). For our highest experimental temperatures, χ_c approaches n(1 - n/2)/T, as expected in the high temperature limit <cit.>. At sufficiently low temperature χ_c is expected to saturate, but we do not reach this limit at our lowest experimental temperature, T/t = 0.3. Our experimental results agree well with DQMC numerics over the full range of experimental temperatures.
 
 




 




We can now use the Nernst-Einstein relation to determine the conductivity from the measured diffusion constant and charge compressibility. We examine the temperature dependence of the resistivity ρ = 1/σ in Fig. <ref>, and observe that it rises without limit, showing no sign of saturation. Assuming the existence of quasiparticles, the maximum resistivity obtained from the Drude relation using the MIR limit is ρ < ρ_max≈√(2 π/n)ħ <cit.>. We find that our resistivity violates this bound for temperatures above T/t ∼ 1.3.  The temperature where ρ exceeds this limit is near the Brinkman-Rice temperature scale, defined by T_BR = (1-n) W, where W = 8t is the bandwidth, which is an estimate of the degeneracy temperature of quasiparticles in a doped Mott insulator. Similar violation of the resistivity bound at T_BR has been observed in DMFT studies <cit.>.

The failure of ρ to saturate at the resistivity bound is similar to behavior observed in bad metals at high temperatures <cit.>. In our system, the violation of the resistivity bound is not associated with the mean free path becoming shorter than the lattice spacing because the diffusion constant does not violate its derived bound, but rather with the temperature dependence of the compressibility <cit.>. This suggests a need for a more careful distinction between the MIR limit on the mean free path and the resistivity bound, despite the presumed equivalence of these concepts in condensed matter experiments.

To further elucidate the temperature dependence of ρ, we fit our results to the form ρ(T) = ρ_o + A T + B T^2. We find the temperature dependence is linear to good approximation as we obtain ρ_o = 1.1(1) ħ, A = 1.55(15) ħ/t, and B = 0.03(3) ħ/t^2. Alternatively, a power law fit to the form ρ(T) = ρ_o + (C T)^α yields ρ_o = 1.2(2) ħ, C = 1.4(2) ħ/t, and α = 1.1(1). Similar fits show the inverse diffusion constant 1/D scales with α = 0.6(1) and the inverse charge compressibility scales with α = 0.85(20). In our temperature range, the linear resistivity is a combined result of the temperature dependence of the diffusivity and compressibility, both of which behave in a non-trivial way. This behavior should be contrasted with the high-temperature regime, T ≫ W, where D saturates to a limiting value and the resistivity inherits its temperature dependence from the compressibility, which scales as χ_c ∝ 1/T  <cit.>. It should also be contrasted with the low-temperature regime usually considered in condensed matter where the compressibility has saturated and the resistivity inherits its temperature dependence from the diffusion constant.

We end with more detailed comparison of resistivity with available theories. At our higher experimental temperatures we compare with FTLM, which is an exact technique, and find reasonable agreement (Fig. <ref>). The experimental resistivity is systematically smaller than the FTLM calculation but within error bars. This may be a result of the uncertainty in determining U/t. At our lowest experimental temperatures, FTLM suffers from finite size effects which become relevant as correlation lengths approach the cluster size. For the 4 × 4 site cluster considered here, these effects limit FTLM resistivity calculations to T/t ≳ 1. 

Because our experiment explores low temperatures which are inaccessible to FTLM, we also compare with an approximate technique, single-site DMFT <cit.> (Fig. <ref>). We find the DMFT tends to overestimate the experimental resistivity at high temperatures. At our highest experimental temperatures, the DMFT resistivity is linear with a positive zero-temperature intercept. This linear scaling crosses over to a second linear scaling with a negative zero-temperature intercept around T/t = 2. This second linear region continues down to about T/t = 0.8 where the resistivity acquires a significant quadratic component. These regimes coincide with two different regimes observed in the DMFT compressibility <cit.>. Previous DMFT studies at stronger interaction strengths have also observed these two linear regimes at intermediate temperatures, finding evidence for resilient quasiparticles in the lower temperature regime <cit.>. We do not observe the change of slope in the resistivity expected near T/t = 2 in either the experimental data (within uncertainties) or the FTLM results. This suggests a need for comparison between more refined DMFT and exact theoretical approaches in the regime where this is possible.





Our experiment paves the way for future studies of the optical conductivity and thermopower, which can be examined near equilibrium using a similar approach. Both of these quantities might be expected to show anomalous scalings, as in the cuprates <cit.>. In line with theoretical work such as <cit.>, searching for direct signatures of resilient quasiparticles using spectroscopic techniques <cit.> would also be very interesting. Further experimental studies will also provide important benchmarks for approximate theoretical methods, as the combination of low temperature, finite-doping, and dynamics is challenging for exact theoretical approaches.





§ ACKNOWLEDGEMENTS



We thank Joseph Thywissen, Martin Zwierlein, and Sean Hartnoll for stimulating discussions. We thank Maxime Charlebois and Patrick Sémon for contributions to the continuous time Monte Carlo impurity solver codes. We thank Dominic Bergeron for assistance with the analytic continuation and comparison of results with the Two-Particle-Self-Consistent (TPSC) approach. This work was supported by the NSF (grant no. DMR-1607277), the David and Lucile Packard Foundation (grant no. 2016-65128), the AFOSR Young Investigator Research Program (grant no. FA9550-16-1-0269), the Canada First Research Excellence Fund, the Natural Sciences and Engineering Research Council of Canada (NSERC) under grant RGPIN-2014-04584, the Research Chair in the Theory of Quantum Materials (AMST), and the Slovenian Research Agency Program P1-0044. Simulations were performed on computers provided by the Canadian Foundation for Innovation, the Ministère de l'Éducation des Loisirs et du Sport (Québec), Calcul Québec, and Compute Canada.  W.S.B. was supported by an Alfred P. Sloan Foundation fellowship. P.T.B. was supported by the DoD through the NDSEG Fellowship Program. Data reported in this paper and the code required to reproduce the data analysis are archived on the Open Science Framework and GitHub <cit.>. 






60
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Coleman(2015)]Coleman2015
  author author P. Coleman, 10.1017/CBO9781139020916 title Introduction to Many-Body Physics (publisher
  Cambridge University Press, year 2015)NoStop
[Ioffe and Regel(1960)]Ioffe1960
  author author A. Ioffe and author A. Regel,  https://elibrary.ru/item.asp?id=21765390 journal
  journal Prog. Semicond. volume 4,  pages 237 (year 1960)NoStop
[Mott(1972)]Mott1972
  author author N. F. Mott, 10.1080/14786437208226973 journal
  journal Philos. Mag. volume 26,  pages 1015 (year 1972)NoStop
[Hussey(2008)]Hussey2008
  author author N. E. Hussey, 10.1088/0953-8984/20/12/123201 journal journal J. Phys.: Condens. Matter volume 20, pages 123201 (year
  2008)NoStop
[Stewart(2001)]Stewart2001
  author author G. R. Stewart, 10.1103/revmodphys.73.797 journal journal Rev. Mod. Phys. volume
  73, pages 797 (year 2001)NoStop
[Gunnarsson et al.(2003)Gunnarsson, Calandra, and Han]Gunnarsson2003
  author author O. Gunnarsson, author M. Calandra,  and author J. E. Han, 10.1103/revmodphys.75.1085 journal
  journal Rev. Mod. Phys. volume 75,  pages 1085 (year 2003)NoStop
[Hartnoll(2014)]Hartnoll2014
  author author S. A. Hartnoll, 10.1038/nphys3174 journal
  journal Nat. Phys. volume 11,  pages 54 (year 2014)NoStop
[Anderson(2008)]Anderson2008
  author author P. W. Anderson, 10.1103/PhysRevB.78.174505 journal journal Phys. Rev. B volume
  78, pages 174505 (year 2008)NoStop
[Varma et al.(1989)Varma,
  Littlewood, Schmitt-Rink, Abrahams, and Ruckenstein]Varma1989
  author author C. M. Varma, author P. B. Littlewood, author S. Schmitt-Rink, author E. Abrahams,  and author A. E. Ruckenstein, 10.1103/physrevlett.63.1996 journal journal Phys. Rev. Lett. volume 63, pages 1996 (year
  1989)NoStop
[Vuččiččevi ćć et al.(2015)Vuččiččevi ćć,
  Tanaskovi ćć, Rozenberg, and Dobrosavljevi ćć]Vucicevic2015
  author author J. Vuččiččevi ćć, author
  D. Tanaskovi ćć, author M. J. Rozenberg,  and author V. Dobrosavljevi ćć, 10.1103/PhysRevLett.114.246402 journal journal
  Phys. Rev. Lett. volume 114, pages
  246402 (year 2015)NoStop
[Hartnoll et al.(2018)Hartnoll, Lucas, and Sachdev]Hartnoll2018
  author author S. A. Hartnoll, author A. Lucas,   and author S. Sachdev, https://www.ebook.de/de/product/31890881/sean_a_hartnoll_andrew_lucas_subir_sachdev_holographic_quantum_matter.html
  title Holographic Quantum Matter (publisher MIT Press, year 2018)NoStop
[Scalapino(2007)]Scalapino2007
  author author D. J. Scalapino, title Handbook of high-temperature
  superconductivity,  (publisher Springer New York, year 2007) Chap. chapter Numerical Studies of the 2D
  Hubbard Model, pp. pages 463–493NoStop
[Jakli čč and  Prelov ššek(2000)]Jaklic2000
  author author J. Jakli čč and author
  P. Prelov ššek, 10.1080/000187300243381 journal journal Adv. Phys. volume
  49, pages 1 (year 2000)NoStop
[Brantut et al.(2012)Brantut, Meineke, Stadler, Krinner, and Esslinger]Brantut2012
  author author J.-P. Brantut, author J. Meineke,
  author D. Stadler, author S. Krinner,  and author T. Esslinger, 10.1126/science.1223175 journal journal
  Science volume 337, pages 1069
  (year 2012)NoStop
[Valtolina et al.(2015)Valtolina, Burchianti, Amico, Neri, Xhani, Seman, Trombettoni, Smerzi, Zaccanti,
  Inguscio, and Roati]Valtolina2015
  author author G. Valtolina, author A. Burchianti, author A. Amico,
  author E. Neri, author
  K. Xhani, author J. A. Seman, author A. Trombettoni, author A. Smerzi, author M. Zaccanti, author M. Inguscio,  and author G. Roati, 10.1126/science.aac9725
  journal journal Science volume 350, pages 1505 (year
  2015)NoStop
[Krinner et al.(2017)Krinner, Esslinger, and Brantut]Krinner2017
  author author S. Krinner, author T. Esslinger,
   and author J.-P. Brantut,  10.1088/1361-648x/aa74a1 journal journal J. Phys.: Condens. Matter volume 29,  pages 343003 (year 2017)NoStop
[Lebrat et al.(2018)Lebrat,
  Gri ššins, Husmann, Häusler, Corman, Giamarchi, Brantut, and Esslinger]Lebrat2017
  author author M. Lebrat, author P. Gri ššins, author D. Husmann, author S. Häusler, author L. Corman, author T. Giamarchi, author J.-P.  Brantut,  and author T. Esslinger, 10.1103/PhysRevX.8.011053
  journal journal Phys. Rev. X volume 8, pages 011053 (year
  2018)NoStop
[Ott et al.(2004)Ott,
  de Mirandes, Ferlaino, Roati,
  Modugno, and Inguscio]Ott2004
  author author H. Ott, author E. de Mirandes,
  author F. Ferlaino, author G. Roati, author
  G. Modugno,  and author
  M. Inguscio, 10.1103/PhysRevLett.92.160601 journal journal
  Phys. Rev. Lett. volume 92, pages
  160601 (year 2004)NoStop
[Strohmaier et al.(2007)Strohmaier, Takasu, Günter, Jördens, Köhl, Moritz, and  Esslinger]Strohmaier2007
  author author N. Strohmaier, author Y. Takasu,
  author K. Günter, author R. Jördens, author
  M. Köhl, author H. Moritz,  and author T. Esslinger, 10.1103/PhysRevLett.99.220601 journal journal
  Phys. Rev. Lett. volume 99, pages
  220601 (year 2007)NoStop
[Schneider et al.(2012)Schneider, Hackermüller, Ronzheimer,
  Will, Braun, Best,
  Bloch, Demler, Mandt,
  Rasch, and Rosch]Schneider2012
  author author U. Schneider, author L. Hackermüller, author J. P. Ronzheimer, author S. Will,
  author S. Braun, author T. Best, author
  I. Bloch, author E. Demler, author S. Mandt, author D. Rasch,  and author A. Rosch, 10.1038/nphys2205 journal journal Nat. Phys. volume
  8, pages 213 (year 2012)NoStop
[Xu et al.()Xu,
  McGehee, Morong, and DeMarco]Xu2016
  author author W. Xu, author W. McGehee,
  author W. Morong,  and author B. DeMarco, @noop    http://arxiv.org/abs/arXiv:1606.06669v4 arXiv:1606.06669v4
  NoStop
[Anderson et al.()Anderson,
  Wang, Xu, Venu, Trotzky, Chevy, and Thywissen]Anderson2017
  author author R. Anderson, author F. Wang,
  author P. Xu, author
  V. Venu, author S. Trotzky, author F. Chevy,  and author J. H.  Thywissen, @noop   http://arxiv.org/abs/arXiv:1712.09965v1 arXiv:1712.09965v1 NoStop
[sup()]supplement
  @noop  note Materials and methods are available as
  supplementary materialsNoStop
[Brown et al.(2017)Brown,
  Mitra, Guardado-Sanchez, Schauß, Kondov, Khatami, Paiva, Trivedi, Huse, and Bakr]Brown2017
  author author P. T. Brown, author D. Mitra,
  author E. Guardado-Sanchez,
  author P. Schauß, author S. S. Kondov, author
  E. Khatami, author T. Paiva, author N. Trivedi, author D. A.  Huse,  and author W. S.  Bakr, 10.1126/science.aam7838 journal journal Science volume
  357, pages 1385 (year 2017)NoStop
[Vidhyadhiraja et al.(2009)Vidhyadhiraja, Macridin,  Şen, Jarrell, and  Ma]Vidhyadhiraja2009
  author author N. S. Vidhyadhiraja, author A. Macridin, author C.  Şen, author M. Jarrell,  and author M. Ma, 10.1103/PhysRevLett.102.206407
  journal journal Phys. Rev. Lett. volume 102, pages 206407 (year
  2009)NoStop
[Khatami and Rigol(2011)]Khatami2011
  author author E. Khatami and author M. Rigol, 10.1103/PhysRevA.84.053611 journal
  journal Phys. Rev. A volume 84,  pages 053611 (year 2011)NoStop
[Hild et al.(2014)Hild,
  Fukuhara, Schauß, Zeiher,
  Knap, Demler, Bloch, and  Gross]Hild2014
  author author S. Hild, author T. Fukuhara,
  author P. Schauß, author J. Zeiher, author
  M. Knap, author E. Demler, author I. Bloch,  and author C. Gross, 10.1103/PhysRevLett.113.147205
  journal journal Phys. Rev. Lett. volume 113, pages 147205 (year
  2014)NoStop
[Pálsson and Kotliar(1998)]Palsson1998
  author author G. Pálsson and author G. Kotliar, 10.1103/physrevlett.80.4775 journal journal Phys. Rev. Lett. volume 80, pages 4775 (year
  1998)NoStop
[Perepelitsky et al.(2016)Perepelitsky, Galatas, Mravlje,
   ŽŽitko, Khatami, Shastry, and Georges]Perepelitsky2016
  author author E. Perepelitsky, author A. Galatas, author J. Mravlje,
  author R.  ŽŽitko, author E. Khatami, author B. S. Shastry,  and author A. Georges, 10.1103/PhysRevB.94.235115 journal journal Phys. Rev. B volume
  94, pages 235115 (year 2016)NoStop
[Kokalj(2017)]Kokalj2017
  author author J. Kokalj, 10.1103/PhysRevB.95.041110 journal journal Phys. Rev. B volume
  95, pages 041110 (year 2017)NoStop
[Mukerjee et al.(2006)Mukerjee, Oganesyan, and Huse]Mukerjee2006
  author author S. Mukerjee, author V. Oganesyan,
   and author D. Huse, 10.1103/PhysRevB.73.035113 journal journal Phys. Rev. B volume 73, pages 035113 (year 2006)NoStop
[Leviatan et al.()Leviatan,
  Pollmann, Bardarson, Huse,  and Altman]Leviatan2017
  author author E. Leviatan, author F. Pollmann,
  author J. H. Bardarson, author D. A. Huse,  and author E. Altman, @noop    http://arxiv.org/abs/arXiv:1702.08894v2 arXiv:1702.08894v2
  NoStop
[Ho and Zhou(2009)]Ho2009
  author author T.-L. Ho and author Q. Zhou,  10.1038/nphys1477 journal journal Nat. Phys. volume 6, pages
  131 (year 2009)NoStop
[Cocchi et al.(2016)Cocchi,
  Miller, Drewes, Koschorreck,
  Pertot, Brennecke, and Köhl]Cocchi2016
  author author E. Cocchi, author L. A. Miller,
  author J. H. Drewes, author M. Koschorreck, author
  D. Pertot, author F. Brennecke,  and author M. Köhl, 10.1103/PhysRevLett.116.175301 journal journal
  Phys. Rev. Lett. volume 116, pages
  175301 (year 2016)NoStop
[Deng et al.(2013)Deng,
  Mravlje,  ŽŽitko, Ferrero, Kotliar, and  Georges]Deng2013
  author author X. Deng, author J. Mravlje,
  author R.  ŽŽitko, author M. Ferrero, author G. Kotliar,   and author A. Georges, 10.1103/PhysRevLett.110.086401 journal journal Phys. Rev. Lett. volume 110, pages 086401 (year 2013)NoStop
[Xu et al.(2013)Xu,
  Haule, and Kotliar]Xu2013
  author author W. Xu, author K. Haule,  and  author G. Kotliar, 10.1103/PhysRevLett.111.036401 journal journal Phys. Rev. Lett. volume 111, pages 036401 (year 2013)NoStop
[Kotliar et al.(2006)Kotliar, Savrasov, Haule, Oudovenko, Parcollet, and Marianetti]Kotliar2006
  author author G. Kotliar, author S. Y. Savrasov, author K. Haule,
  author V. S. Oudovenko, author O. Parcollet,  and author C. A. Marianetti, 10.1103/revmodphys.78.865 journal journal Rev. Mod. Phys. volume 78, pages 865 (year 2006)NoStop
[Stewart et al.(2008)Stewart, Gaebler, and Jin]Stewart2008
  author author J. T. Stewart, author J. P. Gaebler,  and author D. S. Jin, 10.1038/nature07172 journal
  journal Nature volume 454, pages 744 (year 2008)NoStop
[Brown et al.(2018)Brown,
  Mitra, Guardado-Sanchez, Nourafkan, Reymbaut, Hébert,
  Bergeron, Tremblay, Kokalj,
  Huse, Schauß, and Bakr]Brown2018a
  author author P. T. Brown, author D. Mitra,
  author E. Guardado-Sanchez,
  author R. Nourafkan, author A. Reymbaut, author
  C.-D. Hébert, author
  S. Bergeron, author
  A.-M. S. Tremblay, author
  J. Kokalj, author D. A. Huse, author P. Schauß,  and author W. S. Bakr, 10.17605/osf.io/3tf8a
  journal journal Open Science Framework 
  (year 2018), 10.17605/osf.io/3tf8aNoStop
[Kokalj(2018)]Kokalj2018
  author author J. Kokalj,   <https://github.com/jurekokalj/ftlm_hub_cond> (publisher
  GitHub, year 2018)NoStop
[Hebert(2018)]Hebert2018
  author author C.-D. Hébert, <https://github.com/ZGCDDoo/ctint-science>  (publisher
  GitHub, year 2018)NoStop
[Nourafkan(2018)]Nourafkan2018
  author author R. Nourafkan, <https://github.com/RezaNourafkan/Transport-> (publisher
  GitHub, year 2018)NoStop
[Floyd and Steinberg(1976)]Floyd1976
  author author R. W. Floyd and author L. Steinberg, @noop  journal journal
  Proc. Soc. Inf. Displays volume 17, pages 75 (year 1976)NoStop
[Kadanoff and Martin(1963)]Kadanoff1963
  author author L. P. Kadanoff and author P. C. Martin, 10.1016/0003-4916(63)90078-2 journal journal Ann. Phys. volume
  24, pages 419 (year 1963)NoStop
[Shastry(2009)]Shastry2009
  author author B. S. Shastry, http://stacks.iop.org/0034-4885/72/i=1/a=016501
  journal journal Rep. Prog. Phys. volume 72, pages 016501 (year
  2009)NoStop
[Peterson and Shastry(2010)]Peterson2010
  author author M. R. Peterson and author B. S. Shastry, 10.1103/PhysRevB.82.195105 journal journal Phys. Rev. B volume
  82, pages 195105 (year 2010)NoStop
[Arsenault et al.(2013)Arsenault, Shastry, Sémon, and  Tremblay]Arsenault2013
  author author L.-F. Arsenault, author B. S. Shastry, author P. Sémon,   and author A.-M. S. Tremblay,  10.1103/PhysRevB.87.035126 journal journal Phys. Rev. B volume 87, pages 035126 (year 2013)NoStop
[Kokalj and McKenzie(2015)]Kokalj2015
  author author J. Kokalj and author R. H. McKenzie, 10.1103/PhysRevB.91.205121 journal journal Phys. Rev. B volume
  91, pages 205121 (year 2015)NoStop
[Varney et al.(2009)Varney,
  Lee, Bai, Chiesa,
  Jarrell, and Scalettar]Varney2009
  author author C. N. Varney, author C.-R. Lee,
  author Z. J. Bai, author S. Chiesa, author
  M. Jarrell,  and author
  R. T. Scalettar, 10.1103/PhysRevB.80.075116 journal journal Phys.
  Rev. B volume 80, pages 075116
  (year 2009)NoStop
[Kokalj and McKenzie(2013)]Kokalj2013
  author author J. Kokalj and author R. H. McKenzie, 10.1103/PhysRevLett.110.206402 journal journal Phys. Rev. Lett. volume 110, pages 206402 (year
  2013)NoStop
[Poilblanc(1991)]Poilblanc1991
  author author D. Poilblanc, 10.1103/physrevb.44.9562 journal journal Phys. Rev. B volume
  44, pages 9562 (year 1991)NoStop
[Georges and Kotliar(1992)]Georges1992
  author author A. Georges and author G. Kotliar, 10.1103/physrevb.45.6479 journal
  journal Phys. Rev. B volume 45,  pages 6479 (year 1992)NoStop
[Kotliar et al.(2001)Kotliar, Savrasov, Pálsson, and  Biroli]Kotliar2001
  author author G. Kotliar, author S. Y. Savrasov, author G. Pálsson,
   and author G. Biroli, 10.1103/PhysRevLett.87.186401 journal journal Phys. Rev. Lett. volume 87, pages 186401 (year 2001)NoStop
[Gull et al.(2008)Gull,
  Werner, Parcollet, and Troyer]Gull2008
  author author E. Gull, author P. Werner,
  author O. Parcollet,  and  author M. Troyer, 10.1209/0295-5075/82/57003 journal journal EPL volume 82, pages 57003
  (year 2008)NoStop
[Werner et al.(2006)Werner,
  Comanac, de' Medici, Troyer,  and Millis]Werner2006
  author author P. Werner, author A. Comanac,
  author L. de' Medici, author M. Troyer,  and author
  A. J. Millis, 10.1103/PhysRevLett.97.076405 journal journal
  Phys. Rev. Lett. volume 97, pages
  076405 (year 2006)NoStop
[Werner and Millis(2006)]Werner2006a
  author author P. Werner and author A. J. Millis, 10.1103/PhysRevB.74.155107 journal journal Phys. Rev. B volume
  74, pages 155107 (year 2006)NoStop
[Haule(2007)]Haule2007a
  author author K. Haule, 10.1103/PhysRevB.75.155113 journal
  journal Phys. Rev. B volume 75,  pages 155113 (year 2007)NoStop
[Gull et al.(2007)Gull,
  Werner, Millis, and Troyer]Gull2007
  author author E. Gull, author P. Werner,
  author A. Millis,  and author M. Troyer, 10.1103/PhysRevB.76.235123 journal journal Phys.
  Rev. B volume 76, pages 235123
  (year 2007)NoStop
[Sémon et al.(2014)Sémon, Yee, Haule, and Tremblay]Semon2014
  author author P. Sémon, author C.-H. Yee,
  author K. Haule,  and author A.-M. S. Tremblay, 10.1103/PhysRevB.90.075149 journal journal Phys. Rev. B volume 90, pages 075149 (year 2014)NoStop
[Caffarel and Krauth(1994)]Caffarel1994
  author author M. Caffarel and author W. Krauth, 10.1103/PhysRevLett.72.1545 journal journal Phys. Rev. Lett. volume 72, pages 1545 (year
  1994)NoStop
[Arsenault and Tremblay(2013)]Arsenault2013a
  author author L.-F. Arsenault and author A.-M. S. Tremblay, 10.1103/PhysRevB.88.205109 journal journal Phys. Rev. B volume
  88, pages 205109 (year 2013)NoStop
[Maldague(1977)]Maldague1977
  author author P. F. Maldague, 10.1103/physrevb.16.2437 journal journal Phys. Rev. B volume
  16, pages 2437 (year 1977)NoStop
[Calandra and Gunnarsson(2003)]Calandra2003
  author author M. Calandra and author O. Gunnarsson, http://stacks.iop.org/0295-5075/61/i=1/a=088
  journal journal EPL volume
  61, pages 88 (year 2003)NoStop






















§ METHODS

We work with an equal spin mixture of ^6Li hyperfine ground states |1⟩ and |3⟩, numbered up from the lowest energy state, which we label as spin up, |↑⟩, and spin down, |↓⟩, respectively. Our system is well described by the Fermi-Hubbard model with the Hamiltonian,


    ℋ - μ N = -t ∑_𝐢, 𝐣, σ(c^†_𝐢, σ c_𝐣, σ + h.c.) + U ∑_𝐢  n_𝐢,↑ n_𝐢,↓ - μ∑_𝐢,σ n_𝐢,σ,

where c^†_𝐢, σ is the creation operator for a fermion on site 𝐢 = (i_x, i_y) with spin σ∈{↑, ↓}, n_𝐢,σ = c^†_𝐢,σ c_𝐢,σ, t is the hopping rate, U is the on-site interaction, and μ is the chemical potential.

The experimental setup and basic parameters are described in detail in the supplement of ref. <cit.>. After preparing a 2D degenerate Fermi gas, we simultaneously load the optical lattice to a final depth of 6.9(2)E_R and the sinusoidally modulated potential with a 50 intensity ramp. We then turn off the sinusoidal modulation in approximately 10 using the spatial light modulator and observe the decay of the density pattern. 

We work at a field of 616.0(2)G. At this field and lattice depth, we find U/t = 7.4(8) from a band structure calculation, which yields t = 925(10), and spectroscopic measurement of U=7.0(7). 


To prepare clouds of variable temperature, we use two different protocols. To reach temperatures in the range T/t = 0.3 - 2, we hold the cloud in the final trapping configuration for variable time. The gas heats due to technical noise at a rate of 3t per second. To reach even hotter temperatures, we modulate the lattice depth at a frequency of 2. To avoid losses, we perform this modulation at 595G where the interaction is weaker. Finally, we ramp the field to its final value, turn on the DMD potential, and follow the same protocol as before.

For low temperatures, our lattice provides all of the radial confinement. For temperatures hotter than T/t ≈ 3, the compressibility of the gas is reduced and we must provide extra confinement to reach appropriate filling. Therefore we increase our trapping frequency using a 1064 beam. 

To image the density of a single spin state, we freeze the motion of the lattice by ramping the lattice depth to 60  E_R in approximately 100. We checked that this ramp effectively freezes the atomic motion by comparing the measured amplitude modulation without turning off the DMD potential and with turning off the DMD potential and then immediately ramping the lattice depth for our shortest wavelength at our lowest temperatures (where the modulation decays fastest). The modulation depths agreed, indicating that the atomic motion is effectively frozen well before the lattice reaches 60  E_R.



§ DMD CALIBRATION

We engineer our deconfining and sinusoidally modulated potential using up to 15 of 650 coherent light derived from a tapered amplifier fed by a diode laser. Our spatial light modulator is a DLP Discovery 4100 with a DLP7000 digital micromirror device (DMD) in an imaging plane configuration. We image this light onto our atoms using two stages of demagnification. First, we demagnify the DMD image by a factor of 5, then we combine the DMD projection path with our imaging path on a dichroic mirror. Our imaging system demagnifies the light by an additional factor of 30. A single DMD micromirror has a pitch of 13.68, so approximately 8 × 8 mirrors determine the potential at a single lattice site. Our imaging system spatially filters the binary image, resulting in a smooth potential at the atoms.

Before each experiment, we load a series of two images into the DMD memory. The first is the sum of a deconfining Gaussian potential for flattening the atomic density in the central part of the cloud with a sinusoidal modulation pattern. The second is only the deconfining potential. The DMD displays these images successively after receiving a trigger. We use the ALP-4.2 API “uninterruptible binary mode” to keep the image on the DMD until the next trigger. The DMD transitions between images in approximately 10. During this time, all mirrors go to the off state, and then the mirrors needed for the next image are turned to the on state. The motion of the mirrors is underdamped, and we observe the mirrors bouncing by measuring diffracted light on a photodiode.

We produce binary images from continuous potential profiles using the Floyd-Steinberg error diffusion algorithm <cit.>.



§ HYDRODYNAMIC MODEL

Hydrodynamics applies at long wavelengths and low frequencies when there are few conserved quantities, typically only mass, momentum, and energy. In most real materials, electrons cannot be treated hydrodynamically because of couplings to phonons and lattice defects which can absorb energy or momentum. For strongly interacting systems with no external couplings, hydrodynamics is applicable. In lattice systems the momentum is not conserved due to umklapp scattering, and only energy and particle number are conserved. In systems with weak umklapp (momentum relaxation) rate, we can also write down a “momentum conservation” equation including this relaxation rate. For a detailed discussion of when strongly interacting systems can be treated hydrodynamically, see <cit.>.
 
The simplest hydrodynamic theory we can write down accounts for conservation of mass, weak relaxation of momentum, and assumes that energy is decoupled from these two. The two equations describing this are
 
    ∂_t n(𝐫, t)    =    -∇·𝐉(𝐫, t) 
    ∂_t 𝐉(𝐫, t)    =    - Γ( D ∇ n(𝐫, t) + 𝐉(𝐫, t) ),

 where the first equation is the continuity equation, and the second equation reduces to Fick's law, 𝐉(𝐫, t) = - D ∇ n(𝐫, t), when ∂_t 𝐉 = 0. In the limit of strong momentum relaxation, a density modulation instantly creates a current satisfying Fick's law. However, for a finite relaxation rate Γ, the current does not instantly follow the changes in the density, and its time lag is described by eq. <ref>. We can alternatively understand eq. <ref> as a momentum “conservation” equation analogous to the Navier-Stokes equation including a weak momentum relaxation rate Γ, zero viscosity, and neglecting terms of higher order in linear response (e.g. terms proportional to the velocity squared).
 
 We spatially Fourier transform these equations and eliminate J to find,
 
    ∂_t^2 n_k + Γ∂_t n_k + Γ Dk^2 n_k = 0,

 which is the equation that appears in the main text. 
 
 This is the same differential equation that describes a damped harmonic oscillator provided we identify the natural frequency ω_o = √(DΓ) k and the damping rate as γ = Γ /2. This oscillator model crosses over from an underdamped to an overdamped regime at a modulation wavelength λ = 4π√(D/Γ). In a system that can be described using quasiparticles, √(D/Γ) is proportional to the mean free path.
 
 In the overdamped limit, Γ≫√(Γ D) k, we recover diffusive behavior n_k(t) ∝ e^-Dk^2t for finite k. In the underdamped limit, we have sound-waves whose amplitude decays at rate Γ/2. If we take the limit k = 0, the current decays exponentially at rate Γ, identifying this as the momentum relaxation rate. The sound-wave and current relaxation rates differ because the sound wave carries both kinetic and potential energy, shared equally, whereas the uniform (k=0) current excitation carries only kinetic energy. As only the kinetic energy is damped, the sound-wave loses energy at half of the rate of the uniform current excitation.
 


 

§ DETAILS OF FITTING THE HYDRODYNAMIC MODEL

 
 For each temperature, we take images at a series of wavevectors k_i, i=1, ..., N at times t_ij, j=1, ..., j_i. For each (k_i, t_ij) we determine the average 2D density profile. We produce a 1D representation of the modulation by averaging density in the central region of the cloud along the direction perpendicular to the modulation and fit the result to a sinusoidal pattern. 
 For a given wavevector, we first fit the shortest time, t_i1=0, modulation pattern with the sinusoidal amplitude, period, phase, and offset as free parameters. For later times, we fix the phase and the period, leaving only the amplitude and the offset as free parameters. From this procedure we extract a series of amplitudes, a(k_i, t_ij) with uncertainties σ_ij. Because the phase is fixed by the first pattern, a negative amplitude a(k_i, t_ij) < 0 indicates that the modulation pattern at time t_ij is 180^∘ out of phase with the initial pattern.

 To compare our measurements with the hydrodynamic model, we write the solutions to eq. <ref> which satisfy the boundary condition ṅ(t=0) = 0,
 
    n(Γ, D, A, k, t) =  A/2(e^ω̃t + e^-ω̃t)e^-Γ t /2,

 where ω̃= √(Γ^2/4 - Γ D k^2), and A is the amplitude. In the underdamped limit eq. <ref> gives a damped cosine. In the overdamped limit it yields a product of a hyperbolic cosine factor and an exponential factor.
 
 Finally, we perform a non-linear least squares fit which minimizes
 
    ∑_i=1^N ∑_j=1^j_i|a(k_i, t_ij) -  n(Γ, D, A_i, k_i, t_ij) |^2/σ_ij^2,

 with the free parameters Γ, D, and A_i, i=1, ..., N. We determine the uncertainty in the fit parameters with a bootstrapping technique.
 


§ LINEAR RESPONSE THEORY


To connect our hydrodynamic model for the density response of our Fermi-Hubbard system to quantities which can be calculated in theory, we consider the effect of perturbing our system with a time and spatial dependent potential, v_𝐢(t) = F(t) sin(𝐤·𝐫_𝐢). In this experiment we suppose F is turned on slowly starting at t = -∞ and switched off suddenly at t = 0, leading to F(t) = e^η tθ(-t), where η parametrizes the slow turn on. If we suppose that H_o is the Fermi-Hubbard Hamiltonian in the absence of this, perturbation, then the full Hamiltonian is H(t) = H_o + H'(t). We can write the perturbation term as

    H'(t)    =    -∑_𝐢 v_𝐢(t) n_𝐢
       =    - ∑_𝐤 v_𝐤(t) n_-𝐤,

where 𝐫_𝐢 is the position of site 𝐢 and n_k is the spatial Fourier transform of the density.

In linear response theory we think of v as the force which is conjugate to the density response. Given a Hamiltonian of the form in eq. <ref>, we can write the density reaction in terms of a response function Φ

    δ n_𝐤(t)   =   ∫_-∞^t dt'  Φ(𝐤, t - t') v_𝐤(t') 
    Φ(𝐤, t - t')    =    -i/ħΘ(t - t') [n_𝐤(t), n_-𝐤(t')],

where we used translational invariance of the unperturbed system, which ensures that only v_𝐤 contributes to the density response at 𝐤. This equation says that the response of the density to the applied field is given by the density correlations, encapsulated in the retarded Green's function Φ.

Fourier transforming eq. <ref> in time and space leads to,

    δ n_𝐤(ω)   =   χ(𝐤, ω) v_𝐤 (ω)  
    χ(𝐤, ω)    =   ∫_-∞^∞ dt'  e^iω t'Φ(𝐤, t') 
       =    -i/ħ∫_0^∞ dt'  e^i(ω + iη) t'Θ(t') [n_𝐤(t'), n_-𝐤(0) ],

where we added an small positive imaginary part, η, to the frequency so that the integral converges. We refer to χ(𝐤, ω) as the density response function. χ is analytic in the upper-half plane, its real part is symmetric in ω, and its imaginary part is antisymmetric in ω. 



§ NERNST-EINSTEIN EQUATION


In our experiment we have direct access to χ(𝐤, ω) because we control the potential which is the generalized force that couples to the density. We want to measure the charge conductivity, σ, which is the response function for the current. The current is conjugate to the vector potential, which we do not control in this experiment. Fortunately, the continuity equation written using linear response relations connects the density response function with the conductivity,

    σ'(ω) = lim_k→0ω/k^2χ”(𝐤,ω),

where σ'(ω) is the real part of the conductivity and χ”(𝐤, ω) is the imaginary part of the density response function.

To apply this expression, we need to know the form of χ. This is provided by our hydrodynamic model. Adding a force term to eq. <ref> and using the definition of the response function from eq. <ref> leads to the expression

    χ(𝐤, ω) = χ_c/1 - i ω/k^2 D - ω^2/k^2 D Γ,

where χ_c is the charge compressibility, D is the diffusion constant, and Γ is the momentum relaxation rate <cit.>. Inserting this expression for χ in eq. <ref>, we find that the optical conductivity has a Lorentzian profile with half-width half-maximum Γ,

    σ'(ω) = χ_c D/1+(ω/Γ)^2,

and the Nernst-Einstein relationship holds for the DC conductivity,


    σ = χ_c D.


The Nernst-Einstein relation is a consequence of eq. <ref> which holds if the density has a diffusive mode at long times and large wave vectors. The Nernst-Einstein relation does not require the exact form for χ in eq. <ref>. We can, for example, add k^2-dependence to D or Γ.
 


§ LINEARITY




To assess the possibility of non-linear effects which are not included in our hydrodynamic model, we varied the initial amplitude of the density modulation at a fixed wavevector. For each curve, we fit a value for Γ and D, to test how the fitted model parameters change with amplitude. The amplitude versus time curves are shown in Fig. <ref>A for λ≈ 12 sites and temperature T/t = 0.4(1). For each initial amplitude, we fit values for Γ and D using our hydrodynamic model. The fit results are shown in Fig. <ref>B,C. We find that the apparent Γ increases with increasing amplitude, and the apparent D is weakly effected by increasing amplitude. To establish an upper bound on the size of this effect, we perform a linear fit to the hydrodynamic parameters versus amplitude and extrapolate a `zero-amplitude' value. We normalize the curve fit parameters by these values in Fig. <ref>B,C. Based on the statistical error in our fit lines, we find that at a typical experimental amplitude of δ n_↑ = 0.07, Γ is increased by a factor of 1.4(4) and D by a factor of 1.06(10). Our extracted values for Γ appear to increase with amplitude, but the statistical error bar is quite large. This is due to the weak dependence of our model on the value of Γ. In the main text, we are able to obtain smaller error bars on Γ by simultaneously fitting decay curves at different modulation wavelengths. That approach is not feasible here because the degree of linearity may depend on modulation wavelength.

A related but distinct type of non-linearity is dependence of the charge compressibility on density. As the total density approaches half-filling, the compressibility decreases. Therefore, the chemical potential modulation we apply tends to decrease the density at the minimum chemical potential values more than it increases the density at the maximum chemical potential values. This can lead to the density modulation deviating from a sine wave. At T/t = 0.4, the compressibility decreases by ≈ 20 % between n = 0.8 - 0.9. At higher temperature, T/t = 4, the compressibility decreases by ≈ 2 %. We resolve this effect as a shift in the average density between the initial density modulation pattern and the long-time equilibrium density. This effect is largest at the coldest temperatures, and is at most δ n_↑≈ 0.03, which is comparable to the uncertainty in our density.



§ TEMPERATURE FITTING








We determine the temperature of our clouds by fitting to determinantal quantum Monte Carlo (DQMC) results. Spin resolved density n_↑, n_↓, correlation functions n_↑(𝐫)n_↑(0)_c, and singles density n^s data is generated at U/t = 8 on a grid of chemical potentials and temperatures.  The total density and correlation functions are then interpolated on a regular grid of density and temperature points. 

We use these interpolating functions to simultaneously fit the singles density and the single-spin component correlations versus the total density. The only free parameter is the temperature. We apply an imaging fidelity correction of f=0.97 based on our measured hopping and loss rates during imaging. An example fit is shown in Fig. <ref>.



§ COMPRESSIBILITY




We measure the compressibility of our gas in a harmonic trap with no additional potential provided by the DMD. For low temperatures, the lattice beams provide all the radial confinement, leading to ω̅= (2π) 185(10)  Hz. For hotter temperatures, we use a circular beam to provide extra confinement, leading to ω̅= (2π) 280(10)  Hz.  

We determine the harmonic trapping frequencies by fitting the density and nearest-neighbor density correlation profiles of a weakly interacting gas obtained at a field of 568.0(1)G, near the noninteracting point of the |1⟩-|3⟩ mixture, to the expected values for a non-interacting Fermi gas. These are determined from

    n_↑   =   1/N∑_𝐤 f(ϵ_k - μ, T)
    1/4 C^↑(𝐝)    =    - 1/N^2|∑_𝐤 f(ϵ_k - μ,T)e^-i𝐤·𝐝|^2,

where f(ϵ,T) is the Fermi function, and 𝐤 runs over the N allowed lattice momenta. Distance and energy scales are measured in units of the lattice constant and hopping respectively. We assume a harmonic trapping potential, μ(r) = μ_o - 1/2 m ω^2 r^2 and fit our cloud profiles with μ_o, ω, and T as free parameters.

After determining the trapping frequency, we compute the compressibility according to,

    . (∂ n/∂μ) |_T = - 1/mω^2(1/r∂ n/∂ r).


Our cloud is slightly elliptic, with an aspect ratio of ω_x/ω_y ≈ 1.2. Prior to determining the trapping frequency we perform an azimuthal average, which effectively rescales our coordinates (x,y) →( x √(ω_x/ω_y), y √(ω_y/ω_x)). We measure r in these coordinates above, therefore our fitting procedures yields ω̅= √(ω_x ω_y).

We compare our measured compressibility with DQMC, FTLM, and DMFT in Fig. <ref>A. The DQMC and FTLM compressibilities agree well with the experimental data and do not saturate at low temperatures. In contrast, the single-site DMFT compressibility saturates at T/t ≈ 1. The increasing compressibility below this temperature may be associated with short-range correlations <cit.>, which are not accounted for by single-site DMFT. Cellular DMFT results using a 2 × 2 plaquette gives excellent agreement with DQMC, supporting this interpretation.

At high temperatures, the compressibility is expected to scale as 1/T with T χ_c(T) = n(1-n/2), for finite U <cit.>. We plot the compressibility times the temperature, T χ_c(T) in Fig. <ref>B. T χ_c has not yet saturated in the temperature range considered here. We expect that saturation occurs at temperatures much hotter than the bandwidth, T ≫ 8t.



§ THERMOELECTRIC EFFECTS


We prepare our sample in thermal equilibrium, therefore there  are initially no thermal gradients. However, thermal gradients may be generated during the subsequent dynamics. To check this possibility, we looked at a wave vector with underdamped oscillations at our lowest temperature. We measured nearest-neighbor correlations, C^↑(0,1), at the time where the amplitude first crosses zero. Here the density is flat, and any spatial variations in the correlator must be due to thermal gradients. We did not find any evidence for generation of thermal gradients.

Thermoelectric coupling is primarily due to two effects. The first is thermodynamic, and is described by the thermoelectric susceptibility, ζ = -∂^2 Ω/∂μ∂ T = .∂ n/∂ T|_μ = .∂ S/∂μ|_T, where Ω = ϵ - ST - n μ is the grand potential and S is the entropy. This is a static quantity, and can be computed, e.g., by FTLM. In the whole temperature regime accessible by FTLM we find |ζ| ≲ 0.015 t^-1. This is small in the sense that generating a density gradient of 0.01 a^-1 requires a large temperature gradient of ∇ T ≈ 0.8 ta^-1.

The Seebeck coefficient is more difficult to calculate. Using the Mott-Heikes approximation <cit.>   or the Kelvin formula <cit.>   gives a small Seebeck coefficient due to a maximum of entropy which occurs close to n≈ 0.83. This is in agreement with previous observations using different models or techniques <cit.>. A detailed description of particle diffusion in the presence of thermoelectric effects can be found in Ref. <cit.>. In the bad-metallic or high-T regime where the Kelvin formula is a good approximation for the Seebeck coefficient <cit.>, the effect of thermoelectric coupling on particle diffusion is negligible.



§ COMPARISON OF THEORY TECHNIQUES

In the main text, we compare experimental results with three theory techniques, DQMC, FTLM, and DMFT. Each of these have different strengths and weaknesses which makes one or another more suitable for certain comparisons. We provide a broad outline their strengths and weaknesses in this section, and more detailed information on each technique in the following sections.

DQMC is the method of choice for calculating static quantities because it is an exact technique which can access the lowest temperatures we reach in the experiment. On the other hand, dynamical quantities are difficult to extract, as DQMC yields imaginary time Green's functions which must be analytically continued to real time. Analytic continuation is an ill-posed problem, and the statistical uncertainty of DQMC data further complicates matters. In the main text, we use DQMC to compute only static quantities, including the singles density, spin-up density correlator, and compressibility.

FTLM is an exact technique which provides direct access to both static and dynamic correlators but is not capable of reaching temperatures as low as DQMC. The minimum temperature it can access is limited by finite size effects. When correlation lengths exceed the size of the small system used, the results no longer reflect the behavior of the system in the thermodynamic limit. In the main text we use FTLM to compute all dynamic quantities which can be obtained from correlators (momentum relaxation rate and conductivity). We also provide the FTLM compressibility in the supplement to verify that this agrees with the DQMC result.

DMFT is an approximate technique which maps an interacting problem onto a self-consistent quantum impurity problem. It becomes exact in infinite dimensions. The impurity problem can be solved using a variety of techniques. When exact diagonalization is used to solve the impurity problem, dynamical quantities can be obtained directly, without analytic continuation. DMFT calculations can be performed at lower temperatures than attainable by DQMC or DMFT. In the main text we use DMFT to compute all dynamical quantities which can be obtained from correlators. We also provide the DMFT compressibility in the supplement to verify that this agrees with the DQMC result. This comparison provides one test of the DMFT approximations.



§ DETERMINANTAL QUANTUM MONTE CARLO (DQMC)


We perform DQMC calculations using the Quantum Electron Simulation Toolbox (QUEST) <cit.>   on an 8 × 8 homogeneous square lattice. The inverse temperature is split into L=40 imaginary time slices, where L Δτ = β. We perform 5000 warm up sweeps, 50000 measurement sweeps and between 100 and 1000 passes to accumulate adequate statistics. We find that the sign problem at n = 0.83 and U/t = 7.5 becomes important below T/t = 0.5. Reliable results in the range T/t = 0.3-0.5 can be obtained with additional statistics. Below T/t = 0.27, the sign approaches zero.



§ FINITE-TEMPERATURE LANCZOS METHOD (FTLM)


The finite-temperature Lanczos method (FTLM) <cit.>   is an exact diagonalization approach on small clusters (4× 4 in our case). The method employs Lanczos diagonalization, which yields exact extremal eigenstates and effective or approximate eigenstates in the middle of the many-body spectrum.  These states are together with the aid of the sampling over random initial vector used to calculate finite-T properties. Results have unwelcome finite size effects, which are large below some temperature T_fs, but are small at T>T_fs due to shorter correlation lengths at high T. To reduce the finite size effects we employ averaging over twisted boundary condition <cit.>   (with N_Θ=64 different boundary conditions) and summation over all symmetry sectors, e.g., we use the grand canonical ensemble.  We do not show FTLM results for T<T_fs.

Dynamical quantities like the optical conductivity σ(ω) are calculated as correlation functions via evaluation of matrix elements of, e.g., the current operator, between the effective many body eigenstates, which are in this case obtained from two separate Lanczos procedures <cit.>. The spectra, represented as a sum of weighted delta functions δ(ω-ω_i), needs to be further smoothed or broadened by η. Due to the exponential number of many body states and a very dense spectra, particularly at elevated T, the broadening can be relatively small and in our case η∼0.1t ≪Γ. FTLM gives very good results for the static quantities with negligible finite size dependence and no additional broadening needed. Dynamical quantities are more challenging and in particular the optical conductivity at finite frequencies is quite robust while at the lowest frequencies it still shows a small dependence on cluster size or shape and on the broadening used. We estimate the overall uncertainty for the most challenging dc conductivity to be at most 15 percent.

The optical conductivity (Fig. <ref>) for our parameters exhibits a Drude peak at low ω and a separated Hubbard band at ω∼ U. We fit the low-ω Drude peak to a Lorentzian, which is the form predicted by our hydrodynamic model (eq. <ref>), and extract the momentum relaxation rate Γ presented in the main text.





§ DYNAMICAL MEAN-FIELD THEORY (DMFT)


We describe the dynamical mean-field theory self-consistency loop and then the various impurity solvers used in this paper.



 §.§ Self-consistency
 Dynamical mean-field theory is exact in infinite dimension <cit.>. In finite dimension, it approximates the interacting problem by solving a self-consistent quantum impurity problem. A quantum impurity problem is the problem of a single site connected to an infinite bath of non-interacting electrons. The self-consistency is achieved by taking the same self-energy for both the quantum impurity problem and the lattice Green's function, and then requiring that the lattice Green's function projected on a single site equals the impurity Green's function. This takes into account both the localized physics of the atom and the itinerant character of the metal, in competition. When spatial correlations become important in low dimension, a cluster replaces the single site. The latter is known as Cellular Dynamical Mean-Field Theory (CDMFT) <cit.>.

We describe the CDMFT procedure. Single site DMFT is a special case where the cluster is replaced by a single site. For the Hubbard Hamiltonian, defined in eq. <ref>, we write formally an effective action containing an hybridization function Δ̂(τ-τ') that describes the degrees of freedom outside the cluster (the bath) as a time-dependent hopping within the cluster (which is easily pictured from a Feynman path integral point of view),

    S_eff=∫_0^βdτ dτ ^'Ψ _d^†(τ )[ (∂ /∂τ +μ -t̂)δ(τ-τ') - Δ̂(τ-τ') ]Ψ _d(τ ^') 
    
    + U∑_μ∫_0^βdτ n_μ↑n_μ↓.

Hats denote matrices in the cluster degrees of freedom labeled by the Greek letters μ ,ν. Here, t̂ is the hopping of the original Hamiltonian within the cluster. For the case of a 2× 2 plaquette, the spinor is defined by

Ψ _d^†≡ (d_1↑^†,… ,d_4↑^†,d_1↓^†,… ,d_4↓^†).  Physically, this
action corresponds to a cluster embedded in a self-consistently determined medium.

Given the effective action with a starting guess for Δ̂(τ-τ'), the cluster propagator G_c is solved with three methods: Two variants of continuous-time quantum Monte Carlo, and exact diagonalization, on which we comment further below.  Once the cluster Green's function is obtained, we extract the cluster self
energy from Σ̂_c=𝒢̂_0^-1-G
_c^-1 where 𝒢_0^-1 is the quantity in square brackets in the quadratic part of the action, while G_μν ,σ≡ -⟨ Td_μσ(τ )d_νσ^†(0)⟩ is the imaginary-time-ordered Green's function. Using the self-consistency condition in Matsubara frequency,

    iω _n +μ -t̂ - Δ̂(iω _n)=[ N_c/(2π )^2
    ∫ d𝐤̃ G(𝐤,iω _n)]
    ^-1+Σ̂_c(iω _n)

with

    G(𝐤,iω _n)=[ iω _n+μ -t̂(
    𝐤)-Σ̂_c(iω _n)] ^-1,

we recompute the hybridization function Δ̂(iω _n) and iterate till
convergence. Here t̂(𝐤) is the Fourier transform of the
superlattice hopping matrix, N_c is the number of sites within the cluster and the integral over 𝐤 is performed over
the reduced Brillouin zone of the superlattice.



 §.§ Impurity solvers
 The continuous-time quantum Monte Carlo solvers sample observables with a Markov chain defined in the space of Feynman diagrams of all orders. In the continuous-time auxiliary field method (CT-AUX) <cit.>, the action is expanded in powers of the Hubbard interaction. This approach works better when U is less than the bandwidth. Expansion in powers of the hybridization function generates the so-called CT-HYB solver <cit.>, which works better at values of U larger than the bandwidth <cit.>. For CT-HYB, we use a program that contains several improvements for speed <cit.>. The results are obtained from an average over the last 20 converged iterations and typically between 10^8 and 3× 10^9 Monte Carlo updates. High frequency tails are usually of higher quality in the CT-AUX approach.

In the exact-diagonalization approach <cit.>   that we used for single-site DMFT, the impurity problem is represented by an Anderson-like Hamiltonian H_imp with a
discrete number of bath orbitals (here 5 for each spin component) coupled to the impurity 

    H_imp≡∑_mσϵ _mσa_mσ^†a_mσ +∑_m σ(V_m σa_mσ^†c_σ+h.c.)+Un_↑n_↓.

Here, m =1...5 for each spin component such that we have 10 bath energy
levels ϵ _mσ coupled to the impurity via the
bath-cluster hybridization matrix V_m σ. The hybridization function is obtained from

    Δ_σ(iω_n)=∑_mV_mσ^2/iω_n-ϵ_mσ.

The parameters ϵ_m σ and V_m σ 
are determined by imposing the self-consistency condition in Eq. <ref> using a conjugate gradient minimization algorithm with a distance
function

    d=1/N_max∑_n=0^N_max|( 𝒢̂
    _0^' -1(iω _n)-𝒢̂_0^-1(iω
    _n)) | ^2

where N_max is the largest Matsubara frequency index, determined by choosing a high-energy cutoff of about 2000 (energies are given in units of hopping t, and we take ħ =1 and k_B=1). The distance function in Eq.(<ref>)
is computed on the imaginary frequency axis since the hybridization function is a smooth function on that axis. We take a convergence criterion of 10^-5 for the distance. We checked that the compressibility agrees with the continuous-time solvers to three significant digits. 



 §.§ Conductivity calculation


Since in single-site DMFT the vertex corrections vanish, the optical conductivity is calculated from the single-particle spectral weight using

    σ'(ω)=χ^''_jj(ω)/ω    =    π∑_σ∫_-4t^4tdε∫dω' 𝒯(ε) 𝒜(ε, ω') 𝒜(ε, ω' + ω) 
          ×[ f(ω') - f(ω' + ω) ]/ω  ,

where f is the Fermi-Dirac distribution, χ^''_jj(ω) is the imaginary part of the current response function, and 𝒜(𝐤,ω) is the spectral function containing the non-interacting square-lattice dispersion ε_𝐤 and the impurity self-energy, and normalized so that ∫dω 𝒜(𝐤,ω) = 1. Here, the usual integral over wave-vectors has been replaced by an integral over the band energies ε weighted by the longitudinal transport function <cit.> 

    𝒯(ε) = ∑_𝐤( ∂ε_𝐤/∂ k_x)^2 δ(ε - ε_𝐤) = -1/2∫_-4t^ε z N_0(z)  dz
  
containing the non-interacting density of states N_0, normalized so that ∫ N_0(z) dz = 1.

The real part of the conductivity obeys the f-sum rule in the following form <cit.>

    ∫dω/πσ'(ω)=1/N∑_𝐤σ∂^2ϵ_𝐤/∂ k_x^2< n_𝐤σ> = -1/2E_kinetic,

where < n_𝐤σ> is the expectation value of the occupation number in state 𝐤σ and E_kinetic is the expectation value of the kinetic energy for this two-dimensional system with nearest-neighbor hopping only. This means that even in situations where the optical conductivity is dominated by a Drude peak whose width is temperature independent, as in the range 4<T<8 in the inset of Fig. 3, the DC conductivity can decrease with temperature because at high temperature the kinetic energy decreases as 1/T <cit.>
. 

The exact diagonalization method, which is used to obtain the momentum relaxation rate and conductivity presented in the main text, allows one to obtain results directly on the real axis. However, the discrete nature of the bath introduces some uncertainty because the discrete energy levels must be broadened as Lorentzians of width η. When there is a range of η where the results are independent of η, one can be confident of the results. The error bars on the DMFT results in the main text correspond to the difference between η=0.1 and a five times smaller value, η=0.02. The estimate of of the momentum relaxation rate Γ is much less sensitive to η than the DC conductivity. 

We checked that the resistivity obtained in single-site DMFT is independent of the solver used by comparing the exact diagonalization solver result with Padé analytic continuation of the CT-AUX solver. We found good agreement up to T/t=2. At higher temperature we cannot make this comparison because analytic continuation is unreliable. We also verified that Γ agrees for the two DMFT solvers.

