Theory of non-linear diffusion with a physical gapped mode

In a system with one conserved charge the charge diffusion is modified by non-linear self-interactions within an effective field theory (EFT) of diffusive fluctuations. We include the slowest ultraviolet (UV) mode, constructing a UV-regulated EFT. The relaxation time of this UV mode is protected from renormalization, as supported by experimental data in a bad metal system. Furthermore, the retarded density-density Green's function acquires four branch points, eventually increasing the range of applicability. We discuss the fate of long-time-tails as well as implications for the quark gluon plasma.


INTRODUCTION
In a system with a single conservation law, the latetime dynamics is universally described by hydrodynamic diffusion of the associated conserved charge.In the EFT approach to hydrodynamics [1][2][3][4], the linear diffusion equation is derived as the equation of motion from a classical quadratic effective action.If the diffusion constant depends on the transported density, then this equation becomes non-linear [5].This is caused by the (self-)interactions encoded in the corresponding effective action, eventually affecting the correlation functions [6,7].
The onset of diffusion in a system can be characterized by the relaxation time τ , the timescale of the relaxation of that non-conserved quantity which relaxes slowest.In the limiting case τ = 0, the interacting EFT of hydrodynamics predicts long-time-tails [8] and large renormalization effects for the transport coefficients [9].In this limit, which corresponds to the instantaneous response of the current to the source, hydrodynamic equations fail to satisfy the commutation sum rules [10].Therefore a finite relaxation rate 1/τ has to be considered to regulate the theory of diffusion discussed above.Such a relaxation rate has been recently measured in ultra-cold gas of 6 Li in a two-dimensional lattice, an example of a "bad metal", as the testing ground for the Fermi-Hubbard model [11].
The aforementioned relaxation time τ characterizes the timescale of the establishment of local thermal equilibrium in a U(1) static system [12] [13].Since the local equilibration cannot happen arbitrarily quickly, it is expected for τ to be subject to a lower bound [12].It is then important to understand how this bound is affected by fluctuations.This can be investigated via studying the effect of self-interactions in the EFT framework, in particular the renormalization effects.In fact, we need to increase the EFT cutoff of [6,7] to include at least the slowest UV mode.Indeed, this leads us to derive, for the first time, a UV-regulated theory of non-linear diffusion with a physical gapped mode.
Another motivation manifests itself in the long-timetails.It is well known that in a many body system, at late times, the auto-correlation functions feature long-timetails [14,15].Now the question is what happens to these tails when the local thermalization is not instantaneous.Does the thermalization occur exponentially fast [16]?If not, would it possibly lead to a longer tail for correlation functions, thereby delaying the global thermalization?Keeping these questions in mind, in the following, we construct the EFT including fluctuations for a selfinteracting system with a single conserved density together with a single gapped mode; the longest-lived gapped mode, acting as a UV-regulator for the EFT.Our effective action is consistent with the general framework of [1], however, derived through the traditional Martin-Siggia-Rose (MSR) formalism [17].Explicitly calculating the one-loop retarded density-density Green's function, we discuss how the renormalization of the diffusion constant at τ = 0, performed in [6,7], is affected by the gapped mode.In particular we discuss the renormalization of the relaxation time which was not included in [6,7].The renormalized Green's functions then will be used to study the effect of the gapped mode on longtime-tails.We conclude by discussing application of our results to a bad metal [11], expanding quark gluon plasma (QGP) and also to the QGP droplet near the QCD critical point.Natural units, k B = 1 = ℏ, are implied.

