WKB-type-of approximation for rare event statistics in reacting systems

TL;DR

This paper develops a WKB-based approximation method to calculate rare event probabilities in reacting particle systems by transforming the master equation into a Schrödinger-like form and applying saddle point techniques.

Contribution

It introduces a WKB-type approximation for the master equation in reacting systems, extending previous methods to analyze systems with analytically intractable dynamics.

Findings

Derived a WKB approximation for rare event probabilities.

Applied the method to various example systems.

Extended previous analytical approaches to more complex systems.

Abstract

We calculate the probabilities to find systems of reacting particles in states which largely deviate from typical behavior. The rare event statistics is obtained from the master equation which describes the dynamics of the probability distribution of the particle number. We transform the master equation by means of a generating function into a time-dependent "Schr\"odinger equation". Its solution is provided by a separation ansatz and an approximation for the stationary part which is of Wentzel-Kramers-Brillouin (WKB) type employing a small parameter. The solutions of the "classical" equations of motions and a saddle point approximation yield the proper generating function. Our approach extends a method put forward in [V. Elgart and A. Kamenev, Phys. Rev. E 70, 041106 (2004)]. We calculate the rare event statistics for systems where the dynamics cannot be entirely analyzed in an analytical manner. We consider different examples.

Equations

$$\frac{d}{dt}{P}_n(t)=\sum _{i=1}^m\left({\lambda }_i{h}_i(n-{\nu }_i)P_{n-{\nu }_i}(t)-{\lambda }_i{h}_i(n)P_n(t)\right)$$

$$G(\xi ,t)=\sum _{n=0}^{\infty }{\xi }^n{P}_n(t).$$

$$⟨n⟩=\sum _{n=0}^{\infty }n{P}_n(t)={\frac{\partial }{\partial \xi }G(\xi ,t)|}_{\xi =1}$$

$$P_n(t)=\frac{1}{n!}{\frac{{\partial }^n}{\partial {\xi }^n}G(\xi ,t)|}_{\xi =0}.$$

$$P_n(t)=\frac{1}{2\pi i}\oint d\xi G(\xi ,t){\xi }^{-n-1}$$

$$\frac{\partial }{\partial t}G(\xi ,t)=\hat{\mathcal{L}}G(\xi ,t),$$

$$i\lambda \frac{\partial }{\partial it}G(\xi ,t)=\hat{H}G(\xi ,t),$$

$$\hat{\pi }=-i\lambda \frac{\partial }{\partial \xi },$$

$$i\lambda \frac{\partial }{\partial it}\psi (it)$$

$$\hat{H}\phi (\xi )$$

$$\psi (it)=\psi (0)e^{Et/\lambda }.$$

$$\phi (\xi )=A(\xi )\exp (iS(\xi )/\lambda )$$

$$H=H(\xi ,\pi )=E$$

$$\frac{d\xi }{dit}$$

$$\frac{d\pi }{dit}$$

$$\frac{d\xi }{dit}=\tau (\xi ,E),$$

$$\hat{\pi}G(\xi,t)\Big{|}_{\xi=1}=\frac{\lambda}{i}\left.\frac{\partial G(\xi,t)}{\partial\xi}\right|_{\xi=1}=\frac{\lambda}{i}\langle n\rangle=\bar{\pi}(t)\;.$$

$${\pi }_0=\frac{\lambda }{i}{n}_0.$$

$${\int }_{{\xi }_0}^{\xi }\frac{d{\xi }^{\prime }}{\tau ({\xi }^{\prime },E)}=it.$$

$$P_n(t)=\frac{1}{2\pi i}\oint d\xi g(\xi )\exp (f(\xi ,E))$$

$$f^{\prime }(\xi ,E)=\omega (\xi ,E,E^{\prime }).$$

$$0=\frac{\partial }{\partial \xi }{\int }_{{\xi }_0}^{\xi }\frac{d{\xi }^{\prime }}{\tau ({\xi }^{\prime },E(\xi ,t))}.$$

$$E^{\prime }$$

$$\omega ({\xi }_s,{\xi }_0)=0.$$

$${\int }_{{\xi }_0({\xi }_s)}^{{\xi }_s}\frac{d{\xi }^{\prime }}{\tau ({\xi }^{\prime },E_0({\xi }_0({\xi }_s)))}=it.$$

$$f^{\prime \prime }(\xi ,E)=$$

$$E^{\prime \prime }=\frac{{\partial }^{2}E(\xi )}{\partial {\xi }^2}=-(E^{\prime }{)}^2$$

$$P_n(t)$$

$$\frac{d}{dt}{P}_n(t)=\frac{\lambda }{2}\left((n+2)(n+1)P_{n+2}(t)-n(n-1)P_n(t)\right),$$

$$i\lambda \frac{\partial }{\partial it}G(\xi ,t)=\frac{{\lambda }^2}{2}(1-{\xi }^2)\frac{{\partial }^2}{\partial {\xi }^2}G(\xi ,t),$$

