# Corrections to diffusion in interacting quantum systems

Alexios A. Michailidis  Department of Theoretical Physics, Université de
Genève, 1211 Genève, Switzerland Dmitry A. Abanin  Department of Theoretical
Physics, Université de Genève, 1211 Genève, Switzerland Department of Physics,
Princeton University, Princeton, New Jersey 08544, USA  Luca V. Delacrétaz
Department of Theoretical Physics, Université de Genève, 1211 Genève,
Switzerland James Franck Institute, University of Chicago, Chicago, Illinois
60637, USA  Kadanoff Center for Theoretical Physics, University of Chicago,
Chicago, Illinois 60637, USA

###### Abstract

Transport and the approach to equilibrium in interacting classical and quantum
systems is a challenging problem of both theoretical and experimental
interest. One useful organizing principle characterizing equilibration is the
dissipative universality class, the most prevalent one being diffusion. In
this paper, we use the effective field theory (EFT) of diffusion to
systematically obtain universal power-law corrections to diffusion. We then
employ large-scale simulations of classical and quantum systems to explore
their validity. In particular, we find universal scaling functions for the
corrections to the dynamical structure factor $\langle n(x,t)n\rangle$, in the
presence of a single $U(1)$ or $SU(2)$ charge in systems with and without
particle-hole symmetry, and present the framework to generalize the
calculation to multiple charges. Classical simulations show remarkable
agreement with EFT predictions for subleading corrections, pushing precision
tests of effective theories for thermalizing systems to an unprecedented
level. Moving to quantum systems, we perform large-scale tensor-network
simulations in unitary and noisy 1d Floquet systems with conserved
magnetization. We find a qualitative agreement with EFT which becomes
quantitative in the case of noisy systems. Additionally, we show how the
knowledge of EFT corrections allows for fitting methods, which can improve the
estimation of transport parameters at the intermediate times accessible by
simulations and experiments. Finally, we explore non-linear response in
quantum systems and find that EFT provides an accurate prediction for its
behavior. Our results provide a basis for a better understanding of the non-
linear phenomena present in thermalizing systems.

##  I Introduction

One of the main pursuits of condensed matter physics is the understanding of
out-of-equilibrium phenomena in many-body systems. Transport probing the only
slightly out-of-equilibrium, linear response, regime is particularly
accessible experimentally and therefore most important to understand
theoretically, to tie experiments with insight into the fundamental structure
of correlated matter. The experimental accessibility of linear response
observables has allowed to establish some of the most puzzling phenomenology
in condensed matter physics, including the $T$-linear resistivity Gurvitch and
Fiory (1987) of high-$T_{c}$ superconductors and heavy fermion systems,
anomalous Hall angles Chien  _et al._ (1991) and magnetoresistance Hayes  _et
al._ (2016), which have largely eluded explanations despite decades of
activity. More recently, experiments in synthetic quantum matter, such as cold
atoms and superconducting quantum circuits, have offered new tools to explore
quantum transport (e.g., Schneider  _et al._ (2012); Schreiber  _et al._
(2015); Wei  _et al._ (2022); Morvan  _et al._ (2022); Mi  _et al._ (2023)),
further emphasizing the need for a better theoretical understanding of the
landscape of transport phenomena in many-body systems. The theoretical
challenge lies in finding controlled methods to study dynamics in strongly
correlated systems.

Hydrodynamics – broadly understood as the emergent dynamics of conserved
densities in thermalizing systems – offers particularly suitable tools in this
regard, providing a framework to parametrize and understand near-equilibrium
dynamics at late times. In the hydrodynamic limit , the dynamics typically
follow a universal behavior, characterized by a dissipative universality
class, the most prevalent being diffusion. While deriving the hydrodynamics of
diffusion from microscopics is typically a challenging task, experimental and
numerical evidence strongly suggests that it describes the leading order
linear response at late times of thermalizing classical and quantum many-body
systems, across diverse scales and platforms. Yet, the dissipative
universality classes provide information beyond the leading late-time behavior
of linear response: they also include nonlinear response, and universal
scaling corrections to observables. In particular, the leading late time
behavior of simple observables such as the dynamic structure factor $\langle
n(x,t)n\rangle$ can be found from classical hydrodynamic equations Kadanoff
and Martin (1963). However, the understanding of corrections to linear
response and more complicated observables requires a framework for
hydrodynamic fluctuations that systematically treats noise. Several proposals
for doing so exist, including generalizations of the Martin-Siggia-Rose
formalism Martin  _et al._ (1973) to allow for non-Gaussian noise, Fokker-
Planck equations for continuous fields (e.g., An  _et al._ (2021)),
macroscopic fluctuation theory Bertini  _et al._ (2015), and effective field
theories on Schwinger-Keldysh contours Crossley  _et al._ (2017). It is not
clear which of these effective theories – if any – describes thermalizing
many-body systems beyond the leading late time behavior. Furthermore, the
differences between classical and quantum systems in terms of hydrodynamic
fluctuations remain ambiguous, as well as the capacity of these effective
theories to discern them.

Beyond identifying the correct theory of fluctuations, understanding
corrections to observables in thermalizing systems has important experimental
and numerical implications. Starting with numerics, a systematic theory of
scaling corrections is critical for quantum simulations, which can typically
access intermediate times, during which the effects of corrections can be
significant. This can result in inaccurate determination of transport
parameters, or even in an incorrect value of the dynamical exponent $z$, as
illustrated in Figure 1. Additionally, diffusive dynamics are also present in
non-generic systems, such as certain integrable systems (De Nardis  _et al._ ,
2018; Nardis  _et al._ , 2019; Gopalakrishnan  _et al._ , 2018) and non-
interacting systems where diffusion is induced by noise. In these cases, even
if the leading late time behavior is the same, the scaling corrections are
sensitive to the number and type of conserved densities; they therefore offer
precision tests of thermalization, unambiguously distinguishing various
apparently diffusive systems. In experiments, the presence of these
corrections has interesting consequences for the understanding of
thermalization in correlated materials. Power-law corrections to late-time
observables come with time scales related to the local equilibration time –
the time scale at which regular hydrodynamics kicks in. This time scale is
parametrically large in weakly coupled or nearly integrable systems, but
seemingly cannot be made arbitrarily small at strong coupling; this has lead
to the expectation that the local equilibration time is universally bounded by
the ‘Planckian’ time, $\hbar/T$ Sachdev (2007); Zaanen (2004); Bruin  _et al._
(2013); Hartnoll and Mackenzie (2022). We will show that the leading power-law
corrections are in fact entirely fixed by derivatives of diffusivities with
respect to the equilibrium value of the transported density or associated
potential $D^{\prime}\equiv dD/dn$ (e.g. temperature for heat diffusion);
hydrodynamics, therefore, universally ties nonlinear response, scaling
corrections to linear response, and dependence of transport coefficients on
experimental tuning parameters. Since the latter are readily available in
experiments and numerics, this provides a time scale that must be exceeded to
access the asymptotic regime.

A .\begin{overpic}[width=346.89731pt,tics=10,trim=-40.15pt 0.0pt -40.15pt
0.0pt,clip]{nonlinearity_2} \put(10.0,12.0){$\delta n(q)$}
\put(66.0,13.0){$\delta n(q/2)$} \put(66.0,-5.0){$\delta n(q/2)$}
\end{overpic}

B . $\displaystyle\qquad\qquad\quad\langle n(t)n\rangle\ $ $\displaystyle=$
$\displaystyle+$ $\displaystyle+$ $\displaystyle\cdots\qquad\qquad\qquad$
$\displaystyle=$ $\displaystyle\\!\\!\\!\\!\frac{1}{\sqrt{t}}$
$\displaystyle+$ $\displaystyle\\!\\!\frac{\sqrt{\tau}}{t}$ $\displaystyle+$
$\displaystyle\cdots\qquad\qquad\qquad$

C .

Figure 1:  (A) Nonlinear fluctuations of conserved densities are present in
generic many-body systems, and (B) lead to universal corrections to
hydrodynamics at late times. (C) In the case of a single diffusive density,
the leading correction is positive, and can cause a diffusive system (with
autocorrelation function illustrated in green) to appear superdiffusive
(yellow, $z=3/2$ is shown above) at intermediate times. The EFT of diffusion
predicts the coefficient of this correction $\tau=\frac{\chi^{2}D^{\prime
2}}{16\pi D^{4}}$, together with a universal scaling function of
$x/\sqrt{Dt}$, see Eqs. (1) and (2).

In this paper, we use the effective field theory (EFT) of diffusion Crossley
_et al._ (2017) to systematically and quantitatively study the corrections to
observables in generic diffusive systems. We focus on the dynamic structure
factor $\langle n(x,t)n\rangle$ in 1d lattice systems ($>$1d generalizations
are shown in Appendix) where nonlinear corrections are particularly strong,
but also study nonlinear response which offers a complementary verification of
the theory’s predictions. The potential significance of scaling corrections to
correlation functions in quantum systems has been appreciated for some time
Mukerjee  _et al._ (2006); Lux  _et al._ (2014); however, even the corrections
to the considerably simpler autocorrelation function, $\langle
n(0,t)n\rangle$, were only obtained analytically and observed numerically
recently Chen-Lin  _et al._ (2019); Glorioso  _et al._ (2021).

We start by sharpening these results and generalizing the approach in several
directions, to ultimately construct a theory of scaling corrections in
thermalizing systems, providing a framework to make quantitative predictions
systematically in an expansion at late times. We first show that the
coefficient of the leading correction is entirely fixed in terms of transport
parameters of the system and their derivatives with respect to equilibrium
densities – this allows us to establish that the correction is non-negative in
the case of a single diffusing density, making generic diffusive systems
appear superdiffusive at intermediate times (Fig. 1). Next, generalizing to
the dynamic structure factor $\langle n(x,t)n\rangle$, we find the universal
scaling function of $x/\sqrt{Dt}$ accompanying this correction. We also
compute subleading corrections, which arise from higher-order (2-loop)
fluctuation effects, as well as higher derivative terms in the EFT. These are
particularly important in systems with particle-hole symmetry, where the
leading, 1-loop, correction is absent. These new corrections come with their
own universal scaling functions of $x/\sqrt{Dt}$, summarized in Table 1. We
also present the EFT framework required to study study corrections in the
presence of multiple diffusive charges and derive the corrections for the case
of chaotic spin chains with $SU(2)$ symmetry.

Leading order  |  no particle-hole | particle-hole
---|---|---
non-linear  | $\frac{1}{t^{d/2}}F_{1,0}$ | $\frac{1}{t^{d}}\left(F_{2,0}+\log(t)\tilde{F}_{2,0}\right)$
linear  | $\frac{1}{t}F_{0,1}$ | $\frac{1}{t}F_{0,1}$
Table 1: Leading order corrections to the dynamical structure factor of
chaotic diffusive systems from loop corrections (non-linear) or higher
derivative corrections (linear), see Eqs. (1,4). The leading order non-linear
correction (1-loop) vanishes in the presence of particle-hole symmetry, and
therefore the sub-leading correction (2-loop) dominates. These non-linear
corrections are the leading corrections in one dimension ($d=1$).

We then quantitatively test these predictions in numerics. We first consider
classical lattice gases where $D(n)$ is known analytically, so that the theory
prediction can be compared to simulations without requiring any fitting
parameter. We find remarkable agreement for the entire scaling function
accompanying the correction to diffusion, shown in Fig. 4, thereby providing a
test of theories of fluctuating hydrodynamics with unprecedented level of
precision. We next show that the EFT corrections are also present in the
dynamics of interacting quantum spin-chains. In this case, the classical
resources required to accurately capture the dynamics grow rapidly with the
simulation time, and therefore our simulations cannot always reach asymptotic
times. We demonstrate that incorporating the EFT corrections into the fitting
process leads to considerably more accurate transport parameters, such as
diffusivity.

Finally, we discuss nonlinear response. We show that the EFT universally ties
higher-point functions of densities Delacretaz and Mishra (2023) to scaling
corrections to linear response. These observables can therefore be used to
understand which timescale must be exceeded to enter the asymptotic, late time
or low frequency, regime. As controlled experimental probes of nonlinear
response improve Nicoletti and Cavalleri (2016); Chaudhuri  _et al._ (2022),
this offers a quantitative correspondence between these observables and
thermalization. We expect these nonlinear observables, as well as fluctuation
corrections to linear response, to be within reach of current experiments in
cold atoms as well Schweigler  _et al._ (2017); Wei  _et al._ (2022); Mi  _et
al._ (2023). Measuring higher-point functions in numerics can also help
establish unambiguously the dissipative universality class with limited
resources.

The paper is organized as follows. In Section II we present the leading
corrections to the full dynamical structure factor in 1d, for systems with one
abelian local charge. Next we present the corrections for systems which
additionally exhibit particle-hole symmetry, and therefore exhibit a vanishing
leading correction. These corrections originate from both linear and non-
linear fluctuations. However, non-linear fluctuations are logarithmically
stronger in 1d, and are expected to dominate at long times. In Section III we
formulate the EFT formalism and present the main steps towards the calculation
of the corrections. In Section III.A we derive the leading 1-loop corrections,
and in Section III.B we outline the basic steps for the 2-loop calculation
required to get the non-linear corrections in systems with particle-hole
symmetry. Then, in Section III.C we discuss the structure of linear
corrections. Section III.D extends our results to systems with multiple
densities. In particular, we present the result for a single non-abelian
($SU(2)$) charge. We conclude by numerically verifying the leading order
corrections for classical systems, Section III.E. In Section IV we study the
linear response regime of quantum systems (coherent and incoherent) with
magnetization conservation. For incoherent systems we quantitatively verify
the presence of EFT predictions. Coherent dynamics are more complex and they
display longer transient phenomena which persist on all simulated timescales.
Nevertheless, our results qualitatively agree with EFT. In Section V we
explore the nonlinear response through a simple three-point function, which
offers a complementary test for the validity of EFT. In the concluding Section
VI, we consolidate our findings, discuss the relevance of these results for
the field, and outline potential avenues for future research.

##  II Full scaling function–analytical predictions

Whenever diffusion or any other hydrodynamic behavior emerges in a many-body
system, it is inevitably accompanied by scaling corrections that may be
important at intermediate times. These arise from higher-derivative
corrections Landau and Lifshitz (2013) as well as fluctuation (or ‘loop’)
corrections Alder and Wainwright (1970) in the hydrodynamic description. While
these corrections have been appreciated in the context of quantum many-body
systems for some time (see, e.g., Mukerjee  _et al._ (2006)), they are often
ignored. Since accessing late times in quantum simulations is fairly
prohibitive, accounting for these corrections to scaling is crucial even to
correctly capture the dissipative universality class of a given system.

One central result of this work is that these scaling corrections come with
entire universal scaling functions, which can be obtained from the EFT
Crossley  _et al._ (2017). For example, the leading correction to diffusive
correlation functions in one dimension comes from a 1-loop correction and
takes the form

$\langle n(x,t)n\rangle=\frac{\chi}{\sqrt{4\pi
Dt}}\left[F_{0,0}(y)+\frac{1}{\sqrt{t}}F_{1,0}(y)+O\left(\frac{\log
t}{t}\right)\right],$ (1)

where $F_{0,0}$ and $F_{1,0}$ are scaling functions of the scaling variable
$y\equiv x/\sqrt{Dt}$. The leading scaling function $F_{0,0}(y)=e^{-y^{2}/4}$
solves the linearized diffusion equation, whereas the leading correction
$F_{1,0}(y)$ comes from a 1-loop contribution Chen-Lin  _et al._ (2019) which,
as we show below, takes the form

$\displaystyle F_{1,0}(y)$ $\displaystyle=\frac{\chi D^{\prime
2}}{D^{5/2}}\tilde{F}_{1,0}(y)$ (2) $\displaystyle\tilde{F}_{1,0}(y)$
$\displaystyle=\frac{4+y^{2}}{16\sqrt{\pi}}e^{-y^{2}/2}+\frac{y(y^{2}-10)}{32}e^{-y^{2}/4}{\rm
Erf}(y/2)$

We have separated $F_{1,0}$ into a universal dimensionless function
$\tilde{F}_{1,0}$, and a non-universal factor that involves the susceptibility
$\chi$ and the diffusivity $D$ (like the leading order correlator), but also
the derivative of the diffusivity with respect to the background value of the
diffusing density $D^{\prime}\equiv dD(n)/dn$. If this parameter is known,
e.g. by measuring the diffusivity at several densities, the entire functional
form of the $1/\sqrt{t}$ correction to diffusion is fixed. One interesting
feature of this correction is that for $x=0$ it is non-negative
$F_{1,0}(0)=\frac{\chi D^{\prime 2}}{D^{5/2}}\frac{1}{4\sqrt{\pi}}\geq 0$,
which implies that the autocorrelation function approaches its asymptotic
diffusive form from above at late times:

$\langle n(0,t)n\rangle=\frac{\chi}{\sqrt{4\pi Dt}}\left(1+\frac{\chi
D^{\prime 2}}{4\sqrt{\pi}D^{5/2}}\frac{1}{\sqrt{t}}+O\Bigl{(}\frac{\log
t}{t}\Bigr{)}\right)\,.$ (3)

