A giant disparity and a dynamical phase transition in large deviations of the time-averaged size of stochastic populations

TL;DR

This paper analyzes large deviations in stochastic populations, revealing a giant disparity in probabilities of extreme population sizes and identifying a second-order dynamical phase transition through advanced WKB methods.

Contribution

It introduces a novel WKB-based approach combined with large deviation formalism to characterize a giant disparity and a dynamical phase transition in population size fluctuations.

Findings

Giant disparity between small and large population deviations.

Identification of a second-order dynamical phase transition.

Finite N smoothing captured by van-Kampen expansion.

Abstract

We study large deviations of the time-averaged size of stochastic populations described by a continuous-time Markov jump process. When the expected population size $N$ in the steady state is large, the large deviation function (LDF) of the time-averaged population size can be evaluated by using a WKB (after Wentzel, Kramers and Brillouin) method, applied directly to the master equation for the Markov process. For a class of models that we identify, the direct WKB method predicts a giant disparity between the probabilities of observing an unusually small and an unusually large values of the time-averaged population size. The disparity results from a qualitative change in the "optimal" trajectory of the underlying classical mechanics problem. The direct WKB method also predicts, in the limit of $N\to \infty$, a singularity of the LDF, which can be interpreted as a second-order dynamical phase transition. The transition is smoothed at finite $N$, but the giant disparity remains. The smoothing effect is captured by the van-Kampen system size expansion of the exact master equation near the attracting fixed point of the underlying deterministic model. We describe the giant disparity at finite $N$ by developing a different variant of WKB method, which is applied in conjunction with the Donsker-Varadhan large-deviation formalism and involves subleading-order calculations in $1/N$.

Figures

Equations

$${\bar{n}}_T=\frac{1}{T}{\int }_0^{T}n\left(t\right)dt,$$

$$-\ln \mathcal{P}({\bar{n}}_T=Na)\simeq TI\left(a,N\right).$$

$$I\left(a,N\right)\simeq Nf(a),N\gg 1.$$

$${\partial }_t{P}_n\left(t\right)$$

$$-\left[W_{+}(n)+W_{-}(n)\right]P_n\left(t\right),$$

$$W_{\pm }(n)=W_{\pm }(Nq)=N{w}_{\pm }(q)+u_{\pm }(q)+...,$$

$$\frac{dq}{dt}=w_{+}(q)-w_{-}(q).$$

$$P_n\left(t\right)=e^{-Ns\left(\frac{n}{N},t\right)}=e^{-Ns\left(q,t\right)}$$

$$s\left(q\pm \frac{1}{N},t\right)=s\left(q,t\right)\pm \frac{1}{N}\frac{\partial s\left(q,t\right)}{\partial q}+\frac{1}{2N^2}\frac{{\partial }^{2}s\left(q,t\right)}{\partial q^2}+\dots .$$

$$\frac{\partial s}{\partial t}+H_0\left(q,\frac{\partial s}{\partial q}\right)=0,$$

$$H_0(q,p)=w_{+}(q)(e^p-1)+w_{-}(q)\left(e^{-p}-1\right)$$

$$s\left(a,T\right)={\int }_0^{T}p\left(t\right)\dot{q}\left(t\right)dt-ET,$$

$$L(q,\dot{q})=L_0(q,\dot{q})-\mu q,$$

$$H(q,p)=p\dot{q}-L(q,\dot{q})=p\dot{q}-L_0(q,\dot{q})+\mu q,$$

$$H(q,p)=H_0(q,p)+\mu q.$$

$$\dot{q}$$

$$\dot{p}$$

$$\frac{1}{T}{\int }_0^{T}q(t)dt=a,$$

$$\frac{\partial H_0(a,p)}{\partial p}=0,$$

$$s\left(a,T\right)=-ET=-H_0\left[a,p\left(a\right)\right]T.$$

$$f(a)=\frac{s\left(a,T\right)}{T}=-H_0\left[a,p\left(a\right)\right].$$

$$p=\frac{1}{2}\ln \frac{w_{-}(a)}{w_{+}(a)}.$$

$$f(a)={\left[\sqrt{w_{+}(a)}-\sqrt{w_{-}(a)}\right]}^2,$$

$$\varnothing \stackrel{N}{\to }A,A\stackrel{1}{\to }\varnothing$$

$${\pi }_n{W}_{+}(n)={\pi }_{n+1}{W}_{-}(n+1),$$

$${\pi }_n=\frac{e^{-N}{N}^n}{n!},n=0,1,\dots .$$

$$H(q,p,\mu )=e^p-1+q(e^{-p}-1)+\mu q.$$

$$f(a)={\left(1-\sqrt{a}\right)}^2.$$

$$A\stackrel{1}{\to }2A,2A\stackrel{2/N}{\to }A,$$

$$W_{+}(n)=n,W_{-}(n)=\frac{n\left(n-1\right)}{N}.$$

