Master equations and the theory of stochastic path integrals

The theory of continuous-time Markov processes is largely built on two equations: the Fokker–Planck [1–4] and the master equation [4, 5]. Both equations assume that the future of a system depends only on its current state, memories of its past having been wiped out by randomizing forces. This Markov assumption is sufficient to derive either of the two equations. Whereas the Fokker–Planck equation describes systems that evolve continuously from one state to another, the master equation models systems that perform jumps in state space.

Path integral representations of the master equation were first derived around 1980 [6–15], shortly after such representations had been derived for the Fokker–Planck equation [16–19]. Both approaches were heavily influenced by quantum theory, introducing such concepts as the Fock space [20] with its 'bras' and 'kets' [21], coherent states [22–24], and 'normal-ordering' [25] into non-equilibrium theory. Some of these concepts are now well established and the original 'bosonic' path integral representation has been complemented with a 'fermionic' counterpart [26–31]. Nevertheless, we feel that the theory of these 'stochastic' path integrals may benefit from a step back and a closer look at the probabilistic structures behind the integrals. Therefore, the objects imported from quantum theory make place for their counterparts from probability theory in this review. For example, the coherent states give way to the Poisson distribution. Moreover, we use the bras and kets as particular basis functionals and functions whose choice depends on the stochastic process at hand (a functional maps functions to numbers). Upon choosing the basis functions as Poisson distributions, one can thereby recover both a classic path integral representation of averaged observables as well as the Poisson representation of Gardiner and Chaturvedi [32, 33]. The framework presented in this review integrates a variety of different approaches to the master equation. Besides the Poisson representation, these approaches include a spectral method for the computation of stationary probability distributions [34], WKB approximations and other 'semi-classical' methods for the computation of rare event probabilities [35–38], and a variational approach that was proposed in the context of stochastic gene expression [39]. All of these approaches can be treated within a unified framework. Knowledge about this common framework makes it possible to systematically search for new ways of solving the master equation.

Before outlining the organization of this review, let us note that by focusing on the above path integral representations of master and Fokker–Planck equations, we neglect several other 'stochastic' or 'statistical' path integrals that have been developed. These include Edwards path integral approach to turbulence [40, 41], a path integral representation of Haken [42], path integral representations of non-Markov processes [43–58] and of polymers [59–62], and a representation of 'hybrid' processes [63–65]. The dynamics of these stochastic hybrid processes are piecewise-deterministic. Moreover, we do not discuss the application of renormalization group techniques, despite their significant importance. Excellent texts exploring these techniques in the context of non-equilibrium critical phenomena [66–68] are provided by the review of Täuber, Howard, and Vollmayr-Lee [69] as well as the book by Täuber [70]. Our main interest lies in a mathematical framework unifying the different approaches from the previous paragraph and in two path integrals that are based on this framework. Both of these path integrals provide exact representations of the conditional probability distribution solving the master equation. We exemplify the use of the path integrals for elementary processes, which we choose for their pedagogic value. Most of these processes do not involve spatial degrees of freedom but the application of the presented methods to processes on spatial lattices or networks is straightforward. A process with diffusion and linear decay serves as an example of how path integrals can be evaluated perturbatively using Feynman diagrams. The particles' linear decay is treated as a perturbation to their free diffusion. The procedure readily extends to more complex processes. Moreover, we show how the two path integrals can be used for the computation of rare event probabilities. Let us emphasize that we only consider Markov processes obeying the Chapman–Kolmogorov equation and associated master equations [4]. It may be interesting to extend the discussed methods to 'generalized' or 'physical' master equations with memory kernels [71–74].

The organization of this review is summarized in figure 1 and is as follows. In the next section 1.3, we introduce the basic concepts of the theory of continuous-time Markov processes. After discussing the roles of the forward and backward Fokker–Planck equations for processes with continuous sample paths, we turn to processes with discontinuous sample paths. The probability of finding such a 'jump process' in a generic state n at time $\tau >{{t}_{0}}$, given that the process has been in state n₀ at time t₀, is represented by the conditional probability distribution $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$. Whereas the forward master equation evolves this probability distribution forward in time, starting out at time $\tau ={{t}_{0}}$, the backward master equation evolves the distribution $p\left(t,n|\tau,{{n}_{0}}\right)$ backward in time, starting out at time $\tau =t$. Both master equations can be derived from the Chapman–Kolmogorov equation (see left side of figure 1). In section 1.4, we discuss two explicit representations of the conditional probability distribution solving the two master equations. Moreover, we comment on various numerical methods for the approximation of this distribution and for the generation of sample paths. Afterwards, section 1.5 provides a brief historical overview of contributions to the development of stochastic path integrals.

Zoom In Zoom Out Reset image size

The main part of this review begins with section 2. We first exemplify how a generalized probability generating function can be used to determine the stationary probability distribution of an elementary chemical reaction. This example introduces the bra-ket notation used in this review. In section 2.1, we formulate conditions under which a general forward master equation can be transformed into a linear partial differential equation (PDE) obeyed by the generating function. This function is defined as the sum of the conditional probability distribution $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$ over a set of basis functions $\left\{|n\rangle \right\}$, the 'kets' (see middle column of figure 1). The explicit choice of the basis functions depends on the process being studied. We discuss different choices of the basis functions in section 2.2, first for a random walk, afterwards for chemical reactions and for processes whose particles locally exclude one another. Several methods have recently been proposed for the analysis of the PDE obeyed by the generating function. These methods include the variational method of Eyink [83] and of Sasai and Wolynes [39], the WKB approximations [84] and spectral methods of Elgart and Kamenev [35] and of Assaf and Meerson [36–38], and the spectral method of Walczak, Mugler, and Wiggins [34]. We comment on these methods in section 2.3.

In section 3.1, we formulate conditions under which a general backward master equation can be transformed into a novel, backward-time PDE obeyed by a 'marginalized distribution'. This object is defined as the sum of the conditional probability distribution $p\left(t,n|\tau,{{n}_{0}}\right)$ over a set of basis functions $\left\{|{{n}_{0}}\rangle \right\}$ (see middle column of figure 1). If the basis function $|{{n}_{0}}\rangle$ is chosen as a probability distribution, the marginalized distribution also constitutes a true probability distribution. Different choices of the basis function are considered in section 3.2. In section 3.3, the use of the marginalized distribution is exemplified in the calculation of mean extinction times. Afterwards, in section 3.4, we derive yet another linear PDE, which is obeyed by a 'probability generating functional'. This functional is defined as the sum of the conditional probability distribution $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$ over a set of basis functionals $\left\{\langle n|\right\}$, the 'bras'. In section 3.5, we show that the way in which the generating functional 'generates' probabilities generalizes the Poisson representation of Gardiner and Chaturvedi [32, 33].

Sections 4 and 5 share the goal of representing the master equations' solution by path integrals. In section 4.1, we first derive a novel backward path integral representation from the PDE obeyed by the marginalized distribution (see right side of figure 1). Its use is exemplified in sections 4.2, 4.3 and 4.5 in which we solve several elementary processes. Although we do not discuss the application of renormalization group techniques, section 4.5.4 includes a discussion of how the backward path integral representation can be evaluated in terms of a perturbation expansion. The summands of the expansion are expressed by Feynman diagrams. Besides, we derive a path integral representation for Markov processes with continuous state spaces in the 'intermezzo' section 4.4. This representation is obtained by performing a Kramers–Moyal expansion of the backward master equation and it comprises a classic path integral representation [16–19] of the (backward) Fokker–Planck equation as a special case. One can rewrite the representation of the Fokker–Planck equation in terms of an Onsager–Machlup function [79] and, for purely diffusive Brownian motion [82], the representation simplifies to the path integral of Wiener [80, 81]. Moreover, we recover a Feynman–Kac like formula [85], which solves the (backward) Fokker–Planck equation in terms of an average over the paths of an Itô stochastic differential equation [86–88] (or of a Langevin equation [89]).

In section 5, we complement the backward path integral representation with a forward path integral representation. Its derivation in section 5.1 starts out from the PDE obeyed by the generalized generating function (see right side of figure 1). The forward path integral representation can, for example, be used to compute the generating function of generic linear processes as we demonstrate in section 5.2. Besides, we briefly outline how a Kramers–Moyal expansion of the (forward) master equation can be employed to derive a path integral representation for processes with continuous state spaces in section 5.3. This path integral can be expressed in terms of an average over the paths of an SDE proceeding backward in time. Its potential use remains to be explored.

Before proceeding, let us briefly point out some properties of the forward and backward path integral representations. First, the paths along which these path integrals proceed are described by real variables and all integrations are performed over the real line. Grassmann path integrals [26–31] for systems whose particles locally exclude one another are not considered. It is, however, explained in section 2.2.4 how such systems can be treated without the need for Grassmann variables, based on a method recently proposed by van Wijland [90]. Second, transformations of the path integral variables such as the 'Doi-shift' [91] are implemented on the level of the basis functions and functionals. Third, our derivations of the forward and backward path integral representations do not involve coherent states or combinatoric relations for the commutation of exponentiated operators. Last, the path integrals allow for time-dependent rate coefficients of the stochastic processes.

In section 6, we derive a path integral representation of averaged observables (see right side of figure 1). This representation can be derived both from the backward and forward path integral representations (see section 6.1), and by representing the forward master equation in terms of the eigenvectors of creation and annihilation matrices ('coherent states'; see section 6.2). The duality between these two approaches resembles the duality between the wave [92] and matrix [93–95] formulations of quantum mechanics. Let us note that our resulting path integral does not involve a 'second-quantized' or 'normal-ordered' observable [96]. In fact, we show that this object agrees with the average of an observable over a Poisson distribution. In section 6.3, we then explain how the path integral can be evaluated perturbatively using Feynman diagrams. Such an evaluation is demonstrated for the coagulation reaction $2\,A\to A$ in section 6.4, restricting ourselves to the 'tree level' of the diagrams.

In section 7, we review and extend a recent method of Elgart and Kamenev for the computation of rare event probabilities [35]. As explained in section 7.1, this method evaluates a probability distribution by expanding the forward path integral representation from section 5 around 'stationary', or 'extremal', paths. In a first step, one thereby acquires an approximation of the ordinary probability generating function. In a second step, this generating function is transformed back into the underlying probability distribution. The evaluation of this back transformation typically involves an additional saddle-point approximation. In section 7.2, we demonstrate both of the steps for the binary annihilation reaction $2\,A\to \varnothing$, improving an earlier approximation of the process by Elgart and Kamenev [35] by terms of sub-leading order. In section 7.3, we then extend the 'stationary path method' to the backward path integral representation from section 4. The backward path integral provides direct access to a probability distribution without requiring an auxiliary saddle-point approximation. However, the leading order term of its expansion is not normalized. We demonstrate the procedure for the binary annihilation reaction in section 7.4.

Finally, section 8 closes with a summary of the different approaches discussed in this review and outlines open challenges and promising directions for future research.

Our main interest lies in a special class of stochastic processes, namely in the class of continuous-time Markov processes with discontinuous sample paths. These processes are also called 'jump processes'. In the following, we outline the mathematical theory of jump processes and derive the central equations obeyed by them: the forward and the backward master equation. Before going into the mathematical details, let us explain when a system's time evolution can be modelled as a continuous-time Markov process with discontinuous sample paths and what that phrase actually means.

First of all, if the evolution of a system is to be modelled as a continuous-time Markov process, it must be possible to describe the system's state by some variable n. In fact, it must be possible to do so at every point τ in time throughout an observation period [t₀,t]. A variable $n\in \mathbb{Z}$ could, for example, represent the position of a molecular motor along a cytoskeletal filament, or a variable $n\in {{\mathbb{R}}_{\geqslant 0}}$ the price of a stock between the opening and closing times of an exchange. The assumption of a continuous time parameter τ is rather natural and conforms to our everyday experience. Still, a discrete time parameter may sometimes be preferred, for example, to denote individual generations of an evolving population [97]. By allowing τ to take on any value between the initial time t₀ and the final time t, we can choose it to be arbitrarily close to one of those times. Below, this possibility will allow us to describe the evolution of the process in terms of a differential equation.

The (unconditional) probability of finding the system in state n at time τ is represented by the 'single-time' probability distribution $p(\tau,n)$. Upon demanding that the system has visited some state n₀ at an earlier time ${{t}_{0}}<\tau$, the probability of observing state n at time τ is instead encoded by the conditional probability distribution $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$. If the conditional probability distribution is known, the single-time distribution can be inferred from any given initial distribution $p\left({{t}_{0}},{{n}_{0}}\right)$ via $p(\tau,n)={\sum}_{{{n}_{0}}}p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)p\left({{t}_{0}},{{n}_{0}}\right)$. A stochastic process is said to be Markovian if a distribution conditioned on multiple points in time actually depends only on the state that was realized most recently. In other words, a conditional distribution $p\left(t,n|{{\tau}_{k}},{{m}_{k}};\cdots ;{{\tau}_{1}},{{m}_{1}};{{t}_{0}},{{n}_{0}}\right)$ must agree with $p\left(t,n|{{\tau}_{k}},{{m}_{k}}\right)$ whenever $t>{{\tau}_{k}}>{{\tau}_{k-1}},\ldots,{{t}_{0}}$ ¹ . Therefore, a Markov process is fully characterized by the single-time distribution $p\left({{t}_{0}},{{n}_{0}}\right)$ and the conditional distribution $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$. The latter function is commonly referred to as the 'transition probability' for Markov processes [98].

The stochastic dynamics of a system can be modelled in terms of a Markov process if the system has no memory. Let us explain this requirement with the example of a Brownian particle suspended in a fluid [82]. Over a very short time scale, the motion of such a particle is ballistic and its velocity highly auto-correlated [99]. But as the particle collides with molecules of the fluid, that memory fades away. A significant move of the particle due to fluctuations in the isotropy of its molecular bombardment then appears to be completely uncorrelated from previous moves (provided that the observer does not look too closely [100]). Thus, on a sufficiently coarse time scale, the motion of the particle looks diffusive and can be modelled as a Markov process. However, the validity of the Markov assumption does not extend beyond the coarse time scale.

The Brownian particle exemplifies only two of the properties that we are looking for: its position is well-defined at every time τ and its movement is effectively memoryless on the coarse time scale. But the paths of the Brownian particle are continuous, meaning that it does not spontaneously vanish and then reappear at another place. If the friction of the fluid surrounding the Brownian particle is high (over-damped motion), the probability of observing the particle at a particular place can be described by the Smoluchowski equation [101]. This equation coincides with the simple diffusion equation when the particle is not subject to an external force. In the general case, the probability of observing the particle at a particular place with a particular velocity obeys the Klein–Kramers equation [102, 103] (the book of Risken [104] provides a pedagogic introduction to these equations). From a mathematical point of view, all of these equations constitute special cases of the (forward) Fokker–Planck equation [1–4, 104]. For a single random variable $x\in \mathbb{R}$, e.g. the position of the Brownian particle, this equation has the generic form

$$\begin{eqnarray}&&{{\partial}_{\tau}}p\left(\tau,x|{{t}_{0}},{{x}_{0}}\right)=-{{\partial}_{x}}\left[{{\alpha}_{\tau}}(x)p\right]+\frac{1}{2}\partial _{x}^{2}\left[{{\beta}_{\tau}}(x)p\right].\end{eqnarray} \tag{ 1 }$$

The initial condition of this equation is given by the Dirac delta distribution (or generalized function) $p\left({{t}_{0}},x|{{t}_{0}},{{x}_{0}}\right)=\delta \left(x-{{x}_{0}}\right)$. Here we used the letter x for the random variable because the letter n would suggest a discrete state space. The function ${{\alpha}_{\tau}}$ is often called a drift coefficient and ${{\beta}_{\tau}}$ a diffusion coefficient (note, however, that in the context of population genetics, ${{\beta}_{\tau}}$ describes the strength of random genetic 'drift' [105, 106]). For reasons addressed below, the diffusion coefficient must be non-negative at every point in time for every value of x (for a multivariate process, ${{\beta}_{\tau}}$ represents a positive-semidefinite matrix). In the mathematical community, the Fokker–Planck equation is better known as the Kolmogorov forward equation [105], honouring Kolmogorov's fundamental contributions to the theory of continuous-time Markov processes [4]. Whereas the above Fokker–Planck equation evolves the conditional probability distribution forward in time, one can also evolve this distribution backward in time, starting out from the final condition $p\left(t,x|t,{{x}_{0}}\right)=\delta \left(x-{{x}_{0}}\right)$. The corresponding equation is called the Kolmogorov backward or backward Fokker–Planck equation. It has the generic form

$$\begin{eqnarray}&&{{\partial}_{-\tau}}p\left(t,x|\tau,{{x}_{0}}\right)={{\alpha}_{\tau}}\left({{x}_{0}}\right){{\partial}_{{{x}_{0}}}}\,p+\frac{1}{2}{{\beta}_{\tau}}\left({{x}_{0}}\right)\partial _{{{x}_{0}}}^{2}\,p.\end{eqnarray} \tag{ 2 }$$

The forward and backward Fokker–Planck equations provide information about the conditional probability distribution but not about the individual paths of a Brownian particle. The general theory of how partial differential equations connect to the individual sample paths of a stochastic process goes back to works of Feynman and Kac [85, 107] ² . Their theory allows us to write the solution of the backward Fokker–Planck equation (2) in terms of the following Kolmogorov formula, which constitutes a special case of the Feynman–Kac formula [108–110]:

$$\begin{eqnarray}&&p\left(t,x|\tau,{{x}_{0}}\right)={{\langle \langle \delta \left(x-x(t)\right)\rangle \rangle}_{W}}.\end{eqnarray} \tag{ 3 }$$

The brackets ${{\langle \langle \centerdot \rangle \rangle}_{W}}$ represent an average over realizations of a Wiener process W, which evolves through uncorrelated Gaussian increments $\text{d}W$. The Wiener process drives the evolution of the sample path x(s) from $x(\tau )={{x}_{0}}$ to x(t) via the Itô stochastic differential equation (SDE) [86–88]

$$\begin{eqnarray}&&\text{d}x(s)={{\alpha}_{s}}\left(x(s)\right)\text{d}s+\sqrt{{{\beta}_{s}}\left(x(s)\right)}\text{d}W(s).\end{eqnarray} \tag{ 4 }$$

The diffusion coefficient ${{\beta}_{s}}$ must be non-negative because x(s) describes the position of a real particle. Otherwise, the sample path heads off into imaginary space (for a multivariate process, $\sqrt{{{\beta}_{s}}}$ may be chosen as the unique symmetric and positive-semidefinite square root of ${{\beta}_{s}}$ [111]). Algorithms for the numerical solution of SDEs are provided in [108]. In the physical sciences, SDEs are often written as Langevin equations [89] ³ . For a discussion of stochastic differential equations the reader may refer to a recent report on progress [112].

After this brief detour to continuous-time Markov processes with continuous sample paths, let us return to jump processes, whose sample paths are discontinuous. A system that can be modelled as such a process are motor proteins on cytoskeletal filaments [113–115]. The uni-directional walk of a molecular motor such as myosin, kinesin, or dynein along an actin filament or a microtubule is driven by the hydrolysis of adenosine triphosphate (ATP) and is intrinsically stochastic [116]. Once a sufficient amount of energy is available, one of the two 'heads' of the motor unbinds from its current binding site on the filament and moves to the next binding site. Each binding site can only be occupied by a single head. On a coarse-grained level, the state of the system at time τ is therefore characterized by the occupation of its binding sites. With only a single cytoskeletal filament whose binding sites are labelled as $\left\{0,1,2,\ldots \right\}=:\mathbb{L}$, the variable $n\equiv \left({{n}_{0}},{{n}_{1}},{{n}_{2}},\ldots \right)\in {{\left\{0,1\right\}}^{|\mathbb{L}|}}$ can be used to represent the occupied and unoccupied binding sites. Here, $|\mathbb{L}|$ denotes the total number of binding sites along the filament and n_i   =  1 signifies that the ith binding site is occupied. Since the state space of all the binding site configurations is discrete, a change in the binding site configuration involves a 'jump' in state space. Provided that the jumps are uncorrelated from one another (which needs to be verified experimentally), the dynamics of the system can be described by a continuous-time Markov process with discontinuous sample paths. Before addressing further systems for which this is the case, let us derive the fundamental equations obeyed by these processes: the forward and the backward master equation.

In his classic textbook [110], Gardiner presents a succinct derivation of both the (forward) master and the (forward) Fokker–Planck equation by distinguishing between discontinuous and continuous contributions to sample paths. In the following, we are only interested in the master equation, which governs the evolution of systems whose states change discontinuously. To prevent the occurrence of continuous changes, we assume that the state of our system is represented by a discrete variable n and that the space of all states is countable. With the state space being chosen as the set of integers $\mathbb{Z}$, n could, for example, represent the position of a molecular motor along a cytoskeletal filament. On the other hand, a non-negative integer $n\in {{\mathbb{N}}_{0}}$ could represent the number of molecules in a chemical reaction. The minimal jump size is one in both cases. By keeping the explicit role of n unspecified, the following considerations also apply to systems harbouring different kind of molecules (e.g. $n\equiv \left({{n}^{A}},{{n}^{B}},{{n}^{C}}\right)\in \mathbb{N}_{0}^{3}$), and to systems whose molecules perform random walks in a (discrete) space (e.g. $n\equiv {{\left\{{{n}_{i}}\in {{\mathbb{N}}_{0}}\right\}}_{i\in \mathbb{Z}}}$).

To derive the master equation, we start out by marginalizing the joint conditional distribution $p\left(t,n;\tau,m|{{t}_{0}},{{n}_{0}}\right)$ over the state m at the intermediate time τ ($t\geqslant \tau \geqslant {{t}_{0}}$), resulting in

$$\begin{eqnarray}&&p\left(t,n|{{t}_{0}},{{n}_{0}}\right)=\underset{m}{\sum}\,p\left(t,n;\tau,m|{{t}_{0}},{{n}_{0}}\right).\end{eqnarray} \tag{ 5 }$$

Whenever the range of a sum is not specified, it shall cover the whole state space of its summation variable. The above equation holds for arbitrary stochastic processes. But for a Markov process, one can employ the relation between joint and conditional distributions to turn the equation into the Chapman–Kolmogorov equation

$$\begin{eqnarray}&&p\left(t,n|{{t}_{0}},{{n}_{0}}\right)=\underset{m}{\sum}\,p(t,n|\tau,m)p\left(\tau,m|{{t}_{0}},{{n}_{0}}\right).\end{eqnarray} \tag{ 6 }$$

Letting $p\left(t|{{t}_{0}}\right)$ denote the matrix with elements $p\left(t,n|{{t}_{0}},{{n}_{0}}\right)$, the Chapman–Kolmogorov equation can also be written as $p\left(t|{{t}_{0}}\right)=p(t|\tau )p\left(\tau |{{t}_{0}}\right)$ (semigroup property). Note that the matrix notation requires a mapping between the state space of n and n₀ and an appropriate index set $I\subset \mathbb{N}$. However, we also make use of this notation when the state space is countably infinite.

To derive the (forward) master equation from the Chapman–Kolmogorov equation (6), we define

$$\begin{eqnarray}&&{{Q}_{\tau, \Delta t}}(n,m):=\frac{p(\tau + \Delta t,n|\tau,m)-{{\delta}_{n,m}}}{ \Delta t}\end{eqnarray} \tag{ 7 }$$

for all values of n and m and assume the existence and finiteness of the limits

$$\begin{eqnarray}&&{{Q}_{\tau}}(n,m):=\underset{ \Delta t\to 0}{{\lim}}\,{{Q}_{\tau, \Delta t}}(n,m).\end{eqnarray} \tag{ 8 }$$

These are the elements of the transition (rate) matrix ${{Q}_{\tau}}$, which is also called the infinitesimal generator of the Markov process or is simply referred to as the ${{Q}_{\tau}}$-matrix. Its off-diagonal elements ${{w}_{\tau}}(n,m):={{Q}_{\tau}}(n,m)$ denote the rates at which probability flows from a state m to a state $n\ne m$. The 'exit rates' ${{w}_{\tau}}(m):=-{{Q}_{\tau}}(m,m)$, on the other hand, describe the rates at which probability leaves state m. Both ${{w}_{\tau}}(n,m)$ and ${{w}_{\tau}}(m)$ are non-negative for all n and m. All of the processes considered here shall conserve the total probability, requiring that ${\sum}_{n}{{Q}_{\tau}}(n,m)=0$ or, equivalently, ${{w}_{\tau}}(m)={\sum}_{n}{{w}_{\tau}}(n,m)$ (with ${{w}_{\tau}}(m,m):=0$). The finiteness of the exit rate ${{w}_{\tau}}(m)$ and the conservation of total probability imply that we consider a stable and conservative Markov process [117]. In the natural sciences, the master equation is commonly written in terms of ${{w}_{\tau}}$, but most mathematicians prefer ${{Q}_{\tau}}$. These matrices can be converted into one another by employing

$$\begin{eqnarray}&&{{Q}_{\tau}}(n,m)={{w}_{\tau}}(n,m)-{{\delta}_{n,m}}{{w}_{\tau}}(m).\end{eqnarray} \tag{ 9 }$$

We refer to both of the matrices as transition (rate) matrices and to their (identical) off-diagonal elements as transition rates. The transition rates fully specify the stochastic process.

Assuming that the limit in (8) interchanges with a sum over the state m, the (forward) master equation now follows from the Chapman–Kolmogorov equation (6) as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\partial}_{\tau}}p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right) \\ {} & =\underset{ \Delta t\to 0}{{\lim}}\,\frac{p\left(\tau + \Delta t,n|{{t}_{0}},{{n}_{0}}\right)-p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)}{ \Delta t} \\ {} & =\underset{ \Delta t\to 0}{{\lim}}\,\underset{m}{\sum}\,\frac{p(\tau + \Delta t,n|\tau,m)-{{\delta}_{n,m}}}{ \Delta t}p\left(\tau,m|{{t}_{0}},{{n}_{0}}\right) \\ {} & =\underset{m}{\sum}\,{{Q}_{\tau}}(n,m)p\left(\tau,m|{{t}_{0}},{{n}_{0}}\right). \end{array}\end{eqnarray} \tag{ 10 }$$

Thus, the master equation constitutes a set of coupled, linear, first-order ordinary differential equations (ODEs). The time evolution of the distribution starts out from $p\left({{t}_{0}},n|{{t}_{0}},{{n}_{0}}\right)={{\delta}_{n,{{n}_{0}}}}$. In matrix notation, the equation can be written as ${{\partial}_{\tau}}p\left(\tau |{{t}_{0}}\right)={{Q}_{\tau}}p\left(\tau |{{t}_{0}}\right)$. In terms of ${{w}_{\tau}}$, it assumes the intuitive gain-loss form

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\partial}_{\tau}}p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right) \\ {} & =\underset{m}{\sum}\,\left[{{w}_{\tau}}(n,m)p(\tau,m|\centerdot )-{{w}_{\tau}}(m,n)p(\tau,n|\centerdot )\right]. \end{array}\end{eqnarray} \tag{ 11 }$$

The dot inside the probability distribution's argument abbreviates the initial parameters t₀ and n₀, which are of secondary concern right here. That will change below in the derivation of the backward master equation. An omission of the parameters also makes it impossible to distinguish the conditional distribution $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$ from the single-time distribution $p(\tau,n)$. The single-time distribution obeys the master equation as well, as can inferred directly from the relation $p(\tau,n)={\sum}_{{{n}_{0}}}p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)p\left({{t}_{0}},{{n}_{0}}\right)$ or by summing the above master equation over an initial distribution $p\left({{t}_{0}},{{n}_{0}}\right)$. In fact, the single-time distribution would even obey the master equation if the process was not Markovian, but without providing a complete characterization of the process [70, 74]. The master equation (10) or (11) is particularly interesting for transition rates causing an imbalance between forward and backward transitions along closed cycles of states, i.e. for rates violating Kolmogorov's criterion [118] for detailed balance [70]. Such systems are truly out of thermal equilibrium. If detailed balance is instead fulfilled, the system eventually converges to a stationary Boltzmann–Gibbs distribution with vanishing probability currents between states [70]. Whether or not detailed balance is actually fulfilled is, however, not relevant for the methods discussed in this review. Information on the existence and uniqueness of an asymptotic stationary distribution of the master equation can be found in [117].

The name 'master equation' was originally coined by Nordsieck, Lamb, and Uhlenbeck [5] in their study of the Furry model of cosmic rain showers [119]. Shortly before, Feller applied an equation of the same structure to the growth of populations [120] and Delbrück to well-mixed, auto-catalytic chemical reactions [121]. Delbrück's line of research was followed by several others [122–125], most notably by McQuarrie [126–128] (see also the books [98, 110, 129]). In these articles, several elementary chemical reactions are solved by methods that also appear later in this review. When the particles engaging in a reaction can also diffuse in space, their density may exhibit dynamics that are not expected from observations made in well-mixed environments. Hence, reaction–diffusion master equations have been the focus of intense research and have been analysed using path integrals (see, for example, [96, 130–133] and the references in section 1.5). Master equations, and simulations algorithms based on master equations, are now being used in numerous fields of research. They are being applied in the contexts of spin dynamics [134–137], gene regulatory networks [34, 138–143], the spreading of diseases [144–147], epidermal homeostasis [148], nucleosome repositioning [149], ecological [150–163] and bacterial dynamics [158, 164–167], evolutionary game theory [168–179], surface growth [180], and social and economic processes [181–184]. Queuing processes are also often modelled in terms of master equations, but in this context, the equations are typically referred to as Kolmogorov equations [185]. Moreover, master equations and the SSA have helped to understand the formation of traffic jams on highways [186, 187], the walks of molecular motors along cytoskeletal filaments [113–115, 188–190], and the condensation of bosons in driven-dissipative quantum systems [178, 191, 192]. The master equation that was found to describe the coarse-grained dynamics of these bosons coincides with the master equation of the (asymmetric) inclusion process [193–196]. Transport processes are commonly modelled in terms of the (totally) asymmetric simple exclusion process (ASEP or TASEP) [197–200]. The ASEP describes the biased hopping of particles along a one-dimensional lattice, with each lattice site providing space for at most one particle. The ASEP and the TASEP are regarded as paradigmatic models in the field of non-equilibrium statistical mechanics, with many exact mathematical results having been established [201–212]. Some of these results were established by applying the Bethe ansatz to the master equation of the ASEP [205, 209]. The review of Chou, Mallick, and Zia provides a comprehensive account of the ASEP and of its variants [190]. The master equation of the TASEP with Langmuir kinetics was recently used to understand the length regulation of microtubules [213].

Unlike deterministic models, the master equations describing the dynamics of the above systems take into account that discrete and finite populations are prone to 'demographic fluctuations'. The populations of the above systems consist of genes or proteins, infected persons, bacteria or cars and they are typically small, at least compared to the number of molecules in a mole of gas. For example, the copy number of low abundance proteins in Escherichia coli cytosol was found to be in the tens to hundreds [214]. Therefore, the presence or absence of a single protein is much more important than the presence or absence of an individual molecule in a mole of gas. A demographic fluctuation may even be fatal for a system, for example, when the copy number of an auto-catalytic reactant drops to zero. The master equation (10) provides a useful tool to describe such an effect.

Up to this point, we have only considered the forward master equation. But just as the (forward) Fokker–Planck equation (1) is complemented by the backward Fokker–Planck equation (2), the (forward) master equation (10) is complemented by a backward master equation. This equation can be derived from the Chapman–Kolmogorov equation (6) as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\partial}_{-\tau}}p\left(t,n|\tau,{{n}_{0}}\right) \\ {} & =\underset{ \Delta t\to 0}{{\lim}}\,\frac{p\left(t,n|\tau - \Delta t,{{n}_{0}}\right)-p\left(t,n|\tau,{{n}_{0}}\right)}{ \Delta t} \\ {} & =\underset{ \Delta t\to 0}{{\lim}}\,\underset{m}{\sum}\,p(t,n|\tau,m)\frac{p\left(\tau,m|\tau - \Delta t,{{n}_{0}}\right)-{{\delta}_{m,{{n}_{0}}}}}{ \Delta t} \\ {} & =\underset{m}{\sum}\,p(t,n|\tau,m){{Q}_{\tau}}\left(m,{{n}_{0}}\right). \end{array}\end{eqnarray} \tag{ 12 }$$

Here, the transition rate is obtained in the limit ${{\lim}_{ \Delta t\to 0}}{{Q}_{\tau - \Delta t, \Delta t}}\left(m,{{n}_{0}}\right)$ (see (7)). In matrix notation, the backward master equation reads ${{\partial}_{-\tau}}p(t|\tau )=p(t|\tau ){{Q}_{\tau}}$. In terms of ${{w}_{\tau}}$, it assumes the form (see (9))

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\partial}_{-\tau}}p\left(t,n|\tau,{{n}_{0}}\right) \\ {} & =\underset{m}{\sum}\,\left[p(\centerdot |\tau,m)-p\left(\centerdot |\tau,{{n}_{0}}\right)\right]{{w}_{\tau}}\left(m,{{n}_{0}}\right). \end{array}\end{eqnarray} \tag{ 13 }$$

In this equation, the dots abbreviate the final parameters t and n. The backward master equation evolves the conditional probability distribution backward in time, starting out from the final condition $p\left(t,n|t,{{n}_{0}}\right)={{\delta}_{n,{{n}_{0}}}}$. Just as the backward Fokker–Planck equation, the backward master equation proves useful for the computation of mean extinction and first passage times (see [110, 215] and section 3.3). Furthermore, it follows from the backward master equation (12) that the (conditional) average $\langle A\rangle \left(t|\tau,{{n}_{0}}\right):={\sum}_{n}A(n)p\left(t,n|\tau,{{n}_{0}}\right)$ of an observable A fulfils an equation of just the same form, namely

$$\begin{eqnarray}&&{{\partial}_{-\tau}}\langle A\rangle \left(t|\tau,{{n}_{0}}\right)=\underset{m}{\sum}\,\langle A\rangle (t|\tau,m){{Q}_{\tau}}\left(m,{{n}_{0}}\right).\end{eqnarray} \tag{ 14 }$$

The final condition of the equation is given by $\langle A\rangle \left(t|t,{{n}_{0}}\right)=A\left({{n}_{0}}\right)$. The validity of equation (14) is the reason why we later employ a 'backward' path integral to represent the average $\langle A\rangle$ (see section 6).

