key: cord-0729047-0xg2cfry
authors: Bettelheim, Eldad; Smith, Naftali R.; Meerson, Baruch
title: Inverse Scattering Method Solves the Problem of Full Statistics of Nonstationary Heat Transfer in the Kipnis-Marchioro-Presutti Model
date: 2021-12-05
journal: Phys Rev Lett
DOI: 10.1103/physrevlett.128.130602
sha: 5342d12ed9ef6de2d5889903dae65c4a82b24880
doc_id: 729047
cord_uid: 0xg2cfry

We determine the full statistics of nonstationary heat transfer in the Kipnis-Marchioro-Presutti lattice gas model at long times by uncovering and exploiting complete integrability of the underlying equations of the macroscopic fluctuation theory. These equations are closely related to the derivative nonlinear Schr"{o}dinger equation (DNLS), and we solve them by the Zakharov-Shabat inverse scattering method (ISM) adapted by Kaup and Newell (1978) for the DNLS. We obtain explicit results for the exact large deviation function of the transferred heat for an initially localized heat pulse, where we uncover a nontrivial symmetry relation.

Introduction. -Full statistics of currents of matter or energy in macroscopic systems away from thermodynamic equilibrium is a fundamental quantity that has attracted much attention from statistical physicists in the past two decades. Major progress has been achieved in determining this quantity for nonequilibrium steady states in simple models of interacting particles [1] [2] [3] [4] . Nonstationary fluctuations of current, however, proved to be much harder for analysis [5] [6] [7] [8] [9] [10] [11] .

A convenient and widely used family of models for studying the full statistics of currents is stochastic lattice gases [12] [13] [14] [15] . One important example is the Kipnis-Marchioro-Presutti (KMP) model of heat transfer. The KMP model involves immobile particles occupying a whole lattice and carrying continuous amounts of energy. At each random move the total energy of a randomly chosen pair of nearest neighbors is randomly redistributed among them according to uniform distribution. The KMP model originally attracted much interest as the first model for which Fourier's law of heat diffusion at a coarse-grained level was proven rigorously [16] . By now it has become a paradigmatic model of nonequilibrium fluctuations of transport [4, 6-8, 11, 17-27] .

Here we study a full non-stationary heat-transfer statistics in the KMP model on an infinite onedimensional lattice. Suppose that only one particle has a nonzero energy at t = 0. Due to the energy exchange with the neighbors, the energy will start spreading throughout the system. At times much longer than the inverse rate of the energy exchange between the two neighbors (equal to 1/2), and at distances much larger than the lattice constant (equal to 1), the mean coarse-grained temperaturē u(x, t) in the KMP model is governed by the heat diffusion equation [12, 14, 16] ∂ tū (x, t) = ∂ 2 xū (x, t). The initial temperature is a delta-function,ū(x, t = 0) = W δ(x), and so the solution is u(x, t) = (W/ √ 4πt) exp(−x 2 /4t) .

However, in stochastic realizations of the KMP model the coarse-grained temperature will fluctuate around the expected profileū(x, t), see Fig. 1 . To characterize these non-stationary fluctuations, we will consider the total amount of heat W > , observed on the right half line x > 0 at time t = T 1. The expected value of W > is W/2, and we will study the full time-dependent statistics of the heat excess, J = Similar non-stationary large-deviation settings, but with a step-like initial condition for the particle density or temperature, have been recently studied for a whole family of diffusive lattice gases [6] [7] [8] [9] [10] , of which the KMP model is an important particular case. The main working tool of these studies has been the macroscopic fluctuation theory (MFT) [28] : a weak-noise theory, whose starting point is fluctuational hydrodynamics (FH) [12, 14, 29] . The FH is a coarse-grained description of the lattice gas, which is accurate when the characteristic length scale of the problem (here the diffusion length √ T ) and the observation time T are much larger than the lattice constant 1 and the inverse elemental rate 1/2 of the energy exchange, respectively. For diffusive lattice gases with a single conservation law the FH has the form of a single macroscopic Langevin equation, which accounts for the fluctuational contribution to the heat or mass flux. For the KMP model the Langevin equation reads [12, 14] 