$${\pi }_n=\frac{N^n}{\left(e^N-1\right)n!},n=1,2,\dots .$$

Full text

I Introduction

Stochastic population models describe a multitude of processes in nature and society Delbrueck ; Bartlett ; Karlin ; Bailey ; Nisbet ; Andersson ; Allen ; Tsimring ; Allen2015 . Examples include reactions among particles, dynamics of biological populations, spread of diseases, inventions, opinions and rumors, dynamics of financial markets, and many other processes. A convenient mathematical description of stochastic populations is provided by Markov jump processes Kampen ; Gardiner . An important class of these processes has a nontrivial steady state, where the gain and loss processes balance each other on average. A standard characterization of a fluctuating steady state is in terms of its steady-state probability distribution. For a single population in a steady state this distribution, $\pi_{n}$, describes the probabilities to observe $n=0,1,2,\dots$ individuals, which we will call particles.

In recent years there has been a growing interest in a different characterization of fluctuating steady states: by the probability distributions of long-time averages of observables, see Refs. Touchette2009 ; Touchette2018 and references therein. In the context of stochastic populations a natural quantity of this type is the long-time average of the (fluctuating) population size. For many stochastic systems, the probability of observing any value of the time-averaged size, which is different from the expected one, becomes exponentially small as the averaging time $T$ increases.

A similar property, and the ensuing large deviation function (LDF) that encodes it, have received much recent attention from physicists in the context of continuous diffusion processes, described by Langevin equations Touchette2009 ; Touchette2018 ; TouchettePRL2018 ; Lecomte . Recently this characterization has been extended MeersonZilber to one of the best known jump processes in physics: the Ehrenfest Urn Model EUM . The basic version of the Ehrenfest Urn Model involves two urns with $N$ identical balls. The balls hop randomly, one ball at a time, from one urn to the other. For non-interacting balls the LDF of the long-time average of the number of balls in a given urn can be calculated exactly MeersonZilber by using the Donsker-Varadhan (DV) large-deviation formalism DonskerVaradhan , which we briefly review below. If the balls interact, the DV formalism can still be used numerically MeersonZilber , but it does not allow for analytical solution. However, when the total number of balls $N$ is very large, one can circumvent the DV method and employ instead a WKB approximation, directly applied to the master equation Kubo ; Gang ; Dykman ; AssafMeerson . In this way one obtains a controlled large-$N$ approximation of the LDF MeersonZilber .

Here we follow this line of work and identify a class of stochastic populations, for which the direct WKB method predicts a giant disparity between the probabilities of observing an unusually small and an unusually large values of the time-averaged population size. This disparity results from a qualitative change of the character of the optimal (least-action) trajectory of an underlying Hamiltonian classical mechanics problem. The direct WKB method also predicts a jump in the second derivative of the LDF, which can be interpreted as a second-order dynamical phase transition. The sharp transition is smoothed by finite-$N$ effects, unaccounted for by the leading-order WKB approximation, but the giant disparity remains. We show it by developing a different version of WKB theory, which involves subleading-order WKB calculations and is used in conjunction with the DV large-deviation formalism.

Now let us be more specific. Consider a jump process which admits a nontrivial steady state, and denote the population size (the number of particles) in a specific realization of the process at time $t$ by $n\left(t\right)$. We are interested in the probability, $\mathcal{P}(\bar{n}_{T}=Na)$, that the time averaged population size $\bar{n}_{T}$, defined as

[TABLE]

is equal to $Na$, where $N\gg 1$ is the expected population size, and $a$ is a real non-negative number. When $T$ is very large, $\mathcal{P}(\bar{n}_{T}=Na)$ satisfies the large deviations principle, which states that any deviation from the expected value is exponentially unlikely to occur:

[TABLE]

The rate function $I\left(a,N\right)$ is the focus of our interest. It is a non-negative function of $a$ that vanishes at its most probable value which, for $N\gg 1$, is $a\simeq 1$. In the limit $T\rightarrow\infty$, $I\left(a,N\right)$ is independent of the initial condition. An analytical calculation of $I\left(a,N\right)$ is in general impossible. When $N\gg 1$, one can employ a direct WKB approximation Kubo ; Gang ; Dykman ; AssafMeerson , as it has been recently done for the Ehrenfest Urn Model MeersonZilber . In this approximation the $a$- and $N$-dependences of the rate function $I\left(a,N\right)$ factorize MeersonZilber . In the class of models that we will work with here, the rate function $I\left(a,N\right)$, as predicted by the direct WKB approximation, is proportional to $N$:

[TABLE]

It is natural to call the function $f(a)$ the intensive rate function.

