The number statistics and optimal history of non-equilibrium steady states of mortal diffusing particles

Fluctuations of non-equilibrium steady states of driven diffusive lattice gases have attracted a lot of attention in the last two decades [1–7]. Although many of these studies assumed purely diffusive, particle conserving dynamics, lattice gas models with dissipation have also been investigated [8–11]. The absence of detailed balance makes dissipative systems more difficult to handle. There is, however, a steady interest in dissipative models, mostly because of their relevance to experiment in such diverse areas as fluid turbulence [12], granular gases [13, 14] and many other non-equilibrium settings in physics, chemistry, biology and engineering.

The simplest way to characterize a non-equilibrium system is to study its steady state. To maintain a dissipative systems in a steady state one must constantly drive it by injecting energy. The energy of the driven system fluctuates around the mean, and these fluctuations bear the stamp of the non-equilibrium nature of the system. In this work we consider three different driven models of diffusive particles where particles can die: either individually or by annihilating in pairs. In each of these models the system is driven by injecting particles into a half-infinite straight line from a single point-like steady source. The particles diffuse and die so that, at long times, a non-equilibrium steady state sets in. If the particles do not interact, their total number N is Poisson-distributed with mean $\bar{N}$ predicted from a deterministic reaction-diffusion equation for this driven system. Our main interest in this case will be to find the optimal (that is, most probable) density history of the system conditional on observing a given N. We also consider two types of interacting particle models: (i) diffusive lattice gases of mortal particles with constant diffusivity but non-trivial fluctuations (as illustrated, for example, by the simple symmetric exclusion process [1] with mortality), and random walkers on a lattice that can annihilate in pairs. In both these cases the statistics of N is expected to be non-Poissonian, and we calculate the variances describing typical, Gaussian fluctuations of the particle number around the mean. We obtain these results by employing (a dissipative extension of) the Macroscopic Fluctuation Theory (MFT): a coarse-grained low-noise large-deviation theory that employs, as a large parameter, the typical number of particles in the region of interest, see [15] for a recent review. The applicability of the MFT in the driven systems, considered in this work, demands $\bar{N}\gg 1$, and we will work in the parameter regions where this condition is satisfied.

Here is a plan of the remainder of the paper. Section 2 starts with a brief exposition of the expected, or average behavior of a driven system of non-interacting random walkers or Brownian particles that die individually. Then section 2.2 presents the MFT formulation of the problem of particle number statistics for a more general family of diffusive lattice gases of particles that die individually. For the non-interacting particles, we are able to solve, in the same section 2.2, the MFT equations exactly. The solution, via the Hopf–Cole transformation, reproduces (the large-N asymptotic of) the expected Poisson distribution of N. In addition, it gives the previously unavailable optimal density history of the driven system, conditional on observing a given N. Section 3 deals with two examples of interacting particles. The first of them deals with a family of mortal interacting diffusive lattice gases with constant diffusivity but non-trivial fluctuations, as illustrated by the simple symmetric exclusion process with mortal particles. The second example involves random walkers that only interact via pair-wise annihilation. In these examples a full solution of the MFT problem is presently unavailable, and we only calculate the variances of the respective stationary distributions of N. Our main results and their possible extensions are briefly discussed in section 4.

2.1. Expected density behavior

Consider a half-infinite one-dimensional lattice with lattice constant a and suppose that a source of particles at the origin, x = 0, sets a constant particle number density n₀ there. The particles perform random walk at x > 0 with diffusivity D and die individually with rate μ. When D ≫ μa², there is no difference between the discrete random walk and continuous diffusion, and the average particle density ρ(x, t) is governed by the reaction-diffusion equation [1, 16, 17]

$$\begin{equation} \partial_t \rho = -\mu \rho + D\partial_x^2 \rho,\;\;\;\;\;0<x<\infty. \end{equation} \tag{ 1 }$$

At long times the average density profile approaches a steady state, independent of the initial condition:

$$\begin{equation} \bar{\rho}(x) = n_0 {\rm e}^{-\sqrt{\frac{\mu}{D}} x}. \end{equation} \tag{ 2 }$$

Correspondingly, the average steady-state number of particles in this driven system is

$$\begin{equation} \bar{N} =\int_0^{\infty} {\rm d}x\,\bar{\rho}(x) = \sqrt{\frac{D}{\mu}} \,n_0, \end{equation} \tag{ 3 }$$