where u(x, t) is the temperature, and η(x, t) is a delta-correlated Gaussian noise:

The MFT [28] relies on a saddle-point evaluation of the path integral for the stochastic process, described by Eq. (2) . The small parameter of the saddle-point evaluation is again 1/ √ T 1: long times correspond to a weak noise. The saddle-point evaluation of the path integral boils down to a minimization of the action functional [30] , constrained by the specified heat excess J at t = T and obeying the specified initial condition u(x, t = 0). For the statistics of the heat (or mass) excess, the MFT equations and boundary conditions in time were derived in Ref. [6] , and we will present them shortly. For completeness, we also present their derivation in [30] . The solution of the MFT problem describes the optimal path of the process: the most likely time history of the temperature field u(x, t) which dominates the probability distribution P(J, T ) that we are after. The MFT problem, however, has proven to be very hard to solve analytically, especially for quenched (that is, deterministic) initial conditions [31] . In particular, for the KMP model, only small-J [7] and large-J [8] asymptotes have been obtained until now (but for a step-like initial condition). This Letter reports a major advance in this area of statistical mechanics. We present an exact solution to the heat excess statistics problem by uncovering and exploiting complete integrability of the underlying MFT equations. We obtain explicit results for an initially localized heat pulse, u(x, t = 0) = W δ(x), for which we uncover a nontrivial time-reversal mirror symmetry. These are the first exact non-steady-state large-deviation results for the statistics of current in a lattice gas of interacting particles for quenched initial conditions.

Formulation of the MFT problem [6, 30] . -Let us rescale t, x and u by T , √ T and W/ √ T , respectively. The optimal path we are after is described by two coupled Hamilton's equations for the rescaled temperature field u(x, t) and the conjugate "momentum density" field p(x, t) which describes the optimal history of the noise η(x, t), conditioned on the heat excess J.

It is convenient to introduce the (minus) gradient field v(x, t) = −∂ x p(x, t). In the variables u and v, the MFT equations take the form [6, 8, 30] 

The rescaled initial condition is

The condition on the heat excess at t = T becomes

The minimization of the action functional, that enters the constrained path integral, with respect to variations of u(x, t) yields, aside from Eqs. (3) and (4), a second boundary condition in time [6] ,

where λ plays the role of a Lagrange multiplier, to be ultimately fixed by the constraint (6) .

Once u(x, t) and v(x, t) are found, one can calculate the rescaled action, which can be written as [6] [7] [8] 

The action yields the probability density P(J, T, W ) up to a pre-exponent:

Since √ T 1, Eq. (9) has a clear large-deviation structure, and the action s plays the role of a rate function.

A crucial and previously unappreciated observation is that Eqs. (3) and (4) coincides with the derivative nonlinear Schrödinger (DNLS) equation in imaginary time and space [32] . The DNLS equation (with real time and space) describes propagation of nonlinear electromagnetic waves in plasmas and other media [33] . An initial-value problem for the DNLS equation is completely integrable via the Zakharov-Shabat inverse scattering method (ISM) adapted by Kaup and Newell for the DNLS [33] . The MFT formulation presents an difficulty, however, as here one needs to solve a boundaryvalue problem in time, rather than an initial-value problem. Here we overcome this difficulty by (i) making use of a shortcut that allows one to determine the rate function s(j) even without the knowledge of u(x, t) and v(x, t) for all t, and (ii) exploiting a previously unknown symmetry relation [34] , specific to the initial condition (5):

Solution of the MFT problem.-Equations (3) and (4) belong to a class of integrable systems for which a Lax pair exists, i.e., as we explain below, the equations are equivalent to the compatibility condition of a system of two linear differential equations. The latter system defines scattering amplitudes which depend on u and v. The idea behind the approach that we shall use -the ISM -is to consider the time evolution of these scattering amplitudes, which turns out to be very simple, as shown below. By relating these scattering amplitudes, at t = 0 and t = 1, to the fields u and v, the method will enable us to find the heat excess j = j(λ) which suffices for the calculation of s = s(j).

Adapting the derivation of Kaup and Newell [33] to imaginary time and space, we consider the linear system