Therefore, if a dynamic critical exponent $z$ is extracted by fitting the
autocorrelation function as $\langle n(0,t)n\rangle\sim 1/t^{1/z}$ at late
times, a diffusive system will always naively appear to be superdiffusive,
$z<2$. This is illustrated in Fig. 1.

Eq. (1) includes the first two terms in a general expansion in derivatives and
fluctuations, whose structure is shown in Eq. (50). The correction to
diffusion that arises from $\ell$-loop contributions at $n$th order in the
derivative expansion in the EFT scales as $1/t^{n+\ell d/2}$ in $d$ spatial
dimensions and is encoded in a scaling function $F_{\ell,n}(y)$, which is
universal up to one (or a few) non-universal Wilsonian coefficients, similar
to the functions in Eq. (2).

Given that diffusivities generically depend on density, the leading correction
(2) is typically present. However, $D^{\prime}$ may vanish at special values
of the density: this occurs for example if there is a particle-hole (or charge
conjugation) symmetry, which commonly arises in lattice systems at half-
filling. In this case, the leading correction to diffusion takes the form

$\begin{split}&\langle n(x,t)n\rangle=\frac{\chi}{\sqrt{4\pi Dt}}\\\ &\
\times\left[F_{0,0}(y)+\frac{1}{t}\left(F_{0,1}(y)+F_{2,0}(y)\right)+O\left(\frac{\log
t}{t^{2}}\right)\right],\end{split}$ (4)

and comes a higher-derivative contribution $F_{0,1}(y)$, and a two-loop
contribution $F_{2,0}(y)$. The former can be shown to take the form

$\begin{split}F_{0,1}(y)&=\left[c_{1}(y^{2}-2)+c_{2}y^{2}(y^{2}-6)\right]e^{-y^{2}/4}\,,\end{split}$
(5)

where $c_{1},\,c_{2}$ are non-universal transport parameters, while the latter
is obtained in this paper, and is given by:

$\begin{split}F_{2,0}(y)&=\frac{\chi^{2}D^{\prime\prime
2}}{12\sqrt{3\pi}D^{3}}\left[\tilde{F}_{2,0}(y)+\tilde{F}_{2,0}^{\prime}(y)\log
t\right]\\\
\tilde{F}_{2,0}(y)&=\int_{0}^{\infty}\frac{ds}{\pi}\cos(sy)s^{2}\\\
&\\!\\!\times\left[s^{2}e^{-s^{2}}\left(\log\frac{1}{s^{2}}+{\rm
Ei}\left(\tfrac{2s^{2}}{3}\right)\right)-\frac{3}{2}e^{-s^{2}/3}\right]\\\
\tilde{F}_{2,0}^{\prime}(y)&=\frac{y^{4}-12y^{2}+12}{32\sqrt{\pi}}e^{-y^{2}/4}\,,\\\
\end{split}$ (6)

with ${\rm Ei}(z)\equiv-\int_{-z}^{\infty}\frac{du}{u}e^{-u}$. In the first
line, we again separated the scaling function into a non-universal factor,
which now depends on $D^{\prime\prime}\equiv d^{2}D(n)/dn^{2}$, and a
universal scaling function. We have slightly abused notation and allowed for a
factor of $\log t$ in the scaling function instead of writing it explicitly in
(4). Notice that a shift in this logarithm $\log t\to\log(t/\tau)$ can be
absorbed by the higher derivative corrections $c_{1},\,c_{2}$ in (5). At
asymptotically late times, the diffusive autocorrelation function is again
approached from above:

$\langle n(0,t)n\rangle=\frac{\chi}{\sqrt{4\pi
Dt}}\left(1+\frac{\chi^{2}D^{\prime\prime 2}}{32\sqrt{3}\pi D^{3}}\frac{\log
t}{t}+O\Bigl{(}\frac{1}{t}\Bigr{)}\right)\,.$ (7)

We note, however, that this correction has only a $\log t$ enhancement
compared to the non-sign definite $1/t$ corrections from (5). For reader’s
convenience, the universal scaling functions found above are illustrated in
Fig. 7 in the Appendix.

##  III Scaling corrections from the EFT

The universal corrections to diffusion quoted in Eqs. (1) and (4) can be
obtained from the effective field theory (EFT) of diffusion Crossley  _et al._
(2017); Chen-Lin  _et al._ (2019). These arise from thermal fluctuations
(loops) of the hydrodyanmic densities and noise fields. Several qualitative
properties of these loops were already understood shortly after their
discovery in classical numerics Alder and Wainwright (1970), e.g. through
mode-coupling approximations or the Martin-Siggia-Rose approach Martin  _et
al._ (1973), which established first steps towards a general EFT for
fluctuating hydrodynamics. The modern EFT approach completes these
constructions by elevating them into a systematic expansion in derivatives and
fluctuations; we will therefore follow this approach here.

We are interested in studying transport in a system with at least one
conserved quantity, leading to a continuity equation (in the continuum limit)

$\dot{n}+\nabla\cdot j=0\,.$ (8)

The density $n$ could correspond to energy density, charge density,
magnetization density, etc. We will first focus on the situation where a
single density is conserved, and discuss generalizations to multiple densities
in Sec. III.4.

A generating functional for correlation functions of densities and currents in
the thermal state $\rho_{\beta}$ can be written

$Z[A^{1},A^{2}]={\rm Tr}\left(U[A^{1}]\rho_{\beta}U^{\dagger}[A^{2}]\right),$
(9)

where $A^{1},\,A^{2}$ are background gauge fields that couple to the conserved
current in the time-evolution operator

$U[A]=\mathcal{T}\exp\left\\{-i\int_{-\infty}^{\infty}\\!\\!dt\left(H-\int
d^{d}x\,j^{\mu}A_{\mu}(t,x)\right)\right\\}\,,$ (10)

where we have collectively denoted the charge and current density by
$j^{\mu}=(n,j^{i})$. Derivatives of $\log Z$ with respect to $A^{1},\,A^{2}$
can generate correlation functions of $j^{\mu}$ with various time orderings.
If the system thermalizes, one expects the partition function to have a
representation in terms of a local effective Lagrangian of the long-lived
hydrodynamic variables. It is local in space and time because there are no
other long-lived excitations in the thermal state—this is the assumption of
thermalization. In the approach of Ref. Crossley  _et al._ (2017), this
effective Lagrangian is a function of the fluctuating density $n$, and a
conjugate field $\phi_{a}$:

$Z[A^{1},A^{2}]\simeq\int DnD\phi_{a}\,e^{i\int dtd^{d}x\mathcal{L}}\,.$ (11)

What is gained in universality is lost in exactness: while it is not an exact
representation of the microscopic partition function (9), Eq. (11) provides a
systematic expansion for it when background fields $A$ have slow variation in
time (and space) compared to the local equilibration time of the system. We
further motivate this construction in the Supplementary Material, and focus
here on how it is used to obtain universal corrections to diffusion. To
leading order in derivatives, the effective Lagrangian is found to be

$\mathcal{L}=i\sigma(n)\left(\nabla\phi_{a}\right)^{2}-\phi_{a}\left(\dot{n}-\nabla(D(n)\nabla
n)\vphantom{n^{2}_{a}}\right)+\cdots$ (12)

Here, $\sigma(n)$ and $D(n)$ are functions of the density that are not fixed
by the EFT: they correspond to the conductivity and diffusivity of the system.
These functions also play an important role in macroscopic fluctuation theory
Bertini  _et al._ (2015). In the present approach, they are just the leading
terms in a general expansion in derivatives (see for example Jain and Kovtun
(2022) for a discussion of certain terms in the EFT that do not appear in
constitutive relations).

As in most EFTs, it is typically impossible to derive (11) from a microscopic
model of interest. One exception is in the context of strongly interacting
holographic quantum field theories, where progress has been made in deriving
at least the quadratic part of the EFT from microscopics de Boer  _et al._
(2019); Glorioso  _et al._ (2018); Ghosh  _et al._ (2021) (see Policastro  _et
al._ (2002); Herzog and Son (2003); Skenderis and van Rees (2009); Nickel and
Son (2011) for earlier work in this direction); similar derivations may be
possible for lattice systems with large local Hilbert space dimension (e.g.,
Khemani  _et al._ (2018); Chan  _et al._ (2018)), or noisy systems in the
limit of strong noise (e.g., Bauer  _et al._ (2017); Claeys  _et al._ (2022)).

When studying linear response or more general correlation functions, one
expands these functions around the background value of interest for the
density $n=\bar{n}+\delta n$

$D(n)\simeq D+\delta nD^{\prime}+\frac{1}{2}\delta
n^{2}D^{\prime\prime}+\cdots\,,$ (13)

where $D,D^{\prime},D^{\prime\prime},$ etc. on the right-hand side are
evaluated at the background density $\bar{n}$. These are Wilsonian
coefficients of the EFT: the are not fixed by the EFT (and in fact are not
universal), but the EFT instead predicts how they enter in any late time
observable. Since the same coefficients enter in a large number of
observables, the problem is highly overdetermined and the EFT has substantial
predictive power.

In the following subsections, we use the EFT (12) to compute 1-loop and 2-loop
corrections to the retarded Green’s function of the charge density

$G^{R}(\omega,q)=\frac{\sigma q^{2}}{-i\omega+Dq^{2}}+\delta
G^{R}(\omega,q)\,.$ (14)

$G^{R}$ is simply related to the Fourier transform of the dynamical structure
factor through fluctuation-dissipation relations, but has a simpler analytic
structure and is therefore more convenient to work with. In the EFT, it can be
obtained from the mixed correlator (see App. A.2)

$G^{R}(\omega,q)=i\sigma q^{2}\langle n\phi_{a}\rangle(\omega,q)\,.$ (15)

At tree-level, this can be evaluated using the propagators of the fields
obtained from the Gaussian Lagrangian (12):

$\begin{split}\langle
n\phi_{a}\rangle_{0}(\omega,q)&=\frac{1}{\omega+iDq^{2}}\,,\\\ \langle
nn\rangle_{0}(\omega,q)&=\frac{2\sigma
q^{2}}{\omega^{2}+D^{2}q^{4}}\,.\end{split}$ (16)

Using (15), one recovers the leading diffusive behavior in (14). The second
piece $\delta G^{R}(\omega,q)$ comes from loop and higher derivative
corrections, which are studied below.

###  III.1 1-loop

Loop corrections to (14) arise due to nonlinearities in the EFT. For example,
expanding $D(n)$ as in (13) leads to a cubic term

$\mathcal{L}^{(3)}=\frac{1}{2}D^{\prime}\nabla^{2}\phi_{a}n^{2}\,.$ (17)

This produces a cubic vertex which, working perturbatively in these
interactions, will lead to loop corrections to $G^{R}$. Note that the
perturbative expansion is always controlled, because nonlinearities are
irrelevant. Indeed, (16) implies that density fluctuations scale as $\delta
n(t,x)\sim q^{d/2}$; since the cubic nonliniearity is suppressed by an extra
power of $\delta n$, it gives small corrections at late times or long
distances, where $\omega\sim q^{2}\to 0$. This is in contrast to momentum
conserving systems in $d=1$, where nonlinearities are relevant and lead to a
breakdown of diffusion which is replaced by the KPZ universality class Forster
_et al._ (1977). That the perturbative expansion is controlled in the present
situation is a derivation of the EFT, rather than an assumption.

The cubic action also contains a term proportional to $\sigma^{\prime}$. While
this term leads to nonlinear response Delacretaz and Mishra (2023), we show in
the supplementary material that it does not contribute to the 1-loop corrected
two-point function; we therefore ignore it here.

Figure 2: Left: 1-loop correction to diffusion. Right: 2-loop correction to
diffusion at half-filling. The propagators $\langle nn\rangle_{0}(\omega,q)$
and $\langle n\phi_{a}\rangle_{0}(\omega,q)$ from Eq. (16) correspond to the
solid lines, and half-solid half-dashed lines respectively.

The cubic vertex (17) leads to a 1-loop correction to $\langle
n\phi_{a}\rangle$ shown in Fig. 2. Its evaluation, performed in Ref. Chen-Lin
_et al._ (2019), is streamlined here. It is convenient to amputate the
external legs and parametrize the correction as $\delta D(\omega,q)$, namely:

$\delta\langle n\phi_{a}\rangle=-iq^{2}\delta D(\omega,q)\bigl{(}\langle
n\phi_{a}\rangle_{0}\bigr{)}^{2}\,.$ (18)

The one-loop correction then takes the form

$\delta D(\omega,q)=-iD^{\prime 2}\int_{p^{\prime}}q^{\prime 2}\langle
n\phi_{a}\rangle(p^{\prime})\langle nn\rangle(p-p^{\prime})$ (19)

where we have used the short-hand notation $p\equiv\\{\omega,q\\}$, and
$\int_{p}\equiv\int\frac{d\omega d^{d}q}{(2\pi)^{d+1}}$. The loop integrals
can be readily evaluated in any dimension (see supplementary material), and
give

$\delta D(\omega,q)=\frac{\chi D^{\prime
2}}{D^{2}}(-i\omega)\alpha_{d}\left(q^{2}-\frac{2i\omega}{D}\right)\,,$ (20)

where $\chi\equiv\sigma/D$ is the static susceptibility, and with

$\alpha_{d}(z)=\frac{(-z)^{\frac{d}{2}-1}}{(16\pi)^{d/2}\Gamma(\frac{d}{2})}\cdot\begin{cases}i\pi&\hbox{if
$d$ odd}\,,\\\ \log\frac{1}{z}&\hbox{if $d$ even}\,.\end{cases}$ (21)

The general scaling $\delta D/D\sim q^{d}$ agrees with expectations: the cubic
interaction is suppressed by $\delta n\sim q^{d/2}$, and two cubic vertices
were necessary to produce a loop correction. The detailed loop calculation was
necessary to obtain the overall coefficient, as well as the entire dependence
on the dimensionless ratio $Dq^{2}/\omega$. Nevertheless, several aspects of
the result could have been anticipated on general grounds: (i) the fact that
the correction vanishes in the static limit $\lim_{\omega\to 0}\delta
D(\omega,q)=0$ is required by the analyticity of equilibrium thermal
correlators due to the finite thermal correlation length Jain  _et al._
(2021), (ii) the existence of a branch point at $\omega=-\frac{i}{2}Dq^{2}$
follows from a simple cutting argument Chen-Lin  _et al._ (2019).

We are most interested in the case $d=1$:

$\delta D(\omega,q)=\frac{\chi D^{\prime
2}}{D^{2}}(-i\omega)\frac{1}{4}\frac{1}{\sqrt{q^{2}-\frac{2i\omega}{D}}}\,.$
(22)

Fourier transforming this expression, or rather $\delta G^{R}(\omega,q)$, is
straightforward, but presented in the supplementary material for completeness;
this results in a correction to the correlation function shown in Eq. (2).

###  III.2 2-loop–half filling

When $D^{\prime}=0$, e.g. due to particle hole symmetry, there is no one-loop
correction to diffusion. The leading fluctuation corrections come instead from
a two-loop diagram, shown in Fig. 2 arising from the quartic vertex

$\mathcal{L}^{(4)}=\frac{1}{6}D^{\prime\prime}\nabla^{2}\phi_{a}n^{3}\,.$ (23)

Because this interaction scales as
$\mathcal{L}^{(4)}/\mathcal{L}^{(2)}\sim\delta n^{2}\sim q^{d}$, and that two
such vertices will be necessary to give a non-analytic correction to the two-
point function, the two-loop correction will scale as $\delta G^{R}/G^{R}\sim
q^{2d}\sim\omega^{d}$ (up to logarithms). In $d=1$, these corrections are as
large as higher derivative corrections to diffusion, studied in the next
section, which scale as $q^{2}$.

The two-loop correction to $\langle n\phi_{a}\rangle$, with external legs
amputated, is given by

$\displaystyle\delta D(\omega,q)=-\frac{i}{2}D^{\prime\prime 2}$ (24)
$\displaystyle\quad\times\int_{p^{\prime},p^{\prime\prime}}q^{\prime\prime
2}\langle n\phi_{a}\rangle(p^{\prime\prime})\langle
nn\rangle(p^{\prime}-p^{\prime\prime})\langle nn\rangle(p-p^{\prime})\,.$

The integrals are evaluated in the supplementary material. One finds

$\delta D(\omega,q)=\frac{1}{2}\frac{(\chi
D^{\prime\prime})^{2}}{D^{2}}(-i\omega)\beta_{d}\left(q^{2}-\frac{3i\omega}{D}\right)\,,$
(25)

with

$\beta_{d}(z)=\frac{(-z)^{d-1}\log\frac{1}{z}}{(12\sqrt{3}\pi)^{d}\Gamma(d)}\,.$
(26)

This result has the expected $q^{2d}$ scaling, vanishes in the static limit
$\lim_{\omega\to 0}\delta D(\omega,q)=0$, and features the expected 3-diffuson
branch point at $\omega=-\frac{i}{3}Dq^{2}$ Delacretaz (2020).

For $d=1$, this becomes

$\delta D(\omega,q)=\frac{1}{24\sqrt{3}\pi}\frac{(\chi
D^{\prime\prime})^{2}}{D^{2}}(-i\omega)\log\frac{1}{q^{2}-\frac{3i\omega}{D}}\,,$
(27)

leading to a correction to $G^{R}(\omega,q)$ whose Fourier transform is
computed in (84) and shown in (6).

###  III.3 Higher-derivative corrections