For many stochastic populations the intensive rate function $f(a)$ is an analytic function of $a$, and the Ehrenfest Urn Model provides one such example MeersonZilber . As we will see shortly, another such example is provided by the immigration-death process $\emptyset\rightleftarrows A$, a simple and well studied reversible jump process. However, as we find here, there is a class of population models where $f(a)$, as described by the direct WKB method, vanishes on the whole interval $0\leq a\leq 1$. This corresponds to a giant disparity between the probabilities of observing an unusually small and unusually large values of $a$. Furthermore, the predicted $f(a)$ is non-analytic at $a=1$. We will identify this class of models and exemplify it by still another simple reversible process: the branching-coalescence process $A\rightleftarrows 2A$. Then we will obtain more accurate results for this model, by going beyond the direct WKB method.

Here is a plan of the remainder of the paper. In Sec. II we outline the direct WKB method and consider two simple examples: the immigration-death process and the branching-coalescence process. We show that the intensive rate function $f(a)$ is analytic in the former case, and non-analytic and vanishing at $0\leq a\leq 1$ in the latter. In Sec. III we apply WKB approximation, in the leading and subleading orders with respect to $1/N$, to the DV method. In the leading order this calculation reproduces the results of the direct WKB method. In the subleading order the $f(a)=0$ degeneracy is removed due to finite-$1/N$ effects. We also compute the probability distribution $\pi_{n}(a)$ of the population size, conditional on a specified value of $a$. We show that, in the phase-transition region of $a$, the distribution $\pi_{n}(a)$ is bimodal, which reflects the character of the optimal trajectories of the system in this region. We discuss our results and some possible generalizations in Sec. IV.

II Direct WKB method

II.1 General

We start this section by considering a general single-step jump process for a single population. A generalization to multiple-step processes is immediate Dykman ; EscuderoK ; AM2010 . The process is governed by the master equation Kampen ; Gardiner

[TABLE]

for the evolution of the probability $P_{n}\left(t\right)$ of observing $n$ particles at time $t$. The first two terms on the right describe the incoming rates into state $n$ from the states $n-1$ and $n+1$, whereas the last term describes the total outgoing rate from state $n$. We assume throughout this work that the birth and death rates $W_{+}(n)$ and $W_{-}(n)$ are such that Eq. (II.1) admits a nontrivial steady-state solution. We will also assume that the birth and death rates are disparate so that the expected population size $N$ in the steady state is large, and exploit the large parameter $N\gg 1$. Let us expand the process rates $W_{\pm}$ in the powers of $N$. Introducing the rescaled population size $q=n/N$, we obtain EscuderoK ; AM2010

[TABLE]

where $w_{\pm}(q)$ and $u_{\pm}(q)$ are assumed to be $O(1)$.

If one multiplies the master equation (II.1) by $n$, sums it up over all possible values of $n$ and ignores fluctuations Kampen ; Gardiner , one arrives at the deterministic rate equation

[TABLE]

The stationary population, residing in a vicinity of $n=N$, corresponds to an attracting fixed point of Eq. (6) at $q=1$, so that $w_{+}(1)-w_{-}(1)=0$ and $w_{+}^{\prime}(1)-w_{-}^{\prime}(1)<0$, where prime denotes the $q$-derivative.

The direct WKB method for an analysis of fluctuations Kubo ; Gang ; Dykman ; AssafMeerson exploits the large parameter $N\gg 1$ by applying the WKB ansatz

[TABLE]

directly to the master equation. Additionally, we assume $n\gg 1$, treat $q=n/N$ as a continuous variable, and evaluate the action $s$ at $\frac{n\pm 1}{N}=q\pm\frac{1}{N}$ by Taylor expanding it around $q$ Kubo ; Gang ; Dykman :

[TABLE]

Using this expansion of $s$ in the WKB ansatz (7) and plugging the ansatz into the master equation (II.1), we arrive at an equation which involves $q$, $\partial_{q}s$ and $\partial_{t}s$ in different orders of $1/N$, and treat it order by order. In the leading order this procedure yields a Hamilton-Jacobi equation for the “action” $s(q,t)$ Kubo ; Gang ; Dykman :

[TABLE]

where

[TABLE]

is the model-specific Hamiltonian. As one can see, the population size $q$ plays the role of the coordinate of an effective “particle”, whereas $p=\partial_{q}s$ is the momentum. The action along a trajectory $q(t),p(t)$ is given by

[TABLE]

where $E=H_{0}\left[q\left(t\right),p\left(t\right)\right]$ is the (conserved) energy of the “particle”. The constraint $n_{T}=Na\,$ ($a\geq 0$) can be accounted for by adding a term $-\mu\int_{0}^{T}q\left(t\right)dt$, where $\mu$ is a Lagrange multiplier, to the action. In the Lagrangian formulation the constrained Lagrangian is

[TABLE]

where $L_{0}(q,\dot{q})$ is the unconstrained Lagrangian, corresponding to the unconstrained Hamiltonian $H_{0}(q,p)$. The constrained Hamiltonian is equal to LLMechanics

[TABLE]