$$\hat{H}=\frac{{\lambda }^2}{2}(1-{\xi }^2)\frac{{\partial }^2}{\partial {\xi }^2}=-\frac{1}{2}(1-{\xi }^2){\hat{\pi }}^2.$$

Full text

I Introduction

Rare events are important in many situations, especially when their consequences are extreme Mühlbacher and Guhr (2018a). Often, demographic noise or intrinsic fluctuations which are based on stochasticity cause such events. There is a wide variety of applications where stochasticity plays a key role, such as chemical reactions Cardelli et al. (2016); Van Kampen (2007), population dynamics Bartlett (1960); Allen (2010); Nisbet and Gurney (1982); Nåsell (2011); Altland et al. (2011), epidemiology Andersson and Britton (2012); Allen et al. (2008); Keeling and Ross (2008); Black and McKane (2010); Allen (2015) and financial markets Mühlbacher and Guhr (2018b); Gardiner (2009) to name a few examples. Another example for stochasticity is non-demographic, i.e. extrinsic, noise which arises due to interactions between the considered system and a noisy environment Mao et al. (2002); Roberts et al. (2015); Lande (1993); Kamenev et al. (2008). If one ignores the noise, many of such systems can be described by rate equations to obtain a macroscopic, deterministic description of average quantities of the system, particularly mean concentrations. This is referred to as the mean field approximation. However, we are interested in large deviations from a typical system behavior. For example, in population dynamics, this is the case when the average size of the population is large and we are interested in low population numbers. The state of this system is given by the population number at a distinct time. The probabilistic description of a stochastic system in such a state is given by the master equation Gillespie (1992); Weber and Frey (2017). Despite its simplicity it is not possible, apart from a few exceptions Jahnke and Huisinga (2007); McQuarrie (1963); McQuarrie et al. (1964); Laurenzi (2000), to solve the master equation analytically Walczak et al. (2012); Jenkinson and Goutsias (2012); Schnoerr et al. (2017). Hence, other methods are called for. There are many kinds of numerical simulations Li et al. (2008); Pahle (2009); Gillespie et al. (2013); Mauch and Stalzer (2011) to model the stochastic kinetics. The stochastic simulation algorithm Gillespie (1976, 1977, 2007) simulates sample paths of the underlying stochastic process, especially for chemical reactions. Due to limited computational resources simulations are not appropriate for large systems and not for the determination of rare events either. An approximation of the master equation is obtained by a Fokker-Planck equation Risken (1996). To this end one has to apply a van Kampen system size expansion. The Fokker-Planck equation is reliable only for small deviations from the mean field approximation. It fails to give an accurate description of large deviations from the typical evolution of a system Doering et al. (2005); Gaveau et al. (1996). Elgart and Kamenev Elgart and Kamenev (2004) introduced an asymptotic method to calculate the rare event statistics in reaction-diffusion systems by formally relating it to semiclassics. A quantum problem can sometimes be solved by expanding to lowest order in $\hbar$ around the corresponding classical equations of motion. A “Hamiltonian” formulation of reaction-diffusion systems is developed which reformulates the master equation by means of a generating function as a “Schrödinger equation”. The “semiclassical” dynamics of the corresponding Hamiltonian provides all the information necessary for the analysis of rare event statistics. The existence of a small parameter allows the treatment analogous to a Wenzel-Kramers-Brillouin (WKB) approximation Landau and Lifshitz (1977); Miller and Good (1953); Messiah (1967); Griffiths and Schroeter (2018); Bender and Orszag (2013). The WKB approximation is applied to the evolution equation of the generating function. For related studies on the WKB approximation see Assaf and Meerson (2006); Assaf et al. (2010); Assaf and Meerson (2017a); Be’er et al. (2015); Be’er and Assaf (2016); Dykman et al. (1994); Gang (1987); Kubo et al. (1973); Peters et al. (1989) and for reviews see Assaf and Meerson (2017b); Ovaskainen and Meerson (2010); Bressloff (2017). For stochastic population models the WKB approximation allows one also to calculate the mean extinction times and probabilities Assaf and Meerson (2007, 2010); Kessler and Shnerb (2007); Assaf et al. (2009) and for switching rates in multistep reactions see Escudero and Kamenev (2009).