Higher derivative corrections are also captured by the EFT for diffusion (12).
These will either involve extra time derivatives $\partial_{t}$, or two extra
spatial derivatives $\nabla^{2}$ (by reflection symmetry), and therefore give
corrections to the leading behavior that are suppressed by $q^{2}$, or
equivalently $1/t$. One can write the most general such higher derivative
corrections to the EFT (see Crossley  _et al._ (2017)). However, since we are
interested in the two-point function, we can instead directly write the most
general corrections to $G^{R}$. The higher derivative corrections should be
treated perturbatively, as quadratic vertices; the final expression therefore
contains at most two powers of the diffusive propagator. The most general
$O(q^{2})$ correction to the retarded Green’s function is therefore

$\begin{split}G^{R}(\omega,q)&=\frac{\sigma
q^{2}}{-i\omega+Dq^{2}}\left[1+c_{1}(-i\omega)+c_{2}\frac{(-i\omega)^{2}}{-i\omega+Dq^{2}}\right]\\\
&+\tilde{c}_{1}q^{2}+\tilde{c}_{2}(-i\omega)+O(q^{4})\,.\\\\[-10.00002pt]
\end{split}$ (28)

The two coefficients $\tilde{c}_{1},\,\tilde{c}_{2}$ are contact terms and
will not affect the correlation function at separated points: $\tilde{c}_{1}$
has the interpretation of a $q^{2}$ correction to the static susceptibility
$\chi(q)=\chi+\tilde{c}_{1}q^{2}+\cdots$, and $\tilde{c}_{2}$ is in fact
forced to vanish to guarantee $G^{R}(\omega,q\to 0)=0$. Fourier transforming
$\langle
nn\rangle=\frac{2}{1-e^{-\beta\omega}}\operatorname{Im}G^{R}(\omega,q)$ leads
to (5) (we have taken the liberty to redefine the non-universal coefficients
$c_{1},c_{2}$).

Eq. (28) can also be derived through more conventional approaches to
hydrodynamics Kadanoff and Martin (1963); Chaikin  _et al._ (1995): one writes
the linearized constitutive relation for the current in terms of the charge
$\delta n$ up to subleading order in derivatives

$j_{i}(t,x)=-D\partial_{i}\delta n+D_{2}\nabla^{2}\partial_{i}\delta
n+\cdots\,.$ (29)

We have omitted a term $\partial_{t}\partial_{i}\delta n$ which would have the
same scaling as $D_{2}$, because it can be absorbed in $D_{2}$ using the
leading equations of motion $\partial_{t}\delta n=D\nabla^{2}\delta n$. To
obtain response functions, one needs to know the constitutive relation in the
presence of a source $\delta\mu(t,x)$ for charge density. Assuming that the
equilibrium response is given by $\delta n(q)=\chi(q)\delta\mu(q)$, with
$\chi(q)\simeq\chi+\chi_{2}q^{2}+\cdots$ the static susceptibility, the
current constitutive relation in the presence of sources must take the form
(in momentum space)

$\begin{split}j_{i}(t,q)=&-iq_{i}(D+D_{2}q^{2})(\delta
n-(\chi+\chi_{2}q^{2})\delta\mu)\\\ &+\gamma
iq_{i}\partial_{t}\delta\mu+\cdots\,,\end{split}$ (30)

where the combination $\delta n-\chi(q)\delta\mu$ in the first line is
required for the current to vanish in thermal equilibrium. Note however that
this argument allows for terms involving the time derivative of the source, as
in the second line. Inserting the current in the continuity relation
$\partial_{t}n+\partial_{i}j_{i}=0$ and solving for $\delta n$ yields a
retarded Green’s function $G^{R}(\omega,q)\equiv\frac{\delta
n(\omega,q)}{\delta\mu(\omega,q)}$ that matches (28) with
$\chi_{2}=\tilde{c}_{1},\,D_{2}=-D^{2}c_{2},\,\gamma=\tilde{c}_{1}+(c_{1}+c_{2})\sigma$.

###  III.4 Multiple densities

Systems with multiple conserved densities can be studied similarly, by
including all densities in the EFT. The general scaling of loop corrections
remains unchanged; however, mixing of the densities allows for new scaling
functions with qualitatively different features. Indeed, consider for example
systems with conservation of both charge $\dot{n}+\nabla\cdot j=0$ and energy
(or heat) $\dot{\varepsilon}+\nabla\cdot j^{\varepsilon}=0$. Nonlinearities
can now involve both densities, e.g.:

$j_{i}=\cdots+\lambda\delta\varepsilon\partial_{i}\delta n+\cdots\,.$ (31)

The coefficient $\lambda$ arises from a temperature-dependent conductivity
$\partial_{T}\sigma$ (or, equivalently, a density-dependent thermoelectric
conductivity $\partial_{\mu}\alpha$). While it seems similar to the single
density nonlinearity $D^{\prime}$ considered above, this term is qualitatively
different because it is not a total derivative contribution to the current. It
therefore contributions to the $q=0$ optical conductivity $\sigma(\omega)\sim
1+\lambda^{2}|\omega|^{d/2}+\cdots$, as was already recognized in Ref.
Mukerjee  _et al._ (2006).

In order to obtain the universal scaling functions at finite $q$, necessary to
make predictions for the structure function $\langle n(t,x)n\rangle$, the EFT
is generalized to systems with multiple conservation laws in App. B. The
scaling functions are in this case complicated by the fact that there are
several diffusivities, and therefore several natural scaling variables
$y=x/\sqrt{Dt}$. To illustrate the appearance of novel scaling functions with
multiple densities in a simple context, we focus on the hydrodynamics of
densities for a non-abelian internal symmetry, say $SU(2)$. This situation is
simpler because the $SU(2)$ symmetry restricts the susceptibilities to be
diagonal $\chi_{AB}\equiv\frac{dn_{A}}{d\mu^{B}}=\chi\delta_{AB}$, leads to a
single diffusivity $D$, and only allows for one cubic nonlinearity in the EFT
which has a clear similarity with Eq. (31):

$j_{i}^{A}=-D\partial_{i}n^{A}+\lambda\epsilon_{ABC}n^{B}\nabla
n^{C}+\cdots\,.$ (32)

Here $A,B,C,\ldots$ run over the three elements of the $SU(2)$ algebra. The
hydrodynamic description of thermalizing systems with non-abelian internal
symmetries has been studied before Son and Stephanov (2002); Torabian and Yee
(2009); Fernandez-Melgarejo  _et al._ (2017); Glorioso  _et al._ (2021);
Claeys  _et al._ (2022), with the role of the nonlinearity $\lambda$
particularly emphasized in Glorioso  _et al._ (2021); Claeys  _et al._ (2022).
In Sec. B, we show that the 1-loop correction to the density two-point
function is, in one spatial dimension

$\displaystyle G^{R}_{n^{A}n^{B}}(\omega,q)=\delta_{AB}G^{R}(\omega,q)\,,$
(33) $\displaystyle G^{R}(\omega,q)=\frac{\sigma
q^{2}}{Dq^{2}-i\omega}+\frac{\lambda^{2}\chi^{2}}{D}i\omega
q^{2}\frac{\sqrt{q^{2}-\frac{2i\omega}{D}}}{\left(Dq^{2}-i\omega\right)^{2}}+\cdots\,.$

This produces a non-analytic correction at small frequencies to the optical
conductivity

$\begin{split}\sigma(\omega)&=\lim_{q\to
0}\frac{-i\omega}{q^{2}}G^{R}(\omega,q)\\\
&=\sigma-\frac{\lambda^{2}\chi^{2}}{D^{3/2}}(1-i)\sqrt{\omega}+\cdots\,.\end{split}$
(34)

The correction to the density two-point function in spacetime domain can be
found by Fourier transforming (see Sec. B). One finds a correction similar to
(1), with a different universal scaling function

$\displaystyle F_{1,0}^{\rm mult}(y)$
$\displaystyle=\frac{\chi\lambda^{2}}{D^{5/2}}\tilde{F}_{1,0}^{\rm mult}(y)$
(35) $\displaystyle\tilde{F}_{1,0}^{\rm mult}(y)$
$\displaystyle=\frac{4-y^{2}}{4\sqrt{\pi}}e^{-y^{2}/2}+\frac{y(2-y^{2})}{8}e^{-y^{2}/4}{\rm
Erf}(y/2)\,.$

###  III.5 Confirming the EFT with classical numerics

Before turning to quantum simulations, where the limited accessible time
scales make it crucial to account for power-law corrections to diffusion, we
confirm the EFT predictions in classical thermalizing systems. We focus on
classical lattice gases satisfying the ‘gradient condition’, namely where the
current density is a total derivative microscopically. In these situations,
the diffusivity $D(n)$ is known analytically Spohn (2012), making it simple to
perform precision tests of EFT predictions Delacretaz and Mishra (2023).
Indeed, since the loop corrections (1) and (4) only depend on the
susceptibility $\chi$ and derivatives of $D(\rho)$, they are entirely fixed
analytically and can be directly compared to numerics.

Figure 3:  (A): Profile of the dynamical structure factor for the KLS model
with parameters $\delta=\rho=0.9,\epsilon=0$. Different coloured curves denote
different times $t\in\\{200,2000\\}$ with smaller times corresponding to
darker colours. The red dashed curve is the diffusive prediction
$\frac{\chi}{\sqrt{4\pi D}}F_{0,0}(y)$. Inset: Autocorrelation function
($y=0$). Diffusive predictions with (black), and without (red), leading order
corrections, Eq. (3) are shown. (B): Comparison between the correction to
diffusion from simulation data $\Delta n\equiv\sqrt{t}\left[\langle
n(x,t)n\rangle-\frac{\chi}{\sqrt{4\pi Dt}}F_{0,0}(y)\right]$ and the EFT
prediction $F_{1,0}(y)$ (black, dashed), Eq. (1). Inset: Absolute area between
the finite time curves and the analytic expression,
$\mathcal{F}=\int^{3}_{0}dy|\Delta n(t)-F_{1,0}|/\int^{3}_{0}dy|F_{1,0}|$.

As a simple example of a lattice gas satisfying the gradient condition with a
non-trivial $D(\rho)$, we consider the one-dimensional Katz-Lebowitz-Spohn
model Katz  _et al._ (1984); Hager  _et al._ (2001); Baek  _et al._ (2017),
describing a collection of hard core particles hopping on a lattice with rates
depending on the occupation of neighbors:

$\displaystyle 0\,1\,0\,0\quad\xrightarrow{r(1+\delta)}\quad 0\,0\,1\,0$ (36a)
$\displaystyle 1\,1\,0\,1\quad\xrightarrow{r(1-\delta)}\quad 1\,0\,1\,1$ (36b)
$\displaystyle 1\,1\,0\,0\quad\xrightarrow{r(1+\epsilon)}\quad 1\,0\,1\,0$
(36c) $\displaystyle 0\,1\,0\,1\quad\xrightarrow{r(1-\epsilon)}\quad
0\,0\,1\,1$ (36d)

with equal rates for the spatially reversed processes. $\delta$ and $\epsilon$
are two parameters of the model, whereas $r$ defines the unit for time and can
be set to unity. We focus on the model with $\epsilon=0$, corresponding to
infinite temperature dynamics, which allows to use a random initial state as a
thermal state (taking $\epsilon\neq 0$ instead requires prethermalizing the
system, making numerics more costly). In this situation, the susceptibility is
$\chi(\rho)=\rho(1-\rho)$ and diffusivity $D(n)=1+\delta(1-2\rho)$, so that
$D^{\prime}=-2\delta$. This fixes all parameters entering in the leading
correction to diffusion, Eqs. (1), (2). Fig. 3 shows the excellent agreement
between the EFT prediction and numerics. We stress that the entire scaling
function agrees quantitatively with the 1-loop prediction $F_{1,0}(y)$, and
that no fitting parameter is involved in this comparison.

Figure 4:  Tensor network simulations of driven XXZ chain with (A,B)
decoherence, $\gamma=0.1$, and (C,D) staggered field, $g=0.4$. (A,C)
Diffusivity as a function of equilibrium magnetization. Top: different bond
dimensions, $d$. Bottom: different simulated times $T$, using the same bond
dimension $d=256$. Different colours denote fitting methods which take into
account an increasing amount of corrections to leading order diffusion,
I$\rightarrow$ II$\rightarrow$III, see Appendix (C.2). (B,D) Corrections to
diffusion evaluated at different times $t=(40,80,\ldots,400)$, denoted by dark
blue$\rightarrow$yellow colors. Black dashed line denotes linear corrections
F0,1. Red dashed line denotes F1,0 in (D), and F${}_{0,1}+\sqrt{T}$F1,0 with
$T=400$ in (B). The brown dashed line in (D) denotes the combined effect of
2-loop correction and linear corrections F${}_{0,2}+$F1,0, where the
logarithmic in time component of F0,2 is evaluated at $T=400$.

Analytical knowledge of the leading order correction allows for considerably
improved predictions for diffusivity when the available integration time is
short, i.e. when 1-loop effects are strong. Both experiments and simulations,
especially in quantum systems, are usually limited to relatively short
timescales, and uncontrolled extrapolations are therefore employed to get
infinite time properties such as diffusivity. We propose a robust method which
takes into account 1-loop effects by fitting the three-dimensional dataset
$\langle n(x,t)n\rangle_{\rho}$ versus $t,x,\rho$ with Eq. (1). We compare
this method to the fit using just the leading order term (equivalently taking
$F_{1,0}=0$). For example, we aim to approximate $D(\delta=0.9)=1.9-1.8\rho$
around $\rho=0.9$. For that, we simulate a sample of densities
$\rho=(0.85,0.86,\ldots 0.93)$. The diffusivity is parametrized by $D_{\rm
fit}=a-b\rho$ since we know its analytic form. In general, the parametrization
may include additional powers of density as the precise form of diffusivity is
not polynomial. We constrain the time window $t=40-T$, $T=100$ and find the
deviations of the fit from the exact diffusivity:
$\frac{1.9-a}{1.9}\approx(0.135,0.081)$ and
$\frac{-1.8-b}{-1.8}\approx(0.12,0.084)$, where the first/second number in the
parenthesis denote the fit without/with the 1-loop correction. Our results
show a quantitative improvement which increases as $T$ decreases and 1-loop
effects become stronger. In principle, it is possible to perform time
extrapolations to the above fitting method by using time windows of varying
size.

##  IV Quantum Transport

To test and make use of the EFT predictions in a minimal setting, we focus on
quantum-coherent and incoherent chaotic systems with a single conserved
charge. We estimate the dependence of diffusivity on the equilibrium
magnetization using various approaches, and show that the corrections to
diffusion are in agreement with EFT. Moreover, incorporating them in the
fitting methods can significantly improve the diffusivity approximation at
finite times.

###  IV.1 Model and Methods

The conserved charge is chosen to be magnetization,

$N=\sum^{L}_{i=1}\sigma^{z}_{i},$ (37)

where $\sigma^{z}$ is the Pauli-z matrix and $L$ is the system size. Since
magnetization is a sum of local operators, the equilibrium ensemble is a
product state, i.e., $e^{\mu N}=\prod^{L}_{i=1}e^{\mu\sigma^{z}_{i}}$. We did
not consider systems with charges such as energy density which have
equilibrium states with a finite correlation length because this adds an
additional layer of complexity to the simulations, even though the EFT
predictions are the same.

As a minimal chaotic model where magnetization is the only conserved quantity,
we choose the Floquet-XXZ chain with a staggered field, whose stroboscopic
dynamics are generated by a Floquet operator

$\displaystyle U=U_{e}U_{o},\quad
U_{e}=\prod^{L/2}_{i=1}U_{2i,2i-1},\;U_{o}=\prod^{L/2}_{i=1}U_{2i-1,2i-2}$
$\displaystyle
U_{i,i+1}=\text{exp}\left(-i(h^{XXZ}_{i,i+1}+V_{i,i+1})\right).$ (38)

The evolution is performed by first evolving the odd bonds and then the even
bonds, with two-body gates generated by the following operators,

$\begin{split}&h^{XXZ}_{i,i+1}=J(\sigma^{+}_{i}\sigma^{-}_{i+1}+h.c.)+\frac{\Delta}{2}\sigma^{z}_{i}\sigma^{z}_{i+1},\\\
&V_{i,i+1}=g\left((-1)^{i}\sigma^{z}_{i}+(-1)^{i+1}\sigma^{z}_{i+1}\right).\end{split}$
(39)

In the absence of a staggered field $V$ the Floquet-XXZ chain is integrable
Gritsev and Polkovnikov (2017). For our choice of parameters,
$J=\pi/4,\Delta=J-0.2$, the magnetization at $V=0$ displays ballistic
transport. Turning on staggered field leads to integrability breaking, the
system becomes chaotic, and therefore magnetization is expected to diffuse. To
establish our method we alternatively perturb the XXZ-chain with Markovian
noise (dephasing). Dephasing effectively suppresses the generation of operator
entanglement in the simulation, leading to very accurate numerical data. To
simulate the noise-averaged state, we define the dephasing map by the action
of the local channel on the state of a single spin,

