Full Statistics of Nonstationary Heat Transfer in the Kipnis-Marchioro-Presutti Model

We investigate non-stationary heat transfer in the Kipnis-Marchioro-Presutti (KMP) lattice gas model at long times in one dimension when starting from a localized heat distribution. At large scales this initial condition can be described as a delta-function, $u(x,t=0)=W \delta(x)$. We characterize the process by the heat, transferred to the right of a specified point $x=X$ by time $T$, $$ J=\int_X^\infty u(x,t=T)\,dx\,, $$ and study the full probability distribution $\mathcal{P}(J,X,T)$. The particular case of $X=0$ has been recently solved [Bettelheim \textit{et al}. Phys. Rev. Lett. \textbf{128}, 130602 (2022)]. At fixed $J$, the distribution $\mathcal{P}$ as a function of $X$ and $T$ has the same long-time dynamical scaling properties as the position of a tracer in a single-file diffusion. Here we evaluate $\mathcal{P}(J,X,T)$ by exploiting the recently uncovered complete integrability of the equations of the macroscopic fluctuation theory (MFT) for the KMP model and using the Zakharov-Shabat inverse scattering method. We also discuss asymptotics of $\mathcal{P}(J,X,T)$ which we extract from the exact solution, and also obtain by applying two different perturbation methods directly to the MFT equations.

u(x, t = T ) dx , and study the full probability distribution P(J, X, T ).The particular case of X = 0 has been recently solved [Bettelheim et al. Phys. Rev. Lett. 128, 130602 (2022)].At fixed J, the distribution P as a function of X and T has the same long-time dynamical scaling properties as the position of a tracer in a single-file diffusion.Here we evaluate P(J, X, T ) by exploiting the recently uncovered complete integrability of the equations of the macroscopic fluctuation theory (MFT) for the KMP model and using the Zakharov-Shabat inverse scattering method.We also discuss asymptotics of P(J, X, T ) which we extract from the exact solution and also obtain by applying two different perturbation methods directly to the MFT equations.