Here, we consider single-species chemical reactions which can be described by master equations that give the time evolution of the probability to find a distinct number of particles at a given time. One way of dealing with such systems is a spectral formulation and a stationary WKB approximation, see Assaf and Meerson (2006). A different method is the time-dependent semiclassical approximation Elgart and Kamenev (2004). Both methods have different regimes of validity and accuracy. The time-dependent semiclassical approximation is accurate for $1\ll n\leq\langle n\rangle$ which also means for not too long times. In the region $n>\langle n\rangle$ the accuracy of this method breaks down and the spectral formulation and stationary WKB approximation is better suited. We generalize the time-dependent model Elgart and Kamenev (2004) in such a way that reactions with more than two particles of one species can be analyzed. This is not possible in the original method since the corresponding equations cannot be solved analytically. Moreover, we include the calculation of the pre-exponential factor of the distribution. This was disregarded in the original method. Our main focus lies on the probability of rare events, i.e., to find our system in a state far away from the typical behavior. By means of the probabilities we are able to estimate many quantities such as the average extinction time and the lifetime distribution. Even though large deviations from a typical system behavior may be hardly observable, their probabilities are interesting for anyone who has to compensate probable risks which come along with these rare events. In more detail we consider systems which can have, depending on the initial condition, an absorbing state with zero particles left. In general, results for such systems are unavailable in analytical form. However, with our approach, we are able to compute the solution partially in analytical and partially in numerical form. Depending on the reaction scheme, we have to scale the corresponding parameters of the master equation in order to perform the WKB approximation in a proper way. We consider different types of reactions where we combine death or single annihilation, binary annihilation and triple annihilation. As a result we find a very good accordance between the WKB approximation and the exact solution of the master equation especially in the tail of the distribution.

The paper is structured as follows. In section II we go into the details of the approach and generalize it in such a way that more complex reactions like higher order annihilations can be analyzed. In section III we apply the method on a set of examples which can either be solved exactly or by means of the WKB approximation and numerical calculus. We conclude in section IV.

II Deeper look at the method

Consider a system consisting of identical particles which can react with each other according to different reaction schemes $i=1,\dotsc,m$. The chemical reactions can be described by master equations which give the time evolution of the probability to find a given number of particles in the system at a given time. A particular reaction occurs with probability $\lambda_{i}\Delta t$ with $\lambda_{i}\ll 1$ in the time interval $[t,t+\Delta t]$, where $\lambda_{i}$ is the specific probability rate constant. Importantly, all $\lambda_{i}$ are independent of $\Delta t$ or the considered time interval. Since each particle may react with each other the system is fully described by the following master equation

[TABLE]

where $P_{n}(t)$ denotes the probability to find $n$ particles at time $t$, $h_{i}(n)$ is the number of combinations of reacting particles in the system under reaction scheme $i$ when $n$ particles are in the system and $\nu_{i}$ is the change of particle number when the reaction $i$ occurs, see Gillespie (1992). In order to obtain an unique solution of the master equation we have to specify an initial condition. This can be any kind of normalized distribution like the Poisson distribution or for our sake we use a fixed particle number $n_{0}$ and hence $P_{n}(0)=\delta_{n,n_{0}}$, with the Kronecker delta $\delta_{n,n_{0}}$.

We introduce the auxiliary variable $\xi$ which formally plays the rôle as a “position” and the generating function

[TABLE]

The generating function has to fulfill the initial condition $P_{n}(0)$ which gives $G(\xi,0)=\xi^{n}$. Furthermore, from the conservation of probability, we find $G(1,t)=1$. By means of the generating function we are able to calculate the average particle number

[TABLE]

and the probability

[TABLE]

If the particle number $n$ is large it is convenient to use Cauchy’s integral formula

[TABLE]

where the integration is performed over a closed contour encircling $\xi=0$. Multiplying both sides of (1) with $\xi^{n}$ and summing over all $n$ yields the partial differential equation

[TABLE]

where $\hat{\mathcal{L}}$ is a linear differential operator that includes powers of $\partial/\partial\xi$. The requirement of analyticity of $G(\xi,t)$ yields “self-generated” boundary conditions. Besides the universal boundary condition $G(1,t)=1$ the others are specific to the problem at hand. However, there is always the physical initial condition which is determined by the value of $P_{n}(0)$ and hence it is $G(\xi,0)$ that needs to be considered. The partial differential equation (6) can formally be written as a time-dependent “Schrödinger equation” with imaginary time

[TABLE]

where the right hand side is the “Hamilton operator” $\hat{H}$ and we define $\lambda=\lambda_{k}$ for a fixed $k$. Now, we employ the formal analogy of the probability $\lambda$ with Planck’s constant $\hbar$. Moreover, in analogy to quantum mechanics we define the momentum operator

[TABLE]

which will acquire the meaning of a counting operator. By inserting the ansatz $G(\xi,t)=\varphi(\xi)\psi(it)$ into (7) we can separate the variables $\xi$ and $t$ and find the two equations

[TABLE]

We will interpret the constant $E$ as energy or Hamilton function of our system. Up to this point we consider a probabilistic description of our problem, i.e., both variables $\xi$ and $t$ can be chosen independently from each other. However, the separation ansatz leads to the formal problem $G(1,t)=\varphi(1)\psi(it)=1$ where $t$ cannot be chosen arbitrarily. We will solve this problem by a semiclassical approximation in which we use the classical equations of motion that give an explicit dependence of both variables $\xi$ and $t$. The solution of Eq. (9) is

[TABLE]

We solve the stationary “Schrödinger equation” (10) with a WKB approximation. We insert the ansatz

[TABLE]