If the dynamics of a system are restricted to a finite number of states and if its transition rates are independent of time, both the forward master equation ${{\partial}_{\tau}}p\left(\tau |{{t}_{0}}\right)=Qp\left(\tau |{{t}_{0}}\right)$ and the backward master equation ${{\partial}_{-{{\tau}_{0\,}}}}p\left(t|{{\tau}_{0}}\right)=p\left(t|{{\tau}_{0}}\right)Q$ are solved by [216]

$$\begin{eqnarray}&&p\left(\tau |{{\tau}_{0}}\right)={{\text{e}}^{Q\left(\tau -{{\tau}_{0}}\right)}}{\mathbb {1}}\end{eqnarray} \tag{ 15 }$$

(recall that $p\left(\tau |{{\tau}_{0}}\right)$ is the matrix with elements $p\left(\tau,n|{{\tau}_{0}},{{n}_{0}}\right)$). The Chapman–Kolmogorov equation (6) is also solved by the distribution. Although the matrix exponential inside this solution can in principle be evaluated in terms of the (convergent) Taylor expansion ${\sum}_{k=0}^{\infty}\frac{{{\left(\tau -{{\tau}_{0}}\right)}^{k}}}{k!}{{Q}^{k}}$, the actual calculation of this series is typically infeasible for non-trivial processes, both analytically and numerically (a truncation of the Taylor series may induce severe round-off errors and serves as a lower bound on the performance of algorithms in [217]). Consequently, alternative numerical algorithms have been developed to evaluate the matrix exponential. Moler and Van Loan reviewed 'nineteen dubious ways' of computing the exponential in [217, 218]. Algorithms that can deal with very large state spaces are considered in [219, 220]. For time-dependent transition rates, the matrix exponential generalizes to a Magnus expansion [221, 222].

In the previous paragraph, we restricted the dynamics to a finite state space to ensure the existence of the matrix exponential in (15). Provided that the supremum ${{\sup}_{m}}|Q(m,m)|$ of all the exit rates is finite (uniformly bounded Q-matrix), the validity of the above solution extends to state spaces comprising a countable number of states [223]. To see that, we define the left stochastic matrix $P:={\mathbb {1}}+{{\lambda}^{-1}}Q$, with the parameter λ being larger than the above supremum. Writing ${{\text{e}}^{Q \Delta t}}={{\text{e}}^{-\lambda \Delta t}}{{\text{e}}^{\lambda \Delta tP}}$ with $\Delta t:=\tau -{{\tau}_{0}}$, the matrix exponential can be evaluated in terms of the convergent Taylor series

$$\begin{eqnarray}&&p\left(\tau |{{\tau}_{0}}\right)=\underset{k=0}{\overset{\infty}{\sum}}\,\frac{{{\text{e}}^{-\lambda \Delta t}}{{(\lambda \Delta t)}^{k}}}{k!}{{P}^{k}}.\end{eqnarray} \tag{ 16 }$$

Effectively, one has thereby decomposed the continuous-time Markov process with transition matrix Q into a discrete-time Markov chain with transition matrix P, subordinated to a continuous-time Poisson process with rate coefficient λ (the Poisson process acts as a 'clock' with sufficiently high ticking rate λ). Such a decomposition is called a uniformization or randomization and was first proposed by Jensen [224]. The series (16) can be evaluated via numerically stable algorithms and truncation errors can be bounded [225, 226]. Nevertheless, the uniformization method requires the computation of the powers of a matrix having as many rows and columns as the system has states. Consequently, a numerical implementation of the method is only feasible for sufficiently small state spaces. Further information on the method and on its improvements can be found in [224–229].

The mathematical study of the existence and uniqueness of solutions of the forward and backward master equations was pioneered by Feller and Doob in the 1940s [230, 231]. Feller derived an integral recurrence formula [117, 230], which essentially constitutes a single step of the 'path summation representation' that we derive further below. In the following, we assume that the forward and the backward master equations have the same unique solution and we restrict our attention to processes performing only a finite number of jumps during any finite time interval. These conditions, and the conservation of total probability, are, for example, violated by processes that 'explode' after a finite time ⁴ . More information on such processes is provided in [117, 216].

In the following, we complement the above representations of the master equations's solution with a 'path summation representation'. This representation can be derived by examining the steps of the stochastic simulation algorithm (SSA) of Gillespie (its 'direct' version) [75, 76, 78] or by performing a Laplace transformation of the forward master equation. Here we follow the former, qualitative, approach. A formal derivation of the representation via the master equation's Laplace transform is provided in appendix B. Although the basic elements of the SSA had already been known before Gillespie's work [230–236], its popularity largely increased after Gillespie applied it to the study of chemical reactions. As the SSA is restricted to time-independent transition rates, so is the following derivation.

To derive the path summation representation, we prepare a system in state n₀ at time t₀ as illustrated in figure 2. Since the process is homogeneous in time, we may choose t₀  =  0. The total rate of leaving state n₀ is given by the exit rate $w\left({{n}_{0}}\right)={\sum}_{{{n}_{1}}}w\left({{n}_{1}},{{n}_{0}}\right)$. In order to determine how long the system actually stays in state n₀, one may draw a random waiting time ${{\tau}_{0}}$ from the exponential distribution ${{\mathcal{W}}_{{{n}_{0}}}}\left({{\tau}_{0}}\right):=w\left({{n}_{0}}\right){{\text{e}}^{-w\left({{n}_{0}}\right){{\tau}_{0}}}} \Theta \left({{\tau}_{0}}\right)$. The Heaviside step function $\Theta$ prevents the sampling of negative waiting times and is here defined as $\Theta (\tau )=1$ for $\tau \geqslant 0$ and $\Theta (\tau )=0$ for $\tau <0$. Thus far, we only know that the system leaves n₀ but not where it ends up. It could end up in any state n₁ for which the transition rate $w\left({{n}_{1}},{{n}_{0}}\right)$ is positive. The probability that a particular state n₁ is realized is given by $w\left({{n}_{1}},{{n}_{0}}\right)/w\left({{n}_{0}}\right)$. In a numerical implementation of the SSA, the state n₁ is determined by drawing a second (uniformly-distributed) random number. Our goal is to derive an analytic representation of the probability $p\left(t,n|{{t}_{0}},{{n}_{0}}\right)$ of finding the system in state n at time t. Thus, after taking J  −  1 further steps, the sample path ${{\mathcal{P}}_{J}}:=\left\{{{n}_{J}}\leftarrow \cdots \leftarrow {{n}_{1}}\leftarrow {{n}_{0}}\right\}$ should eventually visit state n_J   =  n at some time ${{t}_{J}}\leqslant t$. The total time ${{\tau}_{J-1}}+\cdots +{{\tau}_{0}}$ until the jump to state n occurs is distributed by the convolutions of the individual waiting time distributions, i.e. by $\underset{J-1}{\overset{j=0}{{\star}}}\,{{\mathcal{W}}_{{{n}_{j}}}}$. For example, $\tau :={{\tau}_{1}}+{{\tau}_{0}}$ is distributed by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \left({{\mathcal{W}}_{{{n}_{1}}}} \star {{\mathcal{W}}_{{{n}_{0}}}}\right)(\tau )={\int}_{\mathbb{R}}\text{d}{{\tau}_{0}}\,{{\mathcal{W}}_{{{n}_{1}}}}\left(\tau -{{\tau}_{0}}\right){{\mathcal{W}}_{{{n}_{0}}}}\left({{\tau}_{0}}\right) \end{array}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \qquad \\ \;=\frac{w\left({{n}_{1}}\right)w\left({{n}_{0}}\right)\left({{\text{e}}^{-w\left({{n}_{0}}\right)\tau}}-{{\text{e}}^{-w\left({{n}_{1}}\right)\tau}}\right)}{w\left({{n}_{1}}\right)-w\left({{n}_{0}}\right)} \Theta (\tau ). \end{array}\end{eqnarray} \tag{ 18 }$$

Zoom In Zoom Out Reset image size

The probability that the system still resides in state n_J   =  n at time t is determined by the 'survival probability' ${{\mathcal{S}}_{n}}\left(t-{{t}_{J}}\right)={{\text{e}}^{-w(n)\left(t-{{t}_{J}}\right)}} \Theta \left(t-{{t}_{J}}\right)$. After putting all of these pieces together, we arrive at the following path summation representation of the conditional probability distribution:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} p\left(t,n|{{t}_{0}},{{n}_{0}}\right)=\underset{J=0}{\overset{\infty}{\sum}}\,\underset{\left\{{{\mathcal{P}}_{J}}\right\}}{\sum}\,{{p}_{t-{{t}_{0}}}}\left({{\mathcal{P}}_{J}}\right)~\text{with} \end{array}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{\tau}}\left({{\mathcal{P}}_{J}}\right)=\left({{\mathcal{S}}_{n}}\star \left(\underset{J-1}{\overset{j=0}{{\star}}}\,\frac{w\left({{n}_{j+1}},{{n}_{j}}\right)}{w\left({{n}_{j}}\right)}{{\mathcal{W}}_{{{n}_{j}}}}\right)\right)(\tau ). \end{array}\end{eqnarray} \tag{ 20 }$$

Here, ${\sum}_{\left\{{{\mathcal{P}}_{J}}\right\}}:={\sum}_{{{n}_{1}}}\cdots {\sum}_{{{n}_{J-1}}}$ generates every path with J jumps between n₀ and n_J   =  n. The probability of such a path, integrated over all possible waiting times, is represented by ${{p}_{\tau}}\left({{\mathcal{P}}_{J}}\right)$. By an appropriate choice of integration variables, the probability ${{p}_{\tau}}\left({{\mathcal{P}}_{J}}\right)$ can also be written as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{p}_{\tau}}\left({{\mathcal{P}}_{J}}\right)= & \left(\underset{j=0}{\overset{J-1}{\prod}}\,{\int}_{0}^{\tau}\text{d}{{\tau}_{j}}\,\frac{w\left({{n}_{j+1}},{{n}_{j}}\right)}{w\left({{n}_{j}}\right)}{{\mathcal{W}}_{{{n}_{j}}}}\left({{\tau}_{j}}\right)\right) \\ {} & \centerdot \,\,{{\mathcal{S}}_{n}}\left(\tau -\left({{\tau}_{J-1}}+\cdots +{{\tau}_{0}}\right)\right). \end{array}\end{eqnarray} \tag{ 21 }$$

The survival probability ${{\mathcal{S}}_{n}}$ is included in these integrations. Without the integrations, the expression would represent the probability of a path with J jumps and fixed waiting times. That probability is, for example, used in the master equation formulation of stochastic thermodynamics in associating an entropy to individual paths [237]. In appendix B, we formally derive the path summation representation (19) from the Laplace transform of the forward master equation (11).

The path summation representation (19) does not only form the basis of the SSA but also of some alternative algorithms [238–243]. These algorithms either infer the path probability ${{p}_{\tau}}\left({{\mathcal{P}}_{J}}\right)$ numerically from its Laplace transform or evaluate the convolutions in (20) analytically. The analytic expressions that arise are, however, rather cumbersome generalizations of the convolution in (18) [244, 245]. They simplify only for the most basic processes (e.g. for a random walk or for the Poisson process). Moreover, care has to be taken when the analytic expressions are evaluated numerically because they involve pairwise differences of the exit rate w(n) (see (18)). When these exit rates differ only slightly along a path, a substantial loss of numerical significance may occur due to finite precision arithmetic. Future studies could explore how the convolutions of exponential distributions in (20) can be approximated efficiently (for example, in terms of a Gamma distribution or by analytically determining the Laplace transform of (20), followed by a saddle-point approximation [246] of the corresponding inverse Laplace transformation). In general, both the SSA as well as its competitors suffer from the enormous number of states of non-trivial systems, as well as from the even larger number of paths connecting different states. In [240], these paths were generated using a deterministic depth-first search algorithm, combined with a filter to select the paths that arrive at the right place at the right time. In [242], a single path was first generated using the SSA and then gradually changed into new paths through a Metropolis Monte Carlo scheme. Thus far, the two methods have only been applied to relatively simple systems and their prevalence is low compared to the prevalence of the SSA. Further research is needed to explore how relevant paths can be sampled more efficiently.

The true power of the SSA lies in its generation of sample paths with the correct probability of occurrence. Thus, just a few sample paths generated with the SSA are often sufficient to infer the 'typical' dynamics of a process. A look at individual paths may, for example, reveal that the dynamics of a system are dominated by some spatial pattern, e.g. by spirals [168]. Efficient variations of the above 'direct' version of the SSA are, for example, described in [78, 247–249]. Algorithms for the fast simulation of biochemical networks or processes with spatial degrees of freedom are implemented in the simulation packages [250–257].

The evaluation of the average $\langle A\rangle ={\sum}_{n}A(n)p(t,n|\centerdot )$ of an observable A typically requires the computation of a larger number of sample paths. However, since the occurrence probability of sample paths generated with the SSA is statistically correct, such an average typically converges comparatively fast. Furthermore, each path can be sampled independently of every other path. Therefore, the computation of paths can be distributed to individual processing units, saving real time, albeit no computation time. A distributed computation of the sample paths is most often required, but possibly not even sufficient, if one wishes to compute the full probability distribution $p\left(t,n|{{t}_{0}},{{n}_{0}}\right)$. Vastly more sample paths are required for this purpose, especially if 'rare event probabilities' in the distribution's tails are sought for. In particular, if the probability of finding a system in state n at time t is only 10⁻¹⁰, an average of 10¹⁰ sample paths are needed to observe that event just once. Moreover, the probability of observing any particular state decreases with the size of a system's state space. Thus, the sampling of full distributions becomes less and less feasible as systems become larger. Various other challenges remain open as well; for example, the efficient simulation of processes evolving on multiple time scales. These processes are typically simulated using approximative techniques such as τ-leaping [78, 249, 258–267]. Another challenge is posed by the evaluation of processes with time-dependent transition rates [268–270].

For completeness, let us mention yet another numerical approach to the (forward) master equation. Since the master equation (10) constitutes a set of coupled linear first-order ODEs, it can of course be treated as such and be integrated numerically. The integration is, however, only feasible if the state space is sufficiently small (or appropriately truncated) and if all transitions occur on comparable time scales (otherwise, the master equation is quite probably stiff [271]). Nevertheless, a numerical integration of the master equation may be preferable over the use of the SSA if the full probability distribution is to be computed.

Neither the matrix exponential representation $p\left(t|{{t}_{0}}\right)={{\text{e}}^{Q\left(t-{{t}_{0}}\right)}}{\mathbb {1}}$ of the conditional probability distribution, nor its path summation representation (19) is universally applicable. Moreover, even if the requirements of these solutions are met, the size of the state space or the complexity of the transition matrix may make it infeasible to evaluate them. In the next sections, we formulate conditions under which the conditional probability distribution can be represented in terms of the 'forward' path integral

$$\begin{eqnarray}&&p\left(t,n|{{t}_{0}},{{n}_{0}}\right)={{\langle n{{|}_{t}}{\rlap{-} \int}_{[{{t}_{0}}}^{t)}{{\text{e}}^{-\mathcal{S}}}\,|{{n}_{0}}\rangle}_{{{t}_{0}}}}\end{eqnarray} \tag{ 22 }$$

and in terms of the 'backward' path integral

$$\begin{eqnarray}&&p\left(t,n|{{t}_{0}},{{n}_{0}}\right)={{\langle {{n}_{0}}{{|}_{{{t}_{0}}}}{\rlap{-} \int}_{({{t}_{0}}}^{t]}\,{{\text{e}}^{-{{\mathcal{S}}^{\dagger}}}}\,|n\rangle}_{t}}.\end{eqnarray} \tag{ 23 }$$

The meaning of the integral signs and of the bras $\langle n|$ and kets $|n\rangle$ will become clear over the course of this review. Let us only note that the integrals do not proceed along paths of the discrete variable n, but over the paths of two continuous auxiliary variables that are introduced for this purpose. The relevance of each path is weighed by the exponential factors inside the integrals.

Besides these exact representations of the conditional probability distribution solving the master equations, there exist powerful ways of approximating this distribution and the values of averaged observables. These methods include the Kramers–Moyal [103, 272] and the system-size expansion [98, 273], as well as the derivation of moment equations. Information on these methods can be found in classic text books [98, 110] and in a recent review [274]. Although moment equations encode the complete information about a stochastic process, they typically constitute an infinite hierarchy whose evaluation requires a truncation by some closure scheme [275–283]. On the other hand, the Kramers–Moyal and the system-size expansion approximate the master equation in terms of a Fokker–Planck equation. Both expansions work best if the system under consideration is 'large' (more precisely, they work best if the dynamics are centred around a stable or meta-stable state at a distance $N\gg 1$ from a potentially absorbing state; the standard deviation of its surrounding distribution is then of order $\sqrt{N}$). An extension of the system-size expansion to absorbing boundaries has recently been proposed in [284]. In sections 4.4 and 5.3, we show how Kramers–Moyal expansions of the backward and forward master equations can be used to derive path integral representations of processes with continuous state spaces. When the expansion of the backward master equation stops (or is truncated) at the level of a diffusion approximation, one recovers a classic path integral representation of the (backward) Fokkers–Planck equation [16–19].

The oldest path integral, both in the theory of stochastic processes and beyond [107, 285–287], is presumably Wiener's integral for Brownian motion [80, 81]. The path integrals we consider here were devised somewhat later, namely in the 1970s and 80s: first for the Fokker–Planck (or Langevin) equation [16–19] and soon after for the master equation. For the master equation, the theoretical basis of these 'stochastic' path integrals was developed by Doi [6, 7]. He first expressed the creation and annihilation of molecules in a chemical reaction by the corresponding operators for the quantum harmonic oscillator [288] (modulo a normalization factor), introducing the concept of the Fock space [20] into non-equilibrium theory. Furthermore, he employed coherent states [22–24] to express averaged observables. Similar formalisms as the one of Doi were independently developed by Zel'dovich and Ovchinnikov [8], as well as by Grassberger and Scheunert [10]. Introductions to the Fock space approach, which are in part complementary to our review, are, for example, provided in [70, 91, 289–291]. The review of Mattis and Glasser [289] provides a chronological list of contributions to the field. These contributions include Rose's renormalized kinetic theory [9], Mikhailov's path integral [11, 12, 14], which is based on Zel'dovich's and Ovchinnikov's work, and Goldenfeld's extension of the Fock space algebra to polymer crystal growth [13]. Furthermore, Peliti reviewed the Fock space approach and provided derivations of path integral representations of averaged observables and of the probability generating function [15]. Peliti also expressed the hope that future 'rediscoveries' of the path integral formalism would be unnecessary in the future. However, we believe that the probabilistic structures behind path integral representations of stochastic processes have not yet been clearly exposed. As illustrated in figure 1, we show that the forward and backward master equations admit not only one but two path integral representations: the forward representation (22) and the novel backward representation (23). Although the two path integrals resemble each other, they differ conceptually. While the forward path integral representation provides a probability generating function in an intermediate step, the backward representation provides a distribution that is marginalized over an initial state. Both path integrals can be used to represent averaged observables as shown in section 6. The backward path integral, however, will turn out to be more convenient for this purpose (upon choosing a Poisson distribution as the basis function, i.e. $|n\rangle :=\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}$, the representation is obtained by summing (23) over an observable A(n)). Let us note that even though we adopt some of the notation of quantum theory, our review is guided by the notion that quantum (field) theory is 'totally unnecessary' for the theory of stochastic path integrals. Michael E Fisher once stated the same about the relevance of quantum field theory for renormalization group theory [292] (while acknowledging that in order to do certain types of calculations, familiarity with quantum field theory can be very useful; the same applies to the theory of stochastic path integrals).

Thus far, path integral representations of the master equation have primarily been applied to processes whose transition rates can be decomposed additively into the rates of simple chemical reactions. Simple means that the transition rate of such a reaction is determined by combinatoric counting. Consider, for example, a reaction of the form $k\,A\to l\,A$ in which $k\in {{\mathbb{N}}_{0}}$ molecules of type A are replaced by $l\in {{\mathbb{N}}_{0}}$ molecules of the same type. Assuming that the k reactants are drawn from an urn with a total number of m molecules, the total rate of the reaction should be proportional to the falling factorial ${{(m)}_{k}}:=m(m-1)\cdots (m-k+1)$. The global time scale of reaction is set by the rate coefficient ${{\gamma}_{\tau}}$, which we allow to depend on time. Thus, the rate at which the chemical reaction $k\,A\to l\,A$ induces a transition from state m to state n is ${{w}_{\tau}}(n,m)={{\gamma}_{\tau}}{{(m)}_{k}}{{\delta}_{n,m-k+l}}$. The Kronecker delta inside the transition rate ensures that k molecules are indeed replaced by l new ones. Note that the number of particles in the system can never become negative, provided that the initial number of particles was non-negative. Hence, the state space of n is ${{\mathbb{N}}_{0}}$. Insertion of the above transition rate into the forward master equation (11) results in the 'chemical' master equation

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\partial}_{\tau}}p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right) \\ {} & ={{\gamma}_{\tau}}\left[{{(n-l+k)}_{k}}p(\tau,n-l+k|\centerdot )-{{(n)}_{k}}p(\tau,n|\centerdot )\right]. \end{array}\end{eqnarray} \tag{ 24 }$$

Microphysical arguments for its applicability to chemical reactions can be found in [77]. According to the chemical master equation, the mean particle number $\langle n\rangle ={\sum}_{n}n\,p(\tau,n|\centerdot )$ obeys the equation ${{\partial}_{\tau}}\langle n\rangle ={{\gamma}_{\tau}}(l-k)\langle {{(n)}_{k}}\rangle$. For a system with a large number of particles ($n\gg k$), fluctuations can often be neglected in a first approximation, leading to the deterministic rate equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}\bar{n}={{\gamma}_{\tau}}(l-k){{\bar{n}}^{k}},\end{eqnarray} \tag{ 25 }$$

obeyed by a continuous particle number $\bar{n}(\tau )\in \mathbb{R}$.

Path integral representations of the chemical master equation (24) are sometimes said to be 'bosonic'. First, because an arbitrarily large number of molecules may in principle be present in the system (both globally and, upon extending the system to a spatial domain, locally). Second, because the path integral representations are typically derived with the help of 'creation and annihilation operators' fulfilling a 'bosonic' commutation relation (see section 2.2.2). 'Fermionic' path integrals, on the other hand, have been developed for systems in which the particles exclude one another. Thus, the number of particles in these systems is locally restricted to the values 0 and 1. For systems with excluding particles, the master equation's solution may be represented in terms of a path integral whose underlying creation and annihilation operators fulfil an anti-commutation relation [26–31]. However, van Wijland recently showed that the use of operators fulfilling the bosonic commutation relation is also possible [90]. These approaches are considered in section 2.2.4.

We do not intend to delve further into the historic development and applications of stochastic path integrals at this point. Doing so would require a proper introduction into renormalization group theory, which is of pivotal importance for the evaluation of the path integrals. Readers can find information on the application of renormalization group techniques in the review of Täuber, Howard, and Vollmayr-Lee [69] and in the book of Täuber [70]. Introductory texts are also provided by Cardy's (lecture) notes [91, 293]. Roughly speaking, path integral representations of the chemical master equation (24) have been used to assess how a macroscopic law of mass action changes due to fluctuations, both below [96, 130, 132, 294–305] and above the (upper) critical dimension [295, 306, 307], using either perturbative [14, 70, 130–133, 294–301, 303–305, 308–316] or non-perturbative [306, 307, 317–319] techniques. All of these articles focus on stochastic processes with spatial degrees of freedom for which alternative analytical approaches are scarce. Path integral representations of these processes, combined with renormalization group techniques, have been pivotal in understanding non-equilibrium phase transitions and they contributed significantly to the classification of these transitions in terms of universality classes [67, 70, 320]. Moreover, path integral representations of master equations have recently been employed in such diverse contexts as the study of neural networks [321–323], of ecological populations [156, 157, 162, 163, 324, 325], and of the differentiation of stem-cells [326].

Continuous-time Markov processes with discontinuous sample paths describe a broad range of phenomena, ranging from traffic jams on highways [186] and on cytoskeletal filaments [113–115, 190] to novel forms of condensation in bosonic systems [178, 191, 192]. In the introduction, we laid out the mathematical theory of these processes and derived the fundamental equations governing their evolution: the forward and the backward master equations. Whereas the forward master equation (10) evolves a conditional probability distribution forward in time, the backward master equation (12) evolves the distribution backward in time. In the following main part of this review, we represent the conditional probability distribution solving the master equations in terms of path integrals. The framework upon which these path integrals are based unifies a broad range of approaches to the master equations, including, for example, the spectral method of Walzcak, Mugler, and Wiggins [34] and the Poisson representation of Gardiner and Chaturvedi [32, 33].

We now formulate general conditions under which the forward master equation (10) can be cast into a linear PDE obeyed by the generalized probability generating function

$$\begin{eqnarray}&&|g\left(\tau |{{t}_{0}},{{n}_{0}}\right){{\rangle}_{q}}:=\underset{n}{\sum}\,|n{{\rangle}_{\tau,q}}\,p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right).\end{eqnarray} \tag{ 45 }$$

The basis function $|n{{\rangle}_{\tau,q}}$ is a function of the variable q and possibly of the time variable τ. But unless one of these variables is of direct relevance, its corresponding subscript will be dropped. The explicit form of the basis function depends on the problem at hand and is chosen so that the four conditions (O), (C), (E), and (Q) below are satisfied (the conditions (O) and (C) concern the orthogonality and completeness of the basis, which we already required in the introductory example). The variable n again represents some state from a countable state space. For example, n could describe the position of a molecular motor along a cytoskeletal filament ($n\in \mathbb{Z}$), the copy number of a molecule ($n\in {{\mathbb{N}}_{0}}$), the local copy numbers of the molecule on a lattice ($n\equiv {{\left\{{{n}_{i}}\in {{\mathbb{N}}_{0}}\right\}}_{i\in \mathbb{Z}}}$), or the copy numbers of multiple kinds of molecules ($n\equiv \left({{n}^{A}},{{n}^{B}},{{n}^{C}}\right)\in \mathbb{N}_{0}^{3}$). For the multivariate configurations, the basis function is typically decomposed into a product $|{{n}_{1}}\rangle |{{n}_{2}}\rangle \cdots$ of individual basis functions $|{{n}_{i}}\rangle$, each depending on its own variable q_i . A process with spatial degrees of freedom is considered in section 4.5.2. Besides, we also consider a system of excluding particles in section 2.2.4. There, the (local) number of particles n is restricted to the values 0 and 1. Note that the sum in (45) extends over the whole state space.

The definition of the generating function (45) assumes the existence of a set ${{\{|n\rangle}_{\tau}}\}$ of basis functions for every time $\tau \in \left[{{t}_{0}},t\right]$. In addition, we assume that there exists a set $\left\{\langle m{{|}_{\tau}}\right\}$ of linear basis functionals for every time $\tau \in \left[{{t}_{0}},t\right]$. These bras shall be orthogonal to the kets in the sense that at each time point τ, they act on the kets as

$$\begin{eqnarray}&&{{\langle m{{|}_{\tau}}|n\rangle}_{\tau}}={{\delta}_{m,n}}\quad \quad \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left(\text{O}\right)\end{eqnarray}$$

(for all m and n). Here we note the possible time-dependence of the basis because the (O)rthogonality condition will only be required for equal times of the bras and kets. In addition to orthogonality, the basis shall fulfil the (C)ompleteness condition

$$\begin{eqnarray}&&\underset{n}{{\,\,\,{\sum}^{}}}\,|n{{\rangle}_{\tau}}\langle n{{|}_{\tau}}f=f,\quad \quad \quad \quad \quad \quad \quad \,\,\,\,\,\,\,\quad \quad \quad \left(\text{C}\right)\end{eqnarray}$$

where f represents an appropriate test function. The completeness condition implies that the function f can be decomposed in the basis functions $|n\rangle$ with expansion coefficients $\langle n|\,f$. In the introduction to this section, we introduced various bases fulfilling both the orthogonality and the completeness condition. As in the introductory example, the orthogonality condition allows one to recover the conditional probability distribution via the inverse transformation

$$\begin{eqnarray}&&p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)=\langle n|g\left(\tau |{{t}_{0}},{{n}_{0}}\right)\rangle.\end{eqnarray} \tag{ 48 }$$

Before deriving the flow equation obeyed by the generating function, let us note that the (O)rthogonality condition differs slightly from the corresponding conditions used in most other texts on stochastic path integrals (see, for example, [13, 15] or the 'exclusive scalar product' in [10]). Typically, the orthogonality condition includes an additional factorial $n!$ on its right hand side. The inclusion of this factorial is advantageous in that it accentuates a symmetry between the bases that we consider in sections 2.2.2 and 3.2.2 for the study of chemical reactions. Its inclusion would be rather unusual, however, for the Fourier basis introduced in sections 2.2.1 and 3.2.1. The Fourier basis will be used to solve a simple random walk. Moreover, the factorial obscures a connection between the probability generating functional introduced in section 3.4 and the Poisson representation of Gardiner and Chaturvedi [32, 33]. We discuss this connection in section 3.5.

To derive the flow equation obeyed by the generating function $|g\rangle$, we differentiate its definition (45) with respect to the time variable τ. The resulting time derivative of the conditional probability distribution $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$ can be replaced by the right-hand side of the forward master equation (10). In matrix notation, this equation reads ${{\partial}_{\tau}}p\left(\tau |{{t}_{0}}\right)={{Q}_{\tau}}p\left(\tau |{{t}_{0}}\right)$. Eventually, one finds that

$$\begin{eqnarray}&&{{\partial}_{\tau}}|g\rangle =\underset{n}{\sum}\,\left({{\partial}_{\tau}}|n\rangle +\underset{m}{\sum}\,|m\rangle {{Q}_{\tau}}(m,n)\right)\,p(\tau,n|\centerdot ).\end{eqnarray} \tag{ 49 }$$

Our goal is to turn this expression into a partial differential equation for $|g\rangle$. For this purpose, we require two differential operators. First, we require a differential operator ${{\mathcal{E}}_{\tau}}\left(q,{{\partial}_{q}}\right)$ encoding the time evolution of the basis function. In particular, this operator should fulfil, for all values of n,

$$\begin{eqnarray}&&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\mathcal{E}}_{\tau}}|n\rangle ={{\partial}_{\tau}}|n\rangle.\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \left(\text{E}\right)\end{eqnarray}$$

By the (O)rthogonality condition, one could also define this operator in a 'constructive' way as

$$\begin{eqnarray}&&{{\mathcal{E}}_{\tau}}:=\underset{n}{\sum}\,\left({{\partial}_{\tau}}|n\rangle \right)\langle n|.\end{eqnarray} \tag{ 51 }$$

We call ${{\mathcal{E}}_{\tau}}$ the basis (E)volution operator. In order to arrive at a proper PDE for $|g\rangle$, ${{\mathcal{E}}_{\tau}}\left(q,{{\partial}_{q}}\right)$ should be polynomial in ${{\partial}_{q}}$ (later, in section 2.2.4 we also encounter a case in which it constitutes a power series with infinitely high powers of ${{\partial}_{q}}$). For now, the pre-factors of 1, ${{\partial}_{q}}$, $\partial _{q}^{2}$, may be arbitrary functions of q. Later, in our derivation of a path integral in section 5, we will also require that the pre-factors can be expanded in powers of q. Note that for a multivariate configuration $n\equiv \left({{n}^{A}},{{n}^{B}},\ldots \right)$, the derivative ${{\partial}_{q}}$ represents individual derivatives with respect to $q\equiv \left({{q}^{A}},{{q}^{B}},\ldots \right)$.

The actual dynamics of a jump process are encoded by its transition rate matrix ${{Q}_{\tau}}$ (see section 1.3). The off-diagonal elements of this matrix are the transition rates ${{w}_{\tau}}(n,m)$ from a state m to a state n, and its diagonal elements are the negatives of the exit rates ${{w}_{\tau}}(m)={\sum}_{n}{{w}_{\tau}}(n,m)$ from a state m (with ${{w}_{\tau}}(m,m)=0$; see (9)). We encode the information stored in ${{Q}_{\tau}}$ by a second differential operator called ${{\mathcal{Q}}_{\tau}}\left(q,{{\partial}_{q}}\right)$. This operator should fulfil, for all values of n,

$$\begin{eqnarray}&&{{\mathcal{Q}}_{\tau}}|n\rangle =\underset{m}{\sum}\,|m\rangle {{Q}_{\tau}}(m,n).\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \left(\text{Q}\right)\end{eqnarray}$$

In analogy with the transition (rate) matrix ${{Q}_{\tau}}$, we call ${{\mathcal{Q}}_{\tau}}$ the transition (rate) operator (note that we only speak of 'operators' with respect to differential operators, but not with respect to matrices). Just like the basis (E)volution operator, ${{\mathcal{Q}}_{\tau}}\left(q,{{\partial}_{q}}\right)$ should be polynomial in ${{\partial}_{q}}$. In section 5, it will be assumed that ${{\mathcal{Q}}_{\tau}}\left(q,{{\partial}_{q}}\right)$ can be expanded in powers of both q and ${{\partial}_{q}}$. In a constructive approach, one could also define the operator as

$$\begin{eqnarray}&&{{\mathcal{Q}}_{\tau}}:=\underset{m,n}{\sum}\,|m\rangle {{Q}_{\tau}}(m,n)\langle n|.\end{eqnarray} \tag{ 53 }$$

This constructive definition does not guarantee, however, that ${{\mathcal{Q}}_{\tau}}\left(q,{{\partial}_{q}}\right)$ has the form of a differential operator. This property is, for example, not immediately clear for the Fourier basis function $|n{{\rangle}_{q}}={{\text{e}}^{\text{i}nq}}$ (with $n\in \mathbb{Z}$), which we complemented with the functional $\langle n|\,f={\int}_{-\pi}^{\pi}\frac{\text{d}q}{2\pi}{{\text{e}}^{-\text{i}nq}}f(q)$ in the introduction to this section. Most of the processes that we solve in later sections have polynomial transition rates. Suitable bases and operators for these processes are provided in the next section 2.2. It remains an open problem for the field to find such bases and operators for processes whose transition rates have different functional forms. That is, for example, the case for transition rates that saturate with the number of particles and have the form of a Hill function.

Provided that one has found a transition operator ${{\mathcal{Q}}_{\tau}}$ and a basis (E)volution operator ${{\mathcal{E}}_{\tau}}$ for a (C)omplete and (O)rthogonal basis, it follows from (49) that the generalized generating function $|g\left(\tau |{{t}_{0}},{{n}_{0}}\right)\rangle$ obeys the flow equation ⁶