and we assume that this number is much larger than unity. The actual number of particles N in the steady state fluctuates around $\bar{N}$ because of the shot noises of the diffusion and mortality. We are interested in the steady-state probability distribution of N. For non-interacting random walkers that die individually this probability distribution can be found exactly, by solving the steady-state master equation for the multi-variate probability distribution of observing n₁ particles on site 1, n₂ particles on site 2, etc. The solution has the form of the Poisson product measure with space-dependent parameters (for D ≫ μa², this measure corresponds to the average density profile (2)). This leads to a Poisson distribution of N in the steady state. We will proceed, however, as if we were unaware of these exact results, and employ instead the MFT: a coarse-grained low-noise theory mostly based on the strong inequality $\bar{N}\gg 1$. The purpose is two-fold. First, even for the non-interacting particles, the MFT will give the previously unknown optimal density history of the driven system, conditional on observed N. The optimal density history is well defined only in the limit of $\bar{N}\gg 1$, where it is much more likely than other histories leading to the same N. Second, our main interest is in interacting particle models, where exact microscopic results are usually unavailable. In the next subsection we briefly discuss the basics of the MFT, and formulate the MFT problem for the particle number statistics in a broader context of a family of driven lattice gases of mortal particles.

2.2. Macroscopic fluctuation theory (MFT) of the particle number statistics

2.2.1. Governing equations.

When (i) the length scale of interest is much larger than the lattice constant a, (ii) the time scale of interest is much larger than the inverse rates of the microscopic processes of diffusion and death, and (iii) the typical number of involved particles is much larger than unity, the statistics of large deviations is captured by the MFT [8–10, 15, 18, 19]. The starting point of the derivation of the MFT for diffusing and reacting particles is the exact master equation for the multi-variate probability distribution of observing a certain number of particles on each cite. Here one can either work directly in the physical space, or employ the multi-site probability generating function (that, in the spatially-continuous limit, becomes a probability generating functional). Going over to a path-integral formulation, one then makes a low-noise approximation by evaluating the path integral by the Laplace method that employs the number of particles in the relevant region of space as a large parameter. This procedure yields saddle-point equations (partial differential equations) that can be written in a Hamiltonian form: for the density field q(x, t) and a conjugate field p(x, t) that plays a role of the 'momentum density'. At a qualitative level, the conjugate field p(x, t) describes the magnitude of fluctuations.

If the calculations are performed in the physical space, the saddle point equations (presented here in a form, suitable for a class of diffusive lattice gases of particles that die individually), take the form [8–10, 15, 19]

$$\begin{eqnarray} \partial_t q &=& -\mu q {\rm e}^{-p}+\partial_x\left[D(q) \partial_x q-\sigma(q) \partial_x p\right], \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \partial_t p &=& -\mu ({\rm e}^{-p}-1)- D(q) \partial_x^2 p-\frac{1}{2} \,\sigma^{\prime}(q) (\partial_x p)^2, \end{eqnarray} \tag{ 5 }$$

where D(q) is the gas diffusivity, σ(q) is (twice) the mobility [1], and the prime denotes the derivative with respect to the argument. The terms proportional to μ describe the on-site particle death and its fluctuations. The rest of terms describe diffusive transport and its fluctuations. Equations (4) and (5) are indeed Hamiltonian, as they can be written in terms of variational derivatives:

$$\begin{equation} \partial_t q = \delta H/\delta p\,, \quad \partial_t p = -\delta H/\delta q\,, \end{equation} \tag{ 6 }$$

where

$$\begin{equation} H\{q(x,t),p(x,t)\}= \int_0^{\infty} {\rm d}x\,\mathcal{H} \end{equation} \tag{ 7 }$$

is the Hamiltonian, and

$$\begin{equation} \mathcal{H}(q,p) = \mu q\left({\rm e}^{-p}-1\right) -D(q) \partial_x q \partial_x p +\frac{1}{2}\sigma(q)\!\left(\partial_x p\right)^2 \end{equation} \tag{ 8 }$$

is the Hamiltonian density. Going back to the non-interacting mortal random walkers, we put D(ρ) = D = const and σ(ρ) = 2Dρ, see e.g. [1]. Then equations (4) and (5) become

$$\begin{eqnarray} \partial_t q &=& -\mu q {\rm e}^{-p} + D \partial_x\left(\partial_x q-2 q \partial_x p\right), \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \partial_t p &=& -\mu ({\rm e}^{-p}-1)- D \partial_x^2 p-D (\partial_x p)^2, \end{eqnarray} \tag{ 10 }$$

where 0 < x < ∞. The boundary conditions at the particle source are q(x = 0, t) = n₀ and p(x = 0, t) = 0 [19]. The latter condition is quite intuitive: as we demand a fixed (deterministic) value of the density at x = 0, p must vanish there. Far away from the source there are no particles. This brings the boundary condition q(x = ∞, t) = 0.

Being interested in steady state fluctuations, we can assume that, at t = −∞, the system is at its deterministic steady state [18, 20]: $q(x,t=-\infty)= \bar{\rho}(x)$, see equation (2). We condition the process on observing N particles at some finite moment of time that, without loss of generality, we can set to zero. This imposes an integral constraint on the solution at t = 0:

$$\begin{equation} \int_0^{\infty} {\rm d}x\, q(x,t=0) = N. \end{equation} \tag{ 11 }$$