and separate the resulting equation into its real and imaginary parts. To solve the pair of differential equations we use the standard WKB assumption that all terms of second order or higher in the small parameter, here $\lambda$, can be neglected. The WKB approximation requires that the “quantum” fluctuations are weak, which is true as long as $\langle n(t)\rangle\gg 1$, i.e., for times not too long. In this regime we can apply the condition $\lambda\ll 1$. By putting all solutions of the differential equations together we find the solution for the generating function $G(\xi,t)$. This solution depends on a constant which has to be determined by the initial condition $G(\xi,0)$. Moreover, the generating function also depends on the energy $E$ of our system. To make progress we have to determine the value of $E$.

Up to this point we used a probabilistic interpretation of the system where the variables $\xi$ and $t$ can be chosen independently from each other. To approximate the function $G(\xi,t)$ and determine the value of $E$ we now go into semiclassics. The energy is given by the Hamilton function $H$ which can be inferred from the “Hamilton operator” $\hat{H}$ and reads

[TABLE]

where $\pi$ is the classical “momentum”. The energy is an integral of motion with $dE/dt=0$. The classical equations of motion in imaginary time are the “Hamilton equations”

[TABLE]

Due to the classical equations of motion (14) and (15), the variables $\xi$ and $t$ are no longer independent of each other. From this point, $\xi$ and $\pi$ are viewed as functions of $t$. We solve the energy $E$ for $\pi$ and insert the result into (14). Now, the first equation of motion contains only $\xi$ and the constant energy $E$ and we find

[TABLE]

which gives the explicit classical dependence $\xi=\xi(t)$. Importantly, for the determination of the initial condition of the generating function we find $\xi(0)=\xi_{0}$. In the original method Elgart and Kamenev (2004) the assumption $\xi_{0}=1$ was made. In section III.1 we show in a direct comparison between our and the original method the benefit of taking the variable $\xi_{0}$ into account, rather than making the approximation $\xi_{0}=1$. We do not solve the second equation of motion (15), instead we consider the mean field dynamics. If we are interested in the average particle number $\langle n\rangle$, we have to know the generating function in the vicinity of $\xi=1$, see Eq. (3). This is the constant mean field solution $\bar{\xi}=1$ which solves Eq. (14) because every permissible Hamilton function must satisfy the condition $H(1,\pi)=0$ due to the normalization of probability. By application of the mean field solution the second equation of motion (15) has the solution $\bar{\pi}(t)$ where we specify the initial condition $\bar{\pi}(0)=\pi_{0}$. Acting with the previously defined momentum operator $\hat{\pi}$ on the generating function at the mean field solution gives

[TABLE]

The last equality holds because in the mean field approximation quantum mechanics has to coincide with classical mechanics. Thus, we call the operator $\hat{\pi}$ counting operator which gives the average particle number for a given time. As $\pi_{0}$ has to be constant for all times we can specify its value at time $t=0$ where, due to the initial condition, we find $\langle n\rangle_{t=0}=n_{0}$ and therefore

[TABLE]

We note that $E$ is constant and can also be expressed in terms of $\xi(0)=\xi_{0}$ and $\pi(0)=\pi_{0}=-i\lambda n_{0}$, i.e., $E=E(\xi_{0})=E_{0}$. Eq. (16) can be solved by integration. If we allow for complex reactions the function $\tau(\xi,E)$ will become also complex. In many cases it might not be possible to solve the differential equation (16) analytically. Instead a numerical solution has to be taken into account. Examples therefore are the reaction which combines binary and single annihilation and the combined reaction up to third order annihilation which are discussed in section III.3 and III.4, respectively. This issue is the starting point of our more general approach to calculate the event statistics.

The formal solution of Eq. (16) is

[TABLE]

The initial condition $G(\xi,0)=\xi^{n_{0}}_{0}$ enables us to determine the constant in the generating function $G(\xi,t)$. Before we determine the value $\xi_{0}$ we first calculate the probability $P_{n}(t)$ by means of a saddle point approximation. At this point we have to make an important remark. As already discussed, the starting point of our analysis is a probabilistic description of the process. We proceed by using semiclassical methods in order to find an approximation for the generating function. Therefore, we derive equations of motion that describe the classical trajectories for a constant energy $E$. The classical trajectories give an explicit dependence of the underlying variables which in our case are $\xi$, $\pi$ and $t$. They cannot be chosen independently from each other. The value of the energy determines which trajectory in the phase space describes the relationship between $\xi$ and $t$, i.e., $\xi=\xi(t)$. Eventually, we are interested in probabilities $P_{n}(t)$ again. Therefore, we have to leave the semiclassical description behind and turn back to a probabilistic description. In a probabilistic description the variables $\xi$ and $t$ are independent of each other and can also be chosen independently. Classically this means that we do not move along the trajectories anymore. This only applies if the energy depends on both variables, i.e., $E=E(\xi,t)$. In other words, when we go back from the semiclassical into the probabilistic description the energy does not remain constant. In that sense it is more intuitive to name the quantity $E$ Hamilton function instead of energy. This non-constant Hamilton function can be explained by a short example. In a conservative system the energy is fully determined by the initial conditions. It however differs for different initial conditions. In that sense, the Hamilton function describes, at the same time, a constant of motion and the same system at different energies. Hence, derivatives $E^{\prime}=\partial E/\partial\xi$ have to be taken into account when the saddle point approximation is performed.