where $\dot{q}$ should be expressed through $q$ and $p$ from the relation $p=\partial L/\partial{\dot{q}}$. As a result, the constrained Hamiltonian takes the form MeersonZilber

[TABLE]

As one can see, the constraint is enforced via a simple additive term $\mu q$ in the Hamiltonian. This simplicity occurs because $\partial L/\partial{\dot{q}}$ does not depend on $\mu$, see Eq. (12). The phase-space trajectories $q(t)$ and $p(t)$ satisfy the Hamilton’s equations

[TABLE]

Among all such trajectories, we must choose the one which satisfies appropriate boundary conditions and the condition

[TABLE]

which ultimately sets the value of $\mu$. If there is more than one such trajectory, we must choose the “optimal” one which minimizes the action (11). Having found it, we can evaluate $-\ln\mathcal{P}\left(\bar{n}_{T}=Na\right)\simeq Ns\left(a,T\right)$, where $s\left(a,T\right)$ is its corresponding action.

In the Hamiltonian description the attracting fixed point $q=1$ of the deterministic rate equation (6) becomes a saddle point $(q=1,\,p=0)$ of the unconstrained Hamiltonian $H_{0}$. For the constrained Hamiltonian (14) this saddle point is shifted. For finite $T$ the action $s\left(a,T\right)$ depends on the initial condition $q\left(t=0\right)$, and must be minimized with respect to the final condition $q\left(t=T\right)$. As $T$ increases, the effective “particle” spends more and more time in a close vicinity of the above-mentioned shifted saddle point of the constrained Hamiltonian MeersonZilber . For $T\rightarrow\infty$ this saddle point dominates both the action (11) (where the first term vanishes) and the integral over time in Eq. (17). Then, in view of the Eq. (17), the $q$-value of this fixed point must be equal to $a$. The corresponding value of $p=p(a)$ can be obtained from the algebraic equation

[TABLE]

which follows from Eq. (15). At this stage we are already constraining the process by $a$, so the Lagrange multiplier $\mu$ becomes unnecessary. The optimal action, evaluated at the saddle point, is very simple:

[TABLE]

Finally, combining Eqs. (2), (7) and (19), we obtain the intensive rate function $f(a)$, defined by Eq. (3):

[TABLE]

Equation (20) has quite a general form, and it also holds for multi-step processes. Let us specialize Eq. (20) to the single-step process (II.1). Using Eqs. (10) and (18), we obtain the value of $p$ at the fixed point:

[TABLE]

As a result, Eq. (20) yields

[TABLE]

so one only needs to know the leading-order terms of the birth and death rates (5). Equation (22) reproduces the recent result MeersonZilber for the rate function of the time-averaged number of balls in one of the urns in the Ehrenfest Urn Model.

II.2 $0\rightleftarrows A$

As another simple example, consider the immigration-death process

[TABLE]

with the per-capita immigration rate $N$ and the per-capita death rate $1$. The process rates are $W_{+}(n)=N$ and $W_{-}(n)=n$. This simple process obeys the detailed balance condition

[TABLE]

where the equilibrium distribution $\pi_{n}$ is Poisson:

[TABLE]

Now we assume $N\gg 1$ and $n\gg 1$ and turn to the WKB description. We immediately note that the leading-order rates $w_{+}(q)=1$ and $w_{-}(q)=q$ are exact here. Figure 1 shows the phase portrait of the constrained Hamiltonian,

[TABLE]

The shifted saddle point is at $\left(q,p\right)=\left[\left(1-\mu\right)^{-2}\!,-\ln\left(1-\mu\right)\right]$. It exists for $-\infty<\mu<1$, and its coordinate $q$ can take any value from zero to infinity.

Plugging the rates $w_{+}(q)=1$ and $w_{-}(q)=q$ into Eq. (22), we obtain the intensive rate function of the long-time average population size in this system:

[TABLE]

The function $f(a)$ is strictly convex. It vanishes at the expected value $a=1$ and is positive for all other $a$. It is shown in Fig. 2 alongside with numerical results, obtained by the DV method DonskerVaradhan , which we briefly discuss in the beginning of Sec. III. The perfect agreement, even for a moderate value of $N=10$, is explained by the fact that, for this simple process, the WKB result (27) is exact. The WKB result is also exact for the Ehrenfest Urn Model with non-interacting balls MeersonZilber .

II.3 $A\rightleftarrows 2A$: Giant disparity of probability and dynamical phase transition at $N\to\infty$

Expressions (20) and, for single-step processes, (22) are not always correct. Indeed, the shifted saddle point with coordinate $a$ may not exist, or it may be unaccessible for some initial conditions. But even if it exists and is accessible, there can be a different trajectory with a lesser action, which does not stay in a close vicinity of the shifted saddle point. We will illustrate these general arguments on the example of a well-known branching-coalescence process

[TABLE]