Analogous integral constraints appear in the MFT formulations of the problem of statistics of integrated current in an infinite setting [21] and statistics of particle absorption by an absorber at x = 0 [22]. To account for the integral constraint, we should introduce a Lagrange multiplier λ and minimize the extended action that incorporates the integral constraint. Similarly to [21] and [22], the action minimization does not change the 'bulk' equations (4) and (5) (or (9) and (10)), but yields an additional boundary condition for p(x, t) at t = 0:

$$\begin{equation} p(x,t=0)=\lambda \, \theta(x), \end{equation} \tag{ 12 }$$

where θ(x) is the Heaviside step function, and λ is ultimately set by equation (11) [21, 22].

Note that p(x, t) = 0 is an invariant manifold of equations (4) and (5). The dynamics on this manifold is described by the deterministic equation (1). This is the relaxation path of the system; it solves the problem in the particular case $N=\bar{N}$. For $N\neq \bar{N}$, the solution of the MFT equations describes the optimal activation path: the most likely density history of the driven system conditional on observing N particles. Here p(x, t) ≠ 0. Once q(x, t) and p(x, t) are found, we can evaluate the action S that yields ${\mathcal P} (N)$ up to a pre-exponential factor:

$$\begin{eqnarray} &-&\ln {\mathcal P}(N) \simeq S = \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x\,\left(p\partial_t q-\mathcal{H}\right) \nonumber\\ &=& \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x\,\left[Dq(\partial_x p)^2+\mu q \left(1-{\rm e}^{-p}-p{\rm e}^{-p}\right)\right]. \end{eqnarray} \tag{ 13 }$$

The first term of the integrand comes from the shot noise of diffusion, the second term comes from the shot noise of mortality. Rescaling time μ t → t, the coordinate $\sqrt{\mu/D}\, x \to x$ and the density q/n₀ → q, one can see that $-\ln {\mathcal P}(N)$ obeys a simple scaling relation

$$\begin{equation} -\ln {\mathcal P}(N) = \bar{N} \,f\left(\frac{N}{\bar{N}}\right), \end{equation} \tag{ 14 }$$

where f(z) is the large deviation function of the number of particles. Note that n₀ only enters this scaling relation through $\bar{N}$.

2.2.2. Particle number statistics and optimal path.

We note that equation (10) is decoupled from equation (9). This decoupling only occurs for non-interacting particles, and it greatly simplifies the problem. Let us perform the Hopf–Cole transformation by introducing Q = qe^−p and P = e^p − 1 [18]. The generating functional of this canonical transformation can be chosen to be

$$\begin{equation} F\{q(x,t),Q(x,t)\}=\int_0^\infty {\rm d}x\, \left(q \ln \frac{q}{Q}-q+Q\right). \end{equation} \tag{ 15 }$$

In the new variables Q and P the Hamiltonian is $\tilde{H}\{Q(x,t),P(x,t)\}=\int {\rm d}x\,\tilde{{\mathcal H}}$, where

$$\begin{equation*} \tilde{{\mathcal H}} = -\mu QP -D\, \partial_x Q \,\partial_x P. \end{equation*}$$

As a result, the MFT equations become linear and fully decoupled:

$$\begin{eqnarray} \partial_t Q &=& -\mu Q + D\partial_x^2 Q, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \partial_t P &=& \mu P - D\partial_x^2 P. \end{eqnarray} \tag{ 17 }$$

Note that these equations for Q and P arise immediately, when one employs the Laplace method for the evaluation of the path integral in the formalism of probability generating functional [18].

In the new variables Q and P, the boundary and initial conditions are:

$$\begin{equation} Q(0,t)=n_0 \;\;\; {\rm and}\;\;\;Q(x,-\infty)=\bar{\rho}(x) \end{equation} \tag{ 18 }$$

for Q, and

$$\begin{equation} P(0,t)=0\;\;\; {\rm and}\;\;\;P(x,0)=\left({\rm e}^{{\lambda}}-1\right)\,\theta(x) \end{equation} \tag{ 19 }$$

for P. As a result, Q(x, t) is invariant in time,

$$\begin{equation} Q(x,-\infty<t\leq 0) = \bar{\rho}(x), \end{equation} \tag{ 20 }$$

while

$$\begin{equation} P(x,-\infty<t\leq 0)=\left({\rm e}^{\lambda}-1\right)\,{\rm e}^{\mu t}\,{\rm erf}\,\left(\frac{x}{\sqrt{-4D t}}\right). \end{equation} \tag{ 21 }$$

Now we can determine the optimal path in the original variable q:

$$\begin{eqnarray} &&\fl q(x,t) = Q(x) [1+P(x,t)] = n_0 {\rm e}^{-\sqrt{\frac{\mu}{D}} x} \left[1+\left({\rm e}^{\lambda}-1\right)\,{\rm e}^{\mu t}\,{\rm erf}\,\left(\frac{x}{\sqrt{-4D t}}\right)\right]. \end{eqnarray} \tag{ 22 }$$

Using equation (11), we find $\lambda=\ln (N/\bar{N})$, so