The contour integral can be written as

[TABLE]

with the function $g(\xi)$ and the “free energy” $f(\xi,E)$ that are determined by the generating function $G(\xi,t)$ and the factor $\xi^{-n-1}$ in Cauchy’s integral formula (5). We add the term $\xi^{-n}$ into the exponential because of its factor $n$ which is supposed to be large. We calculate the integral by means of a saddle point approximation. This approximation is justified because we have a small $\lambda$ and assume large $n$. Furthermore, we are interested in times $t$ where $\langle n(t)\rangle$ is sufficiently smaller than $n_{0}$ but $\langle n(t)\rangle\gg 1$ still holds.

The saddle point approximation requires $f^{\prime}(\xi_{s},E)=0$ and we find with a non-constant energy $E(\xi,t)$

[TABLE]

For convenience we drop the arguments of $E(\xi,t)$. Pretending that the analytical solution of (16) is not known, i.e., when we have to evaluate the integral in (19) numerically, we have to determine $E^{\prime}$ in order to solve the saddle point equation $f^{\prime}(\xi_{s},E)=0$. We obtain the derivative of $E$ by deriving the solution of the classical equation of motion (19) with respect to $\xi$ under the condition that the energy depends on $\xi$ and the time $t$ does not

[TABLE]

We remark that the integration is performed over the variable $\xi^{\prime}$ which is due to semiclassics not included in $E(\xi,t)$. By using the Leibniz rule we find

[TABLE]

which again is a function of $\xi$ and $E$. The energy $E$ can be expressed in dependence of $\xi_{0}$, i.e., $E=E(\xi_{0})$. Due to this relation $E^{\prime}$ can also be expressed in terms of $\xi$ and $\xi_{0}$, i.e., $E^{\prime}=\kappa(\xi,\xi_{0})$. Now, we can reduce the derivative $f^{\prime}(\xi,E)=\omega(\xi,\xi_{0})$ on the two variables $\xi$ and $\xi_{0}$ and the saddle point condition becomes

[TABLE]

From its solution we obtain the relationship between $\xi_{s}$ and $\xi_{0}$, i.e., $\xi_{0}=\xi_{0}(\xi_{s})$. Together with Eq. (19) this allows us to calculate a numeric value for $\xi_{s}$, by solving the integral equation

[TABLE]

We note that the integral stems from the semiclassical part of our analysis, hence the energy $E$ has to remain constant during integration. Furthermore we changed the upper limit in (19) to $\xi_{s}$ which is valid, as the constant value of the energy depending on $\xi_{s}$ has to be chosen accordingly. Finding a solution for $\xi_{s}$ can be demanding, depending on the function $\tau(\xi,E)$ and the solution $\xi_{0}(\xi_{s})$. Once, a solution for $\xi_{s}$ is found the value of $\xi_{0}$ can be calculated. The last piece missing for the saddle point approximation is the second derivative of the free energy

[TABLE]

with the second derivative of the energy

[TABLE]

The value of $E$ is determined by $E(\xi_{0})=E_{0}$. Combining all pieces yields the final result for the probability

[TABLE]

if $f^{\prime\prime}(\xi,E)$ is always real. This formalism works for every kind of reaction scheme. However, if the reactions become too complex, e.g. by involving higher order annihilations, it might become impossible to solve the required equations due to technical complications.

III Some examples

III.1 Binary annihilation revisited

For the convenience of the reader we show how the known results of binary annihilation, studied in Elgart and Kamenev (2004), fit into our generalized approach. This model can be solved analytically and exactly McQuarrie et al. (1964).

We consider a system where only the binary annihilation can possibly occur with probability rate $\lambda$. Once, this reaction takes place two particles form an inert aggregate. Since each particle may react with each other the system is fully described by the following master equation

[TABLE]

where $P_{n}(t)$ denotes the probability to find $n$ particles at time $t$. According to Eq. (7) the corresponding time-dependent “Schrödinger equation” with imaginary time is

[TABLE]

where the right hand side can also be written in terms of the “Hamilton operator”

[TABLE]

The separation ansatz $G(\xi,t)=\varphi(\xi)\psi(it)$ yields equation (11) for the time-dependent part and

[TABLE]

for the stationary part. We solve the stationary “Schrödinger equation” with a WKB approximation. We insert the ansatz (12) and separate the resulting equation into its real and imaginary parts which yields

[TABLE]

Eq. (34) can be simplified

[TABLE]

which yields

[TABLE]

with a constant $\tilde{c}$. In Eq. (33) we use the standard WKB assumption that all terms of second order in the small parameter, here $\lambda$, can be neglected. This gives the simplified differential equation

[TABLE]

Its solution is

[TABLE]