with per-capita rates $1$ and $2/N$, respectively. This process exhibits nontrivial dynamics only if it starts from a nonzero number of particles, and we will suppose that this is the case. The branching rate $W_{+}(n)$ is given by the product of the per-capita branching rate $1$ and the current population size $n$. The coalescence rate $W_{-}(n)$ is given by the product of the per-capita coalescence rate $2/N$ and the number of pairs in the population of $n$ particles, $n\left(n-1\right)/2$. As a result,

[TABLE]

This process also obeys the detailed balance condition (24), and the equilibrium distribution is again Poisson:

[TABLE]

However, as we will see now, the large deviation function of the long-time average population size is more interesting here than for the immigration-death process. Applying the direct WKB method and using Eqs. (5) and (29), we identify the leading-order rates $w_{+}(q)=q$ and $w_{-}(q)=q^{2}$. Recklessly using Eq. (22), we would obtain the intensive rate function

[TABLE]

shown in Fig. 3 by the dotted line. This function is obviously non-convex for $a<9/16$. Actually, it is non-convex, according to the supporting line definition of convexity, for all $0<a<1$. This arises suspicion, and indeed Eq. (31) turns out to be incorrect on the whole interval $0<a<1$.

In order to see why, let us explore the phase portrait of the system. The deterministic equation (6) takes the form $\dot{q}=q-q^{2}$. The attracting fixed point is again $q=1$, but now there is also a repelling fixed point at $q=0$. This fact turns out to be crucial, as the effective mechanical system can now stay for a long time at the fixed point $q=0$ without the need to “spend” any action. Let us see how it works. The phase portraits of the unconstrained Hamiltonian,

[TABLE]

is shown in Fig. 4b. The repelling and attracting fixed points of the deterministic dynamics become saddle points $(q,p)=(0,0)$ and $(q,p)=(1,0)$, respectively. In addition, an elliptic fixed point appears at $(q,p)=(1/4,-\ln 2)$.

For the flow described by the constrained Hamiltonian,

[TABLE]

the three fixed points become $\mu$-dependent. The saddle point $(0,0)$ of the unconstrained Hamiltonian becomes $[0,\ln(1-\mu)]$. It exists for $\mu<1$. The nontrivial saddle point exists for $\mu>-1/8$, and the elliptic point exists for $-1/8<\mu<1$. At $\mu=-1/8$ the latter two fixed points merge, and for $\mu<-1/8$ they disappear: compare panels c and d in Fig. 4. As one can check, for all $\mu\geq-1/8$, the value of $q$ of the nontrivial saddle point is greater than or equal to $9/16$; the equality occurs at $\mu=-1/8$. The absence of a nontrivial saddle point with $q<9/16$ (see Fig. 4d) explains why Eq. (31) cannot be correct for $a<9/16$.

What happens for $9/16<a<1$? Here the nontrivial saddle point is accessible for some initial conditions, and inaccessible for others. More importantly, there is a trajectory of a different type, which corresponds to $\mu=0$, provides any specified value of $a\in(0,1)$ and has a lesser action. This trajectory is closely related to the exact one-parameter family of instanton solutions of the unconstrained Hamilton equations:

[TABLE]

These solutions, depicted in Fig. 5, satisfy the boundary conditions $q\left(t=-\infty,C\right)=1$, $q\left(t=\infty,C\right)=0$, $p\left(t=-\infty,C\right)=0$ and $p\left(t=\infty,C\right)=-\infty$. They correspond to the zero-energy phase trajectory $p=\ln q$, indicated by the thick line on Fig. 4b. When $T\gg 1$, they can serve as approximate solutions on the finite interval $0<t<T$, and they are accurate up to corrections which are exponentially small with respect to $T\gg 1$. Furthermore, one can always choose the constant $C$ so that condition (17) is satisfied for any $0<a<1$. The action (11) for each of these trajectories is equal (again, up to exponentially small corrections in $T$) to

[TABLE]

Crucially, this action is independent of $T$. As a result, the intensive rate function $f(a)=s/T$ vanishes at $T\to\infty$ and $0<a<1$, and we arrive at

[TABLE]

Therefore, the probabilities of observing $a<1$ and $a>1$ are dramatically different. Notably, $f(a)$ is continuous with its first derivative $f^{\prime}(a)$ at $a=1$. The second derivative, however, has a jump at $a=1$ which can be interpreted as a dynamical phase transition of second order. The appearance of a giant disparity in $f(a)$, and of a second-order phase transition, in a simple reversible one-population model is quite remarkable.