$$\begin{eqnarray}&&{{\partial}_{\tau}}|g\rangle =\left({{\mathcal{E}}_{\tau}}+{{\mathcal{Q}}_{\tau}}\right)|g\rangle =:{{\tilde{\mathcal{Q}}}_{\tau}}|g\rangle.\end{eqnarray} \tag{ 54 }$$

Its initial condition reads $|g\left({{t}_{0}}|{{t}_{0}},{{n}_{0}}\right)\rangle =\,|{{n}_{0}}\rangle$, with the possibly time-dependent basis function $|{{n}_{0}}\rangle$ being evaluated at time t₀. Although time-dependent bases prove useful in section 7, we mostly work with time-independent bases in the following. The operators ${{\tilde{\mathcal{Q}}}_{\tau}}$ and ${{\mathcal{Q}}_{\tau}}$ then agree because ${{\mathcal{E}}_{\tau}}$ is zero. Therefore, we refer to both ${{\tilde{\mathcal{Q}}}_{\tau}}$ and ${{\mathcal{Q}}_{\tau}}$ as transition (rate) operators.

In our above derivation, we assumed that the ket $|n\rangle$ represents an ordinary function and that the bra $\langle n|$ represents a linear functional. To understand why we chose similar letters for the ${{Q}_{\tau}}$-matrix and the ${{\mathcal{Q}}_{\tau}}$-operator, it is insightful to consider the unit column vectors $|n\rangle :={{\widehat{\boldsymbol{e}}}_{n}}$ and the unit row vectors $\langle n|\,:=\widehat{\boldsymbol{e}}_{n}^{\intercal}$ (with $m,n\in {{\mathbb{N}}_{0}}$). For these vectors, the right hand side of the transition operator (53) simply constitutes a representation of the ${{Q}_{\tau}}$-matrix. Hence, ${{\mathcal{Q}}_{\tau}}$ is not a differential operator in this case but coincides with the ${{Q}_{\tau}}$-matrix. This observation does not come as a surprise because we already noted that the generating function (45) represents the vector $\boldsymbol{p}\left(\tau |{{t}_{0}},{{n}_{0}}\right)$ of all probabilities in this case. Moreover, the corresponding flow equation (54) does not constitute a linear PDE but a vector representation of the forward master equation (10).

Following Doi, the transition operator ${{\tilde{\mathcal{Q}}}_{\tau}}$ could be called a 'time evolution' or 'Liouville' operator [6, 7], or, following Zel'dovich and Ovchinnikov, a 'Hamiltonian' [8]. The latter name is due to the formal resemblance of the flow equation (54) to the Schrödinger equation in quantum mechanics [92] (this resemblance holds for any linear PDE that is first order in ${{\partial}_{\tau}}$). The name 'Hamiltonian' has gained in popularity throughout the recent years, possibly because the path integrals derived later in this review share many formal similarities with the path integrals employed in quantum mechanics [107, 285–287]. Nevertheless, let us point out that ${{\tilde{\mathcal{Q}}}_{\tau}}$ is generally not Hermitian and that the generating function $|g\rangle$ does not represent a wave function and also not a probability (unless one chooses the unit vectors $|n\rangle ={{\widehat{\boldsymbol{e}}}_{n}}$ and $\langle n|\,=\widehat{\boldsymbol{e}}_{n}^{\intercal}$ as basis). In this review, we stick to the name transition (rate) operator for ${{\tilde{\mathcal{Q}}}_{\tau}}$ (and ${{\mathcal{Q}}_{\tau}}$) to emphasize its connection to the transition (rate) matrix ${{Q}_{\tau}}$.

In the previous section, we formulated four conditions under which the master equation (10) can be cast into the linear PDE (54) obeyed by the generalized generating function (45). In the following, we illustrate how these conditions can be met for various stochastic processes.

2.2.1. Random walks.

A simple process that can be solved by the method from the previous section is the one-dimensional random walk. We model this process in terms of a particle sitting at position n of the one-dimensional lattice $\mathbb{L}:=\mathbb{Z}$. The particle may jump to the nearest lattice site on its right with the (possibly time-dependent) rate ${{r}_{\tau}}>0$, and to the site on its left with the rate ${{l}_{\tau}}>0$. Given that the particle was at position n₀ at time t₀, the probability of finding it at position n at time τ obeys the master equation

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{\tau}}p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right) & ={{r}_{\tau}}\left[p(\tau,n-1|\centerdot )-p(\tau,n|\centerdot )\right] \\ {} & \,\,\,\,\,+\,\,{{l}_{\tau}}\left[p(\tau,n+1|\centerdot )-p(\tau,n|\centerdot )\right], \end{array}\end{eqnarray} \tag{ 55 }$$

with initial condition $p\left({{t}_{0}},n|{{t}_{0}},{{n}_{0}}\right)={{\delta}_{n,{{n}_{0}}}}$. One can solve this equation by solving the associated flow equation (54). But for this purpose, we first require an (O)rthogonal and (C)omplete basis, as well as a basis (E)volution operator ${{\mathcal{E}}_{\tau}}$ and a transition operator ${{\mathcal{Q}}_{\tau}}$ (condition (Q)).

An appropriate choice of the orthogonal and complete basis proves to be the time-independent Fourier basis

$$\begin{eqnarray}&&|n{{\rangle}_{q}}:={{\text{e}}^{\text{i}nq}}~\text{and}~\langle n|\,f:={\int}_{-\pi}^{\pi}\frac{\text{d}q}{2\pi}{{\text{e}}^{-\text{i}nq}}f(q)\end{eqnarray} \tag{ 56 }$$

with $n\in \mathbb{Z}$ and test function f. For this basis, the generalized generating function $|g\rangle ={\sum}_{n}{{\text{e}}^{\text{i}nq}}p(\tau,n|\centerdot )=\langle {{\text{e}}^{\text{i}nq}}\rangle$ coincides with the characteristic function. Moreover, the corresponding orthogonality condition $\langle m|n\rangle ={{\delta}_{m,n}}$ agrees with a common integral representation of the Kronecker delta. The completeness of the basis can be shown with the help of a Fourier series representation of the 'Dirac comb'

$$\begin{eqnarray}&&\text{III}(q)=\underset{n\in \mathbb{Z}}{\sum}\,\delta (q-2\pi n)=\frac{1}{2\pi}\underset{n\in \mathbb{Z}}{\sum}\,{{\text{e}}^{\text{i}nq}}.\end{eqnarray} \tag{ 57 }$$

It thereby follows that

$$\begin{eqnarray}&&\underset{n\in \mathbb{Z}}{\sum}\,|n{{\rangle}_{q}}\langle n|\,f={\int}_{-\pi}^{\pi}\text{d}{{q}^{\prime}}\,\text{III}\left(q-{{q}^{\prime}}\right)f\left({{q}^{\prime}}\right)=f(q).\end{eqnarray} \tag{ 58 }$$

Since the Fourier basis function is time-independent, the condition (E) is trivially fulfilled for the evolution operator ${{\mathcal{E}}_{\tau}}:=0$. The only piece still missing is the transition operator ${{\mathcal{Q}}_{\tau}}$. Its defining condition (Q) requires knowledge of the transition matrix ${{Q}_{\tau}}$ whose elements

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{Q}_{\tau}}(m,n) \\ {} & ={{r}_{\tau}}\left({{\delta}_{m,n+1}}-{{\delta}_{m,n}}\right)+{{l}_{\tau}}\left({{\delta}_{m,n-1}}-{{\delta}_{m,n}}\right) \end{array}\end{eqnarray} \tag{ 59 }$$

can be inferred by comparing the master equation (55) with its general form (10). The condition (Q) therefore reads

$$\begin{eqnarray}&&{{\mathcal{Q}}_{\tau}}|n\rangle ={{r}_{\tau}}(|n+1\rangle \,-|n\rangle )+{{l}_{\tau}}(|n-1\rangle \,-|n\rangle ).\end{eqnarray} \tag{ 60 }$$

One can construct a transition operator with this property with the help of the functions $c(q):={{\text{e}}^{\text{i}q}}$ and $a(q):={{\text{e}}^{-\text{i}q}}$. The function c shifts the basis function $|n{{\rangle}_{q}}={{\text{e}}^{\text{i}nq}}$ to the right via $c|n\rangle =\,|n+1\rangle$, and the function a shifts it to the left via $a|n\rangle =\,|n-1\rangle$. Thus, an operator with the property (60) can be defined as

$$\begin{eqnarray}&&{{\mathcal{Q}}_{\tau}}\left(q,{{\partial}_{q}}\right):={{r}_{\tau}}\left[c(q)-1\right]+{{l}_{\tau}}\left[a(q)-1\right].\end{eqnarray} \tag{ 61 }$$

This operator can also be inferred from its constructive definition (53) by making use of the Dirac comb.

After putting the above pieces together, one finds that the generating function $|g\left(\tau |{{t}_{0}},{{n}_{0}}\right){{\rangle}_{q}}$ obeys the flow equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}|g\rangle =\left[{{r}_{\tau}}\left({{\text{e}}^{\text{i}q}}-1\right)+{{l}_{\tau}}\left({{\text{e}}^{-\text{i}q}}-1\right)\right]|g\rangle \end{eqnarray} \tag{ 62 }$$

for ${{t}_{0}}\leqslant \tau \leqslant t$ with the initial value $|g\left({{t}_{0}}|{{t}_{0}},{{n}_{0}}\right)\rangle =|{{n}_{0}}\rangle ={{\text{e}}^{\text{i}{{n}_{0}}q}}$. The flow equation is readily solved for $|g\rangle$ by

$$\begin{eqnarray}&&\exp \left(\left({{\text{e}}^{\text{i}q}}-1\right){\int}_{{{t}_{0}}}^{\tau}\text{d}s\,{{r}_{s}}+\left({{\text{e}}^{-\text{i}q}}-1\right){\int}_{{{t}_{0}}}^{\tau}\text{d}s\,{{l}_{s}}\right)|{{n}_{0}}\rangle.\end{eqnarray}$$

The conditional probability distribution can be recovered from this generating function by performing the inverse Fourier transformation (48), i.e. by evaluating $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)=\langle n|g\rangle$. A series expansion of all of the involved exponentials and some laborious rearrangement of sums eventually result in a Skellam distribution as the solution of the process [329]. We outline the derivation of this distribution in appendix C. The distribution's mean is $\mu ={{n}_{0}}+{\int}_{{{t}_{0}}}^{\tau}\text{d}s\,\left({{r}_{s}}-{{l}_{s}}\right)$ and its variance ${{\sigma}^{2}}={\int}_{{{t}_{0}}}^{\tau}\text{d}s\,\left({{r}_{s}}+{{l}_{s}}\right)$. These moments also follow from the fact that the Skellam distribution describes the difference of two Poisson random variables; one for jumps to the right, the other for jumps to the left. The two moments can, also be obtained more easily by deriving equations for their time evolution from the master equation. Those equations do not couple for the simple random walk. Let us also note that if the two jump rates ${{r}_{\tau}}$ and ${{l}_{\tau}}$ agree, the Skellam distribution tends to a Gaussian for large times.

2.2.2. Chemical reactions.

As a second example, we turn to processes whose transition rates can be decomposed additively into the transition rates of simple chemical reactions. In such a reaction, k₁, k₂,... molecules of types A₁, A₂,... come together to be replaced by l₁, l₂,... molecules of the same types (with ${{k}_{j}},{{l}_{j}}\in {{\mathbb{N}}_{0}}$). Besides reacting with each other, the molecules could also diffuse in space, which can be modelled in terms of hopping processes on a regular, d-dimensional lattice such as $\mathbb{L}:={{\mathbb{Z}}^{d}}$. Upon labelling particles on different lattice sites by their positions, the hopping of a molecule of type A₁ from lattice site $i\in \mathbb{L}$ to lattice site $j\in \mathbb{L}$ could be regarded as the chemical reaction $A_{1}^{(i)}\to A_{1}^{(\,j)}$. We consider the 'chemical' master equation associated to such hopping processes in section 4.5.2.

To demonstrate the generating function approach from section 2.1, we focus on a system with only a single type of molecule A engaging in the 'well-mixed' reaction $k\,A\to l\,A$ ($k,l\in {{\mathbb{N}}_{0}}$). Since the forward master equation (10) is linear in the transition rate ${{Q}_{\tau}}(n,m)$, the following considerations readily extend to networks of multiple reactions, multiple types of molecules, and processes with spatial degrees of freedom. The basis functions and functionals introduced below can, for example, be used to study branching and annihilating random walks, which are commonly modelled in terms of diffusing particles that engage in the binary annihilation reaction $2\,A\to \varnothing$ and the linear growth process $A\to (1+m)A$ [96, 132, 317, 330, 331]. For an odd number of offspring, these walks exhibit an absorbing state phase transition falling into the universality class of directed percolation (according to perturbative calculations in one and two spatial dimensions [96, 132], according to non-perturbative calculations in up to six dimensions [318]). The decomposition of a process into elementary reactions of the form $k\,A\to l\,A$ is also possible in the contexts of growing polymer crystals [13], aggregation phenomena [294], and predator-prey ecosystems [163].

As explained in section 1.5, the transition rate ${{w}_{\tau}}(m,n)={{\gamma}_{\tau}}{{\delta}_{m,n-k+l}}{{(n)}_{k}}$ of the reaction $k\,A\to l\,A$ is determined by combinatorial counting. Its proportionality to the falling factorial ${{(n)}_{k}}=n(n-1)\cdots (n-k+1)$ derives from picking k molecules out of a population of n molecules. The overall time scale of the reaction is set by the rate coefficient ${{\gamma}_{\tau}}$ (this coefficient may also absorb a factorial $k!$ accounting for the indistinguishability of molecules). The above transition rate guarantees that the number of molecules (or 'particles') in the system never drops below zero, provided that it was non-negative initially. Thus, the state space of n is ${{\mathbb{N}}_{0}}$ and we require basis functions $|n\rangle$ and basis functionals $\langle n|$ only for such values.

Before specifying the (C)omplete and (O)rthogonal basis as well as the basis (E)volution operator ${{\mathcal{E}}_{\tau}}$, let us first specify an appropriate transition operator ${{\mathcal{Q}}_{\tau}}$. Its corresponding condition (Q) depends on the transition matrix ${{Q}_{\tau}}$, whose elements ${{Q}_{\tau}}(m,n)={{\gamma}_{\tau}}{{(n)}_{k}}\left({{\delta}_{m,n-k+l}}-{{\delta}_{m,n}}\right)$ can be inferred from the rate ${{w}_{\tau}}(m,n)$ with the help of the relation (9). Consequently, the condition (Q) reads

$$\begin{eqnarray}&&{{\mathcal{Q}}_{\tau}}|n\rangle ={{\gamma}_{\tau}}{{(n)}_{k}}(|n-k+l\rangle -|n\rangle ).\end{eqnarray} \tag{ 63 }$$

This condition is met by the transition operator

$$\begin{eqnarray}&&{{\mathcal{Q}}_{\tau}}(c,a):={{\gamma}_{\tau}}\left({{c}^{l}}-{{c}^{k}}\right){{a}^{k}},\end{eqnarray} \tag{ 64 }$$

provided that there exist a 'creation' operator c and an 'annihilation' operator a acting as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} c|n\rangle =\,|n+1\rangle \,\text{and} \end{array}\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} a|n\rangle =n|n-1\rangle. \end{array}\end{eqnarray} \tag{ 66 }$$

The second of these relations ensures that basis functions with n  <  0 do not appear because $a|0\rangle$ vanishes. Both c and a may depend on time, just as the basis functionals and functions $\langle n|$ and $|n\rangle$ may do. It follows from (65) and (66) that the two operators fulfil the commutation relation

$$\begin{eqnarray}&&\left[{{a}_{\tau}},{{c}_{\tau}}\right]:={{a}_{\tau}}{{c}_{\tau}}-{{c}_{\tau}}{{a}_{\tau}}=1\end{eqnarray} \tag{ 67 }$$

at a fixed time τ. The commutation relation is meant with respect to functions that can be expanded in the basis functions yet to be defined. For brevity, we commonly drop the subscript τ of the creation and annihilation operators. Together with the orthogonality condition, the commutation relation implies that the operators act on the basis functionals as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle n|c=\langle n-1|~\text{and} \end{array}\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle n|a=(n+1)\langle n+1|. \end{array}\end{eqnarray} \tag{ 69 }$$

In quantum mechanics, creation and annihilation operators prove useful in solving the equation of motion of a quantum particle in a quadratic potential (i.e. in solving the Schrödinger equation of the quantum harmonic oscillator) [22, 332]. There, the ket $|n\rangle$ is interpreted as carrying n quanta of energy $\hbar \omega$ in addition to the ground state energy $\hbar \omega /2$ of the oscillator's 'vacuum state' $|0\rangle$ (with Dirac constant $\hbar$ and angular frequency ω). The creation operator then adds a quantum of energy to state $|n\rangle$ via the relation $c|n\rangle =\sqrt{n+1}|n+1\rangle$, and the annihilation operator removes an energy quantum via $a|n\rangle =\sqrt{n}|n-1\rangle$ (a also destroys the vacuum state). These relations differ from the ones in (65) and (66) because in quantum mechanics, the creation and annihilation operators are defined in terms of self-adjoint position and momentum operators, forcing the former operators to be hermitian adjoints of each other (i.e. $c={{a}^{\dagger}}$ and $a={{c}^{\dagger}}$). Consequently, the relations corresponding to (68) and (69) read $\langle n|{{a}^{\dagger}}=\sqrt{n}\langle n-1|$ and $\langle n|a=\sqrt{n+1}\langle n+1|$. Nevertheless, the commutation relation (67) also holds in the quantum world in which its validity ultimately derives from the non-commutativity of the position operator Q and the momentum operator P (i.e. from $[Q,P]=\text{i}\hbar$, which follows from P being the generator of spatial displacements in state space [332]). Besides, let us note that in quantum field theory, the creation operator is interpreted as adding a bosonic particle to an energy state, and the annihilation operator as removing one [333].

One may wonder whether the above interpretations can be transferred to the theory of stochastic processes in which the creation and annihilation operators need not be each other's adjoints. For example, the creation operator in (65) could be interpreted as adding a molecule to a system, and the corresponding annihilation operator (66) as removing a molecule. The commutation relation (67) may then be interpreted as that the addition of a particle to the system (one way to do it) and the removal of a particle (many ways to do it) do not commute. Let us, however, note that these interpretations only apply to the processes discussed in the present section, which can be decomposed into reactions of the form $k\,A\to l\,A$ (with possibly multiple types of particles and spatial degrees of freedom). The interpretations do not apply, for example, to the random walk of a single particle with state space $\mathbb{Z}$, which we discussed in the previous section (there, the 'shift operators' with actions $c|n\rangle =\,|n+1\rangle$ and $a|n\rangle =\,|n-1\rangle$ commute).

Whether creation and annihilation operators with the properties (65) and (66) exist depends on the choice of the basis functions. For the study of chemical reactions, a useful choice proves to be the basis function

$$\begin{eqnarray}&&|n{{\rangle}_{q}}:={{\left(\zeta q+\tilde{q}\right)}^{n}}{{\text{e}}^{-\tilde{x}\left(\zeta q+\tilde{q}\right)}}.\end{eqnarray} \tag{ 70 }$$

Here, $\tilde{q}(\tau )$ and $\tilde{x}(\tau )$ are arbitrary, possibly time-dependent functions and $\zeta \ne 0$ is a free parameter. This parameter only becomes relevant in section 6 in recovering a path integral representation of averaged observables. There, its value is set to $\zeta :=\text{i}$ but typically we choose $\zeta :=1$. For the latter choice, the basis function (70) simplifies to

$$\begin{eqnarray}&&|n{{\rangle}_{q}}={{\left(q+\tilde{q}\right)}^{n}}{{\text{e}}^{-\tilde{x}\left(q+\tilde{q}\right)}}.\end{eqnarray} \tag{ 71 }$$

Alternatively, the parameter ζ could be used to rescale the variable q by a system size parameter N if such a parameter is available.

The two functions $\tilde{q}$ and $\tilde{x}$ may prove helpful in simplifying the flow equation obeyed by the generating function. For example, we chose $\tilde{q}:=1$ and $\tilde{x}:=\gamma /\mu$ in the introduction to section 2 to simplify the flow equation of the reaction $\varnothing \rightleftharpoons A$ (with growth rate coefficient γ and decay rate coefficient μ). Later, in section 7, $\tilde{q}$ and $\tilde{x}$ will act as the stationary paths of a path integral with q and an auxiliary variable x being deviations from them. If both $\tilde{q}$ and $\tilde{x}$ are chosen as zero, the basis function (71) simplifies to the basis function

$$\begin{eqnarray}&&|n{{\rangle}_{q}}={{q}^{n}}\end{eqnarray} \tag{ 72 }$$

of the ordinary probability generating function.

For the general basis function (70), creation and annihilation operators with the properties (65) and (66) can be defined as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} c\left(q,{{\partial}_{q}}\right):=\zeta q+\tilde{q}~\text{and} \end{array}\end{eqnarray} \tag{ 73 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} a\left(q,{{\partial}_{q}}\right):={{\partial}_{\zeta q}}+\tilde{x}. \end{array}\end{eqnarray} \tag{ 74 }$$

Apparently, these operators are not each other's adjoints. For the basis function $|n{{\rangle}_{q}}={{q}^{n}}$ of the ordinary probability generating function, the operators simplify to c  =  q and $a={{\partial}_{q}}$. The corresponding transition operator (64) reads ${{\mathcal{Q}}_{\tau}}\left(q,{{\partial}_{q}}\right)={{\gamma}_{\tau}}\left({{q}^{l}}-{{q}^{k}}\right)\partial _{q}^{k}$, resulting in the flow equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}|g\rangle ={{\gamma}_{\tau}}\left({{q}^{l}}-{{q}^{k}}\right)\partial _{q}^{k}|g\rangle.\end{eqnarray} \tag{ 75 }$$

Thus far, we have not specified the basis functionals. These functionals can be defined using the annihilation operator (74). In particular, we complement the basis function (70) with the functionals

$$\begin{eqnarray}&&\langle n|\,f:=\frac{{{a}^{n}}}{n!}f(q){{|}_{q=-\tilde{{q}}\,/\zeta}},\end{eqnarray} \tag{ 76 }$$

where f represents a test function. The (O)rthogonality of the basis ${{\left\{\langle m|,|n\rangle \right\}}_{m,n\in {{\mathbb{N}}_{0}}}}$ follows directly from the action of the annihilation operator on the basis function $|n\rangle$, and from the vanishing of this function at $q=-\tilde{q}/\zeta$, except for n  =  0. The (C)ompleteness of the basis is also readily established ⁷ .

The only piece still missing is the basis (E)volution operator ${{\mathcal{E}}_{\tau}}$ to encode the time-dependence of the basis function (70). Differentiation of this function with respect to time and use of the creation and annihilation operators (73) and (74) show that this operator is given by

$$\begin{eqnarray}&&{{\mathcal{E}}_{\tau}}(c,a):=\left({{\partial}_{\tau}}\tilde{q}\right)\left(a-\tilde{x}\right)-\left({{\partial}_{\tau}}\tilde{x}\right)c.\end{eqnarray} \tag{ 77 }$$

Let us illustrate the use of the above basis for a process whose transition operator can be reduced to a mere constant by an appropriate choice of $\tilde{q}$ and $\tilde{x}$. This simplification is, however, bought by making the basis functions time-dependent. In particular, we consider the simple growth process $\varnothing \to A$ for a time-independent growth rate coefficient γ. By our previous discussions, one can readily verify that the ordinary probability generating function with basis function $|n{{\rangle}_{q}}={{q}^{n}}$ obeys the flow equation ${{\partial}_{\tau}}|g\rangle =\gamma (q-1)|g\rangle$ for this process. This flow equation can be simplified by redefining the basis function of the generating function as $|n{{\rangle}_{\tau,q}}:={{(q+1)}^{n}}{{\text{e}}^{-\tilde{x}(q+1)}}$, with $\tilde{x}$ solving the rate equation ${{\partial}_{\tau}}\tilde{x}=\gamma$ of the process. Hence, the basis function depends explicitly on time through its dependence on $\tilde{x}(\tau )=\tilde{x}\left({{t}_{0}}\right)+\gamma \left(\tau -{{t}_{0}}\right)$. Upon combining the transition operator (64) with the basis evolution operator (77), one finds that the generating function now obeys the flow equation ${{\partial}_{\tau}}|g\rangle =-\gamma |g\rangle$. The equation is readily solved by $|g\left(\tau |{{t}_{0}},{{n}_{0}}\right)\rangle ={{\text{e}}^{-\gamma \left(\tau -{{t}_{0}}\right)}}|{{n}_{0}}{{\rangle}_{{{t}_{0}}}}$. The conditional probability distribution can be recovered via the inverse transformation (48) as

$$\begin{eqnarray}&&p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)={{\text{e}}^{-\gamma \left(\tau -{{t}_{0}}\right)}}{{\langle n{{|}_{\tau}}|{{n}_{0}}\rangle}_{{{t}_{0}}}}.\end{eqnarray} \tag{ 78 }$$

The coefficient ${{\langle n{{|}_{\tau}}|{{n}_{0}}\rangle}_{{{t}_{0}}}}$ can be evaluated by determining how the functional $\langle n{{|}_{\tau}}$ acts at time t₀. Using ${{a}_{\tau}}-\tilde{x}(\tau )={{a}_{{{t}_{0}}}}-\tilde{x}\left({{t}_{0}}\right)$ and the binomial theorem, one can rewrite the action of this functional as

$$\begin{eqnarray}&&\langle n{{|}_{\tau}}f=\underset{m=0}{\overset{n}{\sum}}\,\frac{{{\left(\gamma \left(\tau -{{t}_{0}}\right)\right)}^{n-m}}}{(n-m)!}\langle m{{|}_{{{t}_{0}}\,}}f.\end{eqnarray} \tag{ 79 }$$

By the (O)rthogonality condition, the conditional probability distribution (78) therefore evaluates to the shifted Poisson distribution

$$\begin{eqnarray}&&p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)=\frac{{{\text{e}}^{-\gamma \left(\tau -{{t}_{0}}\right)}}{{\left(\gamma \left(\tau -{{t}_{0}}\right)\right)}^{n-{{n}_{0}}}}}{\left(n-{{n}_{0}}\right)!}{{ \Theta }_{n-{{n}_{0}}}},\end{eqnarray} \tag{ 80 }$$

where ${{ \Theta }_{n}}:=1$ for $n\geqslant 0$ and ${{ \Theta }_{n}}:=0$ otherwise. The validity of this solution can be verified by solving the corresponding flow equation of the ordinary generating function.

In our above approach, we first established the form of the transition operator because the basis function (70) is not the only possible choice. Ohkubo recently proposed the use of orthogonal polynomials for this purpose [334]. In analogy to the eigenfunctions of the quantum harmonic oscillator [22, 332], one can, for example, choose the basis function as the Hermite polynomial

$$\begin{eqnarray}&&|n{{\rangle}_{q}}:=H{{e}_{n}}(q)={{(-1)}^{n}}{{\text{e}}^{{{q}^{2}}/2}}\partial _{q}^{n}\,{{\text{e}}^{-{{q}^{2}}/2}},\end{eqnarray} \tag{ 81 }$$

with $n\in {{\mathbb{N}}_{0}}$. Hermite polynomials constitute an Appell sequence, i.e. they fulfil ${{\partial}_{q}}H{{e}_{n}}(q)=nH{{e}_{n-1}}(q)$ (18.9.27 in [335]). With $a:={{\partial}_{q}}$, this property coincides with the defining relation $a|n\rangle =n|n-1\rangle$ of the annihilation operator in (66). Furthermore, Hermite polynomials obey the recurrence relation $H{{e}_{n+1}}(q)=qH{{e}_{n}}(q)-nH{{e}_{n-1}}(q)$ (18.9.1 in [335]). Combined with the Appell property, $c:=q-{{\partial}_{q}}$ therefore fulfils the defining relation $c|n\rangle =\,|n+1\rangle$ of a creation operator in (65). After complementing the basis function (81) with the functional

$$\begin{eqnarray}&&\langle m|\,f:=\frac{1}{m!}{\int}_{-\infty}^{\infty}\text{d}q\,\frac{{{\text{e}}^{-{{q}^{2}}/2}}}{\sqrt{2\pi}}H{{e}_{m}}(q)f(q),\end{eqnarray} \tag{ 82 }$$

the (O)rthogonality and (C)ompleteness of the basis are also established (18.3 and 18.18.6 in [335]). Thus, the chemical master equation (24) can be transformed into a flow equation obeyed by a generating function based on Hermite polynomials. To our knowledge, however, no stochastic process has thus far been solved or been approximated along these lines. Besides Hermite polynomials, Ohkubo also proposed the use of Charlier polynomials and mentioned their relation with certain birth-death processes [334].

2.2.3. Intermezzo: the unit vector basis.

One may wonder why we actually bother with explicit representations of the bras $\langle n|$ and kets $|n\rangle$. Often, these objects are introduced only formally as the basis of a 'bosonic Fock space' [6, 7, 91, 289, 336], leaving the impression that the particles under consideration are in fact bosonic quantum particles. Although this impression takes the analogy with quantum theory too far, the analogy has helped in developing new methods for solving the master equations. In the previous sections, we showed how the forward master equation can be cast into a linear PDE obeyed by a probability generating function. Later, in section 5, this PDE is solved in terms of a path integral by which we then recover a path integral representation of averaged observables in section 6. This 'analytic' derivation of the path integral is, however, not the only possible way. In the following, we sketch the mathematical basis of an alternative approach, which employs the unit vectors ${{\widehat{\boldsymbol{e}}}_{n}}$ as the basis. The approach ultimately results in the same path integral representation of averaged observables as we show in section 6.2. The duality between the two approaches resembles the duality between the matrix mechanics formulation of quantum mechanics by Heisenberg, Born, and Jordan [93–95] and its analytic formulation in terms of the Schrödinger equation [92].

The alternative derivation of the path integral also starts out from the forward master equation ${{\partial}_{\tau}}p\left(\tau |{{t}_{0}}\right)=Qp\left(\tau |{{t}_{0}}\right)$. As before, $p\left(\tau |{{t}_{0}}\right)$ represents the matrix of the conditional probabilities $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$, and Q is the transition rate matrix. For simplicity, we assume the transition rate matrix to be independent of time. Moreover, we focus on processes that can be decomposed additively into chemical reactions of the form $k\,A\to l\,A$ with time-independent rate coefficients. The elements of the transition matrix associated to this reaction read $Q(m,n)=\gamma {{(n)}_{k}}\left({{\delta}_{m,n-k+l}}-{{\delta}_{m,n}}\right)$. Since the particle numbers n and n₀ may assume any values in ${{\mathbb{N}}_{0}}$, the probability matrix $p\left(\tau |{{t}_{0}}\right)$ and the transition matrix Q have infinitely many rows and columns.

In the introductory section 1.4, we formulated conditions under which the forward master equation is solved by the matrix exponential $p\left(\tau |{{t}_{0}}\right)={{\text{e}}^{Q\left(\tau -{{t}_{0}}\right)}}{\mathbb {1}}$ (in the sense of the expansion (16), the state space must countable and the exit rates bounded). For the above chemical reaction, those conditions are not necessarily met, and thus the matrix exponential may not exist. Nevertheless, we regard $p\left(\tau |{{t}_{0}}\right)={{\text{e}}^{Q\left(\tau -{{t}_{0}}\right)}}{\mathbb {1}}$ as a 'formal' solution in the following. With the help of certain mathematical 'tricks', this solution will be cast into a path integral in section 6. But before, let us rewrite the transition matrix in the exponential in terms of creation and annihilation matrices. For this purpose, we employ the (infinitely large) unit column vectors $|n\rangle :={{\widehat{\boldsymbol{e}}}_{n}}$ and their orthogonal unit row vectors $\langle n|\,\,:\,=\widehat{\boldsymbol{e}}_{n}^{\intercal}$ as basis (with $n\in {{\mathbb{N}}_{0}}$). Individual probabilities can therefore be inferred from the above solution as

$$\begin{eqnarray}&&p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)=\langle n|{{\text{e}}^{Q\left(\tau -{{t}_{0}}\right)}}|{{n}_{0}}\rangle.\end{eqnarray} \tag{ 83 }$$

For multivariate or spatial processes, the basis vectors can be generalized to tensors (i.e. $|\boldsymbol{n}\rangle =\,|{{n}_{1}}\rangle \otimes |{{n}_{2}}\rangle \otimes \ldots$).

The orthogonality of the basis allows us to rewrite the elements of the transition rate matrix of the reaction $k\,A\to l\,A$ as

$$\begin{eqnarray}&&Q(m,n)=\langle m|\left(\gamma \,{{(n)}_{k}}(|n-k+l\rangle -|n\rangle )\right).\end{eqnarray} \tag{ 84 }$$

A matrix with these elements can be written as

$$\begin{eqnarray}&&Q=\gamma \left({{c}^{l}}-{{c}^{k}}\right){{a}^{k}},\end{eqnarray} \tag{ 85 }$$

with c acting as a 'creation matrix' and a acting as an 'annihilation matrix'. In particular, we define c in terms of its sub-diagonal $(1,1,1,\ldots )$ and a in terms of its super-diagonal $(1,2,3,\ldots )$. All of their other matrix elements are set to zero. It is readily shown that these matrices fulfil $c|n\rangle =\,|n+1\rangle$ and $a|n\rangle =n|n-1\rangle$ with respect to the basis vectors. These relations coincide with our previous relations (65) and (66). Besides, the matrices also fulfil $\langle n|c=\langle n-1|$ and $\langle n|a=(n+1)\langle n+1|$, as well as the commutation relation $[a,c]={\mathbb {1}}$. Further properties of the matrices are addressed in section 6.2. By our previous comments in section 2.1, it is not a surprise that the transition matrix (85) has the same form as the transition operator (64). Both of them are normal-ordered polynomials in c and a. This property will help us in section 6.2 in recovering a path integral representation of averaged observables.

2.2.4. Processes with locally excluding particles.

