Full statistics of nonstationary heat transfer in the Kipnis–Marchioro–Presutti model

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–4]. Nonstationary current fluctuations, however, proved to be much more difficult to capture [5–11].

A convenient and widely used family of models for studying the full statistics of currents is stochastic lattice gases [12–15]. The Kipnis–Marchioro–Presutti (KMP) model [16] is an important member of this family of models [12, 14]. 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 $\bar{u}(x,t)$ obeys the heat diffusion equation

$$\begin{equation}{\partial }_{t}\bar{u}(x,t)={\partial }_{x}^{2}\bar{u}(x,t).\end{equation} \tag{ 1 }$$

The initial coarse-grained temperature profile is a delta-function,

$$\begin{equation}\bar{u}(x,t=0)=W\delta (x),\end{equation} \tag{ 2 }$$

so the mean-field solution for the coarse-grained temperature is

$$\begin{equation}\bar{u}(x,t)=\frac{W}{\sqrt{4\pi t}}\,\mathrm{exp}\left(-\frac{{x}^{2}}{4t}\right).\end{equation} \tag{ 3 }$$

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 figure 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 ${\mathcal{P}}_{0}({J}_{0},T)$ of the total amount of heat J₀, observed at time t = T ≫ 1 on the right half-line:

$$\begin{equation}{J}_{0}={\int }_{0}^{\infty }u(x,t=T)\mathrm{d}x.\end{equation} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

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–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–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 [37].