$\qquad\mathcal{D}_{i}\left(\begin{array}[]{cc}\rho_{1,1}&\rho_{1,0}\\\
\rho_{0,1}&\rho_{0,0}\end{array}\right)=\left(\begin{array}[]{cc}\rho_{1,1}&e^{-\gamma}\rho_{1,0}\\\
e^{-\gamma}\rho_{0,1}&\rho_{0,0}\end{array}\right).$ (40)

The global noise channel is a product of local channels and is applied to the
state following a period of coherent evolution,

$\mathcal{D}=\bigotimes^{N}_{i=1}\mathcal{D}_{i},\quad\rho(t+1)=\mathcal{D}(U\rho(t)U^{{\dagger}}).$
(41)

To study the linear response dynamics, we employ the weak domain wall initial
state proposed by Ljubotina et al. Ljubotina  _et al._ (2019),

$\rho(\mu,\delta,t=0)=\frac{1}{M}e^{\mu
N}\left(\prod^{L/2}_{i=1}e^{-\delta\sigma^{z}_{i}}\bigotimes\prod^{L}_{i=\frac{L}{2}+1}e^{\delta\sigma^{z}_{i}}\right),$
(42)

where $M=\text{tr}\rho$ is the normalization constant and $\delta\rightarrow
0$ generates a weak domain wall perturbation on top of the equilibrium state
characterized by the chemical potential $\mu$. The linear-response regime is
characterized by a quench where the amount of injected magnetization is not
extensive. According to Lieb-Robinson bounds, in the worst case scenario,
$\delta$ will satisfy this criterion provided $\delta\times v_{LR}t_{\rm
max}\ll 1$, where $t_{\rm max}$ is the maximum simulation time and $v_{LR}$ is
the Lieb-Robinson velocity of quantum correlations propagation. We use
$\delta=0.0005$, which is in the linear-response limit for all simulated
times.

The Floquet evolution defined by Eq (IV.1) breaks translation invariance,
since even and odd sites are not equivalent. To simplify analysis, we average
over even and odd sites. In addition, we shift the magnetization by its
equilibrium value,
$\overline{\sigma}=\text{tr}\left(\sigma^{z}_{i}\rho(\mu,0)\right)=\tanh{\mu}$
and normalize its initial magnitude to $1/2$,

$s_{j}=\frac{\text{tr}\left(\sigma^{z}_{2j-1}\rho\right)+\text{tr}\left(\sigma^{z}_{2j}\rho\right)-2\overline{\sigma}}{4|\text{tr}\left(\sigma^{z}_{1}\rho(t=0)\right)|},$
(43)

where $j\in\\{1,2,\ldots,\overline{L}\\}$, and $\overline{L}=L/2$. In this
normalization, the initial state profile is $s_{j\leq\overline{L}/2}=-0.5$ and
$s_{j>\overline{L}/2}=0.5$ for all values of $\mu$ and $\delta$. Since we
doubled the lattice spacing, diffusivity and static susceptibility,
$\chi=\frac{1}{L}\left(\langle N^{2}\rangle-\langle
N\rangle^{2}\right)=(\cosh\mu)^{-2}$, are rescaled accordingly, $D\rightarrow
D/4$, $\chi\rightarrow\chi/2$. We will always present the results of
diffusivity computed using the original lattice spacing. For clarity, we
employ a continuum description of the lattice variables $s_{j}\rightarrow
s(x)$, since the hydrodynamic corrections are defined in continuum.

A simple relation between the domain-wall quench for $\delta\rightarrow 0$ and
the dynamical structure factor was originally derived in Ref. Ljubotina  _et
al._ (2017). In the continuum limit, and under the conventions described in
the previous paragraph, the relation simply reads
$\langle\sigma_{z}(x,t)\sigma_{z}\rangle_{c}=\chi\frac{ds}{dx}$, where the
subscript $c$ stands for the connected part of the correlation function and
the average is performed over the equilibrium state, $\rho(\mu,\delta=0)$.
Distance $x$ is measured from the position of the domain wall at
$\overline{L}/2$. While the two quenches are formally equivalent, we prefer to
transform (simple integral in space) the EFT results to the domain-wall
picture because taking the spatial derivative (or directly calculating the
dynamical structure factor) of the numerical data enhances the errors
generated by the simulations. The corrections of Section II translate to,

$\begin{split}&\mathcal{S}(x,t)=\sum_{m,n}\frac{1}{t^{m/2+n}}\mathbf{F}_{m,n},\\\
&\mathbf{F}_{m,n}=\frac{1}{\sqrt{4\pi Dt}}\int_{0}^{x}dxF_{m,n}.\end{split}$
(44)

The leading order diffusion is given by
$\mathbf{F}_{0,0}=\frac{1}{2}\text{Erf}(y/2)$ and the corrections are
presented in Eq. (A.6). To analyze the numerical data, we propose a
generalized fitting procedure which takes into account the leading corrections
to diffusion, in order to extract the diffusivity as a function of the
equilibrium magnetization, $D(\overline{\sigma})$. This is achieved by
minimizing the difference between the numerical data and the theory on a
three-dimensional surface parametrized by $(\overline{\sigma},x,t)$,

$\text{Min}_{D,\vec{c}}\sum_{i,j,k}|s(\overline{\sigma}_{i},x_{j},t_{k})-\mathcal{S}(\overline{\sigma}_{i},x_{j},t_{k},D(\overline{\sigma}),\vec{c}(\overline{\sigma}))|^{2},$
(45)

where the subscripts $i,j,k$ denote the different samples in the discretized
data-set. The functional minimization over diffusivity profiles is simplified
by employing a Taylor expansion around half-filling,
$D(\overline{\sigma})=\sum^{M}_{i=0}b_{i}(\overline{\sigma})^{2i}$, where
$M=3$ is found to give converged results for the parameter regimes studied in
this work. Only even powers are allowed in the expansion due to the particle-
hole symmetry in the system ($D(\overline{\sigma})=D(-\overline{\sigma})$).
The parameters $\vec{c}=(c_{1},c_{2},\ldots)$ are non-universal parameters
arising from linear-fluctuations described in section III.3. Since we will
only use leading and subleading-order corrections, we just require the two
parameters $\vec{c}=(c_{1},c_{2})$ defined by Eq. (5) and are present in
$\textbf{F}_{0,1}$.

We compare our method to two other approaches which are described in Appendix
C. The fitting methods we use are: (I) A method in which the approximate
diffusivity is extracted by fitting the dynamical structure factor at the
largest time with the leading order profile $\textbf{F}_{0,0}$ . (II) The
method described above in Eq. (45), which takes into account the corrections
$\textbf{F}_{1,0},\,\textbf{F}_{0,1},\,\textbf{F}_{2,0}$ that were explicitly
computed in Sec. III. (III) A method which is based on the minimization of the
deviation of total current of the system from the expected generalized Fick’s
1st law. Many higher derivative and loop corrections cancel in III, so that
this method in fact accounts for more corrections to leading order diffusion
than II, without having to compute them explicitly (namely, all zero-loop and
one-loop higher-derivative corrections $\textbf{F}_{0,n}$ and
$\textbf{F}_{1,n}$, as well as all $\ell$-loop zero-derivative corrections
$\textbf{F}_{\ell,0}$ cancel in III, see Appendix C.2). However, it is a more
expensive method as it requires measuring the full current which is equivalent
to measuring the full structure factor. Method II, on the other hand, can be
equally efficient when few points of the structure factor are sampled.

To probe the strength of the non-linear corrections we calculate the prefactor
of F1,0 defined in Eq. (A.6),

$C_{1}=\frac{\chi D^{\prime 2}}{\sqrt{4\pi D^{5}}},$ (46)

where $D^{\prime}=\frac{dD}{d\overline{\sigma}}$ and the magnetization profile
after subtracting the leading order term,

$\Delta s(y)=s(y)-\textbf{F${}_{0,0}$}.$ (47)

In our simulations we fix the system size to $L=2T$ and maximum time to
$T=400$ in order to avoid finite size effects. We simulate chemical potentials
$\mu=(0,0.01,\ldots)$ which identify the equilibria magnetizations. The
density matrix is represented by a matrix product density operator (MPDO)
Verstraete  _et al._ (2004), and evolved using time evolving block decimation
(TEBD) Vidal (2004). The simulations are performed using the ITensor library
Fishman  _et al._ (2022). Unless stated otherwise the simulations employ bond
dimension $d=256$. As shown in Figure 4 and further explored in Appendix C our
results are not significantly affected by bond dimension. This allows us to
get sufficiently accurate results for $\mu\leq 0.25$ for the staggered
perturbation and for all $\mu$ in the case of dephasing. For, the staggered
field around half filling we have used $d=400$, as $d=256$ is not enough to
accurately extract the weak correction signal.

###  IV.2 Results

Dephasing–To establish our the method we switch off the staggered field
perturbation ($g=0$) and simulate the Floquet-XXZ chain, Eq. (41) in the
presence of dephasing with $\gamma=0.1$. Diffusive transport induced by
dephasing has two distinct features. First, the single particle limit
(equivalently the non-interacting limit $\Delta=0$) in the presence of
dephasing remains diffusive, and therefore, diffusivity is finite for all
magnetizations. This is in contrast to purely interaction-induced diffusion
where the single particle limit is ballistic (free particles). Second, for
increasing strengths of dephasing, magnetization transport becomes less
sensitive on the interactions in the strong noise limit where, to leading
order in $1/\gamma$, $D\propto J^{2}/\gamma$ Bauer  _et al._ (2017) is
independent of magnetization density $\bar{\sigma}$ and hence
$D^{\prime}\equiv\partial_{\bar{\sigma}}D\simeq 0$.

In Figure 4.A we show the dependence of diffusivity on magnetization,
$D_{\text{III}}\sim D_{\text{II}}\approx
5.35+3.03\overline{\sigma}^{2}+0.54\overline{\sigma}^{4}$. For the fit II we
have employed the terms
$(\mathbf{F}_{0,0},\mathbf{F}_{1,0},\mathbf{F}_{0,1})$. In contrast to the
non-interacting ($\Delta=0$) limit, diffusivity does have a magnetization
dependence due to the presence of non-linear corrections. However, the 1-loop
correction $\textbf{F}_{1,0}$ is small at timescales of order $t\sim O(100)$,
since $C_{1}\sim O(10^{-2})$. The methods II,III converge to almost the same
curve (independently of bond dimension), which is indicative of the fast
convergence to the asymptotic behaviour. Method I slightly underestimates the
asymptotic value. The reason is that (I) does not capture the linear
corrections $\textbf{F}_{0,1}$, which dominate at these timescales for all
fillings, despite being suppressed by a factor of $1/\sqrt{t}$ compared to the
1-loop corrections, Figure 4.C. The 1-loop effects only have a visible effect
at the largest simulated times. In all cases, the corrections $\Delta s$ are
accurately captured from our theory. Finally, we observe that the correction
profiles are almost time-independent, which is indicative that besides 1-loop
corrections, higher order corrections are also suppressed.

Unitary dynamics– Shifting to unitary dynamics, in Figure 4.C we show that the
staggered field generates a strong dependence of diffusivity to magnetization,
independently of the fitting method used, $C_{1}\propto O(1)$. However,
diffusivity is considerably affected by the method used as shown by the
$d=256,T=400$ fits: $D_{\text{II}}\approx
1.77+13.7\overline{\sigma}^{2}+220\overline{\sigma}^{4}$,
$D_{\text{III}}\approx 1.73+27\overline{\sigma}^{2}+190\overline{\sigma}^{4}$.
For the fit II we have employed the terms
$(\mathbf{F}_{0,0},\mathbf{F}_{1,0},\mathbf{F}_{0,1})$. This discrepancy is
not due to truncations in the dynamics; instead, we believe it arises from
fitting timescales that are not in the asymptotic diffusive regime, which is
reflected by the dependence of parameters yielded by each method on the
maximum timescale of the simulation. The methods are affected according to the
number of corrections to diffusion they include in the approximation. In that
sense, III is more converged than II, and II is more converged than I, which
employs no additional corrections to asymptotic diffusion. As expected, when
simulation is performed for longer times, different methods tend to show
better agreement.

Due to the discrepancy in the determination of $D$, the form of the
corrections $\Delta s$ depends on the fitting method. Here we have chosen to
use the diffusivity estimated by method II, since by construction it fits
these corrections arising from nonlinearities and higher derivative terms,
while in III the leading such effects are absent. In Appendix C we show that
the different methods result in different quantitative but not qualitative
profiles. Figure 4.D shows that for finite magnetization
$\overline{\sigma}=0.2$, $\Delta s$ scales as $t^{-1/2}$, and has a closely
matching profile to that of the expected 1-loop correction $\mathbf{F}_{1,0}$.
We again note that the 1-loop profile is completely determined by
($D,\partial_{\overline{\sigma}}D)$ and no additional fitting is involved. We
observe that the correction profile shows a significant time dependence,
indicative of higher corrections being still at work at these timescales. In
contrast to finite equilibrium magnetization, $\overline{\sigma}=0$ requires
special attention. First, the correction signal is weaker and requires a bond
dimension $d=400$ to be accurately captured. Additionally, the strength of
non-linear corrections suggests that the leading corrections,
$\mathbf{F}_{2,0}$ and $\mathbf{F}_{0,1}$ will be of similar magnitude at
intermediate timescales. Eventually $\mathbf{F}_{2,0}$, will dominate due to
its logarithmic divergence in time. For this reason we perform a different fit
II around $\overline{\sigma}=0$, by including a few points
$\sigma\in(0,0.01,0.02,0.03)$ and all terms
$(\mathbf{F}_{0,0},\mathbf{F}_{1,0},\mathbf{F}_{0,1},\mathbf{F}_{2,0})$.
Figure 4.D indeed shows that both $\mathbf{F}_{2,0}$ and $\mathbf{F}_{0,1}$
are important at these timescales. However, the available time scales are not
sufficient to see the ultimate dominance of the logarithmic part of
$\mathbf{F}_{2,0}$, which would lead to a profile with an opposite sign around
$y=-0.2,0.2$.

Overall, we have shown that employing EFT corrections to study quantum
transport can significantly improve the estimation of asymptotic transport
parameters such as diffusivity. In addition, these corrections help
understanding the different processes that drive a system towards equilibrium.
For example, dephasing, which is often used to accelerate thermalization,
achieves this at the cost of flattening the diffusivity as a function of
filling or magnetization. We have moreover found that even if noise
$\gamma<\Delta,J$ the system’s behavior is similar to the strong noise limit
$\gamma\gg\Delta,J$, where the equation of motion for the conserved charge can
be perturbatively derived Bauer  _et al._ (2017); Claeys  _et al._ (2022). In
that case, the diffusivity has weak dependence on equilibrium magnetization
$\overline{\sigma}$, leading to smaller loop corrections and faster
thermalization: indeed the leading order non-linear terms in the diffusivity
only appear at 3rd order $O(\Delta^{2}J^{2}/\gamma^{3})$. In our situation,
even if perturbation theory in $1/\gamma$ is not strictly valid , we observe
the same behaviour: Weak non-linear corrections, F1,0, which are almost
invisible at the simulated timescales, and a fast approach to asymptotic times
which is driven by the sub-leading linear corrections F0,1.

Staggered perturbations on the other hand, induce strong non-linear effects,
leading to a slower approach to equilibrium. Additionally, the classical
resources required to simulate the system increase rapidly with the simulation
time, and therefore, the accessible timescales are limited. We have shown that
employing fitting methods which take into account the EFT corrections to
diffusion leads to a significant improvement in the diffusivity estimation for
$\overline{\sigma}>0$. At the same time our numerical data strongly suggest
that the EFT corrections to the dynamical structure factor are present in
interacting quantum systems.

##  V Non-linear response

In the previous section we explored the effects of scaling corrections in the
DSF which is a two-point correlation function. While the effect of
fluctuations on the DSF can be important, it is always subleading at long
times. Here, we go one step further and explore many-body correlations which
would vanish in the absence of nonlinearities in hydrodynamics. In particular
we explore the late-time behavior of the observable,

$s_{3}(x,t)=\frac{\langle\sigma^{z}\left(L/2,t\right)\sigma^{z}\left(L/2+x,t\right)\rangle_{c}}{|2\text{tr}\left(\sigma^{z}_{1}\rho(t=0)\right)|},$
(48)

where the average is performed on the domain wall state defined by Eq. (42).
As we illustrate in Appendix A.6, this observable is, up to a spatial
integration, a special case of the density three-point function. Higher-point
functions generalize full counting statistics in that they allow for operator
insertions at multiple times, and can be obtained from the EFT of diffusion
Delacretaz and Mishra (2023). According to the EFT, the asymptotic behavior of
$s_{3}(x,t)$ in diffusive systems with a single conserved charge is given by a
universal scaling function

$s_{3,\text{EFT}}(y,t)=\frac{1}{\sqrt{t}}\frac{D^{\prime}\chi}{8\sqrt{D^{3}\pi}}\left(e^{-y^{2}/4}+\text{Erf}(y/2)-1\right),$
(49)

where $y=x/\sqrt{Dt}$. As expected for multi-body functions in linear
response, $s_{3,\text{EFT}}(y,t)$ is vanishes as $t\rightarrow\infty$. It also
vanishes in systems which are particle-hole symmetric ($D^{\prime}=0$).