where ψ(x, t, k) is a column vector of dimension 2,

and k is a spectral parameter. As one can check, the com-

is indeed equivalent to Eqs. (3) and (4). Let us define the matrix T (x, y, t, k) as the xpropagator of the system (11), namely, the solution to

with T (x, x, t, k) = I (the identity matrix). At x → ±∞, where the fields u(x, t) and v(x, t) vanish, the matrix U becomes very simple,

Therefore, it is natural to define the full-space propagator G(t, k) as follows:

The entries of the matrix G(t, k) are the scattering amplitudes of the system (11) . The time evolution of G(t, k) is easy to find. Indeed, the matrix T (x, y, t, k) satisfies:

One can check that Eq. (17) is compatible with (14) 

too becomes very simple in the limit x → ±∞,

Plugging (18) into (17), one finds the time evolution of T (x → ∞, y → −∞, t, k) which in turn, using (16) , yields that of G(t, k) :

where we have introduced here a notation for the matrix elements of G(t, k).

Plugging the temporal boundary conditions (5) and (7), we calculate G(0, k) and G(1, k) explicitly by solving Eq. (14), see [30] . Comparing the two solutions and using (19) we obtain

are the Fourier transforms of v(z, 0) restricted to z > 0 and z < 0, respectively:

We solve Eq. (20) in [30] , with the result

where v ± = v(0 ± , 0).

To compute v ± , we demand that Q ± (k) be regular at the origin, corresponding to a vanishing v(z, 0) at infinity. Setting k = 0 in Eq. (22) and using the Sokhotski-Plemelj formula

we obtain after some algebra

Taking the derivative of Eq. (22) with respect to k at k = 0, yields [30] Q Figure 2 shows Re Q + (k) and Im Q + (k) versus k at λ = 1, obtained by plugging Eq. (25) for v ± into Eq. (22) . This figure also shows the same quantities computed by solving Eqs. (3) and (4) numerically with a back-and-forth iteration algorithm [35] . The analytical and numerical curves are almost indistinguishable.

Analytical results for Q+(k), described by Eqs. (22), (23) and (25) Using Eqs. (6), (10) and (21) alongside with the conservation law

(27) Now we use a shortcut which makes the results we have obtained so far sufficient for obtaining the rate function s = s(j). The shortcut comes in the form of the relation ds/dj = λ, which follows from the fact that j and λ are conjugate variables, see e.g. Ref. [36] . It allows one to calculate s(j) bypassing Eq. (8) [which would require the knowledge of the whole optimal path u(x, t)]. We have

Using Eq. (26), we integrate Eq. (28) with respect to λ to get

where Li 2 (z) = ∞ k=1 z k /k 2 is the dilogarithm function, Q + (0) is given by Eq. (26) , and the integration constant was determined from s(λ = 0) = 0. Equations (27) and (29) give the complete rate function s(j) in a parametric form and represent the main result of this work. The exact optimal history of the temperature profile u(x, t) proved difficult to obtain analytically, but it can be computed numerically [30] . Figure 3 shows s(j) alongside with two asymptotes: j → 0 and |j| → 1/2, which correspond to λ → 0 and |λ| → ∞, respectively. Also shown are results of Monte-Carlo simulations. The asymptote λ → 0 can be obtained either from the exact rate function (27) and (29) [30] , or from a perturbative expansion applied directly to the MFT equations [7] . By virtue of the symmetry (10), the latter can be done very easily. Indeed, in the leading order in λ 1 Eqs. (8) , (10) and (1) yield

The shortcut relation ds/dj = λ can be rewritten as (ds/dλ)(dλ/dj) = λ. Combined with Eq. (30) it yields s(j → 0) √ 8πj 2 . Then, from Eq. (9), we see that typical fluctuations of J are normally distributed with variance W 2 /(32πT ) 1/2 . The T −1/2 scaling of the variance should be contrasted with the T 1/2 scaling, obtained for a step-like initial condition [5] [6] [7] . However, the relative magnitude of the fluctuations -the ratio of the standard deviation and the average transferred heathas the same scaling T −1/4 1 in both settings, as to be expected from the law of large numbers. The asymptote of |λ| → ∞ is more subtle [30] . The final result, already in terms of j, is