$$\begin{equation} q(x,t) = n_0 {\rm e}^{-\sqrt{\frac{\mu}{D}} x} \left[1+\left(\frac{N}{\bar{N}}-1\right)\,{\rm e}^{\mu t}\,{\rm erf}\,\left(\frac{x}{\sqrt{-4D t}}\right)\right]. \end{equation} \tag{ 23 }$$

It is easier to calculate the action in the new variables Q and P where, as one can show by a direct calculation [22], the action is equal to the increment of generating functional F from equation (15):

$$\begin{equation} S=F\{q(x,0), Q(x,0)\}-F\{q(x,-\infty), Q(x,-\infty)\}. \end{equation} \tag{ 24 }$$

After some algebra, this gives

$$\begin{eqnarray} - \ln {\mathcal P}(N) \simeq S &=& \int_{0}^{\infty} {\rm d}x\,\left[q(x,0) \ln\frac{q(x,0)}{Q(x,0)}-q(x,0)+Q(x,0)\right] \nonumber\\ &=& N \ln\frac{N}{\bar{N}}-N+\bar{N}, \end{eqnarray} \tag{ 25 }$$

where $\bar{N}$ is given by equation (3). That is, the large deviation function f(z) from equation (14) is equal to f(z) = z lnz − z + 1. The distribution (25) coincides with the N ≫ 1, $\bar{N}\gg 1$ asymptotic of the Poisson distribution with mean $\bar{N}$, as to be expected. In particular, the variance of this distribution coincides with the mean:

$$\begin{equation} V_{{\rm RW}}=\bar{N}=\sqrt{\frac{D}{\mu}}\,n_0. \end{equation} \tag{ 26 }$$

Now let us return to the optimal path (23) that has been previously unknown. Although the statistics of N is time-independent, the optimal path does depend on time. Furthermore, the activation path does not coincide with the time-reversed relaxation path, obtained by solving the deterministic reaction-diffusion equation (1) back in time. This is a clear signature of non-equilibrium. Notice also that, in order to ensure an unusually large or small number of particles at t = 0, the fluctuations create a boundary layer in the density profile at the particle source. This boundary layer becomes a density jump at t = 0,

$$\begin{equation} q(x>0,t=0)= \frac{N}{\bar{N}}\, \bar{\rho}(x), \end{equation} \tag{ 27 }$$

so that the effective boundary condition is $q(x\to 0,t=0)=n_0 N/\bar{N}$, whereas the bulk of the gas particles behaves deterministically. These features can be seen on figure 1 which shows the optimal density histories described by equation (23). The left and right panels correspond to $N=3\bar{N}$ and $N=\bar{N}/3$, respectively. These results are both unexpected and instructive.

Zoom In Zoom Out Reset image size

Now let us consider interacting lattice gases of mortal particles. They may have different σ(q), but for simplicity we will continue to assume a constant diffusivity D = const. For such gases the deterministic equations (1)–(3) continue to hold, while the MFT equations read

$$\begin{eqnarray} &&\partial_t q = -\mu q {\rm e}^{-p} + \partial_x \left[D\partial_x q-\sigma(q) \partial_x p\right], \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} &&\partial_t p = -\mu ({\rm e}^{-p}-1)- D \partial_x^2 p-\frac{1}{2}\sigma^{\prime}(q) (\partial_x p)^2, \end{eqnarray} \tag{ 29 }$$

with the same boundary conditions as before. Once q(x, t) and p(x, t) are known, the probability distribution ${\mathcal P} (N)$ can be evaluated from

$$\begin{eqnarray} \fl &-&\ln {\mathcal P}(N) \simeq S = \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x\,\left[\frac{1}{2}\sigma(q) (\partial_x p)^2+\mu q \left(1-{\rm e}^{-p}-p{\rm e}^{-p}\right)\right]. \end{eqnarray} \tag{ 30 }$$

Rescaling time μ t → t and the coordinate $\sqrt{\mu/D} \,x \to x$, we obtain a scaling relation

$$\begin{equation} -\ln {\mathcal P}(N) = \bar{N} \,f\left(\frac{N}{\bar{N}},n_0\right), \end{equation} \tag{ 31 }$$

where n₀ enters both through $\bar{N}$ and separately¹.

It does not seem possible to solve equations (28) and (29) and determine the large deviation function f(z, n₀) analytically for a general σ(q). Typical, Gaussian fluctuations of the number of particles around $\bar{N}$ are given by a quadratic asymptotic of f(z) at z close to 1. This asymptotic can be found relatively easily via a perturbation theory around the deterministic steady-state solution (2). This theory employs the Lagrange multiplier λ as a small parameter [23]. We set

$$\begin{eqnarray} &&q = \bar{\rho}(x)+ \lambda q_1 + \lambda^2 q_2+\ldots \end{eqnarray} \tag{ 32a }$$

$$\begin{eqnarray} &&p = \lambda p_1 + \lambda^2 p_2+\ldots \end{eqnarray} \tag{ 32b }$$

and plug these expansions into equations (28) and (29). The first-order equations are

