## Simulating transport with cold atoms

Much can be learned about the nature of a solid from how charge and spin propagate through it. Transport experiments can also be performed in quantum simulators such as cold atom systems, in which individual atoms can be imaged using quantum microscopes. Now, two groups have investigated transport in the so-called Fermi-Hubbard model using a two-dimensional optical lattice filled with one fermionic atom per site (see the Perspective by Brantut). Moving away from half-filling to enable charge transport, Brown *et al.* found that the resistivity had a linear temperature dependence, not unlike that seen in the strange metal phase of cuprate superconductors. In a complementary study on spin transport, Nichols *et al.* observed spin diffusion driven by superexchange coupling.

## Abstract

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 lithium-6 in a two-dimensional 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 by 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 developed in this study may be applied to measurements of other transport quantities, including the optical conductivity and thermopower.

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 (*T*)–dependent resistivity ρ(*T*) should scale like *T*^{2} at low temperature (*1*). 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 (*2*, *3*). 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 (*4*) and other anomalous scalings in d- and f-electron materials (*5*). Another is the violation of the resistivity bound ρ < ρ_{max}, which is observed in a wide variety of materials (*6*). 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 (*7*). Approaches introduced to understand these anomalous behaviors include hidden Fermi liquids (*8*), marginal Fermi liquids (*9*), proximity to quantum critical points (*10*) and associated holographic approaches (*11*), and many numerical studies of model systems, most notably the Hubbard model (*12*) and the *t* − *J* model (*13*).

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 because of the finite and isolated nature of these systems. Most fermionic charge transport experiments have focused on studying either mass flow through optically structured mesoscopic devices (*14*–*17*) or bulk transport in lattice systems (*18*–*22*). In this study, we explored 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 is the compressibility. This requires only the assumption of linear response and the absence of thermoelectric coupling (*23*) and does not rest on assumptions concerning quasiparticles.

We realized the two-dimensional Fermi-Hubbard model by using a degenerate spin-balanced mixture of two hyperfine ground states of ^{6}Li in an optical lattice (*24*). 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 by 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. 1, 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. 1, C and D), always keeping the optical lattice at fixed intensity. We measure the density of a single spin component, , by using techniques described in (*24*), giving us access to the total density through .

We work at average total density . This value is close to a conjectured quantum critical point in the Hubbard model (*25*). Our lattice depth is 6.9(2) *E*_{R}, where *E*_{R} is the lattice recoil energy and *E*_{R}/*h* = 14.66 kHz, 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 690 G. 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 (*26*).

We observe the decay of the initial sinusoidal density pattern over a period of a few tunneling times. The short time scale 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. 1, A and C). We fit the average modulation profile to a sinusoid, where the phase and frequency are fixed by the initial pattern (Fig. 2A). The time dependence of the amplitude of the sinusoid quantifies the decay of the density modulation (Fig. 2B). Our experimental technique is analogous to that of Hild *et al*. (*27*), who 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 (*23*). This approach leads to a differential equation for the density decaywhere Γ 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 . 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 by 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. 2B.

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 on the empirical basis that our simple model fits the data and that we have not been able to detect any measurable temperature modulation in the gas (*23*). In addition, theoretical work suggests that the thermopower (Seebeck coefficient) is negligible near our doping (*13*, *28*).

For the remainder of this study, we focus on the temperature dependence of *D* and Γ. 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, , and the nearest-neighbor correlations between spin-up atoms , where **i** = (*i _{x}*,

*i*) is the site index and

_{y}**d**is the separation between sites. We compare these quantities with determinantal quantum Monte Carlo (DQMC) simulations to extract the temperature (

*23*). For temperatures at the low end of the range we can access, 0.3 <

*T*/

*t*< 1, the density correlations are a sensitive thermometer. At higher temperatures, 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. 3A. At high temperatures, *D* is expected to saturate, eventually approaching an infinite temperature-limiting value (*29*). The diffusion constant is closely related to the mean free path, *l*, and is often estimated as , where is the mean quasiparticle velocity (*30*). Therefore, the MIR limit implies a lower bound on the diffusion constant, , where *a* is the lattice constant and *ħ* is Planck’s constant *h* divided by 2π. 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 λ→∞ (*23*), and exact techniques such as diagonalization of finite systems and DQMC are limited to small systems. Even determining the infinite temperature-limiting value is a nontrivial quantum dynamics problem (*31*, *32*).