The sharp transition, however, is observed only in the limit $N\to\infty$. A valuable insight into finite-$N$ effects can be achieved if we go beyond the WKB method and try to interpret the unexpected result $f(a)=0$ at $0\leq a\leq 1$ in the language of the exact microscopic process $A\rightleftarrows 2A$. We can easily do it for one particular value of $n$: $n=1$. What is the probability to observe, for very large $T$, the average population size equal to $1$ fn_nmin ? The only stochastic trajectory that can contribute to this average population size is the one where $n\left(t\right)=1$, while the branching process $A\to 2A$ is suppressed for the whole time $T\gg 1$. The probability of the stochastic trajectory $n\left(t\right)=1$ during the time $T$ is $e^{-T}$, hence $I\left(a=1/N,N\right)=1$. Consequently, in the limit $N\gg 1$, we get $f\left(a=0\right)=0$. On the other hand, $I\left(a,N\right)$ vanishes at $a\simeq 1$ (for $N\gg 1$), so that $f\left(a=1\right)=0$. As we expect the rate function to be convex, it must satisfy $f\left(0\leq a\leq 1\right)=0$, in agreement with the first line of Eq. (36), predicted by the direct WKB method.

The direct-WKB prediction (36) for $f(a)$ is presented in Fig. 3 alongside with numerical results obtained by the DV method DonskerVaradhan for $N=100$. One can see excellent agreement for $a>1$ (which is not surprising, as in this regime the solution comes from the shifted saddle point). The numerical result at $a=1/N$ perfectly agrees with the exact result $I\left(a=1/N,N\right)/N=1/N$, quoted above. It follows immediately from $I\left(a=1/N,N\right)=1$, $I\left(a\simeq 1,N\right)=0$ and the expected convexity of $I(a)$ that

[TABLE]

The numerics indeed indicate a linear behavior of $I\left(a,N\right)/N$ as a function of $a$, with slope $-1/N$, so that the upper bound (37) is in fact very close to the actual value of $I\left(a,N\right)/N$. Our next goal will be to derive a subleading WKB asymptotic which describes the finite $N$ effects observed in the numerics.

III WKB for DV

III.1 DV method

The exact DV method DonskerVaradhan reduces the determination of the rate function $I\left(a,N\right)$ in the $T\rightarrow\infty$ limit to the solution of an eigenvalue problem. In most cases this eigenvalue problem can only be solved numerically. Figures 2 and 3 show such numerical solutions which we obtained for the processes $0\rightleftarrows A$ and $A\rightleftarrows 2A$, respectively. In this section we will show how to solve this eigenvalue problem analytically, albeit approximately, by exploiting the large parameter $N\gg 1$. We will do it by developing a different version of WKB theory, which is applied to the DV eigenvalue problem.

We will start with a brief description of the DV method. Consider a general jump process, governed by a master equation of the form

[TABLE]

For a single-step process (II.1) $\hat{L}^{+}$ is a tridiagonal matrix. The Gärtner-Ellis GartnerEllis theorem states that the rate function $I\left(a,N\right)$ is given by the Legendre-Fenchel transform,

[TABLE]

of the scaled cumulant generating function $\lambda\left(k\right)$, defined as

[TABLE]

where $\braket{...}$ denotes averaging over possible values of $a$. According to DV, $\lambda\left(k,N\right)$ is equal to the largest eigenvalue of an auxiliary operator $\hat{L}^{k}$:

[TABLE]

where

[TABLE]

is a tilted version of the operator $\hat{L}$, which is the Hermitian conjugate of $\hat{L}^{+}$. For any value of $a$, the Legendre-Fenchel transform (39) associates a value of $k$. We will denote it by $k^{*}\left(a,N\right)$, so that

[TABLE]

We used (a truncated version of) this method to perform our numerical calculations.

III.2 $1/N$ expansion and WKB for DV

Employing the small parameter $1/N$, we can seek the rate function perturbatively:

[TABLE]

In the absence of the phase transition in the limit of $N\to\infty$, both $I_{0}(a)$ and $I_{1}(a)$ are of order $1$. For the branching-coalescence process (and other models which exhibits a phase transition of this type) $I_{0}(a)$ and $I_{1}(a)$ are of order $1$ only for $a>1$. For $a<1$ we should expect $I_{0}(a)=0$, so the rate function is

[TABLE]

It is this, more interesting case, that we will focus on lessinteresting . We additionally assume that $a\gg 1/N$, as the WKB approximation does not hold for $n=O(1)$. As we will see, the WKB approximation also breaks down for $a$ very close to $1$, and we will use a different approximation there.

By virtue of Eqs. (43) and (45), $k^{*}\left(a,N\right)$ must be of order $1/N$, and $\lambda\left[k^{*}\left(a,N\right),N\right]$ must be of order $1$, so we should set

[TABLE]

where the unknown functions $k_{1}^{*}(a)$ and $\Lambda(k_{1})$ are of order 1. The last equality in Eq. (47), where we introduce a new function $\Lambda\left[k_{1}^{*}\left(a\right)\right]$, follows from the demand that $\lambda\left[k_{1}^{*}(a)/N,N\right]=O(1)$ does not depend on $N$. Plugging Eqs. (46) and (47) into Eq. (39) yields