As noted in the introductory section 1.3, one can model the movement of molecular motors along a cytoskeletal filament in terms of a master equation [113–115]. In its simplest form, such a model describes the movement of mutually excluding motors as a hopping process on a one-dimensional lattice $\mathbb{L}\subset \mathbb{Z}$, with attachment and detachment of motors at certain boundaries. To respect the mutual exclusion of motors, their local number n_i on a lattice site $i\in \mathbb{L}$ is restricted to 0 and 1. Various other processes can be modelled in similar ways, for example, aggregation processes [26], adsorption processes [337, 338], and directed percolation [29]. In the following, we show how the generating function approach from section 2.1 and its complementary approach from section 2.2.3 can be applied to such processes. To illustrate the mathematics behind these approaches while not burdening ourselves with too many indices, we demonstrate the approaches for the simple, non-spatial telegraph process [339, 340]. This process describes a system that randomly switches between two states that are called the 'on' and 'off' states, or, for brevity, the '1' and '0' states. The rate at which the system switches from state 1 to state 0 is denoted as μ, and the rate of the reverse transition as γ. For simplicity, we assume these rate coefficients to be time-independent. Consequently, the master equation of the telegraph process reads

$$\begin{eqnarray}{{\partial}_{\tau}}\boldsymbol{p}=\left(\begin{array}{c} -\mu & \gamma \\ \mu & -\gamma \end{array}\right)\,\boldsymbol{p},\end{eqnarray} \tag{ 86 }$$

with $\boldsymbol{p}(\tau |\centerdot )={{\left(\,p(\tau,1|\centerdot ),p(\tau,0|\centerdot )\right)}^{\intercal}}$ being the probability vector. Alternatively, the master equation can be written as ${{\partial}_{\tau}}p(\tau,n|\centerdot )={\sum}_{m\in \left\{1,0\right\}}Q(n,m)p(\tau,m|\centerdot )$ with transition matrix elements

$$\begin{eqnarray}&&Q(n,m)=\mu \left({{\delta}_{n,0}}-{{\delta}_{n,1}}\right){{\delta}_{m,1}}+\gamma \left({{\delta}_{n,1}}-{{\delta}_{n,0}}\right){{\delta}_{m,0}}.\,\end{eqnarray} \tag{ 87 }$$

One can solve the above master equation in various ways, for example, by evaluating the matrix exponential in $p\left(\tau |{{\tau}_{0}}\right)={{\text{e}}^{Q\left(\tau -{{\tau}_{0}}\right)}}{\mathbb {1}}$ after diagonalizing the transition matrix. In the following, the telegraph process serves as the simplest representative of processes with excluding particles and it helps in explaining how the methods from the previous sections are applied to such processes. The inclusion of spatial degrees of freedom into the procedure is straightforward.

To apply the generating function technique from section 2.1, we require orthogonal and complete basis functions $\left\{|0\rangle,|1\rangle \right\}$ and basis functionals $\left\{\langle 0|,\langle 1|\right\}$, as well as a transition operator $\mathcal{Q}\left(q,{{\partial}_{q}}\right)$ fulfilling condition (Q), i.e. for $n\in \left\{0,1\right\}$,

$$\begin{eqnarray}&&\mathcal{Q}|n\rangle =\mu (|0\rangle -|1\rangle ){{\delta}_{n,1}}+\gamma (|1\rangle -|0\rangle ){{\delta}_{n,0}}.\end{eqnarray} \tag{ 88 }$$

An operator with this property can be constructed in (at least) two ways: using creation and annihilation operators fulfilling the 'bosonic' commutation relation [a,c]  =  ac  −  ca  =  1, or using operators fulfilling the 'fermionic' anti-commutation relation $\left\{a,c\right\}:=ac+ca=1$.

Let us first outline an approach based on the commutation relation. The approach was recently proposed by van Wijland [90]. To illustrate the approach, we take the 'analytic' point of view with $c\left(q,{{\partial}_{q}}\right)$ and $a\left(q,{{\partial}_{q}}\right)$ being differential operators. However, the following considerations also apply to the approach from the previous section where c and a represented creation and annihilation matrices. To cast a master equation such as (86) into a path integral, van Wijland employed creation and annihilation operators with the same actions as in section 2.2.2, i.e. $c|n\rangle =\,|n+1\rangle$ and $a|n\rangle =n|n-1\rangle$. Besides complying with the commutation relation [a,c]  =  1, these relations imply that the 'number operator' $\mathcal{N}:\,=ca$ fulfils $\mathcal{N}\,|n\rangle =n|n\rangle$. This operator may be used to define the 'Kronecker operator'

$$\begin{eqnarray}&&{{\delta}_{\mathcal{N},m}}:={\int}_{-\pi}^{\pi}\frac{\text{d}u}{2\pi}{{\text{e}}^{\text{i}u\left(\mathcal{N}-m\right)}}\end{eqnarray} \tag{ 89 }$$

in terms of a series expansion of its exponential. The Kronecker operator acts on the basis functions as ${{\delta}_{\mathcal{N},m}}|n\rangle ={{\delta}_{n,m}}|n\rangle$, and thus it can be used to replace the Kronecker deltas in (88). One thereby arrives at the flow equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}|g\rangle =\mathcal{Q}|g\rangle =\left[\mu (a-1){{\delta}_{\mathcal{N},1}}+\gamma (c-1){{\delta}_{\mathcal{N},0}}\right]|g\rangle.\end{eqnarray} \tag{ 90 }$$

On the downside, the flow equation now involves power series with arbitrarily high derivatives with respect to q (upon choosing the operators as c  =  q and $a={{\partial}_{q}}$). In section 5, we show how the solution of a flow equation such as (90) can be expressed in terms of a path integral. The derivation of the path integral requires that the transition operator $\mathcal{Q}\left(q,{{\partial}_{q}}\right)$ is normal-ordered with respect to q and ${{\partial}_{q}}$, i.e. that all the q are to the left of all the ${{\partial}_{q}}$ in every summand. This order can be achieved by the repeated use of $\left[{{\partial}_{q}},q\right]\,f(q)=f(q)$. More information on the procedure can be found in [90]. Van Wijland applied the resulting path integral to the asymmetric diffusion of excluding particles on a one-dimensional lattice. Moreover, Mobilia, Georgiev, and Täuber have employed the method in studying the stochastic Lotka-Volterra model on d-dimensional lattices [153].

An alternative to the above approach lies in the use of operators fulfilling the anti-commutation relation {a, c}  =  ac  +  ca  =  1. This relation is clearly not fulfilled by $c:=q$ and $a:={{\partial}_{q}}$, at least not if q represents an ordinary real variable. However, ${{\partial}_{q}}q+q{{\partial}_{q}}=1$ holds true if q represents a Grassmann variable (see [341, 342] for details). Grassmann variables commute with real and complex numbers (i.e. $[q,\alpha ]=0$ for $\alpha \in \mathbb{C}$), but they anti-commute with themselves and with other Grassmann variables (i.e. $\left\{q,\tilde{q}\right\}=0$). Grassmann variables have, for example, proven useful in calculating the correlation functions of kinetic Ising models [343], even when the system is driven far from equilibrium [344].

Instead of using Grassmann variables for the basis, let us consider a representation of the basis defined by the unit row vectors $\langle 0|\,=(0,1)$ and $\langle 1|\,=(1,0)$, and by the unit column vectors $|0\rangle ={{(0,1)}^{\intercal}}$ and $|1\rangle ={{(1,0)}^{\intercal}}$. Obviously, these vectors fulfil the orthogonality condition $\langle m|n\rangle ={{\delta}_{m,n}}$, and ${\sum}_{n}|n\rangle \langle n|\,={\mathbb {1}}$ is the 2-by-2 unit matrix. The anti-commutation relation {a, c}  =  1 is met by the creation and annihilation matrices

$$\begin{eqnarray}c:={{\sigma}^{+}}=\left(\begin{array}{c} 0 & 1 \\ 0 & 0 \end{array}\right)\,\text{and}\,a:={{\sigma}^{-}}=\left(\begin{array}{c} 0 & 0 \\ 1 & 0 \end{array}\right).\end{eqnarray} \tag{ 91 }$$

The creation matrix acts as $c|0\rangle =\,|1\rangle$ and $c|1\rangle =\mathbf{0}$, and the annihilation matrix as $a|1\rangle =\,|0\rangle$ and $a|0\rangle =\mathbf{0}$ (zero vector). Consequently, the matrix product ca fulfils $ca|0\rangle =\mathbf{0}$ and $ca|1\rangle =|1\rangle$, and thus it serves the same purpose as the Kronecker delta ${{\delta}_{n,1}}$ in (88). Analogously, the operator ac serves the same purpose as ${{\delta}_{n,0}}$. The master equation (86) can therefore be written in a form resembling the flow equation (90), namely as

$$\begin{eqnarray}&&{{\partial}_{\tau}}\boldsymbol{p}=\left[\mu \left(a-{\mathbb {1}}\right)ca+\gamma \left(c-{\mathbb {1}}\right)ac\right]\boldsymbol{p}.\end{eqnarray} \tag{ 92 }$$

In this form, the stochastic process mimics a spin-1/2 problem (especially, if the creation and annihilation matrices are written in terms of the Pauli spin matrices ${{\sigma}_{x}}$, ${{\sigma}_{y}}$, and ${{\sigma}_{z}}$). More information on spin-representations of master equations is provided in [66, 337, 338, 345–348]. A spin-representation has, for example, been employed in the analysis of reaction–diffusion master equations via the density matrix renormalization group [349, 350]. In their study of directed percolation of excluding particles on the one-dimensional lattice $\mathbb{Z}$, Brunel, Oerding, and van Wijland performed a Jordan–Wigner transformation of the spatially extended 'spin' matrices ${{c}_{i}}:=\sigma _{i}^{+}$ and ${{a}_{i}}:=\sigma _{i}^{-}$ ($i\in \mathbb{Z}$) [29]. The transformed operators fulfill anti-commutation relations not only locally but also non-locally, and thus represent the stochastic process in terms of a 'fermionic' (field) theory. While the Jordan–Wigner transformation provides an exact reformulation of the stochastic process, its applicability is largely limited to systems with one spatial dimension. Moreover, the transformation requires that the stochastic process is first rewritten in terms of a spin-1/2 chain (typically possible only for processes with a single species). Based on the (coherent) eigenstates of the resulting fermionic creation and annihilation operators, Brunel et al then derived a path integral representation of averaged observables. The paths of these integrals proceed along the values of Grassmann variables. Further information on the Jordan–Wigner transformation, Grassmann path integrals, and alternative approaches can be found in [26–28, 30, 31, 351–353].

In the following, we outline various approaches that have recently been proposed for the analysis of the generating function's flow equation.

2.3.1. A spectral method for the computation of stationary distributions.

In their study of a linear transcriptional regulatory cascade of genes and proteins, Walczak, Mugler, and Wiggins developed a spectral method for the computation of stationary probability distributions [34]. They described the regulatory cascade on a coarse-grained level in terms of the copy numbers of certain chemical 'species' at its individual steps. By imposing a Markov approximation, the dynamics of the cascade was reduced to a succession of two-species master equations. The solution of each master equation served as input for the next equation downstream. Every of the reduced master equations allowed for the following processes: First, each of the master equation's two species $i\in \left\{1,2\right\}$ is produced in a one-step process whose rate ${{\gamma}_{i}}\left({{n}_{1}}\right)$ depends only on the copy number n₁ of the species coming earlier along the cascade. Second, each of the two species degrades at a constant per-capita rate ${{\mu}_{i}}$. With $\boldsymbol{n}:={{\left({{n}_{1}},{{n}_{2}}\right)}^{\intercal}}\in \mathbb{N}_{0}^{2}$, the corresponding master equation reads

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{\tau}}p\left(\tau,\boldsymbol{n}|\centerdot \right) & ={{\gamma}_{1}}\left({{n}_{1}}-1\right)p\left(\tau,\boldsymbol{n}-{{\widehat{\boldsymbol{e}}}_{1}}|\centerdot \right)-{{\gamma}_{1}}\left({{n}_{1}}\right)p\left(\tau,\boldsymbol{n}|\centerdot \right) \\ {} & \,\,\,\,\,\,+\,\,{{\gamma}_{2}}\left({{n}_{1}}\right)\left[p\left(\tau,\boldsymbol{n}-{{\widehat{\boldsymbol{e}}}_{2}}|\centerdot \right)-p\left(\tau,\boldsymbol{n}|\centerdot \right)\right] \\ {} & \,\,\,\,\,\,+\,\,{{\mu}_{1}}\left[\left({{n}_{1}}+1\right)p\left(\tau,\boldsymbol{n}+{{\widehat{\boldsymbol{e}}}_{1}}|\centerdot \right)-{{n}_{1}}p\left(\tau,\boldsymbol{n}|\centerdot \right)\right] \\ {} & \,\,\,\,\,\,+\,\,{{\mu}_{2}}\left[\left({{n}_{2}}+1\right)p\left(\tau,\boldsymbol{n}+{{\widehat{\boldsymbol{e}}}_{2}}|\centerdot \right)-{{n}_{2}}p\left(\tau,\boldsymbol{n}|\centerdot \right)\right]. \end{array}\end{eqnarray} \tag{ 93 }$$

Here, the unit vector ${{\widehat{\boldsymbol{e}}}_{i}}$ points in the direction of the ith species. The master equation can be cast into a flow equation for the generating function $|g(\tau |\centerdot )\rangle ={\sum}_{\boldsymbol{n}}p\left(\tau,\boldsymbol{n}|\centerdot \right)|\boldsymbol{n}\rangle$ by following the steps in section 2.1. For this purpose, we choose the basis function as a multivariate extension of the basis function from the introduction to section 2, i.e. as $|\boldsymbol{n}{{\rangle}_{\boldsymbol{q}}}:=\,|{{n}_{1}}{{\rangle}_{{{q}_{1}}}}|{{n}_{2}}{{\rangle}_{{{q}_{2}}}}$ with

$$\begin{eqnarray}&&|{{n}_{i}}{{\rangle}_{{{q}_{i}}}}:={{\left({{q}_{i}}+1\right)}^{{{n}_{i}}}}{{\text{e}}^{-\frac{{{{\bar{\gamma}}}_{i}}}{{{\mu}_{i}}}\left({{q}_{i}}+1\right)}}.\end{eqnarray} \tag{ 94 }$$

The values of the auxiliary parameters ${{\bar{\gamma}}_{1}}$ and ${{\bar{\gamma}}_{2}}$ are only specified in a numerical implementation of the method and affect its stability. Information on how their values are chosen is provided in [34]. For our present purposes, the values of the parameters remain unspecified. By differentiating the generating function with respect to time, one finds that the generating function obeys the flow equation ${{\partial}_{\tau}}|g\rangle =\left({{\mathcal{Q}}_{0}}+{{\mathcal{Q}}_{1}}\right)|g\rangle$ with the transition operators

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{Q}}_{0}}:=-\underset{i\in \left\{1,2\right\}}{\sum}\,{{\mu}_{i}}{{q}_{i}}{{\partial}_{{{q}_{i}}}}~\text{and} \end{array}\end{eqnarray} \tag{ 95 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{Q}}_{1}}:=-\underset{i\in \left\{1,2\right\}}{\sum}\,{{q}_{i}}\left({{{\bar{\gamma}}}_{i}}-{{{\hat{\gamma}}}_{i}}\right). \end{array}\end{eqnarray} \tag{ 96 }$$

Here, the two new operators ${{\hat{\gamma}}_{1}}$ and ${{\hat{\gamma}}_{2}}$ are defined in terms of their actions ${{\hat{\gamma}}_{i}}|{{n}_{1}}\rangle ={{\gamma}_{i}}\left({{n}_{1}}\right)|{{n}_{1}}\rangle$ on the basis functions. Surprisingly, the explicit form of these operators is not needed. Thus, the spectral method even allows for non-polynomial growth rates ${{\gamma}_{i}}\left({{n}_{1}}\right)$.

The operator ${{\mathcal{Q}}_{0}}$ has the same form as the transition operator of the bi-directional reaction $\varnothing \rightleftharpoons A$ from the introductory example. Without the perturbation ${{\mathcal{Q}}_{1}}$, the above flow equation could thus be solved by extending the previous ansatz (38) to two species. To accommodate ${{\mathcal{Q}}_{1}}$ as well, it proves useful to generalize that ansatz to

$$\begin{eqnarray}&&|g(\tau |\centerdot )\rangle =\underset{\boldsymbol{k}\in \mathbb{N}_{0}^{2}}{\sum}\,|\boldsymbol{k}\rangle \rangle {{G}_{\boldsymbol{k}}}(\tau |\centerdot ),\end{eqnarray} \tag{ 97 }$$

with yet to be determined expansion coefficients ${{G}_{\boldsymbol{k}}}(\tau |\centerdot )$. The auxiliary ket is defined as $|\boldsymbol{k}\rangle {{\rangle}_{\boldsymbol{q}}}:=\,|{{k}_{1}}\rangle {{\rangle}_{{{q}_{1}}}}|{{k}_{2}}\rangle {{\rangle}_{{{q}_{2}}}}$ with $|{{k}_{i}}\rangle {{\rangle}_{{{q}_{i}}}}=q_{i}^{{{k}_{i}}}$ and is orthogonal to the bra $\langle \langle \boldsymbol{k}|:=\langle \langle {{k}_{1}}|\langle \langle {{k}_{2}}|$ with $\langle \langle {{k}_{i}}|\,f:=\frac{1}{{{k}_{i}}!}\partial _{{{q}_{i}}}^{{{k}_{i}}}\,f\left({{q}_{i}}\right){{|}_{{{q}_{i}}=0}}$. These bras can be used to extract the expansion coefficient via ${{G}_{\boldsymbol{k}}}(\tau |\centerdot )=\langle \langle \boldsymbol{k}|g(\tau |\centerdot )\rangle$. Differentiation of this coefficient with respect to time and imposing stationarity eventually results in a recurrence relation for ${{G}_{\boldsymbol{k}}}$ that can be solved iteratively. The steady-state probability distribution of the stochastic process is then recovered from the generating function via the inverse transformation (48).

In [34], Walczak et al computed the steady-state distribution of the transcriptional regulatory cascade by solving the recurrence relation for $G_{\boldsymbol{k}}^{s}$ numerically. The resulting distribution was compared to distributions acquired via an iterative method and via the stochastic simulation algorithm (SSA) of Gillespie. The spectral method was found to be about 10⁸ times faster than the SSA in achieving the same accuracy. But as we already mentioned in the introductory section 1.4, the SSA generally performs poorly in the estimation of full distributions, especially in the estimation of their tails. As the spectral method has only been applied to cascades with steady-state copy numbers below $n\approx 30$ thus far, it may be challenged by a direct integration of the two-species master equation (after introducing a reasonable cut-off in the copy numbers). The integration of  ∼1000 coupled ODEs does not pose a problem for modern integrators and the integration provides the full temporal dynamics of the process (see [271] for efficient algorithms). In its current formulation, the spectral method is limited to the evaluation of steady-state distributions and to simple one-step birth-death dynamics. It would be interesting to advance the method for the application to more complex processes (possibly with spatial degrees of freedom), and to extend its scope to the temporal evolution of distributions.

2.3.2. WKB approximations and related approaches.

The probability distribution describing the transcriptional regulatory cascade from the previous section approaches a non-trivial stationary shape in the asymptotic time limit $\tau \to \infty$ (see figure 2 in [34]). Often, however, a non-trivial shape of the probability distribution persists only transiently and is said to be quasi-stationary or metastable. The lifetime and shape of such a distribution can often be approximated in terms of a WKB approximation [84]. A WKB approximation of a jump process starts out from an exponential (eikonal) ansatz for the shape of the metastable probability distribution ('real-space' approach) or for the generating function discussed in the previous sections ('momentum-space' approach). Information on the real-space approach can be found in [63, 64, 146, 171, 172, 177, 354–363]. The recent review of Assaf and Meerson provides an in-depth discussion of the applicability of the real- and momentum-space approaches [364].

In the following, we outline the momentum-space WKB approximation for a system in which particles of type A annihilate in the binary reaction $2\,A\to \varnothing$ with rate coefficient μ, and are replenished in the linear reaction $A\to 2\,A$ with rate coefficient $\gamma \gg \mu$. According to a deterministic model of the combined processes with rate equation ${{\partial}_{\tau}}\bar{n}=\gamma \bar{n}-2\mu {{\bar{n}}^{2}}$, the particle number $\bar{n}$ converges to an asymptotic value ${{\bar{n}}_{\infty}}=\gamma /(2\mu )\gg 1$. However, a numerical integration of the (truncated) master equation of the stochastic process shows that the mean particle number stays close to ${{\bar{n}}_{\infty}}$ only for a long but finite time (see figure 3). Asymptotically, all particles become trapped in the 'absorbing' state n  =  0. Consequently, the conditional probability distribution converges to $p\left(\infty,n|{{t}_{0}},{{n}_{0}}\right)={{\delta}_{n,0}}$. Up to a pre-exponential factor, the time after which the absorbing state is reached can be readily estimated using a momentum-space WKB approximation as shown below [35]. The value of the pre-exponential factor was determined by Turner and Malek-Mansour using a recurrence relation [365], by Kessler and Shnerb using a real-space WKB approximation [357], and by Assaf and Meerson upon combining the generating function technique with Sturm-Liouville theory [366] (the latter method was developed in [36, 37]; see also [367]). Using this method, Assaf and Meerson also succeeded in computing the shape of the metastable distribution. Moreover, Assaf et al showed how a momentum-space WKB approximation can be used to determine mean extinction times for processes with time-modulated rate coefficients [38].

Zoom In Zoom Out Reset image size

Upon rescaling time as $\gamma \,\tau \to \tau$, the chemical master equation (24) of the combined processes $A\to 2\,A$ and $2\,A\to \varnothing$ translates into the flow equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}|g\rangle =\left[\left({{q}^{2}}-q\right){{\partial}_{q}}+\frac{1}{2{{{\bar{n}}}_{\infty}}}\left(1-{{q}^{2}}\right)\partial _{q}^{2}\right]|g\rangle \end{eqnarray} \tag{ 98 }$$

of the ordinary generating function $g(\tau ;q|\centerdot )={\sum}_{n}{{q}^{n}}p(\tau,n|\centerdot )$ (see sections 2.1 and 2.2.2). A WKB approximation of the flow equation can be performed by inserting the ansatz

$$\begin{eqnarray}&&g(\tau ;q|\centerdot )={{\text{e}}^{S(\tau,q)}}\end{eqnarray} \tag{ 99 }$$

with the 'action'

$$\begin{eqnarray}&&S(\tau,q)=2\,{{\bar{n}}_{\infty}}\underset{k=0}{\overset{\infty}{\sum}}\,\frac{{{S}_{k}}(\tau,q)}{\bar{n}_{\infty}^{k}}\end{eqnarray} \tag{ 100 }$$

into the equation, followed by a successive analysis of terms that are of the same order with respect to the power of the small parameter $1/{{\bar{n}}_{\infty}}$. Note that the exponential ansatz (99) often includes a minus sign in front of the action, which we neglect for convenience. The pre-factor '2' of the action (100) is also included just for convenience. Upon inserting the ansatz (99) into the flow equation (98), one obtains the equality

$$\begin{eqnarray}&&{{\partial}_{\tau}}{{S}_{0}}+\text{O}\left(1/{{\bar{n}}_{\infty}}\right)=H\left(q,{{\partial}_{q}}{{S}_{0}}\right)+\text{O}\left(1/{{\bar{n}}_{\infty}}\right)\end{eqnarray} \tag{ 101 }$$

with the 'Hamiltonian'

$$\begin{eqnarray}&&H(q,x):=\left({{q}^{2}}-q\right)x+\left(1-{{q}^{2}}\right){{x}^{2}}.\end{eqnarray} \tag{ 102 }$$

Thus, at leading order of $1/{{\bar{n}}_{\infty}}$, we have obtained a closed equation for S₀, which has the form of a Hamilton–Jacobi equation [368]. A Hamilton–Jacobi equation can be solved by the method of characteristics [368], with the characteristic curves q(s) and x(s) obeying Hamilton's equations

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\partial}_{s}}x=\frac{\partial H(q,x)}{\partial q}=(2q-1)x-2q{{x}^{2}}\quad \text{and} \end{array}\end{eqnarray} \tag{ 103 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\partial}_{-s}}q=\frac{\partial H(q,x)}{\partial x}=\left({{q}^{2}}-q\right)+2\left(1-{{q}^{2}}\right)x. \end{array}\end{eqnarray} \tag{ 104 }$$

These equations are, for example, solved by q(s)  =  1 and x(s) being a solution of ${{\partial}_{s}}x=x-2{{x}^{2}}$. Note that this equation corresponds to a rescaled rate equation. The Hamiltonian vanishes along the characteristic curve because H(1, x)  =  0. Further 'zero-energy' lines of the Hamiltonian are given by x  =  0 and by $x(q)=\frac{q}{1+q}$. These lines partition the phase portrait of Hamilton's equations into separate regions as shown in figure 4. The path $\left(q,\frac{q}{1+q}\right)$ from the 'active state' $(q,x)=\left(1,\frac{1}{2}\right)$ to the 'passive state' (0,0) constitutes the 'optimal path to extinction' (see below) [35, 38, 369].

Zoom In Zoom Out Reset image size

In addition to Hamilton's equations, the method of characteristics implies that the action S₀ obeys

$$\begin{eqnarray}&&\frac{\text{d}}{\text{d}s}{{S}_{0}}=-\left(x\,{{\partial}_{-s}}q-H(q,x)\right)\end{eqnarray} \tag{ 105 }$$

along characteristic curves (the minus signs are only included to emphasize a similarity with an action encountered later in section 7.1). According to Elgart and Kamenev, the (negative) value of the action (100) along the optimal path to extinction determines the logarithm of the mean extinction time $\bar{\tau}$ from the metastable state in leading order of ${{\bar{n}}_{\infty}}$ [35]. As this path proceeds along a zero-energy line, the action (100) evaluates to

$$\begin{eqnarray}&&S\approx 2{{\bar{n}}_{\infty}}{{S}_{0}}=2{{\bar{n}}_{\infty}}{\int}_{1}^{0}\frac{q\text{d}q}{1+q}=-2{{\bar{n}}_{\infty}}(1-\ln 2).\end{eqnarray} \tag{ 106 }$$

Upon returning to the original time scale via $\tau \to \gamma \,\tau$, the mean extinction time follows as

$$\begin{eqnarray}&&\bar{\tau}=A\,{{\gamma}^{-1}}{{\text{e}}^{2{{{\bar{n}}}_{\infty}}(1-\ln 2)}}\end{eqnarray} \tag{ 107 }$$

with pre-exponential factor A. The value of this factor has been determined in [357, 365, 366] and reads, in terms of our rate coefficients, $A=2\sqrt{\pi /{{{\bar{n}}}_{\infty}}}$. More information on the above procedure can be found in [35, 38]. The method has also been applied to classic epidemiological models [369] and to a variant of the Verhulst logistic model [370] (exact results on this model with added immigration have been obtained in [371] using the generating function technique).

2.3.3. A variational method.

In the present section, we formulated general conditions under which the forward master equation (10) can be transformed into a linear partial differential equation, a 'flow equation', obeyed by the probability generating function

$$\begin{eqnarray}&&|g\left(\tau |{{t}_{0}},{{n}_{0}}\right)\rangle =\underset{n}{\sum}\,|n\rangle \,p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right).\end{eqnarray} \tag{ 108 }$$

First, the conditions (C) and (O) require that there exists a complete and orthogonal basis comprising a set of basis functions $\left\{|n\rangle \right\}$ and a set of basis functionals $\left\{\langle n|\right\}$. The basis functionals recover the conditional probability distribution via $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)=\langle n|g\left(\tau |{{t}_{0}},{{n}_{0}}\right)\rangle$ and the right choice of the basis functions may help to obtain a simplified flow equation. We introduced different bases for the study of random walks, of chemical reactions, and of processes with locally excluding particles in section 2.2. Moreover, the conditions (E) and (Q) require that there exist two differential operators: a basis evolution operator ${{\mathcal{E}}_{\tau}}$ encoding the possible time-dependence of the basis, and a transition operator ${{\mathcal{Q}}_{\tau}}$ encoding the actual dynamics of the process. The generating function (108) then obeys the flow equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}|g\rangle =\left({{\mathcal{E}}_{\tau}}+{{\mathcal{Q}}_{\tau}}\right)|g\rangle ={{\tilde{\mathcal{Q}}}_{\tau}}|g\rangle.\end{eqnarray} \tag{ 109 }$$

Various methods have recently been proposed for the study of such a flow equation and were outlined in section 2.3. These methods include the variational approach of Eyink [83] and of Sasai and Wolynes [39], WKB approximations and spectral formulations of Elgart and Kamenev and of Assaf and Meerson [35–37], and the spectral method of Walczak, Mugler, and Wiggins [34]. Thus far, most of these methods have only been applied to systems without spatial degrees of freedom and with only one or a few types of particles. Future research is needed to overcome these limitations. Later, in section 5, we show how the solution of the flow equation (109) can be represented by a path integral. The evaluation of the path integral is demonstrated in computing the generating function of general linear processes. Moreover, we explain in section 7 how this path integral connects to a recent method of Elgart and Kamenev for the computation of rare event probabilities [35]. One can also use the path integral to derive a path integral representation of averaged observables. But a simpler route to this representation starts out from a different flow equation, a flow equation obeyed by the 'marginalized distribution'.

The generalized probability generating function (45) was defined by summing the conditional probability distribution $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$ over a set of basis functions, the kets $\left\{|n\rangle \right\}$. In the following, we consider the distribution $p\left(t,n|\tau,{{n}_{0}}\right)$ instead, i.e. we fix the final time t while keeping the initial time τ variable. By summing this distribution over a set of basis functions ${{\{|{{n}_{0}}\rangle}_{\tau}}\}$, one can define the series

$$\begin{eqnarray}&&|\,p(t,n|\tau )\rangle :=\underset{{{n}_{0}}}{\sum}\,p\left(t,n|\tau,{{n}_{0}}\right)|{{n}_{0}}\rangle.\end{eqnarray} \tag{ 110 }$$

As before, the subscript denoting the time-dependence of the basis function is typically dropped. The variables n and n₀ again represent states from some countable state space. Assuming that the basis functions $\left\{|n\rangle \right\}$ and appropriately chosen basis functionals $\left\{\langle n|\right\}$ form a (C)omplete and (O)rthogonal basis, the conditional probability distribution can be recovered from (110) via

$$\begin{eqnarray}&&p\left(t,n|\tau,{{n}_{0}}\right)=\langle {{n}_{0}}|\,p(t,n|\tau )\rangle.\end{eqnarray} \tag{ 111 }$$

We call the function $|\,p(t,n|\tau )\rangle$ a 'marginalized distribution' because it proves most useful when the summation in (110) constitutes a marginalization of the conditional probability distribution $p\left(t,n|\tau,{{n}_{0}}\right)$ over a probability distribution $|{{n}_{0}}\rangle$. The marginalized distribution then represents a 'single-time distribution' with respect to the random variable n in the sense of section 1.3. A basis function to which these considerations apply is the 'Poisson basis function' $|{{n}_{0}}{{\rangle}_{x}}:=\frac{{{x}^{{{n}_{0}}}}{{\text{e}}^{-x}}}{{{n}_{0}}!}$. We make heavy use of this basis function in the study of chemical reactions in section 3.2.2. Since the definition of the marginalized distribution in (110) does not affect the random variable n, the marginalized distribution of course solves the forward master equation ${{\partial}_{t}}|\,p(t,n|\tau )\rangle ={\sum}_{m}{{Q}_{t}}(n,m)|\,p(t,m|\tau )\rangle$ with the initial condition $|\,p(\tau,n|\tau )\rangle =\,|{{n}_{0}}{{\rangle}_{\tau}}$. In the following, we formulate conditions under which $|\,p\rangle$ also obeys a linear PDE evolving backward in time.

Before proceeding, let us briefly note that if the basis function of the marginalized distribution is not chosen as a probability distribution, the name marginalized 'distribution' is somewhat of a misnomer; that is, for example, the case for the Fourier basis $|{{n}_{0}}{{\rangle}_{x}}:={{\text{e}}^{\text{i}{{n}_{0}}x}}$, which we consider in section 3.2.1.

The derivation of the linear PDE obeyed by the marginalized distribution proceeds analogously to the derivation in section 2.1. But instead of employing the forward master equation, we now employ the backward master equation ${{\partial}_{-\tau}}p(t|\tau )=p(t|\tau ){{Q}_{\tau}}$ for this purpose (recall that $p(t|\tau )$ is the matrix with elements $p\left(t,n|\tau,{{n}_{0}}\right)$). Upon differentiating the marginalized distribution (110) with respect to the time parameter τ, one obtains the equation

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\partial}_{-\tau}}|\,p\rangle =\underset{{{n}_{0}}}{\sum}\,p\left(\centerdot |\tau,{{n}_{0}}\right)\left(-{{\partial}_{\tau}}|{{n}_{0}}\rangle +\underset{m}{\sum}\,|m\rangle Q_{\tau}^{\intercal}\left(m,{{n}_{0}}\right)\right). \end{array}\end{eqnarray} \tag{ 112 }$$

The rate $Q_{\tau}^{\intercal}\left(m,{{n}_{0}}\right)={{Q}_{\tau}}\left({{n}_{0}},m\right)$ represents an element of the transposed transition matrix $Q_{\tau}^{\intercal}$.

As in section 2.1, two differential operators are required to turn the above expression into a linear PDE. First, we require a basis evolution operator ${{\mathcal{E}}_{\tau}}\left(x,{{\partial}_{x}}\right)$ fulfilling ${{\mathcal{E}}_{\tau}}|{{n}_{0}}\rangle ={{\partial}_{\tau}}|{{n}_{0}}\rangle$ for all values of n₀. The previous evolution operator (E) serves this purpose. A second differential operator $\mathcal{Q}_{\tau}^{\dagger}\left(x,{{\partial}_{x}}\right)$ is required to encode the information stored in the transition matrix. This operator should be a power series in ${{\partial}_{x}}$ and fulfil, for all n₀,

$$\begin{eqnarray}&&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\mathcal{Q}_{\tau}^{\dagger}|{{n}_{0}}\rangle =\underset{m}{\sum}\,|m\rangle Q_{\tau}^{\intercal}\left(m,{{n}_{0}}\right)\quad \quad \quad \,\,\,\,\,\quad \quad \quad \left({{\text{Q}}^{\intercal}}\right).\end{eqnarray}$$