In a clean system like ours, momentum relaxation can occur only because of 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. 3A, inset). This trend may again be understood as Pauli blocking suppressing momentum relaxation at low temperatures.

We compare the experimental Γ with 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 angular frequency ω ~ *U*, but this does not affect Γ substantially (*23*). 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 the modulation amplitude (*23*). We found that *D* is insensitive to the amplitude in the range explored. Γ shows some amplitude dependence, but because of the large error bars we cannot conclusively say whether this is the source of the discrepancy between experimental and FTLM results (fig. S1).

To extract a resistivity by 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 (*23*, *33*, *34*). The measured compressibility increases with decreasing temperature (Fig. 3B). For our highest experimental temperatures, χ_{c} approaches *n*(1 − *n*/2)/*T*, as expected in the high temperature limit (*29*). 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. 4 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 by using the MIR limit is (*6*, *30*). 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* = 8*t* 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 (*35*, *36*).

The failure of ρ to saturate at the resistivity bound is similar to behavior observed in bad metals at high temperatures (*6*). In our system, the violation of the resistivity bound is associated not 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 (*30*). 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} + *AT* + *BT*^{2}. We find that the temperature dependence is linear to good approximation, as we obtain , , and . Alternatively, a power law fit to the form ρ(*T*) = ρ_{o} + (*CT*)^{α} yields , , and α = 1.1(1). Similar fits show that the inverse diffusion constant 1/*D* scales with α = 0.6(1) and that 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 nontrivial way. This behavior should be contrasted with the high-temperature regime, , where *D* saturates to a limiting value and the resistivity inherits its temperature dependence from the compressibility, which scales as χ_{c} ∝ 1/*T* (*29*). 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. 4). 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 that become relevant as correlation lengths approach the cluster size. For the four-site by four-site cluster considered here, these effects limit FTLM resistivity calculations to .

Because our experiment explores low temperatures that are inaccessible to FTLM, we also compare with an approximate technique, single-site DMFT (*37*) (Fig. 4). We find that 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 substantial quadratic component. These regimes coincide with two different regimes observed in the DMFT compressibility (*23*). Previous DMFT studies at greater interaction strengths have also observed these two linear regimes at intermediate temperatures, finding evidence for resilient quasiparticles in the lower-temperature regime (*35*, *36*). 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 by using a similar approach. Both of these quantities might be expected to show anomalous scalings, as in the cuprates (*4*, *9*). In line with theoretical work such as (*35*, *36*), searching for direct signatures of resilient quasiparticles by using spectroscopic techniques (*38*) 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.

## Supplementary Materials

www.sciencemag.org/content/363/6425/379/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S4

This is an article distributed under the terms of the Science Journals Default License.

## References and Notes

**Acknowledgments:**We thank J. Thywissen, M. Zwierlein, and S. Hartnoll for stimulating discussions. We thank M. Charlebois and P. Sémon for contributions to the continuous-time Monte Carlo impurity solver codes. We thank D. Bergeron for assistance with the analytic continuation and comparison of results with the Two-Particle–Self-Consistent (TPSC) approach. 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.

**Funding:**This work was supported by the NSF (grant DMR-1607277), the David and Lucile Packard Foundation (grant 2016-65128), the AFOSR Young Investigator Research Program (grant 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. 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.

**Author contributions:**W.S.B. and P.S. conceived the study. P.T.B, E.G.-S., and D.M. collected and analyzed the experimental data. R.N., A.R., C.-D.H., S.B., and A.-M.S.T. performed the DMFT simulations. J.K. performed the FTLM simulations. D.A.H. developed the hydrodynamic model and interpreted the results. All coauthors assisted with writing and revising the manuscript.

**Competing interests:**The authors declare no competing financial interests.

**Data and materials availability:**Data reported in this paper and the code required to reproduce the data analysis are archived in Open Science Framework (

*39*). Code associated with the FTLM results from the paper, code associated with DMFT results, and code for the transport calculations are available from GitHub (

*40*

*–*

*42*).