SETUP
Beginning in the absence of interactions, we consider a U (1) conserved charge with a non-relativistic non-conserved current as follows with D being a constant.For finite τ , the combination of the two equations in (1) leads to This equation describes the dynamics of the diffusive density mode together with a gapped mode The theory has two timescales: τ U V ≡ −i/ω 1 is the timescale of the thermalization while τ D ≡ −i/ω 2 is the diffusion time of the mode with momentum k.For the long wavelength modes with τ ≪ (Dk 2 ) −1 : which gives τ D = (Dk 2 ) −1 and τ U V = τ .This can be compered to Fick's law of diffusion, ∂ t n − D∇ 2 n = 0, which encodes one single diffusion mode ω = −iDk 2 .
Since ω 2 is a UV mode for Fick's law, we refer to (2) as the (phenomenological) theory of diffusion, UV-regulated with a physical gapped mode.
Our intent is to make equation ( 2) non-linear.When n does not couple to any dynamical field, this is achieved by including the self-interactions of n [18].These enter by promoting D to be a function of small fluctuations about n = 0: Dropping terms with orders higher than 2, the non-linear version of ( 2) is found to be In general τ could be a function of J in (1).However, since there is no other vector involved in the setup, the scalar τ cannot depend on the vector J at linear order; the first contribution is quadratic in J and turns out not to contribute to any of our results.We take τ constant in what follows.Note also that if n refers to a charge density, then under its changing sign, λ D → −λ D but D and λ ′ D remain the same.Our goal is to study correlations between fluctuations of n.For this, we construct an effective action for n whose equation of motion is (5).We do this in the framework of MSR formalism [17] .The idea is to put a noise term on the right side of (5) and then impose the fluctuationdissipation theorem to fix its strength.Finally, exponentiating this stochastic equation yields the effective action S eff = dtd d xL.As we show in the Supplemental Material, the effective action takes the form where the charge conductivity σ(n is related to the diffusion coefficient by the Einstein relation σ = χD via the charge susceptibility χ.To quartic order in fields it becomes Here n a is an auxiliary field, analogous to the a-field in the Schwinger-Keldysh framework [6].In the limit τ = 0, (7) reduces to the theory of diffusive fluctuations [6].At τ ̸ = 0, (7) defines our novel theory of non-linear diffusion with a physical gapped mode.
In order to study the leading effects caused by the nonlinear terms in (5), we investigate the effect of one-loop corrections on the retarded Green's function in the EFT described by (7).At one-loop order, only cubic interactions in (7) contribute.This is why we dropped higher order terms when expanding D(n) and σ(n).Even λ ′ D and λ ′ σ are not needed.

RESULTS
The retarded Green's function at one-loop can be parameterized by The self-energy is given by where from the first line in (7) we have Considering a hard momentum cutoff, we evaluate the loop integrals.The cutoff-dependent parts of Σ can be absorbed into the bare values of coefficients D and τ .The cutoff-independent parts turn out to be non-analytic, that are given in d = 1, 2, 3 dimensions by ) Here and below f id are analytic functions of ω and k.
The non-analyticity is encoded in with . Having found Σ, the next step is to calculate the loop correction to the numerator of G R, (1) nn .Similarly, we find (see the Supplemental material) The two expressions (10) and (12) fully specify (8).Rewriting it in the form and δτ (ω, k) = 0.The latter answers the first question raised in the Introduction: our theory, Eq. ( 7), predicts that the bound on local thermalization time τ (and thus the onset of diffusion) [12] is protected from renormalization caused by the fluctuations.In the Discussion, we show that this result is supported by the experimental data associated with the bad metallic system of [11].Note that at τ = 0, (14) reduces to the results found in [6,7].Causality implies that G R is analytic in the upper half of the complex frequency plane.However, the structure of F (ω, k) indicates that, due to the interactions, four branch point singularities are induced in the lower half plane (Fig. 1): (15) These branch points correspond to the minimum energy for generating two on-shell excitations with the dispersion (3) in the loops, as explained in the Supplemental Material.The location of branch points will be important when finding the inverse Fourier transform of momentum space correlation functions.
As the last formal result, we find G nn (t, k) = dω 2π G nn (ω, k) e −iωt .Note that G nn (ω, k) can be simply calculated from G R nn (ω, k) via the fluctuation-dissipation theorem.While at τ = 0 one finds G nn (t, k) analytically [19], at τ ̸ = 0 we were not able to represent G nn (t, k) in terms of known special functions.Thus, we find this function numerically.
The results for λ 2 eff ≡ T χ (τ D 5 ) 1/2 λ 2 D = 1 2 are displayed in Fig. 2. For t ≫ τ D , r = τ U V /τ D has no significant effect; the two solid line curves converge similarly.However, the effect of r is significant at t ≪ τ D .As shown in the inset, for r = 1  10 , the correlation function around t = 1 10 τ D decays slower than for r = 1 100 .This is related to the second branch point in Fig. 1, i.e., w22 .As mentioned earlier, we did not find G nn analytically; but one would expect that when evaluating the Fourier integral, among other contributions, picking up the pole w22 would yield a term ∝ e −iω22t ∼ e −t/τ U V .This would give a finite contribution at t ≲ τ U V ∼ rτ D , which justifies the excess shown in the inset.Obviously, when τ U V → 0 [6,7], w22 becomes −i∞, so the effect disappears.We emphasize that this feature is valid in static systems.In a dynamical background, the perturbations may not decay solely by an exponential term [20].
In the limit τ D ≪ t, we find the asymptotic behavior of G (1) nn (t, k) analytically.This "long-time-tail", illustrated in red (r = 1  10 ) and yellow (r = 1 100 ) in Fig. 2, comes from integration over the region surrounded by the dashed semicircle in Fig. 1.For t > 0 and d = 1: This answers the second question raised in the Introduction: the fast thermalization of the non-conserved current J does not affect the fractional power of the long-time-tail.We find t −1/2 , which is identical to the known behavior t −d/2 without UV regulator [8].But the exponential decay associated with diffusive fluctuations changes.When τ = 0, the late-time behavior is given by e − 1 2 Dk 2 t [21].The effect of the gapped mode is to decrease this factor at any k.We emphasize that this exponential decay is specific to static systems near equilibrium.See [20] for a discussion on the existence of late-time universal attractors in a longitudinally expanding QGP.

Bad metal:
Let us check the results of our theory against data from the bad metallic system [11] mentioned in the Introduction.In [11], the charge density n(t, k) is identified with twice the experimentally measured atomic density of one spin component.The data from [11] consists of eight groups of points, each group of points contains ten n(t, k) values, with a fixed value of k.Ref. [11] fits the data using the analytical solution of (2) to determine Γ, D and δn ≡ n(0, k).However, at high temperatures, which is the case in [11], the thermalization length becomes short, on the order of the lattice spacing, and fluctuation corrections are of order one [6].Therefore, instead of fitting data to the linear equation ( 2), as done in [11], we propose the non-linear equation (5).Fig. 3 illustrates the results of this fit (blue) compared to the results of the linear equation fit (orange and red).For the sake of consistency, in both cases we fit data with the numerical solution of the equation.For the linear case we find three parameters as in [11], while in the nonlinear case we also find a fourth parameter: λ D .Our non-linear fit procedure and the dependence of our results on it are detailed step by step in the Supplemental Material.
As shown in the top panel, the linear fit of blue dots (blue line) is coincident with the linear fit of orange dots (orange line).Since each pair of orange-blue dots is associated with a specific k, this coincidence implies δΓ(k) ≈ 0. This may be regarded as the realization of our earlier statement that τ = Γ −1 is protected from the non-linear corrections in our model.
The bottom panel of Fig. 3 indicates that δD(k) is nonzero for almost all values of k in the experiment of [11].In other words, the diffusion constant is significantly renor- versus the amplitude of the solution, namely n(0, k).Orange and red dots are obtained by fitting the data of [11] to Eq. ( 2) while blue dots stem from fitting the data to Eq. ( 5).The orange, red, and blue dashed lines are the corresponding linear fits.Γ and D are normalized by the extrapolated zeroamplitude values Γ0 and D0, respectively.Error bars indicate the standard error of the mean value resulting from our fits.
malized by non-linear effects.This is in agreement with the prediction of our theory.As pointed out in [11], their fit of large-initialamplitude data with the linear equation is merely an attempt to evaluate the possibility of non-linear effects.We now see that the data is qualitatively consistent with our non-linear equation.Due to the large fluctuations existing in the system [6], non-linearities must be taken into account.Therefore our non-linear fit is preferred to [11]'s fit as our fit has at least two quantitative advantages: First, the value of D obtained here is more reliable than that of [11].Second, the value of λ D = dD(n)/dn can only be extracted from the non-linear fit of the data (see Supplemental Material).
Bjorken expansion: A simple dynamical system to study the effect of a UV-regulator on diffusive fluctuations is the Bjorken flow, which is a hydrodynamic model for the longitudinal expansion of the QGP.
The effect of non-linear fluctuations in Bjorken flow has been studied in Ref. [27].Developing a set of hydrokinetic equations, the first fractional power correction to the longitudinal pressure at late times was found as ∝ 1/(τ p T ) 3/2 (τ p is the proper time and 1/τ p the expansion rate of the Bjorken flow).If the expanding flow has a U (1) charge, the fluctuations also result in a fractional power correction of the U (1) current.This is the conse-quence of non-linear mode-coupling between the currents and the hydrodynamic variables [28].In contrast to previous studies, our model also features a UV-regulator, i.e. τ .Without deriving hydro-kinetic equations, we estimate the effect of the UV-regulator on the late-time nonlinear correction to the single charge density, ∆⟨n(τ p )⟩, in "non-fluctuating" Bjorken flow.In contrast to [28], the effect comes from the self-interactions of n.We find (Supplemental Material) The leading correction is similar to Eq. (78a) in [28].In fact this term can be found through the work of [6,7].However, the sub-leading term is found entirely from the theory proposed in this letter, representing the effect of the UV-regulator.It is important but if τ ≪ τ p smaller than the leading non-linear effects.In order to make the above results more accurate in a more realistic scenario, it would be interesting to investigate the effect of the UV-regulator on the hydro-kinetic setup of [28].
Another aspect of our theory of non-linear diffusion, worthy of future exploration, is its application near a critical point, such as for example the critical point expected in the QCD phase diagram.Near such a critical point charge fluctuations, the correlation length, and the relaxation time of our UV-mode all become large [22], such that it has to be included into the long-wavelength description [23], making our non-linear theory of diffusion relevant to the search for the critical point [24][25][26].

CORRELATION FUNCTION IN CLASSICAL AND QUANTUM LIMITS
From the linear equation ( 2),one can easily compute the correlation functions between the fluctuations of the density.Following [1], the equations of one-point functions given in (1) can be promoted to equations for the correlation functions: Note that the above averages are over background thermodynamic fluctuations, the so-called noise fields.In the absence of correlation at non-zero distances, the equal-time correlation functions reduce to those being familiar from thermodynamics.The initial boundary condition is then the vanishing of the equal-time correlation function of the n and J together with Let us define Then by a one-sided Fourier transformation with respect to time and complete transformation with respect to spatial coordinates, equations ( 17) simplify to where ⟨nn⟩ Note that the above first equality is the quantum mechanical version of the relation between the correlation function and its one-sided Fourier transform [1].Since we want to add a UV regulator to the diffusion theory, ω can be of order of T and consequently the presence of such quantum effects are inevitable.In the classical limit ω ≪ T , and then the expression ω 2T coth ω 2T reduces to unity.Deviations from the classical limit may be written as Is the above Taylor expansion convergent?The answer is not, because the left hand side has a set of branch point singularities in the complex frequency plane at ω = ±i2πnT, n ̸ = 0.The domain of convergence of the above expansion is then limited to |ω| < 2πT .Then the next question is how much large the momentum cutoff Λ EFT can be considered in our setup?Based on the above discussion, in order to continue to use the Taylor expansion (21), we may take Λ EFT ≲ 2πT .However, the larger the value of Λ EFT , the slower the series in ( 21) converges.For this reason, we take Λ EFT ≃ 3T .In this way, the first three terms in the series are sufficient to obtain an accurate result in perturbation theory with less than 2% error.On the other hand, the gap mode ω 2 needs to belong to the spectrum; The latter demands τ T ≳ 1/3 which is satisfied by any τ obeying the relaxation time bound in strongly coupled systems, namely τ T ≳ 1. 1 .In summary, in the weak coupling limit, which corresponds to 1/(τ T ) ≪ 1 [3], (21) reduces to unity and quantum effects are unimportant.In the strong coupling limit, below Λ EFT ≲ 2πT , we can safely use (21) to include the quantum effects.Nonetheless, our calculations show that even the first two sub-leading terms do not contribute significantly to the expression of correlation function.For this reason, in the main text we took Q = 1 and represented the results for this case. 1 In addition to the conjectured bound mentioned in footnote 1, there is a stronger bound on τ in a relativistic deformed CFT in 1+1 dimension with finite number of degrees of freedom [2].The latter bound is the consequence of causality and KPZ universality.In the traditional language [1], the correlation function ( 20) is modeled as originating from some microscopic random (noisy) currents, say j.Then the noisy dynamics is governed by In this work we assume that these microscopic currents have Gaussian short-distance correlations2 : To determine the strength of the noise, we combine equations (22): Here θ = −∇ • j is a Gaussian noise field with Now, from the equation ( 23), we can simply find δθ 2 ωk in terms of δn 2 ωk then by using (20), and considering (24), we arrive at with σ = χD.Having specified the correlation function of the noise field, we now want to move on to constructing the effective action that reproduces all the above correlation functions.We will follow the method developed in [5] 3 .

THE EFFECTIVE ACTION
Equation (25) tells us that the average in (20) is in fact over the noise distribution.Then it is reasonable to define the noisy fields, n θ (t, x), with the subscript θ denoting that n θ is the solution to the equation (23).Taking W as the weight of the noise field distribution, we write where J = δ(e.o.m.)/δn is the Jacobian.The idea of finding the effective action is to exponentiate the above delta function and then integrate over the noise field θ.At the linearized level, i.e., when e.o.m. is linear, Jacobian is not field-dependent.By introducing an auxiliary field n a , we exponentiate the delta function as Dn a e i (e.o.m.)na J n(t Considering (23), the integrating over θ gives where From the above quadratic effective action we can read the quadratic Lagrangian; the free propagators then read and nn is given by in complete agreement with (20).Although G nn in this model is well-known from Kadanoff-Martin work [8], however, (32) is the first derivation of this result from an effective action (involving quantum effects).
Let us consider the non-linear equation of motion (5).Repeating the MSR procedure, the Jacobian in ( 27) is no longer field-independent in this case; however, we can exponentiate it by introducing the anticommuting ghost fields ψ and ψ: Then the interacting part of the Lagrangian, up to quartic order, is given by ( 6).Although we assume that the noise is Gaussian (see (24)), but since we expand its magnitude in terms of n, i.e. σ(n) = σ + λ σ δn + 1 2 λ ′ σ δn 2 , we are able to produce the interaction terms nn 2 a and n 2 n 2 a , consistent with the result of Schwinger-Keldysh framework [6].We would also like to point out that we did not apply any constraints from the fluctuation-dissipation theorem to the interaction terms in the action.This is indeed beyond the scope of the MSR formalism.Here we only apply the fluctuation dissipation theorem to the two-point functions in quadraic part of the cation.A detailed analysis of KMS constraints for higher n-point functions, or the so-called generalized fluctuation dissipation theorem, can be found in [9].See also [10] for an interesting discussion of T −reversal symmetry in a system with self-interacting stochastic fields.This symmetry is closely related to the KMS symmetry.
As the last point, note that, as discussed in [5] the ghost terms do not contribute to the following computations.

LOOP CALCULATIONS
In order to study the effect of hydrodynamic interactions on the correlation functions, we parameterize the 1-loop corrections to G nna as the following (p = (ω, k)): with Σ being the self-energy appearing in the retarded Green's function (7).It turns out that there are only two one-loop diagrams contributing to the self-energy In fact, Σ may correct D and τ , as we will discuss below.Then we find with Q(ω) defined in (21).In order to evaluate the frequency integrals, we need to know the analytic structure of the integrand.Each of Green's functions has only two or four simple poles.The expression Q(ω ′ ), however, has infinite number of poles on the imaginary axis in the complex frequency plane.As discussed earlier, as long as the system under study satisfies Λ EFT ≲ 2πT , we can safely expand it about ω ′ = 0. Let us emphasize that in this way we will include all quantum effects in addition to statistical fluctuations.Taking a hard momentum cutoff, performing standard although lengthy calculations, the cutoff-independent part of Σ is found to be given by ( 9) and (10).At d = 1 the functions f 1 and f 2 in (9) are given by and We do not represent the structure of f 1 and f 2 at d = 2, 3 here.As was also mentioned earlier, 2 are suppressed by smaller numerical factors compared to the leading expression.
So far we have found Σ in (7).It is easy to show that the pole singularities of the retarded Green's function (7) are the roots of the following equation4 To find δσ(p), we apply the FDT theorem to (7) and first obtain G nn as follows Here the numerator N (p) contains δσ(p) as follows We then need to calculate N (p).Diagrammatically, we have which leads to (40) Once we find N (p) in (40), we can read δσ(p) from (39).It turns out that to leading order in the expansion of C in (21), f 3 vanishes.The expression of f 4 for d = 1 is given by: Theoretically, it would be interesting to realize the results of this letter in the framework of Schwinger-Keldysh EFT [4,11,12].To our knowledge, inclusion of gapped modes beyond quadratic order [13,14] has not been explored in this context yet.

THRESHOLD SINGULARITIES
Considering the branch points given in Eq. ( 14), ω11 and ω22 correspond to the minimum energy to generate a pair of on-shell ω 1 and ω 2 excitations in the loop, respectively.For instance, in the case of ω11 , we can simply consider the situation in which the two lines of the second loop in (8) carry the excitation ω 1 (k).Conservation of energy and momentum then enforces 3), we find that ω(k) becomes the minimum value at k ′ = k/2, and this value is equal to ω11 in (14).

DISPERSION RELATIONS
From (12), we can also extract the dispersion relation of fluctuations.Defining the dimensionless effective coupling constant and momentum as (Dτ ) d/4 λ D , we have numerically illustrated the spectrum of fluctuations in Fig. 4 for d = 1.For d = 2, 3, the spectrum is qualitatively similar.We find that due to the self-interactions, first, each of the two modes ω 1 and ω 2 splits into two modes.Second, the split dispersion relations avoid a level-crossing; the stronger the self-interactions, the more repulsion between them.
-��� -��� -��� . Spectrum of the excitations in the complex w = ωτ plane.Each colored trajectory, starting with dark red and ending with purple, shows the change of a particular mode when q = (Dk 2 τ ) 1/2 increases from 0 to 0.8.The two modes of the non-interacting theory given by (3) are shown in low opacity.The dashed semicircle is schematically showing the cutoff of Schwinger-Keldysh EFT of diffusive fluctuations discussed in Refs.[6,7].
. Left panel: The dominant part of the analytic structure of Gnn(ω, k) at late times (t >> tD).Right pane: Changing the integral variable from ω to ω and contour deformation.Note that iω11 = 1 − 1 − τ Dk 2 is real-valued.

LONG-TIME-TAILS
In this section we find the long time tail behavior of G R (1) nn and G nn .Let us start by expanding (12).We have Now by using the fluctuation dissipation theorem, G nn = 2T ω ImG R nn , we find To perform the Fourier integral G nn (t, k) = dω 2π G nn (ω, k) e −iωt at t > 0 we consider the analytic structure of G nn (ω, k) in the lower half ω plane.It is shown in Fig. 1.In the asymptotic limit the dominant contribution comes from the branch point closest to the real axis, namely ω11 .Thus to find the long time tail, or equivalently the asymptotic behavior of G nn (t, k), we need to expand G nn (ω, k) about ω11 and keep only the leading term.For this leading term, the branch cut structure is given in the left panel of Fig. 5.This is nothing but the the region surrounded by the dashed semicircle in figure 1.Now taking ω = ω11 + iΩ, we can deform the contour of integration, as depicted in te right panel of Fig. 5.The Fourier integral takes the following form (at d = 1) where g = The integral in the above expression evaluates to 1/ √ πt and we find the result given in (15).

THE EFFECT OF UV REGULATOR ON FLUCTUATIONS IN EXPANDING QGP
Let us recall that the linear response of the system to external sources can be found through the linear response framework.For the system under the consideration in this letter, we may consider a chemical potential source at t < 0 that is turned off at t = 0: µ(t, x) = e ϵt µ(x)θ(−t) [5].Then nn is the regarded Green's function in the absence of non-linear interactions.We find where µ 0 (k) is the Fourier transform of µ(t, x) at t = 0. We would like to use these expressions for the case of Bjorken flow.Since the flow profile only depends on the proper time τ p , we need to have a dynamic quantity that has only a time-dependence.For this reason, we choose to introduce ⟨n(t)⟩ ≡ ⟨n(t, x = 0)⟩.It is given by To find this, two things need to be determined in (47); G (0)R nn (ω, k) and µ 0 (k).Let us emphasize that our final goal is to find the correction to ⟨n(t)⟩ caused by the non-linear interactions.Let us call it ∆⟨n(t)⟩.For this, we evaluate the integral in (47) when , k).The latter is found from our EFT calculations To specify µ 0 (k) we take it into account that in the presence of an external deriving frequency ω ext , there will be an important length scale in the system, the so-called dissipative scale k * Modes with k ≫ k * have been already equilibrated.Modes with k ≪ k * are out of equilibrium and evolve according to linear hydrodynamics.However, for modes k ∼ k * , the equilibrium is ongoing.These modes contribute to the hydrodynamic fluctuations through non-linearities [15].To exclusively estimate the contribution of these modes, we take the chemical potential to have the following form where τ ext ∼ ω −1 ext and µ is the amplitude of µ(t, x) at t = 0.The normaliztion factor has been chosen so that Now we have the necessary ingredients, namely (49) and (51), to calculate (48).We perform the calculations in Bjorken flow.In this case, ω ext is nothing but the expansion rate of flow, i.e.; 1/τ p with τ p being the proper time.The last step is to identify t with τ p and to evaluate the integral in (48).The result is The timescale set by the UV-regulator, namely τ is much smaller than the proper time τ p : τ ≪ τ p .See equation (12) in [15].In our notation, it is expressed as Then it is evident that τ p ≫ τ .To gain more insight into (52), we expand it in powers of τ τp : A more realistic scenario should include not only diffusion.This requires improving the present work by considering the full MIS equations [16][17][18], for recent presentations see [3,19,20].It would be interesting to construct the hydro-kinetic equations for this full set of MIS equations.

FITTING METHOD
To evaluate the possibility of non-linear effects in the bad metallic system examined in [21], the authors of [21] fit the experimental results to δn(t, k ) (6 different values of k) to the analytical solution of equation ( 2).The solution is parameterized by three parameters: Γ, D and the amplitude δn(0, k).For any set of the data points associated with a specific value of k, they find these parameters.The results are given in plots S1-B and S1-C of this reference.
Our goal is to fit the data with a solution to the non-linear equation ( 5) (λ ′ D = 0).However, this equation has no analytical solution.Therefore, we must fit the data to its numerical solution.After trying different fitting Methods in Mathematica, we came to the conclusion that in order to find reliable fitting results, 1. Linear and non-linear differential equations must be solved in the same way.
2. Both linear and non-linear fits must be done using the same method.
The first item above means that both linear and non-linear equations should be solved numerically.This is why the linear fit (orange and red) points in Figure 3 do not precisely coincide with the results from [21].Reference [21] chooses to fit the data to an analytical solution of (2), whereas we do this by fitting a numerical solution of (2).We have used NDSolve[... , Method-> "StiffnessSwitching"] The solution of a linear equation at any k is parameterized by three parameters Γ, D and A ≡ δn(0, k).However, the solution to a nonlinear equation is specified by four parameters: Γ, D, A, and λ.Therefore we call the solutions of linear and nonlinear equations solLinear[GammaL, DL, AL][t] and solNonLinear[GammaN, DN, AN, lambdaN][t], respectively.For the second item above, we found that the appropriate method is NMinimize.For the case of linear fitting, we simply do it as follows  However, for the case of fitting data to the non-linear equation, the situation is more complicated.We were unable to find any non-linear smooth fit to the original data.Instead, we find that at any k it is better to use the solution of the linear equation, i.e. solLinear[GammaL, DL, AL][t] as the benchmark.We then generate a new set of data from it.The number of new data points read from this function depends on k.We define  Results are given in Table .II.
Let us comment on the fitting method and its impact on the results.As mentioned before, in order to advance linear and nonlinear fits consistently, we choose to use a common Method for both fits.Of the existing Method's in Mathematica, we found only one useful: NMinimize.For some other methods, such as Gradient and Newton, we can only fit the data with numerical solutions to "linear" equations.For the LevenbergMarquardt method we cannot even find any fit to the linear equation.).The black dashed curve corresponds to the fitting parameters given in the first row of Table I.
In Fig. 6 we show the results of a linear fit to data with a specific value k (or equivalent λ) found by three different methods.We see that while Mathematica is able to produce convergent results for all three cases.however, this is only NMinimize that can be considered a suitable method.It is worth noting that for some data sets from [21], we find the results of the Newton and Gradient to be divergent functions of time.But the results of NMinimize method, given in Table.I, always converged and performed well for all eight sets of data from [21].With all the above points and observations, we are led to use the NMinimize method.It should be noted that it is not surprising to find that only one of the above methods is suitable for our problem.In any fitting problem related to the solution of differential equations, the appropriate choice of fitting method depends largely on several factors: • The fitted equation and whether the equation has an analytical solution.
• The number of parameters in the equation plus the number of parameters in the boundary conditions.
• The number of data points in any data set and their distribution.
Different fitting methods often do not give the same results for a specific problem.

APPLICATION NEAR A CRITICAL POINT
Fluctuations in the UV-regulated theory of diffusion are also important near a critical point.Let us suppose a simple diffusive charge near the critical point.In the dynamic model H [22], with ξ being the correlation length.It turns out that τ ∼ ξ 2 [23].Near the critical point ξ becomes large, so does the relaxation time τ .As a result, the UV mode of our theory becomes a slow mode and then has to be included even in the long-wavelength limit [24].On the other hand, we find that Σ 3 ∼ ξ 9/2 , becoming large near the critical point.This shows the need to consider non-linearities discussed in this paper.Within a QGP droplet passing by the conjectured QCD critical point, diffusion is coupled to bulk dynamics [25].The model discussed in this paper then needs to be coupled with relativistic hydrodynamic fluctuations.The effect of interactions and self-interactions in non-linear MIS theory on the final proton multiplicity cumulants should then be investigated [26], and is relevant to the search for the critical point [27].

2 ∞Figure 1 .
Figure 1.Contour plot of the phase angle, φ, of the complexvalued α1 = |α1|e iφ .Any discontinuous transition between two distinct colors represents a branch-cut of the retreaded Green's function in the complex w = ω τ plane.For comparison, we have schematically illustrated the range accessible to the diffusive theory of Ref. [6] by a dashed semi-circle in the top panel, while our theory of non-linear diffusion (UVregulated) with a physical gapped mode is valid in the entire range displayed.

Figure 2 .
Figure 2. Solid blue and green curves show G (0)+(1) nn while dashed curves correspond to G (0) nn .With blue and green, we illustrate distinct values of r = τUV /τD: 1/100 and 1/10, respectively.The red and yellow curves show the long-time-tail of G (0)+(1) nn for these two values of r.

Figure 3 .
Figure 3. Variation of fitting parameters Γ (top), D (bottom)versus the amplitude of the solution, namely n(0, k).Orange and red dots are obtained by fitting the data of[11]  to Eq. (2) while blue dots stem from fitting the data to Eq. (5).The orange, red, and blue dashed lines are the corresponding linear fits.Γ and D are normalized by the extrapolated zeroamplitude values Γ0 and D0, respectively.Error bars indicate the standard error of the mean value resulting from our fits.

Figure 6 .
Figure 6.Comparing three different linear fitting methods to the data at a specific k (from Ref.[21]).The black dashed curve corresponds to the fitting parameters given in the first row of TableI.

Table I .
Table.I. Linear fitting results with the NMinimize method.

Table II .
Nonlinear fitting results.