As shown in Figure 5, for both dephasing and staggered perturbations, $s_{3}$
scales according to the EFT prediction, $1/\sqrt{t}$. In the presence of
dephasing, we observe a precise late time agreement of the correction profile
to $s_{3,\text{EFT}}$, which verifies the validity of the EFT prediction. No
fitting parameter was used in this test, as $D$ and $D^{\prime}$ had already
been obtained from the linear response analysis. In the case of staggered
field, we find that while the shape of the profile is qualitatively similar,
there is a quantitative deviation from the EFT prediction (independent of the
diffusivity fitting method). The similarity between the EFT profile (49) and
the numerical data can be illustrated by performing a fit of
$s_{3,\text{EFT}}$ for $y>1$, with now $D$ and $D^{\prime}$ taken to be
independent fitting parameters. The result is illustrated in the inset of
Figure 5.B, the two profiles agree well for $(D\sim
1.85D_{\text{III}},D^{\prime}\sim 0.36D^{\prime}_{\text{III}})$. This
agreement suggests that the staggered profile at $y\gg 1$ has the same
functional form as the EFT prediction and features a Gaussian tail, despite
the apparent disagreement with the values of $D$ and $D^{\prime}$ obtained
from linear response.

While the discrepancy between the EFT prediction and quantum dynamics in
staggered field simulations remains unclear to us, we have verified that this
disagreement is not due to simulation errors. Additionally, the good
convergence with time suggests that subleading corrections are not at play.
Taking into account these observations leaves us with several possibilities:
(i) The diffusivity fits of the previous section may not be accurate for the
staggered field simulations. This could be an artifact of the short accessible
time scale or of large diffusivity fluctuations with magnetization. We have
extensively checked for these artifacts and we did not find any evidence of
them in the data. (ii) This system may simply have a very long local
equilibration time, $\tau_{\rm eq}\gtrsim 200$, so that our numerics never
probes the truly asymptotic regime controlled by the EFT. This could e.g.
arise from the presence of additional long-lived degrees of freedom protected
by approximate symmetries, and could be related to the integrability of the
Floquet-XXZ chain, or to other prethermalization mechanisms Abanin  _et al._
(2017). (iii) Finally, the EFT may fail to capture even the asymptotic
dynamics of coherent many-body Floquet systems. While we do not have a
particular reason to believe that the EFT should fail, there is no proof that
it must emerge in general. We leave this interesting possibility, as well as
further exploration of these phenomena, for future work.

Figure 5:  Three-point correlation at linear response. (A) Noisy system with
bond dimension $d=256$, $\gamma=0.1$ at $\mu=0.65$. (B) Staggered field with
bond dimension $d=600$, $g=0.4$ at $\mu=0.23$. Different coloured curves
denote different times $t\in\\{40,400\\}$ for (A) and $t\in\\{40,200\\}$ for
(B) with smaller times corresponding to darker colours. Red dashed lines
corresponds to EFT predictions with diffusivity estimated by method III.
Crosses in (B) show bond dimension $d=400$ results at $t=200$. Inset: The tail
of the distribution at $t=200$ against the EFT function (black dashed) where
$(D,D^{\prime})$ are taken as fitting parameters.

##  VI Discussion

We have employed the EFT of diffusion to derive the scaling functions of the
leading power law corrections to diffusive transport for thermalizing systems
with one or more conserved local charges. We confirmed these predictions by
numerical simulations in a classical model, finding percent-level agreement of
the entire scaling function without any fitting parameter (see Fig. 3). While
testing subleading EFT predictions in quantum simulations with this same level
of precision is currently beyond reach, due to the rapid growth of required
classical resources, these corrections are expected to be particularly
important there due to the shorter accessible timescales. We showed that
knowledge of these corrections allows for more accurate extraction of
transport parameters, especially when the accessible timescales are very
limited. Furthermore, our results open a number of promising directions for
future research; we list these and other applications below.

#### Precision tests of thermalization—

Our findings can also be used to test possible deviations from standard
diffusion in numerics and experiments. For example, tracking the density (or
temperature) dependence of transport parameters can help estimate power law
corrections to observables. Given that these corrections typically make
diffusive systems appear superdiffusive at intermediate times, it would be
interesting to study them quantitatively in the context of 1d chains showing
apparent anomalous diffusion of superdiffusion De Nardis  _et al._ (2020);
Ljubotina  _et al._ (2023); Chen  _et al._ (2023) (see Glorioso  _et al._
(2021) for preliminary work in this direction, and Wienand  _et al._ (2023)
for related work), as well systems featuring subdiffusion without dipole
conservation Richter and Pal (2022); Gopalakrishnan  _et al._ (2022). Higher-
point functions of local operators offer useful information in this regard.
Indeed, we showed in Sec. V that these are controlled by the same EFT
parameters (in particular, $D^{\prime}$) which lead to power law corrections
to linear response at intermediate times. Measuring higher-point functions of
density (or heat) therefore provides a time-scale that must be exceeded to
access the asymptotic dynamics. The EFT also points to other observables, such
as correlation functions in momentum space $\langle n(q,t)n\rangle$, which
instead are not suitable for precision tests of thermalization because they
receive large fluctuation corrections Delacretaz (2020).

#### Beyond diffusion—

Our theoretical results can be extended in various directions. We have assumed
the dissipative fixed point to be diffusive, however one can similarly study
corrections to subdiffusive or superdiffusive universality classes, or even in
generalized hydrodynamics for integrable models (for the KPZ universality
class, the leading scaling correction was studied in Takeuchi and Sano (2010);
Sasamoto and Spohn (2010); Ferrari and Frings (2011)). These corrections are
also important to incorporate for quantum simulations in higher dimensions
$d>1$, where our ability to numerically study large systems and times is
limited. Our results for the 1-loop and 2-loop corrections, Eqs. (20), (24),
hold in any dimension.

#### Connections to simulation complexity—

We believe that EFT corrections present new hints towards understanding the
hardness of quantum simulations in the linear response regime. We found that
classical resources increase with time faster when non-linear corrections are
stronger, in our case when $\overline{\sigma}$ increases. This obstructed
exploring magnetizations beyond $\overline{\sigma}\approx 0.25$ in the
staggered field simulations. We lack a detailed theory behind this
observation, but believe that this is related to the strong non-linear
corrections, since they enhance multi-point correlation functions such as the
ones explored in Sec. V. This implies that the accurate simulation of a system
with strong non-linear contributions requires keeping more information on
multi-body correlations in the density matrix, which in turn increases the
resources (bond dimension) required by the tensor network simulations.

#### Benchmark for new methods—

Our results on universal corrections to hydrodynamics are also useful to
benchmark theoretical and computational White  _et al._ (2018); Parker  _et
al._ (2019); Claeys  _et al._ (2022); Klein Kvorning  _et al._ (2022);
Rakovszky  _et al._ (2022) approaches to thermalization in many-body systems,
as these will have to reproduce not only the leading diffusive behavior but
the corrections as well. For example, Ref. Parker  _et al._ (2019)
approximates correlation functions based on extrapolations of Lanczos
coefficients, which by design produces a meromorphic $G^{R}(\omega,q)$ that
cannot capture the universal non-analytic corrections (22). Incorporating EFT
results into such constructions is a promising path to ‘bootstrapping’
transport in correlated quantum systems.

#### Distinguishing theories of fluctuating hydrodynamics—

We have also shown that high precision classical stochastic simulations offer
valuable precision tests for theories of fluctuating hydrodynamics, in the
present case confirming the leading and subleading corrections predicted by
the Schwinger-Keldysh EFT approach Crossley  _et al._ (2017); Chen-Lin  _et
al._ (2019). Other approaches for fluctuating hydrodynamics exist, which treat
the noise fields somewhat differently; it would be interesting to further push
these tests, to possibly rule out certain theories and identify the correct
systematic framework. One possible concrete target for the numerics in this
regard are effects arising from non-Gaussianities in noise fields that do not
enter in constitutive relations Jain  _et al._ (2021).

## Acknowledgements

We thank Paolo Glorioso and Zlatko Papic for insightful discussions. This work
was partially supported by the European Research Council via the grant
agreements TANQ 864597 (A.M., D.A.) and by the Swiss National Science
Foundation (D.A.).

## Appendix

##  Appendix A EFT details

###  A.1 General structure of the corrections

The corrections to diffusion shown in Eq. (1) and (4) correspond to the first
few terms arising from an expansion in fluctuations and derivatives in the
EFT. The general form that this expansion takes, in $d$ spatial dimensions,
is:

$\begin{split}\langle n(x,t)n\rangle=\frac{\chi}{(4\pi Dt)^{d/2}}\bigg{[}\ \
&\ F_{0,0}(y)+\frac{1}{t}F_{0,1}(y)+\frac{1}{t^{2}}F_{0,2}(y)+\cdots\\\
+\frac{1}{t^{d/2}}&\left(F_{1,0}(y)+\frac{1}{t}F_{1,1}(y)+\frac{1}{t^{2}}F_{1,2}(y)+\cdots\right)\\\
+\frac{1}{t^{d}}&\left(F_{2,0}(y)+\frac{1}{t}F_{2,1}(y)+\frac{1}{t^{2}}F_{2,2}(y)+\cdots\right)+\cdots\
\bigg{]}\,,\end{split}$ (50)

where $F_{\ell,n}$ are scaling functions of the scaling variable $y\equiv
x/\sqrt{Dt}$. The overall form of the expansion (50) is simple to justify on
general grounds: higher derivative corrections to diffusion come with two
derivatives (assuming reflection or rotation symmetry), and therefore give
corrections suppressed by $\nabla^{2}\sim\frac{1}{x^{2}}\sim\frac{1}{t}$ at
late times. Loop corrections instead come from nonlinearities in the dynamics
of the densities: a single cubic nonlinearity is suppressed by $\delta n\sim
q^{d/2}\sim 1/t^{d/4}$ compared to the linear (Gaussian) dynamics. The first
loop correction requires two insertions of a cubic nonlinearity, and is hence
$1/t^{d/2}$ suppressed. Generalizing, an $\ell$-loop contributions at $n$th
order in the derivative expansion will give a correction to correlation
functions suppressed by $1/t^{n+\ell d/2}$ (up to logarithms); this correction
comes with a dimensionless scaling function $F_{\ell,n}$ and is shown in the
$\ell$th line and $n$th column in Eq. (50). This general structure of
corrections to hydrodynamics applies not only to density two-point functions,
but also to higher point functions Delacretaz and Mishra (2023), as well as to
correlators of arbitrary microscopic operators that have the same quantum
numbers as (composites of) densities Delacretaz (2020); Glorioso  _et al._
(2021).

While the simple scaling argument above predicts the general expansion of
correlation functions at late times in diffusive systems, obtaining the
dimensionless scaling functions $F_{\ell,n}$ in Eq. (50) requires detailed use
of the EFT. The leading diffusive scaling function is well known
$F_{0,0}(y)=e^{-y^{2}/4}$, and captures the density two-point function
universally in any diffusive system. The subsequent $F_{\ell,n}$ capture
scaling corrections to diffusion; they are also universal, up to one or a few
theory-dependent factors.

In this paper, we focus on the first few corrections and explicitly evaluate
$F_{1,0}$ (Eq. (2)), $F_{2,0}$ (Eq. (6)), and $F_{0,1}$ (Eq. (5)). Note that
in (50), we have suppressed certain factors of $\log t$ that can arise from
loop corrections coming with integer powers of $1/t$, see (6) for an example.

###  A.2 Details of the EFT

In this section, we further motivate the EFT representation of the generating
functional (11), and go over several key steps in the construction of the EFT.
Most of the discussion in this section can be found elsewhere, e.g. Refs.
Crossley  _et al._ (2017); Kamenev (2023), but we include it for completeness.

One of the guiding principles in constructing to the Schwinger-Keldysh EFT for
hydrodynamics is to introduce a minimal set of fluctuating degrees of freedom
that will ensure gauge invariance of the generating functional

$Z[A^{1}_{\mu}+\partial_{\mu}\lambda^{1},A^{2}_{\mu}+\partial_{\mu}\lambda^{2}]=Z[A^{1}_{\mu},A^{2}_{\mu}]\,.$
(51)

This is done by introducing phases $\phi^{1,2}$ that always enter in the
combination $A_{\mu}+\partial_{\mu}\phi$ (sometimes called the ‘Stückelberg
trick’):

$Z[A^{1}_{\mu},A^{2}_{\mu}]=\int D\phi^{1}D\phi^{2}\,e^{i\int
dtd^{d}x\mathcal{L}[A_{\mu}^{1}+\partial_{\mu}\phi^{1},A_{\mu}^{2}+\partial_{\mu}\phi^{2}]}\,.$
(52)

This now satisfies (51) for any functional $\mathcal{L}$, because a gauge
transformation can be absorbed through a redefinition of the dynamical fields
$\phi^{1,2}$ that are being integrated over. It is clear that the degrees of
freedom we have introduced are related to the continuous symmetry of the
system. If one had considered instead a system with $N$ separate continuity
relations (8), $2N$ fields would have been introduced. One already notices a
ressemblence with earlier approaches to fluctuating hydrodynamics, where each
continuity relation leads to two degrees of freedom: a density and an
associated noise field. The central assumption in the construction of the EFT
is that $\mathcal{L}$ is a local functional of the fields $\phi^{1,2}$. This
implements the expectation of thermalization: the only long-lived quantities
are associated with symmetries, so that integrating out other degrees of
freedom produces a local EFT with a derivative expansion controlled by the
scale at which the system thermalizes (this scale acts as the UV cutoff of the
hydrodynamic EFT).

It is useful to define the symmetric and antisymmetric combinations of fields
Kamenev (2023)

$\phi_{r}\equiv\frac{\phi_{1}+\phi_{2}}{2}\,,\qquad\phi_{a}\equiv{\phi_{1}-\phi_{2}}\,.$
(53)

One advantage of fields in this basis is that they satisfy the ‘latest time’
property: correlators where the latest time is carried by a $\phi_{a}$ field
vanish due to cyclicity of the trace in Eq. (9)

$\langle O_{1}(t_{1})\cdots O_{n}(t_{n})\phi_{a}(t_{n+1})\rangle=0\,,\quad
t_{n+1}>t_{i}\,,$ (54)

which will lead to simplifications in diagrams below.

The construction bears a resemblance with EFTs for spontaneously broken
phases, where the long-lived degrees of freedom are Goldstone bosons. In fact,
constructing the most general local Lagrangian $\mathcal{L}$ in (52) leads to
an effective description of a thermalizing system in the symmetry broken phase
(a dissipative superfluid). To describe the normal phase, Ref. Crossley  _et
al._ (2017) proposes to forbid the propagating sound mode by imposing an
additional symmetry:

$\phi_{r}(x,t)\to\phi_{r}(x,t)+\lambda(x)\,.$ (55)

See Refs. de Boer  _et al._ (2019); Glorioso  _et al._ (2018) for discussions
on this symmetry in a holographic context. Recently, Ref. Akyuz  _et al._
(2023) proposed a slightly different approach that bypasses the need to impose
this somewhat artificial symmetry, by viewing the density $n_{r}$ rather than
$\phi_{r}$ as the fundamental degree of freedom of the EFT. We expect both of
these approaches to be equivalent.

Otherwise, one simply constructs the most general local functional of the
gauge invariant combinations of fields

$B_{\mu}^{1}\equiv A_{\mu}^{1}+\partial_{\mu}\phi^{1}\,,\qquad
B_{\mu}^{2}\equiv A_{\mu}^{2}+\partial_{\mu}\phi^{2}\,,$ (56)

in an expansion in fields and derivatives. There are few additional
constraints to impose, such as unitarity $Z[A_{1},A_{2}]^{*}=Z[A_{2},A_{1}]$
(which simply follows from the definition (9)), and KMS symmetry; we refer the
reader to Ref. Crossley  _et al._ (2017) for details. To leading order in
derivatives, the action can be expressed as Delacretaz and Mishra (2023)

$\mathcal{L}=\sigma(n)B_{ai}\left(iB_{ai}-\beta
E_{ri}\right)+B_{a0}n-D(n)B_{ai}\partial_{i}n+\cdots\,,$ (57)

where $\beta$ is the inverse temperature,
$E_{ri}=\partial_{0}A_{ri}-\partial_{i}A_{r0}$ is the electric field, and
$B_{a\mu}=A_{a\mu}+\partial_{\mu}\phi_{a}$. We have changed variables from
$\phi_{r}$ to the density $n$. The ellipses denote higher derivative terms
(the most important of which are discussed separately in Sec. III.3), as well
as nonlinear terms that contain higher powers of $\phi_{a}$, which are more
irrelevant than the nonlinearities considered here.

Setting the background fields to zero $A^{1,2}\to 0$ leads to the action (12)
used in the main text. The background fields are however useful to generate
various correlation functions. For example, the retarded Green’s function of
charge density is

$\begin{split}G^{R}(t,x)&\equiv i\theta(t){\rm
Tr}\left(\rho[n(t,x),n]\right)\\\ &=i\langle n_{r}(t,x)n_{a}\rangle\\\
&=i\frac{\delta^{2}\log Z}{\delta(iA_{a0}(t,x))\delta(iA_{r0}(0))}\\\
&=-i\langle
n(t,x)\partial_{i}\left(\sigma(n)\partial_{i}\phi_{a}\right)\rangle\,.\end{split}$
(58)

In particular, when $\sigma(n)=\sigma={\rm const}$, the retarded Green’s
function is simply related to the $\langle n\phi_{a}\rangle$ propagator

