Nonlinear response in diffusive systems

Nonintegrable systems thermalize, leading to the emergence of fluctuating hydrodynamics. Typically, this hydrodynamics is diffusive. We use the effective field theory (EFT) of diffusion to compute higher-point functions of conserved densities. We uncover a simple scaling behavior of correlators at late times, and, focusing on three and four-point functions, derive the asymptotically exact universal scaling functions that characterize nonlinear response in diffusive systems. This allows for precision tests of thermalization beyond linear response in quantum and classical many-body systems. We confirm our predictions in a classical lattice gas.


Introduction
Most interacting systems thermalize at nonzero temperature.Their late time dynamics is then described by the fluctuations of conserved densities, or hydrodynamics, and their global equilibration towards thermal equilibrium.
In the case where a single density is conserved, the hydrodynamic description is typically diffusive [1], with density two-point function taking a universal form ⟨n(x, t)n(0, 0)⟩ = χT (4πD|t|) d/2 e −x 2 /(4D|t|) + O(1/t which is entirely fixed at late times up to two non-universal constants: the charge susceptibility χ ≡ dn dµ , and the diffusivity D.Here T is the temperature and d the spatial dimension.The subleading corrections denote higher derivative [2] and loop contributions [3] respectively -both are irrelevant and lead to small corrections at late times.This two-point function captures linear response in all diffusive systems and organizes it in terms of a few non-universal 'Wilsonian coefficients' of the effective field theory (EFT) of diffusion [4].
Diffusive systems feature nonlinear response as well.For example, a diffusivity D = D(n) or a susceptibility χ = χ(n) that is a nontrivial function of the background density immediately implies nonlinear response.More generally, the EFT of diffusion has interactions that do not have simple interpretations as nonlinearities of transport or thermodynamic parameters (see, e.g., [5]); we will show however that to leading order at late times, nonlinear response is simply controlled by the coefficients in the expansion of D(n) and σ(n) ≡ χ(n)D(n) around the background value of the density The EFT predictions for nonlinear response therefore do not involve, to leading order, any Wilsonian coefficients beyond those measurable within linear response at several values of the background density.
In this paper, we present a simple scaling argument for the late time behavior of equilibrium connected N -point functions in diffusive systems ⟨n((N −1)t) • • • n(2t)n(t)n(0)⟩ ∼ 1/t (N −1)d/2 .Using the EFT of diffusion [4], we confirm this scaling and, focusing on N = 3, obtain the entire three-point function ⟨n(x 2 , t 2 )n(x 1 , t 1 )n(0, 0)⟩ at late times, which is the sum of two universal scaling functions multiplied by the nonuniversal parameters D ′ and σ ′ from (2).The EFT results are asymptotically exact: the only approximation lies in the late time expansion of correlators.Finally, we numerically study a classical lattice gas model and find excellent agreement with our predictions.
One motivation for this work are recent experimental developments in probing the nonlinear response of correlated many-body systems [6,7,8], and particularly progress in studying higher-harmonic generation from low-frequency (single THz) sources [9,10], which may be governed by hydrodynamics when ℏω ≲ T .Our results should apply to diffusive metals, including anomalous metals [11] and strange metals [12,13], as well as to diffusion in correlated insulators [14].We further hope that our results can serve as benchmarks for precision tests of thermalization (or lack thereof) in numerics.
Nonlinear response in thermal states has been studied for some time.Generalized fluctuation-dissipation or KMS relations were found in [15,16].Certain general properties were also studied recently through the lens of the eigenstate thermalization hypothesis in [17,18].Higher order correlation functions were computed in holographic models [19,20], models with relaxational dynamics [19,21], integrable systems [22,23] as well as more general ballistic regimes [24].Nonlinear response is also related to higher cumulants (or full counting statistics) in out-of-equilibrium states (see, e.g., [25,26,27]).

General scaling considerations
A simple scaling argument gives the general form of nonlinear response in equilibrium: from (1), one finds that densities scale as δn ∼ 1/x d/2 and x 2 ∼ Dt.Nonlinear response however requires corrections to diffusive scaling from irrelevant interactions.These leading non-Gaussianities are suppressed by fluctuations δn ∼ 1/x d/2 : a three-point function for example will therefore scale as ⟨nnn⟩ ∼ x −3d/2 x −d/2 ∼ x −2d .A connected N -point function will require N −2 such non-Gaussianities and therefore scale as where g N is a universal dimensionless scaling function (up to a finite number of Wilsonian coefficients) depending on cross-ratios of coordinates.The factor of χT can be obtained from dimensional analysis.We have chosen to parametrize the overall scaling dependence in terms of the geometric mean of the differences in times (assuming that The subleading corrections in (3) are similar to those in (1): higher derivative contributions give corrections with a relative O(1/ t) suppression, and loops give corrections with a relative O(1/ td/2 ) suppression (up to logarithms).
Eq. (1) shows that the two-point function (N = 2) indeed takes the form (3), with the universal scaling function describing linear response in diffusive systems.We study higher-point functions N ≥ 3 below, focusing particularly on N = 3, 4. We will obtain the universal form of the density three-point function in diffusive systems We will find that the cubic scaling function g 3 separates into two universal scaling functions, with non-universal coefficients proportional to D ′ ≡ dD(n)/dn and σ ′ ≡ dσ(n)/dn : with Since both D ′ and σ ′ can be independently extracted from linear response measurements at various values of the density n, the EFT produces a prediction for the three-point function with no fitting parameters.

EFT calculation of higher-point functions
To compute the nonlinear response of diffusive systems, we use the EFT of diffusion developed by Crossley, Glorioso and Liu [4], see [28] for a review and [29,30] for related work.It differs from previous approaches to fluctuating hydrodynamics [31,32] and macroscopic fluctuation theory [26] in that it provides a systematic controlled expansion in fluctuations.
In particular, it captures general nonlinearities in the noise field that are not visible at the level of constitutive relations, and are missed in other approaches.However, these terms correspond to fairly irrelevant operators in the EFT and only give subleading corrections to correlation functions. 1We therefore expect that previous approaches could also be used to obtain nonlinear correlators ⟨n(t )n(0, 0)⟩ to leading order at late times (although this has, to our knowledge, not been done).Nevertheless, in the spirit of using a controlled approach -namely one where corrections to the leading answer can be systematically obtained, order by order -we use the EFT of diffusion in this paper.
Consider the partition function of a microscopic system2 with a conservation law ∂ µ j µ = 0 (µ = 0, 1, . . ., d runs over spacetime indices) where ρ β = e −βH / Tr(e −βH ) is the thermal density matrix, and is the time-evolution operator of the system with density and current coupled to a source A µ .We have introduced two independent background sources A 1 µ , A 2 µ in (8), one on each leg of the Schwinger-Keldysh contour, to generate correlation functions of the density n ≡ j 0 (or current density j i ) with various operator orderings 3 .For example, defining the symmetric Green's function of density is obtained from One can similarly define higher-point functions, for example: with similar expressions for G rra and G raa , see Eq. (A.8).These correspond to various orderings of three point functions of densities, see, e.g., [16].
The EFT of diffusion [4] consists in representing the partition function ( 8) by an effective Lagrangian L of hydrodynamic degrees of freedom where n is the fluctuating density and ϕ a is related to the corresponding noise field.The assumption of hydrodynamics is that the only long-lived excitations in thermalizing systems are conserved densities -the entire nonlocal aspect of Z[A 1 , A 2 ] can therefore be realized through a local effective action of these long-lived degrees of freedom.Appendix A reviews the construction of this effective action.To leading order in the diffusive scaling, the nonlinear action takes the form where The ellipses contain higher derivative terms, which give O(q 2 ) corrections to observables, as well 3 Correlation functions of other microscopic operators can also be obtained from the EFT up to further Wilsonian coefficients, through operator matching equations [33].See [34] for examples in classical spin chains.as further nonlinear terms that are more irrelevant than the leading nonlinear terms.To leading order, all nonlinearities come from expanding σ(n) and D(n) around the background value of the density as in (2).For example, up to cubic order in fields, turning off background fields A µ momentarily, one has The propagators ⟨nϕ a ⟩ and ⟨nn⟩ can be obtained from the first line: where p = {ω, q} denotes frequency and momentum collectively.The terms in the second line of (15) give cubic vertices: these lead to a density three point function (see Fig. 1): where the permutations are obtained by swapping p 1 → p 2 , p 2 → −p 1 − p 2 once, and twice.
Higher-point functions can be similarly computed; the four-point function is obtained in (A.15).

Properties of the three-point function
There are a number of properties that the higher-point functions should satisfy, which serve as useful consistency checks of our calculation.First, notice that (17) vanishes when either q 1 or q 2 → 0, as well as when q 1 → −q 2 ; this must happen because the total charge is conserved and has trivial dynamics.Second, other time orderings of the three-point function G aar and G arr can also be computed from the EFT ( 14): these should be related among each other and with G rrr through nonlinear KMS relations [16].We check in App.A that these relations are satisfied.Third, in the limit where one of the densities is taken to be static, the three-point function reduces to the derivative of a two-point function with respect to chemical potential: where G ar (p) = −iσq 2 ⟨nϕ a ⟩(−p) and G rr (p) = ⟨nn⟩(p).Since d dµ = χ d dn , and G rr , G ar depend on the background density through σ(n), D(n), Eq. ( 18) further elucidates why the three-point function depends on σ ′ and D ′ ; we verify that it holds in App. A.

Fourier transform
The Fourier transform of ( 17) can be evaluated, for direct comparison with numerics or experiments that probe nonlinear response in space and time domain4 .Following Eqs. ( 6) and ( 7), we remove the overall 1/t d scaling and directly extract the universal scaling functions g 3,D ′ and g 3,σ ′ .In this section, we focus on d = 1 spatial dimension for simplicity.The scaling function g 3,σ ′ is most simple and is given by with . The other scaling function g 3,D ′ is much richer and is shown in (A.17).When all three points lie on the same site y 1 = y 21 = 0, it simplifies to Recall that we have taken t 2 > t 1 > 0. In this kinematics, the entire three-point function is given by An interesting kinematic configuration that neatly distinguishes both scaling functions is to take two of the points at equal time, t 1 = t 2 , while keeping x 1 ̸ = x 2 .There one finds from (19) that g 3,σ ′ vanishes exponentially fast as t 2 → t 1 for x 1 ̸ = x 2 and only produces a contact term ∝ δ(x 1 − x 2 ) in this limit.The entire correlator at separated points is then proportional to D ′ : where erf(s) = s 0 du e −u 2 .The fairly long-range correlations at time t in ( 22) are an equilibrium analog of long-range, equal-time correlations that arise in out-of-equilibrium states [26].

OPE in the EFT
Any EFT defined by a path integral satisfies an operator product expansion (OPE) [35] : two operators approaching each other can be approximated by an expansion in local operators of the EFT.This provides a useful organizing principle for nonlinear response in the EFT.For densities, this gives where • • • contains contributions from higher dimension operators such as ∇n, n 2 , etc.
The EFT OPE should be understood as valid in the limit where n(x, t) and n(0, 0) are much closer than any other two EFT operators, but still with a separation greater than the EFT cutoff.We have used the tree level scaling dimension n(x, t) ∼ q d/2 to obtain the scaling of the first term, and furthermore used the fact that second term requires the dimensionful couplings q . However, the next OPE function, describing the n × n → n fusion channel is more interesting.It constrains a limit of the three-point function: The fact that the dependence on x, t and x ′ , t ′ factorizes is an obvious consequence of taking the limit x 2 ∼ Dt ≪ x ′2 ∼ Dt ′ -the nontrivial content of ( 24) is however that the dependence on x ′2 /(Dt ′ ) is entirely fixed.Eq. (19) shows that the σ ′ piece satisfies this property.We show in App.A that the D ′ piece of the three-point function satisfies it as well.0.001 0.010 0.100   1) and ( 21), with no fitting parameters.Inset: dimensionless scaling function (7) g 3 ( t 2 t 1 , x 1,2 = 0) compared to theory (21).(Numerical parameters: δ = 0.9, L = 2 19 , ⟨n⟩ = 0.9, averaged over 5 realizations).

Numerics
We expect the scaling behavior (3) to apply to any diffusive many-body system, with the precise form of the scaling function (17) applying more specifically to systems with a single conserved density.This includes quantum and classical spin chains, random or deterministic unitary circuits and cellular automata with a conserved density, and lattice gases.As a simple test of our predictions, we study nonlinear response numerically in the Katz-Lebowitz-Spohn (KLS) model [36], a classical lattice gas with conserved particle number.One appeal of this model is that the diffusivity D(n) and conductivity σ(n) are known analytically as functions of the background density n [37,38] -all parameters entering our prediction for the three-point function in Eq. ( 17) or ( 21) are therefore fixed.
Let us first verify the universal scaling behavior of correlators from (3): Fig. 2 shows a clear 1/t (N −1)/2 decay of the N -point functions, for N = 2, 3.Moreover, the prefactor of this polynomial behavior agrees with the theory prediction from (21) at late times.As a more refined test of our results, we extract numerically the dimensionless scaling function g 3 (6) and compare it to the prediction ( 19) and (20) (insets of Figs. 2 and 3).In Fig. 2, the predicted scaling function is fairly featureless and difficult to confirm.To enhance its 10 -4 0.001 0.010 0.100 corrections due to loops [3].Dashed lines show the predictions from ( 1) and ( 21) without corrections, and solid lines include 1-loop corrections.For the two-point function, the exact prefactor of the 1-loop correction from [3] was used; for the three-point function, the corresponding computation is not available so the coefficient was fit.Fig. 4 shows the diagrams contributing to this correction.
Figure 4: Loop corrections giving a relative 1/t d/2 correction to the three point function: . In d ≤ 2 this is the leading correction to the three-point function.
features, Eq. ( 21) suggests to study a region of parameters where the D ′ coefficient is large; this is done in Fig. 3, where a good agreement with the predicted scaling function g 3 is found.A larger D ′ is also known to produce larger loop corrections to the correlators [3] -Fig. 3 shows that these corrections are indeed appreciable at intermediate times.

Generalizations
We have focused on nonlinear response in diffusive systems ω ∼ −iq z with z = 2, because this is the most common dissipative universality class observed in thermalizing systems.
However, there are other possibilities: superdiffusion in the KPZ universality class z = 3 2 occurs for sound modes in d = 1 [32], and subdiffusion (e.g., z = 4) can arise in constrained systems with dipole-type symmetries [39,40,41].One can extend a scaling argument similar to Eq. ( 3) for these situations: Assuming that nonlinear static susceptibilities are finite and that the entire N -point function is scaling according to ω ∼ q z in the hydrodynamic regime, one finds The scaling assumption however does not apply to hydrodynamic regimes that include a scale to leading order in dissipation, including balistic modes ω(q) = c s q − iDq z with length scale (D/c s ) 1/(1−z) .The leading nonlinearities in these situations are non-dissipative (see, e.g., [23]), but the subleading diffusive regime can be captured by generalizing the scaling argument above: note that nonlinearities from c s (n) are 1/q z−1 enhanced compared to those from D(n); using N − 2 such vertices then leads to ⟨n To go beyond this general scaling behavior and obtain the universal dimensionless scaling functions for nonlinear response, one would need to use the appropriate EFTs for the different hydrodynamic situations described above.It would be interesting to carry this out to allow for precision tests of thermalization in models exhibiting superdiffusion or subdiffusion.
relevant nonlinear terms -at the cubic level these are: Three new Wilsonian coefficients enter in the cubic action at leading order in derivatives: σ ′ , χ ′ and w.The last one involves the field strength F r,ij which is independent of the dynamical field ϕ r -it will only produce contact term contributions to correlators involving the current, so we ignore it in the following.How suppressed are cubic terms compared to the quadratic action?Fluctuations in ϕ scale as The cubic action is therefore L (3) /L (2) ∼ q d/2 suppressed compared to the quadratic one.
This is what allows for an expansion in fluctuations in the EFT: interactions are irrelevant in the RG sense.For example, they lead to O(q d ) = O(ω d/2 ) corrections to correlation functions [3] (because two cubic vertices are necessary).However, the cubic terms in (A.4) will give the leading contribution to three-point functions of densities or currents.
The density and current operators can be found by taking derivatives with respect to the background fields j µ r ≡ δL/δA aµ .For example, the density n ≡ j 0 r is given by Taking derivatives with respect to µ = A r0 , one sees that χ = dn dµ indeed is the susceptibility, and χ ′ = 1 χ dχ dµ = dχ dn , justifying the notation.Similarly computing the current density j i r shows that σ is the dc conductivity, and σ ′ = dσ dn .To make contact with other approaches to fluctuating hydrodynamics, one can use this equation to trade the degree of freedom ϕ r for n.The action up to cubic order then reads where we defined D(n) ≡ σ(n)/χ(n).This is the action used in the main text, see Eq. ( 15).
One can similarly compute interactions involving more fields -to leading order in derivatives these simply take the form (14).The only Wilsonian coefficients involved are therefore the Taylor expansion coefficients of the diffusivity and conductivity (or susceptibility) around the density of interest (2).This agrees with other existing approaches to hydrodynamics [31,32], including macroscopic fluctuation theory [26] -we therefore expect that the results in this paper on the leading late time density three and four-point functions could be obtained from these other methods (to the best of our knowledge, they have not so far).However, The momenta p 1 , p 2 are carried by the first two arguments of the Green's functions.For example, . One can check that these satisfy nonlinear KMS relations [16], which to the order we are working in read In a certain limit, the three-point functions above should reduce to derivatives of twopoint functions with respect to chemical potential (18).From (A.11), one finds that this is indeed the case i lim In systems with charge conjugation symmetry, such as particle-hole symmetric spin chains at half filling, cubic nonlinearities D ′ , σ ′ vanish.In these situations, the leading non-Gaussianities are quartic: leading to the following four-point function G rrrr (p 3 , p 2 , p 1 ) = 2T σ ′′ (q 1 • q 2 )⟨nn⟩(p 3 )⟨nϕ a ⟩(p 2 )⟨nϕ a ⟩(p 1 )⟨nn⟩(p 4 ) + 5 perm.
The dynamics we consider on this system is the Katz-Lebowitz-Spohn model [36].It consists in a Monte-Carlo type time evolution where at every step, a site i is randomly chosen, and the following '4-site gate' is applied to the charges n i−1 n i n i+1 n i+2 : where the number above the arrow denotes the probability P of the transition happening.
The particle at n i therefore hops to the right with probability depending on the configuration of is neighbors.Processes related to those shown above by reflection occur with equal probability.The model depends on two parameters δ and ϵ.We will set ϵ = 0 below: in this case the rates above clearly satisfy infinite temperature detailed balance, e −βH = 1 (the model with ϵ ̸ = 0 can also be shown to satisfy detailed balance [36,44]).Time will be measured in units of the rate r, which will therefore not appear below.We take a single time step ∆t = 1/r to correspond to L applications of the gate above.
An appealing feature of this model for the purposes of testing the EFT predictions is that D(n) and σ(n) are known analytically as a function of the density measured in units of the lattice constant n = 1 L i n i ∈ [0, 1], and the parameters of the model δ, ϵ [37,38].The coefficients of the nonlinear EFT (15) to leading order in gradients are therefore entirely fixed.For our purposes, it will be sufficient to consider the model with ϵ = 0, in which case these are given by (note that T χ is finite in the limit T → ∞) From these expressions, one can obtain the cubic interactions D ′ and σ ′ appearing in the EFT.
Figs. 2 and 3 show the numerical results for correlation functions, and the corresponding predictions from the EFT, using the values D ′ and σ ′ obtained above.

2 Figure 1 :
Figure 1: Two diagrams give contributions, proportional to σ ′ and D ′ respectively, to the density three-point function G rrr (p 2 , p 1 ).

d/ 2 Λ
to obtain the scaling of the second term.The first scaling function describing the n × n → 1 fusion channel can be simply read off the two point function (1):

Figure 2 :
Figure 2: Log-log plot of the connected two-and three-point functions of density in the KLS model.The gray curves are the theory predictions from (1) and (21), with no fitting

Figure 3 :
Figure 3: Same as Fig. 2, but with a filling ⟨n⟩ = 0.35 chosen such that D ′ is large, leading to a scaling function with clearer features (inset).A large D ′ also leads to enhanced 1/ √ t