$$\begin{eqnarray} &&(\partial_t +\mu -D\partial^2_{x}) \,q_1 = \mu \bar{\rho}(x) p_1-\partial_x[\sigma(\bar{\rho})\partial_x p_1], \end{eqnarray} \tag{ 33a }$$

$$\begin{eqnarray} &&(\partial_t -\mu +D \partial^2_{x}) \,p_1=0. \end{eqnarray} \tag{ 33b }$$

The equation for p₁ is independent of σ(q), and it is decoupled from the equation for q₁. Therefore, we can solve it immediately, with the boundary conditions p₁(0, t) = 0 and p₁(x, 0) = θ(x). The solution is

$$\begin{equation} p_1(x,t\leq 0)={\rm e}^{\mu t} \,{\rm erf}\,\left(\frac{x}{\sqrt{-4Dt}}\right). \end{equation} \tag{ 34 }$$

Now we could plug this expression in equation (33a) and solve for q₁ with the boundary conditions q₁(0, t) = 0 and q₁(x, −∞) = 0. This is unnecessary, however, for the purpose of computing the variance of ${\mathcal P}(N)$, because the latter is independent of q₁. Indeed, we have

$$\begin{eqnarray} &-&\ln {\mathcal P}(N) \simeq S = \frac{\lambda^2}{2} \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x\,\left[\sigma(\bar{\rho}) (\partial_x p_1)^2 +\mu \bar{\rho}p_1^2\right] +{\mathcal O}(\lambda^3). \end{eqnarray} \tag{ 35 }$$

Correspondingly, the variance is equal to

$$\begin{eqnarray} &&V = \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x\,\left[\sigma(\bar{\rho}) (\partial_x p_1)^2 +\mu \bar{\rho}p_1^2\right], \end{eqnarray} \tag{ 36 }$$

with $\bar{\rho}(x)$ from equation (2) and p₁(x, t) from equation (34). Equation (36) does not demand a knowledge of q₁(x, t) and holds for any mortal diffusive lattice gas with D = const.

One well-known example of a gas with constant diffusivity D but non-trivial fluctuations is provided by the SSEP, where each particle can randomly hope to a neighboring lattice site if that site is vacant. If it is occupied by another particle, the move is forbidden. For the SSEP one has σ(ρ) = 2Dρ (1 − ρ) [1]. Here we set the lattice constant a = 1, so that the particle density at the source 0 < n₀ ⩽ 1 is dimensionless. Evaluating the double integral in equation (36) in this case (see appendix A), we obtain

$$\begin{equation} V_{{\rm SSEP}}= \sqrt{\frac{D}{\mu}} \,n_0 \left(1-\frac{2 n_0}{\pi}\right) = \bar{N} \left(1-\frac{2 n_0}{\pi}\right). \end{equation} \tag{ 37 }$$

As $V_{{\rm SSEP}}\neq \bar{N}$, ${\mathcal P}(N)$ is non-Poissonian. As expected on the physical grounds, V_SSEP is smaller than the variance of the total number of non-interacting random walkers with the same n₀, see equation (26), so the distribution is narrower than the Poisson distribution with the same mean. The two variances coincide in the limit of n₀ → 0, where exclusion effects in the SSEP are negligible. That ${\mathcal P}(N)$ is non-Poissonian is not surprising, but even its variance has been previously unknown.

An annihilating random walker (ARW) is immortal when it is alone, but two ARWs on the same lattice site can annihilate, 2A → ∅. Let α be the annihilation rate constant. When the diffusion is sufficiently fast (see the criterion (40) below) the average particle density ρ(x, t) is governed by the continuous reaction-diffusion equation [19]

$$\begin{equation} \partial_t \rho = -\alpha \rho^2 + D\partial_x^2 \rho,\;\;\;\;\;0<x<\infty. \end{equation} \tag{ 38 }$$

The particle source at x = 0 fixes a particle density n₀ of the ARWs there. With this boundary condition, the steady-state average density profile is

$$\begin{equation} \bar{\rho}(x) = n_0 \left(1+\sqrt{\frac{\alpha n_0}{6 D}} \,x\right)^{-2}, \end{equation} \tag{ 39 }$$

it falls off much slower than the exponential profile (2). For the continuous reaction-diffusion equation to be valid, it is necessary that the characteristic length scale ∼(D/α n₀)^1/2 be much larger than the lattice constant a:

$$\begin{equation} \sqrt{\frac{D}{\alpha n_0}}\gg a. \end{equation} \tag{ 40 }$$

The average steady-state number of particles is

$$\begin{equation} \bar{N} =\int_0^{\infty} {\rm d}x\,\bar{\rho}(x)=\left(\frac{6 D n_0}{\alpha}\right)^{1/2}, \end{equation} \tag{ 41 }$$

and we assume $\bar{N}\gg 1$. The MFT equations for this system can be safely derived from the exact master equation for the multi-variate probability distribution by assuming that the typical number of particles on each lattice site is much larger than unity, leading to the strong inequality n₀a ≫ 1 [18, 19]. We believe, however, that it is actually sufficient to require a weaker condition $\bar{N}\gg 1$, alongside with the condition (40). The MFT equations are [19]