$G^{R}(\omega,q)=i\sigma q^{2}\langle n\phi_{a}\rangle(\omega,q)\,.$ (59)

The two-point function can be obtained from $G^{R}$ as usual from a
fluctuation-dissipation relation

$\begin{split}\langle
nn\rangle(\omega,q)&=\frac{2}{1-e^{-\beta\omega}}\operatorname{Im}G^{R}(\omega,q)\,.\end{split}$
(60)

###  A.3 1-loop calculation

The universal leading 1-loop correction to diffusion was computed in Ref.
Chen-Lin  _et al._ (2019) (see also Refs. Jain and Kovtun (2022); Jain  _et
al._ (2021); Winer and Swingle (2022); Sogabe and Yin (2022); Abbasi  _et al._
(2022); Chao and Schaefer (2023) for further studies of loop effects in the
EFT of diffusion). We review the derivation here, and discuss an interesting
cancellation between certain diagrams that simplifies the calculation.

We focus on non-analytic and UV finite corrections to $G^{R}(\omega,q)$. This
correlator also receives UV divergent corrections, which can be absorbed with
local counterterms in the EFT, and renormalize existing transport parameters.
These can be interesting in their own right Kovtun  _et al._ (2011), but do
not lead to power-law corrections to diffusive behavior, which are the focus
of this paper. We can therefore omit diagrams such as the one shown in Fig.
6A, which cannot produce novel singular IR structure. One therefore only needs
to the cubic vertices of the EFT:

$\mathcal{L}^{(3)}=i\sigma^{\prime}n(\nabla\phi_{a})^{2}+\frac{1}{2}D^{\prime}\nabla^{2}\phi_{a}n^{2}\,.$
(61)

Diagramatically, one expects these will generate 1-loop corrections
proportional to $\sigma^{\prime 2}$, $\sigma^{\prime}D^{\prime}$, and
$D^{\prime 2}$. The first is shown in Fig. 6B and can easily be seen to
vanish: indeed, the only such diagram involves a loop where all poles in
frequency poles lie on the same half complex plane. Alternatively, this can be
seen to vanish in time domain using the latest time condition (54). Indeed,
the diagram involves the computation of a correlator
$\langle(n\phi_{a})(t_{1})(n\phi_{a})(t_{2})\rangle$ – since a $\phi_{a}$
field appears at the latest time (whether it is $t_{1}$ or $t_{2}$), the
correlator must vanish. These types of considerations were already well known
to produce diagrammatic simplifications in Schwinger-Keldysh EFTs, see, e.g.,
Ref. Gao  _et al._ (2020). There is in fact a more general argument, showing
that the EFT (12) with $D(n)=\rm const$ but $\sigma(n)$ arbitrary has no loop
corrections to the two-point function; acting with the diffusive kernel on the
two point function, one has

$(\partial_{t}-D\nabla^{2})\langle
n(t,x)n\rangle=-2i\nabla\langle\left(\sigma(n)\nabla\phi_{a}\right)_{t,x}n\rangle\,,$
(62)

