Subdiffusion and heat transport in a tilted 2D Fermi-Hubbard system

Using quantum gas microscopy we study the late-time effective hydrodynamics of an isolated cold-atom Fermi-Hubbard system subject to an external linear potential (a"tilt"). The tilt is along one of the principal directions of the two-dimensional (2D) square lattice and couples mass transport to local heating through energy conservation. We study transport and thermalization in our system by observing the decay of prepared initial density waves as a function of wavelength $\lambda$ and tilt strength and find that the associated decay time $\tau$ crosses over as the tilt strength is increased from characteristically diffusive to subdiffusive with $\tau\propto\lambda^4$. In order to explain the underlying physics we develop a hydrodynamic model that exhibits this crossover. For strong tilts, the subdiffusive transport rate is set by a thermal diffusivity, which we are thus able to measure as a function of tilt in this regime. We further support our understanding by probing the local inverse temperature of the system at strong tilts, finding good agreement with our theoretical predictions. Finally, we discuss the relation of the strongly tilted limit of our system to recently studied 1D models which may exhibit nonergodic dynamics.

Introduction.-While non-interacting particles in a tilted lattice potential have been studied for almost a century [1][2][3][4], the dynamics of strongly tilted and isolated many-body systems with strong interactions have been relatively unexplored. Characterizing the late-time behavior of such closed quantum many-body systems away from equilibrium is a topic of fundamental interest. In a series of recent papers [5][6][7][8][9][10][11] it was shown how irreversible dissipative dynamics can emerge from the unitary evolution of closed quantum systems. Thus generically we expect the transport of conserved quantities in such systems to behave hydrodynamically at late times as long as the system does thermalize. On the experimental front, advances in quantum simulation with cold atoms and other platforms have allowed for unprecedented control of quantum many-body systems, and for the controlled study of their dynamics [12][13][14][15][16][17][18]. For example, in a recent study diffusive charge transport was observed in an isolated strongly-interacting 2D Fermi-Hubbard system [18]. Here we follow that work by observing the dynamics of the same cold-atom Fermi-Hubbard system subject to a strong external linear potential, or "tilt", and find a crossover to qualitatively different subdiffusive behavior at strong tilts.
The dynamics of a weakly tilted 2D Fermi-Hubbard model were studied in Ref. [19] using semiclassical methods. That work formulated an understanding of the dynamics in which regions with positive local temperature heat up and transport charge "up" the tilt, and regions with negative local temperature [20,21] transport charge "down" the tilt as the system approaches global equilibrium. Notably, if a tilted quantum system of particles in a lattice does approach thermal equilibrium within a band, that equilibrium is characterized by an infinite temper- * These authors contributed equally to this work. † Corresponding author. Email: wbakr@princeton.edu ature [19]. In contrast, recent theoretical works [22,23] explored the prospect of a transition to a localized phase in strongly tilted interacting 1D systems. While some evidence for this was found, it was suggested that this was the result of energetically-imposed local kinetic constraints that conserve the center of mass (COM)-a phenomenon later referred to as "Hilbert space fragmentation" [24,25]. This mechanism for nonergodicity at strong tilts depends on factors such as the range of interactions, the dimensionality of the system, and the direction of the tilt. In what follows, we explore a system which does not exhibit such nonergodicity. Thus this work is most directly related to Refs. [18,19], although initial motivation for this study was derived from Refs. [22,23], and investigating any nonergodic aspects of tilted systems is an interesting avenue for future work.
In this work we study the effect of an external tilt on the late-time emergent hydrodynamics of a 2D cold-atom system. This is done by varying the tilt strength and observing the relaxation of prepared initial density waves of various wavelengths λ. We observe a crossover from a diffusive regime at weak tilts, where the relaxation time τ scales like τ ∝ λ 2 , to a subdiffusive regime at stronger tilts, where τ ∝ λ 4 . We then construct a hydrodynamic model that exhibits the same crossover and discuss the underlying physics that leads to the subdiffusive transport. Using the hydrodynamic model we extract the tiltdependent thermal diffusivity of this system. We further verify our understanding of the underlying physics by measuring the local inverse temperature profile of the system, thus confirming a prediction of our theoretical model that this profile should correspond to local equilibrium and be displaced by a quarter wavelength relative to the density profile.
System.-Our system is well-described by the tilted Fermi-Hubbard HamiltonianĤ =Ĥ FH − FN fxCOM whereĤ FH is the conventional Fermi-Hubbard Hamiltonian on a square lattice, F is the tilt strength,N f is the (a) An off-centered beam generates a potential at the atoms that is approximately linear in x and independent of y. Blue-detuned light projected through a spatial light modulator is used to prepare the initial density waves of our experiments, with tunable wavelength in the direction of the tilt and hard walls a distance of 35 a latt apart in the perpendicular direction. The figure is a schematic intended to portray the experimental setup and is not to scale. (b) Density of a single spin component vs. time averaged over ∼ 10 images. The dotted square denotes the region of interest (ROI) in which our measurements were taken. (c) Evolution of the y-averaged density in the ROI of (b) as a function of x. The data corresponds to a system with interaction energy U/t h = 3.9(1), tilt strength F a latt /t h = 0.99 (3), and an initial density modulation of wavelength λ/a latt = 11.46(3). The density profile is shown at times 0 ms (0 ℏ/t h ), 0.5 ms (2.6 ℏ/t h ), and 15 ms (36 ℏ/t h ) from top to bottom.
total number of fermions, andx COM is the x component of the COM. The repulsive on-site interaction energy is denoted by U , and the single-particle hopping energy by t h . We emphasize that the system is tilted in only one of the lattice directions, which we denote with x. We realize our tilted 2D Fermi-Hubbard model by loading a balanced mixture of two hyperfine ground states of 6 Li into an optical lattice [26]. The tilt is generated by an off-centered 1064 nm Gaussian beam of waist ∼ 180 µm, as depicted in Fig. 1(a). The resulting potential is linear to within 10% across a region of length 40 a latt (30 µm), where a latt is the spacing of the optical lattice, and the strength of the potential gradient can be tuned from 0 to ∼ h × 5.5 kHz/a latt [27]. The beam is oriented such that the gradient is aligned with one of the two principal axes of the square lattice. A spatial light modulator (SLM) is used to project sinusoidal potentials of tunable wavelength along the direction of the gradient, and also remove any harmonic confinement from trapping potentials in the region of interest, similar to what was done in [18]. This allows us to prepare initial density modulations of tunable wavelength. We also add "hard walls" in the direction perpendicular to the gradient in order to contain the atoms in that direction and keep the average density constant over the experimental runtime (see Fig. 1(a)).
The atoms are adiabatically loaded into the lattice plus SLM potential at zero gradient (no tilt). The sinusoidal component of the SLM potential is chosen such that the resulting atom-density wave varies spatially with 0.0 ≲ ⟨n i ⟩ ≲ 1.2 (see Fig. 1 We also performed experiments with smaller-amplitude density waves and found no qualitative difference in our results [27]. Once the initial density wave is prepared we suddenly turn off the sinusoidal component of the potential created by the SLM, and turn on the tilt potential, thus initiating the dynamics. We focus on a square region of interest with a size of 35 × 35 lattice sites and measure only the single spin component ⟨n i,↑ ⟩ using fluorescence imaging [26] since in a spin-balanced system ⟨n i ⟩ = 2 ⟨n i,↑ ⟩. We performed all experiments at an optical lattice depth of 7.4(1)E R , where E R /h = 14.66 kHz is the recoil energy and h is Planck's constant. This leads to a hopping rate of t h /h = 820(10) Hz. We work at a magnetic field of 595.29(4) G nearby a Feshbach resonance centered on 690 G. This leads to a scattering length of 472.0(9) a 0 , where a 0 is a Bohr radius, which translates to an interaction energy of U/t h = 3.9(1) in the Fermi-Hubbard Hamiltonian. We tune the tilt strength F to values of up to F a latt /t h ≈ 6 which allows us to explore tilts well above the crossover from diffusive to subdiffusive dynamics.
It is of note that we do not reach tilt strengths so strong that it would be accurate to describe our system over the experimental runtime using an effective Hamiltonian which exactly conserves the COM. Therefore we emphasize that this work does not focus on the physics of fracton-like systems with a strictly conserved dipole moment, nor does it explore the possible nonergodic dynamics in such systems, although these topics are an interesting direction for future research [22][23][24][25][28][29][30]. However, our tilted system does show an emergent conservation of the COM in the long wavelength limit where the potential energy of the tilt dominates the conserved total energy.
Results.-Our experimental protocol consists of preparing initial density waves of various wavelengths in a potential with tilt F and imaging the system's density profile after it has evolved under its own unitary dynamics for some time t. We analyze our data by averaging all measurements from a certain wavelength, tilt, and time, and we also average the density in the direction perpendicular to the tilt. This yields the averaged density profile along the tilted direction as a function of time, as shown in Fig. 1(c). For each wavelength, tilt, and time we fit the density profile to a sinusoid, n(x, t) =n + A(t) cos (ϕ(t) + 2πx/λ), after adjusting for any small amount of atom loss, with the wavelength being fixed by the fit to the initial profile. We extract both the phase ϕ and amplitude A of the sinusoidal fit as a function of time, normalizing the amplitude by its initial value A(0). The main results of this paper are derived from tracking the decay of the amplitude A(t) with time.
Any change in the phase with time is a result of the distance the center of mass "falls down" the tilt as the system heats up in the first band of the lattice potential. More precisely, an initial state with energy density corresponding to a finite temperature in the non-tilted Fermi-Hubbard system will evolve down the gradient of the tilted potential. As this happens the tilt does work ∼ F ∆x COM per particle for a bulk shift of ∆x COM , and this work gets converted locally to kinetic and interaction energy in the system (the t h and U terms) [19]. Since the t h and U terms can only accommodate up to an energy of order ∼ t h + U per particle before reaching infinite temperature, the shift of the COM of the system cannot be more than ∼ (t h + U )/F . We observe phase changes during the dynamics that are consistent with this approximate bound. We corroborate that the atoms are not excited to higher bands using a technique described in [31].
We also use λ/a latt = 7.69(3) for F a latt /t h ≈ 6 as the decay time of the longest-wavelength modulation becomes very large for this tilt. Decay times that we observe vary increasingly with the tilt strength F , from 1−5 ℏ/t h at zero gradient up to 10 3 −10 4 ℏ/t h for F a latt /t h ≈ 6. At each value of the tilt strength we fit a power law of the form τ ∝ λ α to our measured decay times. Diffusive relaxation has a characteristic τ ∝ λ 2 dependence (α = 2), while values of α > 2 indicate slower subdiffusive dynamics. Fig. 2 shows the full analysis for two of the values of F . From the extracted exponents α we observe a crossover from diffusive relaxation at weak tilts, where α ≈ 2, to subdiffusive behavior with an exponent of α ≈ 4 at stronger tilts. This crossover is shown in Fig. 3  theoretical prediction of our hydrodynamic model.
Our observation of diffusive dynamics at weak tilts is consistent with the analysis of Ref. [19], and with the diffusive transport observed in previous experiments on the same system at F = 0 [18], albeit at lower temperatures. The crossover to subdiffusion with α ≈ 4 at strong tilts was, until now, previously unobserved, and its observation and explanation is the main result of this work. Below, and more completely in the Supplement, we construct a hydrodynamic model of our system to help explain these observations. We also further test our Diffusive to subdiffusive crossover. Extracted scaling exponent α for τ ∝ λ α from datasets at different tilts (orange circles). As the tilt is increased from F a latt /t h = 0 to F a latt /t h ≈ 6 the relaxation of initial density waves crosses over from characteristically diffusive (α = 2) to subdiffusive with α ≈ 4. The shaded curve is a prediction of our hydrodynamic model which is derived in detail in the Supplement [27].
understanding of the mechanism behind the subdiffusive transport by experimentally verifying our model's predictions for the local temperature profile.
Hydrodynamic model.-We denote the non-tilt energy density due to t h and U terms by e(x, t), and the number density of fermions by n(x, t). Our system is, on average, uniform along the y direction, so e and n are assumed to only depend on x and t. n is a conserved density and so is ϵ = e − F xn, the total energy density including the tilted potential.
For nonzero tilt, our system heats up to near infinite temperature within the lowest band, where the thermodynamic properties are readily calculated using the hightemperature expansion. There are then three unknown transport coefficients in the most general formulation of our model: diffusivities for each of the two conserved densities and a thermopower coefficient which might be significant for this system since the energy and atom transport are strongly coupled by the tilt. Our data does not have enough detail to allow us to estimate all three of these transport parameters. However, in the stronger-tilt regime where τ ∼ λ 4 , a tilt-dependent thermal diffusivity is the only transport coefficient that enters in the relaxation, and thus this one parameter can be determined from our measurements. We therefore present a less general version of our model for this strong-tilt regime here, leaving the more general model to the Supplement.
Let us first consider the infinite temperature equilibrium that our system thermalizes to at late times. This is a limit of zero inverse temperature (β → 0) and infinite chemical potential (µ → ∞), with a finite spatially uniform βµ; we call this equilibrium valueβµ. This uniform equilibrium has atom number densityn = 2eβ µ /(1+eβ µ ) per site. It is convenient when separating the energy into tilt and nontilt terms to choose the interaction term at each site to be U (n ↑ − (n/2))(n ↓ − (n/2)). This choice amounts to changing the total energy and potential V (x) by constants, so it does not change the physics. With this choice, the equilibrium nontilt energy density vanishes: The density profile at finite time has an additional sinusoidal component: n(x, t) =n + A 0 e −t/τ cos kx with k = 2π/λ (choosing the origin so there is no added phase in the argument of the cosine). In the strong tilt, small k regime we are considering now, this density profile is at local equilibrium with a time-dependent and spatially nonuniform inverse temperature β(x, t). We assume the system is also near global equilibrium, so we work to lowest order in A 0 and β. Near position x, if we have local equilibrium in the tilted potential V (x) = −F x in this high temperature limit, the density is given by n(x) = 2e β(µ+F x) /(1 + e β(µ+F x) ). So, in the long wavelength limit we are considering here, the density gradient is dn/dx = F n(1 − (n/2))β(x). Thus to leading order the temperature profile is given by −A 0 ke −t/τ sin kx = Fn(1 − (n/2))β(x, t). Using this result along with a high temperature expansion to write e as a function of β to leading order, we obtain the nontilt energy profile (1) at local equilibrium to lowest order in A 0 and k. Now that we have determined the profiles of n and e assuming local equilibrium, next we consider the dynamics and use energy and number conservation to determine the relaxation time τ . In the regime we are now considering, the rate-limiting bottleneck is the transport of nontilt energy (heat) through the system. This limits the rate at which tilt energy can be converted to heat and dissipated to the rest of the system, and thus the rate at which the whole system relaxes. The relaxation of the number density implies, via the continuity equation for atom number, an atom number current density of This current flows along the tilt direction, locally converting tilt energy to nontilt energy. In addition, there is a heat current j h (x, t) = −D th ∇e(x, t) flowing due to the temperature gradients, where D th (F ) is a tilt-dependent thermal diffusivity. Conservation of energy is theṅ showing the contribution of heat diffusion and the conversion of energy from tilt to nontilt due to the atom current j n . In the strong tilt regime we are considering, the two terms on the RHS of Eqn. 3 are each much larger in magnitude than the LHS: the motion of the atoms converts tilt energy to nontilt energy and this is dissipated by thermal transport, while the amplitude of the inhomogeneities decays slowly (D th k 2 τ ≫ 1). In this strong tilt regime, the decay rate is and the condition for the validity of this regime is We use Eqn. 4 to extract the thermal diffusivity D th as a function of tilt strength in the regime consistent with τ ∝ λ 4 and plot the result in Fig. 4. From the validity condition of Eqn. 5 we can also estimate the location of the crossover shown in Fig. 3. Plugging in the experimental values of U/t h = 4 andn = 0.6, and any value of k from the experimental range ka latt ∈ [2π/24, 2π/12], we get the condition that α ≈ 4 when F a latt /t h ≫ 1, which is consistent with the data shown in Fig. 3. A more complete model is detailed in the Supplement, and this model is used to derive the superimposed curve of Fig. 3 which agrees quantitatively with our experimental results. This more detailed model also gives the thermal diffusivity D th in terms of all of the transport coefficients, including the thermopower. We therefore conclude that our hydrodynamic model captures the essential physics leading to the main observation of this paper: the crossover from diffusive to subdiffusive relaxation with τ ∝ λ 4 as the tilt becomes strong.
The picture we have laid out in this section is one where, at strong tilts and long wavelengths, the system quickly achieves local equilibrium, locking the local inverse temperature to the density profile. As the density profile decays, local number density currents flow, and by conservation of energy this necessitates the flow of nontilt energy in the system. It is this flow of nontilt energy that we claim bottlenecks the relaxation in the large F regime, and thus D th sets the relaxation rate of the system. A prediction of this understanding is local equilibrium between β(x, t) and n(x, t). We verify this prediction by measuring the single component density and singlon occupancy profiles in our system and solving for the inverse temperature in the atomic limit, which is an effective method of thermometry at such high temperatures. In Fig. 5(a) and (b) we show both the density and local inverse temperature profiles, the decay of both of their amplitudes, and the phase difference between them in time (inset). From this we see that the β profile is at local equilibrium, locked at a quarter wavelength phase shift from the density profile, and both profiles decay together in time, as predicted by our understanding of the subdiffusive regime of this system. Summary and outlook.-We studied a new regime of thermalization in a square-lattice cold-atom Fermi-Hubbard system subject to an external linear potential. Our system was effectively closed and evolved under its own unitary dynamics starting from prepared initial density waves of various wavelengths λ. By observing how the amplitude of these initial density modulations evolved in time we found two qualitatively different hydrodynamic regimes and a crossover between them: At weak tilts the system relaxes diffusively, in accordance with previous theory [19] and experiments [18]. At strong tilts, we found a new regime where the system relaxes subdiffusively with a decay time τ that scales as τ ∝ λ 4 . We argued that this subdiffusive behavior is a result of having to "drain" the large reservoir of tilt energy via the bottleneck of heat transport en route to global equilibrium, and is captured effectively by a hydrodynamic description with the system remaining near local equilibrium. To test this understanding we measured the local temperature profile and do indeed find that the system remains near local equilibrium as it relaxes in this subdiffusive regime. In the Supplement, we also develop and present a more complete and detailed hydrodynamic model that quantitatively captures the crossover between the diffusive and subdiffusive regimes (Fig. 3).
In the strongly tilted regime we used our model to extract the tilt-strength-dependent thermal diffusivity that bottlenecks the relaxation of the system. One perspective on why this novel subdiffusive regime appears is that in the strong-tilt and long-wavelength limit the center-of- mass potential energy is the dominant part of the total energy, so energy conservation becomes an emergent almost-conservation of the center of mass. In contrast to recent theoretical studies of potential ergodicity breaking in tilted 1D systems [22,23], in this work we focused on the novel effects of a tilt on the approach to equilibrium in an isolated system that does indeed thermalize. This thermalization was robust because our system had a tilt potential along only one of the two principal axes of the lattice, and the resulting unconstrained motion of atoms in the perpendicular direction produced good thermal baths in each such row of the lattice. To arrest this thermalization more microscop-ically, one avenue of future exploration will be to apply tilt potentials along both axes of the lattice to suppress such local thermalization.

SUPPLEMENTAL MATERIAL
The experimental setup and basic parameters are already described in detail in the supplement of Ref. [26]. The spatial light modulator calibration is also explained in the supplement of Ref [18].

Tilt potential calibration
To calibrate the gradient and characterize its homogeneity across the region of interest, we used the SLM to prepare an initial state consisting of three thin stripes of width ∼ 1 a latt and a separation of ∼ 20 a latt , with their long direction oriented orthogonal to the tilt direction. Each stripe consists of a spin-polarized gas of the lowest hyperfine ground state of 6 Li.
For weak tilts, we are able to directly measure Bloch oscillations of these non-interacting particles. We do so by fitting a Gaussian profile to the density profile integrated along the direction perpendicular to the tilt which is used to quantify the "breathing" oscillation of the width of the stripes. This is similar to what was done in [32]. From the theory of Bloch oscillations, we expect the width of each stripe to oscillate with a maximal half-width of A = 4t h /F and a period of T = h/F a latt . Thus, by fitting a sinusoid to the evolution of the width of each stripe, we can extract the tilt strength at their respective positions. Fig. S1(a) shows an example of such oscillations.
For stronger tilts, directly measuring the Bloch oscillations becomes challenging due to their small amplitude. Instead we use a modulation technique analogous to what was done in [33]. We modulate the lattice potential at frequencies on the order of the tilt strength. This brings lattice sites that were decoupled due to the tilt into resonance which results in photon-assisted tunneling. We again measure the width of the thin stripes versus modulation frequency and observe a broadening of the stripes at resonance. Fig. S1(b) shows an example of such a measurement.
We corroborated that for the same potential strength at intermediate tilts, the gradient extracted using the two techniques agrees.  (7) kHz with a maximal difference of 7.5% between stripes.

Linear Response
Our hydrodynamic model assumes linearity in the amplitude of the initial inhomogeneities. In this experiment, we worked with relatively large amplitude density modulations. In a previous study ( [18]), we worked with very small amplitude modulations and fit to a linear hydrodynamic model we developed. In the "tilted" system studied in this work, we are no longer working close to a ground state, and as such, the strength of the modulation is not expected to be as important. Fig. S2 shows a comparison between the decay of strong and weak density modulations in a tilted potential. We observe that when we normalize the sinusoid amplitude and look at its decay, there is no measurable difference between the decays within the errorbars. This justifies working with strong modulations in this work to reduce the statistical error in the measurements for a fixed number of repetitions.

Complete hydrodynamic model
We write our hydrodynamic theory in terms of the particle number density n(x, t) and nontilt energy density e(x, t), as well as their corresponding currents j n (x, t) and j e (x, t). Total particle number and total energy are conserved, and these conservation laws can be written asṅ + ∇ · j n = 0 (S1) e + ∇ · j e − F j n = 0. (S2) The entropic "force" laws that describe how currents are driven in this system are of the form j e = M e χ e + M ne χ n and j n = M n χ n + M en χ e , where χ e and χ n are the entropic forces determined by the profiles of e and n, and we insist on writing the forces in a "canonical basis" for which Onsager's reciprocal relations take the simple form M en = M ne . The M coefficients are dynamical coefficients that are, in general, difficult to determine from the microscopic model. The off-diagonal coefficient M ne is associated with thermopower-type effects in our system, thus this model is quite general aside from its assumption of linearity, which is well-supported by our experimental measurements. In this canonical basis, the force χ e is determined by the local change of entropy when an infinitesimal current of nontilt energy flows but no particle current flows. The force χ n is determined in a similar fashion, with an infinitesimal particle current and no nontilt energy current, but note that if an infinitesimal particle current flows, then due to energy conservation there must be a production (or depletion) of nontilt energy. Thus the forces in this system take the form where s is the entropy density of the Fermi-Hubbard model, and we have expanded this entropy density near infinitetemperature equilibrium with n =n and e =ē(n). The coefficients s ee , s nn , and s ne come from this high-temperature expansion: and we emphasize that these coefficients are known functions of U , t h andn for the Fermi-Hubbard model. Now that we have specified our model, we proceed in determining its eigenmodes and respective relaxation rates, with a particular focus on the slowest mode, which is representative of the late-time behavior that we analyze in the main text. To do this, we first organize our model into a matrix eigenvalue problem. The eigenmodes are functions of definite wavelength λ, and they decay to equilibrium at a rate τ −1 , i.e. e(x, t) −ē(n) = e −Γt (a cos kx + b sin kx) and n(x, t) −n = e −Γt (c cos kx + d sin kx), where k = 2π/λ and Γ = 1/τ . We therefore write the above deviations of e and n from global equilibrium as ρ = (a b c d) T in the basis {e −Γt cos kx, e −Γt sin kx, e −Γt cos kx, e −Γt sin kx}. In this language the currents are driven according to j = M U ρ and the conservation laws are written as −Γρ = −∇j, where Thus our model is solved via the eigenvalue problem Γρ Γ =∇M U ρ Γ . There are two solutions for Γ, each with a multiplicity of two corresponding to pure cos and sin waves for n(x, t) −n. The only solution we need to consider at late times is the slow mode with Γ = Γ − and n(x, t) −n ∝ cos kx. This eigenmode is representative of the dynamics of all monochromatic initial conditions at late times. In what follows we discuss some important features of the slowest eigenmode.
In the limit of small F (and/or large k) the slowest mode is diffusive, i.e. Γ ∝ k 2 . In the limit of large F (and/or small k) the slowest decay rate is where D th = M e − (M 2 ne /M n ) s ee is the thermal diffusivity, and we will discuss why we identify it as such below. Thus we see that our model crosses over from diffusive to subdiffusive with τ ∝ λ 4 as 1/F λ becomes small.
If we assume a scaling of the form Γ − ∝ k α we can estimate the exponent α by α = d log Γ− d log k k=ke evaluated at some k in the experimental range k ∈ [2π/24, 2π/12] denoted k e . The general expression for α evaluated this way depends on the dynamical coefficients M , but in the limit where Me Mn , Mne Mn ≪ snn see this dependence drops out and we get a parameter-free estimate of α as a function of F . In this limit α(F ) = 2 + 2 and this is the theoretical estimate of α(F ) that we use to compare to experimental results in the main text (Fig. 3). Now we examine the structure of the slowest eigenmode itself and explain why we identify D th as mentioned above. At small k/F , to leading order, the slowest eigenmode has (S11) We see that in this mode a modulation of number density with amplitude O(1) comes with a slow subdiffusive number density current j n ∝ k 3 that is "out of phase" by a quarter wavelength. This number current converts tilt energy to nontilt energy and this generates a small out of phase, nontilt energy profile with amplitude ∝ k/F . That nontilt energy diffuses and we see that the ratio of amplitudes of the resulting "in phase" (with n(x, t)) energy current to the energy profile it is depleting is |j e |/|e| = D th k with D th = M e − (M 2 ne /M n ) s ee as mentioned earlier. This is why we identify D th as such. The process of diffusing the nontilt energy that is generated by the particle current that is relaxing the density profile is the bottleneck process and obeys a diffusion equation with diffusivity D th . That is why this diffusivity shows up as the one unknown coefficient in Γ − , and thus we use our data to determine it in the regime where τ ∝ λ 4 where this mechanism is valid. Now let's address the β profile in this mode. The "in phase" component of e −ē(n) shown in Eqn. S11, which is larger than the out of phase component by a factor of F/k, is due to the difference betweenē(n) andē(n), and not due to a nonzero β component that is in phase. Since β is proportional to e −ē(n) at high temperatures, to leading order in the high-temperature limit β(x, t) is set by the out of phase component of e −ē(n) which is the same as e −ē(n) for that component because the out of phase component of n −n is zero by definition (since "in phase" and "out of phase" are defined relative to the n(x, t) profile here). Thus this model predicts an out of phase local β modulation with amplitude amp(β(x, t)) = −s ee   s nn s ee − s 2 ne s 2 ee   k F (S12) where we have used the high temperature expressions for s ee , s ne , and s nn . Indeed in the main text we show measurements of the local β that are consistent with this prediction (Fig. 5). We call the nontilt energy current that results from this local β profile the "heat current" j h . Thus D th is the diffusivity corresponding to the heat current that is being driven by the nontilt energy that is generated by the relaxation of the particle number distribution.

High temperature expansion
We compute the grand partition function of the Fermi-Hubbard model in the high temperature expansion to second order in β and evaluate the second partial derivatives of the entropy density with respect to n, the particle number density, and e, the energy density due to t h and U terms, in order to compute the coefficients s ee , s nn , and s ne . The results are s ee = 16 n(2 −n)(32t 2 h +n(2 −n)U 2 ) (S14) s nn = 64t 2 h + 2n(2 +n)U 2 n(2 −n)(32t 2 h +n(2 −n)U 2 ) (S15) s ne = 8nŪ n(2 −n)(32t 2 h +n(2 −n)U 2 ) . (S16) For the experimental parametersn = 0.6 and U/t h = 3.9 these coefficients take the values s nn ≈ 2.96, s ee ≈ 0.43, s ne ≈ 0.50 in units where t h = 1.

Simultaneous fitting of model
As explained in the previous sections, there is a fast and a slow exponential decay solution to our hydrodynamic model. In the strong tilt regime, Eqn. S9 shows that the slow decay depends only on the thermal diffusivity D th .
We perform a simultaneous fit to all wavelengths at a given tilt strength as explained in the supplement of [18]. The fitting function is and it is fitted only to the late-time decay. Here, A 0 is a fitting parameter that can vary for each wavelength while D th is fitted globally to all wavelengths. The parameters F and k are fixed according to our experimentally measured values. The results of fitting this model to measurements in the strong tilt regime are shown in Fig. S3.  19.33(7) a latt (purple), and 23.3(2) a latt (pink) at different tilts. The lines are simultaneous fits of the hydrodynamic model to the long-time decay after the initial average heating (phase change). We are able to extract the thermal diffusivity through this fitting method.