$$\begin{eqnarray} \partial_t q &=& - \alpha q^2 {\rm e}^{-2p} + D \partial_x\left(\partial_x q-2 q \partial_x p\right), \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} \partial_t p &=& -\alpha q ({\rm e}^{-2p}-1)- D \partial_x^2 p-D (\partial_x p)^2, \end{eqnarray} \tag{ 43 }$$

whereas

$$\begin{equation} \mathcal{H}(q,p) = \frac{1}{2}\alpha q^2\left({\rm e}^{-2p}-1\right) -D \partial_x q \partial_x p +D q\!\left(\partial_x p\right)^2 \end{equation} \tag{ 44 }$$

is the Hamiltonian density [18, 19]. The first term comes from the on-site annihilations, the second and third terms come from the diffusion. Correspondingly,

$$\begin{eqnarray} \fl &-&\ln {\mathcal P}(N) \simeq S = \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x\,\left[D q (\partial_x p)^2+\frac{\alpha}{2} q^2 \left(1-{\rm e}^{-2p}-2p{\rm e}^{-2p}\right)\right]. \end{eqnarray} \tag{ 45 }$$

As one can check, by performing rescalings described by equation (48) below and additional rescaling q/n₀ → q,

$$\begin{equation} -\ln {\mathcal P}(N) = \bar{N} \,f\left(\frac{N}{\bar{N}}\right). \end{equation} \tag{ 46 }$$

Here n₀ only enters through $\bar{N}$, as for the random walkers who die individually, see equation (14).

As in section 3, we can calculate analytically the variance of the total number of ARWs in the steady state. We make the ansatz (32a) and (32b) in equations (42) and (43) and obtain, in the first order in λ ≪ 1,

$$\begin{eqnarray} &\left[\partial_t +2 \alpha \bar{\rho}(x) -D\partial^2_{x}\right] \,q_1 = 2 \alpha \bar{\rho}^2(x) p_1-2 D \partial_x[\bar{\rho}(x)\partial_x p_1], \end{eqnarray} \tag{ 47a }$$

$$\begin{eqnarray} &\left[\partial_t -2 \alpha \bar{\rho}(x) +D \partial^2_{x}\right] \,p_1=0. \end{eqnarray} \tag{ 47b }$$

As for the mortal SSEP, equation (47b) for p₁ is decoupled from that for q₁, and its solution suffices for computing the variance we are after. Let us reverse and rescale time and transform the coordinate:

$$\begin{equation} \tau= -\frac{\alpha n_0 t}{6} \;\;\;\;\;\;{\rm and}\;\;\;\; y= 1+\sqrt{\frac{\alpha n_0}{6 D}} \,x, \end{equation} \tag{ 48 }$$

so that equation (47b) becomes

$$\begin{equation} \partial_{\tau} p_1 +\frac{12 p_1}{y^2} = \partial^2_{y} p_1. \end{equation} \tag{ 49 }$$

We need to solve it for 1 < y < ∞ and 0 < τ < ∞ subject to the boundary condition p₁(y = 1, τ) = 0 and initial condition

$$\begin{equation} p_1(y,\tau=0)=\theta(y-1) . \end{equation} \tag{ 50 }$$

The rescaled problem for p₁(y, τ) is parameter-free, and we solve it in appendix B. Once p₁(y, τ) is found, we can calculate the variance of N:

$$\begin{eqnarray} V_{{\rm ARW}} &=& \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x\,\left[2 D \bar{\rho} (\partial_x p_1)^2 +2 \alpha \bar{\rho}^2 p_1^2\right] \nonumber\\ &=&2 \bar{N} \int_0^{\infty} {\rm d}\tau\, \int_1^{\infty} \frac{{\rm d}y}{y^2} \left[(\partial_{y} p_1)^2+\frac{6 p_1^2}{y^2}\right]. \end{eqnarray} \tag{ 51 }$$

As expected from equation (46), V_ARW is proportional to $\bar{N}$. Using equation (B10) of appendix B for p(y, τ), we can represent the proportionality coefficient (a dimensionless number of order unity) as a quadruple integral. The integration over τ is elementary. The integration over y is very tedious, but can be performed explicitly with 'Mathematica'. We evaluated the remaining double integral numerically, leading to $V_{{\rm ARW}}\simeq 0.78 \bar{N}$.

As a check, we also solved equation (49) numerically in the region 1 < y < L and 0 < τ < T with the boundary conditions p₁(y = 1, t) = 0 and ∂_yp₁(y = L, t) = 0 and initial condition (50), taking L and T sufficiently large. Then we used the numerical solution to compute the double integral in the second line of equation (51). The result comes quite close, $V_{{\rm ARW}}\simeq 0.77 \bar{N}$. As $V_{{\rm ARW}} <\bar{N}$, the distribution ${\mathcal P}(N)$ is narrower than a Poisson distribution with the same mean.