where ∆ ≡ 1/2 − |j| 1, and W −1 (. . . ) is the proper branch of the product log (Lambert W ) function [37]. At j = 1/2, s diverges and P vanishes, as to be expected. Nested-log large-current asymptotes similar to Eq. (31) appear to be typical for the KMP model and other models of the hyperbolic universality class [8, 9, 11] .

Discussion.-By combining the MFT and the ISM, we calculated exactly the rate function s(j), see Eqs. (27) and (29) , which describes the full long-time statistics of nonstationary heat transfer in the KMP model for an initially localized heat pulse. This is the first exact non-steady-state large-deviation result for the statistics of current in a lattice gas of interacting particles for quenched initial conditions. It opens the way to extensions of the ISM to additional fluctuating quantities of the KMP model. Another challenging goal is to apply the ISM to the simple symmetric exclusion process (SSEP) [12] [13] [14] [15] -a lattice-gas model with quite different properties [9] . Encouragingly, the MFT equations for the SSEP (see e.g. [6] ) can be mapped to Eqs. (3) and (4) via a canonical transformation. This transformation, however, complicates the boundary conditions in time.

From a more general perspective, the MFT of lattice gases is a particular case of the weak-noise theory, or optimal fluctuation method (OFM): a highly versatile framework which captures a broad class of large deviations in macroscopic systems. For non-stationary processes the OFM equations -coupled nonlinear partial differential equations for the optimal path -are usually very hard to solve exactly. One class of problems of this type, which has received much recent attention, deals with the complete one-point height statistics of an interface whose dynamics is described by the Kardar-Parisi-Zhang (KPZ) equation [38] . The OFM captures the complete KPZ height statistics at short times [39] [40] [41] [42] [43] [44] .

Here too, a previous analytical progress in the solution of the OFM equations was limited to asymptotics of very large or very small interface height. But very recently these OFM equations -which coincide with the Nonlinear Schrödinger equation (NLS) (not the derivative one) [42] -have been solved exactly [45, 46] by the ISM for several "standard" initial conditions. The two integrable systems, the NLS and DNLS, are closely related, so our approach can be compared with that of Refs. [45, 46] . We used only standard techniques of the ISM which do not rely on additional tools, such as Fredholm determinants used in Refs. [45, 46] . Because of its relative simplicity our approach appears to be more readily adaptable to solving additional large-deviation problems [47] .

Here we give some technical details of the calculations described in the main text of the Letter.

We begin by defining an auxiliary potential

Integrating Eq. (2) of the main text with respect to x and using the boundary condition u(x → −∞, t) → 0, we obtain a Langevin equation for ψ(x, t):

Now we use a path-integral approach, by writing the probability density functional P [η] of the white Gaussian noise term η(x, t):

We now express η through ψ via Eq. (S2), enabling us to write the probability density for a given history of the system:

where S [ψ] is the action functional. The presence of the large parameter √ T 1 enables us to calculate P(J, T ) via a saddle-point evaluation of the path integral. Within this framework, − ln P (J, T ) is given by minimum of the action functional S, constrained on the initial condition u(x, 0) = W δ(x) and on a given value of heat excess J. The heat excess, which we recall is J = ∞ 0 u(x, t = T ) dx − W/2, is conveniently rewritten as J = W/2 − ψ(0, T ) (where we used the conservation law ∞ −∞ u(x, t = T ) dx = W ). To take into account the latter constraint on J, we add the term −Λψ (0, T ) to the action, where Λ is a Lagrange multiplier, i.e., we minimize the constrained functional

The subsequent procedure is pretty standard. Consider a variation ψ(x, t) → ψ(x, t) + δψ(x, t). This leads to a variation of S Λ [ψ], which, to first order in δψ is given by

where we have defined the momentum density gradient

Integrating by parts in Eq. (S6), we obtain