with some initial value $\xi_{0}$. Now, we can put all terms together and obtain the generating function

[TABLE]

where all constants have been put into $c$. Before determining this constant by the initial condition of the generating function, we first go into semiclassics. The Hamilton function is given by

[TABLE]

The classical equations of motion in imaginary time are

[TABLE]

We solve Eq. (40) for $\pi$ and insert the solution into Eq. (41). The ensuing differential equation

[TABLE]

has the solution

[TABLE]

which is analogous to Eq. (19). The mean field dynamics gives us the average particle number. By application of the mean field solution $\bar{\xi}=1$ the second equation of motion (42) has the solution

[TABLE]

The average particle number is

[TABLE]

for $n_{0}\gg 1$. The energy can be written as

[TABLE]

Now, we are able to determine the constant $c$ in Eq. (39) with the initial condition $G(\xi,0)=\xi^{n}_{0}$

[TABLE]

Obviously, the normalization $G(1,t)=1$ remains preserved. We calculate the probability $P_{n}(t)$ by means of the contour integral (5) and a saddle point approximation. The contour integral now reads

[TABLE]

with

[TABLE]

The saddle point approximation requires $f^{\prime}(\xi_{s},E)=0$ and we find with a non-constant Hamilton function $E(\xi,t)$

[TABLE]

For convenience we drop the arguments of $E(\xi,t)$. We see that inserting the solution (44) of the classical equation of motion deletes the first bracket and yields

[TABLE]

Interestingly, the $E^{\prime}$ term has been removed and with Eq. (47) we immediately obtain the relation $\xi_{0}=\xi_{0}(\xi_{s})$ according to Eq. (24)

[TABLE]

Combining Eqs. (54), (47) and (44) yields an implicit equation for $\xi_{s}$ which has to be solved by numerical methods. Having found the solution for $\xi_{s}$, a value for $\xi_{0}$ can be obtained by Eq. (54). In order to obtain the second derivative of the free energy we need $E^{\prime}$ which is according to Eq. (23)

[TABLE]

This result can directly be obtained by deriving Eq. (44). The second derivative of $f(\xi)$ at $\xi_{s}$ can be simplified to

[TABLE]

which is always real for real $\xi_{s}$. Combining all pieces yields the final result for the probability

[TABLE]

In Fig. 1 we compare the approximation (57) with the exact solution. We set $n_{0}=300$, $\lambda=0.01$ and $t=0.5$. The exact solution is calculated by numerical integration of the master equation. The maximum probability is around the average particle number which according to Eq. (46) is $\langle n(0.5)\rangle=120$. Overall the approximation coincides over many orders of magnitude very well with the exact solution.

We compare the new result (57) with the original one. In fact, we do not use the original result of Elgart and Kamenev (2004), where the prefactor was included manually, we use the time-dependent solution of Assaf and Meerson (2006) that includes the prefactor by calculation. Fig. 2 shows on a log-scale the ratios of the exact result and the semiclassical approximation for both, the new result (57) and the original result. Clearly, in the region $n\leq\langle n\rangle$ we see the high accuracy of the new result which outperforms the original result. For large $n$ the accuracy of the time-dependent semiclassical approximation deteriorates.

III.2 Single annihilation

The simplest reaction is the single annihilation, where one particle forms an inert aggregate with probability rate $\lambda$. In terms of population dynamics this reaction is called death process. This process can be solved analytically and it is well known in literature, see, e.g. Chang (1968); Jahnke and Huisinga (2007). Nevertheless, we apply the method to give a different representation how the model can be solved. Its master equation reads

[TABLE]

This master equation can be solved exactly by different methods, see Jahnke and Huisinga (2007). By means of the generating function (2) Eq. (58) can be written as “Schrödinger equation”

[TABLE]

with the momentum operator $\hat{\pi}$. By applying the separation ansatz, we find the same formal solution for the time-dependent part as for the binary annihilation, see Eq. (11). The $\xi$-dependent part of the separation ansatz is more interesting

[TABLE]

it has the exact solution

[TABLE]

The classical energy is $E=i\lambda(1-\xi)\pi$ and the classical equations of motion are

[TABLE]

Both equations are easily solvable. The second equation gives the mean field dynamics resulting in the average particle number which is $\langle n\rangle=n_{0}e^{-\lambda t}$. The solution for (62)

[TABLE]

can be inserted into (61) which yields $\varphi(\xi)=\varphi(\xi_{0})e^{-Et/\lambda}$. Hence, we find the generating function

[TABLE]

In the last step we used the initial condition that at time $t=0$ there are $n_{0}$ particles in the system. Solving Eq. (64) for $\xi_{0}$ and inserting into the generating function yields

[TABLE]

The probability to find $n$ particles at time $t$ can be calculated by means of Eq. (4). We find

[TABLE]

III.3 Binary and single annihilation

We consider a system where two reactions can occur. The binary annihilation occurs with probability rate $\lambda$ and the single annihilation occurs with probability rate $\sigma$. Both probabilities are of the same order of magnitude. The master equation reads