This work addressed the statistics of the total number of particles N that are present at any chosen time in the steady state of a driven lattice gas composed of mortal diffusing particles. The formalism we used is that of the dissipative Macroscopic Fluctuation Theory (MFT). For non-interacting random walkers who die individually, the MFT formulation of the problem is exactly soluble and yields the expected Poissonian statistics of N with the mean predicted by the simple reaction-diffusion equation (1). It also provides a fascinating and instructive visualization of large deviations of N in the form of the optimal density history of the driven system conditional on N. For interacting diffusing particles we calculated the variance of the distribution of N, and found that the distribution is narrower than a Poisson distribution with the same mean.

The variance calculations that we showed here is a first step towards studying the complete statistics of N. Extending our perturbation theory for the MFT to higher orders in λ, one should be able to compute several higher moments of ${\mathcal P}(N)$, as it has been recently done in the problem of melting of an Ising quadrant [24]. We also note that it is possible to compute the distribution of N numerically by solving the full MFT equations with the proper boundary with the Chernykh-Stepanov iteration algorithm [25]. This algorithm was originally developed for evaluating the probability distribution of large negative velocity gradients in the Burgers turbulence. Later on it was used in studies of different types of large deviations in diffusive lattice gases, with and without on-site reactions [18], [19], [23], [26–29]. This algorithm is much more computationally efficient than microscopic stochastic simulations.

It would be very interesting, and challenging, to directly probe the tails of ${\mathcal P}(N)$, that are beyond the reach of the small-λ perturbation theory. For the SSEP involving immortal particles, μ = 0, the limit of very large transferred mass, in an infinite system, can be described by neglecting the term −D∂_xq ∂_xp in the Hamiltonian density (8). The ensuing reduced MFT equations turn out to be exactly soluble [28, 30]. Whether a similar reduction is possible in the problem of extreme statistics of the total number of interacting mortal particles is an open question.

On a more general note, understanding non-equilibrium systems requires, among other things, intuition which one acquires by learning from examples. The prototypical dissipative systems, considered in this work, are helpful in gaining such an intuition.

This research was supported by grant No. 2012145 from the United States—Israel Binational Science Foundation (BSF).

Here we evaluate the integral in equation (36) for the SSEP. We start with calculating the integral

$$\begin{eqnarray} I_1 &=& \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x\,\sigma(\bar{\rho}) (\partial_x p_1)^2 \nonumber\\ &=& \int_0^{\infty} {\rm d}x \, \int_{-\infty}^0 {\rm d}t\, 2D n_0 \,{\rm e}^{-\sqrt{\frac{\mu}{D}} x} \left(1-n_0\,{\rm e}^{-\sqrt{\frac{\mu}{D}} x}\right) \left(-\frac{{\rm e}^{-\frac{x^2}{2Dt}+2\mu t}}{\pi D t}\right) \nonumber\\ &=& \frac{4n_0}{\pi} \int_0^{\infty} {\rm d}x \,{\rm e}^{-\sqrt{\frac{\mu}{D}} x} \left(1-n_0\,{\rm e}^{-\sqrt{\frac{\mu}{D}} x}\right) {\rm K}_0\left( \sqrt{\frac{4\mu}{D}}x\right), \end{eqnarray} \tag{ A1 }$$

where K₀(...) is the modified Bessel function of the second kind. Evaluating the remaining integral, we obtain

$$\begin{equation} I_1=\sqrt{\frac{D}{\mu}} \,n_0 \left(\frac{4\sqrt{3}}{9}-\frac{2n_0}{\pi}\right). \end{equation} \tag{ A2 }$$

Now we evaluate the integral

$$\begin{eqnarray} I_2 &=& \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x\,\mu \bar{\rho}p_1^2 \nonumber\\ &=& \int_{-\infty}^0 {\rm d}t\, \int_0^{\infty} {\rm d}x \, \mu n_0 \,{\rm e}^{-\sqrt{\frac{\mu}{D}} x+2 \mu t} \,{\rm erf}^2\,\left(\frac{x}{\sqrt{-4Dt}}\right)\nonumber\\ &=& \sqrt{\frac{4D}{\mu}} \,n_0 \int_0^{\infty} {\rm d}z \sqrt{z}\, {\rm e}^{-2z} \int_0^{\infty} {\rm d}u\,{\rm e}^{-2\sqrt{z} \,u} \,{\rm erf}^2\,u. \end{eqnarray} \tag{ A3 }$$

With some patience, this double integral can be evaluated using the book of integrals [31], and the result is²

$$\begin{equation} I_2= \sqrt{\frac{D}{\mu}} \,n_0 \left(1-\frac{4\sqrt{3}}{9}\right). \end{equation} \tag{ A4 }$$

Summing up I₁ and I₂ we obtain equation (37).

To find the spectrum of the problem and the proper eigenfunctions, we make the ansatz p₁(y, τ) = ψ(y, Γ) e^{−Γ2 τ} and arrive at the equation