By the (C)ompleteness of the basis, one could also define this operator constructively as

$$\begin{eqnarray}&&\mathcal{Q}_{\tau}^{\dagger}:=\underset{m,{{n}_{0}}}{\sum}\,|m\rangle Q_{\tau}^{\intercal}\left(m,{{n}_{0}}\right)\langle {{n}_{0}}|.\end{eqnarray} \tag{ 114 }$$

We wrote these expressions in terms of the transposed transition matrix because for the unit column vectors $|m\rangle ={{\widehat{\boldsymbol{e}}}_{m}}$ and the unit row vectors $\langle n|=\widehat{\boldsymbol{e}}_{n}^{\intercal}$, $\mathcal{Q}_{\tau}^{\dagger}$ and $Q_{\tau}^{\intercal}$ coincide (with $m,n\in {{\mathbb{N}}_{0}}$). As we mostly consider bases whose kets and bras represent functions and functionals, we call $\mathcal{Q}_{\tau}^{\dagger}$ the 'adjoint' transition operator in the following. Often, an operator ${{O}^{\dagger}}\left(x,{{\partial}_{x}}\right)$ is said to be the adjoint of an operator $O\left(x,{{\partial}_{x}}\right)$ if the following relation holds with respect to two test functions f and g:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {\int}^{}\text{d}x\,\left[{{O}^{\dagger}}\left(x,{{\partial}_{x}}\right)f(x)\right]g(x) \\ ={\int}^{}\text{d}x\,f(x)\left[O\left(x,{{\partial}_{x}}\right)g(x)\right]. \end{array}\end{eqnarray} \tag{ 115 }$$

However, whether an operator ${{\mathcal{Q}}_{\tau}}$ complementing the above $\mathcal{Q}_{\tau}^{\dagger}$ actually exists will not be important in the following (except in our discussion of the Poisson representation in section 3.5).

Provided that both a basis evolution operator ${{\mathcal{E}}_{\tau}}$ and an adjoint transition operator $\mathcal{Q}_{\tau}^{\dagger}$ are found for a particular process, it follows from the backward-time equation (112) that the marginalized distribution $|\,p(t,n|\tau )\rangle$ obeys the flow equation ⁸

$$\begin{eqnarray}&&{{\partial}_{-\tau}}|\,p\rangle =\left(-{{\mathcal{E}}_{\tau}}+\mathcal{Q}_{\tau}^{\dagger}\right)|\,p\rangle =:\tilde{\mathcal{Q}}_{\tau}^{\dagger}|\,p\rangle.\end{eqnarray} \tag{ 116 }$$

The evolution of this equation proceeds backward in time, starting out from the final condition $|\,p(t,n|t)\rangle =\,|n\rangle$.

In section 4, we show how the flow equation (116) can be solved in terms of a 'backward' path integral. Provided that the basis function $|n\rangle$ is chosen as a probability distribution, this path integral represents a true probability distribution: the marginalized distribution. The fact that the backward path integral represents a probability distribution distinguishes it from the 'forward' path integral in section 5. The forward path integral represents the probability generating function (45). Both the forward path integral and the backward path integral can be used to derive a path integral representation of averaged observables, as we show in section 6. The derivation of this representation from the backward path integral, however, is significantly easier. In fact, the path integral representation of the average of an observable A will follow directly by summing the backward path integral representation of the marginalized distribution over A(n), i.e. via $\langle A\rangle ={\sum}_{n}A(n)|\,p(t,n|\tau )\rangle$. Note that this average also obeys the flow equation (116). In section 3.3, we demonstrate how the flow equation can be used to compute mean extinction times.

Before introducing bases for the analysis of different stochastic processes, let us briefly note that if a process is homogeneous in time, its marginalized distribution depends only on the difference $t-\tau$. The above flow equation can then be rewritten so that it evolves $|\,p(\tau,n|0)\rangle$ forward in time τ, starting out from the initial condition $|\,p(0,n|0)\rangle =\,|n\rangle$. In section 3.3, we make use of this property to compute mean extinction times.

To demonstrate the application of the marginalized distribution (110), let us reconsider the random walk from section 2.2.1 and the chemical reaction from section 2.2.2. The 'Poisson basis function' introduced in the latter section will be employed in the computation of mean extinction times in section 3.3. Moreover, it will allow us to recover the Poisson representation of Gardiner and Chaturvedi [32, 33] in section 3.5.

3.2.1. Random walks.

The following solution of the random walk largely parallels the previous derivation in section 2.2.1. In particular, we again use the orthogonal and complete Fourier basis $|n{{\rangle}_{x}}={{\text{e}}^{\text{i}nx}}$ and $\langle n|\,f={\int}_{-\pi}^{\pi}\frac{\text{d}x}{2\pi}{{\text{e}}^{-\text{i}nx}}f(x)$ with $n\in \mathbb{Z}$. Due to the time-independence of the basis, the (E)volution operator ${{\mathcal{E}}_{\tau}}$ is zero. The condition (${{\text{Q}}^{\intercal}}$) on the adjoint transition operator $\mathcal{Q}_{\tau}^{\dagger}$ is specified by the transition matrix (59) of the process and reads

$$\begin{eqnarray}&&\mathcal{Q}_{\tau}^{\dagger}|{{n}_{0}}\rangle ={{r}_{\tau}}\left(|{{n}_{0}}-1\rangle \,-\,|{{n}_{0}}\rangle \right)\,+\,{{l}_{\tau}}\left(|{{n}_{0}}+1\rangle \,-|{{n}_{0}}\rangle \right).\end{eqnarray} \tag{ 117 }$$

As before, we employ the operators $c(x)={{\text{e}}^{\text{i}x}}$ and $a(x)={{\text{e}}^{-\text{i}x}}$, which act on the basis functions as $c|n\rangle =\,|n+1\rangle$ and $a|n\rangle =\,|n-1\rangle$, respectively. Therefore, the adjoint transition operator with the above property can be defined as

$$\begin{eqnarray}&&\mathcal{Q}_{\tau}^{\dagger}\left(x,{{\partial}_{x}}\right):={{r}_{\tau}}\left(a(x)-1\right)+{{l}_{\tau}}\left(c(x)-1\right).\end{eqnarray} \tag{ 118 }$$

As this operator does not contain any derivatives, it is self-adjoint in the sense of (115). Due to a mismatch in signs, it is, however, not the adjoint of the previous operator ${{\mathcal{Q}}_{\tau}}$ in (61). This mismatch could be corrected by redefining the above basis function as $|n{{\rangle}_{x}}:={{\text{e}}^{-\text{i}nx}}$. Ignoring this circumstance, the flow equation of the marginalized distribution follows as

$$\begin{eqnarray}&&{{\partial}_{-\tau}}|\,p\rangle =\left[{{r}_{\tau}}\left({{\text{e}}^{-\text{i}x}}-1\right)+{{l}_{\tau}}\left({{\text{e}}^{\text{i}x}}-1\right)\right]|\,p\rangle,\end{eqnarray} \tag{ 119 }$$

and is solved by

$$\begin{eqnarray}&&|\,p\rangle =\exp \left(\left({{\text{e}}^{-\text{i}x}}-1\right){\int}_{\tau}^{t}\text{d}s\,{{r}_{s}}+\left({{\text{e}}^{\text{i}x}}-1\right){\int}_{\tau}^{t}\text{d}s\,{{l}_{s}}\right)|n\rangle.\end{eqnarray} \tag{ 120 }$$

The conditional probability distribution is recovered via the inverse Fourier transformation $p\left(t,n|\tau,{{n}_{0}}\right)=\langle {{n}_{0}}|\,p\rangle$. Upon inserting the explicit representation of the basis, the derivation proceeds as in appendix C (with the substitutions $x\to -q$, $\tau \to {{t}_{0}}$ and $t\to \tau$). Eventually, one recovers a Skellam distribution as the solution of the process.

3.2.2. Chemical reactions.

To prepare the computation of mean extinction times in the next section as well as the derivation of the Poisson representation in section 3.5, we now reconsider processes that can be decomposed additively into chemical reactions of the form $k\,A\to l\,A$. Our later derivation of a path integral representation of averaged observables is also restricted to such processes. The state space of the number of molecules is again ${{\mathbb{N}}_{0}}$. Analogous to section 2.2.2, we require an (O)rthogonal and (C)omplete basis, as well as an (E)volution operator ${{\mathcal{E}}_{\tau}}$ and an adjoint transition operator $\mathcal{Q}_{\tau}^{\dagger}$ (condition (${{\text{Q}}^{\intercal}}$)) to specify the flow equation of the marginalized distribution.

As discussed in section 2.2.2, the elements of the transition rate matrix of the reaction $k\,A\to l\,A$ with rate coefficient ${{\gamma}_{\tau}}$ are given by ${{Q}_{\tau}}(m,n)={{\gamma}_{\tau}}{{(n)}_{k}}\left({{\delta}_{m,n-k+l}}-{{\delta}_{m,n}}\right)$. The condition (${{\text{Q}}^{\intercal}}$) on the adjoint transition operator therefore reads

$$\begin{eqnarray}&&\mathcal{Q}_{\tau}^{\dagger}|{{n}_{0}}\rangle ={{\gamma}_{\tau}}\left({{\left({{n}_{0}}-l+k\right)}_{k}}|{{n}_{0}}-l+k\rangle -{{\left({{n}_{0}}\right)}_{k}}|{{n}_{0}}\rangle \right).\end{eqnarray} \tag{ 121 }$$

This condition is met by

$$\begin{eqnarray}&&\mathcal{Q}_{\tau}^{\dagger}(c,a):={{\gamma}_{\tau}}{{c}^{k}}\left({{a}^{l}}-{{a}^{k}}\right),\end{eqnarray} \tag{ 122 }$$

provided that there exist operators c and a fulfilling

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} c|n\rangle =(n+1)|n+1\rangle ~\text{and} \end{array}\end{eqnarray} \tag{ 123 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} a|n\rangle =\,|n-1\rangle, \end{array}\end{eqnarray} \tag{ 124 }$$

respectively. We again call c the creation and a the annihilation operator, even though the pre-factors in the above relations differ from the ones in (65) and (66). Similar relations also hold with respect to the basis functionals, namely $\langle n|c=n\langle n-1|$ and $\langle n|a=\langle n+1|$ (assuming the orthogonality of the basis). The operators also fulfil the commutation relation [a, c]  =  1.

The actions of the creation and annihilation operators on the basis functions hint at how the basis functions and functionals from section 2.2.2 can be adapted to meet the present requirements. In particular, the relations (123) and (124) can be fulfilled by moving the factorial from the basis functional (76) to the basis function (70) so that

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} |n{{\rangle}_{x}}:=\frac{{{\left(\zeta x+\tilde{{x}}\,\right)}^{n}}{{\text{e}}^{-\tilde{{q}}\,\left(\zeta x+\tilde{{x}}\,\right)}}}{n!}\,\text{and} \end{array}\end{eqnarray} \tag{ 125 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle n|\,f:={{a}^{n}}f(x){{|}_{x=-\tilde{x}/\zeta}}. \end{array}\end{eqnarray} \tag{ 126 }$$

Let us briefly note that we changed the argument of the basis function from q to x as compared to section 2.2.2 because both the generating function approach and the marginalized distribution approach thereby result in the same path integral representation of averaged observables (see section 6). Apart from this notational change, the basis evolution operator

$$\begin{eqnarray}&&{{\mathcal{E}}_{\tau}}(c,a):=\left({{\partial}_{\tau}}\tilde{x}\right)\left(a-\tilde{q}\right)-\left({{\partial}_{\tau}}\tilde{q}\right)c\end{eqnarray} \tag{ 127 }$$

keeps the form it had in (77). The creation and annihilation operators also keep their previous forms in (73) and (74), i.e.

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} c\left(x,{{\partial}_{x}}\right):=\zeta x+\tilde{x}\,\text{and} \end{array}\end{eqnarray} \tag{ 128 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} a\left(x,{{\partial}_{x}}\right):={{\partial}_{\zeta x}}\,+\,\tilde{{q}}\,. \end{array}\end{eqnarray} \tag{ 129 }$$

Despite their similar appearance, the operator $\mathcal{Q}_{\tau}^{\dagger}$ in (122) is not the adjoint of the operator ${{\mathcal{Q}}_{\tau}}$ in (64). Nevertheless, the two operators fulfil ${{\mathcal{Q}}_{\tau}}(q,x)=\mathcal{Q}_{\tau}^{\dagger}(x,q)$ for scalar arguments. This relation is essentially the reason why we interchanged the letters x and q as compared to section 2.2.2. Both the generating function approach and the marginalized distribution approach thereby lead to the same path integral representation of averaged observables, as will be shown in section 6.1.

The choice of the parameters $\zeta \ne 0$, $\tilde{x}(\tau )$, and $\tilde{q}(\tau )$ in the basis function (125) depends on the problem at hand. In section 7, $\tilde{x}$ and $\tilde{q}$ will act as 'stationary' or 'extremal' paths, with x and an auxiliary variable q being deviations from them. For $\zeta :=\tilde{q}:=1$ and $\tilde{x}:=0$, the basis function instead simplifies to the Poisson distribution

$$\begin{eqnarray}&&|n{{\rangle}_{x}}=\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}.\end{eqnarray} \tag{ 130 }$$

In the following, we make heavy use of this 'Poisson basis function'. It will play a crucial role in the formulation of a path integral representation of averaged observables in section 6 and in recovering the Poisson representation in section 3.5.

Let us demonstrate the use of the Poisson basis function for the linear decay process $A\to \varnothing$ with rate coefficient ${{\mu}_{\tau}}$. For the above choice of ζ, $\tilde{x}$, and $\tilde{q}$, the creation operator (128) reads c  =  x and the annihilation operator (129) reads $a={{\partial}_{x}}\,+\,1$. With the transition operator (122), the flow equation (116) obeyed by the marginalized distribution $|\,p(t,n|\tau ){{\rangle}_{x}}$ follows as

$$\begin{eqnarray}&&{{\partial}_{-\tau}}|\,p\rangle =-{{\mu}_{\tau}}x{{\partial}_{x}}|\,p\rangle.\end{eqnarray} \tag{ 131 }$$

This equation is solved by the Poisson distribution

$$\begin{eqnarray}&&|\,p(t,n|\tau ){{\rangle}_{x}}=\frac{{{\left({{\alpha}_{t,\tau}}x\right)}^{n}}{{\text{e}}^{-{{\alpha}_{t,\tau}}x}}}{n!}\end{eqnarray} \tag{ 132 }$$

whose mean ${{\alpha}_{t,\tau}}\,x$ decays proportionally to ${{\alpha}_{t,\tau}}:={{\text{e}}^{-{\int}_{\tau}^{t}\text{d}s\,{{\mu}_{s}}}}$. To interpret this solution, let us recall that the sum in the definition of the marginalized distribution (110) does not affect the particle number n. Therefore, upon relabelling the time parameters, the above solution (132) also solves the forward master equation of the process, namely

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\partial}_{\tau}}|\,p\left(\tau,n|{{t}_{0}}\right){{\rangle}_{x}} \\ {} & ={{\mu}_{\tau}}{{((n+1)|\,p\left(\tau,n+1|{{t}_{0}}\right)\rangle}_{x}}-n|\,p\left(\tau,n|{{t}_{0}}\right){{\rangle}_{x}}). \end{array}\end{eqnarray} \tag{ 133 }$$

Unlike the conditional distribution $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)$, however, the marginalized distribution $|\,p\left(\tau,n|{{t}_{0}}\right){{\rangle}_{x}}$ describes the dynamics of a population whose particle number at time t₀ is Poisson distributed with mean x. The conditional distribution is recovered from it via the inverse transformation $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)=\langle {{n}_{0}}|\,p\left(\tau,n|{{t}_{0}}\right)\rangle$ with $\langle {{n}_{0}}|\,f={{\left({{\partial}_{x}}+1\right)}^{{{n}_{0}}}}f(x){{|}_{x=0}}$. This transformation results in the Binomial distribution

$$\begin{eqnarray}&&p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)=\left(\begin{array}{c} {{n}_{0}} \\ n \end{array}\right){{\left({{\alpha}_{\tau,{{t}_{0}}}}\right)}^{n}}{{\left(1-{{\alpha}_{\tau,{{t}_{0}}}}\right)}^{{{n}_{0}}-n}}.\end{eqnarray} \tag{ 134 }$$

Both the mean value ${{\text{e}}^{-{\int}_{{{t}_{0}}}^{\tau}\text{d}s\,{{\mu}_{s}}}}{{n}_{0}}$ and the variance $\left(1-{{\text{e}}^{-{\int}_{{{t}_{0}}}^{\tau}\text{d}s\,{{\mu}_{s}}}}\right)\,{{\text{e}}^{-{\int}_{{{t}_{0}}}^{\tau}\text{d}s\,{{\mu}_{s}}}}{{n}_{0}}$ of this distribution decay exponentially for large times, provided that ${{\mu}_{s}}>0$.

If at most two reactants and two products are involved in a reaction, the flow equation (116) of the marginalized distribution has the mathematical form of a backward Fokker–Planck equation. That is, for example, the case for the coagulation reaction $2\,A\to A$. For the Poisson basis function, the marginalized distribution $|\,p(t,n|\tau ){{\rangle}_{x}}$ of this process obeys the flow equation

$$\begin{eqnarray}&&{{\partial}_{-\tau}}|\,p\rangle ={{\alpha}_{\tau}}(x){{\partial}_{x}}|\,p\rangle +\frac{1}{2}{{\beta}_{\tau}}(x)\partial _{x}^{2}|\,p\rangle,\end{eqnarray} \tag{ 135 }$$

with the drift coefficient ${{\alpha}_{\tau}}(x):=-{{\mu}_{\tau}}{{x}^{2}}$ and the diffusion coefficient ${{\beta}_{\tau}}(x):=-2{{\mu}_{\tau}}{{x}^{2}}$. Unlike the diffusion coefficient of the 'true' backward Fokker–Planck equation (2), this diffusion coefficient may be negative (e.g. for $x\in \mathbb{R}\backslash \{0\}$). The final condition $|\,p(t,n|t){{\rangle}_{x}}=\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}$ and the Feynman–Kac formula (A.3) imply that the above flow equation is solved by

$$\begin{eqnarray}&&|\,p(t,n|\tau ){{\rangle}_{x}}={{\langle \langle \frac{x{{(t)}^{n}}{{\text{e}}^{-x(t)}}}{n!}\rangle \rangle}_{W}}.\end{eqnarray} \tag{ 136 }$$

Here, x(s) obeys the Itô SDE

$$\begin{eqnarray}&&\text{d}x(s)={{\alpha}_{s}}\left(x(s)\right)\text{d}s+\sqrt{{{\beta}_{s}}\left(x(s)\right)}\text{d}W(s),\end{eqnarray} \tag{ 137 }$$

which evolves x(s) from $x(\tau )=x$ to x(t). The symbol ${{\langle \langle \centerdot \rangle \rangle}_{W}}$ represents an average over realizations of the Wiener process W. Since the diffusion coefficient ${{\beta}_{\tau}}(x)=-2{{\mu}_{\tau}}{{x}^{2}}$ of the coagulation process can take on negative values, the sample paths of the SDE may acquire imaginary components. This circumstance does not prevent the use of (136) for the calculation of $|\,p(t,n|\tau )\rangle$, although it may complicate numerical evaluations. Recently, Wiese attempted the evaluation of the average in (136) via the generation of sample paths [336]. Over a short time interval $[\tau,t]$, he found a good agreement between the resulting distribution and a distribution effectively acquired via the stochastic simulation algorithm (SSA) in section 1.3. Over larger time intervals, however, the integration of the SDE encountered problems regarding its numerical convergence. Future research is needed to overcome this limitation. Moreover, it remains an open challenge to specify the boundary conditions of the PDE (135) to enable its direct numerical integration (an analogous problem is encountered for the flow equation of the generating function, see [366]).

One often wishes to know the mean time at which a process first hits some target in state space. Such a target could, for example, be a state in which no more particles are left in the system. If the particles only replenish through auto-catalysis, the process will then come to a halt. The mean time after which that happens is called the mean extinction time. For Markov processes with continuous sample paths, mean extinction times and, more generally, first-passage times, are commonly calculated with the help of the backward Fokker–Planck equation (2). For jump processes, one can use the backward master equation for this purpose. The calculation, however, is typically feasible only for one-step processes and involves the solution of a recurrence relation [110, 215]. In [374], Drummond et al recently showed how the calculation can be simplified using the Poisson representation of Gardiner and Chaturvedi [32, 33]. We introduce this representation in section 3.5. In the following, we outline how mean extinction times can instead be inferred in an analogous way using the marginalized distribution from the previous sections. For the Poisson basis function $|{{n}_{0}}{{\rangle}_{x}}=\frac{{{x}^{{{n}_{0}}}}{{\text{e}}^{-x}}}{{{n}_{0}}!}$, the marginalized distribution is a true (single-time) probability distribution and is more easily interpreted than the integral kernel of the Poisson representation.

We again consider a process that can be decomposed additively into chemical reactions of the form $k\,A\to l\,A$. In addition, the transition matrix $\mathcal{Q}_{\tau}^{\dagger}={{\mathcal{Q}}^{\dagger}}$ of the process shall now be time-independent and the particles in the system shall be Poisson distributed with mean x at time $\tau =0$. Consequently, the probability of finding n particles in the system at time $\tau \geqslant 0$ is described by the marginalized distribution $|\,p(\tau,n|0){{\rangle}_{x}}$, provided that its basis function is chosen as $|{{n}_{0}}{{\rangle}_{x}}:=\frac{{{x}^{{{n}_{0}}}}{{\text{e}}^{-x}}}{{{n}_{0}}!}$ (see the definition (110) of the marginalized distribution). Since the Poisson basis function is independent of time and the process under consideration homogeneous in time, the marginalized distribution $|\,p(\tau,n|0){{\rangle}_{x}}$ obeys the forward-time flow equation (see (116))

$$\begin{eqnarray}&&{{\partial}_{\tau}}|\,p\rangle ={{\mathcal{Q}}^{\dagger}}|\,p\rangle.\end{eqnarray} \tag{ 138 }$$

We define the probability of finding the system in an 'active' state with n  >  0 particles at time $\tau \geqslant 0$ as

$$\begin{eqnarray}&&\alpha (\tau ;x):=\underset{n=1}{\overset{\infty}{\sum}}\,|\,p(\tau,n|0){{\rangle}_{x}}=1\,-\,|\,p(\tau,0|0){{\rangle}_{x}}.\end{eqnarray} \tag{ 139 }$$

Over time, this probability flows into the 'absorbing' state n  =  0 at rate $f(\tau,x):=-{{\partial}_{\tau}}\alpha (\tau ;x)$. Since $f(\tau,x)\text{d}\tau$ is the probability of becoming absorbed during the time interval $[\tau,\tau + \Delta t]$, one can define the mean extinction time as ${{\langle \tau \rangle}_{x}}:={\int}_{0}^{\infty}\text{d}\tau \,\tau \,f(\tau,x)$. An integration by parts transforms this average into

$$\begin{eqnarray}&&{{\langle \tau \rangle}_{x}}=-{\int}_{0}^{\infty}\text{d}\tau \,\tau {{\partial}_{\tau}}\alpha (\tau ;x)={\int}_{0}^{\infty}\text{d}\tau \,\alpha (\tau ;x).\end{eqnarray} \tag{ 140 }$$

The boundary terms of the integration by parts vanished because we assume that all particles are eventually absorbed (in particular, we assume ${{\lim}_{\tau \to \infty}}\tau \alpha (\tau ;x)=0$). The definition of the probability $\alpha (\tau ;x)$ in (139) implies that it fulfils the same forward-time flow equation as the marginalized distribution $|\,p(\tau,n|0){{\rangle}_{x}}$. Since ${{\lim}_{\tau \to \infty}}\alpha (\tau ;x)=0$, the mean extinction time (140) therefore obeys

$$\begin{eqnarray}&&-{{\mathcal{Q}}^{\dagger}}\left(x,{{\partial}_{x}}\right){{\langle \tau \rangle}_{x}}=\alpha (0;x)=1-{{\text{e}}^{-x}}.\end{eqnarray} \tag{ 141 }$$

The last equality follows from the fact that a certain fraction of all particles, namely ${{\text{e}}^{-x}}$, has already been in the absorbing state initially. The above derivation readily extends to higher moments of the mean extinction time.

For the linear decay process $A\to \varnothing$ with decay rate coefficient μ, the equation (141) for the mean extinction time reads

$$\begin{eqnarray}&&\mu {{\partial}_{x}}{{\langle \tau \rangle}_{x}}=\left(1-{{\text{e}}^{-x}}\right)/x.\end{eqnarray} \tag{ 142 }$$

This equation implies that the mean extinction time $\mu {{\langle \tau \rangle}_{x}}$ increases logarithmically with the particles' mean initial distance x from the absorbing state (see figure 5). The explicit solution of the equation is $\mu {{\langle \tau \rangle}_{x}}=\ln x+\gamma -\text{Ei}(-x)$ with Euler's constant γ and the exponential integral $\text{Ei}$ (6.2.6 in [335]; the exponential integral ensures that ${{\langle \tau \rangle}_{0}}=0$). The application of the functional $\langle {{n}_{0}}|\,f={{\left({{\partial}_{x}}\,+\,1\right)}^{{{n}_{0}}}}f(x){{|}_{x=0}}$ to the solution returns the mean extinction time of particles whose initial number is not Poisson distributed but that is fixed to some value ${{n}_{0}}\geqslant 0$. The corresponding mean extinction time $\mu {{\langle \tau \rangle}_{{{n}_{0}}}}$ is given by the harmonic number ${{H}_{{{n}_{0}}}}:={\sum}_{i=1}^{{{n}_{0}}}\frac{1}{i}$. One can also infer this result directly from the backward master equation using the methods discussed in [110, 215].

Zoom In Zoom Out Reset image size

Drummond et al have extended the above computation to chemical reactions with at most two reactants and two products [374]. For these reactions, the flow equation (138) has the mathematical form of a (forward) Fokker–Planck equation.

Both the probability generating function (45) and the marginalized distribution (110) were defined by summing the conditional probability distribution over a set of basis functions, either at the final or at the initial time. In addition, one can define the 'probability generating functional'

$$\begin{eqnarray}&&\langle g\left(\tau |{{t}_{0}},{{n}_{0}}\right)|\,:=\underset{n}{\sum}\,\langle n|\,p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)\end{eqnarray} \tag{ 143 }$$

by summing the conditional distribution over the set $\left\{\langle n|\right\}$ of basis functionals. The definition of the generating functional is meant with respect to test functions that can be expanded in the basis functions $\left\{|n\rangle \right\}$. As before, we assume that the basis functions and functionals constitute a (C)omplete and (O)rthogonal basis. The probability generating functional 'generates' probabilities in the sense that

$$\begin{eqnarray}&&p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)=\langle g\left(\tau |{{t}_{0}},{{n}_{0}}\right)|n\rangle.\end{eqnarray} \tag{ 144 }$$

In the next section, we show how this inverse transformation reduces to the Poisson representation of Gardiner and Chaturvedi [32, 33] upon choosing the basis function $|n\rangle$ as a Poisson distribution.

The derivation of the (functional) flow equation of $\langle g|$ proceeds analogously to the derivations in sections 2.1 and 3.1. Differentiation of its definition (143) with respect to τ and use of the forward master equation ${{\partial}_{\tau}}p\left(\tau |{{t}_{0}}\right)={{Q}_{\tau}}p\left(\tau |{{t}_{0}}\right)$ result in

$$\begin{eqnarray}&&{{\partial}_{\tau}}\langle g|\,=\underset{n}{\sum}\,\left({{\partial}_{\tau}}\langle n|\,+\underset{m}{\sum}\,Q_{\tau}^{\intercal}(n,m)\langle m|\right)\,p(\tau,n|\centerdot ).\end{eqnarray} \tag{ 145 }$$

This equation can be cast into a linear PDE by using the basis (E)volution operator ${{\mathcal{E}}_{\tau}}$ and the adjoint transition operator $\mathcal{Q}_{\tau}^{\dagger}$ (condition (${{\text{Q}}^{\intercal}}$)). In particular, since the Kronecker delta $\langle n|m\rangle ={{\delta}_{n,m}}$ is independent of time, it holds that $\langle n|{{\partial}_{\tau}}|m\rangle =\left(-{{\partial}_{\tau}}\langle n|\right)|m\rangle$. Consequently, the basis evolution operator in condition (E) fulfils

$$\begin{eqnarray}&&\langle n|{{\mathcal{E}}_{\tau}}=-{{\partial}_{\tau}}\langle n|\end{eqnarray} \tag{ 146 }$$

with respect to functions that can be expanded in the basis functions. Moreover, the (O)rthogonality and (C)ompleteness of the basis imply that the adjoint transition operator in condition (${{\text{Q}}^{\intercal}}$) acts on basis functionals as

$$\begin{eqnarray}&&\langle n|\mathcal{Q}_{\tau}^{\dagger}=\underset{m}{\sum}\,Q_{\tau}^{\intercal}(n,m)\langle m|.\end{eqnarray} \tag{ 147 }$$

Both of the above expressions hold for all values of n. If the two differential operators ${{\mathcal{E}}_{\tau}}$ and $\mathcal{Q}_{\tau}^{\dagger}$ exist, the generating functional $\langle g\left(\tau |{{t}_{0}},{{n}_{0}}\right)|$ obeys the functional flow equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}\langle g|\,=\langle g|\left(-{{\mathcal{E}}_{\tau}}+\mathcal{Q}_{\tau}^{\dagger}\right)=\langle g|\tilde{{\mathcal{Q}}}\,_{\tau}^{\dagger},\end{eqnarray} \tag{ 148 }$$

with initial condition $\langle g\left({{t}_{0}}|{{t}_{0}},{{n}_{0}}\right)|=\langle {{n}_{0}}|$. Note that this flow equation employs the same operator as the flow equation (116) of the marginalized distribution. Both equations can be used to derive the 'backward' path integral representation considered in section 4.

As a side note, let us remark that the flow equation ${{\partial}_{\tau}}|g\rangle =\left({{\mathcal{E}}_{\tau}}+{{\mathcal{Q}}_{\tau}}\right)|g\rangle$ of the generating function in (54) also admits a functional counterpart. In particular, the series

$$\begin{eqnarray}&&\langle \,p(t,n|\tau )|\,:=\underset{{{n}_{0}}}{\sum}\,p\left(t,n|\tau,{{n}_{0}}\right)\langle {{n}_{0}}|\end{eqnarray} \tag{ 149 }$$

obeys the flow equation

$$\begin{eqnarray}&&{{\partial}_{-\tau}}\langle \,p|=\langle \,p|\left({{\mathcal{E}}_{\tau}}+{{\mathcal{Q}}_{\tau}}\right)={{\tilde{\mathcal{Q}}}_{\tau}}\langle \,p|,\end{eqnarray} \tag{ 150 }$$

with final value $\langle \,p(t,n|t)|\,=\langle n|$. The corresponding inverse transformation reads $p\left(t,n|\tau,{{n}_{0}}\right)=\langle \,p(t,n|\tau )|{{n}_{0}}\rangle$. Both the flow equation obeyed by the generating function and the above flow equation can be used to derive the 'forward' path integral representation in section 5. Further uses of the series (149) remain to be explored.

Assuming that the action of the generating functional $\langle g|$ on a function f can be expressed in terms of an integral kernel (also called g) as

$$\begin{eqnarray}&&\langle g\left(\tau |{{t}_{0}},{{n}_{0}}\right)|\,f={\int}_{-\infty}^{\infty}\text{d}x\,g\left(\tau ;x|{{t}_{0}},{{n}_{0}}\right)f(x),\end{eqnarray} \tag{ 151 }$$

the insertion of the Poisson basis function $|n{{\rangle}_{x}}=\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}$ into the inverse transformation (144) results in

$$\begin{eqnarray}&&p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)={\int}_{-\infty}^{\infty}\text{d}x\,g\left(\tau ;x|{{t}_{0}},{{n}_{0}}\right)\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}.\end{eqnarray} \tag{ 152 }$$

A representation of the probability distribution of this form is called a 'Poisson representation' and was first proposed by Gardiner and Chaturvedi [32, 33]. Since the integration in (152) proceeds along the real line, the above representation is referred to as a 'real' Poisson representation [375]. Although the use of a real variable may seem convenient, its use typically results in the kernel being a 'generalized function', i.e. a distribution. For example, the initial condition $p\left({{t}_{0}},n|{{t}_{0}},{{n}_{0}}\right)={{\delta}_{n,{{n}_{0}}}}$ is recovered for the integral kernel $g\left({{t}_{0}};x|{{t}_{0}},{{n}_{0}}\right)=\delta (x){{\left({{\partial}_{x}}+1\right)}^{{{n}_{0}}}}$. By the definition

$$\begin{eqnarray}&&{\int}_{-\infty}^{\infty}\text{d}x\,\left(\partial _{x}^{j}\delta (x)\right)f(x):={\int}_{-\infty}^{\infty}\text{d}x\,\delta (x){{\left(-{{\partial}_{x}}\right)}^{j}}f(x)\end{eqnarray} \tag{ 153 }$$

of distributional derivatives (1.16.12 in [335], $j\in {{\mathbb{N}}_{0}}$), this kernel can also be written as $\left[{{\left(1-{{\partial}_{x}}\right)}^{{{n}_{0}}}}\delta (x)\right]$. The integral kernel $\delta \left(x-{{x}_{0}}\right)$ instead results for a Poisson distribution with mean x₀. In [32], Gardiner and Chaturvedi used the real Poisson representation to calculate steady-state probability distributions of various elementary reactions. Furthermore, the real Poisson representation was employed by Elderfield [376] to derive a stochastic path integral representation. Droz and McKane [377] later argued that this representation is equivalent to a path integral representation based on Doi's Fock space algebra. Our discussion of the backward and forward path integral representations in sections 4 and 5 clarifies the similarities and differences between these two approaches.

To circumvent the use of generalized functions, various alternatives to the real Poisson representation have been proposed: the 'complex' and 'positive' Poisson representations [110, 375] as well as the 'gauge' Poisson representation [378]. The former two representations as well as the real representation are discussed in the book of Gardiner [110]. The complex Poisson representation is obtained by continuing the Poisson basis function in (152) into the complex domain and performing the integration around a closed path $\mathcal{C}$ around 0 (once in counter-clockwise direction). Upon using Cauchy's differentiation formula to redefine the basis functional as

