Corrections to diffusion in interacting quantum systems

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 ⟨ n ( x, t ) n ⟩ , 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 ther-malizing 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 [1] of high-T c superconductors and heavy fermion systems, anomalous Hall angles [2] and magnetoresistance [3], 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., [4][5][6][7][8]), 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 nearequilibrium 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 ⟨n(x, t)n⟩ can be found from classical hydrodynamic equations [9].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 [10] to allow for non-Gaussian noise, Fokker-Planck equations for continuous fields (e.g., [11]), macroscopic fluctuation theory [12], and effective field theories on Schwinger-Keldysh contours [13].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 [14][15][16] and noninteracting 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, ℏ/T [17][18][19][20].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 ′ ≡ 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.
In this paper, we use the effective field theory (EFT) of diffusion [13] to systematically and quantitatively study the corrections to observables in generic diffusive systems.The EFT relies on two mild assumptions: (1) the locality of the generator of the dynamics, and (2) the thermalization of the system, i.e. the only collective excitations that survive at late times are conserved densities and associated noise fields.Therefore, it is expected to apply to a broad range of quantum and classical systems.We focus on the dynamic structure factor ⟨n(x, t)n⟩ in 1d lattice systems (>1d generalizations are shown in Appendix A 3) where nonlinear corrections are particularly strong, but also study nonlinear response (in Appendix D), 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 [21,22]; however, even the corrections to the considerably simpler autocorrelation function, ⟨n(0, t)n⟩, were only obtained analytically and observed numerically recently [23,24].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 ⟨n(x, t)n⟩, we find the universal scaling function of x/ √ Dt accompanying this correction.We also compute subleading corrections, which arise from Leading order no particle-hole particle-hole non-linear Table I: Leading order corrections to the dynamical structure factor of chaotic diffusive systems from loop corrections (nonlinear) 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 subleading correction (2-loop) dominates.These non-linear corrections are the leading corrections in one dimension (d = 1).
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/ √ Dt, summarized in Table I.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.
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 [25] 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 [26,27], 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 [6,8,28].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 lo-cal 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 1loop corrections, and in Section III.B we outline the basic steps for the 2-loop calculation required to get the nonlinear 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 the concluding Section V, we consolidate our findings, discuss the relevance of these results for the field, and outline potential avenues for future research.In Appendix D we explore the nonlinear response through a simple threepoint function, which offers a complementary test for the validity of EFT.

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 [29] as well as fluctuation (or 'loop') corrections [30] in the hydrodynamic description.While these corrections have been appreciated in the context of quantum many-body systems for some time (see, e.g., [21]), 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 [13].For example, the leading correction to diffusive correlation functions in one dimension comes from a 1-loop correction and takes the form where F 0,0 and F 1,0 are scaling functions of the scaling variable y ≡ x/ √ 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 1loop contribution [23] which, as we show below, takes the form F1,0 (y) = 4 + y 2 16 √ π e −y 2 /2 + y(y 2 − 10) 32 e −y 2 /4 Erf(y/2) We have separated F 1,0 into a universal dimensionless function F1,0 , and a non-universal factor that involves the susceptibility χ 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 ′ ≡ dD(n)/dn.If this parameter is known, e.g. by measuring the diffusivity at several densities, the entire functional form of the 1/ √ t correction to diffusion is fixed.One interesting feature of this correction is that for π ≥ 0, which implies that the autocorrelation function approaches its asymptotic diffusive form from above at late times: (3) Therefore, if a dynamic critical exponent z is extracted by fitting the autocorrelation function as ⟨n(0, t)n⟩ ∼ 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. (A1).The correction to diffusion that arises from ℓ-loop contributions at nth order in the derivative expansion in the EFT scales as 1/t n+ℓd/2 in d spatial dimensions and is encoded in a scaling function F ℓ,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 ′ 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 and comes a higher-derivative contribution F 0,1 (y), and a two-loop contribution F 2,0 (y) + F ′ 2,0 (y) log t.The former can be shown to take the form where c 1 , c 2 are non-universal transport parameters, while the latter is obtained in this paper, and is given by: with Ei(z) ≡ − ∞ −z du u e −u .In the first line, we again separated the scaling function into a non-universal factor, which now depends on D ′′ ≡ d 2 D(n)/dn 2 , and a universal scaling function.Notice that a shift in the logarithm log t → log(t/τ ) 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: 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. 6 in Appendix A.

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 [13,23].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 [30], e.g. through mode-coupling approximations or the Martin-Siggia-Rose approach [10], 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) 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 D.
A generating functional for correlation functions of densities and currents in the thermal state ρ β can be written where A 1 , A 2 are background gauge fields that couple to the conserved current in the time-evolution operator where we have collectively denoted the charge and current density by j µ = (n, j i ).Derivatives of log Z with respect to A 1 , A 2 can generate correlation functions of j µ 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. [13], this effective Lagrangian is a function of the fluctuating density n, and a conjugate field ϕ a : 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 Appendix A 2, 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 Here, σ(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 [12].In the present approach, they are just the leading terms in a general expansion in derivatives (see for example [31] 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 [32][33][34] (see [35][36][37][38] for earlier work in this direction); similar derivations may be possible for lattice systems with large local Hilbert space dimension (e.g., [39,40]), or noisy systems in the limit of strong noise (e.g., [41,42]).
When studying linear response or more general correlation functions, one expands these functions around the background value of interest for the density n = n + δn where D, D ′ , D ′′ , etc. on the right-hand side are evaluated at the background density 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 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) At tree-level, this can be evaluated using the propagators of the fields obtained from the Gaussian Lagrangian (12): Using (15), one recovers the leading diffusive behavior in (14).The second piece δG R (ω, q) comes from loop and higher derivative corrections, which are studied below.

A. 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 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 δn(t, x) ∼ q d/2 ; since the cubic nonliniearity is suppressed by an extra power of δn, it gives small corrections at late times or long distances, where ω ∼ q 2 → 0. This is in contrast The propagators ⟨nn⟩0(ω, q) and ⟨nϕa⟩0(ω, q) from Eq. ( 16) correspond to the solid lines, and half-solid half-dashed lines respectively.
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 [43].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 σ ′ .While this term leads to nonlinear response [25], in Appendix A 3 we show that it does not contribute to the 1-loop corrected two-point function; we therefore ignore it here.
The cubic vertex (17) leads to a 1-loop correction to ⟨nϕ a ⟩ shown in Fig. 2. Its evaluation, performed in Ref. [23], is streamlined here.It is convenient to amputate the external legs and parametrize the correction as δD(ω, q), namely: The one-loop correction then takes the form where we have used the short-hand notation p ≡ {ω, q}, and p ≡ dωd d q (2π) d+1 .The loop integrals can be readily evaluated in any dimension (see Appendix A 3), and give where χ ≡ σ/D is the static susceptibility, and with The general scaling δD/D ∼ q d agrees with expectations: the cubic interaction is suppressed by δn ∼ 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 /ω.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 ω→0 δD(ω, q) = 0 is required by the analyticity of equilibrium thermal correlators due to the finite thermal correlation length [44], (ii) the existence of a branch point at ω = − i 2 Dq 2 follows from a simple cutting argument [23].
We are most interested in the case d = 1: Fourier transforming this expression, or rather δG R (ω, q), is straightforward, but presented in Appendix A 5 for completeness; this results in a correction to the correlation function shown in Eq. (2).

B. 2-loop-half filling
When D ′ = 0, there is no one-loop correction to diffusion.This situation naturally arises in particle-hole symmetric systems at half filling, because D ′ ≡ dD (n)  dδn is odd under particle-hole (or charge conjugation) symmetry δn → −δn and must therefore vanish.The leading fluctuation corrections come instead from a two-loop diagram, shown in Fig. 2 arising from the quartic vertex Because this interaction scales as L (4) /L (2) ∼ δn 2 ∼ q d , and that two such vertices will be necessary to give a nonanalytic correction to the two-point function, the twoloop correction will scale as δG R /G R ∼ q 2d ∼ ω 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 ⟨nϕ a ⟩, with external legs amputated, is given by The integrals are evaluated in Appendix A 4. One finds δD(ω, q) = 1 2 with This result has the expected q 2d scaling, vanishes in the static limit lim ω→0 δD(ω, q) = 0, and features the expected 3-diffuson branch point at ω = − i 3 Dq 2 [45].For d = 1, this becomes leading to a correction to G R (ω, q) whose Fourier transform is computed in (A35) and shown in (6).

C. Higher-derivative corrections
Higher derivative corrections are also captured by the EFT for diffusion (12).These will either involve extra time derivatives ∂ t , or two extra spatial derivatives ∇ 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 [13]).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 The two coefficients c1 , c2 are contact terms and will not affect the correlation function at separated points: c1 has the interpretation of a q 2 correction to the static susceptibility χ(q) = χ + c1 q 2 + • • • , and c2 is in fact forced to vanish to guarantee G R (ω, q → 0) = 0. Fourier transforming ⟨nn⟩ = 2 1−e −βω Im G R (ω, 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 [9,46]: one writes the linearized constitutive relation for the current in terms of the charge δn up to subleading order in derivatives We have omitted a term ∂ t ∂ i δn which would have the same scaling as D 2 , because it can be absorbed in D 2 using the leading equations of motion ∂ t δn = D∇ 2 δn.To obtain response functions, one needs to know the constitutive relation in the presence of a source δµ(t, x) for charge density.Assuming that the equilibrium response is given by δn(q) = χ(q)δµ(q), with χ(q) ≃ χ+χ 2 q 2 +• • • the static susceptibility, the current constitutive relation in the presence of sources must take the form (in momentum space) where the combination δn − χ(q)δµ 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 ∂ t n + ∂ i j i = 0 and solving for δn yields a retarded Green's function G R (ω, q) ≡ δn(ω,q) δµ(ω,q) that matches (28) with

D. 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 ṅ + ∇ • j = 0 and energy (or heat) ε + ∇ • j ε = 0. Nonlinearities can now involve both densities, e.g.: The coefficient λ arises from a temperature-dependent conductivity ∂ T σ (or, equivalently, a density-dependent thermoelectric conductivity ∂ µ α).While it seems similar to the single density nonlinearity D ′ 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 σ(ω) ∼ 1 + λ 2 |ω| d/2 + • • • , as was already recognized in Ref. [21].
In order to obtain the universal scaling functions at finite q, necessary to make predictions for the structure function ⟨n(t, x)n⟩, 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/ √ 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 χ AB ≡ dn A dµ B = χδ 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): Here A, B, C, . . .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 [24,42,[47][48][49], with the role of the nonlinearity λ particularly emphasized in [24,42].In Appendix B, we show that the 1-loop correction to the density two-point function is, in one spatial dimension This produces a non-analytic correction at small frequencies to the optical conductivity 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 Erf(y/2) .

E. 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 [50], making it simple to perform precision tests of EFT predictions [25].Indeed, since the loop corrections ( 1) and ( 4) only depend on the susceptibility χ and derivatives of D(ρ), they are entirely fixed analytically and can be directly compared to numerics.
As a simple example of a lattice gas satisfying the gradient condition with a non-trivial D(ρ), we consider the one-dimensional Katz-Lebowitz-Spohn model [51][52][53], describing a collection of hard core particles hopping on a lattice with rates depending on the occupation of neighbors: with equal rates for the spatially reversed processes.δ and ϵ 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 ϵ = 0, corresponding to infinite temperature dynamics, which allows to use a random initial state as a thermal state (taking ϵ ̸ = 0 instead requires prethermalizing the system, making numerics more costly).In this situation, the susceptibility is 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.Analytical knowledge of the leading order correction allows for considerably improved predictions for diffusivity when the available integration time is short, i.e. when < l a t e x i t s h a 1 _ b a s e 6 4 = " z r R q w C e e 8 + 7 < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 L w V v 8 A 8 8 j f g g b B J P F p n a q 3 p f 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 ⟨n(x, t)n⟩ ρ versus t, x, ρ 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(δ = 0.9) = 1.9 − 1.8ρ around ρ = 0.9.For that, we simulate a sample of densities ρ = (0.85, 0.86, . . .0.93).The diffusivity is parametrized by D fit = a − bρ 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: 1.9−a 1.9 ≈ (0.135, 0.081) and  < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 L w V v 8 A 8 8 j f g g b B J P F p n a q 3 p f 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.

A. Model and Methods
The conserved charge is chosen to be magnetization, where σ 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 µN = L i=1 e µσ 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 The evolution is performed by first evolving the odd bonds and then the even bonds, with two-body gates generated by the following operators, In the absence of a staggered field V the Floquet-XXZ chain is integrable [54].For our choice of parameters, J = π/4, ∆ = 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 XXZchain 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, The global noise channel is a product of local channels and is applied to the state following a period of coherent evolution, To study the linear response dynamics, we employ the weak domain wall initial state proposed by Ljubotina et al. [55], where M = trρ is the normalization constant and δ → 0 generates a weak domain wall perturbation on top of the equilibrium state characterized by the chemical potential µ.The linear response regime is characterized by a quench where the amount of injected magnetization is not extensive.Inspired by the leading non-linear correction, Eq. 1, a natural condition for the linear response regime is δ ≪ χD ′ √ tmaxD , where t max is the maximum simulation time.We use δ = 0.0005, which satisfies the condition and is also numerically checked to be in the linear response limit for all simulated times.
The Floquet evolution defined by Eq (38) breaks translation invariance, since even and odd sites are not equivalent.To simplify the analysis, we average over even and odd sites.In addition, we shift the magnetization by its equilibrium value, σ = tr (σ z i ρ(µ, 0)) = tanh µ and normalize its initial magnitude to 1/2, where j ∈ {1, 2, . . ., L}, and L = L/2.In this normalization, the initial state profile is s j≤L/2 = −0.5 and s j>L/2 = 0.5 for all values of µ and δ.Since we doubled the lattice spacing, diffusivity and static susceptibility, , are rescaled accordingly, D → D/4, χ → χ/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 → s(x), since the hydrodynamic corrections are defined in continuum.
A simple relation between the domain-wall quench for δ → 0 and the dynamical structure factor was originally derived in Ref. [56].In the continuum limit, and under the conventions described in the previous paragraph, the relation simply reads ⟨σ z (x, t)σ z ⟩ c = χ ds dx , where the subscript c stands for the connected part of the correlation function and the average is performed over the equilibrium state, ρ(µ, δ = 0).Distance x is measured from the position of the domain wall at 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.Therefore, the magnetization profile in the domain-wall quench takes form, where F m,n are the functions described in Section II.The leading order diffusion is given by F 0,0 = 1 2 Erf(y/2) and the corrections are presented in Eq. (A37).
The diffusive corrections are numerically explored by simulating the dynamics of large system sizes, using tensor network techniques.We employ the matrix product density operator (MPDO) [57,58] representation of the state, and evolve it with time evolving block decimation (TEBD) [59] algorithm.The simulations are performed using the ITensor library [60].We simulate the dynamics from equlibria states with chemical potentials µ = (0, 0.01, . ..).In our simulations we evolve up to a time T = 400 and fix the system size to L = 2T in order to avoid finite size effects.The numerical results are shown to be convergent for the bond dimensions employed.
To probe the strength of the non-linear corrections we calculate the prefactor of F 1,0 defined in Eq. (A37), where D ′ = dD dσ and the magnetization profile after subtracting the leading order term, ∆s(y) = s(y) − F 0,0 .
The estimation of diffusivity D(σ), is achieved by three different fitting schemes, labelled (I, II, III).In the following, we present a brief overview of each scheme and further elaborate on the novel method II, which for reasons that will become clear later, is more accurate than leading order diffusion fitting and can be applied to any diffusive system efficiently.Methods I and III are elaborated in Appendix C. I is a scheme in which the approximate diffusivity is extracted by fitting the dynamical structure factor at the largest time with the leading order profile F 0,0 .II, additionally takes into account the EFT corrections F 1,0 , F 0,1 , F 2,0 that were explicitly computed in Sec.III.III 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, and in fact account for more corrections to leading order diffusion than II, without having to compute them explicitly.Namely, all zeroloop and one-loop higher-derivative corrections F 0,n and F 1,n , as well as all ℓ-loop zero-derivative corrections F ℓ,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.Additionally, III is less accurate than II in systems with multiple conserved charges as explained in Appendix C.
We now further elaborate on II, which is a general fitting scheme for the estimation of D(σ).II is inspired by our classical simulations (see Section III E), where we found that diffusivity estimation by finite time simulations is more accurate when the corrections to the leading order diffusion are taken into account.The fitting is performed as follows: We simulate the dynamics for different equilibria magnetizations and store the local magnetization values for different sites at different times, s(σ i , x j , t k ), where the subscripts i, j, k denote different samples in the discretized dataset.The diffusivity is estimated by fitting the numerical dataset with the the EFT function Eq. ( 44) using a simple least-squares method, The functional minimization over D(σ) is simplified by employing a Taylor expansion around half-filling, 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 our system (D(σ) = D(−σ)).The parameters ⃗ c = (c 1 , c 2 , . ..) are non-universal parameters arising from linear-fluctuations described in section III C. Since we will only use leading and subleadingorder corrections, we just require the two parameters ⃗ c = (c 1 , c 2 ) defined by Eq. ( 5) and are present in F 0,1 .

B. Results
Dephasing-To establish the efficiency of the fitting methods and the accuracy of EFT predictions, we switch off the staggered field perturbation (g = 0) and simulate the Floquet-XXZ chain, Eq. ( 41) in the presence of dephasing with γ = 0.1.Diffusive transport induced by dephasing has two distinct features.First, the single particle limit (equivalently the non-interacting limit ∆ = 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/γ, D ∝ J 2 /γ [41] is independent of magnetization density σ and hence D ′ ≡ ∂ σ D ≃ 0.
In Figure 4.A we show the dependence of diffusivity on magnetization, D III ∼ D II ≈ 5.35 + 3.03σ 2 + 0.54σ 4 .For the fit II we have employed the terms (F 0,0 , F 1,0 , F 0,1 ).In contrast to the non-interacting (∆ = 0) limit, diffu-sivity does have a magnetization dependence due to the presence of non-linear corrections.However, the 1-loop correction F 1,0 is small at timescales of order t ∼ O(100), since C 1 ∼ 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 F 0,1 , which dominate at these timescales for all fillings, despite being suppressed by a factor of 1/ √ 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 ∆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 ∝ O(1).Unsurprisingly, the fitted diffusivity is then considerably affected by the method used as shown by the d = 256, T = 400 fits: D II ≈ 1.77 + 13.7σ 2 + 220σ 4 , D III ≈ 1.73 + 27σ 2 + 190σ 4 .For the fit II we have employed the terms (F 0,0 , F 1,0 , 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 simulations are 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 ∆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 to corrections with similar profiles.Figure 4.D shows that for finite magnetization σ = 0.2, ∆s scales as t −1/2 , and has a closely matching profile to that of the expected 1loop correction F 1,0 .We again note that the 1-loop profile is completely determined by (D, ∂ σ 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, σ = 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 nonlinear corrections suggests that the leading corrections, F 2,0 and F 0,1 will be of similar magnitude at intermediate timescales.Eventually F 2,0 , will dominate due to its logarithmic divergence in time.For this reason we perform a different fit II around σ = 0, by including a few points σ ∈ (0, 0.01, 0.02, 0.03) and all terms (F 0,0 , F 1,0 , F 0,1 , F 2,0 ). Figure 4.D indeed shows that both F 2,0 and 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 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 γ < ∆, J the system's behavior is similar to the strong noise limit γ ≫ ∆, J, where the equation of motion for the conserved charge can be perturbatively derived [41,42].In that case, the diffusivity has weak dependence on equilibrium magnetization σ, leading to smaller loop corrections and faster thermalization: indeed the leading order non-linear terms in the diffusivity only appear at 3rd order O(∆ 2 J 2 /γ 3 ).In our situation, even if perturbation theory in 1/γ is not strictly valid, we observe the same behaviour: Weak non-linear corrections, F 1,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 F 0,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 σ > 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. 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 [61][62][63] (see [24] for preliminary work in this direction, and [64] for related work), as well systems featuring subdiffusion without dipole conservation [65,66].Higher-point functions of local operators offer useful information in this regard.Indeed, we show in Appendix D, that these are controlled by the same EFT parameters (in particular, D ′ ) 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 ⟨n(q, t)n⟩, which instead are not suitable for precision tests of thermalization because they receive large fluctuation corrections [45].
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 [67][68][69]).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 σ increases.This obstructed exploring magnetizations beyond σ ≈ 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 multipoint correlation functions such as the ones explored in Appendix D. 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 [42,[70][71][72][73] 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. [71] approximates correlation functions based on extrapolations of Lanczos coefficients, which by design produces a meromorphic G R (ω, 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 [13,23].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 [44].

APPENDIX Appendix A: EFT details 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: where F ℓ,n are scaling functions of the scaling variable y ≡ x/ √ Dt.The overall form of the expansion (A1) 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 ∇ 2 ∼ 1 x 2 ∼ 1 t at late times.Loop corrections instead come from nonlinearities in the dynamics of the densities: a single cubic nonlinearity is suppressed by δn ∼ q d/2 ∼ 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 ℓ-loop contributions at nth order in the derivative expansion will give a correction to correlation functions suppressed by 1/t n+ℓd/2 (up to logarithms); this correction comes with a dimensionless scaling function F ℓ,n and is shown in the ℓth line and nth column in Eq. (A1).This general structure of corrections to hydrodynamics applies not only to density two-point functions, but also to higher point functions [25], as well as to correlators of arbitrary microscopic operators that have the same quantum numbers as (composites of) densities [24,45].
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 ℓ,n in Eq. (A1) 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 ℓ,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 (A1), 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.

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.[13,74], 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 This is achieved by introducing phases ϕ 1,2 that always enter in the combination A µ + ∂ µ ϕ (sometimes called the 'Stückelberg trick'): (A3) This now satisfies (A2) for any functional L, because a gauge transformation can be absorbed through a redefinition of the dynamical fields ϕ 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 L is a local functional of the fields ϕ 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 [74] One advantage of fields in this basis is that they satisfy the 'latest time' property: correlators where the latest time is carried by a ϕ a field vanish due to cyclicity of the trace in Eq. ( 9) 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 L in (A3) leads to an effective description of a thermalizing system in the symmetry broken phase (a dissipative superfluid).
To describe the normal phase, Ref. [13] proposes to forbid the propagating sound mode by imposing an additional symmetry: See Refs.[32,33] for discussions on this symmetry in a holographic context.Recently, Ref. [75] proposed a slightly different approach that bypasses the need to impose this somewhat artificial symmetry, by viewing the density n r rather than ϕ 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 in an expansion in fields and derivatives.There are few additional constraints to impose, such as unitarity ] (which simply follows from the definition ( 9)), and KMS symmetry; we refer the reader to Ref. [13] for details.To leading order in derivatives, the action can be expressed as [25] is the electric field, and B aµ = A aµ + ∂ µ ϕ a .We have changed variables from ϕ r to the density n.The ellipses denote higher derivative terms (the most important of which are discussed separately in Sec.III C), as well as nonlinear terms that contain higher powers of ϕ a , which are more irrelevant than the nonlinearities considered here.
Setting the background fields to zero A 1,2 → 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 In particular, when σ(n) = σ = const, the retarded Green's function is simply related to the ⟨nϕ a ⟩ propagator The two-point function can be obtained from G R as usual from a fluctuation-dissipation relation (A11)

1-loop calculation
The universal leading 1-loop correction to diffusion was computed in Ref. [23] (see also Refs.[31,44,[76][77][78][79] 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 (ω, 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 [80], 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. 5A, which cannot produce novel singular IR structure.One therefore only needs to the cubic vertices of the EFT: Diagramatically, one expects these will generate 1-loop corrections proportional to σ ′2 , σ ′ D ′ , and D ′2 .The first is shown in Fig. 5B and can easily be seen to vanish: indeed, the only such diagram involves a loop where all poles in frequency lie on the same half complex plane.Alternatively, this can be seen to vanish in time domain using the latest time condition (A5).Indeed, the diagram involves the computation of a correlator ⟨(nϕ a )(t 1 )(nϕ a )(t 2 )⟩ -since a ϕ 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. [81].There is in fact a more general argument, showing that the EFT (12) with D(n) = const but σ(n) arbitrary has no loop corrections to the two-point function; acting with the diffusive kernel on the two point function, one has where we used the equation of motion δL/δϕ a = 0 (this equation holds up to contact terms ∝ δ(t)δ(x).The righthand side vanishes because a ϕ 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 This is known to occur in certain lattice gas models [50].
Here we have temporarily ignored higher-derivative corrections, which will enter as usual through F 0,n as in (A1).
We have established that the σ ′2 contribution to ⟨nn⟩(ω, q) vanishes.One can in fact show that the σ ′ D ′ contribution vanishes as well, although the argument is slightly more subtle.Diagramatically, the two diagrams that give corrections to ⟨nn⟩(ω, q) are shown in Fig. 5C.They can be shown to cancel by explicit calculationhowever, 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 ′ vertex, the remaining object to be computed is a two-point function between n and the normal-ordered composite operator n 2 : Crucially, this object is to be computed in the theory with σ ′ as its only cubic interaction (since the D ′ 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 (A13), one finds again that the result vanishes, which implies that this correlator must be proportional to the diffusive two-point function The diagram therefore at most renormalizes χ or D, without producing new non-analytic structures.
We are left with the D ′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 ⟨nϕ a ⟩(ω, q) by (A10) (given that we have shown that one can set σ(n) = σ = const in the action).This is done in Sec.III, where it is shown that the loop can be expressed as a correction D → D + δD(ω, q) with ( 19) where in the second line we inserted the propagators (16) and evaluated the integral over frequencies dω 2π .Here χ ≡ σ/D.Changing integration variables to q ′ → k ≡ q ′ − 1 2 q and defining this becomes The term in the first line is a UV-divergent contribution to the diffusivity, and can be absorbed with a counterterm D → D + δ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 ≥ 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 One therefore finds Eq. (20).

2-loop calculation
In systems with charge conjugation symmetry, the EFT must by invariant under This forbids cubic terms in the EFT: in particular, D ′ , σ ′ = 0.The leading nonlinearities are instead quartic, and can be found again by expanding (12): The leading fluctuation correction to transport then comes from 2-loop diagrams involving the two quartic vertices, such as those in Figs. 2 and 5D.These have not been computed before -we evalute them below.The contribution proportional to σ ′′2 vanishes due to the latest time condition (A5).The contributions proportional to σ ′′ D ′′ can also be shown to vanish, following the same argument as in the previous section: they involve computing ⟨n(t, x)n 3 (0, 0)⟩ in the theory with only the quartic interaction σ ′′ .The only remaining contribution is the one proportional to D ′′2 .It is studied in Sec.III and leads to a correction (see Eq. ( 24)) where in the third line we used the result of the 1-loop calculation (A17), and in the last one we performed the ω ′ integral by residue (note that α d (q ′2 − 2iω ′ D ) is analytic in the upper-half ω ′ plane).Changing integration variables to where we defined Let us first focus on d odd, where one can write with ) .The corrections then takes the form with This final integral has several power-law UV divergences: these are analytic in ω, 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 Note that the c d,n contribution is always analytic, so that it can be ignored.The coefficients of interest are found to be This second equality implies that the q 2 z d−1 log 1 z correction vanishes, so that δD(ω, q) is proportional to (−iω)z d−1 log 1 z -this guarantees that static correlators ω → 0 are analytic as expected [44].We are then left with One can repeat this calculation in d even, where instead of (A26) one has α d (z) = ãd z The result in general dimension therefore takes the form (25).

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 where ⟨nn⟩(ω, q) can be obtained from the retarded Green's function using (A11) where we have expanded to linear order in δ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 Let us now turn to the Fourier transform of the correction.Introducing the dimensionless variables w = ω/(Dq 2 ) and τ = tDq 2 , one has (A35) 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 The leading order Gaussian spreading is F0,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.
where we defined z = 1 2 −iw and deformed the contour to pick up the discontinuity across the two-diffuson branch cut, Discf (z) = f (z + i0 + ) − f (z − i0 + ).Evaluating the integral gives (A36) with Erfi(z) ≡ Erf(iz)/i.To perform the final Fourier transform dq 2π 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 (A35), and evaluates the discontinuity across the 3-diffuson branch cut to obtain where Ei(z) ≡ − ∞ −z du u e −u (ExpIntegralEi[z] in Mathematica).We were not able to express the final Fourier transform dq 2π 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. 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, 4 (−6 + y 2 )Erf(y/2) 16 The parameters in the term F 0,1 are defined in order to absorb the numerical constant.Equation (D2) follows from a simple manipulation of the three-point function derived in [25], 4 (Erf(y 2 /2) + sign(y 1 − y 2 )) + where y 1 = x 1 / √ Dt, y 2 = x 2 / √ Dt, and 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, Appendix B: EFT for multiple diffusing densities

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 A µI , where I = 1, 2 denotes the SK contour, and the index A = 1, . . ., N .The action will be made up of the gauge invariant combinations Imposing diagonal shift symmetry and KMS as before, the quadratic action to leading order in derivatives is where χ AB , σ AB are matrices, and dots denote matrix multiplication.Both χ and σ have to be symmetric by time-reversal symmetry (Onsager relation).The cubic action to leading order in derivatives is The cubic interactions arise as before as dependence of transport or thermodynamic parameters on potentials: .In terms of these variables, the full action up to cubic order is The derivatives of transport parameters are now taken with respect to densities, e.g.
and the diffusion matrix has been defined as 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 ω = −iD A q 2 .While the first two lines in L (3) (B5) are analogous to the N = 1 case, the term in the third line is qualitatively new.It is a contribution to j ir = δL/δ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.

Non-abelian densities
The EFT can also be straightforwardly generalized to non-abelian Lie groups [24,75].We will focus here on SU (2) for concreteness.Since 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 δ AB , ϵ ABC ).The nonlinear transformation of ϕ A a requires more attention: instead of using B µa = ∂ µ ϕ a + A µa , one should use the Maurer-Cartan form 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 The term in the last line of (B5) produced the cubic term above: we have written its coefficient as ∂ A σ BC ≡ χλϵ ABC .Note that there are also cubic terms in the first line, coming from expanding (B8): however, one can show using the leading order equation of motion that they will not contribute to the 1-loop correction studied below.

1-loop correction
Let us study the 1-loop correction to ⟨n r n a ⟩.There are two contributions: the first comes from the nonlinear piece in r n B a ⟩ = δ AB ⟨n r n a ⟩, this leads to the following correction to ⟨n r n a ⟩: The other contribution comes from two insertions of the cubic interaction S (3) : Summing these two diagrams one finds (in Like before, there are several consistency checks that this result satisfies: it vanishes when q → 0 (as it must by current conservation), and when ω → 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 [42] -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.

Fourier transformation
We'd like to compute the fourier transform of Fourier transforming dω 2π e −iωt by picking up the cut as usual, one finds The final Fourier transform dq 2π 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.

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 [82].Operator entanglement is a generated by quantum dynamics and is expected to increase linearly with time S op ∝ t, corresponding to an exponential scaling of the required classical resources (bond dimension of local tensors) [83,84].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 [70,85].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, γ, in the system, coherences are destroyed at a timescale t γ ∼ 1/γ.This leads to a saturation of S op .Therefore, if the bond dimension is enough to produce accurate results up to t γ , it will also < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 L w V v 8 A 8 8 j f g g b B J P F p n a q 3 p f   .At µ = 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 µ = 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 µ = 0.3 the profiles agree well, however deviations in (B), become significant for t ≥ 200, and therefore, we don't fit these data in the main text.
be accurate for t ≫ t γ .We confirm this for the results shown in the main text (γ = 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 7.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 7.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.

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 E suggest that when dealing with limited resources, making use of our knowledge of the general struc-ture 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, where an explicit time dependence on diffusivity α σ,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. (C1) 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 σ on the spatial profile of magnetization s(x) defined by Eq. (43).One can then estimate diffusivity from the longest simulation time t = T .How-ever, the true diffusivity is defined as D = lim t→∞ α σ,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(σ) is extracted, the non-linear corrections are either calculated by discrete derivatives, or analytically, by first performing a low order polynomial fit, where only even powers appear due to the particle-hole symmetry in the system.In Figure 7.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(σ, 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. (C2) and similar trial functions for the linear corrections, c 1/2 = M i=0 b 1/2 i (σ) 2i .We note that c 1/2 are taken constant in the coherent simulations since the maximum reliable chemical potential is ≈ 0.25 and the dependence of the constants with chemical potential is weak.Figure 7.D shows the correction to the leading order diffusive profile at T = 150.We observe that for sufficiently low µ ≤ 0.25 the correction is more than one order of magnitude stronger than the bond dimension difference (Figure 7.B), indicative of a quantitatively accurate result.For the incoherent simulation where we fit µ = 0 − 1 it is important to allow for the constants to depend on σ.
The third method (III) studies the relaxation of the total current J = dx j(x).On general grounds, the constitutive relation for the current density is j = −D(n)∂ x n + higher derivatives. (C3) 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, where in the second line we expanded the density around the equilibrium magnetization n(x) = σ + s(x), and D(σ) is approximated by Eq. (C2).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) ≃ χ(σ)δ.Let us now turn to the higherderivative terms in Eq. (C3).Because of the argument above, ∂ 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 ⊃ ∂ 2n+1 x s, which otherwise would have lead to corrections of the form F 0,n .This also includes certain nonlinear higher derivative terms: for example, all terms involving two diffusive fluctuations are total derivatives s∂ 2n+1 x s = ∂ x (. ..), so that 1-loop corrections with any number of derivatives F 1,n do not contribute.The first correction to (C4) 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 ⊃ s 2 ∂ 3 x s [21].This will lead to a ℓ = 2 loop correction at order n = 2 in derivatives (F 2,2 ) to the current decay J(t) scaling as 1/t n+ℓ/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 [24]).Because the convergence to the late time value of the current lim t→∞ 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 F 1,0 , scaling as 1/t n+ℓ/2 = 1/ √ t.The diffusivity D(σ) is therefore obtained from the late time current using (C4), and fit as a function of magnetization using (C2), 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.(C3).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 8, however method I clearly underestimates diffusivity, leading to larger deviation in the corrections.DSF can be important, it is always subleading at long times.Here, we go one step further and explore manybody correlations which would vanish in the absence of nonlinearities in hydrodynamics.In particular we explore the late-time behavior of the observable, 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 [25].According to the EFT, the asymptotic behavior of s 3 (x, t) in diffu-sive systems with a single conserved charge is given by a universal scaling function As expected for multi-body functions in linear response, s 3,EFT (y, t) is vanishes as t → ∞.It also vanishes in systems which are particle-hole symmetric (D ′ = 0).
As shown in Figure 9, for both dephasing and staggered perturbations, s 3 scales according to the EFT prediction, 1/ √ t.In the presence of dephasing, we observe a precise late time agreement of the correction profile to s 3,EFT , which verifies the validity of the EFT prediction.No fitting parameter was used in this test, as D and D ′ 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 (D2) and the numerical data can be illustrated by performing a fit of s 3,EFT for y > 1, with now D and D ′ taken to be independent fitting parameters.The result is illustrated in the inset of Figure 9.B, the two profiles agree well for (D ∼ 1.85D III , D ′ ∼ 0.36D ′ III ).This agreement suggests that the staggered profile at y ≫ 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 ′ 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, τ eq ≳ 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 [86].(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 1 :
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 τ = χ 2 D ′2 16πD 4 , together with a universal scaling function of x/ √ Dt, see Eqs. (1) and (2).
r H A e r J e r b d x 6 5 w 1 m d l B f 2 C 9 / w A 5 I p R N < / l a t e x i t > (B)
1.8−b 1.8 ≈ (0.12, 0.084), where the first/second number in the parenthesis denote the fit without/with the 1-loop < l a t e x i t s h a 1 _ b a s e 6 4 = " E J I W w 6 6 M c K v D 7 C 2 r t w k u n j 5 / V X I = " > A A A B 8 n i c b V D L S g N B E J z 1 G e M r 6 t H L Y A h E h L A r + M B T I B f B S w T z g G Q J s 5 P Z Z M j s 7 D L T K 4 Y l P y B 4 9 O J B E a / + j N 7 8 C b / B 2 S Q H T S x o K K q 6 6 e 7 y I s E 1 2 P a X t b C 4 t L y y m l n L r m 9 s b m 3 n d n b r O r H A e r J e r b d x 6 5 w 1 m d l B f 2 C 9 / w A 5 I p R N < / l a t e x i t > (B)

Figure 4 :
Figure 4: Tensor network simulations of driven XXZ chain with (A,B) decoherence, γ = 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→ II→III, see Appendix (C 2).(B,D) Corrections to diffusion for d = 256, evaluated at different times t = (40, 80, . . ., 400), denoted by dark blue→yellow colors.Black dashed line denotes linear corrections F0,1.Red dashed line denotes F1,0 in (D), and F0,1 + √ 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 F0,2+F1,0, where the logarithmic in time component of F0,2 is evaluated at T = 400.

Figure 5 :
Figure 5: 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 σ ′2 vanishes due to the latest time condition (A5).(c) The two 1-loop contributions proportional to σ ′ D ′ cancel.(d) A similar cancellation happens at 2-loop for the σ ′′ D ′′ contribution.

d 2 − 1
log 1 z .One finds the same results as above, apart for a sign b d,n → −b d,n .

Figure 6 :
Figure 6: Universal scaling functions describing non-linear corrections to the diffusive structure factor as a function of the hydrodynamic variable y = x/ √ Dt.The leading order Gaussian spreading is F0,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.
r H A e r J e r b d x 6 5 w 1 m d l B f 2 C 9 / w A 3 m 5 R M < / l a t e x i t > (A) < l a t e x i t s h a 1 _ b a s e 6 4 = " b B r H A e r J e r b d x 6 5 w 1 m d l B f 2 C 9 / w A 5 I p R N < / l a t e x i t > (B) < l a t e x i t s h a 1 _ b a s e 6 4 = " E J I W w 6 6 M c K v D 7 C 2 r t w k u n j 5 / V X I = " > A A A B 8 n i c b V D L S g N B E J z 1 G e M r 6 t H L Y A h E h L A r + M B T I B f B S w T z g G Q J s 5 P Z Z M j s 7 D L T K 4 Y l P y B 4 9 O J B E a / + j N 7 8 C b / B 2 S Q H T S x o K K q 6 6 e 7 y I s E 1 2 P a X t b C 4 t L y y m l n L r m 9 s b m 3 n d n b r O o w

Figure 7 :
Figure 7: Estimation of stability of corrections for different fillings µ = 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 µ = 0.3.(B): Difference between the profiles at times t = 100, 200, 400, denoted by full,dashed and dotted lines respectively.For µ = 0.3 the difference is almost an order of magnitude larger than µ = 0, 0.2, indicating the enhancement of simulation error at larger fillings.(C): Predictions of diffusivity by employing the fitting method (I).At µ = 0 diffusivity decreases with time indicating that the system is not yet at asymptotic diffusion.At larger µ we observe a monotonic increase over time, compatible with the EFT predictions.(D): Correction to diffusion, ∆s = s − F0,0 using method II and times up to T = 150 for d = 256, 600.At µ = 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 µ = 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 µ = 0.3 the profiles agree well, however deviations in (B), become significant for t ≥ 200, and therefore, we don't fit these data in the main text.
F ℓ,0 (the first two ℓ = 1, 2, coming from D ′ , D ′′ , were computed in Sec.III).All of these corrections do not contribute to correlators of the total current: defining C(n) such that C ′ (n) = D(n) and integrating from x = 1 to x = L, we haveJ ≃ C(n(L)) − C(n(1)) ≃ D(σ) (s(L) − s(1)) ,(C4) Appendix D: Non-linear responseIn the previous sections, we explored the effects of scaling corrections in the DSF which is a two-point correlation function.While the effect of fluctuations on the

Figure 8 :
Figure8: Sub-leading correction estimation for diffusivity fits using different fitting methods for µ = 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(σ).
One can change variables to the density as before ϕ A r → n A ≡ δL