The evaluation of the full distribution ${\mathcal{P}}_{0}({J}_{0},T)$ for the quenched (that is deterministic) initial condition (2) was achieved in [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 (ISM) [38, 39]. In this work we extend the approach of [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 $\mathcal{P}(J,X,T)$ of the total amount of heat,

$$\begin{equation}J={\int }_{X}^{\infty }u(x,t=T)\mathrm{d}x\end{equation} \tag{ 5 }$$

observed at time t = T ≫ 1 to the right of x = X. Because of the energy conservation in the KMP model, the distribution $\mathcal{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):

$$\begin{equation}\bar{J}=\bar{J}(X,T)={\int }_{X}^{\infty }\bar{u}(x,T)\mathrm{d}x=\frac{W}{2}\text{erfc}\left(\frac{X}{\sqrt{4T}}\right),\end{equation} \tag{ 6 }$$

where erfc z = 1 − erf z, and erf z is the error function. For X = 0 one obtains $\bar{J}=W/2$. The main objective of this work is to determine the full statistics of fluctuations of J around the mean value $\bar{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 $\mathcal{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 [40–42]. The OFM has been also applied to turbulence [43], stochastic reactions [44, 45], non-equilibrium surface growth and related models [46–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/\sqrt{T}\ll 1$ [27]. For a general one-component diffusive lattice gas, the fluctuational hydrodynamics has the form of the stochastic partial differential equation [12, 14]

$$\begin{equation}{\partial }_{t}u={\partial }_{x}\left[D(u){\partial }_{x}u+\sqrt{\sigma (u)}\eta \right],\end{equation} \tag{ 7 }$$

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, $\left\langle \eta (x,t)\right\rangle =0$ and $\quad \left\langle \eta (x,t)\eta ({x}^{\prime },{t}^{\prime })\right\rangle =\delta (x-{x}^{\prime })\delta (t-{t}^{\prime })$. For the KMP model, where u(x, t) is the coarse-grained temperature, one obtains D = 1 and σ(u) = 2u², and equation (7) becomes

$$\begin{equation}{\partial }_{t}u={\partial }_{x}\left({\partial }_{x}u+\sqrt{2}u\eta \right).\end{equation} \tag{ 8 }$$

Without the noise term equation (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 equation (7). In our problem, the small parameter of the saddle-point evaluation is again $1/\sqrt{T}$. The saddle-point 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 $\mathcal{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. Bettelheim et al [27] reported a first major advance in this area of statistical mechanics by finding an exact solution to the MFT problem for the distribution $\mathcal{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 ISM. A second major advance in this area has been recently reported in [37]. The authors employed the MFT for determining the full statistics of J defined in equation (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 [6]. The authors of [37] reproduced this asymptotic by combining a generalization of the canonical Cole–Hopf transformation with the ISM.

The recent works [27, 37] (and [50–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 ISM for obtaining exact solutions for different large-deviation statistics. The present work is a natural next step in exploring this potential.

In section 2 we present the formulation of the MFT problem for $\mathcal{P}(J,X,T)$ [6, 8] and briefly discuss the scaling behavior of $\mathcal{P}(J,X,T)$ and its relation with the universal scaling behavior of the tracer statistics in single-file processes. In section 3 we solve the MFT problem by the ISM. In section 4 we discuss some asymptotic limits of $\mathcal{P}(J,X,T)$. Our results are briefly summarized in section 5. Some of the technical details of the calculations are given in the appendices.

We rescale the variables t, x and u by T, $\sqrt{T}$ and $W/\sqrt{T}$, respectively. As a result, X and J become rescaled by $\sqrt{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

$$\begin{equation}{\partial }_{t}u={\partial }_{x}({\partial }_{x}u+2{u}^{2}v),\end{equation} \tag{ 9 }$$

$$\begin{equation}{\partial }_{t}v={\partial }_{x}(-{\partial }_{x}v+2u{v}^{2}).\end{equation} \tag{ 10 }$$

To make this paper more self-contained, a derivation of equations (9) and (10) is given in appendix A. The initial condition, after the rescaling, is

$$\begin{equation}u(x,t=0)=\delta (x).\end{equation} \tag{ 11 }$$

The rescaled transferred heat at t = T is constrained by the equation

$$\begin{equation}{\int }_{X}^{\infty }u(x,t=1)\mathrm{d}x=J.\end{equation} \tag{ 12 }$$

The particular case X = 0 was solved in [27]. Being interested in X ≠ 0 we can, without loss of generality, set X > 0. The constrained optimization problem yields, aside from equations (9) and (10), a second boundary condition in time [6],

$$\begin{equation}v(x,t=1)=-\lambda \delta (x-X),\end{equation} \tag{ 13 }$$

where λ is a Lagrange multiplier, to be ultimately determined from equation (12). This MFT formulation possesses a nontrivial symmetry

$$\begin{equation}v(x,t)=-\lambda u\left(X-x,1-t\right),\end{equation} \tag{ 14 }$$

which generalizes the symmetry found for the particular case X = 0 in [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–8]:

$$\begin{equation}s={\int }_{0}^{1}\mathrm{d}t{\int }_{-\infty }^{\infty }\mathrm{d}x\,{u}^{2}{v}^{2}.\end{equation} \tag{ 15 }$$

This action determines the probability density $\mathcal{P}(J,X,T)$ up to a pre-exponential factor. In the original, dependent variables we have

$$\begin{equation}\mathrm{ln}\,\mathcal{P}(J,X,T)\simeq -\sqrt{T}s\left(\frac{J}{W},\frac{X}{\sqrt{T}}\right),\end{equation} \tag{ 16 }$$

which describes the scaling behavior of $\mathcal{P}(J,X,T)$. As we will see in the following, at fixed J/W and small $X/\sqrt{T}$, the rate function s is quadratic in $X/\sqrt{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 [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/\sqrt{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 equations (9) and (10) coincide with the derivative nonlinear Schrö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 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. Bettelheim et al [27] overcame this difficulty in the particular case X = 0 by (i) making use of a shortcut (see, e.g. [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 [27] to an arbitrary X.

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 equations (9) and (10) is the existence of a corresponding Lax pair. This means than equations (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).

3.1. Scattering amplitudes and their dynamics

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

$$\begin{equation}\begin{cases}{\partial }_{x}\boldsymbol{\psi }(x,t,k)=U(x,t,k)\boldsymbol{\psi }(x,t,k),\quad \hfill \\ {\partial }_{t}\boldsymbol{\psi }(x,t,k)=V(x,t,k)\boldsymbol{\psi }(x,t,k),\quad \hfill \end{cases}\end{equation} \tag{ 17 }$$

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

$$\begin{equation}U(x,t,k)=\left(\begin{matrix}\hfill -\text{i}k/2\hfill & \hfill -\text{i}v\sqrt{\text{i}k}\hfill \\ \hfill -\text{i}u\sqrt{ik}\hfill & \hfill \text{i}k/2\hfill \end{matrix}\right),\end{equation} \tag{ 18 }$$

$$\begin{equation}V(x,t,k)=\left(\begin{matrix}\hfill {k}^{2}/2-\text{i}kuv\hfill & \hfill -\text{i}{(\sqrt{\text{i}k})}^{3}v+\text{i}\sqrt{\text{i}k}{\partial }_{x}v-\text{i}\sqrt{\text{i}k}2{v}^{2}u\hfill \\ \hfill -\text{i}{(\sqrt{\text{i}k})}^{3}u+\text{i}\sqrt{\text{i}k}{\partial }_{x}u-\text{i}\sqrt{\text{i}k}2{u}^{2}v\hfill & \hfill -{k}^{2}/2+\text{i}kuv\hfill \end{matrix}\right)\end{equation} \tag{ 19 }$$

constitute the Lax pair, and k is a spectral parameter. It is now straightforward to check that the compatibility condition ∂_t ∂_x ψ = ∂_x ∂_t ψ , which corresponds to

$$\begin{equation}{\partial }_{t}U-{\partial }_{x}V+\left[U,V\right]=0,\end{equation} \tag{ 20 }$$

is indeed equivalent to equations (9) and (10).

It is useful to define the matrix $\mathcal{T}(x,y,t,k)$ as the x-propagator of the system (17), namely, the solution to the equation

$$\begin{equation}{\partial }_{x}\mathcal{T}(x,y,t,k)=U(x,t,k)\mathcal{T}(x,y,t,k)\end{equation} \tag{ 21 }$$

with the initial condition $\mathcal{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,

$$\begin{equation}U(x\to \pm \infty ,t,k)=\left(\begin{matrix}\hfill -\text{i}k/2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \text{i}k/2\hfill \end{matrix}\right).\end{equation} \tag{ 22 }$$

As a result, the entries of $\mathcal{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:

$$\begin{equation}G(t,k)=\underset{\underset{y\to \,-\infty }{x\to \infty }}{\mathrm{lim}}\left(\begin{matrix}\hfill {\text{e}}^{\text{i}kx/2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\text{e}}^{-\text{i}kx/2}\hfill \end{matrix}\right)\mathcal{T}(x,y,t,k)\left(\begin{matrix}\hfill {\text{e}}^{-\text{i}ky/2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\text{e}}^{\text{i}ky/2}\hfill \end{matrix}\right).\end{equation} \tag{ 23 }$$

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 $\mathcal{T}(x,y,t,k)$ satisfies:

$$\begin{equation}{\partial }_{t}\mathcal{T}(x,y,t,k)=V(x,t,k)\mathcal{T}(x,y,t,k)-\mathcal{T}(x,y,t,k)V(y,t,k).\end{equation} \tag{ 24 }$$

As one can verify directly, equation (24) is compatible with (21) (i.e. ${\partial }_{t}{\partial }_{x}\mathcal{T}={\partial }_{x}{\partial }_{t}\mathcal{T}$) due to equation (20). The matrix V(x, t, k) too becomes very simple in the limit x → ±∞,

$$\begin{equation}V(x\to \pm \infty ,t,k)=\frac{{k}^{2}}{2}\left(\begin{matrix}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill \end{matrix}\right).\end{equation} \tag{ 25 }$$

Plugging equation (25) into (24), one finds the time evolution of $\mathcal{T}(x\to \infty ,y\to -\infty ,t,k)$ which in turn, using (23), yields that of G(t, k):

$$\begin{equation}G(t,k)=\left(\begin{matrix}\hfill a(t,k)\hfill & \hfill \tilde{b}(t,k)\hfill \\ \hfill b(t,k)\hfill & \hfill \tilde{a}(t,k)\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill a(0,k)\hfill & \hfill \tilde{b}(0,k){\text{e}}^{{k}^{2}t}\hfill \\ \hfill b(0,k){\text{e}}^{-{k}^{2}t}\hfill & \hfill \tilde{a}(0,k)\hfill \end{matrix}\right)\end{equation} \tag{ 26 }$$

In equation (26) we have introduced a notation for the matrix elements of G(t, k).

3.2. Solving the scattering problem at t = 0 and t = 1

Let us find the matrix $\mathcal{T}(x,y,0,k)$ at t = 0. This is achieved by solving equation (21) at t = 0, which reads (using u(x, 0) = δ(x))

$$\begin{equation}{\partial }_{x}\mathcal{T}\left(x,y,0,k\right)=\left(\begin{matrix}\hfill -\frac{\text{i}k}{2}\hfill & \hfill -\text{i}v\left(x,0\right)\sqrt{\text{i}k}\hfill \\ \hfill -\text{i}\delta \left(x\right)\sqrt{\text{i}k}\hfill & \hfill \frac{\text{i}k}{2}\hfill \end{matrix}\right)\mathcal{T}\left(x,y,0,k\right).\end{equation} \tag{ 27 }$$

The general solution to equation (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 $\mathcal{T}(x,x,t,k)=I$, we obtain (see appendix B for details)

$$\begin{equation}\mathcal{T}(x,y,0,k)=\begin{cases}\left(\begin{matrix}\hfill {\text{e}}^{-\text{i}k(x-y)/2}\hfill & \hfill -\text{i}\sqrt{\text{i}k}\,{\text{e}}^{-\text{i}k(x-y)/2}{I}_{v}(x,y)\hfill \\ \hfill 0\hfill & \hfill {\text{e}}^{\text{i}k(x-y)/2}\hfill \end{matrix}\right),\quad \hfill & xy > 0,\hfill \\ \left(\begin{matrix}\hfill {\text{e}}^{\text{i}k(y-x)/2}\left[1\pm ik{I}_{v}(x,0)\right]\hfill & \hfill -\text{i}\sqrt{\text{i}k}\,{\text{e}}^{\text{i}k(y-x)/2}\left[{I}_{v}(x,y)\pm \text{i}k{I}_{v}(0,y){I}_{v}(x,0)\right]\hfill \\ \hfill \pm \text{i}\sqrt{\text{i}k}{\text{e}}^{\text{ik}(x+y)/2}\hfill & \hfill {\text{e}}^{\text{i}k(x-y)/2}\pm \text{i}k\,{\text{e}}^{\text{i}k(x+y)/2}{I}_{v}(0,y)\hfill \end{matrix}\right),\quad \hfill & xy< 0,\hfill \end{cases}\end{equation} \tag{ 28 }$$

where

$$\begin{equation}{I}_{v}(x,y)={\int }_{y}^{x}v(z,0){\text{e}}^{\text{i}k(z-y)}\mathrm{d}z,\end{equation} \tag{ 29 }$$

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 equation (28) into (23), with the result:

$$\begin{equation}G(0,k)=\left(\begin{matrix}\hfill 1-\text{i}k{Q}_{+}(k)\hfill & \hfill -\text{i}\sqrt{\text{i}k}\left[Q(k)-\text{i}k{Q}_{-}(k){Q}_{+}(k)\right]\hfill \\ \hfill -\text{i}\sqrt{\text{i}k}\hfill & \hfill 1-\text{i}k{Q}_{-}(k)\hfill \end{matrix}\right)\end{equation} \tag{ 30 }$$

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:

$$\begin{equation}{Q}_{-}(k)={\int }_{-\infty }^{0}\!v(z,0){\text{e}}^{\text{i}kz}\,\mathrm{d}z,\quad {Q}_{+}(k)={\int }_{0}^{\infty }\!v(z,0){\text{e}}^{\text{i}kz}\,\mathrm{d}z,\quad Q(k)={Q}_{+}(k)+{Q}_{-}(k).\end{equation} \tag{ 31 }$$

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 equation (21) at t = 1. The solution is very similar to that at t = 0, and one gets

$$\begin{align}\hfill & \mathcal{T}(x,y,1,k)\hfill \\ \hfill & \qquad =\begin{cases}\left(\begin{matrix}\hfill {\text{e}}^{-\text{i}k(x-y)/2}\hfill & \hfill 0\hfill \\ \hfill -\text{i}\sqrt{\text{i}k}\,{\text{e}}^{\text{i}k(x-y)/2}{I}_{u}(x,y)\hfill & \hfill {\text{e}}^{\text{i}k(x-y)/2}\hfill \end{matrix}\right),\quad \hfill & \left(x-X\right)\left(y-X\right) > 0,\hfill \\ \left(\begin{matrix}\hfill {\text{e}}^{-\text{i}k(x-y)/2}\pm \lambda \text{i}k\,{\text{e}}^{\text{i}k(x+y-2X)/2}{I}_{u}(X,y)\hfill & \hfill \pm \text{i}\lambda \sqrt{\text{i}k}\,{\text{e}}^{-\text{i}k(x+y-2X)/2}\hfill \\ \hfill -\text{i}\sqrt{\text{i}k}\,{\text{e}}^{\text{i}k(x-y)/2}\left[{I}_{u}(x,y)\pm \lambda \text{i}k{I}_{u}(X,y){I}_{u}(x,X)\right]\hfill & \hfill {\text{e}}^{\text{i}k(x-y)/2}\left[1\pm \lambda \text{i}k{I}_{u}(x,X)\right]\hfill \end{matrix}\right),\quad \hfill & \left(x-X\right)\left(y-X\right)< 0,\hfill \end{cases}\hfill \end{align} \tag{ 32 }$$

where

$$\begin{equation}{I}_{u}(x,y)={\int }_{y}^{x}u(z,1){\text{e}}^{-\text{i}k(z-y)}\mathrm{d}z,\end{equation} \tag{ 33 }$$

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

$$\begin{equation}G(1,k)=\left(\begin{matrix}\hfill 1+\lambda \text{i}k{R}_{+}(k)\hfill & \hfill +\text{i}\lambda \sqrt{\text{i}k}\,{\text{e}}^{\text{i}kX}\hfill \\ \hfill -\text{i}\sqrt{\text{i}k}\,{\text{e}}^{-\text{i}kX}\left[R(k)+\lambda \text{i}k{R}_{-}(k){R}_{+}(k)\right]\hfill & \hfill 1+\lambda \text{i}k{R}_{-}(k)\hfill \end{matrix}\right)\end{equation} \tag{ 34 }$$

where

$$\begin{equation}\begin{gathered}{R}_{+}(k)={\int }_{-\infty }^{X}u(z,1){\text{e}}^{-\text{i}kz}\,\mathrm{d}z,\quad {R}_{-}(k)={\int }_{X}^{\infty }u(z,1){\text{e}}^{-\text{i}kz}\,\mathrm{d}z,\\ R(k)={R}_{+}(k)+{R}_{-}(k).\end{gathered}\end{equation} \tag{ 35 }$$

3.3. 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< \sqrt{8}$. (iii) For $X > \sqrt{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–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 > \sqrt{8}$.

3.3.1. Full solution for X = 0, and for X > 0 with sufficiently small λ

Comparing the upper-right elements of G(0, k) from equation (30) and of G(1, k) from equation (34), using equation (26), we obtain

$$\begin{equation}\text{i}kQ(k)-\text{i}k{Q}_{-}(k)\times \text{i}k{Q}_{+}(k)=-\lambda \text{i}k\,{\text{e}}^{\text{i}kX-{k}^{2}}.\end{equation} \tag{ 36 }$$

Equation (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 equation (36). We begin by completing the squares in equation (36), by writing

$$\begin{equation}\text{i}k{Q}_{\pm }(k)=1-{\text{e}}^{{{\Phi}}_{\pm }(k)-{{\Phi}}_{\pm }(0)},\end{equation} \tag{ 37 }$$

where the constant −Φ_±(0) in the exponent is found by requiring that Q_±(k) is regular at the origin. This equation turns equation (36) into

$$\begin{equation}{\text{e}}^{{{\Phi}}_{+}(k)+{{\Phi}}_{-}(k)-{{\Phi}}_{+}(0)-{{\Phi}}_{-}(0)}=1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}},\end{equation} \tag{ 38 }$$

which has the solution

$$\begin{equation}{{\Phi}}_{\pm }(k)=\pm {\int }_{-\infty }^{\infty }\frac{\mathrm{ln}\left(1+\text{i}\lambda {k}^{\prime }\,{\text{e}}^{\text{i}{k}^{\prime }X-{k}^{\prime 2}}\right)}{{k}^{\prime }-k\mp \text{i}{0}^{+}}\frac{\mathrm{d}{k}^{\prime }}{2\pi \text{i}}.\end{equation} \tag{ 39 }$$

Equation (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-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 half-planes, f_±(k), respectively, given by

$$\begin{equation}{f}_{\pm }(k)=\int \frac{f({k}^{\prime })}{{k}^{\prime }-k\mp \text{i}{0}^{+}}\frac{\mathrm{d}{k}^{\prime }}{2\pi \text{i}}.\end{equation} \tag{ 40 }$$

The decomposition (40) is applied to the logarithm of equation (38), leading to equation (39).

Taking the derivative of equation (37) with respect to k, setting k to 0 and dividing by i we obtain:

$$\begin{equation}{Q}_{\pm }(0)=\pm {\int }_{-\infty }^{\infty }\frac{\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)}{{(k\mp \text{i}{0}^{+})}^{2}}\frac{\mathrm{d}k}{2\pi }\end{equation} \tag{ 41 }$$

which, using the Sokhotski–Plemelj formula, can be also written in terms of a principal value integral,

$$\begin{equation}{Q}_{\pm }(0)=\pm {\fint}_{-\infty }^{\infty }\frac{\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)}{{k}^{2}}\frac{\mathrm{d}k}{2\pi }\mp \frac{\lambda }{2}.\end{equation} \tag{ 42 }$$

Using equations (12), (14) and (31) alongside with the conservation law ${\int }_{-\infty }^{\infty }u(x,t)\mathrm{d}x=1$, we determine J = J(λ, X):

$$\begin{align}\hfill J\left(\lambda ,X\right)& ={\int }_{X}^{\infty }u\left(x,t=1\right)\mathrm{d}x=1-{\int }_{-\infty }^{X}u\left(x,t=1\right)\mathrm{d}x=1+\frac{1}{\lambda }{\int }_{0}^{\infty }v\left(x,0\right)\mathrm{d}x\hfill \\ \hfill & =1+\frac{{Q}_{+}\left(0\right)}{\lambda }=1+{\int }_{-\infty }^{\infty }\frac{\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)}{{(k-\text{i}{0}^{+})}^{2}}\frac{\mathrm{d}k}{2\pi \lambda }\hfill \end{align} \tag{ 43 }$$

which can also be written as

$$\begin{equation}J\left(\lambda ,X\right)={\fint}_{-\infty }^{\infty }\frac{\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)}{{k}^{2}}\frac{\mathrm{d}k}{2\pi \lambda }-\frac{1}{2}.\end{equation} \tag{ 44 }$$

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. [56]. This enables us to calculate s(J, X) bypassing equation (15) (which would require us to know the entire optimal path u(x, t)). Using $J\left(\lambda ,X\right)=1+{Q}_{+}\left(0\right)/\lambda$ from equation (43), we have

$$\begin{equation}\frac{\partial s}{\partial \lambda }=\frac{\partial s}{\partial J}\frac{\partial J}{\partial \lambda }=\lambda \frac{\partial J}{\partial \lambda }=\frac{\partial {Q}_{+}\left(0\right)}{\partial \lambda }-\frac{{Q}_{+}\left(0\right)}{\lambda }.\end{equation} \tag{ 45 }$$

Using equation (42), we integrate equation (45) with respect to λ to get

$$\begin{equation}s(\lambda ,X)={Q}_{+}(0)+{\int }_{-\infty }^{\infty }\frac{{\text{Li}}_{2}\left(-\text{i}k\lambda \,{\text{e}}^{\text{i}kX-{k}^{2}}\right)}{2\pi {k}^{2}}\,\mathrm{d}k\,+\frac{\lambda }{2}.\end{equation} \tag{ 46 }$$

where ${\text{Li}}_{2}\left(z\right)={\sum }_{k=1}^{\infty }{z}^{k}/{k}^{2}$ is the dilogarithm function, Q₊(0) is given by equation (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.

Figure 2 shows s(J, X) alongside with the asymptotic $J\to \bar{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 equations (9) and (10) numerically with a back-and-forth iteration algorithm [63]. The analytical and numerical curves are in good agreement.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now deal with the question of the uniqueness of the solution. Let us denote ${a}_{+}\left(k\right)=a\left(k,0\right)$ and ${a}_{-}\left(k\right)=\tilde{a}\left(k,0\right)$ (using the notation introduced in (26)). From equations (30) and (34), we have

$$\begin{equation}{a}_{\pm }\left(k\right)=1-\text{i}k{Q}_{\pm }\left(k\right)=1+\lambda \text{i}k{R}_{\pm }\left(k\right).\end{equation} \tag{ 47 }$$

Now we turn to equation (36) and to its solution, given by equation (37). In terms of a_± these two equations read

$$\begin{equation}1-{a}_{+}\left(k\right){a}_{-}\left(k\right)=-\lambda \text{i}k\,{\text{e}}^{\text{i}kX-{k}^{2}},\end{equation} \tag{ 48 }$$

$$\begin{equation}{a}_{\pm }\left(k\right)={\text{e}}^{{{\Phi}}_{\pm }(k)-{{\Phi}}_{\pm }(0)},\end{equation} \tag{ 49 }$$

respectively. However, there are more solutions to equation (48). Given any solution ${a}_{\pm }\left(k\right)$ to equation (48), it is easily seen that, for any ${k}_{1},{k}_{2}\in \mathbb{C}$, where k₁ is in the upper half plane and k₂ in the lower half plane,

$$\begin{equation}{\hat{a}}_{\pm }\left(k\right)={\left(\frac{{k}_{2}}{{k}_{1}}\frac{k-{k}_{1}}{k-{k}_{2}}\right)}^{\pm 1}{a}_{\pm }\left(k\right)\end{equation} \tag{ 50 }$$

is also a solution, since only the product ${a}_{+}\left(k\right){a}_{-}\left(k\right)$ 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₁ and k₂ solve the equation

$$\begin{equation}1+\text{i}\lambda {k}_{i}{\text{e}}^{\text{i}{k}_{i}X-{k}_{i}^{2}}=0.\end{equation} \tag{ 51 }$$

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.

3.3.2. Solution for $X< \sqrt{8}$ at all λ

In the analysis up to now, we tacitly assumed that in the term $\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)$ 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\in \mathbb{R}$,

$$\begin{equation}1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\in {\mathbb{R}}_{-}.\end{equation} \tag{ 52 }$$

The values of k and λ at which these branch cuts cross the real axis can be found by solving the equation $1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}=0$ for real k and λ. The solutions are k = ±k_R,n , where

$$\begin{equation}\qquad {k}_{R,n}=\frac{\left(n+\frac{1}{2}\right)\pi }{X},\qquad n=0,1,2,\dots ,\end{equation} \tag{ 53 }$$

$$\begin{equation}{\lambda }_{c}\left(X,n\right)=\frac{{\left(-1\right)}^{n}X}{\left(n+\frac{1}{2}\right)\pi }\mathrm{exp}\left[\frac{{\left(n+\frac{1}{2}\right)}^{2}{\pi }^{2}}{{X}^{2}}\right].\end{equation} \tag{ 54 }$$

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 equation (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 equation (39), and in equations (41), (43) and (46) stemming from it, such as to avoid crossing any branch cut. Possible contours are shown in figure 4, but any contour that does not cross the branch cuts would do.

Zoom In Zoom Out Reset image size

3.3.3. Solution for $X > \sqrt{8}$ at all λ

At $X > \sqrt{8}$, it is impossible to avoid the branch cuts altogether at sufficiently large λ, see below. In fact, the analytical continuation of the formula of equation (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< \sqrt{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< \sqrt{8}$, which consisted of changing the contour of integration). Let us first define, for the case $X< \sqrt{8}$, three critical values of λ, which will be important in the following:

$$\begin{equation}{\lambda }_{\text{c}0}={\lambda }_{\text{c}}\left(X,0\right)< {\lambda }_{\text{c}1}=\frac{4\,{\text{e}}^{\frac{1}{2}+\frac{{X}^{2}-X\sqrt{{X}^{2}-8}}{8}}}{X-\sqrt{{X}^{2}-8}}< {\lambda }_{\text{c}2}=\frac{4\,{\text{e}}^{\frac{1}{2}+\frac{{X}^{2}+X\sqrt{{X}^{2}-8}}{8}}}{X+\sqrt{{X}^{2}-8}}.\end{equation} \tag{ 55 }$$

At λ = λ_c1 and λ = λ_c2, the number of purely imaginary solutions to the equation (51) changes.

Let us assume for now

$$\begin{equation}{\lambda }_{\text{c}2}< \underset{n=2,4,6,\dots }{\mathrm{min}}{\lambda }_{\text{c}}\left(X,n\right),\end{equation} \tag{ 56 }$$

as this makes the analysis a little simpler. At λ = λ_c0 two branch cuts cross the real axis. The branch points are denoted by κ₁, κ₂ and the points at which the branch cuts cross the real axis are denoted by ±k_R0. The branch points κ₁ and κ₂, which satisfy $1+\text{i}\lambda {\kappa }_{i}\,{\text{e}}^{\text{i}{\kappa }_{i}X-{\kappa }_{i}^{2}}=0$, are reflections of each other with respect to the imaginary axis, ${\kappa }_{1}=-{\kappa }_{2}^{\ast }$. As λ increases beyond λ_c0, rather than changing the contour of integration, one could instead represent the solution as

$$\begin{equation}{\hat{a}}_{\pm }^{(1)}\left(k\right)={\left(\frac{{k}_{R}^{2}}{{\kappa }_{1}{\kappa }_{2}}\frac{(k-{\kappa }_{1})(k-{\kappa }_{2})}{{k}_{R}^{2}-{k}^{2}}\right)}^{\pm 1}{a}_{\pm }\left(k\right),\end{equation} \tag{ 57 }$$

while keeping the contour of integration in the definition of a_±, equation (49), through Φ, equation (39), to be along the real axis. ${\hat{a}}_{\pm }^{(1)}\left(k\right)$ is a legitimate solution to equation (48), as it is reached through a repeated use of the transformation described in equation (50). ⁴ It can be easily shown 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:

$$\begin{equation}{\hat{a}}_{\pm }^{(2)}\left(k\right)={\left(\frac{{k}_{R}^{2}}{{\kappa }_{3}{\kappa }_{2}}\frac{(k-{\kappa }_{3})(k-{\kappa }_{2})}{{k}_{R}^{2}-{k}^{2}}\right)}^{\pm 1}{a}_{\pm }\left(k\right),\end{equation} \tag{ 58 }$$

$$\begin{equation}{\hat{a}}_{\pm }^{(3)}\left(k\right)={\left(\frac{{k}_{R}^{2}}{{\kappa }_{3}{\kappa }_{1}}\frac{(k-{\kappa }_{3})(k-{\kappa }_{1})}{{k}_{R}^{2}-{k}^{2}}\right)}^{\pm 1}{a}_{\pm }\left(k\right).\end{equation} \tag{ 59 }$$

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 ${\hat{a}}^{(1)}$ respects this symmetry for λ < λ_c1, where κ₁ and κ₂ are mirror images of each other. When λ increases beyond λ_c2 the branch point κ₂ remains imaginary, while κ₁ and κ₃ are mirror images of each other. This means that beyond λ_c2 only ${\hat{a}}^{(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 equations (49) and (39) in equation (59). This figure also shows the same quantities computed by solving equations (9) and (10) numerically with a back-and-forth iteration algorithm [63]. The analytical and numerical curves are in good agreement.

Zoom In Zoom Out Reset image size

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 $\mathrm{R}\mathrm{e}\,\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}=0$ and passes through the point $\text{i}\frac{X+\sqrt{{X}^{2}-8}}{4}$, see figure 4(b). For λ > λ_c2, this curve contains the branch points κ₁ and κ₃, 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:

$$\begin{equation}{a}_{\pm }\to {\left(\frac{{k}_{R}^{2}}{{\kappa }_{1}{\kappa }_{2}}\frac{(k-{\kappa }_{1})(k-{\kappa }_{2})}{{k}_{R}^{2}-{k}^{2}}\right)}^{\pm 1}{a}_{\pm },\end{equation} \tag{ 60 }$$

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):

$$\begin{equation}{a}_{\pm }\left(k\right)=\begin{cases}{\text{e}}^{{{\Phi}}_{\pm }(k)-{{\Phi}}_{\pm }(0)},\quad \hfill & \lambda < {\lambda }_{\text{c}2},\hfill \\ {\left(\frac{{\kappa }_{1}}{{\kappa }_{3}}\frac{k-{\kappa }_{3}}{k-{\kappa }_{1}}\right)}^{\pm 1}\,{\text{e}}^{{{\Phi}}_{\pm }(k)-{{\Phi}}_{\pm }(0)},\quad \hfill & {\lambda }_{\text{c}1}< \lambda < {\lambda }_{\text{c}2},\hfill \\ {\left(\frac{{\kappa }_{2}}{{\kappa }_{3}}\frac{k-{\kappa }_{3}}{k-{\kappa }_{2}}\right)}^{\pm 1}\,{\text{e}}^{{{\Phi}}_{\pm }(k)-{{\Phi}}_{\pm }(0)},\quad \hfill & \lambda > {\lambda }_{\text{c}1}.\hfill \end{cases}\end{equation} \tag{ 61 }$$

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₁. As J grows above J₁, λ decreases down to a value λ_c1 which corresponds to some J = J₂ > J₁. 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₁ and J₂ are the inflection points, at which ${\partial }_{J}^{2}s=0$. These behaviors are schematically plotted in figure 6. We emphasize that equation (61) holds independently of the assumption (56), so it describes the solution at all $X > \sqrt{8}$ and at all λ > 0.

Zoom In Zoom Out Reset image size

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< \sqrt{8}$, we obtain for $J\left(\lambda ,X\right)=1+{Q}_{+}\left(0\right)/\lambda$ the following expression: ⁵

$$\begin{equation}J\left(\lambda ,X\right)=\int \frac{\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)}{{k}^{2}}\frac{\mathrm{d}k}{2\pi \lambda }+\begin{cases}0,\quad \hfill & J< {J}_{1},\hfill \\ \frac{\text{i}}{\lambda }\left(\frac{1}{{\kappa }_{1}}-\frac{1}{{\kappa }_{3}}\right),\quad \hfill & {J}_{1}< J< {J}_{2},\hfill \\ \frac{\text{i}}{\lambda }\left(\frac{1}{{\kappa }_{2}}-\frac{1}{{\kappa }_{3}}\right),\quad \hfill & J > {J}_{2},\hfill \end{cases}\end{equation} \tag{ 62 }$$

and by integrating again equation (45) with s(λ = 0, X) = 0, we find for the action:

$$\begin{align}\hfill s\left(\lambda ,X\right)& =\int \!\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)\frac{\mathrm{d}k}{2\pi {k}^{2}}+\!\int \!{\text{Li}}_{2}\left(-\text{i}k\lambda \,{\text{e}}^{\text{i}kX-{k}^{2}}\right)\frac{\mathrm{d}k}{2\pi {k}^{2}}\hfill \\ \hfill & \quad +\begin{cases}0,\quad \hfill & J< {J}_{1},\hfill \\ 2\text{i}\left({\kappa }_{3}-{\kappa }_{1}\right)+X\,\mathrm{ln}\left(\frac{{\kappa }_{3}}{{\kappa }_{1}}\right),\quad \hfill & {J}_{1}< J< {J}_{2},\hfill \\ 2\text{i}\left({\kappa }_{3}-{\kappa }_{2}\right)+X\,\mathrm{ln}\left(\frac{{\kappa }_{3}}{{\kappa }_{2}}\right),\quad \hfill & J > {J}_{2},\hfill \end{cases}\hfill \end{align} \tag{ 63 }$$

a function which we plot in figure 2 as a function of J for fixed X = 3. Note that the integration contour here is similar to the one drawn in figure 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 equation (46). To arrive at equation (63), one has to evaluate the integral

$$\begin{equation}{\int }^{\lambda }\frac{\mathrm{d}\,\mathrm{ln}{\lambda }^{\prime }}{{\kappa }_{i}({\lambda }^{\prime })}=\frac{\mathrm{ln}(\lambda )}{{\kappa }_{i}(\lambda )}+{\int }^{{\kappa }_{i}(\lambda )}\frac{\mathrm{ln}\,\lambda ({\kappa }_{i}^{\prime })}{{\kappa }_{i}^{\prime 2}}\,\mathrm{d}{\kappa }_{i}^{\prime }.\end{equation} \tag{ 64 }$$

By definition the root ${\kappa }_{i}^{\prime }$ satisfies the following equation:

$$\begin{equation}\mathrm{ln}(\lambda ({\kappa }_{i}^{\prime }))=\frac{\text{i}\pi }{2}-\mathrm{ln}\,{\kappa }_{i}^{\prime }-\text{i}{\kappa }_{i}^{\prime }X+{\kappa }_{i}^{\prime 2},\end{equation} \tag{ 65 }$$

so that the integral on the right-hand side of equation (64) can be easily performed:

$$\begin{equation}{\int }^{\lambda }\frac{\mathrm{d}\,\mathrm{ln}{\lambda }^{\prime }}{{\kappa }_{i}({\lambda }^{\prime })}=-\text{i}X\,\mathrm{ln}{\kappa }_{i}(\lambda )+\frac{1}{{\kappa }_{i}(\lambda )}+2{\kappa }_{i}(\lambda )\end{equation} \tag{ 66 }$$

allowing us to write equation (63).

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.

4.1.  λ ≪ 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 equation (15), we obtain

$$\begin{equation}s(\lambda ,X)\simeq {\lambda }^{2}{\int }_{0}^{1}\mathrm{d}t{\int }_{-\infty }^{\infty }\mathrm{d}x\,{\bar{u}}^{2}(x,t){\bar{u}}^{2}(X-x,1-t)=\frac{{\lambda }^{2}\,{\text{e}}^{-\frac{{X}^{2}}{2}}}{8\sqrt{2\pi }},\end{equation} \tag{ 67 }$$

where we have plugged the zeroth-order solution (3), rescaled by W. The same result (67) follows from the exact rate function, described by equations (43) and (46), via a straightforward expansion in λ.

Rewriting the shortcut relation ∂s/∂J = λ in the form of (∂s/∂λ)(∂λ/∂J) = λ, and combining it with equation (67), we can express λ in terms of J and X:

$$\begin{equation}\lambda \simeq 4\sqrt{2\pi }\,{\text{e}}^{\frac{{X}^{2}}{2}}(J-\bar{J}).\end{equation} \tag{ 68 }$$

Equations (67) and (68) yield (still in the rescaled variables) a quadratic in $J-\bar{J}$ asymptotic which strongly depends on X:

$$\begin{equation}s(J\to \bar{J},X)\simeq 2\sqrt{2\pi }\,{\text{e}}^{\frac{{X}^{2}}{2}}{(J-\bar{J})}^{2}.\end{equation} \tag{ 69 }$$

The corresponding probability distribution $\mathcal{P}(J,X,T)$ describes typical, small Gaussian fluctuations of the transferred heat around its mean value $\bar{J}$. In the original variables

$$\begin{equation}\mathrm{ln}{\mathcal{P}}_{\text{Gauss}}\simeq -2\sqrt{2\pi T}\,{\text{e}}^{\frac{{X}^{2}}{2T}}{\left[\frac{J}{W}-\frac{1}{2}\,\text{erfc}\left(\frac{X}{\sqrt{4T}}\right)\right]}^{2}.\end{equation} \tag{ 70 }$$

In the particular case X = 0 this Gaussian asymptotic coincides with the one obtained in [27]. Note that equation (69) is applicable only for |λ| ≪ 1, that is for $\vert J-\bar{J}\vert \ll {\text{e}}^{-{X}^{2}/2}$. This condition becomes restrictive at X ≫ 1, where the expected value of the transferred heat is extremely small.

Evaluating equations (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 $\lambda \simeq 2\sqrt{2}{X}_{\text{m}}$, and (in the rescaled variables)

$$\begin{equation}s\simeq \frac{{X}_{\text{m}}^{2}}{\sqrt{2\pi }}.\end{equation} \tag{ 71 }$$

This corresponds to a Gaussian asymptotic of the probability distribution ρ(X_m, T) of the median. In the original variables

$$\begin{equation}\mathrm{ln}\,\rho ({X}_{\text{m}},T)\simeq -\frac{{X}_{\text{m}}^{2}}{\sqrt{2\pi T}}.\end{equation} \tag{ 72 }$$

The standard deviation of the energy median position, ${\langle {X}_{\text{m}}^{2}\rangle }^{1/2}={(\pi /2)}^{1/4}{T}^{1/4}$, exhibits the familiar one-dimensional tracer scaling with time, T^1/4.

4.2.  λ ≫ 1, X ≫ 1

4.2.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 [8], where a step-like initial u-profile were considered, and X was zero.

The adiabatic perturbation theory of [8] exploits the strong inequality λ ≫ 1 which guarantees a large length-scale 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₁ − p₂ (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 $\mathcal{H}(t)$, corresponding to the MFT equations (9) and (10) rewritten in the Hamiltonian form (that is, in the canonical variables u and p). The quantities M(t), V(t) and $\mathcal{H}(t)$ slowly depend on time. They are formally defined as

$$\begin{equation}\left\{M(t),V(t),\mathcal{H}(t)\right\}={\int }_{{X}_{\mathrm{d}}(t)-\varepsilon }^{{X}_{\mathrm{d}}(t)+\varepsilon }\left\{u,-v,h\right\}\mathrm{d}x,\end{equation} \tag{ 73 }$$

where h = v ∂_x u + u² v² is the density of the Hamiltonian of equations (9) and (10) rewritten in terms of the canonical variables u and p, and ɛ > 0 can be fixed at an arbitrary small value as X → ∞ and λ → ∞. It was shown by direct calculations in [8] that

$$\begin{equation}M(t)\simeq \frac{1}{{\tilde{v}}_{1}(t)},\quad V(t)\simeq \frac{1}{{u}_{2}(t)},\quad \text{and}\quad \mathcal{H}(t)\simeq c(t),\end{equation} \tag{ 74 }$$

where ${\tilde{v}}_{1}(t)$ is the value of v in front of the doublon, and u₂(t) is the value of u behind the doublon. As we will see shortly, the doublon velocity c(t) is actually constant. Using equations (57)–(59) and (62)–(64) of [8] and the relation ${\dot{X}}_{\mathrm{d}}(t)=c(t)$, we arrive at the equations

$$\begin{equation}\dot{M}=-c{u}_{2},\end{equation} \tag{ 75 }$$

$$\begin{equation}\dot{V}=c{\tilde{v}}_{1},\end{equation} \tag{ 76 }$$

$$\begin{equation}\dot{\mathcal{H}}=0.\end{equation} \tag{ 77 }$$

Then, using equation (74), we obtain

$$\begin{equation}\dot{M}=-\frac{c}{V},\end{equation} \tag{ 78 }$$

$$\begin{equation}\dot{V}=\frac{c}{M},\end{equation} \tag{ 79 }$$

$$\begin{equation}c=\text{const}.\end{equation} \tag{ 80 }$$

The doublon must arrive at t = 1 at the point x = X, which gives c = X. Equations (78) and (79) have the first integral

$$\begin{equation}M(t)V(t)={k}_{1}(X),\end{equation} \tag{ 81 }$$

where the function k₁(X) is a priori unknown. Using this first integral in equation (78), we obtain

$$\begin{equation}M(t)=\mathrm{exp}\left[{k}_{2}(X)-\frac{Xt}{{k}_{1}(X)}\right].\end{equation} \tag{ 82 }$$

The a priori unknown function k₂(X) must be set to zero in view of the condition M(t = 0) = 1, and we obtain

$$\begin{equation}M(t)=\mathrm{exp}\left[-\frac{Xt}{{k}_{1}(X)}\right].\end{equation} \tag{ 83 }$$

Then equation (81) yields

$$\begin{equation}V(t)={k}_{1}(X)\mathrm{exp}\left[\frac{Xt}{{k}_{1}(X)}\right].\end{equation} \tag{ 84 }$$

Using the two remaining boundary conditions [64],

$$\begin{equation}M(t=1)=j\quad \text{and}\quad V(t=1)=\lambda ,\end{equation} \tag{ 85 }$$

we determine k₁(X) and λ:

$$\begin{equation}{k}_{1}(X)=\frac{X}{\mathrm{ln}\frac{1}{j}}\quad \text{and}\quad \lambda =\frac{X}{j\,\mathrm{ln}\,\frac{1}{j}}.\end{equation} \tag{ 86 }$$

Using the first of these relations in equation (83), we see that the energy M(t) of the doublon decays exponentially in time with rate ln(1/j):

$$\begin{equation}M(t)=\mathrm{exp}\left(-t\mathrm{ln}\frac{1}{j}\right).\end{equation} \tag{ 87 }$$

In its turn, the doublon's momentum V(t) grows exponentially in time with the same rate:

$$\begin{equation}V(t)=\frac{X}{\mathrm{ln}\frac{1}{j}}\mathrm{exp}\left(t\,\mathrm{ln}\,\frac{1}{j}\right),\end{equation} \tag{ 88 }$$

which completes the task of determining the optimal path of the doublon in the adiabatic approximation. Note that M(t) and V(t) from equations (87) and (88) respectively, together with λ from equation (86), satisfy the symmetry relation (14) integrated over the doublon region (X_d(t) − ɛ, X_d(t) + ɛ).

The rescaled action, according to [8], is

$$\begin{equation}s(j,X)\simeq c\,\mathrm{ln}(cMV)=X\,\mathrm{ln}\frac{{X}^{2}}{\mathrm{ln}\,\frac{1}{j}},\quad X\gg 1.\end{equation} \tag{ 89 }$$

The corresponding large-X asymptotic of the probability distribution is

$$\begin{equation}\mathrm{ln}\,\mathcal{P}(J,X,T)\simeq -X\,\mathrm{ln}\left(\frac{{X}^{2}}{T\,\mathrm{ln}\,\frac{W}{J}}\right),\quad X\gg \sqrt{T}.\end{equation} \tag{ 90 }$$

This expression describes a faster-than-exponential tail of $\mathcal{P}$ as a function of X. This tail depends only logarithmically on T, and extremely weakly on J/W.

4.2.2. From exact rate function

For large X we are automatically in the regime of $X > \sqrt{8}$. As the contour of integration we can use the curve $\mathrm{arg}\,\text{i}k\,{\text{e}}^{\text{i}kX-{k}^{2}}=0$, which at Re k → ±∞ merges with the line $\mathrm{I}\mathrm{m}\,k=\frac{X}{2}$. In fact, one can ascertain that for large X the entire curve $\mathrm{arg}\,\text{i}k\,{\text{e}}^{\text{i}kX-{k}^{2}}=0$ has $\mathrm{I}\mathrm{m}\,k\simeq \frac{X}{2}$. Then, in view of equation (63) for the action, it can be shown that the integrand is vanishingly small unless $\lambda \sim {\text{e}}^{\frac{{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 ${\text{e}}^{\frac{{X}^{2}}{4}}$. In fact, we obtain J₁ = 0 and J₂ = e⁻¹. Given all this, the action in this region can be written as:

$$\begin{equation}s\left(\lambda ,X\right)\simeq \begin{cases}2\text{i}\left({\kappa }_{3}(\lambda )-{\kappa }_{1}(\lambda )\right)+X\,\mathrm{ln}\left(\frac{{\kappa }_{3}(\lambda )}{{\kappa }_{1}(\lambda )}\right),\quad \hfill & {J}_{1}< J< {J}_{2},\hfill \\ 2\text{i}\left({\kappa }_{3}(\lambda )-{\kappa }_{2}(\lambda )\right)+X\,\mathrm{ln}\left(\frac{{\kappa }_{3}(\lambda )}{{\kappa }_{2}(\lambda )}\right),\quad \hfill & J > {J}_{2}.\hfill \end{cases}\end{equation} \tag{ 91 }$$

The expression for the current, equation (62), also benefits from the vanishing of the integral, and we can write approximately:

$$\begin{equation}J\left(\lambda ,X\right)\simeq \begin{cases}\frac{\text{i}}{\lambda }\left(\frac{1}{{\kappa }_{1}}-\frac{1}{{\kappa }_{3}}\right),\quad \hfill & {J}_{1}< J< {J}_{2},\hfill \\ \frac{\text{i}}{\lambda }\left(\frac{1}{{\kappa }_{2}}-\frac{1}{{\kappa }_{3}}\right),\quad \hfill & J > {J}_{2}.\hfill \end{cases}\end{equation} \tag{ 92 }$$

As we shall also see, for J = O(1), that λ ≃ X and ${\kappa }_{1,2}\sim \frac{1}{\lambda }$. Therefore we can approximate equation (51) for κ_j , j ∈ {1, 2} as follows:

$$\begin{equation}-\frac{1}{\lambda }=\text{i}{\kappa }_{j}\,{\text{e}}^{\text{i}{\kappa }_{j}X}.\end{equation} \tag{ 93 }$$

This algebraic equation for κ_j can be solved in terms of the Lambert function W_j−2:

$$\begin{equation}{\kappa }_{j}=-\frac{\text{i}}{X}{W}_{j-2}\left(-\frac{X}{\lambda }\right),\quad j\in \left\{1,2\right\}.\end{equation} \tag{ 94 }$$

As for κ₃, in the large X limit, where λ ∼ X, one finds that κ₃ is of order X and has the following form:

$$\begin{equation}{\kappa }_{3}=\text{i}X-\text{i}\frac{\mathrm{ln}(\lambda X)}{X}.\end{equation} \tag{ 95 }$$

Then J approximates to

$$\begin{equation}J\sim \frac{\text{i}}{\lambda {\kappa }_{j}}=-\frac{X}{\lambda }{W}_{j-2}^{-1}\left(-\frac{X}{\lambda }\right),\quad j\in \left\{1,2\right\}\end{equation} \tag{ 96 }$$

Since this equation can give any J from 0 to 1, we thus confirm that J₁ = 0, and only the second and third branches are relevant in equations (62) and (63). We also see that J₂ is obtained at the branch point of the Lambert function, −e⁻¹. Namely ${J}_{2}=\frac{-{\text{e}}^{-1}}{{W}_{0}(-{\text{e}}^{-1})}={\text{e}}^{-1}$.

To find s(J) from equation (91) for s(λ), we first write $J=\frac{\text{i}}{\lambda {\kappa }_{j}}={\text{e}}^{\text{i}{\kappa }_{i}X}$, the last equality being a consequence of the definition of κ_i . This gives ${\kappa }_{i}=-\text{i}\frac{\mathrm{ln}\,J}{X}$. Then, using again the equality $J=\frac{\text{i}}{\lambda {\kappa }_{j}}$, we obtain the desired relation for λ in terms of X and J,

$$\begin{equation}\quad \lambda =\frac{X}{J\,\mathrm{ln}\,\frac{1}{J}},\end{equation} \tag{ 97 }$$

which coincides with our second perturbative relation in equation (86). Using this expression, we can rewrite κ_j in terms of J only:

$$\begin{equation}{\kappa }_{j}=\begin{cases}\frac{\text{i}\,\mathrm{ln}(1/J)}{X},\quad \hfill & j\in \left\{1,2\right\},\hfill \\ \text{i}X-\frac{\text{i}}{X}\,\mathrm{ln}\left(\frac{{X}^{2}}{J\,\mathrm{ln}(1/J)}\right),\quad \hfill & j=3,\hfill \end{cases}\end{equation} \tag{ 98 }$$

Then one can substitute this into the action equation (91) and arrive at

$$\begin{equation}s(J,X)=X\,\mathrm{ln}\left(\frac{{X}^{2}}{{\text{e}}^{2}\,\mathrm{ln}\,\frac{1}{J}}\right)+O(1).\end{equation} \tag{ 99 }$$

A comparison of equation (99) with the perturbative result (89) shows that the latter misses an important subleading factor e⁻² under the external logarithm.

4.3.  λ ≫ 1 and X held constant: J → 1

For large λ we are automatically on branch 3. We assume here that $X > \sqrt{8}$, but the same result can be obtained for $X< \sqrt{8}$ by largely the same means, so we do not expound on the case $X< \sqrt{8}$ here. The integration contour for the integrals in question can be performed on the curve $\mathrm{arg}\,\text{i}k\,{\text{e}}^{\text{i}kX-{k}^{2}}=0$, which at Re k → ±∞ merges with the line Im k = X/2, as shown in figure 4(b).

We shall need the asymptotic expression for κ_i in this limit. These are easily computed from equation (51) to be

$$\begin{equation}{\kappa }_{2}^{-1}=-\text{i}\lambda +\text{i}X-\frac{\text{i}}{\lambda }+O\left(\frac{1}{{\lambda }^{2}}\right)\end{equation} \tag{ 100 }$$

$$\begin{equation}{\kappa }_{3}=-{\kappa }_{1}^{\ast }=\frac{\text{i}X}{2}+\sqrt{\mathrm{ln}\,\lambda }-\frac{\text{i}\pi +\frac{{X}^{2}}{2}-\mathrm{ln}\,\mathrm{ln}\,\lambda }{4\sqrt{\mathrm{ln}\,\lambda }}+O\left(\frac{1}{\mathrm{ln}\,\lambda }\right).\end{equation} \tag{ 101 }$$

Computing the first integral in equation (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 figure 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 κ₃ + δ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 κ₃ running from κ₃ + δA to κ₃ − δA. On this part of the contour the expression inside the logarithm can be expanded around κ₃. One gets:

$$\begin{equation*}\int \frac{\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)}{{k}^{2}}\frac{\mathrm{d}k}{2\pi }=\frac{\delta {A}^{2}}{2\pi {\kappa }_{3}}+O\left(\frac{1}{{(\mathrm{ln}\,\lambda )}^{3/2}}\right).\end{equation*}$$

For the third part of the contour, we define A = κ₃ + δA and integrate over that part of the contour running from A to the imaginary axis along the contour of the right panel of figure 4. In fact, it is easier to integrate from A to −A*, which is ultimately what is needed. Given that, we write:

$$\begin{align}\hfill {\int }_{-{A}^{\ast }}^{A}\frac{\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)}{{k}^{2}}\frac{\mathrm{d}k}{2\pi }& ={\int }_{-{A}^{\ast }}^{A}\frac{\mathrm{ln}\left(\text{i}\lambda k\right)+\text{i}kX-{k}^{2}}{{k}^{2}}\frac{\mathrm{d}k}{2\pi }+O\left(\frac{1}{\lambda }\right)\hfill \\ \hfill & =-{\left.\frac{\mathrm{ln}\left(\text{i}\lambda k\right)+1}{2\pi k}+\frac{\text{i}X\mathrm{ln}\,k-k}{2\pi }\right\vert }_{-{A}^{\ast }}^{A}+O\left(\frac{1}{\lambda }\right).\hfill \end{align} \tag{ 102 }$$

We now expand the last expression taken at the upper limit near κ₃ and use

$$\begin{equation}\mathrm{ln}(\text{i}\lambda {\kappa }_{3})=-\text{i}\pi -\text{i}{\kappa }_{3}X+{\kappa }_{3}^{2}\end{equation} \tag{ 103 }$$

to simplify the combined results of the calculation. Applying the symmetry k → −k*, we find for J, which in turn is given by equation (62), the following:

$$\begin{align}\hfill J(\lambda ,X)& =\int \frac{\mathrm{ln}\left(1+\text{i}\lambda k\,{\text{e}}^{\text{i}kX-{k}^{2}}\right)}{\lambda {k}^{2}}\frac{\mathrm{d}k}{2\pi }-\frac{\text{i}}{\lambda {\kappa }_{3}}+\frac{\text{i}}{\lambda {\kappa }_{2}}\hfill \\ \hfill & \simeq \frac{{Q}_{+}(0)}{\lambda }\simeq 1-\frac{2{(\mathrm{ln}\,\lambda )}^{1/2}}{\pi \lambda }-\frac{{(\mathrm{ln}\,\lambda )}^{-1/2}}{\pi \lambda }\left(1+\frac{{X}^{2}}{4}+\frac{1}{2}\,\mathrm{ln}\,\mathrm{ln}\,\lambda \right),\hfill \end{align} \tag{ 104 }$$

from which one gets:

$$\begin{equation}-\frac{{\pi }^{2}}{2}{(1-J)}^{2}=-\frac{2\,\mathrm{ln}\,\lambda }{{\lambda }^{2}}\left[1+\frac{1+\frac{{X}^{2}}{4}+\frac{1}{2}\,\mathrm{ln}\,\mathrm{ln}\,\lambda }{\mathrm{ln}\lambda }+O\left(\frac{1}{{\mathrm{ln}}^{2}\,\lambda }\right)\right]\end{equation} \tag{ 105 }$$

Thus to this order:

$$\begin{equation}\mathrm{ln}\,\lambda =-\frac{1}{2}{W}_{-1}\left(-\frac{{\pi }^{2}}{2}{(1-J)}^{2}\right)+1+\frac{{X}^{2}}{4}+\frac{1}{2}\,\mathrm{ln}\left[-\frac{1}{2}{W}_{-1}\left(-\frac{{\pi }^{2}}{2}{(1-J)}^{2}\right)\right].\end{equation} \tag{ 106 }$$

To compute the action we use again equation (45) to write:

$$\begin{equation}s(\lambda ,X)={Q}_{+}(0)-\int {Q}_{+}(0)\frac{\mathrm{d}\lambda }{\lambda }\simeq \frac{4{(\mathrm{ln}\,\lambda )}^{3/2}}{3\pi }+\frac{{X}^{2}}{\pi }{(\mathrm{ln}\,\lambda )}^{1/2}.\end{equation} \tag{ 107 }$$

Re-writing this in terms of J instead of λ using equation (106), we finally get, in the limit J → 1 at fixed X,

$$\begin{align}\hfill s(J\to 1,X)& \simeq \frac{4}{3\pi }{\left[-\frac{1}{2}{W}_{-1}\left(-\frac{{\pi }^{2}}{2}{(1-J)}^{2}\right)\right]}^{3/2}+\frac{2}{\pi }{\left[-\frac{1}{2}{W}_{-1}\left(-\frac{{\pi }^{2}}{2}{(1-J)}^{2}\right)\right]}^{1/2}\hfill \\ \hfill & \quad \times \left[1+\frac{3{X}^{2}}{4}+\frac{1}{2}\,\mathrm{ln}\left(-\frac{1}{2}{W}_{-1}\left(-\frac{{\pi }^{2}}{2}{(1-J)}^{2}\right)\right)\right].\hfill \end{align} \tag{ 108 }$$

The leading term, that is the first term on the right-hand side of equation (108), coincides with the result obtained in [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 equation (108).

Here we studied the full distribution $\mathcal{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 > \sqrt{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 $\mathcal{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–8].

The MFT of lattice gases can be viewed as a particular case of the 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–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. Krajenbrink et al [65] was posted shortly before the present paper was submitted. The authors of [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.

The research of E B and B M is supported by the Israel Science Foundation (Grant No. 1466/15 and 1499/20, respectively).

Let us first explain why the effective amplitude of the noise becomes small in the long-time limit T ≫ 1. We rescale time and space, $\tilde{t}=t/T$, $\tilde{x}=x/\sqrt{T}$, and define $\tilde{u}\left(\tilde{x},\tilde{t}\right)\equiv u\left(\tilde{x}\sqrt{T},\tilde{t}T\right)$, and the dimensionless noise term $\tilde{\eta }\left(\tilde{x},\tilde{t}\right)={T}^{3/4}\eta \left(\tilde{x}\sqrt{T},\tilde{t}T\right)$. Under these rescalings, equation (8) of the main text becomes

$$\begin{equation}{\partial \,}_{\tilde{t}}\tilde{u}={\partial }_{\tilde{x}}\left({\partial }_{\tilde{x}}\tilde{u}+{T}^{-1/4}\sqrt{2}\tilde{u}\tilde{\eta }\right)\end{equation} \tag{ A1 }$$

from which it is clearly observed that the noise is effectively weak in the limit T ≫ 1. Note, however, that the long-time limit must be taken with constant $X/\sqrt{T}$.

We now derive of the MFT equations, equations (9) and (10) of the main text. We first define a 'potential'

$$\begin{equation}\psi \left(x,t\right)={\int }_{-\infty }^{x}u\left(y,t\right)\mathrm{d}y.\end{equation} \tag{ A2 }$$

Integrating equation (8) 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):

$$\begin{equation}{\partial }_{t}\psi ={\partial }_{x}^{2}\psi +\sqrt{2}{\partial }_{x}\psi \eta .\end{equation} \tag{ A3 }$$

An equivalent description of the dynamics (A3) is provided by the path-integral representation. The probability density functional P[η] of the white Gaussian noise term η(x, t) is

$$\begin{equation}P\left[\eta \right]\sim \mathrm{exp}\left(-{\int }_{0}^{T}\mathrm{d}t{\int }_{-\infty }^{\infty }\mathrm{d}x\,\frac{{\eta }^{2}}{2}\right).\end{equation} \tag{ A4 }$$

The advantage of using the ψ representation is that the noise term in equation (A3) appears without a spatial derivative, and one can therefore express η through ψ via equation (A3). Using this in equation (A4), one finds that the probability density for a given realization of ψ(x, t) is:

$$\begin{equation}-\mathrm{ln}\,P\left[\psi \right]\simeq S\left[\psi \right]\equiv \frac{1}{4}{\int }_{0}^{T}\mathrm{d}t{\int }_{-\infty }^{\infty }\mathrm{d}x{\left(\frac{{\partial }_{t}\psi -{\partial }_{x}^{2}\psi }{{\partial }_{x}\psi }\right)}^{2},\end{equation} \tag{ A5 }$$

where $S\left[\psi \right]$ is the action functional. As explained above, in the limit $\sqrt{T}\gg 1$, the noise is effectively weak. This enables us to calculate $\mathcal{P}(J,T)$ via a saddle-point evaluation of the path integral. Within this framework, $-\mathrm{ln}\,\mathcal{P}\left(J,T\right)$ is given by the minimum of the action functional S with respect to all realizations ψ(x, t), constrained on the initial condition u(x, 0) = Wδ(x) and on a given value of the transferred heat J = W − ψ(X, T) (where we used the conservation law ${\int }_{-\infty }^{\infty }u(x,t=T)\mathrm{d}x=W$). The latter constraint on J is taken into account by adding the term $-{\Lambda}\psi \left(0,T\right)$ to the action, where Λ is a Lagrange multiplier, i.e. we minimize the constrained functional

$$\begin{equation}{S}_{{\Lambda}}\left[\psi \right]=S\left[\psi \right]-{\Lambda}\psi \left(0,T\right)=\frac{1}{4}{\int }_{0}^{T}\mathrm{d}t{\int }_{-\infty }^{\infty }\mathrm{d}x{\left(\frac{{\partial }_{t}\psi -{\partial }_{x}^{2}\psi }{{\partial }_{x}\psi }\right)}^{2}-{\Lambda}\psi \left(X,T\right).\end{equation} \tag{ A6 }$$

To minimize S_Λ, we must require its linear variation with respect to ψ(x, t) to vanish. A given variation ψ(x, t) → ψ(x, t) + δψ(x, t) causes a variation of ${S}_{{\Lambda}}\left[\psi \right]$ which, to first order in δψ, is given by

$$\begin{align}\hfill \delta {S}_{{\Lambda}}& ={S}_{{\Lambda}}\left[\psi +\delta \psi \right]-{S}_{{\Lambda}}\left[\psi \right]\hfill \\ \hfill & =\frac{1}{2}{\int }_{0}^{T}\mathrm{d}t{\int }_{-\infty }^{\infty }\mathrm{d}x\left(\frac{{\partial }_{x}^{2}\psi -{\partial }_{t}\psi }{{\partial }_{x}\psi }\right)\left[\frac{{\partial }_{x}^{2}\delta \psi -{\partial }_{t}\delta \psi }{{\partial }_{x}\psi }-\frac{{\partial }_{x}^{2}\psi -{\partial }_{t}\psi }{{\left({\partial }_{x}\psi \right)}^{2}}{\partial }_{x}\delta \psi \right]\,-{\Lambda}\delta \psi \left(X,T\right)\hfill \\ \hfill & ={\int }_{0}^{T}\mathrm{d}t{\int }_{-\infty }^{\infty }\mathrm{d}x\,{\partial }_{x}p\left({\partial }_{x}^{2}\delta \psi -{\partial }_{t}\delta \psi -2{\partial }_{x}p{\partial }_{x}\psi {\partial }_{x}\delta \psi \right)-{\Lambda}\delta \psi \left(X,T\right),\hfill \end{align} \tag{ A7 }$$

where we have defined the momentum density gradient

$$\begin{equation}{\partial }_{x}p=\frac{{\partial }_{x}^{2}\psi -{\partial }_{t}\psi }{2{\left({\partial }_{x}\psi \right)}^{2}}.\end{equation} \tag{ A8 }$$

Integrating by parts in equation (A7), we obtain

$$\begin{align}\hfill \delta {S}_{{\Lambda}}& ={\int }_{0}^{T}\mathrm{d}t{\int }_{-\infty }^{\infty }\mathrm{d}x\left\{{\partial }_{x}^{3}p+{\partial }_{xt}p+2{\partial }_{x}\left[{\left({\partial }_{x}p\right)}^{2}{\partial }_{x}\psi \right]\right\}\delta \psi \hfill \\ \hfill & \quad +{\int }_{-\infty }^{\infty }\mathrm{d}x\left[{\partial }_{x}p\left(x,0\right)\delta \psi \left(x,0\right)-{\partial }_{x}p\left(x,T\right)\delta \psi \left(x,T\right)-{\Lambda}\delta \left(x-X\right)\delta \psi \left(x,T\right)\right],\hfill \end{align} \tag{ A9 }$$

where the first two terms in the single integral are the boundary terms originating from the integration by parts in time. Since ψ(x, 0) is given (i.e. the initial condition is quenched), the variation $\delta \psi \left(x,0\right)$ vanishes.

By requiring the double integral in (A9) to vanish for arbitrary δψ, one obtains the second MFT equation, equation (10) in the main text (recalling that v = −∂_x p and u = ∂_x ψ). The first MFT equation, equation (9) in the main text, follows from equation (A8) 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 and X should be rescaled by 1/W and $\sqrt{T}$, respectively. The MFT equations invariant under these rescalings. In addition, we must require the single integral in equation (A9) to vanish for arbitrary $\delta \psi \left(x,T\right)$. This yields the boundary condition

$$\begin{equation}{\partial }_{x}p\left(x,T\right)=-{\Lambda}\delta (x-X).\end{equation} \tag{ A10 }$$

After the rescaling, this becomes equation (13) in the main text, where λ = WΛ is a rescaled Lagrange multiplier, and X is rescaled as well. Finally, using equation (A8) in (A5) one finds that the action can be rewritten as $S={\int }_{0}^{T}\mathrm{d}t{\int }_{-\infty }^{\infty }\mathrm{d}x\,{u}^{2}{v}^{2}$ which, after the rescaling, leads to $S=\sqrt{T}s$ where the rescaled action s is given by equation (15) in the main text. The large parameter $S\sim \sqrt{T}\gg 1$ justifies a posteriori the saddle-point approximation that we used.

In this appendix, we solve equation (27) in the main text, obtaining equation (28). Let us introduce the notation

$$\begin{equation}\mathcal{T}\left(x,y,0,k\right)=\left(\begin{matrix}\hfill {\mathcal{T}}_{11}\hfill & \hfill {\mathcal{T}}_{12}\hfill \\ \hfill {\mathcal{T}}_{21}\hfill & \hfill {\mathcal{T}}_{21}\hfill \end{matrix}\right)\end{equation} \tag{ B1 }$$

for the entries of the matrix, with the x, y, k dependence suppressed for brevity. We solve equation (27) separately for the two columns of the matrix. For the left column, equation (27) reads

$$\begin{equation}{\partial }_{x}{\mathcal{T}}_{11}=-\frac{\text{i}k}{2}{\mathcal{T}}_{11}-\text{i}v\left(x,0\right)\sqrt{\text{i}k}{\mathcal{T}}_{21},\end{equation} \tag{ B2 }$$

$$\begin{equation}{\partial }_{x}{\mathcal{T}}_{21}=-\text{i}\delta \left(x\right)\sqrt{\text{i}k}{\mathcal{T}}_{11}+\frac{\text{i}k}{2}{\mathcal{T}}_{21}.\end{equation} \tag{ B3 }$$

It is convenient to solve these equations separately in the two regions x < 0 and x > 0. The general solution reads

$$\begin{equation}{\mathcal{T}}_{11}=\begin{cases}{A}_{-}\,{\text{e}}^{-\text{i}kx/2}-\text{i}\sqrt{\text{i}k}{B}_{-}\,{\text{e}}^{-\text{i}kx/2}\,{\text{e}}^{\text{i}ky}{I}_{v}\left(x,y\right),\quad \hfill & x< 0,\hfill \\ {A}_{+}\,{\text{e}}^{-\text{i}kx/2}-\text{i}\sqrt{\text{i}k}{B}_{+}\,{\text{e}}^{-\text{i}kx/2}\,{\text{e}}^{\text{i}ky}{I}_{v}\left(x,y\right),\quad \hfill & x > 0,\hfill \end{cases}\end{equation} \tag{ B4 }$$

$$\begin{equation}{\mathcal{T}}_{21}=\begin{cases}{B}_{-}\,{\text{e}}^{\text{i}kx/2},\quad \hfill & x< 0,\hfill \\ {B}_{+}\,{\text{e}}^{\text{i}kx/2},\quad \hfill & x > 0.\hfill \end{cases}\end{equation} \tag{ B5 }$$

In order to find the four coefficients A_± and B_± (that depend on y and k), we use the following four conditions: (i) ${\mathcal{T}}_{11}$ must be continuous at x = 0. (ii) ${\mathcal{T}}_{21}$ has a jump discontinuity at x = 0, whose magnitude is found by integrating equation (B3) over an infinitesimally small region around x = 0:

$$\begin{equation}{\left.{\mathcal{T}}_{21}\right\vert }_{x={0}^{+}}-{\left.{\mathcal{T}}_{21}\right\vert }_{x={0}^{-}}=-\text{i}\sqrt{\text{i}k}{\left.{\mathcal{T}}_{11}\right\vert }_{x=0}.\end{equation} \tag{ B6 }$$

(iii) ${\left.{\mathcal{T}}_{11}\right\vert }_{x=y}=1$. (iv) ${\left.{\mathcal{T}}_{21}\right\vert }_{x=y}=0$. The latter two conditions come from the boundary condition $\mathcal{T}(x,x,t,k)=I$. These conditions yield four linear equations on the coefficients, whose solution then gives the left column of the matrix $\mathcal{T}\left(x,y,0,k\right)$ as given in equation (28). For the right column of the matrix, the solution is very similar (replacing ${\mathcal{T}}_{11}\to {\mathcal{T}}_{12}$ and ${\mathcal{T}}_{21}\to {\mathcal{T}}_{22}$), the difference being that conditions (iii) and (iv) give way to ${\left.{\mathcal{T}}_{12}\right\vert }_{x=y}=0$ and ${\left.{\mathcal{T}}_{22}\right\vert }_{x=y}=1$, respectively. Solving the resulting equations for the coefficients then gives the right column of the matrix $\mathcal{T}\left(x,y,0,k\right)$ as given in equation (28).