$$\begin{eqnarray}&&\langle n|\,f:=\partial _{x}^{n}{{\text{e}}^{x}}f(x){{|}_{x=0}}={{{\oint}^{}}_{\mathcal{C}}}\text{d}x\frac{n!}{2\pi \text{i}}\frac{{{\text{e}}^{x}}}{{{x}^{n+1}}}f(x),\end{eqnarray} \tag{ 154 }$$

it becomes apparent that the initial condition $p\left({{t}_{0}},n|{{t}_{0}},{{n}_{0}}\right)={{\delta}_{n,{{n}_{0}}}}$ is then recovered for the kernel $g\left({{t}_{0}};x|{{t}_{0}},{{n}_{0}}\right)=\frac{{{n}_{0}}!}{2\pi \text{i}}\frac{{{\text{e}}^{x}}}{{{x}^{{{n}_{0}}+1}}}$. As the interpretation of this kernel is also not straightforward, we refrain from calling it a 'quasi-probability distribution' [32].

Let us exemplify the flow equation obeyed by the integral kernel for the reaction $k\,A\to l\,A$ with rate coefficient ${{\gamma}_{\tau}}$. The flow equation can be inferred from the corresponding flow equation (148) of the generating functional and the kernel's definition in (151). Since we employ the Poisson basis function, the results from section 3.2.2 imply that the adjoint transition operator (122) reads

$$\begin{eqnarray}&&\mathcal{Q}_{\tau}^{\dagger}={{\gamma}_{\tau}}{{x}^{k}}\left({{\left({{\partial}_{x}}+1\right)}^{l}}-{{\left({{\partial}_{x}}+1\right)}^{k}}\right).\end{eqnarray} \tag{ 155 }$$

Consequently, one finds that the integral kernel $g(\tau ;x|\centerdot )$ obeys the flow equation

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{\tau}}g & ={{\mathcal{Q}}_{\tau}}\left(x,{{\partial}_{x}}\right)g \end{array}\end{eqnarray} \tag{ 156 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & ={{\gamma}_{\tau}}\left({{\left(1-{{\partial}_{x}}\right)}^{l}}-{{\left(1-{{\partial}_{x}}\right)}^{k}}\right){{x}^{k}}g. \end{array}\end{eqnarray} \tag{ 157 }$$

To arrive at this equation, we performed repeated integrations by parts while ignoring any potential boundary terms (see the definition of the adjoint operator in (115)). The importance of boundary terms is discussed in [378].

Thus far, most studies employing the Poisson representation have focused on networks of bimolecular reactions ${\sum}_{j}{{k}_{j}}\,{{A}_{j}}\to {\sum}_{j}{{l}_{j}}\,{{A}_{j}}$ with ${\sum}_{j}{{k}_{j}}\leqslant 2$ and ${\sum}_{j}{{l}_{j}}\leqslant 2$. For these networks, the flow equation (156) assumes the mathematical form of a forward Fokker–Planck equation with the derivatives being of at most second order (the corresponding diffusion coefficient may be negative). It has been attempted to map the resulting equation to an Itô SDE [32], but it should be explored whether this procedure is supported by the Feynman–Kac formula in appendix A. The numerical integration of an SDE with potentially negative diffusion coefficient was attempted for the bi-directional reaction $2\,A\leftrightharpoons \varnothing$ in [379]. The value of the integration has remained inconclusive. Recently, an exponential ansatz for the integral kernel $g(\tau ;x|\centerdot )$ has been considered to approximate its flow equation in the limit of weak noise [380] (see section 2.3.2). Moreover, a gauge Poisson representation was recently employed in a study of the coagulation reaction $2\,A\to A$ [381].

In the previous section, we formulated general conditions under which the master equation can be transformed into a partial differential equation obeyed by the probability generating function (45). In the present section, we complemented this flow equation by a backward-time flow equation obeyed by the marginalized distribution (110) and by a functional flow equation obeyed by the probability generating functional (143). Whereas the marginalized distribution

$$\begin{eqnarray}&&|\,p(t,n|\tau )\rangle =\underset{{{n}_{0}}}{\sum}\,p\left(t,n|\tau,{{n}_{0}}\right)|{{n}_{0}}\rangle \end{eqnarray} \tag{ 158 }$$

was defined as the sum of the conditional probability distribution over a set of basis functions, the generating functional

$$\begin{eqnarray}&&\langle g\left(\tau |{{t}_{0}},{{n}_{0}}\right)|\,=\underset{n}{\sum}\,\langle n|\,p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)\end{eqnarray} \tag{ 159 }$$

was defined as the sum of the distribution over a set of basis functionals. In section 3.2, we introduced (O)rthogonal and (C)omplete basis functions and functionals for the study of different stochastic processes, including the Poisson basis for the study of chemical reactions ($\,|n{{\rangle}_{x}}=\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}$ and $\langle m|\,f={{\left({{\partial}_{x}}+1\right)}^{m}}f(x){{|}_{x=0}}$). Provided that there also exist a basis (E)volution operator ${{\mathcal{E}}_{\tau}}$ and an adjoint transition operator $\mathcal{Q}_{\tau}^{\dagger}$ (condition (${{\text{Q}}^{\intercal}}$)), the marginalized distribution obeys the backward-time flow equation

$$\begin{eqnarray}&&{{\partial}_{-\tau}}|\,p\rangle =\left(-{{\mathcal{E}}_{\tau}}+\mathcal{Q}_{\tau}^{\dagger}\right)|\,p\rangle =\tilde{\mathcal{Q}}_{\tau}^{\dagger}|\,p\rangle \end{eqnarray} \tag{ 160 }$$

with final condition $|\,p(t,n|t)\rangle =\,|n\rangle$. Moreover, the probability generating functional (143) then obeys the functional flow equation (148). The inverse transformation $p\left(\tau,n|{{t}_{0}},{{n}_{0}}\right)=\langle g\left(\tau |{{t}_{0}},{{n}_{0}}\right)|n\rangle$ of this functional generalizes the Poisson representation of Gardiner and Chaturvedi [32, 33] as shown in section 3.5. In section 3.3, we showed how the flow equation (160) obeyed by the marginalized distribution can be used to compute mean extinction times. Furthermore, the equation will prove useful in the derivation of path integral representations of the master equation and of averaged observables in sections 4 and 6. Future studies could explore whether the flow equation obeyed by the marginalized distribution can be evaluated in terms of WKB approximations or spectral methods.

We first derive the backward path integral representation (162) because its application to actual processes is more intuitive than the application of the forward path integral representation. For this purpose, we consider the flow equation ${{\partial}_{-\tau}}|\,p(\centerdot |\tau )\rangle =\tilde{\mathcal{Q}}_{\tau}^{\dagger}|\,p(\centerdot |\tau )\rangle$ in (116) obeyed by the marginalized distribution. Its final condition reads $|\,p(t,n|t)\rangle =\,|n\rangle$. The conditional probability distribution is recovered from the marginalized distribution via $p\left(t,n|{{t}_{0}},{{n}_{0}}\right)=\langle {{n}_{0}}|\,p\left(t,n|{{t}_{0}}\right)\rangle$. As the first step, we split the time interval [t₀,t] into N pieces ${{t}_{0}}\leqslant {{t}_{1}}\leqslant \cdots \leqslant {{t}_{N}}:=t$ of length $\Delta t:=\left(t-{{t}_{0}}\right)/N\ll 1$. Over the time interval $\left[{{t}_{0}},{{t}_{1}}\right]$, the flow equation is then solved by ⁹

$$\begin{eqnarray}&&|\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{{{x}_{0}}}}=\mathcal{L}_{{{t}_{0}}}^{\dagger}\left({{x}_{0}},{{\partial}_{{{x}_{0}}}}\right)|\,p\left(t,n|{{t}_{1}}\right){{\rangle}_{{{x}_{0}}}}\end{eqnarray} \tag{ 163 }$$

with the generator $\mathcal{L}_{\tau}^{\dagger}:=1+\tilde{\mathcal{Q}}_{\tau}^{\dagger} \Delta t+\text{O}\left({{(\Delta t)}^{2}}\right)$. Alternatively, the following derivation could be performed by solving the flow equation (148) obeyed by the generating functional over the time interval [t_N−1, t] as (see figure 6)

$$\begin{eqnarray}&&\langle g\left(t|{{t}_{0}},{{n}_{0}}\right)|=\langle g\left({{t}_{N-1}}|{{t}_{0}},{{n}_{0}}\right)|\mathcal{L}_{{{t}_{N-1}}}^{\dagger}\left({{x}_{N-1}},{{\partial}_{{{x}_{N-1}}}}\right).\,\end{eqnarray} \tag{ 164 }$$

Zoom In Zoom Out Reset image size

Here, the generating functional $\langle g|$ acts on functions f(x_N−1), meets the initial condition $\langle g\left({{t}_{0}}|{{t}_{0}},{{n}_{0}}\right)|\,=\langle {{n}_{0}}|$, and is transformed back into the conditional probability distribution via $p\left(t,n|{{t}_{0}},{{n}_{0}}\right)=\langle g\left({{t}_{0}}|{{t}_{0}},{{n}_{0}}\right)|n\rangle$. Both of the above approaches result in the same path integral representation. In the following, we use the equation (163) for this purpose.

As the next step of the derivation, we insert the integral representation of a Dirac delta between the generator ${{\mathcal{L}}^{\dagger}}$ and the marginalized distribution $|\,p\rangle$ in (163), turning the right-hand side of this expression into

$$\begin{eqnarray}&&\mathcal{L}_{{{t}_{0}}}^{\dagger}\left({{x}_{0}},{{\partial}_{{{x}_{0}}}}\right){\int}_{{{\mathbb{R}}^{2}}}\frac{\text{d}{{x}_{1}}\text{d}{{q}_{1}}}{2\pi}\,{{\text{e}}^{-\text{i}{{q}_{1}}\left({{x}_{1}}-{{x}_{0}}\right)}}|\,p\left(t,n|{{t}_{1}}\right){{\rangle}_{{{x}_{1}}}}.\end{eqnarray} \tag{ 165 }$$

The integrations over x₁ and q₁ are both performed along the real line from $-\infty$ to $+\infty$. If the adjoint transition operator $\tilde{\mathcal{Q}}_{\tau}^{\dagger}$, and thus also $\mathcal{L}_{\tau}^{\dagger}$, are normal-ordered power series, we can replace ${{\partial}_{{{x}_{0}}}}$ by $\text{i}{{q}_{1}}$ in the above expression and pull $\mathcal{L}_{{{t}_{0}}}^{\dagger}$ to the right of the exponential. Making use of the final condition $|\,p(t,n|t){{\rangle}_{{{x}_{N}}}}=\,|n{{\rangle}_{t,{{x}_{N}}}}$, the above steps can be repeated until (165) reads

$$\begin{eqnarray}&&{\int}^{}-\,_{1}^{N}\,\left(\underset{j=1}{\overset{N}{\prod}}\,{{\text{e}}^{-\text{i}{{q}_{j}}\left({{x}_{j}}-{{x}_{j-1}}\right)}}\mathcal{L}_{{{t}_{j-1}}}^{\dagger}\left({{x}_{j-1}},\text{i}{{q}_{j}}\right)\right)|n{{\rangle}_{t,{{x}_{N}}}},\end{eqnarray} \tag{ 166 }$$

with the abbreviation

$$\begin{eqnarray}&&{\int}^{}-_{k}^{\,l}:=\underset{j=k}{\overset{l}{\prod}}\,{{{\int}^{}}_{{{\mathbb{R}}^{2}}}}\frac{\text{d}{{x}_{j}}\text{d}{{q}_{j}}}{2\pi}.\end{eqnarray} \tag{ 167 }$$

To proceed, we now replace $\mathcal{L}_{\tau}^{\dagger}=1+\tilde{\mathcal{Q}}_{\tau}^{\dagger} \Delta t+\text{O}\left({{(\Delta t)}^{2}}\right)$ by the exponential ${{\text{e}}^{\tilde{\mathcal{Q}}_{\tau}^{\dagger} \Delta t}}$ (for brevity, we drop the correction term in the following). Note that the exponential ${{\text{e}}^{\tilde{\mathcal{Q}}_{\tau}^{\dagger} \Delta t}}$ does not involve differential operators because those were all replaced by the new variables $\text{i}{{q}_{j}}$. Whether the exponentiation can be made mathematically rigorous is an open question and may possibly be answered positively only for a restricted class of stochastic processes. Upon performing the exponentiation, one obtains the following discrete-time path integral representation of the marginalized distribution:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} |\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{{{x}_{0}}}}={\rlap{-} \int}_{1}^{\,N}\,{{\text{e}}^{-\mathcal{S}_{N}^{\dagger}}}\,|n{{\rangle}_{t,{{x}_{N}}}}\,\text{with} \end{array}\end{eqnarray} \tag{ 168 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{S}_{N}^{\dagger}:=\underset{j=1}{\overset{N}{\sum}}\, \Delta t\left(\text{i}{{q}_{j}}\frac{{{x}_{j}}-{{x}_{j-1}}}{ \Delta t}-\tilde{\mathcal{Q}}_{{{t}_{j-1}}}^{\dagger}\left({{x}_{j-1}},\text{i}{{q}_{j}}\right)\right). \end{array}\end{eqnarray} \tag{ 169 }$$

The final condition $|\,p(t,n|t){{\rangle}_{x}}=\,|n{{\rangle}_{t,x}}$ is trivially fulfilled because only N  =  0 time slices fit between t and t. The object $\mathcal{S}_{N}^{\dagger}$ is called an 'action'. The exponential factor ${{\text{e}}^{-\mathcal{S}_{N}^{\dagger}}}$ weighs the contribution of every path $\left({{x}_{1}},{{q}_{1}}\right)\to \left({{x}_{2}},{{q}_{2}}\right)\to \cdots \to \left({{x}_{N}},{{q}_{N}}\right)$.

As the final step of the derivation, we take the continuous-time limit $N\to \infty$, or $\Delta t\to 0$, at least formally. The following continuous-time expressions effectively serve as abbreviations of the above discretization scheme with the identifications $x\left({{t}_{0}}+j \Delta t\right):={{x}_{j}}$ and $q\left({{t}_{0}}+j \Delta t\right):={{q}_{j}}$. Further comments on the discretization scheme are provided in section 4.4. Combined with the inverse transformation (111), we thus arrive at the following backward path integral representation of the conditional probability distribution:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} p\left(t,n|{{t}_{0}},{{n}_{0}}\right)={{\langle {{n}_{0}}{{|}_{{{t}_{0}},x\left({{t}_{0}}\right)}}|\,p\left(t,n|{{t}_{0}}\right)\rangle}_{x\left({{t}_{0}}\right)}} \end{array}\end{eqnarray} \tag{ 170 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{with}\,|\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{x\left({{t}_{0}}\right)}}={\rlap{-} \int}_{({{t}_{0}}}^{t]}\,{{\text{e}}^{-{{\mathcal{S}}^{\dagger}}}}\,|n{{\rangle}_{t,x(t)}} \end{array}\end{eqnarray} \tag{ 171 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{and}~{{\mathcal{S}}^{\dagger}}:={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\left[\text{i}q{{\partial}_{\tau}}x-\tilde{\mathcal{Q}}_{\tau}^{\dagger}\left(x,\text{i}q\right)\right]. \end{array}\end{eqnarray} \tag{ 172 }$$

Here we included x(t₀) as an argument of the functional $\langle {{n}_{0}}{{|}_{{{t}_{0}},x\left({{t}_{0}}\right)}}$ to express that this functional acts on $x(\tau )$ only for $\tau ={{t}_{0}}$. Moreover, we defined ${\rlap{-} \int}_{({{t}_{0}}}^{t]}:={{\lim}_{N\to \infty}}{\int}^{}-\,_{1}^{N}$ to indicate that the path integral involves integrations over x(t) and q(t), but not over x(t₀) and q(t₀).

The evaluation of the above path integral representation for an explicit process involves two steps. First, the path integral (171) provides the marginalized distribution $|\,p\left(t,n|{{t}_{0}}\right)\rangle$. Second, this marginalized distribution is mapped to the conditional probability distribution $p\left(t,n|{{t}_{0}},{{n}_{0}}\right)$ by the action of the functional $\langle {{n}_{0}}|$.

Let us demonstrate the above two steps for the bi-directional reaction $\varnothing \rightleftharpoons A$. The stationary distribution of this process was already derived in the introduction to section 2. Both the growth rate coefficient γ and the decay rate coefficient μ shall be time-independent, but this assumption is easily relaxed. As the first step, the backward path integral (170) requires us to choose an appropriate basis. The time-independent Poisson basis with the basis function $|n{{\rangle}_{x}}=\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}$ proves to be convenient for this purpose (see section 3.2.2). The adjoint transition operator follows from (122), (128), and (129) as ${{\tilde{\mathcal{Q}}}^{\dagger}}\left(x,{{\partial}_{x}}\right)=(\gamma -\mu x){{\partial}_{x}}$. Insertion of this operator into the action (172) results in

$$\begin{eqnarray}&&{{\mathcal{S}}^{\dagger}}={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\text{i}q\,\left[{{\partial}_{\tau}}x-(\gamma -\mu x)\right]\end{eqnarray} \tag{ 173 }$$

The integration over the path $q(\tau )$ from $\tau ={{t}_{0}}$ to $\tau =t$ is performed most easily in the discrete-time approximation. The marginalized distribution $|\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{x\left({{t}_{0}}\right)}}$ with $x\left({{t}_{0}}\right)\equiv {{x}_{0}}$ thereby follows as the following product of Dirac deltas:

$$\begin{eqnarray}&&\left(\underset{j=1}{\overset{N}{\prod}}\,{{{\int}^{}}_{\mathbb{R}}}\text{d}{{x}_{j}}\,\delta \left({{x}_{j}}-{{x}_{j-1}}-\left(\gamma -\mu \,{{x}_{j-1}}\right) \Delta t\right)\right)|n{{\rangle}_{t,{{x}_{N}}}}.\end{eqnarray}$$

The function $|n{{\rangle}_{t,{{x}_{N}}}}$ is also integrated over. Upon taking the continuous-time $\Delta t\to 0$ and performing the integration over the path $x(\tau )$, one finds that the marginalized distribution is given by the Poisson distribution

$$\begin{eqnarray}&&|\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{x\left({{t}_{0}}\right)}}=\,|n{{\rangle}_{t,x(t)}}=\frac{x{{(t)}^{n}}{{\text{e}}^{-x(t)}}}{n!},\end{eqnarray} \tag{ 174 }$$

with $x(\tau )$ solving ${{\partial}_{\tau}}x=\gamma -\mu x$. This equation coincides with the rate equation of the process $\varnothing \rightleftharpoons A$. The solution $x(\tau )={{\text{e}}^{-\mu \left(\tau -{{t}_{0}}\right)}}\left(x\left({{t}_{0}}\right)-\frac{\gamma}{\mu}\right)+\frac{\gamma}{\mu}$ of the rate equation specifies the mean and the variance of the above Poisson distribution. For $\mu >0$, both of them converge to the asymptotic value $\gamma /\mu$. In the next section and in section 4.5, we generalize these results to more general processes.

The marginalized distribution (174) does not only solve the master equation of the reaction $\varnothing \rightleftharpoons A$, but it also establishes the link between the path integral variable x and the moments of the particle number n. In particular, the mean particle number evaluates to $\langle n\rangle (t)=x(t)$, while higher order moments can determined via $\langle {{n}^{k}}\rangle (t)={\sum}_{l=0}^{k}\left\{\begin{array}{c} k \\ l \end{array}\right\}x{{(t)}^{l}}$. The curly braces denote a Stirling number of the second kind (see section 26.8 in [335]).

The conditional probability distribution can be calculated from the marginalized distribution by applying the functional $\langle {{n}_{0}}|\,f={{\left({{\partial}_{x\left({{t}_{0}}\right)}}+1\right)}^{{{n}_{0}}}}f\left(x\left({{t}_{0}}\right)\right){{|}_{x\left({{t}_{0}}\right)=0}}$ to the latter (see (126) and (129)). In the limit of a vanishing decay rate ($\mu \to 0$), for which $x(\tau )=x\left({{t}_{0}}\right)+\gamma \left(\tau -{{t}_{0}}\right)$, we thereby recover the shifted Poisson distribution (80), i.e.

$$\begin{eqnarray}&&p\left(t,n|{{t}_{0}},{{n}_{0}}\right)=\frac{{{\text{e}}^{-\gamma \left(t-{{t}_{0}}\right)}}{{\left(\gamma \left(t-{{t}_{0}}\right)\right)}^{n-{{n}_{0}}}}}{\left(n-{{n}_{0}}\right)!}{{ \Theta }_{n-{{n}_{0}}}}.\end{eqnarray} \tag{ 175 }$$

Thus, the mean and variance of the conditional probability distribution grow linearly with time. If the particles decay but are not replenished ($\mu >0$ but $\gamma =0$), we instead recover the Binomial distribution (134).

We now show how the backward path integral representation (171) of the marginalized distribution can be expressed in terms of an average over the paths of an Itô stochastic differential equation (SDE). The resulting expression bears similarities with the Feynman–Kac formula (3), especially when the adjoint transition operator $\tilde{\mathcal{Q}}_{\tau}^{\dagger}\left(x,{{\partial}_{x}}\right)$ of the stochastic process under consideration is quadratic in ${{\partial}_{x}}$. In the general case, however, functional derivatives act on the average over paths. These derivatives can, for example, be evaluated in terms of perturbation expansions as we demonstrate in section 4.5.4. The procedure outlined below serves as a general starting point for the exact or approximate evaluation of the backward path integral (171).

In the following, we consider a stochastic process whose adjoint transition operator has the generic form

$$\begin{eqnarray}&&\tilde{{\mathcal{Q}}}\,_{\tau}^{\dagger}\left(x,{{\partial}_{x}}\right)={{\alpha}_{\tau}}(x){{\partial}_{x}}\,+\,\frac{1}{2}{{\beta}_{\tau}}(x)\partial _{x}^{2}\,+\,\mathcal{P}_{\tau}^{\dagger}\left(x,{{\partial}_{x}}\right).\end{eqnarray} \tag{ 176 }$$

As before, we call ${{\alpha}_{\tau}}$ a drift and ${{\beta}_{\tau}}$ a diffusion coefficient. The object $\mathcal{P}_{\tau}^{\dagger}$ is referred to as the (adjoint) perturbation operator, or simply as the perturbation. The perturbation operator absorbs all the terms of higher order in ${{\partial}_{x}}$ and possibly also terms of lower order. Thus, the above form of $\tilde{\mathcal{Q}}_{\tau}^{\dagger}$ is not unique and ${{\alpha}_{\tau}}$, ${{\beta}_{\tau}}$, and $\mathcal{P}_{\tau}^{\dagger}$ should be chosen so that the evaluation of the expressions below becomes as simple as possible. If the perturbation operator $\mathcal{P}_{\tau}^{\dagger}$ is zero, those expressions simplify considerably.

As the first step, let us rewrite the backward path integral representation (171) of the marginalized distribution as

$$\begin{eqnarray}&&|\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{x\left({{t}_{0}}\right)}}={\int}_{{{\mathbb{R}}^{2}}}\frac{\text{d}{{x}_{N}}\text{d}{{q}_{N}}}{2\pi}{{\text{e}}^{-\text{i}{{q}_{N}}{{x}_{N}}}}{{\mathcal{Z}}_{0,0}}|n{{\rangle}_{t,{{x}_{N}}}},\end{eqnarray} \tag{ 177 }$$

where ${{\mathcal{Z}}_{Q=0,X=0}}$ represents a value of the (Q,X)-generating functional

$$\begin{eqnarray}&&{{\mathcal{Z}}_{Q,X}}:={\rlap{-} \int}_{({{t}_{0}}}^{\,t)}{{\text{e}}^{\text{i}{{q}_{N}}{{x}_{N}}-{{\mathcal{S}}^{\dagger}}+{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\left[Q(\tau )x(\tau )+X(\tau )\text{i}q(\tau )\right]}}.\end{eqnarray} \tag{ 178 }$$

Hence, we have singled out the integrations over x_N and q_N before performing the continuous-time limit (see (168) and (171)). We call ${{\mathcal{Z}}_{Q,X}}$ a generating functional because a functional differentiation of ${{\mathcal{Z}}_{Q,X}}$ with respect to $Q(\tau )$ generates a factor $x(\tau )$, and a functional differentiation of ${{\mathcal{Z}}_{Q,X}}$ with respect to $X(\tau )$ generates a factor $\text{i}q(\tau )$.

Given the action (172) with the adjoint transition operator (176), the above properties of the (Q,X)-generating functional can now be used to rewrite this function as (see appendix D)

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{Z}}_{Q,X}}={{\text{e}}^{\text{i}{{q}_{N}}\frac{\delta}{\delta Q(t)}+{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\mathcal{P}_{\tau}^{\dagger}\left(\frac{\delta}{\delta Q(\tau )},\frac{\delta}{\delta X(\tau )}\right)}}\mathcal{Z}_{Q,X}^{0} \end{array}\end{eqnarray} \tag{ 179 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{with}\,\mathcal{Z}_{Q,X}^{0}:={{\langle \langle {{\text{e}}^{{\int}_{{{t}_{0}}}^{t}\text{d}\tau Q(\tau )x(\tau )}}\rangle \rangle}_{W}}. \end{array}\end{eqnarray} \tag{ 180 }$$

Here, ${{\langle \langle \centerdot \rangle \rangle}_{W}}$ represents the average over realizations of a Wiener process $W(\tau )$. This Wiener process influences the evolution of the path $x(\tau )$ through the Itô SDE

$$\begin{eqnarray}&&\text{d}x(\tau )=\left[{{\alpha}_{\tau}}(x)+X(\tau )\right]\text{d}\tau +\sqrt{{{\beta}_{\tau}}(x)}\text{d}W(\tau ).\end{eqnarray} \tag{ 181 }$$

The temporal evolution of $x(\tau )$ starts out from x(t₀), which is determined by the argument of the marginalized distribution (177). In section 4.5.4, we demonstrate how the marginalized distribution can be evaluated in terms of a perturbation expansion of the (Q,X)-generating functional (179) for a process of diffusing particles that are also decaying. In order to perform the perturbation expansion, let us already note that the value $x(\tau )$ of the path at time τ depends on the 'source' $X\left({{\tau}^{\prime}}\right)$ only for times ${{\tau}^{\prime}}<\tau$, and that $\frac{\delta}{\delta Q(\tau )}{\int}_{{{t}_{0}}}^{t}\text{d}{{\tau}^{\prime}}Q\left({{\tau}^{\prime}}\right)f\left({{\tau}^{\prime}}\right)=f(\tau )$ holds for all $\tau \in \left({{t}_{0}},t\right]$. Moreover, let us note that the (Q,X)-generating functional depends on q_N and x(t₀) but not on x_N . These properties follow from the derivation of the representation (179), which we outline in appendix D.

An important special case of the above representation is constituted by processes whose perturbation operator $\mathcal{P}_{\tau}^{\dagger}$ is zero. For Q  =  X  =  0, the generating functional (179) then simplifies to ${{\mathcal{Z}}_{0,0}}={{\text{e}}^{\text{i}{{q}_{N}}x(t)}}$ and the marginalized distribution (177) follows in terms of the Feynman–Kac like fomula

$$\begin{eqnarray}&&|\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{x\left({{t}_{0}}\right)}}=\langle {{\langle |n\rangle}_{t,x(t)}}\rangle {{\rangle}_{W}}.\end{eqnarray} \tag{ 182 }$$

Here, $x(\tau )$ solves the Itô SDE

$$\begin{eqnarray}&&\text{d}x(\tau )={{\alpha}_{\tau}}(x)\text{d}\tau +\sqrt{{{\beta}_{\tau}}(x)}\text{d}W(\tau )\end{eqnarray} \tag{ 183 }$$

with the initial value x(t₀). For the Poisson basis function $|n{{\rangle}_{x}}=\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}$, the above representation of the marginalized distribution coincides with our earlier result (136), which we encountered in the discussion of the coagulation reaction $2\,A\to A$. For the bi-directional reaction $\varnothing \rightleftharpoons A$ from the previous section with adjoint transition operator ${{\tilde{\mathcal{Q}}}^{\dagger}}\left(x,{{\partial}_{x}}\right)=(\gamma -\mu x){{\partial}_{x}}$, the Itô SDE (183) simplifies to the deterministic rate equation ${{\partial}_{\tau}}x=\gamma -\mu x$. Thus, no averaging is required to evaluate the marginalized distribution (182) and one recovers our earlier solution (174).

Our above discussion only applies to processes in well-mixed environments with a single type of particles. A generalization of the results to processes with multiple types of particles or with spatial degrees of freedom is straightforward. A spatial process will be discussed in section 4.5. In a multivariate generalization of the above procedure, the adjoint transition operator (176) includes a vector-valued drift coefficient $\boldsymbol{\alpha }_{{\tau}}\left(\boldsymbol{x}\right)$ and a diffusion matrix ${{\beta}_{\tau}}\left(\boldsymbol{x}\right)$. Moreover, the derivation of the corresponding generating functional (179) requires that there exists a matrix $\sqrt{{{\beta}_{\tau}}}:={{\gamma}_{\tau}}$ fulfilling ${{\gamma}_{\tau}}\gamma _{\tau}^{\intercal}={{\beta}_{\tau}}$. If the diffusion matrix ${{\beta}_{\tau}}$ is positive-semidefinite, its positive-semidefinite and symmetric square root $\sqrt{{{\beta}_{\tau}}}$ can be determined via diagonalization [111].

In section 4.5, we exemplify the use of our above results for various well-mixed and spatial processes. But before, let us show how the results from the previous sections can be used to derive a path integral representation for processes with continuous state spaces.

The transition rate ${{w}_{\tau}}(m,n)$ denotes the rate at which probability flows from a state n to a state m, or, considering an individual sample path, the rate at which particles jump from state n to state m. Thus, ${{\kappa}_{\tau}}(\Delta n,n):={{w}_{\tau}}(n+ \Delta n,n)$ denotes the rate at which the state n is left via jumps of size $\Delta n$ (with $n\in {{\mathbb{N}}_{0}}$ and $\Delta n\in \mathbb{Z}$). Thus far, we only considered jumps between the states of a discrete state space. But in the following, we derive a path integral representation for processes having a continuous state space.

4.4.1. Processes with continuous state spaces.

We consider a process whose state is characterized by a continuous variable $x\in {{\mathbb{R}}_{\geqslant 0}}$. The change from the letter n to the letter x is purely notational and emphasizes that the state space is now continuous (the change also highlights a formal similarity between the linear PDEs discussed so far and the ones derived below). The conditional probability distribution $p\left(\tau,x|{{t}_{0}},{{x}_{0}}\right)$ describing the system shall be normalized as ${{{\int}^{}}_{\mathbb{R}}}\text{d}x\,p(\tau,x|\centerdot )=1$ and obey the master equation

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{\tau}}p(\tau,x|\centerdot )={{{\int}^{}}_{\mathbb{R}}}\text{d} \Delta x & \left[{{\kappa}_{\tau}}(\Delta x,x- \Delta x)p(\tau,x- \Delta x|\centerdot ) \right. \\ {} & \left. -{{\kappa}_{\tau}}(\Delta x,x)p(\tau,x|\centerdot )\right] \end{array}\end{eqnarray} \tag{ 184 }$$

with the initial condition $p\left({{t}_{0}},x|{{t}_{0}},{{x}_{0}}\right)=\delta \left(x-{{x}_{0}}\right)$. The structure of this master equation is equivalent to the structure of the master equation (11). Provided that the product ${{\kappa}_{\tau}}(\Delta x,x)p(\tau,x|\centerdot )$ is analytic in x, one can perform a Taylor expansion of the above master equation to obtain the (forward) Kramers–Moyal expansion [1, 103, 272]

$$\begin{eqnarray}&&{{\partial}_{\tau}}p(\tau,x|\centerdot )=\underset{m=1}{\overset{\infty}{\sum}}\,\frac{{{(-1)}^{m}}}{m!}\partial _{x}^{m}\left[\mathcal{M}_{\tau}^{(m)}(x)p(\tau,x|\centerdot )\right].\end{eqnarray} \tag{ 185 }$$

The 'jump moments' $\mathcal{M}_{\tau}^{(m)}$ are defined as

$$\begin{eqnarray}&&\mathcal{M}_{\tau}^{(m)}(x):={{{\int}^{}}_{\mathbb{R}}}\text{d} \Delta x\,{{(\Delta x)}^{m}}{{\kappa}_{\tau}}(\Delta x,x).\end{eqnarray} \tag{ 186 }$$

By Pawula's theorem [104, 383, 384], the positivity of the conditional probability distribution requires that the Kramers–Moyal expansion either stops at its first or second summand, or that it does not stop at all. If the expansion stops at its second summand, it assumes the form of a (forward) Fokker–Planck equation. The drift coefficient of this Fokker–Planck equation is given by $\mathcal{M}_{\tau}^{(1)}(x)$ and its non-negative diffusion coefficient by $\mathcal{M}_{\tau}^{(2)}(x)$. The sample paths of the Fokker–Planck equation are continuous, however, contradicting our earlier assumption of the process making discontinuous jumps in state space. For a jump process, the Kramers–Moyal expansion cannot stop. Nevertheless, a truncation of the Kramers–Moyal expansion at the level of a Fokker–Planck equation often provides a decent approximation of a process, provided that fluctuations cause only small relative changes of its state x.

The backward analogue of the master equation (184) reads

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\partial}_{-\tau}}p\left(t,x|\tau,{{x}_{0}}\right) \\ {} & ={\int}_{\mathbb{R}}\text{d} \Delta x\,\left[p\left(\centerdot |\tau,{{x}_{0}}+ \Delta x\right)-p\left(\centerdot |\tau,{{x}_{0}}\right)\right]{{\kappa}_{\tau}}\left(\Delta x,{{x}_{0}}\right), \end{array}\end{eqnarray} \tag{ 187 }$$

with the final condition $p\left(t,x|t,{{x}_{0}}\right)=\delta \left(x-{{x}_{0}}\right)$. Given the analyticity of $p\left(\centerdot |\tau,{{x}_{0}}\right)$ in x₀, the backward master equation can be rewritten in terms of the backward Kramers–Moyal expansion