I. Introduction 2 II. Formulation of the MFT problem 5
III. Solution of the MFT problem 6 A. Scattering amplitudes and their dynamics 6 B. Solving the scattering problem at t = 0 and t = 1 7 C. Recovering the rate function s(J, X) from the scattering amplitudes 8 1. Full solution for X = 0, and for X > 0 with sufficiently small λ 8 2. Solution  One fundamental problem of statistical mechanics deals with full statistics of currents of matter or energy in classical many-body systems away from thermodynamic equilibrium.This challenging problem has continued to attract attention of workers in the past two decades.There has been a remarkable progress in determining the current statistics for nonequilibrium steady states in simple models of interacting particles [1][2][3][4].Nonstationary current fluctuations, however, proved to be much more difficult to capture [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].The Kipnis-Marchioro-Presutti (KMP) model [16] is an important member of this family of models.The KMP model is considered a prototypical model of energy transfer by diffusion: not in the least because, for this model, the Fourier's law of heat diffusion at a coarse-grained level was proven rigorously [16].The model's definition is quite simple: There are immobile agents occupying the whole lattice and carrying continuous amounts of a scalar quantity which can be viewed as energy.At each (exponentially distributed) random move the total energy of a randomly chosen pair of nearest neighbors is redistributed among them randomly according to uniform distribution, and the process continues.Because of its simplicity the KMP model has become a paradigmatic model in the studies of fluctuations, including large deviations, of heat transfer in classical macroscopic nonequilibrium systems [4, 6-8, 11, 17-27].
To set up our ideas, let us consider the KMP model on a one-dimensional lattice.Suppose that only one agent has energy W (or a few neighboring agents have total energy W ) at t = 0. Due to the interactions with the neighbors, the energy will spread in a random fashion throughout the system.Looking for universality, we will focus on times much longer than the inverse rate of the elemental energy exchange between the two neighbors (equal to 1/2), and on distances much larger than the lattice constant (equal to 1).As it was shown already in the original paper of KMP [16], the mean, or expected, coarse-grained temperature ū(x, t) obeys the heat diffusion equation The initial coarse-grained temperature profile is a delta-function, ū(x, t = 0) = W δ(x), so the mean-field solution for the coarse-grained temperature is However, in individual realizations of the stochastic process, the coarse-grained temperature u(x, t) will of course fluctuate around the mean-field solution.These fluctuations are exemplified by Fig. 1 which shows a snapshot, and a coarse-grained version of it, of a Monte Carlo simulation of the microscopic KMP model in this setting.To characterize fluctuations of the heat transfer, we have recently evaluated [27] the full probability distribution P 0 (J 0 , T ) of the total amount of heat J 0 , observed at time t = T 1 on the right half-line: Our work [27] continued a line of studies of large deviations of non-stationary mass or energy transfer in different diffusive lattice gas models [5][6][7][8][9][10].In parallel, there has been much interest in the large deviations of the position of a tracer (or a tagged particle) in single-file systems such as the simple symmetric exclusion process (SSEP) [26,[28][29][30][31][32][33][34].(The SSEP is another well-known diffusive lattice gas model [12,14,35,36], where each particle can hop with equal probabilities to one of the adjacent lattice sites, but a hop to a site occupied by another particle is forbidden.)However, prior to our work [27], the full statistics of the corresponding quantities were determined only for annealed initial conditions, which assume a "pre-thermalization" of the lattice gas (or of the regions of x < 0 and x > 0 thereof at different densities or temperatures) [5,6,26,31,32], see also the more recent Ref. [37].
The evaluation of the full distribution P 0 (J 0 , T ) for the quenched (that is deterministic) initial condition (2) was achieved in Ref. [27] by employing the macroscopic fluctuation theory (MFT) [22], which will be introduced shortly, and by uncovering and exploiting complete integrability of the MFT equations for the KMP model by the Zakharov-Shabat inverse scattering method [38,39].In this work we extend the approach of Ref. [27] to a more general characterization of the stochastic heat transfer in the KMP model.We specify an arbitrary point x = X and evaluate the full time-dependent probability distribution P(J, X, T ) of the total amount of heat, observed at time t = T 1 to the right of x = X.Because of the energy conservation in the KMP model, the distribution P(J, X, T ) must have compact support: 0 < J < W .The mean, or expected value of J at fixed X readily follows from the deterministic solution (3): where erfc z = 1 − erf z, and erf z is the error function.For X = 0 one obtains J = W/2.The main objective of this work is to determine the full statistics of fluctuations of J around the mean value J.The presence of a new parameter 0 < X < ∞ allows for a more detailed characterization of the system, as we show here.Importantly, when the position x = X is varied at fixed J, it becomes an analog of the position of a tracer in single-file systems such as the SSEP [26,31,32].Indeed, in a single-file system, the total number of particles to the right of the tracer position x = X does not change with time, and this can be mimicked by keeping J = const in the KMP model.As a result, as we discuss below, the probability distribution P(J, X, T ) as a function of X at J = const has the same dynamical scaling properties as the probability distribution of the position x = X of a tracer in single-file systems, such as the SSEP.
Like many previous studies of the statistics of energy transfer [6-10, 27, 37] and of the tracer position [29,30], we will employ the MFT [22]: a weak-noise theory based on a saddle-point evaluation of the exact path integral for the fluctuational hydrodynamics of diffusive lattice gases.The MFT is a variant of the optimal fluctuation method (OFM) which goes back to Refs.[40][41][42].The OFM has been also applied to turbulence [43], stochastic reactions [44,45], non-equilibrium surface growth and related models [46][47][48][49] and to many other systems.A full list of references on different applications of the OFM would exceed a hundred.
Fluctuational hydrodynamics of diffusive lattice gases [12,14] is a coarse-grained description of the lattice gas on length scales much larger than the lattice constant 1, and time scales much larger than the inverse rate 1/2 of the energy exchange between the neighboring agents.The corresponding small parameter in our problem is 1/ √ T 1 [27].For a general one-component diffusive lattice gas, the fluctuational hydrodynamics has the form of the stochastic partial differential equation [12,14] where u(x, t) is the coarse-grained density or temperature, D(u) is the diffusion (or heat diffusion) coefficient, σ(u) is the gas compressibility, and η(x, t) is a delta-correlated Gaussian noise, η(x, t) = 0 and η(x, t)η(x , t ) = δ(x − x )δ(t − t ).For the KMP model, where u(x, t) is the coarse-grained temperature, one obtains D = 1 and σ(u) = 2u 2 , and Eq. ( 7) becomes Without the noise term Eq. ( 8) becomes the heat diffusion equation (1).As mentioned above, the MFT relies on a saddle-point evaluation of the path integral for the stochastic process, described by Eq. (7).In our problem, the small parameter of the saddle-point evaluation is again 1/ √ T .The saddlepoint evaluation leads to a constrained optimization problem, where the (integral) constraint comes from conditioning the process on the specified transferred heat (5) subject to the initial condition (2).The integral constraint (5) is conveniently taken into account via a Lagrange multiplier [6].We will present a complete Hamiltonian formulation of the MFT problem shortly.The solution of this problem describes the optimal path of the process: the most likely time histories of the coarse-grained temperature field and of the noise field which dominate the probability distribution P(J, X, T ) that we are after.
The MFT, being stripped of unnecessary microscopic details, is an appealing tool for studying large deviations in lattice gases on large scales and at long times.However, the MFT equations -two coupled nonlinear partial differential equations -are usually very hard to solve exactly, without reliance on additional small or large parameters.Ref. [27] reported a first major advance in this area of statistical mechanics by finding an exact solution to the MFT problem for the distribution P(J, X, T ) in the particular case of X = 0.As already mentioned above, the exact solution came from uncovering and exploiting complete integrability of the MFT equations by the Zakharov-Shabat inverse scattering method.A second major advance in this area has been recently reported in Ref. [37].The authors employed the MFT for determining the full statistics of J defined in Eq. ( 5) for the SSEP with a step-like annealed initial condition.Earlier this problem was exactly solved (by a Bethe ansatz) in its microscopic formulation, and the long-time asymptotic was extracted, in Ref. [6].The authors of Ref. [37] reproduced this asymptotic by combining a generalization of the canonical Cole-Hopf transformation with the inverse scattering method.
The recent works [27] and [37] (and Refs.[50][51][52] where complete integrability was uncovered [50] and exploited [51,52] for the determination of full short-time statistics of the interface height as described by the Kardar-Parisi-Zhang (KPZ) equation [53]), have shown a great potential for the combination of optimal fluctuation methods, such as the MFT, with the inverse scattering method for obtaining exact solutions for different large-deviation statistics.The present work is a natural next step in exploring this potential.
In Sec.II we present the formulation of the MFT problem for P(J, X, T ) [6,8] and briefly discuss the scaling behavior of P(J, X, T ) and its relation with the universal scaling behavior of the tracer statistics in single-file processes.In Sec.III we solve the MFT problem by the inverse scattering method.In Sec.IV we discuss some asymptotic limits of P(J, X, T ).Our results are briefly summarized in Sec.V. Some of the technical details of the calculations are given in the Appendices.

II. FORMULATION OF THE MFT PROBLEM
We rescale the variables t, x and u by T , √ T and W/ √ T , respectively.As a result, X and J become rescaled by √ T and by W , respectively.The optimal path of the process is described by two coupled Hamilton's equations: partial differential equations for the rescaled temperature field u(x, t) and the conjugate "momentum density" field p(x, t) (the latter describes the optimal history of the noise η(x, t) conditioned on the transferred heat J at given X).
We also introduce the (minus) momentum density gradient field v(x, t) = −∂ x p(x, t).As was shown earlier [6,8,27], the MFT equations can be written as To make this paper more self-contained, a derivation of Eqs. ( 9) and ( 10) is given in Appendix A. The initial condition, after the rescaling, is The rescaled transferred heat at t = T is constrained by the equation The particular case X = 0 was solved in Ref. [27].Being interested in X = 0 we can, without loss of generality, set X > 0. The constrained optimization problem yields, aside from Eqs. ( 9) and ( 10), a second boundary condition in time [6], where λ is a Lagrange multiplier, to be ultimately determined from Eq. ( 12).This MFT formulation possesses a nontrivial symmetry which generalizes the symmetry found for the particular case X = 0 in Ref. [27].The symmetry ( 14) is special to the delta-function initial condition (11).Once u(x, t) and v(x, t) are determined, the rescaled mechanical action can be evaluated [6][7][8]: This action determines the probability density P(J, X, T ) up to a pre-exponential factor.In the original, dependent variables we have ln which describes the scaling behavior of P(J, X, T ).As we will see in the following, at fixed J/W and small X/ √ T , the rate function s is quadratic in X/ √ T , describing Gaussian fluctuations, where the standard deviation of X scales with time as T 1/4 , that is sub-diffusively.This scaling coincides with the one universally observed for the tracer position in single file diffusion processes, see Ref. [29,30] and references therein.This coincidence is of course not accidental.Indeed, the mathematical formulation of the MFT problem that we have just presented, is very similar to the MFT formulation for the statistics of the position of a tracer [29,30].The differences are in the particular u-dependence of the multiplicative noise term in the Langevin equation (8), and in the delta-function initial condition (11) (rather than a step-like one) that we adopted here.Neither of them alters the subdiffusive scaling T 1/4 of the standard deviation of X. Finite or large X/ √ T corresponds to non-Gaussian large deviations which are already non-universal, and which we will explore.At fixed X and small J/W , the rate function s is also quadratic and therefore describes Gaussian fluctuations around the mean, which are universal for a whole family of lattice gases.For a finite J/W (large deviations) the distribution is non-Gaussian.
It was observed [27] that Eqs. ( 9) and ( 10) coincide with the derivative nonlinear Shrödinger equation (DNLSE) in imaginary time and space [54].With real time and space, the DNLSE describes nonlinear electromagnetic wave propagation in plasmas and other systems [55].The initial-value problem for the DNLSE is completely integrable by using the Zakharov-Shabat inverse scattering method (ISM) adapted by Kaup and Newell for the DNLSE [55].The heat transfer statistics presents an additional difficulty, as it involves a boundary-value problem in time.Ref. [27] overcame this difficulty in the particular case X = 0 by (i) making use of a shortcut (see, e.g.Ref. [56]) which makes it possible to determine the rate function s(J, X) even in the absence of the complete solutions u(x, t) and v(x, t), and (ii) exploiting the symmetry relation (14).Here we extend the approach of Ref. [27] to an arbitrary X.

III. SOLUTION OF THE MFT PROBLEM
Let us outline the key ideas behind the approach that we shall use to solve the MFT problem: the ISM [38,39].A key property related to the integrability of Eqs. ( 9) and ( 10) is the existence of a corresponding Lax pair.This means than Eqs.( 9) and ( 10) are equivalent to the compatibility condition of a system of two linear differential equations defined by a pair of matrices.For the latter system, one can define scattering amplitudes which depend on u and v.
As explained below, one can then solve for the time evolution of the scattering amplitudes, which turns out to be very simple.By relating the fields u and v to the scattering amplitudes, at t = 0 and t = 1, the ISM will enable us to find the transferred heat J = J(λ, X) which suffices for the calculation of s = s(J, X).

A. Scattering amplitudes and their dynamics
Adapting the derivation of Kaup and Newell [55] to imaginary time and space, we consider the linear system where ψ(x, t, k) is a column vector of dimension 2, the matrices constitute the Lax pair, and k is a spectral parameter.It is now straightforward to check that the compatibility condition is indeed equivalent to Eqs. ( 9) and (10).
It is useful to define the matrix T (x, y, t, k) as the x-propagator of the system (17), namely, the solution to the equation with the initial condition T (x, x, t, k) = I (the identity matrix).At x → ±∞, the fields u(x, t) and v(x, t) vanish, and the matrix U takes the simple form, As a result, the entries of T (x, y, t, k) oscillate as a function of x and y with wave numbers whose magnitude is k/2 at x → ±∞, y → ±∞.Therefore, it is natural to define the full-space propagator G(t, k) as follows: The matrix G(t, k) plays an important role in the ISM, because its (u-and v-dependent) entries are the scattering amplitudes of the system (17).The time evolution of G(t, k) is remarkably simple.Indeed, one finds that the time evolution of the matrix T (x, y, t, k) satisfies: As one can verify directly, Eq. ( 24) is compatible with (21) (i.e.∂ t ∂ x T = ∂ x ∂ t T ) due to Eq. ( 20).The matrix V (x, t, k) too becomes very simple in the limit x → ±∞, Plugging Eq. ( 25) into (24), one finds the time evolution of T (x → ∞, y → −∞, t, k) which in turn, using (23), yields that of G(t, k) : In Eq. ( 26) we have introduced a notation for the matrix elements of G(t, k).
B. Solving the scattering problem at t = 0 and t = 1 Let us find the matrix T (x, y, 0, k) at t = 0.This is achieved by solving Eq. ( 21) at t = 0, which reads (using The general solution to Eq. ( 27) is found by solving it in the regions x < 0 and x > 0 and then matching the solutions at x = 0.Then, taking into account the initial condition T (x, x, t, k) = I, we obtain (see Appendix B for details) where and in the second case in (28), the sign ± is the same as the sign of y.We now compute G(0, k), by plugging Eq. ( 28) into (23), with the result: in terms of Q ± (k) which are the Fourier transforms of v(z, 0) restricted to z > 0 and z < 0, respectively, and of their sum: It is useful to compare this result to the one obtained at t = 1.Here we have v(x, 1) = −λδ(x − X), which enables us to solve Eq. ( 21) at t = 1.The solution is very similar to that at t = 0, and one gets where and in the second case in (32), the sign ± is to be taken as the opposite of the sign of y − X. Similarly to the t = 0 case, we now compute where C. Recovering the rate function s(J, X) from the scattering amplitudes In this section we calculate the rate function s(J, X) by relating the scattering amplitudes at t = 0 and t = 1.We find that several cases must be considered.
(i) In the simplest case, in which λ is sufficiently small, the derivation closely follows that of [27].This case completely covers the case X = 0. (ii) For X = 0, and larger values of λ, technical problems are encountered when evaluating some of the integrals that appear in the solution in case (i).These are fairly simple to circumvent, by deforming the integration contours in the complex plane.This completely covers the case X < √ 8. (iii) For X > √ 8 a third regime, of very large λ, exists.This regime is rather interesting as it involves multiple solutions at a given λ.That is, there is a regime of λ's for which s and J are multi-valued functions of λ.Nevertheless, s remains a single-valued, analytic function of J, unlike in many other systems in which additional branches of the rate function lead to singularities which have the character of dynamical phase transitions [11,23,50,[57][58][59][60][61][62].As we explain below, this rather unusual behavior comes about concurrently with a non-convexity of s(J, X) as a function of J for X > √ 8.
1. Full solution for X = 0, and for X > 0 with sufficiently small λ Comparing the upper-right elements of G(0, k) from Eq. ( 30) and of G(1, k) from Eq. ( 34), using Eq. ( 26), we obtain Eq. ( 36) provides a key step towards the solution to our problem because, after solving it, J = J(λ, X) is obtained from Q + (0) with the use of the symmetry (14), and this, as shown below, suffices to determine of the rate function s(J, X).We now turn to the solution to Eq. (36).We begin by completing the squares in Eq. (36), by writing where the constant −Φ ± (0) in the exponent is found by requiring that Q ± (k) is regular at the origin.This equation turns Eq. ( 36) into which has the solution Eq. ( 39) is derived by noting that Q ± (k) are analytic in the upper and lower half complex k plane, respectively, and are well-behaved when k is allowed to reach infinity through the respective half-planes1 .We then use the wellknown decomposition f (k) = f + (k) + f − (k) of a general function f (k) into functions analytic in the upper and lower half-planes, f ± (k), respectively, given by The decomposition ( 40) is applied to the logarithm of Eq. ( 38), leading to Eq. (39).
Taking the derivative of Eq. ( 37) with respect to k, setting k to 0 and dividing by i we obtain: which, using the Sokhotski-Plemelj formula, can be also written in terms of a principal value integral, Using Eqs. ( 12), ( 14) and ( 31) alongside with the conservation law which can also be written as In order to proceed, we utilize a useful shortcut, which enables us to obtain the rate function s = s(J, X) given that we have calculated J = J(λ, X).The shortcut makes use of the relation ∂ J s = λ, a property that follows from the fact that J and λ are conjugate variables, see e.g.Ref. [56].This enables us to calculate s(J, X) bypassing Eq. (15) [which would require us to know the entire optimal path u(x, t)].Using J (λ, X) = 1 + Q + (0) /λ from Eq. ( 43), we have Using Eq. ( 42), we integrate Eq. ( 45) with respect to λ to get where Li is given by Eq. ( 42), and the integration constant was determined from s(λ = 0, X) = 0. Equations ( 43) and (46) give the rate function s(J, X) in a parametric form.37) and (39) (solid lines), versus numerical results (dashed lines), for λ = X = 1, corresonding to J = 0.31 . . . . Figure 2 shows s(J, X) alongside with the asymptotic J → J, which corresponds to λ → 0. This asymptotic can be obtained either from the exact rate function (43) and (46), or from a small-λ perturbative expansion applied directly to the MFT equations [7,27].Figure 3 shows Re Q + (k) and Im Q + (k) versus k at X = 1.This figure also shows the same quantities computed by solving Eqs. ( 9) and ( 10) numerically with a back-and-forth iteration algorithm [63].The analytical and numerical curves are in good agreement.
We now deal with the question of the uniqueness of the solution.Let us denote a + (k) = a (k, 0) and a − (k) = ã (k, 0) [using the notation introduced in ( 26)].From Eqs. ( 30) and ( 34), we have Now we turn to Eq. ( 36) and to its solution, given by Eq. (37).In terms of a ± these two equations read respectively.However, there are more solutions to Eq. (48).Given any solution a ± (k) to Eq. ( 48), it is easily seen that, for any k 1 , k 2 ∈ C, where k 1 is in the upper half plane and k 2 in the lower half plane, is also a solution, since only the product a + (k) a − (k) enters in (48).This transformation can be applied any number of times.Which of these solutions is the proper solution to the MFT equations, i.e., the solution which describes histories of u and v that are real, nonnegative and minimize the action s constrained on a given J? One useful test is based on the fact that the solutions for u and v must exhibit a Gaussian decay at x → ±∞.This is because u and v are small there, and therefore the nonlinear terms in the MFT equations ( 9) and ( 10) are negligible, so these equations become the diffusion and anti-diffusion equations respectively.For localized initial conditions, the solutions are Gaussian at x → ±∞ at all times.On the other hand, the factors which appear in the transformation (50) lead, in general, to an exponential rather than a Gaussian decay of u(x, 1) or v(x, 0) at x → ±∞ due to the poles introduced by the Fourier transforms of these functions.The exponential decay is avoided, if k 1 and k 2 solve the equation In this case the pole cancels out, and the inclusion of the extra factor in a ± (k) does not alter the Gaussian decay at infinity.Another possibility is that k i is real (or infinitesimally below or above the real axis).These observations will be useful below.
2. Solution for X < √ 8 at all λ In the analysis up to now, we tacitly assumed that in the term ln 1 + iλke ikX−k 2 branch cuts of the logarithm are never encountered in the integrals given above.This is indeed the case for X = 0 and also at sufficiently small λ for X = 0.However, at larger λ, branch cuts are encountered at X = 0, i.e., for some k ∈ R, The values of k and λ at which these branch cuts cross the real axis can be found by solving the equation 1 + iλke ikX−k 2 = 0 for real k and λ.The solutions are k = ±k R,n , where It follows that for each n for which λ > λ c (X, n), there is an additional branch cut crossing the real axis at k = ±k R,n .The branch cut crossing the real axis would cause Φ ± to be non-analytic as a function of λ if the real axis would be retained as the contour of integration in Eq. ( 39).This would lead to non-analyticity in all the physical quantities (J, S, etc.).To retain analyticity, one must change the contour of integration in Eq. (39), and in Eqs. ( 41), ( 43), (46) stemming from it, such as to avoid crossing any branch cut.Possible contours are shown in Fig. 4, but any contour that does not cross the branch cuts would do.

Solution for X >
√ 8 at all λ At X > √ 8, it is impossible to avoid the branch cuts altogether at sufficiently large λ, see below.In fact, the analytical continuation of the formula of Eq. ( 39) leads to a complex J.We show these statements and obtain the proper solutions for any value λ > 0 below (at λ < 0, the solution for the case X < √ 8 is valid here as well).It is convenient at this point to use a different method to analytically continue the solutions from small λ (as opposed to the one used for X < √ 8, which consisted of changing the contour of integration).Let us first define, for the case X < √ 8, three critical values of λ, which will be important in the following: At λ = λ c1 and λ = λ c2 , the number of purely imaginary solutions to the equation ( 51) changes.Let us assume for now .A convenient choice for integration contours are the red lines while the points at which the branch cuts cross the real axis are denoted by ±kR0 in accordance with Eq. ( 53).In panel (b) the branch points κi are denoted, while in panel (a) three branch points are arbitrarily denoted as ki, i = 0, 1, 3.
as this makes the analysis a little simpler.At λ = λ c0 two branch cuts cross the real axis.The branch points are denoted by κ 1 , κ 2 and the points at which the branch cuts cross the real axis are denoted by ±k R0 .The branch points κ 1 and κ 2 , which satisfy 1 + iλκ i e iκiX−κ2 i = 0, are reflections of each other with respect to the imaginary axis, As λ increases beyond λ c0 , rather than changing the contour of integration, one could instead represent the solution as â( 1) while keeping the contour of integration in the definition of a ± , Eq. ( 49), through Φ, Eq. ( 39), to be along the real axis.â(1) ± (k) is a legitimate solution to Eq. ( 48), as it is reached through a repeated use of the transformation described in Eq. ( 50) 2 .It can be easily shown that that for λ < λ c2 this is equivalent to just changing the contour of integration defining a ± to avoid the branch cuts, while for λ > λ c2 it is impossible to find such a contour.This allows us to define for λ c1 < λ < λ c2 two more solutions: â( 3) Symmetry to reflection with respect to the imaginary axis, which is present since all κ i are purely imaginary at λ c1 < λ < λ c2 , ensures that all three solutions produce real physical fields u(x, 1) and v(x, 0).Only â(1) respects this symmetry for λ < λ c1 , where κ 1 and κ 2 are mirror images of each other.When λ increases beyond λ c2 the branch point κ 2 remains imaginary, while κ 1 and κ 3 are mirror images of each other.This means that beyond λ c2 only â(3) yields real physical fields.Figure 5 shows Re Q + (k) and Im Q + (k) versus k at X = 3.5, λ = 20 > λ c2 , obtained by substituting Eqs. ( 49) and (39) in Eq. ( 59).This figure also shows the same quantities computed by solving Eqs. ( 9) and ( 10) numerically with a back-and-forth iteration algorithm [63].The analytical and numerical curves are in good agreement.
To avoid having to write additional terms in the equation for a ± each time additional branch cuts cross the real axis, it is now convenient, after establishing which solutions are real, to change the contour of integration such that it does not cross any of these other branch cuts.One may choose the curve that satisfies Reλke ikX−k 2 = 0 and passes through the point i X+ √ X 2 −8 4 , see Figure 4 (b).For λ > λ c2 , this curve contains the branch points κ 1 and κ 3 , and in computing the integral, the branch cut of the integrand must be drawn so as not to cross this curve.This change of integration contour from the real axis to this contour affects the change: for all λ > λ c0 , and leaves a ± unchanged for 0 < λ < λ c0 .In this new representation (namely, for this choice of contour of integration), we have the following formulas (note that there are three branches of the solution at λ c1 < λ < λ c2 , whose origin will be explained shortly): As J is increased, λ grows up to the value λ c2 along the branch corresponding to the first line in (61), up to a value J = J 1 .As J grows above J 1 , λ decreases down to a value λ c1 which corresponds to some J = J 2 > J 1 .As J is increased even further, λ increases again: from λ c1 to arbitrary large λ, so that λ → ∞ corresponds to J → 1.The ensuing non-monotonic dependence of λ on J implies that s is a non-convex function of J (this is easy to see from ∂ J s = λ), which is quite an unusual situation in large-deviation theory.J 1 and J 2 are the inflection points, at which ∂ 2 J s = 0.These behaviors are schematically plotted in Fig. 6.We emphasize that Eq. ( 61) holds independently of the assumption (56), so it describes the solution at all X > √ 8 and at all λ > 0. We now calculate Q ± (0), and from it we calculate J = J(λ, X) and s = s(λ, X).By following a similar route to the one above for X < √ 8 , we obtain for J (λ, X) = 1 + Q + (0) /λ the following expression:3 FIG. 6.Schematic plots of (a) J(λ, X), (b) s(λ, X) and (c) s(J, X) for X > √ 8 and λ > 0. Interestingly, J and s are not single-valued functions of λ, and s is not a convex function of J.
and by integrating again Eq. ( 45) with s(λ = 0, X) = 0, we find for the action: a function which we plot in Fig. 2 as a function of J for fixed X = 3.Note that the integration contour here is similar to the one drawn in Fig. 4 (b).As such it is away from the real axis, where the integrands have a pole, and therefore no principal value integral is necessary.Consequently, the Sokhotski-Plemelj formula is also not employed here in contrast, e. g., to Eq. (46).To arrive at Eq. ( 63), one has to evaluate the integral By definition the root κ i satisfies the following equation: so that the integral on the right hand side of Eq. ( 64) can be easily performed: allowing us to write Eq. ( 63).

IV. ASYMPTOTIC LIMITS
Our exact results greatly simplify in three asymptotic limits that we present here.Two of the three asymptotic formulas can be obtained either from the exact results, or directly from the MFT equations by using two different perturbation methods.
A. λ 1: linear theory In the limit of small λ the optimal path can be found by a straightforward perturbative expansion around the zeroth-order solution (3) [7].Here, due to the symmetry relation (14), the leading order calculation in λ 1 is especially simple [27].Using Eq. ( 15), we obtain where we have plugged the zeroth-order solution (3), rescaled by W .The same result (67) follows from the exact rate function, described by Eqs. ( 43) and ( 46), via a straightforward expansion in λ.
Rewriting the shortcut relation ∂s/∂J = λ in the form of (∂s/∂λ)(∂λ/∂J) = λ, and combining it with Eq. (67), we can express λ in terms of J and X: Equations ( 67) and ( 68) yield (still in the rescaled variables) a quadratic in J − J asymptotic which strongly depends on X: The corresponding probability distribution P(J, X, T ) describes typical, small Gaussian fluctuations of the transferred heat around its mean value J.In the original variables In the particular case X = 0 this Gaussian asymptotic coincides with the one obtained in Ref. [27].Note that Eq. ( 69) is applicable only for |λ| 1, that is for |J − J| e −X 2 /2 .This condition becomes restrictive at X 1, where the expected value of the transferred heat is extremely small.Evaluating Eqs. ( 68) and (69) at J = 1/2, one can probe small fluctuations of the energy median x = X m , which is an analog of the tracer position in the SSEP and other single-file systems.At fixed J the applicability condition of the perturbation expansion λ 1 demands that |X m | 1.Using the small-argument asymptotic of erfc(...), we obtain in the leading order λ 2 √ 2X m , and (in the rescaled variables) This corresponds to a Gaussian asymptotic of the probability distribution ρ(X m , T ) of the median.In the original variables The standard deviation of the energy median position, X 2 m 1/2 = (π/2) 1/4 T 1/4 , exhibits the familiar one-dimensional tracer scaling with time, 1. From adiabatic perturbation theory The large-deviation limit λ 1 and X 1 is especially interesting.The reason for that is a unique way of the most probable transfer of a significant amount of energy to large distances, specific for lattice-gas models (such as the KMP) for which the second derivative of the compressibility σ(u) with respect to u [8] is positive.In this case the optimal path has the form of a traveling doublon: a pair of coupled large-amplitude solitary pulses of u and v [8,11].The amplitudes of the u and v pulses slowly (adiabatically) vary with time [8].An adiabatic perturbation theory, which describes this nontrivial process in our present setting, is a slightly modified version of the perturbation theory of Ref. [8], where a step-like initial u-profile were considered, and X was zero.
The adiabatic perturbation theory of Ref. [8] exploits the strong inequality λ 1 which guarantees a large lengthscale separation between the relatively small interior region of the doublon and the large-scale exterior regions in front of the doublon and behind it.In the exterior regions one can neglect the diffusion terms in the MFT equations ( 9) and (10).The diffusion, however, is important for the description of the doublon itself (the interior region).Importantly, when viewed from the exterior regions, the doublon can be effectively described as a single point-like singularity located at a point x = X d (t) which travels in space.The singularity carries a finite energy M (t), a finite jump of the "momentum" V (t) = p 1 − p 2 (where the subscripts 1 and 2 refer to "in front of the singularity" and "behind the singularity", respectively), and a finite value of the Hamiltonian H(t), corresponding to the MFT equations (9) and In its turn, the doublon's momentum V (t) grows exponentially in time with the same rate: which completes the task of determining the optimal path of the doublon in the adiabatic approximation.Note that M (t) and V (t) from Eqs. ( 87) and (88) respectively, together with λ from Eq. ( 86), satisfy the symmetry relation ( 14) integrated over the doublon region (X d (t) − ε, X d (t) + ε).
The rescaled action, according to Ref. [8], is The corresponding large-X asymptotic of the probability distribution is ln P(J, X, T ) −X ln This expression describes a faster-than-exponential tail of P as a function of X.This tail depends only logarithmically on T , and extremely weakly on J/W .

From exact rate function
For large X we are automatically in the regime of X > √ 8.As the contour of integration we can use the curve arg ike ikX−k 2 = 0, which at Rek → ±∞ merges with the line Imk = X 2 .In fact, one can ascertain that for large X the entire curve arg ike ikX−k 2 = 0 has Imk X 2 .Then, in view of Eq. ( 63) for the action, it can be shown that the integrand is vanishingly small unless λ ∼ e X 2 4 .In addition, we shall see a posteriori that we obtain all possible values of J = O(1) by keeping λ much smaller than e X 2 4 .In fact, we obtain J 1 = 0 and J 2 = e −1 .Given all this, the action in this region can be written as: (91) The expression for the current, Eq. ( 62), also benefits from the vanishing of the integral, and we can write approximately: As we shall also see, for J = O(1), that λ X and κ 1,2 ∼ 1 λ .Therefore we can approximate Eq. ( 51) for κ j , j ∈ {1, 2} as follows: This algebraic equation for κ j can be solved in terms of the Lambert function W j−2 : As for κ 3 , in the large X limit, where λ ∼ X, one finds that κ 3 is of order X and has the following form: Then J approximates to Since this equation can give any J from 0 to 1, we thus confirm that J 1 = 0, and only the second and third branches are relevant in Eqs. ( 62) and (63).We also see that J 2 is obtained at the branch point of the Lambert function, −e −1 .Namely J 2 = −e −1 W0(−e −1 ) = e −1 .To find s(J) from Eq. (91) for s(λ), we first write J = i λκj = e iκiX , the last equality being a consequence of the definition of κ i .This gives κ i = −i ln J X .Then, using again the equality J = i λκj , we obtain the desired relation for λ in terms of X and J, which coincides with our second perturbative relation in Eq. ( 86).Using this expression, we can rewrite κ j in terms of J only: Then one can substitute this into the action Eq. ( 91) and arrive at A comparison of Eq. ( 99) with the perturbative result (89) shows that the latter misses an important subleading factor e −2 under the external logarithm.
C. λ 1 and X held constant: J → 1 For large λ we are automatically on branch 3. We assume here that X > √ 8, but the same result can be obtained for X < √ 8 by largely the same means, so we do not expound on the case X < √ 8 here.The integration contour for the integrals in question can be performed on the curve arg ike ikX−k 2 = 0, which at Re k → ±∞ merges with the line Im k = X/2, as shown in Fig. 4 (b).
We shall need the asymptotic expression for κ i in this limit.These are easily computed from Eq. ( 51) to be Computing the first integral in Eq. ( 63), we may first note the symmetry k → −k * , under which the integral simply gets complex conjugated.This means that one may first compute only the integral on that part of the contour in Fig. 4 (b) that runs from ∞ to the imaginary axis, and then take the real value.We decompose this contour into three parts, the first runs from ∞ to κ 3 + δA, where δA ∼ (ln λ) −1/4 .The integrand in this region is vanishingly small and can be dropped.The second part of the contour runs on a small circular arc of radius δA and center κ 3 runnig from κ 3 + δA to κ 3 − δA.On this part of the contour the expression inside the logarithm can be expanded around κ 3 .One gets: For the third part of the contour, we define A = κ 3 + δA and integrate over that part of the contour running from A to the imaginary axis along the contour of the right panel of Fig. 4. In fact, it is easier to integrate from A to −A * , which is ultimately what is needed.Given that, we write: We now expand the last expression taken at the upper limit near κ 3 and use to simplify the combined results of the calculation.Applying the symmetry k → −k * , we find for J, which in turn is given by Eq. ( 62), the following: from which one gets: Thus to this order: To compute the action we use again Eq. ( 45) to write: Re-writing this in terms of J instead of λ using Eq. ( 106), we finally get, in the limit J → 1 at fixed X, s(J → 1, X) The leading term, that is the first term on the right-hand side of Eq. ( 108), coincides with the result obtained in Ref. [27] for the particular case X = 0.As one can see, a finite X affects the J → 1 asymptotic only in a subleading order, described by the second term on the right-hand side of Eq. (108).

V. DISCUSSION
Here we studied the full distribution P(J, X, T ) of the total heat J transferred to the right of a point x = X in the KMP model at long times for an initially localized heat pulse, thereby extending our previous work [27] which was limited to X = 0. We calculated exactly the rescaled rate function s(J, X) that describes this distribution, by combining the MFT and the ISM.In particular we found that for X > √ 8 (in properly rescaled variables) s is not convex as a function of J, a rather unusual situation in large-deviation problems.Furthermore, we studied three different asymptotics of s(J, X, T ) by extracting them from the exact solution and also obtaining two of them by applying two different perturbation methods directly to the MFT equations.The perturbative solutions also yield the optimal histories of the system conditioned on J.This gives an important additional physical insight, at present unaccessible for the ISM.The distribution P(J, X, T ) as a function of X at fixed J has the same long-time dynamical scaling as the distribution of a tracer particle in a single-file transport of interacting particles, as in the SSEP.It would be interesting to try to apply the ISM to other large-deviation problems for the KMP model, in which non-stationary heat transfer statistics problems were previously solved analytically only in some limiting cases [6][7][8].
The MFT of lattice gases can be viewed as a particular case of the optimal fluctuation method (OFM), a versatile tool which can be used to study large deviations in a broad class of macroscopic systems that also includes paradigmatic surface growth models (such as the KPZ equation), reaction-diffusion systems, and many more [40][41][42][43][44][45][46][47][48][49].Mathematically, the OFM formulation involves a set of coupled nonlinear partial differential equations whose solution gives the optimal history of the system, conditioned on the rare event in question.It is generally very hard to solve these equations exactly.Therefore, uncovering and exploiting exact integrability of these equations appears to be a way to make a significant progress in some of these systems.Indeed, the present work is one of several very recent studies [27,37,51,52,65] which use the ISM in order to obtain exact solutions for full statistics of different physical quantities in nonlinear spatio-temporal stochastic systems.We believe that this line of work has a great potential for the future.
Note added.Ref. [65] was posted shortly before the present paper was submitted.The authors of Ref. [65] focused on a crossover between different large-deviation regimes for a model which is mathematically identical to the KMP model.It appears that their work is closely related to ours, and it will be interesting to make a detailed comparison of the results.
It is convenient to solve these equations separately in the two regions x < 0 and x > 0. The general solution reads A − e −ikx/2 − i √ ikB − e −ikx/2 e iky I v (x, y) , x < 0, A + e −ikx/2 − i √ ikB + e −ikx/2 e iky I v (x, y) , x > 0, (B4) B − e ikx/2 , x < 0, B + e ikx/2 , x > 0. (B5) In order to find the four coefficients A ± and B ± (that depend on y and k), we use the following four conditions: (i) T 11 must be continuous at x = 0. (ii) T 21 has a jump discontinuity at x = 0, whose magnitude is found by integrating Eq. (B3) over an infinitesimally small region around x = 0: These conditions yield four linear equations on the coefficients, whose solution then gives the left column of the matrix T (x, y, 0, k) as given in Eq. (28).For the right column of the matrix, the solution is very similar (replacing T 11 → T 12 and T 21 → T 22 ), the difference being that conditions (iii) and (iv) give way to T 12 | x=y = 0 and T 22 | x=y = 1, respectively.Solving the resulting equations for the coefficients then gives the right column of the matrix T (x, y, 0, k) as given in Eq. ( 28).

FIG. 1 .
FIG.1.A Monte-Carlo simulation of the microscopic KMP model.Initially only one agent, at x = 0, had a nonzero energy W = 1.Plotted is the simulated energy profile u as a function of x at time t = 3 × 10 4 (bars), its spatial average over each 50 consecutive lattice sites (solid line), and the theoretical mean-field profile (3) (dashed line).

FIG. 4 .
FIG. 4. Contour plot of arg ike ikX−k 2 for λ = 10 and X = 2 (a) and X = 3 (b), in the complex k plane.Dotted and solid lines correspond to ike ikX−k 2 ∈ R± respectively.The solid lines are branch cuts of ln 1 + iλke ikX−k 2. A convenient choice for integration contours are the red lines while the points at which the branch cuts cross the real axis are denoted by ±kR0 in accordance with Eq. (53).In panel (b) the branch points κi are denoted, while in panel (a) three branch points are arbitrarily denoted as ki, i = 0, 1, 3.