$$\begin{equation} \frac{{\rm d}^2 \psi(y,\Gamma)}{{\rm d}y^2} +\left(\Gamma^2-\frac{12}{y^2}\right) \psi(y,\Gamma)=0 \end{equation} \tag{ B1 }$$

that we need to solve with the boundary condition ψ(y = 1, Γ) = 0. Equation (B1) is the Shrödinger equation for a quantum particle with energy Γ² in the potential

$$\begin{equation} V(y) = \left\{\begin{array}{cc} \frac{12}{y^2}, & y>1,\\ \infty, & y\leq 1. \end{array}\right. \end{equation} \tag{ B2 }$$

The spectrum of the problem is continuous, 0 < Γ < ∞. Two linearly independent solutions of equation (B1) can be chosen as

$$\begin{equation} \psi_1(y,\Gamma) =\frac{\Gamma y \left(\Gamma ^2 y^2-15\right) \sin (\Gamma y)+3 \left(2 \Gamma ^2 y^2-5\right) \cos (\Gamma y)}{\Gamma ^3 y^3} \end{equation} \tag{ B3 }$$

and

$$\begin{equation} \psi_2 (y,\Gamma) =\frac{\Gamma y \left(\Gamma ^2 y^2-15\right) \cos (\Gamma y)-3 \left(2 \Gamma ^2 y^2-5\right) \sin (\Gamma y)}{\Gamma ^3 y^3}. \end{equation} \tag{ B4 }$$

At Γ → ∞ (or at y → ∞) these solutions become sin (Γy) and cos (Γy) as expected. All the eigenfunctions, vanishing at y = 1, can be written as

$$\begin{equation} \psi(y,\Gamma)=a(\Gamma) \phi(y,\Gamma), \end{equation} \tag{ B5 }$$

where

$$\begin{equation} \phi(y,\Gamma)= \psi_1 (y,\Gamma) \psi_2(1,\Gamma)- \psi_1 (1,\Gamma) \psi_2 (y,\Gamma) \end{equation} \tag{ B6 }$$

and a(Γ) is a yet undetermined amplitude. The eigenfunctions (B5) are orthogonal, and they can be normalized as follows:

$$\begin{equation} \int_1^{\infty} {\rm d}y\, \psi (y,\Gamma) \,\psi (y, \Gamma^{\prime})=\delta (\Gamma-\Gamma^{\prime}), \end{equation} \tag{ B7 }$$

where δ is Dirac's delta function. As a result,

$$\begin{equation} a^{-2}(\Gamma)=\int_0^{\infty} {\rm d}\Gamma^{\prime} \int_1^{\infty} {\rm d}y \,\phi(y,\Gamma)\,\phi(y,\Gamma^{\prime}). \end{equation} \tag{ B8 }$$

Evaluating this double integral with a help of 'Mathematica', we obtain

$$\begin{equation} a(\Gamma)=\sqrt{\frac{2}{\pi}}\,\frac{\Gamma^3}{\sqrt{\Gamma^6+6\Gamma^4+45 \Gamma^2+225}}. \end{equation} \tag{ B9 }$$

The solution for p₁(y, τ) can be written as

$$\begin{equation} p_1(y,\tau)=\int_0^{\infty} {\rm d} \Gamma A(\Gamma) \,\psi (y,\Gamma) \,{\rm e}^{-\Gamma^2 \tau}, \end{equation} \tag{ B10 }$$

where A(Γ) is the projection of the initial condition (50) on the normalized eigenfunctions (B5). That is, $A(\Gamma)=\int_1^{\infty} {\rm d}y\, \psi(y,\Gamma)$. This integral can be also evaluated with 'Mathematica', resulting in a tedious formula

$$\begin{eqnarray*}&&\fl A(\Gamma) = \bigg[4 \Gamma^3+90 \Gamma +6 {\rm Ci}(\Gamma) \left[3 \left(5-2 \Gamma ^2\right) \sin \Gamma +\Gamma \left(\Gamma ^2-15\right) \cos \Gamma \right]\nonumber\\ &&-3 \Gamma \left(\Gamma ^2-15\right) [\pi -2 {\rm Si}(\Gamma )] \sin \Gamma -9 \left(2 \Gamma ^2-5\right) [\pi -2 {\rm Si}(\Gamma )] \cos \Gamma\bigg]\nonumber\\&&\times \bigg[2 \sqrt{2 \pi} \, \Gamma \, \sqrt{\Gamma ^6+6 \Gamma ^4+45 \Gamma ^2+225}\bigg]^{-1}, \end{eqnarray*}$$

where

$$\begin{equation*} {\rm Si}(z) =\int_0^z \frac{\sin y}{y}\,{\rm d}y\;\;\;{\rm and}\;\;\;{\rm Ci}(z) =-\int_z^{\infty} \frac{\cos y}{y}\,{\rm d}y \end{equation*}$$

are the sine and cosine integrals, respectively. Now p₁(y, τ) in equation (B10) is fully determined in terms of a single integral over Γ.