$$\begin{eqnarray}&&{{\partial}_{-\tau}}p\left(\centerdot |\tau,{{x}_{0}}\right)=\tilde{\mathcal{Q}}_{\tau}^{\dagger}\left({{x}_{0}},{{\partial}_{{{x}_{0}}}}\right)p\left(\centerdot |\tau,{{x}_{0}}\right)\end{eqnarray} \tag{ 188 }$$

with the adjoint transition operator

$$\begin{eqnarray}&&\tilde{{\mathcal{Q}}}\,_{\tau}^{\dagger}\left({{x}_{0}},{{\partial}_{{{x}_{0}}}}\right):=\underset{m=1}{\overset{\infty}{\sum}}\,\frac{1}{m!}\mathcal{M}_{\tau}^{(m)}\left({{x}_{0}}\right)\partial _{{{x}_{0}}}^{m}.\end{eqnarray} \tag{ 189 }$$

This operator is the adjoint of the operator in the forward Kramers–Moyal expansion (185) in the sense that ${\int}^{}\text{d}x\,\left[\tilde{\mathcal{Q}}_{\tau}^{\dagger}\left(x,{{\partial}_{x}}\right)f(x)\right]g(x)={\int}^{}\text{d}xf(x)\left[{{\tilde{\mathcal{Q}}}_{\tau}}\left(x,{{\partial}_{x}}\right)g(x)\right]$ (provided that all boundary terms in the integrations by parts vanish).

4.4.2. Path integral representation of the backward Kramers–Moyal expansion.

The adjoint transition operator $\tilde{\mathcal{Q}}_{\tau}^{\dagger}\left({{x}_{0}},{{\partial}_{{{x}_{0}}}}\right)$ of the backward Kramers–Moyal expansion (189) is normal-ordered with respect to x₀ and ${{\partial}_{{{x}_{0}}}}$. Moreover, the backward Kramers–Moyal expansion (188) has the same form as the flow equation (116) obeyed by the marginalized distribution. One can therefore follow the steps in section 4.1 to represent the backward Kramers–Moyal expansion by the path integral

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} p\left(t,x|{{t}_{0}},{{x}_{0}}\right)={\rlap{-} \int}_{({{t}_{0}}}^{t]}\,{{\text{e}}^{-{{\mathcal{S}}^{\dagger}}}}\,\delta \left(x-x(t)\right){{|}_{x\left({{t}_{0}}\right)={{x}_{0}}}} \end{array}\end{eqnarray} \tag{ 190 }$$

with the action

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{S}}^{\dagger}}:={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\left[\text{i}q(\tau ){{\partial}_{\tau}}x(\tau )-\tilde{\mathcal{Q}}_{\tau}^{\dagger}\left(x(\tau ),\text{i}q(\tau )\right)\right]. \end{array}\end{eqnarray} \tag{ 191 }$$

The integral sign in (190) is again defined as the continuous-time limit of (167) and traces out all paths of $x(\tau )$ and $q(\tau )$ for $\tau \in \left({{t}_{0}},t\right]$ (note that $x(\tau )$ differs from x and x₀, which are fixed parameters; for brevity, however, $x(\tau )$ is occasionally abbreviated as x below). A diagrammatic computation of multi-time correlation functions based on a path integral representation equivalent to (190) has recently been considered in [385].

Let us specify the path integral representation (190) for a model of the chemical reaction $k\,A\to l\,A$ with rate coefficient ${{\gamma}_{\tau}}$ (and $k,l\in {{\mathbb{N}}_{0}}$). As we assume the particle 'number' x to be continuous, the model is only reasonable in an approximate sense for large values of x. In defining the transition rate of the reaction, one has to ensure the non-negativity of x and the conservation of probability. A possible choice of the transition rate is

$$\begin{eqnarray}&&{{\kappa}_{\tau}}(\Delta x,x):={{\gamma}_{\tau}}{{x}^{k}} \Theta (x-k)\delta \left(\Delta x-(l-k)\right).\end{eqnarray} \tag{ 192 }$$

Here, the Heaviside step function $\Theta (x-k)$ ensures that a sufficient number of particles are present to engage in a reaction (and it prevents the loss of probability to negative values of x).

Assuming a smooth approximation of the Heaviside step function that vanishes for x  <  0, the jump moments (186) follow as

$$\begin{eqnarray}&&\mathcal{M}_{\tau}^{(m)}(x)={{\gamma}_{\tau}}{{(l-k)}^{m}}{{x}^{k}} \Theta (x-k).\end{eqnarray} \tag{ 193 }$$

The corresponding adjoint transition operator evaluates to

$$\begin{eqnarray}&&\tilde{\mathcal{Q}}_{\tau}^{\dagger}\left(x,\text{i}q\right)={{\gamma}_{\tau}}{{\left(x{{\text{e}}^{-\text{i}q}}\right)}^{k}}\left[{{\left({{\text{e}}^{\text{i}q}}\right)}^{l}}-{{\left({{\text{e}}^{\text{i}q}}\right)}^{k}}\right] \Theta (x-k).\end{eqnarray} \tag{ 194 }$$

Path integrals with such an operator have been noted in [386–388], but their potential use remains to be fully explored. Curiously, upon ignoring the step function, the operator (194) has the same structure as the transition operator (122) upon identifying $x{{\text{e}}^{-\text{i}q}}$ with the creation operator c, and ${{\text{e}}^{\text{i}q}}$ with the annihilation operator a (Cole-Hopf transformation [386–390]). Note, however, that the stochastic processes associated to the two transition operators differ from each other (discrete versus continuous state space; different transition rates). The connection between the two associated path integrals (171) and (190) should be further explored.

4.4.3. Path integral representation of the backward Fokker–Planck equation.

With the drift coefficient ${{\alpha}_{\tau}}(x):=\mathcal{M}_{\tau}^{(1)}(x)$ and the non-negative diffusion coefficient ${{\beta}_{\tau}}(x):=\mathcal{M}_{\tau}^{(2)}(x)$, a truncation of the adjoint transition operator (189) at its second summand reads

$$\begin{eqnarray}&&\tilde{{\mathcal{Q}}}\,_{\tau}^{\dagger}\left({{x}_{0}},{{\partial}_{{{x}_{0}}}}\right)={{\alpha}_{\tau}}\left({{x}_{0}}\right){{\partial}_{{{x}_{0}}}}\,+\,\frac{1}{2}{{\beta}_{\tau}}\left({{x}_{0}}\right)\partial _{{{x}_{0}}}^{2}.\end{eqnarray} \tag{ 195 }$$

Thus, the corresponding Kramers–Moyal expansion (188) recovers the backward Fokker–Planck equation (2). Since the action (191) evaluates to

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{S}}^{\dagger}}:={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\left(\text{i}q\left[{{\partial}_{\tau}}x-{{\alpha}_{\tau}}(x)\right]+\frac{1}{2}{{\beta}_{\tau}}(x){{q}^{2}}\right), \end{array}\end{eqnarray} \tag{ 196 }$$

the corresponding path integral (190) coincides with a classic path integral representation of the (backward) Fokker–Planck equation. The original development of this representation goes back to works of Martin, Siggia, and Rose [16], de Dominicis [17], Janssen [18, 19], and Bausch, Janssen, and Wagner [19]. The application of this path integral to stochastic processes is, for example, discussed in the book of Täuber [70].

The transition operator (195) of the backward Fokker–Planck equation has the same form as the transition operator (176) in section 4.3, but it does not involve a perturbation operator $\mathcal{P}_{\tau}^{\dagger}$. Therefore, one can follow the steps in that section to evaluate the path integral (190) in terms of an average over the paths of an Itô SDE. This procedure shows that the backward Fokker–Planck equation is solved by

$$\begin{eqnarray}&&p\left(t,x|{{t}_{0}},{{x}_{0}}\right)={{\langle \langle \delta \left(x-x(t)\right)\rangle \rangle}_{W}},\end{eqnarray} \tag{ 197 }$$

with $x(\tau )$ solving the Itô SDE

$$\begin{eqnarray}&&\text{d}x(\tau )={{\alpha}_{\tau}}\left(x(\tau )\right)\text{d}\tau +\sqrt{{{\beta}_{\tau}}\left(x(\tau )\right)}\text{d}W(\tau )\end{eqnarray} \tag{ 198 }$$

with initial value $x\left({{t}_{0}}\right)={{x}_{0}}$. Hence, we have recovered the Feynman–Kac formula (3) (apart from a notational change in the time parameters).

4.4.4. The Onsager–Machlup function.

The connection between the above path integral representation of the (backward) Fokker–Planck equation and the work of Onsager and Machlup [79] becomes apparent upon the completion of a square (as in a Hubbard-Stratonovich transformation [391, 392]). For this purpose, the diffusion coefficient ${{\beta}_{\tau}}$ in the transition operator (195) must not only be non-negative but positive(-definite). One can then complete the square in the variable q and perform the path integration over this variable. Returning to the discrete-time approximation (169) of the action, one obtains the following representation of the conditional probability distribution $p\left(t,x|{{t}_{0}},{{x}_{0}}\right)$ in terms of convolutions of Gaussian distributions:

$$\begin{eqnarray}&&\underset{N\to \infty}{{\lim}}\,\left(\underset{j=0}{\overset{N-1}{\prod}}\,{\int}_{\mathbb{R}}\text{d}{{x}_{j+1}}\,{{\mathcal{G}}_{{{\mu}_{j}},\sigma _{j}^{2}}}\left({{x}_{j+1}}-{{x}_{j}}\right)\right)\delta \left(x-{{x}_{N}}\right).\end{eqnarray} \tag{ 199 }$$

The Dirac delta function is included in the integrations. Moreover, ${{\mu}_{j}}:={{\alpha}_{{{t}_{j}}}}\left({{x}_{j}}\right) \Delta t$ acts as the mean and $\sigma _{j}^{2}:={{\beta}_{{{t}_{j}}}}\left({{x}_{j}}\right) \Delta t$ as the variance of the Gaussian distribution

$$\begin{eqnarray}&&{{\mathcal{G}}_{\mu,{{\sigma}^{2}}}}(x):=\frac{{{\text{e}}^{-{{(x-\mu )}^{2}}/\left(2{{\sigma}^{2}}\right)}}}{\sqrt{2\pi {{\sigma}^{2}}}}.\end{eqnarray} \tag{ 200 }$$

As before, $\Delta t=\left(t-{{t}_{0}}\right)/N$ denotes the time intervals between ${{t}_{0}}\leqslant {{t}_{1}}\leqslant \cdots \leqslant {{t}_{N}}=t$.

Upon identifying $x\left({{t}_{0}}+j \Delta t\right)$ with x_j , the representation (199) can be rewritten as

$$\begin{eqnarray}&&p\left(t,x|{{t}_{0}},{{x}_{0}}\right)={\rlap{-} \int}_{({{t}_{0}}}^{t]}{{\text{e}}^{-{{\mathcal{S}}^{\dagger}}}}\,\delta \left(x-x(t)\right){{|}_{x\left({{t}_{0}}\right)={{x}_{0}}}}.\end{eqnarray} \tag{ 201 }$$

This integral proceeds only over paths $x(\tau )$ with $\tau \in \left({{t}_{0}},t\right]$ because the q-variables have already been integrated over. Paths are weighed by the exponential factor ${{\text{e}}^{-{{\mathcal{S}}^{\dagger}}}}$ with the action

$$\begin{eqnarray}&&{{\mathcal{S}}^{\dagger}}:=\underset{N\to \infty}{{\lim}}\,\underset{j=0}{\overset{N-1}{\sum}}\,\frac{{{\left[{{x}_{j+1}}-{{x}_{j}}-{{\alpha}_{{{t}_{j}}}}\left({{x}_{j}}\right) \Delta t\right]}^{2}}}{2{{\beta}_{{{t}_{j}}}}\left({{x}_{j}}\right) \Delta t}.\end{eqnarray} \tag{ 202 }$$

One may abbreviate the continuous-time limit of this action by the integral

$$\begin{eqnarray}&&{{\mathcal{S}}^{\dagger}}={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\frac{{{\left[{{\partial}_{\tau}}x-{{\alpha}_{\tau}}(x)\right]}^{2}}}{2{{\beta}_{\tau}}(x)}.\end{eqnarray} \tag{ 203 }$$

The integrand of this action is called an Onsager–Machlup function [79, 393] and the representation (201) may, consequently, be called an Onsager–Machlup representation (or 'functional' in the sense of functional integration). Note that the Onsager–Machlup representation involves the limit

$$\begin{eqnarray}&&{\rlap{-} \int}_{({{t}_{0}}}^{t]}:=\underset{N\to \infty}{{\lim}}\,\underset{j=0}{\overset{N-1}{\prod}}\,{{{\int}^{}}_{\mathbb{R}}}\frac{\text{d}{{x}_{j+1}}}{\sqrt{2\pi {{\beta}_{{{t}_{j}}}}\left({{x}_{j}}\right) \Delta t}}.\end{eqnarray} \tag{ 204 }$$

Mathematically rigorous formulations of the Onsager–Machlup representation have been attempted in [394–396]. For the other path integrals discussed in this review, such attempts have not yet been made to our knowledge.

4.4.5. Alternative discretization schemes.

Up to this point, we have discretized time in such a way that the evaluation of the path integrals eventually proceeds via the solution of an Itô stochastic differential equation (see section 4.3). Alternative discretization schemes have been proposed as well and are, for example, employed in stochastic thermodynamics [237]. To illustrate these schemes, let us, for simplicity, assume that the drift coefficient $\alpha (x)$ does not depend on time and that the diffusion coefficient $D:={{\beta}_{\tau}}$ is constant. A general discretization scheme—called the α-scheme but here we use a κ—consists of shifting the argument x_j of the drift coefficient in the action (202) to ${{\bar{x}}_{j}}:=\kappa \,{{x}_{j+1}}+(1-\kappa ){{x}_{j}}$ by writing

$$\begin{eqnarray}&&{{x}_{j}}={{\bar{x}}_{j}}-\kappa \left({{x}_{j+1}}-{{x}_{j}}\right)\end{eqnarray} \tag{ 205 }$$

(with $\kappa \in [0,1]$). Afterwards, the drift coefficient is expanded in powers of ${{x}_{j+1}}-{{x}_{j}}$, which is of order $\sqrt{ \Delta t}$ along relevant paths. Following [397], it suffices to keep only the first two terms of the expansion of $\alpha \left({{x}_{j}}\right)$ so that

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{x}_{j+1}}-{{x}_{j}}-\alpha \left({{x}_{j}}\right) \Delta t \\ {} & \approx \left[1+\kappa \,{{\alpha}^{\prime}}\left({{{\bar{x}}}_{j}}\right) \Delta t\right]\left({{x}_{j+1}}-{{x}_{j}}\right)-\alpha \left({{{\bar{x}}}_{j}}\right) \Delta t. \end{array}\end{eqnarray} \tag{ 206 }$$

As the next step, the integration variables x₁,..., x_N are transformed according to

$$\begin{eqnarray}&&\left[1+\kappa \,{{\alpha}^{\prime}}\left({{\bar{x}}_{j}}\right) \Delta t\right]\left({{x}_{j+1}}-{{x}_{j}}\right)\mapsto {{x}_{j+1}}-{{x}_{j}}.\end{eqnarray} \tag{ 207 }$$

The Jacobian of this transformation vanishes everywhere except on its diagonal and sub-diagonal. Its determinant therefore contributes an additional factor to the path integral and requires a redefinition of the action (202) as

$$\begin{eqnarray}&&{{\mathcal{S}}^{\dagger}}:=\underset{N\to \infty}{{\lim}}\,\underset{j=0}{\overset{N-1}{\sum}}\,\left(\frac{{{\left[{{x}_{j+1}}-{{x}_{j}}-\alpha \left({{{\bar{x}}}_{j}}\right) \Delta t\right]}^{2}}}{2D \Delta t}+\kappa \,{{\alpha}^{\prime}}\left({{\bar{x}}_{j}}\right) \Delta t\right).\end{eqnarray}$$

One may abbreviate this limit by

$$\begin{eqnarray}&&{{\mathcal{S}}^{\dagger}}={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\left(\frac{{{\left[{{\partial}_{\tau}}x-\alpha (x)\right]}^{2}}}{2D}+\kappa \,{{\alpha}^{\prime}}(x)\right).\end{eqnarray} \tag{ 208 }$$

The Itô version of this action with $\kappa =0$ has proved to be convenient in perturbation expansions of path integrals because so called 'closed response loops' can be omitted right from the start (see section 4.5 in [70]). Moreover, this prescription does not require a change of variables and makes the connection to the Feynman–Kac formula in appendix A most apparent. The Stratonovich version of the action with $\kappa =\frac{1}{2}$ is, for example, employed in Seifert's review on stochastic thermodynamics [237]. For a recent, more thorough discussion of the above discretization schemes, as well as of conflicting approaches, we refer the reader to [398] (the above action is discussed in the appendices).

4.4.6. Wiener's path integral.

After this detour to processes with continuous state spaces, let us return to the evaluation of the backward path integral representation (171). In the following, we show how the method introduced in section 4.3 can be applied to several elementary jump processes. We already applied a simplified version of the method in section 4.2 to solve the bi-directional reaction $\varnothing \rightleftharpoons A$. We now consider a generic reaction $k\,A\to l\,A$ with rate coefficient ${{\gamma}_{\tau}}$. Using the Poisson basis function $|n{{\rangle}_{x}}=\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}$, the adjoint transition operator of this reaction can be written as (see section 3.2.2)

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \tilde{\mathcal{Q}}_{\tau}^{\dagger}\left(x,{{\partial}_{x}}\right)={{\mathcal{Q}}^{\dagger}}(c,a)={{\gamma}_{\tau}}{{c}^{k}}\left({{a}^{l}}-{{a}^{k}}\right) \end{array}\end{eqnarray} \tag{ 209 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & =\left[{{\gamma}_{\tau}}(l-k){{x}^{k}}\right]{{\partial}_{x}}+\left[{{\gamma}_{\tau}}\left(l(l-1)-k(k-1)\right){{x}^{k}}\right]\frac{\partial _{x}^{2}}{2} \\ {} & \,\,\,\,\,\,+\,\,\mathcal{P}_{\tau}^{\dagger}\left(x,{{\partial}_{x}}\right). \end{array}\end{eqnarray} \tag{ 210 }$$

Here we employed the creation operator c  =  x and the annihilation operator $a={{\partial}_{x}}+1$, and we performed a series expansion with respect to ${{\partial}_{x}}$. Terms of third and higher order in ${{\partial}_{x}}$ were shoved into the perturbation operator $\mathcal{P}_{\tau}^{\dagger}\left(x,{{\partial}_{x}}\right)$. In the following, we approximate the marginalized distribution of the reaction $k\,A\to l\,A$ by first dropping both the diffusion coefficient and the perturbation operator from (210). Afterwards, we reintroduce the diffusion coefficient and show how the marginalized distribution follows as the average of a Poisson distribution over the paths of an Itô SDE. For the pair generation process $\varnothing \to 2\,A$, this representation is exact because the perturbation operator associated to this process is zero (see section 4.5.1). Later, in section 4.5.4, we solve a process with a non-vanishing perturbation operator $\mathcal{P}_{\tau}^{\dagger}$.

As a first approximation of the reaction $k\,A\to l\,A$, let us drop all the terms of the adjoint transition operator (210) except for the drift coefficient ${{\alpha}_{\tau}}(x)={{\gamma}_{\tau}}(l-k){{x}^{k}}$. The SDE (181) then simplifies to the deterministic rate equation ${{\partial}_{\tau}}x={{\gamma}_{\tau}}(l-k){{x}^{k}}$ of the reaction. Its solution $x(\tau )$ acts as the mean of the marginalized distribution, which, according to (177)–(180), is again given by the Poisson distribution $|\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{x\left({{t}_{0}}\right)}}=\,|n{{\rangle}_{t,x(t)}}$ in (174).

Going one step further, one may keep the diffusion coefficient in (210). The marginalized distribution then reads

$$\begin{eqnarray}&&|\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{x\left({{t}_{0}}\right)}}={{\langle \langle \frac{x{{(t)}^{n}}{{\text{e}}^{-x(t)}}}{n!}\rangle \rangle}_{W}}\end{eqnarray} \tag{ 211 }$$

with $x(\tau )$ solving the Itô SDE

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{d}x= & {{\gamma}_{\tau}}(l-k){{x}^{k}}\text{d}\tau \\ {} & +\sqrt{{{\gamma}_{\tau}}\left[l(l-1)-k(k-1)\right]{{x}^{k}}}\,\text{d}W(\tau ). \end{array}\end{eqnarray} \tag{ 212 }$$

Hence, the Poisson distribution in (211) is averaged over all possible sample paths of the SDE (212), whose evolution starts out from the initial value x(t₀) (see (136)). The above expression has also recently been derived by Wiese [336], although based on the forward path integral that we will discuss in section 5. Apparently, the expression under the square root in (212) is non-negative only if $k\in \left\{0,1\right\}$ or if $l\geqslant k\geqslant 2$. If that is not the case, $x(\tau )$ strays off into the complex domain. SDEs with imaginary noise (or the corresponding Langevin equations) have been studied in several recent articles, most notably for the binary annihilation reaction $2\,A\to \varnothing$ and for the coagulation process $2\,A\to A$ [336, 400, 401]. The numerical evaluation of such SDEs, however, often encountered severe convergence problems [336, 379]. The appearance of imaginary noise has been linked to the anti-correlation of particles in spatial systems [400], but as (212) shows, no spatial degrees of freedoms are actually required for its emergence. As Wiese pointed out, imaginary noise generally appears when, over time, the marginalized distribution (211) becomes narrower than a Poisson distribution [336].

4.5.1. Pair generation.

For the pair generation process $\varnothing \to 2\,A$ with growth rate coefficient ${{\gamma}_{\tau}}$, the drift and the (squared) diffusion coefficients agree: ${{\alpha}_{\tau}}={{\beta}_{\tau}}=2{{\gamma}_{\tau}}$. As neither of them depends on $x(\tau )$, the Itô SDE is readily solved by $x(t)=x\left({{t}_{0}}\right)+{\int}_{{{t}_{0}}}^{t}\text{d}\tau 2{{\gamma}_{\tau}}+{\int}_{{{t}_{0}}}^{t}\text{d}W(\tau )\sqrt{2{{\gamma}_{\tau}}}$. After introducing the (rather daunting) parameter

$$\begin{eqnarray}&&{{\eta}_{k}}:=\langle \langle {{\text{e}}^{-\left({\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\gamma}_{\tau}}+{\int}_{{{t}_{0}}}^{t}\text{d}W(\tau )\sqrt{2{{\gamma}_{\tau}}}\right)}}\\ &&\quad \centerdot \,\,{{\left({{\left({\int}_{{{t}_{0}}}^{t}\text{d}\tau \,2{{\gamma}_{\tau}}\right)}^{1/2}}+\frac{{\int}_{{{t}_{0}}}^{t}\text{d}W(\tau )\sqrt{2{{\gamma}_{\tau}}}}{{{\left({\int}_{{{t}_{0}}}^{t}\text{d}\tau \,2{{\gamma}_{\tau}}\right)}^{1/2}}}\right)}^{k}}\rangle {{\rangle}_{W}},\end{eqnarray} \tag{ 213 }$$

one can employ the binomial theorem to rewrite the right-hand side of the marginalized distribution (211) as

$$\begin{eqnarray}&&\frac{{{\text{e}}^{-\left(x\left({{t}_{0}}\right)+{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\gamma}_{\tau}}\right)}}}{n!}\underset{k=0}{\overset{n}{\sum}}\,\left(\begin{array}{c} n \\ k \end{array}\right){{\left({\int}_{{{t}_{0}}}^{t}\text{d}\tau \,2{{\gamma}_{\tau}}\right)}^{k/2}}{{\eta}_{k}}\,x{{\left({{t}_{0}}\right)}^{n-k}}.\end{eqnarray}$$

If the rate coefficient γ is independent of time, one can show that the parameter ${{\eta}_{k}}$ coincides with the kth moment of a Gaussian distribution with zero mean and unit variance (i.e. ${{\eta}_{k}}=0$ for odd k and ${{\eta}_{k}}=(k-1)!!$ for even k; note that the sum of two Gaussian random variables is again a Gaussian random variable, with its mean and variance following additively). The marginalized distribution $|\,p\left(t,n|{{t}_{0}}\right){{\rangle}_{x\left({{t}_{0}}\right)}}$ therefore follows as

$$\begin{eqnarray}&&{{\text{e}}^{-\left(x\left({{t}_{0}}\right)+{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\gamma}_{\tau}}\right)}}\underset{k=0}{\overset{\lfloor n/2\rfloor}{\sum}}\,\frac{{{\left({\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\gamma}_{\tau}}\right)}^{k}}}{k!}\,\frac{{{\left(x\left({{t}_{0}}\right)\right)}^{n-2k}}}{(n-2k)!}.\end{eqnarray} \tag{ 214 }$$

Here, $\lfloor n/2\rfloor$ represents the integral part of n/2. The marginalized distribution (214) solves the master equation of the pair generation process and is initially of Poisson shape. For x(t₀)  =  0, the distribution effectively keeps that shape, although only on the set of all even numbers (the reaction $\varnothing \to 2\,A$ cannot create an odd number of particles when starting out from zero particles). Using Mathematica by Wolfram Research, we verified that the distribution also applies when ${{\gamma}_{\tau}}$ depends on time. The evolution of the distribution is shown in figure 7 for the rate coefficient ${{\gamma}_{\tau}}=1/\left(1+\sqrt{\tau}\right)$. We refrain from computing the conditional distribution from the marginalized distribution because the computation is unwieldy and does not shed any more light on the path integral approach.

Zoom In Zoom Out Reset image size

4.5.2. Diffusion on networks.

The backward path integral can also be used to solve spatial processes. In fact, stochastic path integrals have been most useful in the study of such processes. If particles engaging in a chemical reaction can also diffuse in space, their density may evolve in ways that are not expected from the well-mixed, non-spatial limit. That is, for example, the case for particles that annihilate one another in the reaction $2\,A\to \varnothing$ while diffusing along a one-dimensional line. From the well-mixed limit, one would expect that the particle density decays asymptotically as t⁻¹ with time (see section 7), but instead it decays as t^−1/2 [130, 294, 402]. The reason behind this surprising decay law is the rapid condensation of the system into a state in which isolated particles are separated by large voids. From that point on, further decay of the particle density requires that two particles first find each other through diffusion. Therefore, the process is 'diffusion-limited' and exhibits a slower decay of the particles (see e.g. [403] on the related 'first-passage' problem). Path integral representations of the master equation, combined with renormalization group techniques, have been successfully applied to the computation of decay laws in spatially extended systems, both regarding the (universal) exponents and the pre-factors of these laws [130–132, 294, 306, 404]. For a broader discussion of systems exhibiting a transition into an absorbing state see [68]. Before turning to the particle density and, more generally, to observables of the particle number, we now show how the backward path integral helps in computing the full probability distribution of a spatial process.

As a first step, we consider a pure diffusion process with particles hopping between neighbouring nodes of a network $\mathbb{L}$. For now, the topology of the network may be arbitrary, being either random or regular. In the next section, the network topology will be chosen as a regular lattice. In the limit of a vanishing lattice spacing (and an infinite number of lattice sites), a 'field theory' will be obtained. In section 4.5.4, the particles will also be allowed to decay.

The configuration of particles on the network may be represented by the vector $\boldsymbol{n}\in \mathbb{N}_{0}^{|\mathbb{L}|}$, with $|\mathbb{L}|$ being the total number of network nodes. The configuration changes whenever a particle hops from some node $i\in \mathbb{L}$ to a neighbouring node $j\in {{{\mathsf {N}}}_{i}}\subset \mathbb{L}$. The probability $p\left(\tau,\boldsymbol{n}|{{t}_{0}},\boldsymbol{n}_{{0}}\right)$ of finding the system in configuration $\boldsymbol{n}$ then obeys the master equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}p\left(\tau,\boldsymbol{n}|\centerdot \right)=\underset{i\in \mathbb{L}}{\sum}\,{{\varepsilon}_{\tau,i}}\underset{j\in {{{\mathsf {N}}}_{i}}}{\sum}\,\left[\left({{n}_{i}}+1\right)p\left(\tau,\boldsymbol{n}+{{\widehat{\boldsymbol{e}}}_{i}}-{{\widehat{\boldsymbol{e}}}_{j}}|\centerdot \right)\right. \\ &&\left. \qquad \qquad \ \;-{{n}_{i}}\,p\left(\tau,\boldsymbol{n}|\centerdot \right)\right].\end{eqnarray} \tag{ 215 }$$

Here, ${{\widehat{\boldsymbol{e}}}_{i}}$ represents a unit vector that points in direction $i\in \mathbb{L}$. Moreover, ${{\varepsilon}_{\tau,i}}$ acts as a hopping rate and may depend both on time and on the node from which a particle departs. Apparently, the above master equation has the same structure as the chemical master equation (24). It may therefore be cast into a linear PDE by extending the operators and basis functions from section 3.2.2 to multiple variables. In particular, we extend the Poisson basis function to

$$\begin{eqnarray}&&|\boldsymbol{n}{{\rangle}_{\boldsymbol{x}}}:=\underset{i\in \mathbb{L}}{\prod}\,\frac{x_{i}^{{{n}_{i}}}{{\text{e}}^{-{{x}_{i}}}}}{{{n}_{i}}!}\end{eqnarray} \tag{ 216 }$$

and employ the creation and annihilation operators

$$\begin{eqnarray}&&\boldsymbol{c}:=\boldsymbol{x}~\text{and}~\boldsymbol{a}:=\nabla +\mathbf{1}.\end{eqnarray} \tag{ 217 }$$

According to the flow equation (116), the marginalized distribution evolves via ${{\partial}_{-\tau}}|\,p(t,n|\tau )\rangle =\tilde{\mathcal{Q}}_{\tau}^{\dagger}|\,p(t,n|\tau )\rangle$ with the adjoint transition operator

$$\begin{eqnarray}&&\tilde{\mathcal{Q}}_{\tau}^{\dagger}\left(\boldsymbol{x},\nabla \right)=\underset{i\in \mathbb{L}}{\sum}\,\left({\underline \Delta }{{\varepsilon}_{\tau,i}}{{x}_{i}}\right){{\partial}_{{{x}_{i}}}}.\end{eqnarray} \tag{ 218 }$$

Here we introduced the discrete Laplace operator ${\underline \Delta }{{f}_{i}}:={\sum}_{j\in {{{\mathsf {N}}}_{i}}}\left(\,{{f}_{j}}-{{f}_{i}}\right)$, which acts both on ${{\varepsilon}_{\tau,i}}$ and on x_i .

By following the steps in section 4.3, one finds that the marginalized distribution is given by the multivariate Poisson distribution

$$\begin{eqnarray}&&|\,p\left(t,\boldsymbol{n}|{{t}_{0}}\right){{\rangle}_{\boldsymbol{x}\left({{t}_{0}}\right)}}=\underset{i\in \mathbb{L}}{\prod}\,\frac{{{x}_{i}}{{(t)}^{{{n}_{i}}}}{{\text{e}}^{-{{x}_{i}}(t)}}}{{{n}_{i}}!}.\end{eqnarray} \tag{ 219 }$$

Its mean $\boldsymbol{x}(t)$ solves the discrete diffusion equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}{{x}_{i}}=\underset{j\in {{{\mathsf {N}}}_{i}}}{\sum}\,\left({{\varepsilon}_{\tau,j}}{{x}_{j}}-{{\varepsilon}_{\tau,i}}{{x}_{i}}\right)={\underline \Delta }{{\varepsilon}_{\tau,i}}{{x}_{i}},\end{eqnarray} \tag{ 220 }$$

with the initial condition $\boldsymbol{x}\left({{t}_{0}}\right)$. Let us note that the marginalized distribution solves the master equation (215), but with the initial number of particles being Poisson distributed locally.

The conditional distribution $p\left(t,\boldsymbol{n}|{{t}_{0}},\boldsymbol{n}_{{0}}\right)$ follows from the marginalized distribution (219) by applying the functional

$$\begin{eqnarray}&&\langle \boldsymbol{n}_{{\mathbf{0}}}{{|}_{\boldsymbol{x}\left({{t}_{0}}\right)}}\,f=\left[\underset{i\in \mathbb{L}}{\prod}\,{{\left({{\partial}_{{{x}_{i}}\left({{t}_{0}}\right)}}+1\right)}^{{{n}_{0,i}}}}\right]\,f\left(\boldsymbol{x}\left({{t}_{0}}\right)\right){{|}_{\boldsymbol{x}\left({{t}_{0}}\right)=\mathbf{0}}}\end{eqnarray} \tag{ 221 }$$

to it (see (126) with (129)). The evaluation requires a prior solution of the discrete diffusion equation (220). This equation can be written in matrix form as ${{\partial}_{\tau}}\boldsymbol{x}={{M}_{\tau}}\boldsymbol{x}$ with ${{M}_{\tau,ij}}=\left({{A}_{ij}}\,-\,|{{{\mathsf {N}}}_{i}}|{{\delta}_{ij}}\right){{\varepsilon}_{\tau,j}}$. Here, A represents the symmetric adjacency matrix of the network and $|{{{\mathsf {N}}}_{i}}|$ represents the number of neighbours of node $i\in \mathbb{L}$. The matrix equation can in principle be solved through a Magnus expansion. The solution has the generic form $\boldsymbol{x}(\tau )=G\left(\tau |{{t}_{0}}\right)\boldsymbol{x}\left({{t}_{0}}\right)$ with the 'propagator' G solving ${{\partial}_{\tau}}G\left(\tau |{{t}_{0}}\right)={{M}_{\tau}}G\left(\tau |{{t}_{0}}\right)$. We write the elements of the propagator as $G\left(\tau,i|{{t}_{0}},\,j\right)$. Its flow starts out from $G\left({{t}_{0}}|{{t}_{0}}\right)={\mathbb {1}}$. For later purposes, let us note the time-reversal property in terms of the matrix inverse $G{{(t|\tau )}^{-1}}=G(\tau |t)$ and also the conservation law ${{\mathbf{1}}^{\intercal}}G\left(\tau |{{t}_{0}}\right)={{\mathbf{1}}^{\intercal}}$.