[TABLE]

Its corresponding “Schrödinger equation” for the generating function is

[TABLE]

where we use the momentum operator $\hat{\pi}=-i\lambda\partial/\partial\xi$ in the second line. The separation ansatz yields to the known result (11) for the time-dependent part and for the stationary part we find

[TABLE]

In the manner of section III.1 we use the ansatz $\varphi(\xi)=A(\xi)\exp(iS(\xi)/\lambda)$, sort the resulting equation for its real and imaginary part. Neglecting all terms with prefactor $\lambda^{2}$ or $\lambda\sigma$ we find

[TABLE]

The solution of (71) is (38) which is already discussed. We can insert the derivatives of $S(\xi)$ and insert them into (72) to obtain a differential equation for $A(\xi)$ with the solution

[TABLE]

with a constant $c$. Hence, the generating function, fulfilling the initial condition $G(\xi_{0},0)=\xi^{n_{0}}_{0}$, is

[TABLE]

The classical energy is $E=\tfrac{1}{2}(\xi^{2}-1)\pi^{2}-i\sigma(\xi-1)\pi$ and the classical equations of motion are

[TABLE]

First, we evaluate the average particle number by means of the mean field dynamics

[TABLE]

The solution gives the average particle number

[TABLE]

In order to make analytical progress we simplify Eq. (75) furthermore by dropping the $\sigma^{2}$ term. This is fully consistent with the WKB method since by Taylor approximation the leading order of Eq. (75) is $\sigma^{2}$. Due to this approximation we find the same classical dynamics as in the case of binary annihilation, see equations (41) and (42). The solution for the corresponding classical equation of motion, determining $\xi$, is given in (44). Hence, the generating function has the same arguments in the exponential function as the generating function of the binary annihilation (39). It solely differs from (39) by the additional prefactor $\left(\tfrac{\xi_{0}+1}{\xi+1}\right)^{\sigma/\lambda}$. We calculate the probability $P_{n}(t)$ by means of the contour integral with a saddle point approximation. The saddle point condition (53) combined with the solution of the approximated classical equation of motion (44) gives the dependence

[TABLE]

Recombining this equation with (44) yields an equation for $\xi_{s}$ which can be solved numerically. Finally the probability reads

[TABLE]

We notice that the probability has the same form as the probability for the binary annihilation. The differences are the additional prefactor and the value for $\xi_{s}$ at the saddle point. The additional prefactor contributes only little to the value of (80) compared with the exponential. Furthermore, the value of the saddle point does not vary significantly from that of the binary annihilation. Hence, the probability $P_{n}(t)$ for the single and binary annihilation does not vary considerable from that for the pure binary annihilation. This means that the binary annihilation dominates the single annihilation considerably. In Fig. 3 we compare the approximation with the exact probabilities.

We find some structure of the approximated probability near the maximum of the distribution. This is due to the prefactor which approaches zero in the denominator. Usually, in the saddle point approximation we consider the exponent only and do not consider the prefactor. If we do so, we would have to adapt an artificial prefactor in such a way that the normalization condition is still fulfilled. In that case the structure would vanish. Overall the approximation coincides over many orders of magnitude well with the exact solution.

III.4 Third order reaction

We study a reaction which combines the single and binary annihilation from the last example with the trimolecular annihilation. The master equation reads

[TABLE]

We introduce the rate constant $\mu$ for the trimolecular reaction, which we assume to be of order $\lambda^{3}$. The rate constant $\sigma$ of the single annihilation is of order $\lambda$. The corresponding “Schrödinger equation” is

[TABLE]

with the momentum operator $\hat{\pi}=-i\lambda\partial/\partial\xi$. Applying the separation ansatz $G(\xi,t)=\psi(t)\varphi(\xi)$ yields to the result (11) for the time dependent part. The stationary “Schrödinger equation” can be written as

[TABLE]

We insert the ansatz $\varphi(\xi)=A(\xi)\exp(iS(\xi)/\lambda)$ to solve this equation. We sort for the real and imaginary parts and drop all terms of order $\lambda^{2}$ and above which gives

[TABLE]

For convenience of the reader, we dropped the arguments of $A(\xi)$ and $S(\xi)$. Once again we remark the scaling behavior $\sigma\sim\lambda$ and $\mu\sim\lambda^{3}$. Interestingly, we find Eq. (84) which we also found in the last example and also for the binary annihilation. Its solution is given in Eq. (38). Furthermore, each single term of Eq. (85) is of order $\lambda$. This means that we formally can cancel the small parameter $\lambda$ and so its order of magnitude does not appear in both equations. With $S^{\prime}(\xi)$ from (84) we can solve Eq. (85) and find

[TABLE]

With the initial condition $P_{n}(0)=\delta_{n,n_{0}}$ we are able to determine the constants. The probability now reads

[TABLE]

with

[TABLE]

To make progress we specify the energy

[TABLE]

in terms of $\xi$. Therefore we have to solve the equations of motion