where we used the equation of motion $\delta\mathcal{L}/\delta\phi_{a}=0$
(this equation holds up to contact terms $\propto\delta(t)\delta(x)$. The
right-hand side vanishes because a $\phi_{a}$ field is at the latest time
(assuming without loss of generality that $t>0$), showing that the two-point
function is unaffected by nonlinearities in this theory, and is equal to

$\langle n(t,x)n\rangle=\frac{\chi}{(4\pi Dt)^{d/2}}e^{-x^{2}/(4Dt)}\,.$ (63)

This is known to occur in certain lattice gas models Spohn (2012). Here we
have temporarily ignored higher-derivative corrections, which will enter as
usual through $F_{0,n}$ as in (50).

A

B

C

D

Figure 6:  Diagrams not contributing to transport corrections to diffusion:
(a) Diagrams where the external momentum does not flow through a loop cannot
produce new IR singularities; they only renormalize tree-level transport
parameters. (b) The 1-loop contribution proportional to $\sigma^{\prime 2}$
vanishes due to the latest time condition (54). (c) The two 1-loop
contributions proportional to $\sigma^{\prime}D^{\prime}$ cancel. (d) A
similar cancellation happens at 2-loop for the
$\sigma^{\prime\prime}D^{\prime\prime}$ contribution.

We have established that the $\sigma^{\prime 2}$ contribution to $\langle
nn\rangle(\omega,q)$ vanishes. One can in fact show that the
$\sigma^{\prime}D^{\prime}$ contribution vanishes as well, although the
argument is slightly more subtle. Diagramatically, the two diagrams that give
corrections to $\langle nn\rangle(\omega,q)$ are shown in Fig. 6C. They can be
shown to cancel by explicit calculation—however, the cancellation only happens
after performing the integral over frequency, and dropping a UV divergence in
the integral over momenta. Note that, after amputating one external leg on the
$D^{\prime}$ vertex, the remaining object to be computed is a two-point
function between $n$ and the normal-ordered composite operator $n^{2}$:

$\langle n(t,x)n^{2}(0,0)\rangle\,.$ (64)

Crucially, this object is to be computed in the theory with $\sigma^{\prime}$
as its only cubic interaction (since the $D^{\prime}$ vertex has already been
used). By time-reversal symmetry, one can take $t>0$. Acting with the
diffusive kernel and using the equation of motion as in (62), one finds again
that the result vanishes, which implies that this correlator must be
proportional to the diffusive two-point function

$\langle n(t,x)n^{2}(0,0)\rangle\propto\frac{1}{t^{d/2}}e^{-x^{2}/(4Dt)}\,.$
(65)

The diagram therefore at most renormalizes $\chi$ or $D$, without producing
new non-analytic structures.

We are left with the $D^{\prime 2}$ contribution to the two-point function,
shown in Fig. 2. It is simplest to study the 1-loop correction to the retarded
Green’s function $G^{R}$, which is simply related to $\langle
n\phi_{a}\rangle(\omega,q)$ by (59) (given that we have shown that one can set
$\sigma(n)=\sigma=$ const in the action). This is done in Sec. III, where it
is shown that the loop can be expressed as a correction $D\to D+\delta
D(\omega,q)$ with (19)

$\displaystyle\delta D(\omega,q)$ $\displaystyle=-iD^{\prime
2}\int_{p^{\prime}}q^{\prime 2}\langle n\phi_{a}\rangle(p^{\prime})\langle
nn\rangle(p-p^{\prime})$ (66)
$\displaystyle=-i\chi(D^{\prime})^{2}\int\frac{d^{d}q^{\prime}}{(2\pi)^{d}}\frac{q^{\prime
2}}{\omega+iD\left[q^{\prime 2}+(q-q^{\prime})^{2}\right]}$

where in the second line we inserted the propagators (16) and evaluated the
integral over frequencies $\int\frac{d\omega}{2\pi}$. Here
$\chi\equiv\sigma/D$. Changing integration variables to $q^{\prime}\to k\equiv
q^{\prime}-\frac{1}{2}q$ and defining

$z\equiv q^{2}-\frac{2i\omega}{D}\,,$ (67)

this becomes

$\begin{split}\delta D(\omega,q)&=-\frac{\chi D^{\prime
2}}{2D}\int\frac{d^{d}k}{(2\pi)^{d}}\\\ &-i\omega\frac{\chi D^{\prime
2}}{D^{2}}\int\frac{d^{d}k}{(2\pi)^{d}}\frac{1}{z+4k^{2}}\,.\end{split}$ (68)

The term in the first line is a UV-divergent contribution to the diffusivity,
and can be absorbed with a counterterm $D\to D+\delta D$ in the EFT. The term
in the second line instead has interesting non-analytic IR behavior. It is
entirely UV finite in $d=1$. In $d\geq 2$ it produces additional UV
divergences that are analytic in $z$, and can be absorbed by higher-derivative
counterterms in the EFT. The UV finite part is given by

$\begin{split}\alpha_{d}(z)&\equiv\int\frac{d^{d}k}{(2\pi)^{d}}\frac{1}{z+4k^{2}}-{\rm
UV}\\\
&=\frac{(-z)^{\frac{d}{2}-1}}{(16\pi)^{d/2}\Gamma(\frac{d}{2})}\cdot\begin{cases}i\pi&\hbox{if
$d$ odd}\,,\\\ \log\frac{1}{z}&\hbox{if $d$ even}\,.\end{cases}\end{split}$
(69)

One therefore finds Eq. (20).

###  A.4 2-loop calculation

In systems with charge conjugation symmetry, the EFT must by invariant under

$n\to-n\,,\qquad\phi_{a}\to-\phi_{a}\,.$ (70)

This forbids cubic terms in the EFT: in particular,
$D^{\prime},\,\sigma^{\prime}=0$. The leading nonlinearities are instead
quartic, and can be found again by expanding (12):

$\mathcal{L}^{(4)}=\frac{i}{2}\sigma^{\prime\prime}n^{2}(\nabla\phi_{a})^{2}+\frac{1}{6}D^{\prime\prime}\nabla^{2}\phi_{a}n^{3}\,.$
(71)

The leading fluctuation correction to transport then comes from 2-loop
diagrams involving the two quartic vertices, such as those in Figs. 2 and 6D.
These have not been computed before – we evalute them below. The contribution
proportional to $\sigma^{\prime\prime 2}$ vanishes due to the latest time
condition (54). The contributions proportional to
$\sigma^{\prime\prime}D^{\prime\prime}$ can also be shown to vanish, following
the same argument as in the previous section: they involve computing $\langle
n(t,x)n^{3}(0,0)\rangle$ in the theory with only the quartic interaction
$\sigma^{\prime\prime}$. The only remaining contribution is the one
proportional to $D^{\prime\prime 2}$. It is studied in Sec. III and leads to a
correction (see Eq. (24))

$\displaystyle\delta D(\omega,q)$ $\displaystyle=-\frac{i}{2}D^{\prime\prime
2}\\!\\!\int_{p^{\prime},p^{\prime\prime}}\\!\\!\\!q^{\prime\prime 2}\langle
n\phi_{a}\rangle(p^{\prime\prime})\langle
nn\rangle(p^{\prime}-p^{\prime\prime})\langle nn\rangle(p-p^{\prime})$
$\displaystyle=\frac{\chi D^{\prime\prime
2}}{2D^{2}}\int_{p^{\prime}}(-i\omega^{\prime})\alpha_{d}\left(q^{\prime
2}-\frac{2i\omega^{\prime}}{D}\right)\langle nn\rangle(p-p^{\prime})$ (72)
$\displaystyle=\frac{(\chi D^{\prime\prime})^{2}}{2D^{2}}\times$
$\displaystyle\int_{q^{\prime}}\\!\\!\left(-i\omega+D(q-q^{\prime})^{2}\right)\alpha_{d}\left(q^{\prime
2}+2(q-q^{\prime})^{2}-\frac{2i\omega}{D}\right)\,,$

where in the third line we used the result of the 1-loop calculation (66), and
in the last one we performed the $\omega^{\prime}$ integral by residue (note
that $\alpha_{d}(q^{\prime 2}-\frac{2i\omega^{\prime}}{D})$ is analytic in the
upper-half $\omega^{\prime}$ plane). Changing integration variables to
$q^{\prime}\to k\equiv\frac{3}{\sqrt{2}}(q^{\prime}-\frac{2}{3}q)$ leads to

$\begin{split}\delta D&(\omega,q)=\frac{(\chi
D^{\prime\prime})^{2}}{2D^{2}}\left(\tfrac{\sqrt{2}}{3}\right)^{d}\times\\\
&\int_{k}\\!\\!\left(-i\omega+\frac{1}{9}D(q^{2}+2k^{2})\right)\alpha_{d}\left(\frac{2}{3}(k^{2}+z)\right)\,,\end{split}$
(73)

where we defined

$z\equiv q^{2}-\frac{3i\omega}{D}\,.$ (74)

Let us first focus on $d$ odd, where one can write

$\alpha_{d}(z)=a_{d}z^{\frac{d}{2}-1}\,,$ (75)

with $a_{d}=\frac{(-1)^{\frac{d-1}{2}}\pi}{(16\pi)^{d/2}\Gamma(d/2)}$. The
corrections then takes the form

$\begin{split}\delta D(\omega,q)&=\frac{3(\chi
D^{\prime\prime})^{2}}{4D^{2}}\times\\\
&\left[\left(-i\omega+\frac{D}{9}q^{2}\right)J_{d,0}(z)+\frac{2D}{9}J_{d,2}(z)\right]\,,\end{split}$
(76)

with

$J_{d,n}(z)=\left(\tfrac{\sqrt{2}}{3}\right)^{d}a_{d}\int\frac{d^{d}k}{(2\pi)^{d}}k^{n}(k^{2}+z)^{\frac{d}{2}-1}\,.$
(77)

This final integral has several power-law UV divergences: these are analytic
in $\omega,q$ and can therefore be absorbed with counterterms in higher-
derivative corrections to hydrodynamics. We thus focus on the UV finite (or UV
log-divergent), non-analytic part. By dimensional analysis, it has the form

$J_{d,n}(z)=z^{d-1+\frac{n}{2}}\left[b_{d,n}\log\frac{\Lambda^{2}}{z}+c_{d,n}\right]+\hbox{UV}\,.$
(78)

Note that the $c_{d,n}$ contribution is always analytic, so that it can be
ignored. The coefficients of interest are found to be

$b_{d,0}=\frac{1}{(12\sqrt{3}\pi)^{d}}\frac{1}{\Gamma(d)}\,,\qquad
b_{d,2}=-\frac{1}{2}b_{d,0}\,.$ (79)

This second equality implies that the $q^{2}z^{d-1}\log\frac{1}{z}$ correction
vanishes, so that $\delta D(\omega,q)$ is proportional to
$(-i\omega)z^{d-1}\log\frac{1}{z}$ – this guarantees that static correlators
$\omega\to 0$ are analytic as expected Jain  _et al._ (2021). We are then left
with

$\delta D(\omega,q)=\frac{1}{2}b_{d,0}\frac{(\chi
TD^{\prime\prime})^{2}}{D^{2}}(-i\omega)z^{d-1}\log\frac{1}{z}\,.$ (80)

One can repeat this calculation in $d$ even, where instead of (75) one has
$\alpha_{d}(z)=\tilde{a}_{d}z^{\frac{d}{2}-1}\log\frac{1}{z}$. One finds the
same results as above, apart for a sign $b_{d,n}\to-b_{d,n}$. The result in
general dimension therefore takes the form (25).

###  A.5 Fourier transformation of corrections to diffusion

In this section, we detail the computation of the inverse Fourier transform
that produces the two-point function, focusing on $d=1$ for simplicity

$\langle n(t,x)n\rangle\equiv\int\frac{d\omega dq}{(2\pi)^{2}}e^{-i\omega
t+iqx}\langle nn\rangle(\omega,q)\,,$ (81)

where $\langle nn\rangle(\omega,q)$ can be obtained from the retarded Green’s
function using (60)

$\begin{split}\langle
nn\rangle(\omega,q)&\simeq\frac{2}{\omega}\operatorname{Im}\frac{\sigma
q^{2}}{-i\omega+(D+\delta D(\omega,q))q^{2}}\\\ &\simeq\frac{2\sigma
q^{2}}{\omega^{2}+D^{2}q^{4}}-\frac{2}{\omega}\operatorname{Im}\frac{\sigma
q^{4}\delta D(\omega,q)}{(-i\omega+Dq^{2})^{2}}\,,\end{split}$ (82)

where we have expanded to linear order in $\delta D$ (given by (22) or (27))
because we are only interested in the leading correction.

The Fourier transformation of the first term is straightforward, and given by

$\langle n(t,x)n\rangle_{0}=\frac{\chi}{\sqrt{4\pi Dt}}e^{-y^{2}/4}\,,\qquad
y\equiv\frac{x}{\sqrt{Dt}}\,.$ (83)

Let us now turn to the Fourier transform of the correction. Introducing the
dimensionless variables $w=\omega/(Dq^{2})$ and $\tau=tDq^{2}$, one has

$\begin{split}\delta\langle
nn\rangle(t,q)&=\int\frac{d\omega}{2\pi}e^{-i\omega t}\langle
nn\rangle(\omega,q)\\\ &={\chi
q^{2}}\int\frac{dw}{2\pi}e^{-iw\tau}\frac{\delta
D/(-i\omega)}{(1-iw)^{2}}\,.\end{split}$ (84)

We have taken $t>0$, implying that one can disregard non-analyticities in the
upper-half $w$ plane. Inserting the 1-loop expression (22) leads to

$\delta\langle nn\rangle(t,q)=\frac{i\chi^{2}D^{\prime
2}|q|e^{-\tau/2}}{4\sqrt{2}D^{2}}\\!\int_{-\infty}^{0}\\!\frac{dz}{2\pi}{\rm
Disc}\frac{e^{z\tau}}{\sqrt{z}(z+\frac{1}{2})^{2}}\,,$

where we defined $z=\frac{1}{2}-iw$ and deformed the contour to pick up the
discontinuity across the two-diffuson branch cut, ${\rm
Disc}f(z)=f(z+i0^{+})-f(z-i0^{+})$. Evaluating the integral gives

$\delta\langle nn\rangle(t,q)=\frac{\chi^{2}D^{\prime
2}|q|e^{-\tau/2}}{4D^{2}}\left[e^{-\tau/2}(1+\tau){\rm
Erfi}\\!\left(\tfrac{\sqrt{\tau}}{\sqrt{2}}\right)-\frac{\sqrt{2\tau}}{\sqrt{\pi}}\right]$
(85)

with ${\rm Erfi}(z)\equiv{\rm Erf}(iz)/i$. To perform the final Fourier
transform $\int\frac{dq}{2\pi}e^{iqx}$, one can express the integrand as a
product of two Fourier transforms, and evaluate their convolution. This gives
Eq. (1).

The 2-loop contribution can be obtained similarly: one inserts (27) in (84),
and evaluates the discontinuity across the 3-diffuson branch cut to obtain

$\delta\langle nn\rangle(t,q)=\frac{\chi^{3}D^{\prime\prime
2}q^{2}}{24\sqrt{3}\pi D^{2}}\left[\tau
e^{-\tau}\left(\log\tfrac{1}{q^{2}}+{\rm
Ei}(\tfrac{2\tau}{3})\right)-\frac{3}{2}e^{-\frac{\tau}{3}}\right]$

where ${\rm Ei}(z)\equiv-\int_{-z}^{\infty}\frac{du}{u}e^{-u}$
(ExpIntegralEi[z] in Mathematica). We were not able to express the final
Fourier transform $\int\frac{dq}{2\pi}e^{iqx}$ in terms of known special
functions – the resulting integral is shown in Eq. (6). The 1-loop and 2-loop
scaling functions are shown in Fig. 7.

Figure 7:  Universal scaling functions describing non-linear corrections to
the diffusive structure factor as a function of the hydrodynamic variable
$y=x/\sqrt{Dt}$. The leading order Gaussian spreading is
$F_{0,0}=e^{-y^{2}/4}$ and the leading order non-linear fluctuations are given
by Eq. (2) for general systems and by Eq. (6) in the presence of particle-hole
symmetry.

###  A.6 EFT predictions for the domain wall quench

In this section we present the EFT predictions for the domain-wall initial
condition, Eq. (42). As explained in the main text, the conversion from the
dynamical structure factor to the domain wall picture is a simple integration
in space, Eq. (44). Here we present the equations employed in the fitting
processes,

$\displaystyle\textbf{F}_{0,0}=\frac{1}{2}\text{Erf}(y/2)$
$\displaystyle\textbf{F}_{1,0}=\frac{\chi D^{\prime 2}}{\sqrt{4\pi
D^{5}}}\left(-\frac{4e^{-\frac{y^{2}}{2}}y}{8\sqrt{\pi}}-\frac{e^{-\frac{y^{2}}{4}}(-6+y^{2})\text{Erf}(y/2)}{16}\right)$
$\displaystyle\textbf{F}_{0,1}=e^{-\frac{y^{2}}{4}}y(c_{1}+c_{2}y^{2})$
$\displaystyle\textbf{F}_{2,0}=\frac{\chi^{2}D^{\prime\prime
2}}{24\pi\sqrt{3}D^{3}}\left(\tilde{\textbf{F}}_{2,0}+\textbf{F}_{2,0}^{\prime}\log
t\right)$ (86)
$\displaystyle\tilde{\textbf{F}}_{2,0}=\int_{0}^{\infty}\frac{ds}{\pi}\sin(sy)s$
$\displaystyle\\!\\!\times\left[s^{2}e^{-s^{2}}\left(\log\frac{1}{s^{2}}+{\rm
Ei}\left(\tfrac{2s^{2}}{3}\right)\right)-\frac{3}{2}e^{-s^{2}/3}\right]$
$\displaystyle\textbf{F}_{2,0}^{\prime}=-\frac{e^{-\frac{y^{2}}{4}}y(y^{2}-6)}{16\sqrt{\pi}}.$

The parameters in the term $\textbf{F}_{0,1}$ are defined in order to absorb
the numerical constant.

Equation (49) follows from a simple manipulation of the three-point function
derived in Delacretaz and Mishra (2023),

$\begin{split}&f(\bar{x},x,t)=\langle
n(x_{1},t)n(x_{2},t)n(0,0)\rangle=\frac{\chi^{2}D^{\prime}}{8\pi
D^{2}t}\times\\\
&e^{-\frac{y^{2}_{1}+y^{2}_{2}}{4}}\left(1-\frac{\sqrt{\pi}y_{1}}{2}e^{\frac{y^{2}_{2}}{4}}\left(\text{Erf}(y_{2}/2)+\text{sign}(y_{1}-y_{2})\right)\right)+\\\
&y_{2}\leftrightarrow y_{1},\end{split}$ (87)

where $y_{1}=x_{1}/\sqrt{Dt}$, $y_{2}=x_{2}/\sqrt{Dt}$, and
$\bar{x}=(x_{1}+x_{2})/2$, $x=x_{1}-x_{2}$ correspond to the centre of mass
coordinates. In the domain-wall initial state the centre of mass will be
summed over all positions, while the second operator is fixed at $x_{2}=0$.
This corresponds to a spatial integration with respect to the center of mass
coordinate,

$s_{3,\text{EFT}}=\frac{1}{\chi}\int^{0}_{\infty}f(\bar{x},x,t)d\bar{x}.$ (88)

##  Appendix B EFT for multiple diffusing densities

###  B.1 General construction of the EFT

The EFT approach can be generalized to account for multiple continuity
relations giving rise to conserved densities. Consider a thermalizing system
with abelian $U(1)^{N}$ symmetry and couple it to background fields $A_{\mu
I}^{A}$, where $I=1,2$ denotes the SK contour, and the index $A=1,\ldots,N$.
The action will be made up of the gauge invariant combinations

$B_{\mu I}^{A}\equiv A_{\mu I}^{A}+\partial_{\mu}\phi_{I}^{A}\,.$ (89)

Imposing diagonal shift symmetry

$B_{ir}^{A}\to B_{ir}^{A}+\partial_{i}\lambda^{A}(\vec{x})\,,$ (90)

and KMS as before, the quadratic action to leading order in derivatives is

$\mathcal{L}^{(2)}=B_{0a}\cdot\chi\cdot
B_{0r}+B_{ia}\cdot\sigma\cdot\left(iTB_{ia}-\dot{B}_{ir}\right)\,,$ (91)

where $\chi_{AB},\,\sigma_{AB}$ are matrices, and dots denote matrix
multiplication. Both $\chi$ and $\sigma$ have to be symmetric by time-reversal
symmetry (Onsager relation). The cubic action to leading order in derivatives
is

$\mathcal{L}^{(3)}=\frac{1}{2}B_{0a}\cdot\partial_{A}\chi\cdot
B_{0r}B_{0r}^{A}+B_{ia}\cdot\partial_{A}\sigma\cdot\left(iTB_{ia}-\dot{B}_{ir}\right)B_{0r}^{A}\,.$
(92)

The cubic interactions arise as before as dependence of transport or
thermodynamic parameters on potentials: $\partial_{A}\chi\equiv
d\chi/d\mu^{A},\,\partial_{A}\sigma\equiv d\sigma/d\mu^{A}$. One can change
variables to the density as before $\phi_{r}^{A}\to
n^{A}\equiv\frac{\delta\mathcal{L}}{\delta A_{a0}^{A}}$. In terms of these
variables, the full action up to cubic order is

$\displaystyle\mathcal{L}^{(2)}=$ $\displaystyle\,B_{0a}\cdot
n+B_{ia}\cdot\sigma\cdot(iTB_{ia}-F_{0i,r})-B_{ia}\cdot
D\cdot\partial_{i}n\,,$ $\displaystyle\mathcal{L}^{(3)}=$
$\displaystyle\,-\frac{1}{2}B_{ia}\cdot\partial^{A}D\cdot\partial_{i}\left(nn_{A}\right)$
$\displaystyle\,+B_{ia}\cdot\partial^{A}\sigma\cdot(iTB_{ia}-F_{0i,r})n_{A}$
(93)
$\displaystyle\,+\frac{1}{2}B_{ia}\cdot\partial^{A}\sigma\cdot\chi^{-1}\cdot\left(n\partial_{i}n_{A}-\partial_{i}nn_{A}\right)\,.$

The derivatives of transport parameters are now taken with respect to
densities, e.g.

$\partial^{A}\sigma=\frac{d}{dn_{A}}\sigma=(\chi^{-1})^{AB}\partial_{B}\sigma\,,$
(94)

and the diffusion matrix has been defined as

$D\cdot\chi=\sigma\,,\qquad\hbox{or}\qquad
D_{A}{}^{B}=\sigma_{AC}(\chi^{-1})^{CB}\,.$ (95)

As the product of two symmetric positive matrices, $D$ can be diagonalized –
its eigenvalues correspond to the location of poles of the density two-point
function $\omega=-iD_{A}q^{2}$.

While the first two lines in $\mathcal{L}^{(3)}$ (B.1) are analogous to the
$N=1$ case, the term in the third line is qualitatively new. It is a
contribution to $j_{ir}=\delta\mathcal{L}/\delta A_{ia}$, as expected (see Eq.
(31)). These lead to different corrections to diffusion. We study them below,
in the more constrained situation of non-abelian densities.

###  B.2 Non-abelian densities

The EFT can also be straightforwardly generalized to non-abelian Lie groups
Glorioso  _et al._ (2021); Akyuz  _et al._ (2023). We will focus here on
$SU(2)$ for concreteness. Because the densities $n^{A}$ transform linearly (in
the adjoint representation) under group action, one can implement the symmetry
by making sure that they are contracted with group covariant tensors (for
$SU(2)$, these are $\delta_{AB},\,\epsilon_{ABC}$). The nonlinear
transformation of $\phi_{a}^{A}$ requires more attention: instead of using
$B_{\mu a}=\partial_{\mu}\phi_{a}+A_{\mu a}$, one should use the Maurer-Cartan
form

$iB_{\mu
a}=e^{-i\phi_{a}^{A}T_{A}}\left(\partial_{\mu}+iA_{\mu}^{B}T_{B}\right)e^{i\phi_{a}^{C}T_{C}}$
(96)

which also transforms in the adjoint of $SU(2)$. Here $T_{A}$ are the
generators of the algebra. The cubic $SU(2)$ invariant action is therefore

$\displaystyle\mathcal{L}^{(2)}=$ $\displaystyle\,B_{0a}\cdot n+\sigma
B_{ia}\cdot(iTB_{ia}-F_{0i,r})-DB_{ia}\cdot\partial_{i}n\,,$
$\displaystyle\mathcal{L}^{(3)}=$
$\displaystyle\,\lambda\epsilon_{ABC}B_{ia}^{A}n^{B}\left(\chi
F_{0i}^{C}+\partial_{i}n^{C}\right)\,.$ (97)

The term in the last line of (B.1) produced the cubic term above: we have
written its coefficient as
$\partial_{A}\sigma_{BC}\equiv\chi\lambda\epsilon_{ABC}$. Note that there are
also cubic terms in the first line, coming from expanding (96): however, one
can show using the leading order equation of motion that they will not
contribute to the 1-loop correction studied below.

###  B.3 1-loop correction

Let us study the 1-loop correction to $\langle n_{r}n_{a}\rangle$. There are
two contributions: the first comes from the nonlinear piece in

$n_{a}^{A}\equiv\frac{\partial\mathcal{L}}{\partial
A_{0r}}=-\sigma\nabla^{2}\phi_{a}^{A}+\chi\lambda\epsilon^{ABC}\partial_{i}\left(\partial_{i}\phi_{a}^{B}n^{C}\right)+\cdots\,.$
(98)

Writing $\langle n_{r}^{A}n_{a}^{B}\rangle=\delta^{AB}\langle
n_{r}n_{a}\rangle$, this leads to the following correction to $\langle
n_{r}n_{a}\rangle$:

$\begin{split}2i\chi\lambda^{2}\langle
n\phi\rangle(p)q_{i}q_{j}\int_{p^{\prime}}\left[q^{\prime}_{j}(2q^{\prime}-q)_{i}\right]\langle
n\phi\rangle(p^{\prime})\langle nn\rangle(p-p^{\prime})\,.\end{split}$ (99)

The other contribution comes from two insertions of the cubic interaction
$S^{(3)}$:

$\begin{split}&2\sigma\lambda^{2}[\langle
n\phi\rangle(p)]^{2}q^{2}q_{i}\times\\\
&\int_{p^{\prime}}\left[q^{\prime}_{j}(2q^{\prime}-q)_{i}(2q-q^{\prime})_{j}\right]\langle
nn\rangle(p-p^{\prime})\langle n\phi_{a}\rangle(p^{\prime})\,.\end{split}$
(100)

Summing these two diagrams one finds (in $d=1$ dimensions)

$\displaystyle G^{R}_{nn}(\omega,q)=i\langle n_{r}n_{a}\rangle(\omega,q)$
(101) $\displaystyle\simeq\frac{\sigma
q^{2}}{-i\omega+Dq^{2}}+\frac{\lambda^{2}T\chi^{2}}{D}q^{2}(i\omega)\frac{\sqrt{q^{2}-\frac{2i\omega}{D}}}{\left(Dq^{2}-i\omega\right)^{2}}+\cdots\,.$

Like before, there are several consistency checks that this result satisfies:
it vanishes when $q\to 0$ (as it must by current conservation), and when
$\omega\to 0$ (as it must by analyticity of static correlators). The two
diagrams above do not satisfy the latter check individually, only their sum
does. Our result slightly differs from one obtained in a strong noise
expansion in Claeys  _et al._ (2022) – because their result does not become
analytic in the static limit, we suspect that they may have missed a
contribution to the 1-loop correction.

Figure 8:  Estimation of stability of corrections for different fillings
$\mu=0,0.2,0.3$ and bond dimensions $d=$ $256$-red, $400$-blue, $600$-green,
for the staggered XXZ chain. (A): Profiles of magnetizations at 3 different
times $t=100,200,400$ (three visibly different sets of curves). Bond dimension
$d=600$ is only shown for $t=100$ the maximum integrated time is $t=150$. A
slight difference is only visible for $\mu=0.3$. (B): Difference between the
profiles at times $t=100,200,400$, denoted by full,dashed and dotted lines
respectively. For $\mu=0.3$ the difference is almost an order of magnitude
larger than $\mu=0,0.2$, indicating the enhancement of simulation error at
larger fillings. (C): Predictions of diffusivity by employing the fitting
method (I). At $\mu=0$ diffusivity decreases with time indicating that the
system is not yet at asymptotic diffusion. At larger $\mu$ we observe a
monotonic increase over time, compatible with the EFT predictions. (D):
Correction to diffusion, $\Delta s=s-\textbf{F}_{0,0}$ using method II and
times up to $T=150$ for $d=256,600$. At $\mu=0$ the difference between this
plot and the main text is due to the smaller timescale used here in order to
compare with $d=600$. For $\mu=0.2$, the profiles are more than one order of
magnitude larger than what the error is estimated from (B), suggesting that
the quantitative structure of the correction is not considerably affected by
errors. For $\mu=0.3$ the profiles agree well, however deviations in (B),
become significant for $t\geq 200$, and therefore, we don’t fit these data in
the main text.

###  B.4 Fourier transformation

We’d like to compute the fourier transform of

$\displaystyle\delta\langle nn\rangle(\omega,q)$
$\displaystyle=\frac{2T}{\omega}\operatorname{Im}\delta G^{R}_{nn}(\omega,q)$
(102)
$\displaystyle=\frac{\lambda^{2}(T\chi)^{2}}{D}q^{2}\frac{\sqrt{q^{2}-\frac{2i\omega}{D}}}{\left(Dq^{2}-i\omega\right)^{2}}+\hbox{c.c.}\,.$

Fourier transforming $\int\frac{d\omega}{2\pi}e^{-i\omega t}$ by picking up
the cut as usual, one finds

$\begin{split}\delta\langle
nn\rangle(t,q)&=-\frac{\lambda^{2}(T\chi)^{2}}{D^{2}}\frac{2}{\sqrt{\pi}}|q|e^{-\tau/2}\times\\\
&\left[-\sqrt{\tau/2}+(\tau-1)F\left(\sqrt{\tau/2}\right)\right]\,.\end{split}$
(103)

The final Fourier transform $\int\frac{dq}{2\pi}e^{iqx}$ can be performed by
convoluting the Fourier transforms of both products above. The result is shown
in Eq. (35).

##  Appendix C Details on Quantum Transport

In this section we (1) identify the effects of the truncation of the bond
dimension in the DMPO dynamics and (2) compare different fitting approaches
for the extraction of diffusivity. We find that even when truncation is only
weakly affecting the simulation, the fact the system is not yet at the
asymptotic regime can lead to different fitting results depending on the
method.

###  C.1 Effects of information truncation in the dynamics

The main source of error in tensor network approaches to quantum dynamics is
the truncation of information, e.g. operator entanglement in the case of
density matrix evolution Zanardi (2001). Operator entanglement is a generated
by quantum dynamics and is expected to increase linearly with time
$S_{op}\propto t$, corresponding to an exponential scaling of the required
classical resources (bond dimension of local tensors) Hartmann  _et al._
(2009); Žnidarič  _et al._ (2008). In practice, one sets a maximum dimension
for the local tensors effectively bounding the amount of operator entanglement
in the state. The effects of this truncation to the long-time dynamics in
linear-response quenches is an active area of research White  _et al._ (2018);
von Keyserlingk  _et al._ (2022). Since there is no theory for the effects of
truncation to the dynamics we simply change the bond dimension and compare the
results. If there is agreement between bond dimensions we assume that the
truncation is weak.

In the presence of dephasing, $\gamma$, in the system, coherences are
destroyed at a timescale $t_{\gamma}\sim 1/\gamma$. This leads to a saturation
of $S_{op}$. Therefore, if the bond dimension is enough to produce accurate
results up to $t_{\gamma}$, it will also be accurate for $t\gg t_{\gamma}$. We
confirm this for the results shown in the main text ($\gamma=0.1$) by
comparing bond dimensions $d=128,256$ (not shown).

Coherent simulations on the other hand are much more demanding as the amount
of resources increases exponentially with time. In Figure 8.A we show that
different bond dimensions agree well at different times at the leading order
level. To get a better understanding of the accuracy we calculate the
distribution differences between different bond dimensions, Figure 8.B. The
values of these differences will be employed to estimate the accuracy of sub-
leading effects. These results show that as we deviate from half-filling the
simulations become more demanding, leading to a decrease of accuracy for the
same resources. While we don’t fully understand this phenomenon, it’s likely
related to the increased strength of the non-linear corrections.

###  C.2 Fitting methods

Following the raw data comparison, we perform a consistency check between
different fitting approaches to the simulation data. The classical numerics
presented in Ref. III.5 suggest that when dealing with limited resources,
making use of our knowledge of the general structure of corrections to
diffusion improves the precision in fitting transport parameters.
Consequently, we employ three qualitatively different fitting methods that are
designed to take an increasing amount of corrections to diffusion into
account. The first two methods (I,II) perform fits on the dynamical structure
factor. The third method (III) is based on Fick’s first law and the total
current in the system.

Method I assumes knowledge of only the leading order diffusion,

$\mathcal{S}_{I}(x,t)=\mathbf{F}_{0,0}=\frac{1}{2}\text{Erf}(x/\sqrt{4\alpha_{\overline{\sigma},t}t}),$
(104)

where an explicit time dependence on diffusivity
$\alpha_{\overline{\sigma},t}$ is assumed in order to encapsulate the finite
time corrections. This method is simple to implement and commonly used.
However it is entirely phenomenological, and, strictly speaking, incorrect –
indeed, the scaling corrections discussed in Sec. III imply that the
autocorrelation function does not take the form Eq. (104) at intermediate
times. We nevertheless study this method to compare it to other methods that
are consistent with EFT predictions. The fit is performed for each time $t$
and equilibrium magnetization $\overline{\sigma}$ on the spatial profile of
magnetization $s(x)$ defined by Eq. (43). One can then estimate diffusivity
from the longest simulation time $t=T$. However, the true diffusivity is
defined as $D=\lim_{t\rightarrow\infty}\alpha_{\overline{\sigma},t}$. In
principle, it is possible to extrapolate diffusivity with some appropriate
function of inverse time. However, in simulations which have increasing errors
with time, such extrapolations can capture artifacts, so we will not be
performing them. Once diffusivity $D(\overline{\sigma})$ is extracted, the
non-linear corrections are either calculated by discrete derivatives, or
analytically, by first performing a low order polynomial fit,

$D(\overline{\sigma})=\sum^{M}_{i=0}b_{i}(\overline{\sigma})^{2i},$ (105)

where only even powers appear due to the particle-hole symmetry in the system.
In Figure 8.C we show how this method performs on the staggered field
simulations for different bond-dimensions. Compared to methods II and III, I
gives a time-dependent illustration of the systems behaviour, by showing how
the asymptotic limit is approached.

Method II, which is explained in the main text is based on a full fit of the
hypersurface $s(\overline{\sigma},x,t)$. This is done by either employing the
leading order diffusion or more elaborate assumption for the fitting function,
Eq. (44) based on our knowledge of the leading corrections to diffusion. We
use the same trial function for diffusivity as before, Eq. (105) and similar
trial functions for the linear corrections,
$c_{1/2}=\sum^{M}_{i=0}b^{1/2}_{i}(\overline{\sigma})^{2i}$. We note that
$c_{1/2}$ are taken constant in the coherent simulations since the maximum
reliable chemical potential is $\approx 0.25$ and the dependence of the
constants with chemical potential is weak. Figure 8.D shows the correction to
the leading order diffusive profile at $T=150$. We observe that for
sufficiently low $\mu\leq 0.25$ the correction is more than one order of
magnitude stronger than the bond dimension difference (Figure 8.B), indicative
of a quantitatively accurate result. For the incoherent simulation where we
fit $\mu=0-1$ it is important to allow for the constants to depend on
$\bar{\sigma}$.

The third method (III) studies the relaxation of the total current $J=\int
dx\,j(x)$. On general grounds, the constitutive relation for the current
density is

$j=-D(n)\partial_{x}n+\hbox{higher derivatives}.$ (106)

We will ignore the higher derivative terms for now, and come back to them
below. The first term is responsible for all higher-loop corrections to the
dynamic structure factor at leading order in gradients, $\textbf{F}_{\ell,0}$
(the first two $\ell=1,2$, coming from $D^{\prime},D^{\prime\prime}$, were
computed in Sec. III). All of these corrections do not contribute to
correlators of the total current: defining $C(n)$ such that
$C^{\prime}(n)=D(n)$ and integrating from $x=1$ to $x=L$, we have

$\begin{split}J&\simeq C(n(L))-C(n(1))\\\ &\simeq
D(\bar{\sigma})\left(s(L)-s(1)\right),\end{split}$ (107)

where in the second line we expanded the density around the equilibrium
magnetization $n(x)=\bar{\sigma}+s(x)$, and $D(\bar{\sigma})$ is approximated
by Eq. (105). It is justified to drop higher order terms in $s(x)$ (which
could otherwise lead to fluctuating corrections), because the dynamics has not
affected the magnetization sufficiently far from the domain wall position at
any time. This is true since for all times integrated, the Lieb-Robinson’s
light-cone has not reached the system’s boundary, meaning that
$s(L)=-s(1)\simeq\chi(\bar{\sigma})\delta$. Let us now turn to the higher-
derivative terms in Eq. (106). Because of the argument above,
$\partial^{n}_{x}s(x)=0$ at the boundaries $x=1,L$ so that other contributions
to the current that are total derivatives vanish. This includes linear higher
derivative terms $j\supset\partial_{x}^{2n+1}s$, which otherwise would have
lead to corrections of the form $\textbf{F}_{0,n}$. This also includes certain
nonlinear higher derivative terms: for example, all terms involving two
diffusive fluctuations are total derivatives
$s\partial_{x}^{2n+1}s=\partial_{x}(\ldots)$, so that 1-loop corrections with
any number of derivatives $\textbf{F}_{1,n}$ do not contribute. The first
correction to (107) comes from the leading EFT operator that is parity-odd,
and not a total derivative. For the case of a single diffusive density, this
is $j\supset s^{2}\partial_{x}^{3}s$. This will lead to a $\ell=2$ loop
correction at order $n=2$ in derivatives ($\textbf{F}_{2,2}$) to the current
decay $J(t)$ scaling as $1/t^{n+\ell/2}=1/t^{3}$ at late times. While a
precise exponent is difficult to extract from the numerics, we observe fast
decay of the total current, consistent with the observation that many of the
leading EFT corrections vanish in this observable (fast polynomial decay due
to high order hydrodynamic tails is known to be difficult to observe
quantitatively Glorioso  _et al._ (2021)). Because the convergence to the late
time value of the current $\lim_{t\to\infty}J(t)$ is therefore fast, there is
no need to extrapolate in time to obtain accurate results. We note that this
method becomes less powerful in the presence of multiple conserved charges,
since there is a larger number of terms which are not total derivatives, see
Eq. (31). In this case, the decay of the current is due to a 1-loop correction
$\textbf{F}_{1,0}$, scaling as $1/t^{n+\ell/2}=1/\sqrt{t}$.

Figure 9:  Sub-leading correction estimation for diffusivity fits using
different fitting methods for $\mu=0.2$. Faster increasing diffusivity results
to suppressed edge corrections, while slower increasing diffusivity leads to
suppressed corrections around $y=0$. The analytical result (red dashed line)
denotes the EFT result F1,0, which relies only on the fitted diffusivity
$D(\overline{\sigma})$.

The diffusivity $D(\bar{\sigma})$ is therefore obtained from the late time
current using (107), and fit as a function of magnetization using (105),
limiting ourselves to $M=2$. The parameters $b_{0},b_{1},b_{2}$ are estimated
by minimizing the distance between the measured currents at different fillings
and Eq. (106). To avoid over-parametrization, we fit a number of equilibrium
magnetization much larger than the number of free parameters (3 in this case).

The sub-leading correction is qualitatively similar for fitting methods II and
III despite the small variations in diffusivity, Figure 9, however method I
clearly underestimates diffusivity, leading to larger deviation in the
corrections.

## References

  * Gurvitch and Fiory (1987) M. Gurvitch and A. T. Fiory, “Resistivity of ${\mathrm{la}}_{1.825}$${\mathrm{sr}}_{0.175}$${\mathrm{cuo}}_{4}$ and ${\mathrm{yba}}_{2}$${\mathrm{cu}}_{3}$${\mathrm{o}}_{7}$ to 1100 k: Absence of saturation and its implications,” Phys. Rev. Lett. 59, 1337–1340 (1987).
  * Chien  _et al._ (1991) T. R. Chien, Z. Z. Wang, and N. P. Ong, “Effect of zn impurities on the normal-state hall angle in single-crystal ${\mathrm{yba}}_{2}$${\mathrm{cu}}_{3\mathrm{-}\mathit{x}}$${\mathrm{zn}}_{\mathit{x}}$${\mathrm{o}}_{7\mathrm{-}\mathrm{\delta}}$,” Phys. Rev. Lett. 67, 2088–2091 (1991).
  * Hayes  _et al._ (2016) Ian M Hayes, Ross D McDonald, Nicholas P Breznay, Toni Helm, Philip JW Moll, Mark Wartenbe, Arkady Shekhter, and James G Analytis, “Scaling between magnetic field and temperature in the high-temperature superconductor bafe2 (as1- x p x) 2,” Nature Physics 12, 916–919 (2016).
  * Schneider  _et al._ (2012) Ulrich Schneider, Lucia Hackermüller, Jens Philipp Ronzheimer, Sebastian Will, Simon Braun, Thorsten Best, Immanuel Bloch, Eugene Demler, Stephan Mandt, David Rasch, and Achim Rosch, “Fermionic transport and out-of-equilibrium dynamics in a homogeneous hubbard model with ultracold atoms,” Nature Physics 8, 213–218 (2012).
  * Schreiber  _et al._ (2015) Michael Schreiber, Sean S Hodgman, Pranjal Bordia, Henrik P Lüschen, Mark H Fischer, Ronen Vosk, Ehud Altman, Ulrich Schneider, and Immanuel Bloch, “Observation of many-body localization of interacting fermions in a quasirandom optical lattice,” Science 349, 842–845 (2015).
  * Wei  _et al._ (2022) David Wei, Antonio Rubio-Abadal, Bingtian Ye, Francisco Machado, Jack Kemp, Kritsana Srakaew, Simon Hollerith, Jun Rui, Sarang Gopalakrishnan, Norman Y Yao, _et al._ , “Quantum gas microscopy of kardar-parisi-zhang superdiffusion,” Science 376, 716–720 (2022).
  * Morvan  _et al._ (2022) Alexis Morvan, TI Andersen, Xiao Mi, Charles Neill, Andre Petukhov, Kostyantyn Kechedzhi, DA Abanin, Alexios Michailidis, R Acharya, F Arute, _et al._ , “Formation of robust bound states of interacting microwave photons,” Nature 612, 240–245 (2022).
  * Mi  _et al._ (2023) X Mi, AA Michailidis, S Shabani, KC Miao, PV Klimov, J Lloyd, E Rosenberg, R Acharya, I Aleiner, TI Andersen, _et al._ , “Stable quantum-correlated many body states via engineered dissipation,” arXiv preprint arXiv:2304.13878 (2023).
  * Kadanoff and Martin (1963) Leo P Kadanoff and Paul C Martin, “Hydrodynamic equations and correlation functions,” Annals of Physics 24, 419–469 (1963).
  * Martin  _et al._ (1973) P. C. Martin, E. D. Siggia, and H. A. Rose, “Statistical Dynamics of Classical Systems,” Phys. Rev. A 8, 423–437 (1973).
  * An  _et al._ (2021) Xin An, Gök çe Başar, Mikhail Stephanov, and Ho-Ung Yee, “Evolution of non-gaussian hydrodynamic fluctuations,” Phys. Rev. Lett. 127, 072301 (2021).
  * Bertini  _et al._ (2015) Lorenzo Bertini, Alberto De Sole, Davide Gabrielli, Giovanni Jona-Lasinio, and Claudio Landim, “Macroscopic fluctuation theory,” Rev. Mod. Phys. 87, 593–636 (2015), arXiv:1404.6466 [cond-mat.stat-mech] .
  * Crossley  _et al._ (2017) Michael Crossley, Paolo Glorioso, and Hong Liu, “Effective field theory of dissipative fluids,” JHEP 09, 095 (2017), arXiv:1511.03646 [hep-th] .
  * De Nardis  _et al._ (2018) Jacopo De Nardis, Denis Bernard, and Benjamin Doyon, “Hydrodynamic diffusion in integrable systems,” Phys. Rev. Lett. 121, 160603 (2018).
  * Nardis  _et al._ (2019) Jacopo De Nardis, Denis Bernard, and Benjamin Doyon, “Diffusion in generalized hydrodynamics and quasiparticle scattering,” SciPost Phys. 6, 049 (2019).
  * Gopalakrishnan  _et al._ (2018) Sarang Gopalakrishnan, David A. Huse, Vedika Khemani, and Romain Vasseur, “Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems,” Phys. Rev. B 98, 220303 (2018).
  * Sachdev (2007) Subir Sachdev, “Quantum phase transitions,” in  _Handbook of Magnetism and Advanced Magnetic Materials_ (Wiley Online Library, 2007).
  * Zaanen (2004) Jan Zaanen, “Why the temperature is high,” Nature 430, 512–513 (2004).
  * Bruin  _et al._ (2013) J. A. N. Bruin, H. Sakai, R. S. Perry, and A. P. Mackenzie, “Similarity of scattering rates in metals showing ¡i¿t¡/i¿-linear resistivity,” Science 339, 804–807 (2013), https://www.science.org/doi/pdf/10.1126/science.1227612 .
  * Hartnoll and Mackenzie (2022) Sean A. Hartnoll and Andrew P. Mackenzie, “Colloquium: Planckian dissipation in metals,” Rev. Mod. Phys. 94, 041002 (2022), arXiv:2107.07802 [cond-mat.str-el] .
  * Mukerjee  _et al._ (2006) Subroto Mukerjee, Vadim Oganesyan, and David Huse, “Statistical theory of transport by strongly interacting lattice fermions,” Phys. Rev. B 73, 035113 (2006), arXiv:cond-mat/0503177 .
  * Lux  _et al._ (2014) Jonathan Lux, Jan Müller, Aditi Mitra, and Achim Rosch, “Hydrodynamic long-time tails after a quantum quench,” Phys. Rev. A 89, 053608 (2014).
  * Chen-Lin  _et al._ (2019) Xinyi Chen-Lin, Luca V. Delacrétaz, and Sean A. Hartnoll, “Theory of diffusive fluctuations,” Phys. Rev. Lett. 122, 091602 (2019), arXiv:1811.12540 [hep-th] .
  * Glorioso  _et al._ (2021) Paolo Glorioso, Luca V. Delacrétaz, Xiao Chen, Rahul M. Nandkishore, and Andrew Lucas, “Hydrodynamics in lattice models with continuous non-Abelian symmetries,” SciPost Phys. 10, 015 (2021), arXiv:2007.13753 [cond-mat.stat-mech] .
  * Delacretaz and Mishra (2023) Luca V. Delacretaz and Ruchira Mishra, “Nonlinear Response in Diffusive Systems,” (2023), arXiv:2304.03236 [cond-mat.str-el] .
  * Nicoletti and Cavalleri (2016) Daniele Nicoletti and Andrea Cavalleri, “Nonlinear light–matter interaction at terahertz frequencies,” Adv. Opt. Photon. 8, 401–464 (2016), arXiv:1608.05611 [cond-mat.str-el] .
  * Chaudhuri  _et al._ (2022) Dipanjan Chaudhuri, David Barbalas, Ralph Romero III, Fahad Mahmood, Jiahao Liang, John Jesudasan, Pratap Raychaudhuri, and N. P. Armitage, “Anomalous high-temperature THz nonlinearity in superconductors near the metal-insulator transition,” (2022), arXiv:2204.04203 [cond-mat.supr-con] .
  * Schweigler  _et al._ (2017) Thomas Schweigler, Valentin Kasper, Sebastian Erne, Igor Mazets, Bernhard Rauer, Federica Cataldini, Tim Langen, Thomas Gasenzer, Jürgen Berges, and Jörg Schmiedmayer, “Experimental characterization of a quantum many-body system via higher-order correlations,” Nature 545, 323–326 (2017).