It proves insightful to consider the master equation (215) of the multi-particle hopping process for the random walk of a single particle on the one-dimensional lattice $\mathbb{L}=\mathbb{Z}$, which we already considered in sections 2.2.1 and 3.2.1 (with symmetric hopping rates ${{\varepsilon}_{\tau,i}}={{l}_{\tau}}={{r}_{\tau}}$). The presence of only a single particle can be enforced by choosing $\boldsymbol{n}_{{0}}$ as n_0,k   =  1 for one $k\in \mathbb{Z}$ and n_0,j   =  0 for all $j\ne k$. By following the above steps, one eventually finds that the master equation (215) is solved by the conditional probability distribution

$$\begin{eqnarray}&&p\left(t,\boldsymbol{n}|{{t}_{0}},\boldsymbol{n}_{{0}}\right)=\underset{i\in \mathbb{Z}}{\sum}\,\left(\underset{j\ne i}{\prod}\,{{\delta}_{{{n}_{j}},0}}\right){{\delta}_{{{n}_{i}},1}}G\left(t,i|{{t}_{0}},k\right),\end{eqnarray} \tag{ 222 }$$

with the propagator G solving the master equation (55) of the simple random walk. Hence, the propagator evaluates to a Skellam distribution.

4.5.3. Diffusion in continuous space.

To make the transition to a field theory, one may specify the network as the the d-dimensional lattice $\mathbb{L}={{\left(l\mathbb{Z}\right)}^{d}}$ with the lattice spacing l  >  0 going to zero. In order to take this limit, we define the variable $x\left(\tau,\boldsymbol{r}\right):={{x}_{\boldsymbol{r}}}(\tau )$, the rescaled Laplace operator $\Delta :={\underline \Delta }/{{l}^{2}}$, and the rescaled diffusion coefficient ${{D}_{\tau}}\left(\boldsymbol{r}\right):={{l}^{2}}{{\varepsilon}_{\tau,\boldsymbol{r}}}$ for $\boldsymbol{r}\in \mathbb{L}$. Assuming that these definitions can be continued to $\boldsymbol{r}\in {{\mathbb{R}}^{d}}$, the discrete diffusion equation (220) becomes a PDE for the 'field' $x\left(\tau,\boldsymbol{r}\right)$ in the limit $l\to 0$, namely ${{\partial}_{\tau}}x\left(\tau,\boldsymbol{r}\right)= \Delta \left({{D}_{\tau}}\left(\boldsymbol{r}\right)x\left(\tau,\boldsymbol{r}\right)\right)$. Here, $\Delta$ represents the ordinary Laplace operator ¹⁰ . The solution of the PDE acts as the mean of the multivariate Poisson distribution (219) whose extension to $\boldsymbol{r}\in \mathbb{R}$ is, however, not quite obvious. If the diffusion coefficient is homogeneous in space, the PDE is solved by

$$\begin{eqnarray}&&x\left(\tau,\boldsymbol{r}\right)={{{\int}^{}}_{{{\mathbb{R}}^{d}}}}\text{d}\boldsymbol{r}_{{0}}\,G\left(\tau,\boldsymbol{r}|{{t}_{0}},\boldsymbol{r}_{{0}}\right)x\left({{t}_{0}},\boldsymbol{r}_{{0}}\right),\end{eqnarray} \tag{ 224 }$$

with the Gaussian kernel

$$\begin{eqnarray}&&G\left(\tau,\boldsymbol{r}|{{\tau}^{\prime}},\boldsymbol{r}^{{\prime}}\right)=\frac{{{\text{e}}^{-{{\left(\boldsymbol{r}-\boldsymbol{r}^{{\prime}}\right)}^{2}}/4{\int}_{{{\tau}^{\prime}}}^{\tau}\text{d}s\,{{D}_{s}}}}}{\sqrt{4\pi {\int}_{{{\tau}^{\prime}}}^{\tau}\text{d}s\,{{D}_{s}}}}.\end{eqnarray} \tag{ 225 }$$

The action (172) can also be extended into continuous space. For that purpose, one may rescale the second path integral variable as $q\left(\tau,\boldsymbol{r}\right):={{l}^{-d}}{{q}_{\boldsymbol{r}}}(\tau )$ so that in the limit $l\to 0$,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{S}}^{\dagger}} & ={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\underset{\boldsymbol{r}\in \mathbb{L}}{\sum}\,\left(\text{i}{{q}_{\boldsymbol{r}}}\left({{\partial}_{\tau}}{{x}_{\boldsymbol{r}}}-{\underline \Delta }{{\varepsilon}_{\tau,\boldsymbol{r}}}{{x}_{\boldsymbol{r}}}\right)\right) \end{array}\end{eqnarray} \tag{ 226 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \to {\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{\int}_{{{\mathbb{R}}^{d}}}\text{d}\boldsymbol{r}\,\text{i}q\left({{\partial}_{\tau}}x- \Delta {{D}_{\tau}}x\right). \end{array}\end{eqnarray} \tag{ 227 }$$

The above rescaling of x and q is not unique and depends on the problem at hand. Instead of dividing ${{q}_{\boldsymbol{r}}}$ by the volume factor l^d , this factor is sometimes employed to cast ${{x}_{\boldsymbol{r}}}$ and the particle number ${{n}_{\boldsymbol{r}}}$ into densities [130]. In the study of branching and annihilating random walks with an odd number of offspring, yet another kind of rescaling may bring the action into a Reggeon field theory [70, 405–407] like form [96]. The critical behaviour of these random walks falls into the universality class of directed percolation (DP) [96, 132]. Information on this universality class can be found in the book [68], as well as in the original articles of Janssen and Grassberger, which established the extensive scope of the DP class [408, 409].

4.5.4. Diffusion and decay.

Let us exemplify how one can accommodate a non-vanishing perturbation operator $\mathcal{P}_{\tau}^{\dagger}$ in the adjoint transition operator (176). As in the section before the last, we consider a system of particles that are hopping between the nodes of an arbitrary network $\mathbb{L}$. The corresponding diffusive transition operator (218) shall specify the drift coefficient ${{\alpha}_{\tau,i}}\left(\boldsymbol{x}\right)={\underline \Delta }{{\varepsilon}_{\tau,i}}{{x}_{i}}$ in the multivariate extension of (176). The coefficient ${{\beta}_{\tau,ij}}\left(\boldsymbol{x}\right)$ is zero. Besides allowing for diffusion, we now allow the particles to decay in the linear reaction $A\to \varnothing$. For the sake of brevity, we assume that the corresponding decay rate coefficient ${{\mu}_{\tau}}$ is spatially homogeneous. This assumption can be relaxed but the equations that result cannot be written in matrix form and involve many indices. We treat the adjoint transition operator of the decay process as the perturbation, which reads, according to the flow equation (131),

$$\begin{eqnarray}&&\mathcal{P}_{\tau}^{\dagger}\left(\boldsymbol{x},\nabla \right)=-{{\mu}_{\tau}}\underset{i\in \mathbb{L}}{\sum}\,{{x}_{i}}{{\partial}_{{{x}_{i}}}}.\end{eqnarray} \tag{ 228 }$$

The combined process could also be solved directly by treating the decay process as part of the drift coefficient ${{\alpha}_{\tau,i}}\left(\boldsymbol{x}\right)$. For pedagogic reasons, however, we wish to outline a perturbative solution using Feynman diagrams.

The first step of the derivation is to solve the differential equation (181) for $\boldsymbol{x}$, which now reads ${{\partial}_{\tau}}{{x}_{i}}={\underline \Delta }{{\varepsilon}_{\tau,i}}{{x}_{i}}+{{X}_{i}}(\tau )$. The homogeneous solution of this equation is given by $\boldsymbol{x}_{{h}}(\tau ):=G\left(\tau |{{t}_{0}}\right)\boldsymbol{x}\left({{t}_{0}}\right)$, with G being the propagator from the end of section 4.5.2. The solution of the full equation can be written as

$$\begin{eqnarray}&&\boldsymbol{x}(\tau )=\boldsymbol{x}_{{h}}(\tau )+{\int}_{{{t}_{0}}}^{\tau}\text{d}{{\tau}^{\prime}}G\left(\tau |{{\tau}^{\prime}}\right)\boldsymbol{X}\left({{\tau}^{\prime}}\right).\end{eqnarray} \tag{ 229 }$$

As the next step, the evaluation of the marginalized distribution (177) requires us to compute the $\left(\boldsymbol{Q},\boldsymbol{X}\right)$-generating functional (179) for $\boldsymbol{Q}=\boldsymbol{X}=\mathbf{0}$. By performing a series expansion of its leading exponential, this function can be written as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{Z}}_{\mathbf{0},\mathbf{0}}}= & \underset{k=0}{\overset{\infty}{\sum}}\,\underset{l=0}{\overset{k}{\sum}}\,\frac{1}{l!}{{\left(\text{i}\boldsymbol{q}_{{N}}\centerdot \frac{\delta}{\delta \boldsymbol{Q}(t)}\right)}^{l}} \\ {} & \centerdot {\int}_{{{t}_{0}}}^{t}\text{d}{{\tau}_{k-l}}\,\mathcal{P}_{{{\tau}_{k-l}}}^{\dagger}\cdots {\int}_{{{t}_{0}}}^{{{\tau}_{2}}}\text{d}{{\tau}_{1}}\,\mathcal{P}_{{{\tau}_{1}}}^{\dagger}\,\mathcal{Z}_{\boldsymbol{Q},\boldsymbol{X}}^{0}{{|}_{\boldsymbol{Q}=\boldsymbol{X}=\mathbf{0}}}, \end{array}\end{eqnarray} \tag{ 230 }$$

with $\mathcal{P}_{\tau}^{\dagger}=\frac{\delta}{\delta \boldsymbol{X}(\tau )}\centerdot \left(-{{\mu}_{\tau}}\right)\frac{\delta}{\delta \boldsymbol{Q}(\tau )}$. To evaluate ${{\mathcal{Z}}_{\mathbf{0},\mathbf{0}}}$, it helps to note that $\ln \mathcal{Z}_{\boldsymbol{Q},\boldsymbol{X}}^{0}={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\boldsymbol{Q}\centerdot \boldsymbol{x}$, from which it follows that for ${{t}_{0}}\leqslant \tau \leqslant t$:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\delta \ln \mathcal{Z}_{\boldsymbol{Q},\boldsymbol{X}}^{0}}{\delta \boldsymbol{Q}(\tau )}{{|}_{\boldsymbol{Q}=\boldsymbol{X}=\mathbf{0}}}=\boldsymbol{x}_{{h}}(\tau )\,\,\,\text{and} \end{array}\end{eqnarray} \tag{ 231 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{\delta}^{2}}\ln \mathcal{Z}_{\boldsymbol{Q},\boldsymbol{X}}^{0}}{\delta {{Q}_{i}}(\tau )\delta {{X}_{j}}\left({{\tau}^{\prime}}\right)}{{|}_{\boldsymbol{Q}=\boldsymbol{X}=\mathbf{0}}}=G\left(\tau,i|{{\tau}^{\prime}},j\right) \Theta \left(\tau -{{\tau}^{\prime}}\right). \end{array}\end{eqnarray} \tag{ 232 }$$

The Heaviside step function is used with $\Theta (0):=0$ to take into account that $\boldsymbol{x}_{{h}}(\tau )$ depends on $\boldsymbol{X}\left({{\tau}^{\prime}}\right)$ only for ${{\tau}^{\prime}}<\tau$ (see appendix D). All other derivatives of $\ln \mathcal{Z}_{\boldsymbol{Q},\boldsymbol{X}}^{0}$, as well as itself, vanish for $\boldsymbol{Q}=\boldsymbol{X}=\mathbf{0}$.

Upon inserting the perturbation (228) into (230), one observes that every summand of the expansion can be written in terms of a combination of the derivatives (231) and (232), and a terminal factor $\text{i}\boldsymbol{q}_{{N}}$. For example, the summand with k  =  2 and l  =  1 reads

$$\begin{eqnarray}&&\text{i}\boldsymbol{q}_{{N}}\centerdot \frac{\delta}{\delta \boldsymbol{Q}(t)}{\int}_{{{t}_{0}}}^{t}\text{d}{{\tau}_{1}}\,\frac{\delta}{\delta \boldsymbol{X}\left({{\tau}_{1}}\right)}\centerdot \left(-{{\mu}_{{{\tau}_{1}}}}\right)\frac{\delta}{\delta \boldsymbol{Q}\left({{\tau}_{1}}\right)}\,\mathcal{Z}_{\boldsymbol{Q},\boldsymbol{X}}^{0},\end{eqnarray} \tag{ 233 }$$

with $\boldsymbol{Q}=\boldsymbol{X}=\mathbf{0}$ being taken after the evaluation of the functional derivatives. The evaluation of these derivatives results in

$$\begin{eqnarray}&&\text{i}\boldsymbol{q}_{{N}}\centerdot {\int}_{{{t}_{0}}}^{t}\text{d}{{\tau}_{1}}\,G\left(t|{{\tau}_{1}}\right)\left(-{{\mu}_{{{\tau}_{1}}}}\right)\boldsymbol{x}_{{h}}\left({{\tau}_{1}}\right).\end{eqnarray} \tag{ 234 }$$

One can represent this expression graphically by a Feynman diagram according to the following rules. First, every diagram ends in a sink , which contributes the factor $\text{i}\boldsymbol{q}_{{N}}$. The incoming line, or 'leg', of the sink represents the derivative $\frac{\delta}{\delta \boldsymbol{Q}(t)}$. This leg may either be left dangling, resulting in a factor $\boldsymbol{x}_{{h}}(t)$ according to (231), or it may be 'contracted' with the outgoing line of a vertex . The two legs of this vertex reflect the two derivatives in the perturbation $\mathcal{P}_{\tau}^{\dagger}=\frac{\delta}{\delta \boldsymbol{X}(\tau )}\centerdot \left(-{{\mu}_{\tau}}\right)\frac{\delta}{\delta \boldsymbol{Q}(\tau )}$. According to (232), the contraction results in a propagator $G(t,i|\tau,j)$ and the vertex itself contributes a factor $-{{\mu}_{\tau}}$. The incoming leg of the vertex may again be left dangling or it may be connected to another vertex. Therefore, each Feynman diagram is a straight line for the linear decay process. The expression (234) can therefore be represented graphically by

Note that dangling outgoing lines are not permitted because the derivative $\delta \ln \mathcal{Z}_{\boldsymbol{Q},\boldsymbol{X}}^{0}/\delta \boldsymbol{X}(\tau )$ vanishes for $\boldsymbol{Q}=\boldsymbol{X}=\mathbf{0}$.

For more complex, non-linear processes, the Feynman diagrams may contain multiple kinds of vertices, each representing an individual summand of $\mathcal{P}_{\tau}^{\dagger}$. If there exist vertices with more than two legs, the diagrams may exhibit internal loops (see section 6.4). It is then usually impossible to evaluate the full perturbation series and it needs to be truncated at a certain order in the number of loops. Further information about these techniques, and about how renormalization group theory comes into play, is provided, for example, by the book of Täuber [70]. Information on a non-perturbative renormalization group technique can be found in [410].

For the simple diffusion process with decay, all the summands of (230) are readily cast into Feynman diagrams. However, it turns out that individual summands of the expansion may be associated to multiple diagrams that are not connected to one another. Furthermore, diagrams that represent summands of lower order keep reappearing as unconnected components of summands of higher order. Hence, there appears to be redundant information involved. This redundancy is removed by a classical theorem from diagrammatic analysis. This theorem states that the logarithm $\ln {{\mathcal{Z}}_{\mathbf{0},\mathbf{0}}}$ is given by the sum of only the connected diagrams [411], i.e. by

This sum can be evaluated with the help of $G(t|\tau )G\left(\tau |{{t}_{0}}\right)=G\left(t|{{t}_{0}}\right)$ as

$$\begin{eqnarray}\begin{array}{*{35}{l}} \ln {{\mathcal{Z}}_{\mathbf{0},\mathbf{0}}} & =\underset{k=0}{\overset{\infty}{\sum}}\,\text{i}\boldsymbol{q}_{{N}}\centerdot {\int}_{{{t}_{0}}}^{t}\text{d}{{\tau}_{k}}\,G\left(t|{{\tau}_{k}}\right)\left(-{{\mu}_{{{\tau}_{k}}}}\right) \\ \,\,\,\,\, & \,\,\,\,\,\cdots {\int}_{{{t}_{0}}}^{{{\tau}_{2}}}\text{d}{{\tau}_{1}}\,G\left({{\tau}_{2}}|{{\tau}_{1}}\right)\left(-{{\mu}_{{{\tau}_{1}}}}\right)\boldsymbol{x}_{{h}}\left({{\tau}_{1}}\right) \end{array}\end{eqnarray} \tag{ 237 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & =\text{i}\boldsymbol{q}_{{N}}\centerdot \boldsymbol{x}_{{h}}(t){{\text{e}}^{-{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\mu}_{\tau}}}}. \end{array}\end{eqnarray} \tag{ 238 }$$

After inserting this expression into the marginalized distribution (177), one recovers the multivariate Poisson distribution (219). Its mean, however, has now acquired the pre-factor ${{\text{e}}^{-{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\mu}_{\tau}}}}$, reflecting the decay of the particles.

In the present section, we introduced the novel backward path integral representation

$$\begin{eqnarray}&&|\,p\left(t,n|{{t}_{0}}\right)\rangle ={\rlap{-} \int}_{({{t}_{0}}}^{t]}\,{{\text{e}}^{-{{\mathcal{S}}^{\dagger}}}}\,|n{{\rangle}_{t}}\end{eqnarray} \tag{ 239 }$$

of the marginalized distribution (see section 3)

$$\begin{eqnarray}&&|\,p\left(t,n|{{t}_{0}}\right)\rangle =\underset{{{n}_{0}}}{\sum}\,p\left(t,n|{{t}_{0}},{{n}_{0}}\right)|{{n}_{0}}{{\rangle}_{{{t}_{0}}}}.\end{eqnarray} \tag{ 240 }$$

When the basis function $|n{{\rangle}_{x}}$ is chosen as a probability distribution (e.g. $|n{{\rangle}_{x}}=\frac{{{x}^{n}}{{\text{e}}^{-x}}}{n!}$ for all $n\in {{\mathbb{N}}_{0}}$), the backward path integral represents a true probability distribution: the marginalized distribution $|\,p\left(t,n|{{t}_{0}}\right)\rangle$. This distribution solves the forward master equation (10) for the initial condition $|\,p\left({{t}_{0}},n|{{t}_{0}}\right)\rangle =|n\rangle$ and it transforms into the conditional probability distribution as $p\left(t,n|\tau,{{n}_{0}}\right)=\langle {{n}_{0}}|\,p(t,n|\tau )\rangle$. In section 4.3, we showed how the backward path integral (239) can be expressed in terms of an average over the paths of an Itô stochastic differential equation. This method provided the exact solutions of various elementary stochastic processes, including the simple growth, the linear decay, and the pair generation processes. Moreover, we showed how the path integral can be evaluated perturbatively using Feynman diagrams for a process with diffusion and linear particle decay. We hope that the backward path integral (239) will prove useful in the study of reaction–diffusion processes. Thus far, the critical behaviour of such processes could only be approached via path integral representations of averaged observables or of the generating function [69, 412]. In section 6, we show how the former representation readily follows from the backward path integral (239) upon summing the marginalized distribution over an observables A(n). The corresponding representation is commonly applied in the study of diffusion-limited reactions (e.g. [69, 96, 130–132], but we here show how it can be freed of some of its quantum mechanical ballast (such as 'second-quantized' or 'normal-ordered' observables and coherent states). Besides, we showed in section 4.4 how one can derive a path integral representation of processes with continuous state spaces whose (backward) master equations admit a Kramers–Moyal expansion. Provided that this expansion stops at the level of a diffusion approximation, one recovers a classic path integral representation of the (backward) Fokker–Planck equation and also the Feynman–Kac formula (3). Moreover, the representation can be rewritten in terms of an Onsager–Machlup function and, for diffusive Brownian motion, it simplifies to the path integral of Wiener.

The forward path integral can be derived both from the flow equation (54) obeyed by the generating function or from the flow equation (150) obeyed by the series (149) (see figure 8). As it is more convenient to work with functions than with functionals, we use the flow equation of the generating function for this purpose. The derivation parallels a derivation of Elgart and Kamenev [35]. As in section 4.1, we first split the time interval [t₀, t] into N pieces ${{t}_{0}}\leqslant {{t}_{1}}\leqslant \cdots \leqslant {{t}_{N}}:=t$ of length $\Delta t$. Over a sufficiently small interval $\Delta t$, the flow equation (54) is then solved by

$$\begin{eqnarray}&&|g\left(t|{{t}_{0}},{{n}_{0}}\right){{\rangle}_{{{q}_{N}}}}={{\mathcal{L}}_{{{t}_{N-1}}}}\left({{q}_{N}},{{\partial}_{{{q}_{N}}}}\right)|g\left({{t}_{N-1}}|{{t}_{0}},{{n}_{0}}\right){{\rangle}_{{{q}_{N}}}},\end{eqnarray} \tag{ 241 }$$

with the generator ${{\mathcal{L}}_{\tau}}:=1+{{\tilde{\mathcal{Q}}}_{\tau}} \Delta t+\text{O}\left({{(\Delta t)}^{2}}\right)$. After inserting the integral form of a Dirac delta between $\mathcal{L}$ and $|g\rangle$, the right-hand side of the equation reads

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\mathcal{L}}_{{{t}_{N-1}}}}\left({{q}_{N}},{{\partial}_{{{q}_{N}}}}\right){{{\int}^{}}_{{{\mathbb{R}}^{2}}}}\frac{\text{d}{{q}_{N-1}}\text{d}{{x}_{N-1}}}{2\pi}{{\text{e}}^{-\text{i}{{x}_{N-1}}\left({{q}_{N-1}}-{{q}_{N}}\right)}}\,\,|g\left({{t}_{N-1}}|\centerdot \right){{\rangle}_{{{q}_{N-1}}}}. \end{array}\end{eqnarray} \tag{ 242 }$$

Zoom In Zoom Out Reset image size

Assuming that the transition operator ${{\tilde{\mathcal{Q}}}_{\tau}}$, and therefore also ${{\mathcal{L}}_{\tau}}$, are normal-ordered with respect to q and ${{\partial}_{q}}$, we may replace ${{\partial}_{{{q}_{N}}}}$ by $\text{i}{{x}_{N-1}}$ and interchange ${{\mathcal{L}}_{{{t}_{N-1}}}}$ with the exponential. This procedure is repeated N times before invoking the exponentiation ${{\mathcal{L}}_{\tau}}=\exp \left({{\tilde{\mathcal{Q}}}_{\tau}} \Delta t\right)$. Using the inverse transformation (48) and the initial condition $|g\left({{t}_{0}}|{{t}_{0}},{{n}_{0}}\right)\rangle =\,|{{n}_{0}}\rangle$, one obtains the discrete-time path integral representation

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & p\left(t,n|{{t}_{0}},{{n}_{0}}\right)={{\langle n{{|}_{t,{{q}_{N}}}}|g\left(t|{{t}_{0}},{{n}_{0}}\right)\rangle}_{{{q}_{N}}}} \end{array}\end{eqnarray} \tag{ 243 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \text{with}\,|g\left(t|{{t}_{0}},{{n}_{0}}\right){{\rangle}_{{{q}_{N}}}}={\rlap{-} \int}_{0}^{\,N-1}\,{{\text{e}}^{-{{\mathcal{S}}_{N}}}}\,|{{n}_{0}}{{\rangle}_{{{t}_{0}},{{q}_{0}}}}\,\text{and} \end{array}\end{eqnarray} \tag{ 244 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\mathcal{S}}_{N}}:=\underset{j=0}{\overset{N-1}{\sum}}\, \Delta t\left(\text{i}{{x}_{j}}\frac{{{q}_{j}}-{{q}_{j+1}}}{ \Delta t}-{{{\tilde{\mathcal{Q}}}}_{{{t}_{j}}}}\left({{q}_{j+1}},\text{i}{{x}_{j}}\right)\right). \end{array}\end{eqnarray} \tag{ 245 }$$

Here we again employed the abbreviation (167), i.e.

$$\begin{eqnarray}&&{\rlap{-} \int}_{k}^{l}=\underset{j=k}{\overset{l}{\prod}}\,{{{\int}^{}}_{{{\mathbb{R}}^{2}}}}\frac{\text{d}{{q}_{j}}\text{d}{{x}_{j}}}{2\pi}.\end{eqnarray} \tag{ 246 }$$

Moreover, the initial condition $p\left({{t}_{0}},n|{{t}_{0}},{{n}_{0}}\right)={{\delta}_{n,{{n}_{0}}}}$ is again trivially fulfilled for N  =  0. Upon taking the continuous-time limit $N\to \infty$, so that $\Delta t\to 0$, the forward path integral representation of the master equation follows as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & p\left(t,n|{{t}_{0}},{{n}_{0}}\right)={{\langle n{{|}_{t,q(t)}}|g\left(t|{{t}_{0}},{{n}_{0}}\right)\rangle}_{q(t)}} \end{array}\end{eqnarray} \tag{ 247 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}}\text{with}\,|g\left(t|{{t}_{0}},{{n}_{0}}\right){{\rangle}_{q(t)}}={\rlap{-} \int}_{[{{t}_{0}}}^{t)}{{\text{e}}^{-\mathcal{S}}}\,|{{n}_{0}}{{\rangle}_{{{t}_{0}},q\left({{t}_{0}}\right)}} \end{array}\end{eqnarray} \tag{ 248 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \text{and}~\mathcal{S}:={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\left[\text{i}x{{\partial}_{-\tau}}q-{{{\tilde{\mathcal{Q}}}}_{\tau}}\left(q,\text{i}x\right)\right]. \end{array}\end{eqnarray} \tag{ 249 }$$

The limit ${\rlap{-} \int}_{[{{t}_{0}}}^{t)}:={{\lim}_{N\to \infty}}{\rlap{-} \int}_{0}^{N-1}$ now involves integrations over x(t₀) and q(t₀), but not over x(t) and q(t).

The forward path integral (248) can, for example, be used to derive the generating function of the generic linear process $A\to l\,A$ with rate coefficient ${{\mu}_{\tau}}$ (and $l\in {{\mathbb{N}}_{0}}$). For that purpose, we choose the basis function as $|n{{\rangle}_{\tau,q}}={{q}^{n}}$ so that $|g\rangle$ coincides with the ordinary generating function. Since the basis function does not depend on time, (64) alone specifies the transition operator

$$\begin{eqnarray}&&{{\tilde{\mathcal{Q}}}_{\tau}}\left(q,{{\partial}_{q}}\right)={{\mu}_{\tau}}\left({{q}^{l}}-q\right){{\partial}_{q}}\end{eqnarray} \tag{ 250 }$$

of the flow equation ${{\partial}_{\tau}}|g(\tau |\centerdot )\rangle ={{\tilde{\mathcal{Q}}}_{\tau}}\left(q,{{\partial}_{q}}\right)|g(\tau |\centerdot )\rangle$. Consequently, the action

$$\begin{eqnarray}&&\mathcal{S}={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\text{i}x\left[{{\partial}_{-\tau}}q-{{\mu}_{\tau}}\left({{q}^{l}}-q\right)\right]\end{eqnarray} \tag{ 251 }$$

is linear in $\text{i}x$. To evaluate the path integral (248), it helps to reconsider the discrete-time approximation

$$\begin{eqnarray}&&{{\mathcal{S}}_{N}}=\underset{j=0}{\overset{N-1}{\sum}}\,\text{i}{{x}_{j}}\left({{q}_{j}}-\left({{q}_{j+1}}+{{\mu}_{{{t}_{j}}}}\left(q_{j+1}^{l}-{{q}_{j+1}}\right) \Delta t\right)\right).\end{eqnarray} \tag{ 252 }$$

of the action. Upon integrating over all the q_j -variables and taking the limit $\Delta t\to 0$, one obtains the generating function

$$\begin{eqnarray}&&|g\left(t|{{t}_{0}},{{n}_{0}}\right){{\rangle}_{q(t)}}=q{{\left({{t}_{0}}\right)}^{{{n}_{0}}}}~\end{eqnarray} \tag{ 253 }$$

with $q(\tau )$ solving ${{\partial}_{-\tau}}q={{\mu}_{\tau}}\left({{q}^{l}}-q\right)$. The unique real solution of this equation with final value q(t) reads

$$\begin{eqnarray}&&q(\tau )=\frac{q(t)}{{{\left[q{{(t)}^{l-1}}+{{\text{e}}^{(l-1){\int}_{\tau}^{t}\text{d}s\,{{\mu}_{s}}}}\left(1-q{{(t)}^{l-1}}\right)\right]}^{1/(l-1)}}}.\,\end{eqnarray} \tag{ 254 }$$

For the linear growth, or Yule-Furry [119, 413], process $A\to 2\,A$ (l  =  2), the inverse transformation (247) casts the generating function into

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & p\left(t,n|{{t}_{0}},{{n}_{0}}\right)=\left({{\frac{1}{n!}\partial _{q(t)}^{n}|g\left(t|{{t}_{0}},{{n}_{0}}\right)\rangle}_{q(t)}}\right){{|}_{q(t)=0}} \end{array}\end{eqnarray} \tag{ 255 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & =\left(\begin{array}{c} n-1 \\ n-{{n}_{0}} \end{array}\right){{\left({{\text{e}}^{-{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\mu}_{\tau}}}}\right)}^{{{n}_{0}}}}{{\left(1-{{\text{e}}^{-{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\mu}_{\tau}}}}\right)}^{n-{{n}_{0}}}} \end{array}\end{eqnarray} \tag{ 256 }$$

for n  >  0 and into ${{\delta}_{0,{{n}_{0}}}}$ for n  =  0. The above distribution has the form of a negative Binomial distribution with probability of success ${{\text{e}}^{-{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\mu}_{\tau}}}}$, number of failures n  −  n₀, and number of successes n₀ [414]. The mean value ${{\text{e}}^{{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\mu}_{\tau}}}}{{n}_{0}}$ of the marginalized distribution (256) grows exponentially with time as long as ${{\mu}_{\tau}}>0$, and so does its variance $\left({{\text{e}}^{{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\mu}_{\tau}}}}-1\right)\,{{\text{e}}^{{\int}_{{{t}_{0}}}^{t}\text{d}\tau \,{{\mu}_{\tau}}}}{{n}_{0}}$.

For the linear decay process $A\to \varnothing$ (l  =  0), the inverse transformation in (255) instead recovers the Binomial distribution (134).

The above procedure can be generalized to processes whose transition operator is of the form

$$\begin{eqnarray}&&{{\tilde{{\mathcal{Q}}}\,}_{\tau}}\left(q,{{\partial}_{q}}\right)={{\alpha}_{\tau}}(q){{\partial}_{q}}\,+\,\frac{1}{2}{{\beta}_{\tau}}(q)\partial _{q}^{2}\,+\,{{\mathcal{P}}_{\tau}}\left(q,{{\partial}_{q}}\right).\end{eqnarray} \tag{ 257 }$$

The derivation proceeds analogously to section 4.3 and appendix D. The probability distribution is thereby expressed as an average over the paths of a stochastic differential equation proceeding backward in time. The merit of such a representation remains to be explored. In the following, we briefly outline how the procedure is applied to processes with continuous sample paths.

In section 4.4, we explained how the forward master equation (184) of a process with a continuous state space can be written in terms of the (forward) Kramers–Moyal expansion

$$\begin{eqnarray}&&{{\partial}_{\tau}}p(\tau,q|\centerdot )=\underset{m=1}{\overset{\infty}{\sum}}\,\frac{{{(-1)}^{m}}}{m!}\partial _{q}^{m}\left[\mathcal{M}_{\tau}^{(m)}(q)p(\tau,q|\centerdot )\right]\end{eqnarray} \tag{ 258 }$$

with initial condition $p\left(\tau,q|{{t}_{0}},{{q}_{0}}\right)=\delta \left(q-{{q}_{0}}\right)$. Here we changed the letter from x to q to keep in line with the notation used in sections 5.1 and 5.2. If the Kramers–Moyal expansion stops with its second summand, it coincides with the (forward) Fokker–Planck equation

$$\begin{eqnarray}&&{{\partial}_{\tau}}p(\tau,q|\centerdot )=-{{\partial}_{q}}\left[\mathcal{M}_{\tau}^{(1)}(q)p\right]+\frac{1}{2}\partial _{q}^{2}\left[\mathcal{M}_{\tau}^{(2)}(q)p\right].\end{eqnarray} \tag{ 259 }$$

To apply the procedure from the previous section to this equation, it needs to be brought into the form ${{\partial}_{\tau}}p(\tau,q|\centerdot )={{\tilde{\mathcal{Q}}}_{\tau}}\left(q,{{\partial}_{q}}\right)p(\tau,q|\centerdot )$ with a normal-ordered transition operator ${{\tilde{\mathcal{Q}}}_{\tau}}\left(q,{{\partial}_{q}}\right)$. It is readily established that this operator can be expressed in the form of (257) with the coefficients

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\alpha}_{\tau}}(q):=-\mathcal{M}_{\tau}^{(1)}+{{\partial}_{q}}\mathcal{M}_{\tau}^{(2)}, \end{array}\end{eqnarray} \tag{ 260 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\beta}_{\tau}}(q):=\mathcal{M}_{\tau}^{(2)},\quad \text{and} \end{array}\end{eqnarray} \tag{ 261 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{P}}_{\tau}}(q):=-{{\partial}_{q}}\mathcal{M}_{\tau}^{(1)}+\frac{1}{2}\partial _{q}^{2}\mathcal{M}_{\tau}^{(2)}. \end{array}\end{eqnarray} \tag{ 262 }$$

The Fokker–Planck equation now has the same form as the flow equation obeyed by the generating function in section 5.1. Therefore, we can follow the steps in that section to represent the solution of the Fokker–Planck equation by the path integral

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} p\left(t,q|{{t}_{0}},{{q}_{0}}\right)={\rlap{-} \int}_{[{{t}_{0}}}^{t)}\,{{\text{e}}^{-\mathcal{S}}}\,\delta \left({{q}_{0}}-q\left({{t}_{0}}\right)\right){{|}_{q(t)=q}} \end{array}\end{eqnarray} \tag{ 263 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{with}~\mathcal{S}:={\int}_{{{t}_{0}}}^{t}\text{d}\tau \,\left[\text{i}x{{\partial}_{-\tau}}q-{{{\tilde{\mathcal{Q}}}}_{\tau}}\left(q,\text{i}x\right)\right]. \end{array}\end{eqnarray} \tag{ 264 }$$