[TABLE]

Let us determine the function $\Lambda(k_{1})$, specializing our calculations to the branching-coalescence model. The auxiliary operator $\hat{L}^{k}$ of the DV method is

[TABLE]

where we have set $k\simeq k_{1}/N$, see Eq. (46). Let us consider the eigenvalue problem for the right eigenvector $r_{n}$,

[TABLE]

and find the maximal eigenvalue $\xi=N\xi_{0}+\xi_{1}+\dots$ by using WKB approximation. As in the direct WKB method, we define $q=n/N$, make a WKB ansatz

[TABLE]

and treat $q$ as a continuous variable, so that

[TABLE]

When we plug this ansatz into Eq. (50), the factor $e^{NS(q)}$ will cancel out, so we will ignore it from now on. Before making the ansatz, we expand $S(q)$ in the powers of $1/N$:

[TABLE]

so that

[TABLE]

Now we make the ansatz and obtain

[TABLE]

In the $O(N)$ order we obtain a time-independent Hamilton-Jacobi equation

[TABLE]

with the Hamiltonian $H_{0}\left(q,p_{0}\right)$, defined in Eq. (32). The quantity $\xi_{0}$ plays the role of the (conserved) energy. We immediately notice that $H_{0}\left(q,p_{0}\right)$, and therefore $\xi_{0}$, do not depend on $k$, so that we can evaluate them at $k=0$. But for $k=0$, the right eigenvector, which corresponds to the maximal eigenvalue, is a constant which can be set to be $r_{k=0}=1$ Touchette2018 . As a result, $p_{0}(q)\equiv S_{0}^{\prime}\left(q\right)=0$, and Eq. (56) yields $\xi_{0}=\xi_{0}^{(k=0)}=0$. This result coincides with the one obtained by the direct WKB method, where we found that the energy of the optimal trajectory is zero for any $0<a\leq 1$ fn_eigenvector .

The sub-leading correction emerges from the terms $O(1)$ of Eq. (III.2). Using the equalities $p_{0}(q)=p_{0}^{\prime}(q)=0$, we obtain

[TABLE]

This is again a time-independent Hamilton-Jacobi equation, with the Hamiltonian $H_{1}\left(q,p_{1}\right)=q\left(1-q\right)p_{1}+k_{1}q$. Solving Eq. (57) for $S_{1}^{\prime}(q)$, we see that $S_{1}^{\prime}(q)$ diverges at $q=1$ unless we set $\xi_{1}=k_{1}$. As a result, $S_{1}(q)=k_{1}\ln q$ fn_subleading . The divergence of $S_{1}(q)$ at $q=0$ is not dangerous, because at $q=O(1/N)$, that is at $n=O(1)$, our continuous WKB approximation for $r_{n}$ breaks down anyway.

By virtue of $\xi_{1}=k_{1}$ and Eqs. (41) and (47), we obtain

[TABLE]

Before we use this result in Eq. (48), we will show that $k_{1}$ is bounded from below by $-1$. Let us return to the exact eigenvalue problem (50):

[TABLE]

In particular, for $n=1$ we obtain

[TABLE]

All components of the eigenvector, corresponding to the maximum eigenvalue, and in particular $r_{1}$ and $r_{2}$, must have the same sign Touchette2015 ; Touchette2018 . (A familiar analog is a ground-state wave function which does not have nodes.) Therefore, the right hand side of Eq. (60) must be positive. Using the scalings $k\simeq k_{1}/N$ and $\xi\simeq\xi_{1}$ [see Eqs. (41), (46) and (47)], we obtain

[TABLE]

To avoid excess of accuracy, we ignore the term $-k_{1}/N$ and obtain the inequality $\xi_{1}>-1$. Applying now the WKB result (58), we see that $k_{1}>-1$. To remind the reader, this result is valid only in the WKB regime. Using this bound and Eq. (58) in Eq. (48), we obtain

[TABLE]

that is the maximum is reached at the boundary. The nonzero rate function (62) removes the degeneracy $I\left(a,N\right)=0$ at $0\leq a\leq 1$, that we previously observed at $N\to\infty$.

The continuous WKB theory breaks down at $a=O\left(1/N\right)$, so Eq. (62) is inapplicable there. It is also inapplicable too close to $a=1$, where the WKB action is small. Fortunately, in a close vicinity of $a=1$ one can approximate AssafMeerson ; MeersonZilber the exact master equation of the jump process by a continuous Fokker-Planck equation, using the Van Kampen system-size expansion Kampen . This leads to an effective Ornstein-Uhlenbeck process, and we obtain (see e.g. Eq. (49) in Ref. Touchette2018 )

[TABLE]