where the first two terms in the single integral are the boundary terms originating from the integration by parts in time. For the quenched (deterministic) initial condition, the variation δψ (x, 0) vanishes. The second MFT equation, Eq. (4) in the main text is now obtained by requiring the double integral in (S8) to vanish for arbitrary δψ (recalling that v = −∂ x p and u = ∂ x ψ). The first MFT equation, Eq. (3) in the main text, follows from Eq. (S7) after multiplying by the denominator and then taking a spatial derivative. Note that, under the rescalings of x, t and u described in the text, v should be rescaled by 1/W . These rescalings leave the MFT equations invariant. Requiring the single integral in Eq. (S8) to vanish for arbitrary δψ (x, T ), we obtain the boundary condition

After the rescaling, this becomes Eq. (7) in the main text, where λ = W Λ is a rescaled Lagrange multiplier. Finally, using Eq. (S7) in (S4) we find that the action can be rewritten as S = Let us find the matrix T (x, y, 0, k) at t = 0. By solving Eq. (14) of the main text at t = 0, using u(x, 0) = δ(x), one gets

ike ik(x+y)/2 e ik(x−y)/2 ± ike ik(x+y)/2 I u (0, y) , xy < 0 ,

where

and in the second case in (S10), the sign ± is to be taken as the sign of y. Plugging (S10) into Eq. (16) of the main text, we compute G(0, k):

in terms of Q ± (k) which are defined in Eq. (21) in the main text, with Q(k) = Q + (k) + Q − (k). It is useful to compare this result to the one obtained at t = 1. Here we have v(x, 1) = −λδ(x). Similarly to the t = 0 case, one gets

where in the second case, the sign ± is to be taken as the opposite of the sign of y. Now compute G(1, k):

Comparing the upper-right elements of G(0, k) from Eq. (S12) and of G(1, k) from Eq. (S14), using Eq. (19) 

where v ± = v(0 ± , 0). This turns Eq. (20) of the main text into

which has the solution

provided the condition

is satisfied. Eq. (S18) is derived by noting that Q ± are analytic in the upper and lower half plane respectively and are well-behaved when k is allowed to reach infinity through the respective half-planes. We then use the well-known decomposition f (k) = f + (k) + f − (k) of a general function f (k) into functions analytic in the upper and lower halfplanes, f ± (k), respectively, given by

. This decomposition is applied to the logarithm of Eq. (S17). Plugging Eq. (S18) into Eq. (S16), we obtain the solution given in Eqs. (22) and (23) of the main text.

Taking the derivative of Eq. (22) with respect to k at k = 0, yields the equation:

The imaginary part of the integrand is an odd function of k and therefore does not contribute to the integral. Keeping only the real part and simplifying it, we obtain (26) of the main text, where we replaced the principal value integral by a regular integral since the integrand is regular at k = 0.

At small λ Eq. (26) of the main text yields

therefore j(λ 1) λ/(4 √ 2π). Now, using the |z| 1 asymptotic Li 2 (−z) −z in the integrand of Eq. (29) of the main text, we obtain after a simple algebra: s(λ 1)

). This yields the asymptotic behavior s(j 1) 2 √ 2π j 2 given in the main text. Now we consider the |λ| 1 asymptote and start from Eq. (26) of the main text. Due to the symmetry j(−λ) = −j(λ), s(−λ) = s(λ) we consider only λ > 0. Let us denote the integrand in Eq. (26) (including the factor 1/4π) by

We can recast it as

where lambda, k) = ln 1+λ 2 k 2 e −2k 2

The integral over the first term of Eq. (S22) yields λ/2, and we now focus on the integral of lambda, k). For λ 1, lambda, k) as a function of k behaves as follows (see Fig. 4 ). At 1 λk < √ ln λ, lambda, k) is approximately constant and equal to −1/(2π). This asymptote is obtained when neglecting 1 in the numerator and in the denominator of the fraction inside the logarithm in Eq. (S23). For k √ ln λ, Φ(λ, k) behaves as