[TABLE]

by inserting the solution of the energy for $\pi$ into Eq. (92) which gives $\partial\xi/\partial it=\tau(\xi,E)$. The resulting equation cannot be solved analytically. Hence, we have to solve this equation by evaluating the integral (19) numerically. However, this can only be done if we know the functional coherence of $\xi$ and $\xi_{0}$, which is obtained by the saddle point approximation, i.e., $f^{\prime}(\xi_{s},E)=0$. We find

[TABLE]

with $E^{\prime}$ from Eq. (23). Again, this equation has to be solved numerically. A value for $\xi_{s}$ can be obtained by solving Eq. (25). This procedure requires the numerical solution of two equations combined with a numerical integration. We remark that there are three possible solutions for $\pi$ when solving Eq. (90). This is relevant for the numerics as we find $\xi_{s}<1$ for $n<\langle n\rangle$ and $\xi_{s}>1$ for $n>\langle n\rangle$ where $\langle n\rangle$ is the average particle number. By performing the numerics we have to use different solutions for $\pi$ when we calculate the values for $\xi_{s}$. This is depending on the value of $n$ being above or below $\langle n\rangle$. The mean field solution $\bar{\xi}=1$ gives us the dynamics of the average particle number

[TABLE]

with $\bar{\pi}=-i\lambda\langle n\rangle$. We show the average particle number $\langle n(t)\rangle$ in Fig. 4. We specify the parameters $\lambda=\tfrac{1}{100}$, $\sigma=\lambda$, $\mu=\lambda^{3}$ and $n_{0}=300$. The average particle number at $t=1$ is $\langle n(1)\rangle\approx 74.15$. Hence, the transition $\xi_{s}>1$ arises when $n\geq 75$ for $t=1$ if we only consider integer particle numbers.

Having found the values for $\xi_{s}$ we can now perform the saddle point approximation which gives us

[TABLE]

with the energy

[TABLE]

and the second derivative

[TABLE]

which is always real. In Fig. 5 we show the probability distribution $P_{n}(1)$ for the same parameter specification as in Fig. 4. Again we find a good agreement over many orders of magnitude between the semiclassical approximation and the exact solution, especially for the left tail of the distribution.

IV Conclusions

Extending the approach by Elgart and Kamenev, we put forward a model which is able to describe the probability of rare events in reaction-diffusion systems described by master equations. The systems consist of single-species particles with infinite-range interaction and we assume that spatial degrees of freedom are irrelevant. By means of a generating function we transform the master equation into a time-dependent “Schrödinger equation” in imaginary time. As the master equation gives a probabilistic description of the system, so does the evolution equation for the generating function. In short, the master equation is equivalent to the time-dependent “Schrödinger equation” derived by means of the generating function. We separate the time-dependent part of the “Schrödinger equation” and the stationary part and thereby introduce a new constant which we will interpret as the energy or Hamilton function of the system. We apply a WKB approximation to solve the stationary part of the “Schrödinger equation”, we identify a small parameter which is the analogue to Planck’s constant $\hbar$. In order to find the rare event statistics we are interested in large deviations from a typical system behavior. A typical behavior can be calculated by the mean field approximation which allows us to calculate quantities such as the average particle number. Within the WKB approximation we can derive the classical equations of motion which determine the phase portrait of the system. The trajectories of the phase portrait are determined by the energy which is constant in the semiclassical description. These classical equations of motion need to be solved either analytically or numerically in order to find the solution for the generating function. Finally, the probability to find $n$ particles at time $t$ is calculated by Cauchy’s integral formula on which we apply a saddle point approximation. At this point we make a transition from the semiclassical description in which the dynamics is determined by the phase portrait to a probabilistic description. Hence, when performing the saddle point approximation the energy does not remain constant and becomes a function of the phase space parameters. This dependence is determined by the classical equations of motion. We remark that the accuracy of the method breaks down in the vicinity of the average particle number and for small times when the condition $\langle n(t)\rangle\ll n_{0}$ is not satisfied.

We studied some systems. The single annihilation is exactly solvable. The binary annihilation process and the process which combines binary and single annihilation is analytically solved by means of the WKB approximation. For both systems we find good agreement over many orders of magnitude with the exact solution of the master equation. Finally, we analyzed a process that combines single, binary and triple annihilation. In general, the leading order of the stationary “Schrödinger equation” is obtained by the order of the highest annihilation process. The third order differential equation can be solved analytically when we require a certain scaling of its parameters and use the WKB approximation. However, the semiclassical equations of motion have to be solved numerically. Once more, in the tails of the distribution, we find a very good agreement with the exact solution over many orders of magnitude.

Rare events define the tail of the distribution. We are not interested in typical fluctuations near the maximum of the distribution where we find oscillations for some examples. Nontheless, the oscillations can easily be removed by dropping the prefactor of the distribution and replacing it with a normalization constant.

Acknowledgements.

One of us (A.M.) acknowledges support from Studienstiftung des deutschen Volkes.