corresponding to the Gaussian asymptotic of $\mathcal{P}\left(\bar{n}_{T}=a\right)$ as a function of $a$ near the expected value $a=1$. Equation (63) also follows from a Taylor expansion around $a=1$ of the $a>1$ branch of $I(a)$ from Eq. (36). Notably, the smooth Gaussian asymptotic clearly implies that the second-order phase transition is an attribute of the leading-order WKB approximation. Using Eqs. (62) and (63), we can establish the applicability condition of Eq. (62) at $a$ close to $1$: $1-a\gg 1/N$. Putting everything together, we obtain

[TABLE]

The asymptotics (66) and (67) coincide in their joint region $0\leq a-1\ll 1$.

Our numerics, see Fig. 6, show that the exact rate function is strictly convex. This feature is most clearly observed at not too large $N$ or at small $a$. The approximate expression (64) is not strictly convex, but it is natural to assume that strict convexity will appear in the next order of the perturbation procedure minus_infinity . Figure 6 compares the asymptotics (64)-(67) with numerical results, obtained by the DV method DonskerVaradhan , for $N=20$ and $100$. As one can see, the WKB approximation improves for larger $N$.

Figure 7 shows our numerical results for the stationary distribution $\pi_{n}(a)$ of the branching-annihilation process, conditioned on a given value of $a$ (or $k$). This distribution is given by the properly normalized product of the right and left eigenvectors, $\pi_{n}=r_{n}l_{n}$ for a fixed value of $a$ (or $k$) Chetrite2013 ; Touchette2015 . Remarkably (but a posteriori not surprisingly), for $1/N<a<1$ the conditional distribution $\pi_{n}$ is bimodal: one of its two peaks is at $n\simeq N$, the other is at $n=1$. As $N$ increases, the two peaks become well separated, and the area under the peak at $n=N$ approaches $a$, whereas the area under the peak at $n=1$ approaches $1-a$. This is in a remarkable agreement with the prediction of our direct WKB analysis of Sec. II.3. Indeed, for $a<1$ the optimal trajectory is such that the population has size $N$ for time $aT$ and then, after a short transient, size $1/N$ (which, in the direct WKB method, is indistinguishable from zero) for time $(1-a)T$.

IV Discussion

Time-averaged quantities provide a useful characterization of stochastic populations. Here, using a WKB approximation, we observed, at $N\to\infty$, a dynamical phase transition and a giant disparity in the large deviation function (LDF) of the time-averaged population size of the branching-coalescence process. These effects result from a qualitative change in the “optimal” trajectory of the underlying Hamiltonian classical mechanics problem. The phase transition is in fact smoothed at finite $N$, but the probabilities of observing unusually small and unusually large values of the time-averaged population size remain dramatically different. In order to resolve the degeneracy $f(a)=0$ at $0\leq a\leq 1$, obtained in the leading WKB order, we developed a different version of WKB method, which is applied in conjunction with the Donsker-Varadhan large-deviation formalism and goes beyond the leading-order calculations in $1/N$. A similar in spirit methodology has been recently implemented for a continuous diffusion process Derrida , where a (very different) singularity of a large deviation function, and its regularization, have been studied.

From a more general perspective, our work provides an additional example of a singularity in a LDF of a time-integrated fluctuating quantity which appears in the weak-noise limit $N\to\infty$. A somewhat similar in spirit singularity (also of second order) appears in the LDF of the time- and space-integrated current in certain diffusive lattice gases in systems with periodic boundaries Bertini2006 ; Bodineau2007 . The details, however, are quite different.

Our results, obtained for the branching-annihilation model, can be immediately generalized to a whole family of single-species stochastic population models, which have a nontrivial (that is, non-empty) steady state, corresponding to a non-trivial attracting fixed point of deterministic theory. A sufficient condition for the emergence of a giant probability disparity at large $N$, and a second-order dynamical phase transition in the limit of $N\to\infty$, is the presence of a repelling fixed point at $q=0$. We can even relax the demand of the non-empty steady state by conditioning the process on the population survival by the time $T$ over which the long-time average is calculated.

It would be interesting to complete our calculations of $I\left(a,N\right)$ in the regions $a=O\left(1/N\right)$ and $1-a=O\left(1/N\right)$, where WKB theory breaks down. It would be also interesting to study statistics of the time-averaged population size for stochastic populations which exhibit switches between two metastable states. In the deterministic description, these populations have two attracting fixed points, and a repelling point in between. A host of additional interesting questions concerns multi-population systems which describe coexistence of species.

Finally, a giant disparity and a second-order dynamical phase transition should also appear, in the weak-noise limit, in LDFs of time-averaged quantities in continuous systems with multiplicative noise. In analogy to our work, this can happen in the presence of an attracting fixed point $q\neq 0$ and a repelling fixed point $q=0$ in the deterministic limit.

ACKNOWLEDGMENTS

We are grateful to Hugo Touchette and Tal Agranov for useful discussions, and acknowledge support from the Israel Science Foundation (Grant No. 807/16). N.R.S. was supported by the Clore Foundation.