This asymptote is obtained when neglecting the second term in the numerator and 1 in the denominator of the fraction inside the logarithm in Eq. (S23). Importantly, the transition from the asymptote −1/(2π) to the asymptote (S24) occurs in a narrow boundary layer around k = √ ln λ, whose width goes to zero as λ → ∞. Furthermore, the region of k 1/λ contributes a term O(1/λ) to the integral which, as we shall see a posteriori, is negligible. Therefore, we can divide the integration region into two subregions: 0 < k < √ ln λ and √ ln λ < k < ∞, and use the asymptote Φ(λ, k) −1/(2π) in the former subregion, and Eq. (S24) in the latter one. Keeping the leading and two subleading terms in the result and multiplying it by 2 to account for k < 0, we obtain Using the relation (27) of the main text between j and Q + (0), we obtain in the leading order:

To calculate the large-λ asymptote of s(λ), we plug (S26) into the first equality in (28) in the main text to get

Integrating this over λ we obtain

(Note that, in the limit λ 1 that we consider here, the integration constant is not important.) Solving Eq. (S26) for λ, we obtain

where ∆ ≡ 1/2 − j 1, and W −1 (. . . ) is the proper branch of the product log (Lambert W ) function [1] . Plugging Eq. (S29) into Eq. (S28), we obtain the large-j asymptote (31) of the main text.

Our exact solution gives the rate function s(j) but it does not give the optimal path u(x, t) at all times. To calculate the latter analytically, one would need to solve Eq. (14) of the main text at arbitrary times 0 ≤ t ≤ 1, and this is far more challenging than solving it only at t = 0 and t = 1, as we did. The optimal path can, however, be computed numerically, using the back-and-forth iteration algorithm due to Chernykh and Stepanov [2] . The results of one such calculation are shown in Fig. 5 . The optimal temperature profile u(x, t) for λ = 10 (corresponding to j 0.38) at times 1/4, 1/2, 3/4 and 1. Noticeable is a shock-like singularity of u at x = 0 and t = 1.

Here we describe how the numerical data in Fig. 3 of the main text was generated. We performed Monte-Carlo (MC) simulations of the microscopic KMP model. The dynamical object in the simulations is the vector of energies at 2L + 1 lattice sites, u −L , . . . , u −1 , u 0 , u 1 , . . . u L , where the site 0 represents the origin, x = 0. Initially, all of the energy is at the origin, u i = δ i,0 (this corresponds to W = 1). At every Monte-Carlo step, we randomly choose one of the 2L adjacent pairs of lattice sites (i, i + 1), and randomly redistribute their total energy between them, i.e., u i , u i+1 →ũ i ,ũ i+1 whereũ i is sampled from a uniform distribution between 0 and u i +u i+1 , andũ i+1 = u i +u i+1 −ũ i . Since (in our convention) the microscopic rate of this elementary process is 2, the total number of steps in the simulation is randomly sampled from a Poisson distribution with mean 4T L. The excess heat was measured at t = T as J = u 0 /2 + 

where the probability density function P(J, T ) was computed from the histogram of the simulated J's, and V = 1/ √ 32πT is our prediction for the variance [see the paragraph following Eq. (30) of the main text]. The term √ 2πV in Eq. (S30) comes from the normalization of the central part of the distribution, rather than from a systematic calculation of the pre-exponential factor in P(J, T ) which is beyond the leading-order MFT. The normalization approximation of the pre-exponential factor introduces a small systematic error at moderate J which is discernible in Fig. 3 of the main text. The relative error, introduced by the pre-exponential factor, must go to zero as √ T goes to infinity.

As to be expected, direct MC simulations become prohibitively long for larger J and T . Special methods of large deviation sampling should be used for these purposes.

Large Scale Dynamics of Interacting Particles

Stochastic Interacting Systems: Contact, Voter, and Exclusion Processes

Scaling Limits of Interacting Particle Systems

A Kinetic View of Statistical Physics

Statistical Physics

See Supplemental Material

For an annealed step-like initial condition (a stepfunction with equilibrium density distributions with given average densities u− and u+ at x < 0 and x > 0), the full heat transfer statistics for the KMP model was found by Derrida and Gerschenfeld [6]. They achieved it by uncovering an exact mapping between the KMP model and the simple symmetric exclusion process

Encyclopedia of Integrable Systems (L. D. Landau Institute for Theoretical Physics

10) leaves Eqs. (3) and (4) and the boundary conditions (5) and (7) invariant