One simple question we can ask is what is the distribution of the population in the stationary limit $t \to \infty$. Such distribution must be invariant over time translations, thus it is obtained by setting Eq. (5.5) to zero. We obtain

$$0 = j_{\infty}(1, 0) = j_{\infty}(2, 1) = \dots \qquad (5.6)$$

leading to the principle of detailed balance, which prescribes

$$k_{+}p_{\infty}(X) = k_{-}(X + 1)p_{\infty}(X + 1) \qquad (5.7)$$

yielding

$$p_{\infty}(X) = \alpha \frac{\left(\frac{k_{+}}{k_{-}}\right)^{X}}{X!}, \qquad (5.8)$$

where $\alpha = p(0)$ is fixed by normalization $\sum_{X} p_{\infty}(X) = 1$. Noting that $x_{\infty} \equiv x(\infty) = k_{+}/k_{-}$ is the unique fixed point $\lim_{t \to \infty} x(t)$ of the corresponding deterministic rate equation (see Lecture 4)

$$\frac{dx(t)}{dt} = k_{+} - k_{-}x(t) \qquad (5.9)$$

we find that $\alpha = e^{-x_{\infty}}$. Therefore Eq. (5.8) is just the the Poisson distribution with parameter $x_{\infty}$.

We can plot the stationary distribution and compare it to the average time spent at state X via a Monte-Carlo Markov Chain / Doob-Gillespie algorithm (see Lecture 1). A simple imperative approach with basic Python is in Fig. 5.1 Its output Fig. 5.2 compares the total time spent at states vs. the total number of times a state is visited against the true solution Eq. (5.8), showing that only the former is correct (confusion between these two observables is a frequent early mistake in such kind of coding).

<figure><img src="image_44.png" /><figcaption>Figure 5.2: The output of algorithm 5.1. Top histogram: time spent at a state, normalized; Bottom histogram: number of times a state is visited, normalized; Bullets: values of the Poisson distribution.</figcaption></figure>
The Doob-Gillespie algorithm produces the exact statistics for such kind of observables. However, notice that the estimation of $p_{t \to \infty}$ should be done on $n$ independent samples at large enough $t$. Here instead we are integrating along a single realization, hoping that sample averages equal time averages (ergodic principle). This assumption is not trivial and it lends itself to some errors that in principle should be evaluated or dealt with, in particular [73]:

- Relaxation error: we started counting the time spent at a state at $t = 0$, but in principle we should wait long-enough to let the system relax, and then we should start counting for long-enough. But how long is long?

- Autocorrelation error: the probability of being back to some state right after having been there is larger than being at a more distant state. This means that by counting all of the states visited we are favouring more probable ones ("the rich get richer"). To avoid this we should put a long-enough time between one sample and the next. But again, how long is long?

Of course all such petty issues are washed away by machine learning$^{17}$.

**Additional.** Equation (5.5) is the backward Kolmogorov (Fokker-Panck) equation. One interesting question we can ask is: is there an analog of Brownian motion (Langevin equation) for chemical processes? That is, an equation directly expressed in terms of $\hat{X}(t)$?

Notice that we can write

$$ \hat{X}(t) = \hat{X}_0 + \hat{N}_+(t) - \hat{N}_-(t) \qquad (5.10) $$

where $\hat{X}_0$ is sampled with probability $p_0$ and $\hat{N}_\pm(t)$ are respectively the number of times reaction $\pm$ occurs.

Now let $\hat{N}(t)$ be a unit Poisson process with distribution $p(\hat{N}(t) \equiv \mathcal{N}) = e^{-t}t^\mathcal{N}/\mathcal{N}!$. Clearly the forward reaction is a Poisson process

$$ \hat{N}_+(t) = \hat{N}(k_+t) = \hat{N} \left( \int_0^t k_+ ds \right) = \hat{N} \left( \int_0^t k_+ [\emptyset](s) ds \right) \qquad (5.11) $$

where we made some dull manipulations. The backward reaction is not a homogeneous Poisson process because its rate depends on the number of molecules present in the reactor. However, inspired by the latter expression we are tempted to write

$$ \hat{N}_-(t) = \hat{N} \left( \int_0^t k_- \hat{X}(s) ds \right). \qquad (5.12) $$

This given, a wrong argument is that, given that the mean of a Poissonian with parameter $\lambda$ is $\lambda$, then the mean of Eq. (5.10) is

$$ \langle \hat{X}(t) \rangle = \langle \hat{X}_0 \rangle + k_+ t - \int_0^t k_- \langle \hat{X}(s) \rangle ds \qquad (5.13) $$

<footnote>$^{17}$I am being overtly ironic about the way the machine learning hype is propagating into scientific narratives. Inference based on high-dimensional nonlinear functions is very fine-tuned to specific tasks, and in my opinion it is not poised as comes to the creativeness and rigour of scientific inquiry. However, there might exist interesting specific applications of such algorithms for CRNs, see passages of Ref. [74] for a review.</footnote>
where $\langle \cdot \rangle$ is the average over samples. This equation in fact is just the integral solution of Eq. (5.9), which on the other hand can be obtained by multiplying by $X$ and summing over $X$ in the master equation Eq. (5.5).

This latter line of derivations is only true for linear systems, but it can be made more correct in terms of the mode in place of the mean in the large-scale limit, see e.g. Ref. [75] for a thorough explanation of random-time-change Poisson processes and their large-scale limits, that we will briefly introduce in Sec. 5.6, by means of a scaling parameter $\Omega$. Physicists often talk about “mean-field” behaviour by replacing

$$ \langle X_t^2 \rangle \rightarrow \langle X_t \rangle^2 \qquad (5.14) $$

which does reproduce the deterministic dynamics. However, this is a strong assumption that implies that the mean behaviour is representative of the process, in some sense. In fact, CRNs can display a very different behaviour among their stochastic and their deterministic counterparts. For example in Ref. [76] it has been shown that a so-called strongly endotactic CRN (whose deterministic dynamics falls into a closed domain) can have stochastic escape routes to infinity (and therefore no stationary distribution), while the deterministic dynamics stays confined in a compact set. This is a case where the mean and mode of the process deviate wildly (although it requires irreversible reactions).

## 5.3 Fishing for Poisson-like

Now consider

$$ \emptyset \rightleftharpoons 2X \qquad (5.15) $$

creating a complex $C_2 = 2X$ made of two interacting copies of species X out of the empty complex $C_1 = \emptyset$. The population space is

$$ 0 \rightarrow 1 \rightarrow 2 \rightarrow 3 \rightarrow 4 \rightarrow \dots \qquad (5.16) $$

There are two disconnected components $\mathcal{C}_{even} = 2N$ and $\mathcal{C}_{odd} = 2N + 1$ (even and odd integers), called stoichiometric compatibility classes. Correspondingly, the parity of $\hat{X}(t)$ is conserved, given an initial even/odd number of molecules. Notice that this (broken) symmetry is not present in the corresponding deterministic rate equation $dx(t)/dt = k_+-k_-x(t)^2$ (in Quantum Field Theory similar things are called anomalies). The stationary distribution now reads

$$ p_{\infty}(X) = \alpha \frac{x_{\infty}^{X}}{X!} [p_{\text{even}} \mathbf{1}_{\mathcal{C}_{\text{even}}}(X) + p_{\text{odd}} \mathbf{1}_{\mathcal{C}_{\text{odd}}}(X)] \qquad (5.17) $$

where $\mathbf{1}_{\mathcal{C}}(a)$ is the indicator function<fnref n="18" /> and $p_{\text{even/odd}}$ are the normalized probabilities that the initial state is either even or odd, and $\alpha$ is again fixed by normalization yielding

$$ \alpha = [p_{\text{even}} \cosh x_{\infty}) + p_{\text{odd}} \sinh x_{\infty}]^{-1}. \qquad (5.18) $$

Notice that Eq. (5.17) looks Poissonian, but it is not because for a stochastic variable to be Poissonian it has to have the same expression over all the integers. Here they are split in two, and therefore we talk of Poisson-like or Poisson-form distribution.

<footnote n="18">¹⁸Namely, 1 if $a \in \mathcal{A}$ else 0. Quite conveniently, in the Doob-Gillespie algorithm shown in Fig. 5.1 the indicator function has not to be implemented explicitly because negative populations are never reached, given that the reaction $0 \rightarrow -1$ always has zero rate</footnote>
## 5.4 Moieties and product-form

Now consider

$$X_2 + X_3 \rightleftharpoons 2X_1. \qquad (5.19)$$

Each single realization of the reaction is a random event due to the collision of two molecules, e.g. of X₂ and X₃ forward, or of X₁ with itself backwards. MAK prescribes

$$\begin{aligned} r_+(X_1, X_2, X_3) &= k_+ X_2 X_3 \\ r_-(X_1, X_2, X_3) &= k_- X_1 (X_1 - 1) \end{aligned} \qquad (5.20)$$

where the “-1” is due to the fact that there is one less molecule to bounce into (this correction is lost at the deterministic level).

The stoichiometric matrix

$$S = \begin{pmatrix} +2 \\ -1 \\ -1 \end{pmatrix} \qquad (5.21)$$

has a two-dimensional left-null space (obviously). As a meaningful basis of null co-vectors such that $\ell^\top S = 0$ we prefer to choose ones that have small positive integer entries, for example $\ell_{\text{mass}}^\top = (1, 1, 1)$ (total mass conservation) and $\ell_{\text{moie}}^\top = (1, 2, 0)$. The former corresponds to mass conservation. Chemists call the latter a moiety, that is, a part of a molecule that has a name because it is also part of other molecules, in this case the two “disks” in X₂ and the one “disk” in X₁:

<figure><img src="image_45.png" /></figure>

(mass is the ultimate moiety given that it is part of all molecules. Notice that we could have also chosen $\ell_{\text{moie}}^\top = (1, 0, 2)$, in which case the moiety would be the disk with a cavity). In general finding right or left null vectors of a stoichiometric matrix of chemical meaning is not just a purely theoretical exercise, but it entails some understanding of the chemistry behind the formalism. See Ref. [74] for some interesting recent advancements on the interpretation of moieties.

The associated deterministic system reads

$$\begin{aligned} \frac{d}{dt}x_2(t) &= k_-x_1(t)^2 - k_+x_2(t)x_3(t) \\ \lambda_{\text{mass}} &= x_1(t) + x_2(t) + x_3(t) \\ \lambda_{\text{moie}} &= x_1(t) + 2x_2(t) \end{aligned} \qquad (5.22)$$

where $\lambda_{\text{mass}}$, $\lambda_{\text{moie}}$ are the deterministically conserved quantities. Given these values, the above system has a unique fixed point $x_\infty = x(\infty) = \lim_{t \to \infty} x(t)$. Notice that by the first equation the ratio

$$\frac{k_+}{k_-} = \frac{x_1(\infty)^2}{x_2(\infty)x_3(\infty)} \qquad (5.23)$$

does not depend explicitly on the values of the conserved quantities. However, in general each initial population will converge to a different fixed point parametrized by $\lambda_{\text{mass}}$, $\lambda_{\text{moie}}$.
The stochastic population space is composed of compatibility classes analogous to the even/odd classes in the previous example according to the values

$$L_{\text{mass}} = X_1 + X_2 + X_3 \\ L_{\text{moie}} = X_1 + 2X_2. \qquad (5.24)$$

Visually, for $L_{\text{mass}} = 0$ (red dot), $L_{\text{mass}} = 1$ (green dots), $L_{\text{mass}} = 2$ (blue dots):

<figure><img src="image_46.png" /></figure>

For higher masses and moieties each compatibility class is at most a finite interval

$$(\bar{X}_1, \bar{X}_2, \bar{X}_3) = (\bar{X}_1 + 2, \bar{X}_2 - 1, \bar{X}_3 - 1) = \cdots = (\bar{X}_1 + 2N, \bar{X}_2 - N, \bar{X}_3 - N)$$

ranging from either $\bar{X}_1 = 0, 1$ for some $N$ such that either $\bar{X}_2 = N$ or $\bar{X}_3 = N$. On each such subspace we can again use the principle of detailed balance

$$k_-(X_1 + 2)(X_1 + 1)p_\infty(X_1 + 2, X_2 - 1, X_3 - 1) = k_+X_2X_3p_\infty(X_1, X_2, X_3) \qquad (5.25)$$

to propagate the solution to any population compatible with $(\bar{X}_1, \bar{X}_2, \bar{X}_3)$. We find for example

$$\begin{aligned} p_\infty(\bar{X}_1 + 4, \bar{X}_2 - 2, \bar{X}_3 - 2) &= \frac{k_+(\bar{X}_2 - 1)(\bar{X}_3 - 1)}{k_-(\bar{X}_1 + 4)(\bar{X}_1 + 3)} p_\infty(\bar{X}_1 + 2, \bar{X}_2 - 1, \bar{X}_3 - 1) \\ &= \left(\frac{k_+}{k_-}\right)^2 \frac{(\bar{X}_2 - 1)\bar{X}_2(\bar{X}_3 - 1)\bar{X}_3}{(\bar{X}_1 + 4)(\bar{X}_1 + 3)(\bar{X}_1 + 2)(\bar{X}_1 + 1)} p_\infty(\bar{X}_1, \bar{X}_2, \bar{X}_3) \\ &= \left[\frac{x_1(\infty)^2}{x_2(\infty)x_3(\infty)}\right]^2 \frac{\bar{X}_1!\bar{X}_2!\bar{X}_3!}{(\bar{X}_1 + 4)!(\bar{X}_2 - 2)!(\bar{X}_3 - 2)!} p_\infty(\bar{X}_1, \bar{X}_2, \bar{X}_3) \end{aligned} \qquad (5.26)$$

where in the second passage we used Eq. (5.23) and played a trick with factorials. We find

$$p_\infty(\bar{X}_1 + 4, \bar{X}_2 - 2, \bar{X}_3 - 2) = \alpha \frac{x_1(\infty)^{\bar{X}_1+4} x_2(\infty)^{\bar{X}_2-2} x_3(\infty)^{\bar{X}_3-2}}{(\bar{X}_1+4)! (\bar{X}_2-2)! (\bar{X}_3-2)!} \qquad (5.27)$$

where we played yet another trick by multiplying and dividing by $x_1(\infty)^{\bar{X}_1} x_2(\infty)^{\bar{X}_2} x_3(\infty)^{\bar{X}_3}$ and defining

$$\alpha = \frac{\bar{X}_1! \bar{X}_2! \bar{X}_3!}{x_1(\infty)^{\bar{X}_1} x_2(\infty)^{\bar{X}_2} x_3(\infty)^{\bar{X}_3}} p_\infty(\bar{X}_1, \bar{X}_2, \bar{X}_3). \qquad (5.28)$$
Proceeding like above once again to the next population one would find

$$p_{\infty}(\bar{X}_1 + 6, \bar{X}_2 - 3, \bar{X}_3 - 3) = \alpha \frac{x_1(\infty)^{\bar{X}_1+6}}{(\bar{X}_1+6)!} \frac{x_2(\infty)^{\bar{X}_2-3}}{(\bar{X}_2-3)!} \frac{x_3(\infty)^{\bar{X}_3-3}}{(\bar{X}_3-3)!}, \quad (5.29)$$

and so on. This kind of distribution is called product-form Poisson-like. Notice, quite interestingly, that while a solution of the related deterministic system enters the distribution, we could express this latter just in terms of $k_+/k_-$, as exemplified by the first passage of the derivation of Eq. (5.26). Thus the distribution does not depend on the deterministic conserved quantities $\lambda_{\text{mass}}$, $\lambda_{\text{moie}}$, and one can plug in any solution of the deterministic system. The computation of $\alpha$, containing all of the correlation between the species, is far from trivial.

## 5.5 One trivial cycle

Consider

$$\begin{align}
\emptyset &\stackrel{1}{\rightleftharpoons} X \tag{5.30} \\
\emptyset &\stackrel{2}{\rightleftharpoons} X
\end{align}$$

where we assume that we can distinguish the two reactions. The population space is

$$0 \mathop{\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\dots} 1 \mathop{\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\dots} 2 \mathop{\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\dots} 3 \dots \quad (5.31)$$

By lumping together the two reaction rates into an effective $r_{\pm}(X) = r_{\pm 1}(X) + r_{\pm 2}(X)$ we re-obtain the first example in Sec. 5.2 with reaction rate constants $k_{\pm} = k_{\pm 1} + k_{\pm 2}$, thus from a kinetic point of view this example has nothing to offer and the stationary population is the Poissonian with parameter $(k_{\pm 1} + k_{\pm 2})/(k_{-1} + k_{-2})$.

Here we are rather interested in dynamical behaviour. Because we can distinguish two different mechanisms, when the system has relaxed to the stationary distribution there will be a net stationary circulation of currents

$$0 \mathop{\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\dots} 1 \mathop{\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\dots} 2 \mathop{\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\leftarrow\dots} 3 \dots \quad (5.32)$$

In Lecture 2 the cycle affinities were defined as the log-product of ratio of forward and backward rates of a process along any cycle. Here we would have to consider all of the infinite cycles, but we have that populations cancel out and the affinity is the same for all cycles

$$\log \frac{k_{+1}Xk_{-2}(X+1)}{k_{+2}(X+1)k_{-1}X} = \log \frac{k_{+1}k_{-2}}{k_{+2}k_{-1}} =: \mathcal{A}. \quad (5.33)$$

For a chemical interpretation, we can think of $\emptyset$ to actually be species that are chemostatted, e.g.

$$\begin{align}
Y_1 + Y_2 &\stackrel{1}{\rightleftharpoons} X \tag{5.34} \\
Y_3 &\stackrel{2}{\rightleftharpoons} X
\end{align}$$

where all $[Y_1]$, $[Y_2]$, $[Y_3]$ are assumed to be kept fixed in time (see Lecture 4). Notice that this only amounts to a rewriting the rate constants $k_{+1} = k'_{+1}[Y_1][Y_2]$ and
$k_{+2} = k'_{+2}[Y_3]$ on the assumption that when the chemostatted species are held to 1 (or to some standard chemical potential, on which we do not delve into), detailed balance holds $k'_{+1}/k'_{-1} = k'_{+2}/k'_{-2}$. In this case the affinity reads

$$ \mathcal{A} = \log[Y_1] + \log[Y_2] - \log[Y_3] \qquad (5.35) $$

and it dictates the direction of transport of the chemostatted chemicals through the environment as a cycle is performed

$$ Y_1 + Y_2 \rightarrow Y_3, \qquad (5.36) $$

which is reminiscent of the school's logo.

Now, defining the time-integrated current $\hat{J}(t) = \hat{N}_+(t) - \hat{N}_-(t)$ as the total number of reactions +1 minus the total number of reactions -1 occuring up to time $t$, that is, the number of all the transitions of the kind

$$ 0 \underset{1}{\overset{-1}{\rightleftharpoons}} 1 \underset{2}{\overset{-1}{\rightleftharpoons}} 2 \underset{3}{\overset{-1}{\rightleftharpoons}} \dots, \qquad (5.37) $$

minus that of transitions

$$ 0 \underset{1}{\overset{-1}{\rightleftharpoons}} 1 \underset{2}{\overset{-1}{\rightleftharpoons}} 2 \underset{3}{\overset{-1}{\rightleftharpoons}} \dots \qquad (5.38) $$

<figure><img src="image_47.png" /><figcaption>Figure 5.3: For $k_{+1} = 1, k_{-1} = 2, k_{+2} = 3, k_{-2} = 4$, initial population $\hat{X}(0) = 0$, number of samples 50000, for three values of the final time, histograms of the currents are their Gaussian fit, and log-ratios of positive to negative currents' probabilities.</figcaption></figure>
one question we can ask is whether the following well-known fluctuation relation (see e.g. [77]) is satisfied:

$$ \log \frac{p(\hat{J}(t) \equiv \mathcal{J})}{p(\hat{J}(t) \equiv -\mathcal{J})} \stackrel{?}{=} \mathcal{A}\mathcal{J}. \qquad (5.39) $$

We plot the left-hand side quantity in the bottom panel of Fig. 5.3, for various values of time. First we notice that at short times there is a systematic bias due to dependence on the initial population. So we need to wait long enough to claim that a fluctuation relation is satisfied. But even then, is this relation of statistical significance<fnref n="19" />? Notice in fact that, from the first panel in Fig. 5.3, the distribution of the currents appears to be reasonably Gaussian, as we expect by the onset of the Central Limit Theorem. Thus the linearity here appears to be just a statistical artifact, and it is hard to claim a truly nonequilibrium phenomenon.

One way out of this conundrum may be to sample from an initial distribution such that the fluctuation relation is satisfied at all times, if it exists. One might be tempted to speculate that such initial distribution should be the stationary one. The top panel of Fig. 5.4 hints that the initial distribution that yields the fluctuation relation at all times is the one obtained by first letting the system relax to the equilibrium distribution of reaction 2 only, such that the population space is

<figure><img src="image_48.png" /><figcaption>Figure 5.4: Top panel: Same rates as above, number of samples 5000000, final time $t = .4$, the log-ratio of forward to backward probabilities for the initial population $\hat{X}_0$: (blue) equal to 0 with certainty; (yellow) sampled with Poissonian with parameter $k_{+2}/k_{-2}$; (green) sampled with Poissonian with parameter $(k_{+2} + k_{+1})/(k_{-1} + k_{-2})$. Bottom panel: Same as above, but with stopping time the total number of reactions $\hat{N}_{+1} + \hat{N}_{-1} = 8$.</figcaption></figure>

<footnote n="19">¹⁹How can we claim to observe a nonequilibrium behaviour? For this to be, we would need to develop some sound statistical test of the fluctuation relation not being explained by a null hypothesis representing the linear regime.</footnote>
$$ \underbrace{0}_{2} \overbrace{1}_{2} \underbrace{2}_{2} \overbrace{3}_{2} \dots, \qquad (5.40) $$

and then at time $t = 0$ turn on the other transitions and count the currents. The general proof for generic graphs is in [78], but this does not generally translate to CRNs so well.

A second way out is to consider fluctuation relations at stopping time, as in Ref. [3]. For example, in the bottom panel of Fig. 5.4 we stop the process exactly after a fixed number reactions $\hat{N}_{+1} + \hat{N}_{-1}$ have occurred, and it appears that convergence of the fluctuation relation is slightly faster.

## 5.6 One less trivial cycle: Schlögl model

Take

$$ \begin{gathered} 0 \underset{2}{\overset{1}{\rightleftharpoons}} X \\ 2X \underset{2}{\overset{2}{\rightleftharpoons}} 3X \end{gathered} \qquad (5.41) $$

with asymmetric population space

$$ \underbrace{0}_{2} \overbrace{1}_{2} \underbrace{1}_{2} \overbrace{2}_{2} \overbrace{3}_{2} \dots \qquad (5.42) $$

To compute the stationary distribution we can lump together reactions to obtain the effective rates

$$ r_{+}(X) = k_{+1}\Omega + k_{+2}\frac{X(X-1)}{\Omega}, \quad \text{for } X \ge 1 \qquad (5.43) $$

$$ r_{-}(X) = k_{-1}X + k_{-2}\frac{X(X-1)(X-2)}{\Omega^2}, \quad \text{for } X \ge 2 $$

where we introduced a scaling parameter $\Omega$ to approach the deterministic limit (the reasoning behind the scaling is that by rescaling $X = \Omega x$ we want to reproduce a deterministic equation such as Eq. (5.9), therefore rate constants need to scale consistently, see also Lecture 8). Notice that these latter lumped reaction rates are not MAK.

State $X = 2$ is pivotal so we can compute the stationary distribution in the same way as above by propagating from $\alpha = p_{\infty}(2)$:

$$ p_{\infty}(0) = \alpha \left( \frac{k_{-1}}{k_{+1}\Omega} \right)^2 \qquad (5.44) $$

$$ p_{\infty}(1) = \alpha \frac{k_{-1}}{k_{+1}\Omega} \qquad (5.45) $$

$$ p_{\infty}(X) = \alpha \prod_{X>2} \frac{k_{+1}\Omega + k_{+2}X(X-1)/\Omega}{k_{-1}X + k_{-2}X(X-1)(X-2)/\Omega^2}, \quad \text{for } X > 2. \qquad (5.46) $$

This distribution is definitely not Poisson-like (see plots in Fig. 5.5). In fact, this is a well-studied toy model (see e.g. [77,79,80]) due to Schlögl to study bistability and critical points. The Schlögl model is a special case of a the broader class of networks that have non-zero deficiency, a concept that we will explain below in greater detail. For now, just be aware that in this case the deficiency is $\delta = 1$ because the system obviously has a cycle, but the network of complexes depicted in Eq. (5.41) has no “visible” cycle. Eventually we can make the cycle visible by turning from a representation in terms of a graph to one in terms of a hypergraph (aka Petri net):
<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /><figcaption>Figure 5.5: The stationary distribution for various values of the volume, for $k_{+1} = 0.5$, $k_{-1} = 3$, $k_{+2} = 4.6$, $k_{-2} = 1$.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Figure 5.6: The mean (blue more regular curve) and most probable (yellow more zig-zagged curve) entropy production rate for the simple model (on the left) with rescaled rates $k_+ \rightarrow k_+\Omega$, $k_- \rightarrow k_-\Omega$ and the Schlögl model (on the right), $k_{+1} = 0.5$, $k_{-1} = 3$, $k_{+2} = 1$ as a function of $k_{-2}$, for four different values of the volume $\Omega = 4, 4^2, 4^3, 4^4$ (top to bottom), $\dot{X}(0) = 0$, final time $t = 10.0$ and $n = 10$ samples.</figcaption></figure>
Notice that this space is very symmetrical, and that it is basically an infinite copy-paste of the above network of complexes 5.48. This is not always the case, as already observed in the SchlGLYPH<suppress>l model.

The stoichiometric matrix reads

$$
\mathbf{S} = \begin{array}{c} \begin{array}{c} \begin{array}{ccc} 1 & 2 & 3 \end{array} \\ \begin{array}{ccc} X_1 & X_2 & X_3 \end{array} \end{array} \begin{pmatrix} +2 & -2 & 0 \\ -1 & +1 & 0 \\ -1 & 0 & +1 \end{pmatrix} = \mathbf{KD} \end{array} \qquad (5.50)
$$

where

$$
\mathbf{D} = \begin{array}{c} \begin{array}{c} \begin{array}{ccc} 1 & 2 & 3 \end{array} \\ \begin{array}{ccc} C_1 & C_2 & C_3 \end{array} \end{array} \begin{pmatrix} -1 & 0 & +1 \\ +1 & -1 & 0 \\ 0 & +1 & -1 \end{pmatrix} \end{array} \qquad (5.51)
$$

is the incidence matrix of the complete graph whose vertices are the complexes and

$$
\mathbf{K} = \begin{array}{c} \begin{array}{c} \begin{array}{ccc} C_1 & C_2 & C_3 \end{array} \\ \begin{array}{ccc} X_1 & X_2 & X_3 \end{array} \end{array} \begin{pmatrix} 0 & 1 & 1 \\ 2 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix} \end{array} \qquad (5.52)
$$

is sometimes referred to as Kirchhoff matrix, telling which species belong to which com-
plexes. The quantity

$$
\delta = \dim \ker \mathbf{S} - \dim \ker \mathbf{D} \qquad (5.53)
$$

is the deficiency briefly mentioned above, where for a graph dim ker S is the number of
independent cycles of the graph (see Lecture 2). Then, roughly, the deficiency is the
number of independent stoichiometric cycles that cannot be visualized as cycles in the
graph of complexes (but can eventually be visualized as cycles in the hypergraph).

In this case the deficiency turns out to be zero. Then several important results apply.
In particular [82], the corresponding deterministic system has a unique fixed point x(∞)
subject to the deterministic conservation law ℓ = x₁(∞) + 2x₂(∞) found by solving the
continuity equations

$$
\overbrace{k_{+1}x_2(\infty)x_3(\infty) - k_{-1}x_1(\infty)^2}^{\nu_1} = \overbrace{k_{+3}x_2(\infty) - k_{-3}x_2(\infty)x_3(\infty)}^{\nu_2} = \overbrace{k_{+2}x_1(\infty)^2 - k_{-2}x_2(\infty)}^{\nu_3} \quad (5.54)
$$

where the overbraces define the stationary deterministic currents. Notice that once again
there is no deterministic analogue of the conservation of parity. A proof that such fixed
points are globally attractive is under review since seven years [83].

The stationary distribution of the master equation is found by solving the continuity
equation at each node of the population network

$$
\sum_{\rho \in \mathcal{R}} j_{\infty}(\rho, X_1, X_2, X_3) = 0, \qquad (5.55)
$$
where $\rho$ labels reactions and $\mathcal{R} = \{\pm 1, \pm 2, \pm 3\}$. Let us write this explicitly:

$$0 = + k_{+1}X_2X_3p_{\infty}(X_1, X_2, X_3) - k_{-1}(X_1 + 2)(X_1 + 1)p_{\infty}(X_1 + 2, X_2 - 1, X_3 - 1) \\ + k_{-1}X_1(X_1 - 1)p_{\infty}(X_1, X_2, X_3) - k_{+1}(X_2 + 1)(X_3 + 1)p_{\infty}(X_1 - 2, X_2 + 1, X_3 + 1) \\ + k_{+2}X_1(X_1 - 1)p_{\infty}(X_1, X_2, X_3) - k_{-2}(X_2 + 1)p_{\infty}(X_1 - 2, X_2 + 1, X_3) \\ + \text{etc.} \qquad (5.56)$$

Clearly, solving this by hand is not feasible, but the important ACK theorem [84] stipulates that, since this system has zero deficiency, its stationary distribution is a product-form Poisson-like distribution with density

$$\alpha \frac{x_1(\infty)^{X_1}}{X_1!} \frac{x_2(\infty)^{X_2}}{X_2!} \frac{x_3(\infty)^{X_3}}{X_3!} \qquad (5.57)$$

on any stoichiometric compatibility class. Let us check that this works by plugging into the above Eq. (5.56). After some tedious work one obtains:

$$0 = [k_{+1}x_2(\infty)x_3(\infty) - k_{-1}x_1(\infty)^2] \frac{x_1(\infty)^{X_1}x_2(\infty)^{X_2-1}x_3(\infty)^{(X_3-1)}}{X_1!(X_2-1)!(X_3-1)!} \\ + \underbrace{[k_{-1}x_1(\infty)^2 - k_{+1}x_2(\infty)x_3(\infty)]}_{-\iota_1} \frac{x_1(\infty)^{X_1-2}x_2(\infty)^{X_2}x_3(\infty)^{X_3}}{(X_1-2)!X_2!X_3!} \\ + \underbrace{[k_{+2}x_1(\infty)^2 - k_{-2}x_2(\infty)]}_{\iota_2} \frac{x_1(\infty)^{X_1-2}x_2(\infty)^{X_2}x_3(\infty)^{X_3}}{(X_1-2)!X_2!X_3!} \\ + \text{etc.} \qquad (5.58)$$

where we recognized the deterministic stationary currents defined in Eq. (5.54). But then notice that since $\iota_1 = \iota_2 = \iota_3$ and the prefactor is the same, the second and third terms cancel out. The same can be proven for all the other terms, which cancel out cyclically.

Notice that, as anticipated, the stationary distribution does not depend on the deterministic conserved quantities, and that we did not need to specify the stoichiometric compatibility class (that would be necessary to compute the normalization). For networks with non-zero deficiency, if one tries to plug in the multi-Poisson-like distribution in the generator of the Chemical Master Equation one obtains that the prefactors do not match as nicely, and terms do not cancel out. The above procedure is, in fact, the blueprint of the original proof of the ACK theorem.

## 5.8 Deficiency and response

First we start with the simple linear network

$$\begin{align*} X_1 &\stackrel{1}{\rightleftharpoons} X_2 \\ X_2 &\stackrel{2}{\rightleftharpoons} X_3 \\ X_3 &\stackrel{3}{\rightleftharpoons} X_1 \end{align*} \qquad (5.59) $$

(clockwise orientation for rate constants $k_+$, counterclockwise for $k_-$). The unique deterministic fixed points $\boldsymbol{x}(\infty)$ such that $\boldsymbol{1} \cdot \boldsymbol{x}(\infty) = \ell_{\text{mass}}$ are found by applying the Markov Chain matrix-tree theorem (see Lecture 1). Having deficiency $\delta = 0$ the stationary distribution is product-form Poisson-like with stoichiometric subspaces labelled by
1 · X = L<sub>mass</sub>. As mentioned above, there is no relation between the choice of ℓ<sub>mass</sub> and L<sub>mass</sub>. The system is nonequilibrium with cycle affinity

$$
\mathcal{A} = \log \frac{k_{+1}k_{+2}k_{+3}}{k_{-1}k_{-2}k_{-3}}. \qquad (5.60)
$$

We assume the above reactions to be hidden from the observer. We want to study how
they affect the behaviour of an additional observable reaction

$$
X_1 + 2X_2 \stackrel{4}{\rightleftharpoons} 3X_3. \tag{5.61}
$$

Overall, the system has stoichiometric matrix

$$
\mathbf{S} = \begin{pmatrix} -1 & 0 & +1 & -1 \\ +1 & -1 & 0 & -2 \\ 0 & +1 & -1 & +3 \end{pmatrix}. \qquad (5.62)
$$

The observable reaction respects the conservation law, and thus by a similar reasoning as
in § 5.2 overall the system has two chemical cycles (right-null vectors). The graph in the
space of complexes is

$$
\left.
\begin{array}{c}
C_1 \\
2 \\
C_3
\end{array}
\right\} \begin{array}{c}
\frac{1}{3} \\
\diagup \\
\diagdown \\
\frac{1}{2}
\end{array}
C_2
C_4 \xrightarrow{4} C_5
\qquad (5.63)
$$

with incidence matrix

$$
\mathbf{D} = \begin{pmatrix} -1 & 0 & +1 & 0 \\ +1 & -1 & 0 & 0 \\ 0 & +1 & -1 & 0 \\ 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & +1 \end{pmatrix}. \qquad (5.64)
$$

The hypergraph instead is

<figure><img src="image_52.png" /></figure>

and there are two independent linear combinations of hyperedges such that an equal
number of arrows go into and out of every node, e.g.

<figure><img src="image_53.png" /></figure>
corresponding to the null eigenvectors of the stoichiometric matrix

$$ \begin{pmatrix} +1 \\ +1 \\ +1 \\ 0 \end{pmatrix}, \quad \begin{pmatrix} 0 \\ -2 \\ +1 \\ +1 \end{pmatrix}. \qquad (5.65) $$

Notice that, differing from the first, the second is not a null vector of the incidence matrix, while it is of the stoichiometric matrix. Therefore the deficiency is $\delta = 1$. While graph cycles admit a simple physical interpretation in terms of stuff that is conserved throughout, hypercycles do not allow such simple visuals. An attempt to construct a dictionary and geometric interpretation is in Ref. [20].

Now let $a = \log k_{+4}/k_{-4}$ be a parameter regulating the intensity of the observable reaction, let $\hat{J}(t)$ be the its time-integrated current, and let

$$ j(a) = \lim_{t \to \infty} \frac{\langle \hat{J}(t) \rangle}{t} \qquad (5.66) $$

be the stationary mean current (notice that in general this is different from the deterministic current $\iota$ because the system has nonzero deficiency). We are interested in the relationship between the integrated current's scaled variance (as a function of $a$)

$$ \kappa(a) := \lim_{t \to \infty} \frac{\langle \hat{J}(t)^2 \rangle - \langle \hat{J}(t) \rangle^2}{t} \qquad (5.67) $$

and the current's response to a perturbation of $a$

$$ \rho(a) := \frac{\partial j(a)}{\partial a}. \qquad (5.68) $$

In the first panel of Fig. 5.7 we plot the mean current $j(a)$ and in the second the response-to-twice-variance $\kappa(a)/2\rho(a)$ for several choices of rate constants and $L_{\text{mass}} = 6$. Different data series correspond to different values of $a$. Different data points along a series correspond to different values of the affinity $\mathcal{A}$. The picture suggests that for the value $a = a^*$ such that the mean current is $j(a^*) = 0$ the following fluctuation-dissipation relation is satisfied

$$ \kappa(a^*) = 2\rho(a^*). \qquad (5.69) $$

This latter is a milestone of close-to-equilibrium statistical mechanics. However, here it would be displayed far from equilibrium ($\mathcal{A} \neq 0$) in states characterized by vanishing observable currents – called stalling instead of equilibrium – whose observable currents vanish, but that have non-observable flows.

**Additional:** The fluctuation-dissipation relation above can be easily derived if for some appropriate initial distribution the following integral fluctuation relation holds:

$$ 1 = \langle e^{(a^*-a)\hat{J}} \rangle (a) \qquad (5.70) $$

for all values of $a$, in which case the exponent is a candidate for a physically meaningful measurement of (partial) entropy production. Notice that in general the average itself depends on $a$. In fact, by taking the first derivative with respect to $a$ and evaluating at
$a = a^*$ we find

$$ 
\begin{aligned}
0 &= \left. \frac{d}{da} \langle e^{(a^*-a)\hat{J}} \rangle (a) \right|_{a=a^*} \\
&= \left[ \left. \frac{\partial}{\partial a} \langle e^{(a^*-a)\hat{J}} \rangle (a) - \langle \hat{J} e^{(a^*-a)\hat{J}} \rangle (a) \right] \right|_{a=a^*} \\
&= \left[ \left. \frac{\partial}{\partial a} 1 - \langle \hat{J} e^{(a^*-a)\hat{J}} \rangle (a) \right] \right|_{a=a^*} \\
&= -\langle \hat{J} \rangle (a^*)
\end{aligned} 
\quad (5.71) $$

where in the second passage we used the integral fluctuation relation above. This confirms that the average current at $a = a^*$ stalls. By taking the second total derivative with

<figure><img src="image_54.png" /><figcaption>Figure 5.7: Top panel: the average current for trajectories stopping at time $t = 50$, 10000 samples, values of the hidden rates $k_{+1} = k_{+2} = k_{+3} = \exp A/6$ and $k_{-1} = k_{-2} = k_{-3} = \exp -A/6$ and of the observable rates $k_{+4} = \exp a/2$, $k_{-4} = \exp -a/2$ as a function of the affinity $A$, for different values of $a$. Bottom panel: under the same conditions, ratio of the response over twice the variance of the current. For this subset of rates the stalling value is $a^* = 0$, which follows from the following argument. First notice that the stationary distribution of the hidden network is independent of the $k$'s. Hence the value of the rate constants $k_{\pm 4}$ for which the observable current vanishes is independent of the other rate constants. So we can choose hidden symmetric rates $k_{+1} = k_{-1}$ etc. such that $A = 0$ for which stalling is equilibrium, and then detailed balance prescribes $k_{+4} = k_{-4}$.</figcaption></figure>
respect to a we find

$$ 
\begin{aligned}
0 &= \frac{d^2}{da^2} \langle e^{(a^*-a)\hat{J}} \rangle (a)|_{a=a^*} \\
&= \frac{d}{da} \left[ \frac{\partial}{\partial a} \langle e^{(a^*-a)\hat{J}} \rangle (a) - \langle \hat{J}e^{(a^*-a)\hat{J}} \rangle (a) \right] \Big|_{a=a^*} \\
&= \left[ \frac{\partial^2}{\partial a^2} \langle e^{(a^*-a)\hat{J}} \rangle (a^*) - 2\frac{\partial}{\partial a} \langle \hat{J}e^{(a^*-a)\hat{J}} \rangle (a) + \langle \hat{J}^2e^{(a^*-a)\hat{J}} \rangle (a) \right] \Big|_{a=a^*} \\
&= -2\frac{\partial}{\partial a} \langle \hat{J} \rangle (a^*) + \langle \hat{J}^2 \rangle (a^*)
\end{aligned} 
\quad (5.72) $$

which is, upon taking the long-time limit, Eq. (5.69).

So, why didn't we just show computationally the validity of the integral fluctuation relation, instead of considering first-order response? In Fig. 5.8 we show that an estimator of the right-hand side of Eq. (5.70) becomes extremely noisy outside of the linear regime. Adding samples tends to linearize the expression around the stalling value $a^* = 0$, but this is a particularly slow process (in terms of samples). In particular, there is a systematic bias towards small values further away from stalling due to the fact that most probably the current is directed in the direction of the effective affinity, making the exponential very small. However, very rare events of a large current opposite to the effective affinity can occur, producing the spikes above the expected value of 1. This clash between most typical behaviour and fat tails makes convergence very problematic (see e.g. [85]).

<figure><img src="image_55.png" /><figcaption>Figure 5.8: For initial population (2,2,2), time $t = 50$, internal rates $k_{+1} = k_{+2} = k_{+3} = 2$, $k_{-1} = k_{-2} = k_{-3} = 1$ and observable rates $k_{+4} = \exp a/2$, $k_{-4} = \exp -a/2$, the integral fluctuation relation estimator as a function of the stalling parameter $a$, and for several values of the total number of samples.</figcaption></figure>
5.9 Considerations

While preparing this lecture a couple of interesting questions arose that may be worth
investigating in the future.

In Sec. 5.5 it was shown that the fluctuation relation can be recovered at all times if an appropriate initial probability distribution is prepared: the question remains open in what sense and under which conditions this distribution has physical sense (that is, it is operationally realizable) or is just a mathematical artifact. There, there is only one current in the system. For more currents, as in Sec. 5.8, there is preliminary work in progress by the Author, extending the formalism of Ref. [86], indicating that in certain circumstances the integral fluctuation relation can be recovered.

Here and there it was suggested that properties of the CRN were linked to some sort
of symmetry of the population space. The language of symmetries being that of groups
(e.g. the population space of $\emptyset \rightleftharpoons 2\mathbb{X}$ must be related to $\mathbb{Z}_2$, the cyclic group of order two)
I wonder whether there could be a systematic characterization of deficiency in terms of
group theory.

It is informally known that the experimental observability of the (integral) fluctuation
relation is problematic, in particular as it comes to deciding whether a system is far from
equilibrium. There exist algorithms that fix the problem computationally (e.g. by such
techniques as cloning, adaptive sampling, etc.), but in my view there are foundational
issues when it comes to experiments with very little sample size. Statistics (delineation
of a null hypothesis, p-value estimation, power study etc.) could help in the choice and
design. Here it was suggested that response out of stalling states may be the proper ground
for making statistically significant claims. Another option is to resort to stopping times
different than the external clock time $t$.
# 6 Dynamical Large Deviations: at long times

William D. Piñeros and Vivien Lecomte.<fnref n="20" /> Rare events represent anomalous realizations in the dynamics of stochastic systems. Due to their potentially significant impact, investigating and understanding rare events is crucial across various disciplines, including physics, living systems and climatology. This lecture serves as an introduction to the core principles of large deviation theory (LDT), offering an easily comprehensible outline of its formal structure and its applications. We further motivate its use through the specific case of diffusive systems in the large-size and long-time limits where we explore current fluctuations. We also introduce and explain the use of a population dynamics algorithm via a discrete, Markovian picture, and show how it allows us to explicitly evaluate relevant related LDT functions (so-called rate functions).

## 6.1 Introduction: Why large deviations?

Rare events are by definition anomalous situations that fall outside typical behavior. These events can encompass anything from harmless incidents such as a fortunate winning streak in a game of dice to significantly influential occurrences, such as the abrupt formation of crystals due to nucleation in a supercooled liquid, a sliding layer of snow triggering an avalanche, or large-scale natural disasters causing irreversible changes to a landscape. These motivate one to find and employ some systematic and practical framework for estimating their probability of occurrence. In terms of probability distributions, the focus lies not only on determining the most probable state, which typically centers around the mean, but also on examining its behavior up to its tails. The question consists in assessing the likelihood of a stochastic process fluctuating from its usual state to a particular rare value.

The theory of large deviations (LDT) answers precisely this question and provides us with a framework that allows to systematically estimate the probability of observing a particular rare event relative to its most likely value. In particular, LDT can provide an asymptotic estimation of such probabilities in the limit of very large observation windows. For instance, consider a coin toss game where we count the number of heads after $n$ coin tosses. LDT then allows us to estimate the likelihood of achieving any (possibly large) number of heads, in the asymptotics of large $n$ coin tosses. Certainly, in this particular example, one can explicitly calculate probabilities for all values of $n$ (as we will see). The strength of LDT, however, lies in scenarios where such explicit computations are generally unfeasible but become viable in the limit large number of occurrences.

LDT has expanded its range of applications beyond its original mathematical framework, serving as a relevant tool to characterize rare events across various fields. For example, in statistical physics, LDT allows one to investigate dynamical phase transitions (DPTs), which occur between phases that are characterized by different values of dynamical observables (such as flows or activities). Consider for instance a lattice one-dimensional system of particles with the ability to move either left or right (provided their target site is empty): values of particle current significantly lower than the average are characterized by a “jammed” phase with the formation of a cluster of particles. Similar DPTs have been described through the application of LDT to the distribution of entropy production in systems such as active Brownian particles, with a phase transition marked by a collective particle alignment [87]. Furthermore, LDT may aid in uncovering DPTs in other systems, such as chemical reaction networks or conformational changes in energy-driven biomolecules, by identifying rare dynamical regimes of their dynamics. At larger

<footnote n="20">²⁰VL was the lecturer; WP wrote the chapter.</footnote>
<figure><img src="image_56.png" /><figcaption>Figure 6.1: a. Temperature deviation in the northern hemisphere over a 90 days window. Large deviations in temperature were obtained from an LDT algorithm biased towards achieving a heat wave (rare event). Figure from Ref. [88]. b. A one dimensional lattice hopping model where particles tend to jump right (figure adapted Ref. [96]). Conditioning the system to present a low time-integrated particle currents, induces the emergence of density fronts (red arrows).</figcaption></figure>

scales, in the area of climatology and geophysical modeling, LDT has gained prominence in evaluating climate anomalies, including sustained, long-term heat waves across spatially extended regions (see Fig. 6.1 and Ref. [88]).

These lecture notes are intended as an extremely simple introduction to the notion of large-deviation scaling. The interested reader may continue with textbook references on the subject, for instance Refs. [89–95].

## 6.2 Basics and formalism

Having introduced some of the motivation behind LDT we now present its formal settings. Consider a stochastic process and a time-additive observable A defined on some observation window duration t as:

$$A = \int_{0}^{t} dt' \alpha(t') \qquad (6.1)$$

where α(t') depends on an underlying stochastic dynamics at time t'. Examples of such quantities can be obtained from any noisy phenomena such as a temperature series in a climate model over some period of observation, and arise naturally in physical systems when calculating quantities like work (integrated power) and other forms of integrated currents (e. g. the entropy production in a non-equilibrium system).

Of course, given the stochastic nature of the system, values of A are generally not fixed but fluctuate around a steady value that represents its more likely outcome. One may therefore want to characterize the probability distribution P(A, t) of A at time t that describes the rare events of interest. The LDT then tells that, in the limit of large t, P(A, t) scales as

$$P(A = at, t) \asymp e^{-tI(a)}, \qquad (6.2)$$

where ≍ is shorthand to denote a the logarithmic equivalent. In other words:

$$\lim_{t \to \infty} \frac{1}{t} \log P(A = at, t) = -I(a), \qquad (6.3)$$
and $I(a)$ is a so-called rate function. If the rate function $I(a)$ reaches its minimum in $\bar{a}$, then $\bar{a}$ is the typical value of $A/t$. As its name implies, it quantifies the rate at which the probability of observing a value $A$ deviates (exponentially rarely) from that of its typical value $\bar{a}t$. Fig. 6.2a provides a generic example for the behavior of such rate functions.

To better illustrate these concepts and help highlight the origin of the scaling form of Eq. (6.2) for the rate function we now turn to two basic examples.

### 6.2.1 A coin tossing game

Consider a fair coin toss game where a player wins a prize based on the number $F$ of heads obtained after $N$ trials. Thus $F$ in this case is the additive observable of interest and $N$ is the duration of our window of observation. This means that we are interested in the quantity

$$F = \sum_{i=1}^{N} \delta_{H,x_i} \qquad (6.4)$$

where $\delta$ is the Kronecker delta and $x_i \in \{T, H\}$ indicates the value tails or heads respectively for the random value $x_i$ of the process at instance $i$.

Since the coin is fair, the average of $F$ is $N/2$; its second cumulant is also proportional to $N$. This implies that small deviations of $F$ from its average scale as $\sqrt{N}$; in fact, the central limit theorem states that the distribution of $(F - N/2)/\sqrt{N}$ is Gaussian; in other words:

$$\lim_{N \to \infty} P(F = \frac{N}{2} + \sqrt{N} \delta \hat{F}) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{1}{2} \frac{\delta \hat{F}^2}{\sigma^2}} \qquad (6.5)$$

for some variance $\sigma > 0$ independent of $N$.

The large-deviation approach focuses on deviations of $F$ that are far more rare than those normal deviations. Namely, we are interested in the probability $P(F/N = f)$ that the time average $F/N$ achieves some value $f$. In our coin toss problem, this can be solved by representing the number of heads as per a binomial distribution. Explicitly

$$P(F/N = f) = \frac{N!}{2^N (fN)!(1-f)N!} \qquad (6.6)$$

In the large-$N$ limit, taking the log on both sides of Eq. (6.6) and apply Stirling's expansion $\ln N! = N \ln N - N + O(\ln N)$, one finds

$$\ln P(F/N = f) \approx -N \ln 2 - N f \ln f - N(1-f) \ln(1-f) \qquad (6.7)$$

$$P(F/N = f) \approx e^{-NI(f)},$$

where we identified the rate function of this process as $I(f) \equiv \ln 2 + f \ln f + (1-f) \ln(1-f)$. This scaling is an example of the LDT scaling of Eq. (6.2). It is represented in Fig. 6.2b. As expected, the most likely fraction of heads is $\bar{f} = 1/2$, which is the point where $I(f)$ is zero and reaches its minimum. Importantly, the rate function $I(f)$ describes fluctuations of $F/N$ that go beyond the Gaussian fluctuations around its average value.

### 6.2.2 Gaussian Sums

Consider now a sum $A$ of $N$ independent and identically distributed (i.i.d.) Gaussian random variables $\alpha$ with individual mean values $\mu$ and standard deviations $\sigma$:

$$A = \sum_{i=1}^{N} \alpha_i, \qquad (6.8)$$
<figure><img src="image_57.png" /><figcaption>Figure 6.2: a. A schematic rate function $I(a)$ illustrating the typical event (blue dashed line) and a rare event of an observable $a$ (orange dashed line). The magnitude of $I(a)$ (circle markers) dictates the likelihood, via $P \approx -tI(a)$, of observing a rare event relative to the typical. Hence larger magnitudes indicate lower probability. b. Rate function in Eq. (6.8) for the binomial example. The most typical value falls on $f = 0.5$ where $I(f) = 0$ c. Rate function for the Gaussian process example which yields a simple parabola.</figcaption></figure>

where $\alpha_i$ represents the individual realization of the random variable. As before, $N$ is the duration of the window of observation and $A$ the additive observable. We are interested in finding values of $a = A/N$ which represent possibly atypical values of $A$ away from its most likely value $N\mu$. Using that the sum of $N$ Gaussian i.i.d. variables yields another Gaussian variable with effective mean $\mu' = N\mu$ and square standard deviation $\sigma'^2 = N\sigma^2$, one has

$$P(A/N = a) = \sqrt{\frac{N}{2\pi\sigma^2}} e^{-\frac{N(a-\mu)}{2\sigma^2}}. \qquad (6.9)$$

$$\frac{1}{N} \ln P(A/N = a) = -\frac{(a-\mu)^2}{2\sigma^2} - \frac{1}{2N} \ln(2\pi\sigma^2) + \frac{1}{2N} \ln N,$$

and taking the large-$N$ limit one finds

$$\lim_{N \to \infty} \frac{1}{N} \ln P(A/N = a) = -\frac{(a-\mu)^2}{2\sigma^2} \qquad (6.10) \\ = I(a).$$

Note that terms sublinear on $N$ drop out and we end up with the rate function $I(a)$ which in this case is just quadratic. In other words we have shown that $P(A/N = a)$ presents again the LDT scaling $\asymp e^{-NI(a)} = e^{-N\frac{(a-\mu)^2}{2\sigma^2}}$. The most likely value $\bar{a} = \mu$ is the location of the minimum of the rate function.

In this specific example, the central-limit theorem and the LDT yield the same Gaussian expression for the large-$N$ behavior of the distribution of $A$, but this is in general not the case, as we have seen in the coin toss example above.

## 6.3 Application: Large Deviations in Diffusive Systems

We now move on to examples of physical relevance to model stochastic phenomena, and consider as a paradigmatic application the case of diffusive dynamics. Chiefly, we are
concerned with a physical system – discrete, particle based, or continuum – coupled to a thermal bath, and/or whose internal micro-dynamics gives rise to an inherent noise that drives fluctuations in the observables of the system. For instance, in the classic example of a 1D overdamped Brownian particle coupled to a thermal bath, one may model the dynamics of its position X as

$$ \dot{X} = \mu F + \sqrt{2D}\xi, \qquad (6.11) $$

where F represents a deterministic forcing in the system (possibly depending on X), and ξ is a Gaussian white noise with zero mean and correlation $\langle\xi(t)\xi(t')\rangle = \delta(t' - t)$. This noise represents for instance fluctuations coming from a thermal bath. The probability density ρ(x, t) of its configurational density is governed by the Fokker-Planck equation

$$ \begin{aligned} \partial_t \rho &= -\partial_x (\mu F \rho - D \partial_x \rho) \\ &= -\partial_x j, \end{aligned} \qquad (6.12) $$

where we identified the probability current j ≡ μFρ − D∂$_{x}$ρ. One is then typically interested in quantifying the fluctuations of time-additive observables A such as the work or the total entropy production. They depend on the level of non-equilibrium drive in the system, measured by j. This implies that knowledge of the fluctuations of j informs us on the fluctuations of such time-additive observables, and therefore gives access to the associated distributions of P(A, t). Identifying the distribution of j, via the dynamics of ρ, is therefore of great interest, but also necessitates the solution of the associated Fokker-Planck equation, which may be typically unfeasible in large and interacting systems.

We can alternatively attempt to coarse-grain over small details and consider a macroscopic picture in a steady-state where dynamics is determined by a very large number of particles evolving in time over diffusive time scales. This procedure allows one to transform a discrete (or particle) representation into a continuum limit, for large systems. Their dynamics may be externally driven or coupled to thermodynamic gradients and can describe in or out of equilibrium regimes. The usefulness of this approach, called Macroscopic Fluctuation Theory (MFT, see Ref. [97] for a review), is then that it provides us with a more general description of the non-equilibrium dynamics through the joint fluctuations of the densities and currents in the system. These latter may in part be captured by a joint-probability distribution function P(ρ, j) and therefore neatly encapsulate all relevant information about rare events.

Concretely, consider as an example a one-dimensional system coupled to two distinct thermal reservoirs and without external driving. This may represent some substance, (material, particles) confined in one dimension but subject to a thermal gradient which generates a density current. As before, we are interested in characterizing the system activity and its fluctuations through a current observable Q integrated over a given time interval. In particular we are interested in estimating P(Q, t) which, in the large-size limit N ≫ 1 and long-time limit t → ∞ LDT allows one to estimate as P(Q, t) ≍ e − tNI ( Q ), where I ( Q ) is the rate function associated with fluctuations in Q. As we will see, it is possible to estimate I ( Q ) from knowledge of the fluctuations in ρ and j which represent the underlying dynamics.

We first define the time-integrated current as

$$ Q = \int_{0}^{t} dt' \int_{0}^{1} dx \, j(x, t'), \qquad (6.13) $$

where j represents the change of density in the system as given by the following dynamics

$$ \partial_t \rho = -\partial_x j \qquad (6.14) $$

$$ j = -D \partial_x \rho + \sqrt{2\sigma\eta}. \qquad (6.15) $$
Here $\eta$ is a white noise term i.e. $\langle \eta(x,t)\eta(x',t') \rangle = \frac{1}{N}\delta(x-x')\delta(t-t')$ and $\sigma$ is its associated diffusion coefficient. Following the MFT, this noise term is considered to capture the underlying microscopic dynamics in this large time and particle limit. Furthermore, the Gaussian form of the noise allows one to immediately estimate the probability of observing a given trajectory of noise $\eta$ as

$$P(\eta) = e^{-\frac{N}{2} \int_0^t dt' \int_0^1 dx \, \eta(x,t')} . \qquad (6.16)$$

Then, rearranging from Eq. (6.15) we see that we can rewrite $\eta = \frac{j+D\partial_x\rho}{\sqrt{2\sigma}}$ and thus obtain a relation for the joint probability distribution of $\rho$ and $j$ as

$$P(\rho, j) = e^{-\frac{N}{2} \int_0^t dt' \int_0^1 dx \frac{(j+D\partial_x\rho)^2}{\sigma}} \qquad (6.17) \\ = e^{-tN I(\rho, j)},$$

where we identify $I \equiv \frac{1}{t} \int_0^t dt' \int_0^1 dx \frac{(j+D\partial_x\rho)^2}{\sigma}$ as a joint rate function of $\rho$ and $j$ given the large limit assumptions in $N$ and $t$, and where $j$ and $\rho$ are constrained to verify the continuity equation $\partial_t\rho + \partial_xj = 0$.

Next, we wish to estimate the likelihood of observing a given fluctuation of $Q$ through knowledge of $I(Q)$. Here, one expects that knowing the joint rate function of $\rho$ and $j$, and through the dependence of $Q$ on these variables, that we will be able to estimate fluctuations on $Q$ as well. Indeed, LDT allows one to map, or transform from a rate function $I(a)$ to another $I(b)$ in a mathematical procedure known as a contraction if there exists a function, or constraint that relates $a$ to $b$. Briefly, if $b = f(a)$ where $f$ is some function, then the contraction is defined as

$$I(b) = \min_{a: b=f(a)} I(a), \qquad (6.18)$$

where by an abuse of notation we denote by the same symbol the rate functions of the random variables $a$ and $b$. This result follows from a so-called saddle point analysis which states that for integrals of the form $\int dx e^{-rI(x)} \asymp e^{-r \min I(x)}$ as $r \to \infty$. Thus, for an exponentially weighted probability distribution, one seeks to pick the most likely event fulfilling some constraint which amounts to finding values that minimize the argument. Its namesake follows in the case where the linking function is many-to-one so that it leads to a net reduction, or contraction in the number of variables. More rigorous derivations of these results go beyond the scope of this lecture but are available from reviews in the literature.

Using this procedure we then finally obtain our desired estimates of $P(Q, t)$ via $I(Q)$. Namely (denoting for short $\int dt' \int dx = \int_0^t dt' \int_0^1 dx$),

$$I(Q) = \min_{j: Q = \int dt' \int dx \, j} \min_{\rho: \partial_t \rho + \partial_x j = 0} I(\rho, j) \qquad (6.19) \\ = \min_{j: Q = \int dt' \int dx \, j} \min_{\rho: \partial_t \rho + \partial_x j = 0} \int_0^t dt' \int_0^1 dx \frac{(j + D \partial_x \rho)^2}{2\sigma},$$

where the double minimization arises due to the bivariate of $I(\rho, j)$ with respect to the densities and currents, and the additional expressions underneath represent constraints as per Eq. (6.13) and Eq. (6.14). Thus, fluctuations in $Q$ are driven by the most likely fluctuations in $\rho$ and $j$ that coincide in achieving the given value of $Q$.

We can therefore appreciate how the use of LDT allows us to forego sampling through all possible realizations of a dynamics in a system, and instead provides a direct means to characterize observable fluctuations of a quantity of interest via knowledge of its corresponding rate functions. Our original goal of estimating the likelihood of rare events is then reduced to finding rate functions which we now address next.
## 6.4 Biased dynamics: making rare dynamics typical

Rare events are, by definition, rare and hence not necessarily observable within some feasible window of time in numerical approaches. Instead we may condition, or choose relevant trajectories, such that a particular value of an observable $A/t = a$ is obtained upon averaging. This is equivalent to working in a conditioned probability distribution of the form

$$P(x; A/t = a) = \int \mathcal{D}(X) p(X = x) \delta\left(\frac{1}{t}A(X) - a\right), \qquad (6.20)$$

where $\mathcal{D}(X)$ represents the path integral over all realizations of trajectories $X$. Such conditioned probabilities might seem familiar as they are the dynamical analogy to a microcanonical ensemble in equilibrium statistical systems, where an arbitrary $A$, rather than energy, now plays the constraining variable. Naturally, just as in the microcanonical ensemble, such conditioned ensemble can also be used to obtain other observables $\mathcal{O}$ of interest under the condition that $A/t = a$. Thus one may obtain averages of $\mathcal{O}$ as

$$\begin{aligned} \bar{\mathcal{O}}(a) &= \frac{\int \mathcal{D}(x) P(x; A/t = a) \mathcal{O}(x)}{\int \mathcal{D}(x) P(x; A/t = a)} \qquad (6.21) \\ &= \frac{\langle \mathcal{O} \delta\left(\frac{1}{t}A(X) - a\right) \rangle}{\langle \delta\left(\frac{1}{t}A(X) - a\right) \rangle} \end{aligned}$$

where $\mathcal{D}(x)$ denotes a path integral over all conditioned paths. Here, $\bar{\mathcal{O}}(a)$ represents the average of observable $\mathcal{O}$ over trajectories conditioned to $A/t = a$.

However, much like the microcanonical ensemble, rather than working with ensembles conditioned on $A/t = a$, it is often easier and more practical to instead have a modified ensemble where trajectory probabilities are biased by a weight $e^{\lambda A}$ where $\lambda$ is a fixed parameter. In fact, if $\mathcal{O}$ does not behave exponentially in time, it is possible to show that such procedure respects expectation values so that this ‘microcanonical’ average $\bar{\mathcal{O}}$ may now be re-expressed as

$$\bar{\mathcal{O}}(a) = \frac{\langle \mathcal{O}(X) e^{\lambda A(X)} \rangle}{\langle e^{\lambda A(X)} \rangle}, \qquad (6.22)$$

for a well-chosen $\lambda$, conjugated to $a$. Thus we can reinterpret this biased ensemble as a ‘canonical’ dynamical ensemble where the sum in the denominator may be interpreted as analogous to a dynamical partition function, and $\lambda$ plays the role of a conjugate variable to $A$ akin to the conjugate relation of temperature and energy in thermodynamics. In fact, the dynamical partition function $\langle e^{\lambda A(X)} \rangle$ implicitly encodes for all relevant statistical measurements of the system by recognizing: it is also the moment generating function (MGF) of the original probability density $P(X)$. That is, by noting that

$$\langle e^{\lambda A} \rangle = \int \mathcal{D}(X) P(X) e^{\lambda A} \qquad (6.23)$$

it follows that,

$$\frac{\partial^m}{\partial \lambda^m} \int \mathcal{D}(X) P(X) e^{\lambda A} \Big|_{\lambda=0} = \int \mathcal{D}(X) A^m P(X) e^{\lambda A} \Big|_{\lambda=0} = \langle A^m \rangle, \qquad (6.24)$$

so that the average, variances and higher order statistical measures of the biased ensemble are indeed readily accessible from the derivatives of the MGF in $\lambda = 0$. Furthermore, if we now consider the log of Eq. (10.34), known as the cumulant generating function (CGF), we
may then obtain mean-centered moments $\langle(A - \langle A \rangle)^2\rangle$ of the distribution, which therefore provides more relevant information of the deviations of a quantity from its average value.

In fact, it turns out that knowledge of the CGF allows us to compute the rate function $I(a)$ directly in the long-time limit as follows. In particular, it can be shown that

$$ \lim_{t \to \infty} \frac{1}{t} \ln \langle e^{\lambda A} \rangle = \psi(\lambda) \qquad (6.25) $$

where $\psi(\lambda, A)$ is known as the scaled cumulant generating function (SCGF). Then, following a result in LDT known as the Gärtner-Ellis theorem, we can now obtain the rate function as

$$ I(a) = \max_{\lambda} \lambda a - \psi(\lambda) . \qquad (6.26) $$

Notice that this extremalization principle gives precisely the value $\lambda$ conjugated to $a$ that was needed in Eq. (6.22). We can understand this result briefly by supposing that $P(A/t = a) \asymp e^{-tI(a)}$ and evaluating the MGF as

$$ \begin{aligned} \langle e^{\lambda A} \rangle &= \int da P(A = at) e^{\lambda at} \\ &= \int da e^{t(\lambda a - I(a))} \\ &\asymp e^{\max_{\lambda} t(\lambda a - I(a))}, \end{aligned} \qquad (6.27) $$

where the last line follows from the saddle-point evaluation for $t \to \infty$. The result in Eq. (6.25) then follows from taking the log on both sides of Eq. (6.28) and letting $t \to \infty$. Thus the likelihood of fluctuations of $A$ can be obtained from the limiting calculation of $\psi(\lambda)$ alone.

In practice, different analytical and numerical methods are available in obtaining $\psi$ for actual systems. However, the general challenge lies in obtaining accurate averages of the observables for increasing values of the bias $\lambda$. In what follows we will motivate a numerical approach to calculate $\psi$ in Markovian systems as inspired in a population dynamics algorithm.

### 6.4.1 Population dynamics algorithm: a Markov jump process motivation

Consider a Markov process defined by a finite number of configurations $\{x\}$ and whose dynamics are determined by transition rates $r(x|x')$ to go from state $x'$ to $x$. For example, in a discrete graph system these configurations could correspond to the number of nodes and the transition rate the probability of jumping from one node to another. The dynamics of the probability distribution of states $p(x, t)$ is then given by the master equation

$$ \partial_t p(x, t) = \sum_{x' \neq x} \left[ r(x|x')p(x', t) - r(x'|x)p(x, t) \right], \qquad (6.28) $$

and which we can interpret as the difference between a gain term, from states $x'$ going to $x$, and a loss term from states leaving $x$.

Now consider an additive observable $A$ which can be written over a time history of $K$ jumps as

$$ A = \sum_{0 \le k \le K-1} \alpha_{x_k x_{k+1}}, \qquad (6.29) $$
where $\alpha_{x_k x_{k+1}}$ are individual contributions along the history path of successively visited states $\{x_k\}_{0 \le k \le K}$. The quantity $\alpha_{x_{x'}}$ represents the amount by which $A$ is increased during a jump from state $x$ to state $x'$. For instance, in a 1D particle system, taking $\alpha = \pm 1$ depending on whether a particle jumps towards right/left will yield for $A$ a time-integrated current. Our goal is to compute the likelihood of observing fluctuations of this observable through knowledge of its corresponding SCGF $\psi(\lambda)$.

We may achieve this by first considering a joint probability distribution $p(x, A, t)$ of the observable and the states as

$$ \partial_t p(x, A, t) = \sum_{x' \neq x} \left[ r(x|x')p(x', A - \alpha_{x'x}, t) - r(x'|x)p(x, A, t) \right]. \quad (6.30) $$

where the $-\alpha_{x'x}$ term accounts for the fact that we are coming from the previous value of $A$ at the $x'$ state. Next we are interested in biasing dynamics with respect to specific values $A$ which, as discussed in the previous section, can be achieved by exponentially weighing the corresponding probability distribution $p(x, A, t)$. This then corresponds to the biased ‘canonical’ ensemble expressed as

$$ \hat{p}(x, \lambda, t) = \sum_{A} e^{\lambda A} p(x, A, t), \quad (6.31) $$

where $\lambda$ is the conjugate variable with $A$ and determines the ‘strength’ of the biasing. We can then multiply by the corresponding weights $e^{\lambda A}$ in the Eq. (6.30), sum over all $A$ and find the desired biased, joint dynamics

$$ \partial_t \hat{p}(x, \lambda, t) = \sum_{x' \neq x} \left[ e^{-\lambda \alpha_{x',x}} r(x|x') \hat{p}(x', \lambda, t) - r(x'|x) \hat{p}(x, \lambda, t) \right]. \quad (6.32) $$

Finally using the result in Eq. (6.25) we recognize that

$$ \sum_x \hat{p}(x, \lambda, t) = \langle e^{\lambda A} \rangle \underset{t \to \infty}{\asymp} e^{t\psi(\lambda)}, \quad (6.33) $$

so that our task involves evaluating the biased dynamics of Eq. (6.32) and obtaining the expectation values with respect to $p(x, \lambda, t)$ and hence $\psi(\lambda)$. Note further that the exponential modification of the rates means the probability is no longer conserved in Eq. (6.32) and instead we are now dealing with an open system.

Analytically, such calculation is equivalent to solving a maximum eigenvalue problem of the form $R_\lambda |\hat{p}\rangle = e^{t\psi(\lambda)}|\hat{p}\rangle$, where $|\hat{p}\rangle$ represents the eigenvector and $R_\lambda$ an evolution operator with matrix entries $R_{x,y} = e^{\lambda \alpha_{yx}} r(x|y) - \delta_{x,y} r(x)$ (with $r(x) = \sum_y r(y|x)$ is the escape rate). However, for more general cases where closed solutions are not tractable, we can instead consider a numerical approach that, as we will now demonstrate, amounts to a population dynamics algorithm with respect to $\lambda$ and $x$.

First let us then re-write Eq. (6.32) as

$$ \partial_t \hat{p}(x, \lambda, t) = \sum_{x' \neq x} r_{\lambda}(x|x') \hat{p}(x', \lambda, t) - r_{\lambda}(x) \hat{p}(x, \lambda, t) + \delta r_{\lambda}(x) \hat{p}(x, \lambda, t), \quad (6.34) $$

where

$$ r_{\lambda}(x'|x) = e^{\lambda \alpha_{x'|x}} r(x'|x) \quad (6.35) $$

$$ r_{\lambda}(x) = \sum_{x' \neq x} r_{\lambda}(x'|x) \quad (6.36) $$

$$ \delta r_{\lambda}(x) = r_{\lambda}(x) - r(x), \quad (6.37) $$
<figure><img src="image_58.png" /><figcaption>Figure 6.3: Schematic of a population dynamics evolution with a biased selection rule for a system with $N_c = 4$. Clones displaying larger biases towards a desired event are multiplied, while those with average or below performance are deleted. The reproduction success i.e. population deletion/multiplication is controlled by the bias parameter $\lambda$. Surviving trajectories (highlighted here in red) are those representative of the biased ensemble.</figcaption></figure>

and are, respectively, the biased transition rates from $x'$ to $x$, the net biased escape rate from $x$, and the difference between the original and modified escape rates in the biased ensemble. Importantly, we can recognize the first terms in Eq. (6.34) as corresponding to a *conserved* dynamics of rates $r_\lambda$, plus an additional term $\delta r_\lambda$ representing *gain* or *loss* in the probability depending on sign.

Furthermore, the Feynman-Kac formula tells that Eq. (6.34) implies

$$ \langle e^{\lambda A} \rangle = \langle e^{\int_0^t dt' \delta r_\lambda(x(t'))} \rangle_\lambda \quad (6.38) $$

where the subscript in $\lambda$ indicates the biased ensemble (i.e. of transition rates $r_\lambda$), so that calculation of $\psi(t)$, from its long time limit relation in Eq. (6.33), may be now be reinterpreted as a history average over $\delta r_\lambda$. Given that this term is precisely the injection/loss of probability in the dynamics, we can interpret this procedure as a population dynamics algorithm where population members, all evolving in parallel in the biased dynamics, play the role of probability outcomes to be removed (probability loss) or replicated (probability gain) as they evolve in time.

The idea to represent probability loss or gain comes from similar quantum-mechanical problems (finding the ground-state energy of an operator using diffusion Monte-Carlo methods [98]) and was applied to varied problems in statistical mechanics [99] and mathematics [100]. Its application to compute SCGFs was put forward in [96, 101].

Concretely, we can now formulate a population dynamics algorithm as follows. First, define $Y \equiv e^{\int_0^t dt' \delta r_\lambda(x(t'))}$ and discretize the integral in Eq. (6.38) as $Y = \prod_{t_k=0}^t Y(x_k)^{\Delta t'}$ where $Y(x_k)^{\Delta t'} \equiv e^{\delta t'(r_\lambda(x_k) - r(x_k))}$, and $x_k$ represent configurations at a given time $t_k$. Then, starting from a very large population of $N_0$ clones all evolving in parallel under the biased dynamics:

* Compute $Y(x_k)^{\Delta t}$ during an interval $\Delta t$

* Clone $n_i(x_k)$ is pruned or replicated with rate $Y(x_k)^{\Delta t}$. For instance, given random
a variable $\epsilon$ uniformly distributed on $[0, 1]$, make $y = \lfloor Y(x_k) \Delta t' + \epsilon \rfloor$ copies of the clone. This then alters the population by a factor $F_t = \frac{N+y-1}{N}$

* Evolve new modified population for the next interval $\Delta t$ and repeat process up to $t$.

A schematic of this population dynamics with selection is illustrated in Fig. 6.3. The SCGF is then evaluated from the change in population as

$$ \psi(\lambda, t) \approx \frac{1}{t} \ln\langle Y \rangle_{\lambda} \qquad (6.39) \\ = \frac{1}{t} \ln \prod_{t} F_{t}. $$

Therefore, $\psi(t)$ can ultimately be interpreted as the exponential rate of growth or decrease in a population. However, in practice such cloning procedure may lead to diverging or vanishing populations which could be undesirable from a numerical point of view. Instead, one may also choose to maintain a constant population throughout the evolution of the system by different selection rules that prune/copy clones from a statistically weighted chance with respect to the entire population. Specific examples of these rules, as well as other choice of cloning regimes, are for instance described in the review [102].

Altogether, through this cloning algorithm we are therefore able to sample for biased values of an observable $A$ for any Markovian system of choice. Such procedure is formally exact in the limit of large number of clones, and simulation times, and has been used successfully to probe for phase transitions in both discrete and particle based systems. However, achieving large parameter limits can be challenging in practice, and finite-size and finite-time scaling analysis can be helpful [103, 104].
# 7 Dynamical Large Deviations: at large size

**Alexandre Lazarescu.** In the second part of the lecture, we focus on large deviations at large size and finite time. Using examples of population models, and particularly chemical networks, we will see how far we can push the large deviations principle to obtain information on not only time averages, but also precise time evolutions of large systems, using path integrals.

## 7.1 Dynamical large deviations

In the first part of this lecture, we saw how the large deviations principle can be applied to describe fluctuations of time-averaged observables in a Markovian process, when time becomes large. In this second part, we focus on a different scaling: we do not put any restrictions on time, but we assume the system to be very large. This large-size limit is the same as the one we take in equilibrium statistical mechanics in order to reach the thermodynamic limit; here, we merely consider the distributions to be time-dependent.

If we focus on interacting particle models, there are essentially two ways to reach that limit: either look at particles on a fixed network, and increase the number of particles until there is a macroscopic number of particles on each site, or consider particles on a lattice or in continuous space, and rescale space until the density of particles at every point is finite. The first class of systems is generally called *population models*, and includes chemical networks and certain ecological models, and the second is often referred to as *diffusion models* (see Section 6.3), even if the resulting dynamics is sometimes ballistic. In this lecture, we will focus on populations, but we will provide a few formulas relating to diffusions when appropriate. In short, we will explain where formula (6.17) comes from, and what we can do with it.

Our approach will be the following: we first set up a few explicit population models where we can take a large size limit; we then present the fundamental building bloc of dynamical large deviations: the instantaneous displacement distribution, and its generating function; we use said building blocks to construct path integrals; from those path integrals, we derive equations of motion, both in a Lagrangian and a Hamiltonian formalism; we make a brief aside to present the similarities and differences between those formalisms within statistical mechanics, quantum mechanics, and classical mechanics; and, finally, we describe the long time behaviour of those path integrals.

### 7.1.1 Formalism

The majority of objects and relations that we will present in this lecture are quite general: they will involve generic Markov processes, with some static observables (densities) and some dynamical observables (currents), and a relation between the two (continuity equation). This abstraction makes our framework quite powerful, but also rather formal. It is then useful to first consider a few concrete cases, where explicit calculations can be made, and pictures can be drawn.

#### Population models

As is customary, we start with independent particles. Consider $N$ independent random walkers governed by a process $R$ (as defined in Section 1.2). They live on a fixed finite network with sites $x$, and at every time $t$ we denote the occupancies of every site as $n_x(t)$, and the corresponding densities as $\rho_x \equiv n_x/N$.
The jump process R is Markovian and performs transitions between sites ($x \to y$) with rates $r_{yx}$, which can be translated into transitions between occupancies, where occupancy vectors transform as {$n_z$} $\to$ {$n_z - \delta_{z,x} + \delta_{z,y}$} with rate

$$R_{yx}^{(N)} = n_x r_{yx} = N \rho_x r_{yx} \equiv N w_{yx}(\rho). \qquad (7.1)$$

We can count the number of jumps between each pair of neighbours at any time (which we will denote with a generic # symbol for "number"), over a short period, i.e. between times $t$ and $t + \delta t$. This yields a current observable

$$\lambda_{yx} = \frac{\#[x \to y]_t^{t+\delta t}}{N \delta t} \qquad (7.2)$$

which is a random variable taking a time-dependent value for each realisation. Note that the typical value of this random variable is related to the semi-currents defined in Section 1.4, in the case of independent particles (there is simply a factor $N$ between them).

This current observable completely determines the evolution of the density observable: the change in density at any site is the difference between all incoming and outgoing currents. This can be formally written using a graph divergence $\nabla$ (which is the opposite of the incidence matrix $D$ defined in Section 2.2) such that we have the usual continuity equation $\dot{\rho} = -\nabla \cdot \lambda$. This matrix can be written as

$$\nabla_{z,(y,x)} = \delta_{z,x} - \delta_{z,y} \qquad (7.3)$$

and acts on edges $(y, x)$ to produce a vector on sites $z$.

We can make things slightly more complex by changing the type of transitions that are allowed in terms of occupancies, by removing or adding various numbers of particles from any number of sites. If we rename sites as species, what we get is a chemical network.

Consider thus a chemical system of volume $V$ with species $x$. We collect the number of particles of each species at time $t$ into a composition vector $n_x(t)$, and we define the corresponding concentrations as $\rho_x = n_x/V$. In order to transform groups of particles into other groups of particles (i.e. to perform chemical reactions), it is useful to first determine which groups can appear: we define complexes $\gamma$ with stoichiometry $\nu_x^\gamma$ (the number of particles of species $x$ which appear in complex $\gamma$). Note that those stoichiometric coefficients are sometimes indexed by the label and direction of a reaction instead (as in Section 4.2).

Given this, we can cast reactions as being transformations between complexes ($\gamma \to \gamma'$), with the corresponding change of occupancy {$n_z$} $\to$ {$n_z - \nu_z^\gamma + \nu_z^{\gamma'}$}. We call $\gamma$ the reactants, and $\gamma'$ the products, as is customary. Note that we treat forward and backward reactions as different transformations, because we don't assume micro-reversibility. The simplest transition rates we can then consider are called mass action rates, which contain a reaction constant $K$ which we scale as $k$ times a suitable power of $V$, and a combinatorial factor counting the number of ways to choose the reactants:

$$R_{\gamma'\gamma}^{(N)} = \left( \prod_x \frac{n_x!}{[n_x - \nu_x^\gamma]^!} \right) K_{\gamma'\gamma} = V \rho^{\nu^\gamma} k_{\gamma'\gamma} \equiv V w_{\gamma'\gamma}(\rho). \qquad (7.4)$$

Note the similarity in structure with the previous simpler case (7.1), which was also made of a combinatorial factor $n_x$ multiplied by a constant term tied to the type of transition $r_{yx}$. This case can be recovered if all reactions involve only one reactant and one product, like for enzyme conformation changes, for instance.
As for independent particles, we can define a (chemical) current observable by counting the number of times a transformation is made between two instants:

$$ \lambda_{\gamma'\gamma} \equiv \frac{\#[\gamma \to \gamma'_{t}^{t+\delta t}]}{V \delta t}. \qquad (7.5) $$

And just as before, this current observable completely determines the evolution of the concentration of each species in the system: the change in concentration is the difference between the number of particles created and the number of particles destroyed, where each number can be determined by the number of occurrences of each reaction, multiplied by the stoichiometric coefficient of the species in that reaction. This can be written in terms of a continuity equation $\dot{\rho} = -\nabla \cdot \lambda$ which involves a chemical divergence $\nabla$ which is essentially twice the stoichiometric matrix 4.4:

$$ \nabla_{z,(\gamma',\gamma)} = \nu_z^{\gamma} - \nu_z^{\gamma'}. \qquad (7.6) $$

This structure can be generalised, and will be the basis of our abstract analysis: we will always require our systems to have a well-defined volume $V$, states $\rho$, currents $\lambda$, state-dependent rates $R(\rho) \sim Vw(\rho)$, and a divergence operator $\nabla$ such that the continuity equation $\dot{\rho} = -\nabla \cdot \lambda$ translates currents into changes of state. The appropriate state and current observables are often obvious from the definition of the model. If the currents do not completely determine $\dot{\rho}$, it often means that some parts of the process were forgotten (either from missing transitions or an inappropriate bunching of transitions). A suitable set of currents will be called *fundamental currents*, and we should note that it is not unique: a transition can always be duplicated and its rate split in halves, without changing the process.

Note that, due to common usage, we will use the same symbol $\nabla$ for gradients, acting on state vectors and producing edge vectors. Also due to usage, the gradient will be defined as minus the transpose of the divergence, i.e. for instance

$$ \nabla_{(y,x),z} = -\delta_{z,x} + \delta_{z,y}. \qquad (7.7) $$

The distinction between the two symbols will be made through the fact that products in edge-space are denoted by a dot, as in $\nabla \cdot \lambda$, whereas products in state-space are not, such as $\nabla \rho$.

### Example

Let us conclude this part with a concrete example: the replicator, involving a single species A with occupancy $n$ and three reactions

$$ A \to 2A \quad \text{with rate } an \qquad \text{such that } w_{21}(\rho) = a\rho \qquad (7.8) $$

$$ 2A \to A \quad \text{with rate } bn(n-1)/V \qquad \text{such that } w_{12}(\rho) = b\rho^2 \qquad (7.9) $$

$$ A \to \emptyset \quad \text{with rate } cn \qquad \text{such that } w_{01}(\rho) = c\rho. \qquad (7.10) $$

We can define three complexes $\gamma_1 = A$, $\gamma_2 = 2A$, $\gamma_0 = \emptyset$, with stoichiometries $\nu^{(1)} = 1$, $\nu^{(2)} = 2$, $\nu^{(0)} = 0$.

In this case, the divergence is a one-line matrix $\nabla = [1-2, 2-1, 1-0] = [-1, 1, 1]$ in basis $\{A \to 2A, 2A \to A, A \to \emptyset\}$.

This example will reappear regularly as illustration in the rest of the lecture.
**Displacement distribution at large volume** Having defined the necessary mathematical objects, we can now look at the dynamics of our process. We know the rates of every possible transition, each associated to a current observable. Looking at the joint statistics of all currents $\lambda(t)$ starting from a given state $\rho(t)$ is in general complicated: the probabilities of that collection of jumps depends on the order in which they are performed (the individual processes of the transition matrix do not commute). However, if the size $V$ is large and all occupancies $n_x$ scale with $V$, a small displacement $\delta\rho$ will not affect the rates much: the processes become effectively independent. This is the most crucial approximation we make, and it relies on a specific scaling of the dynamics: typical occupancies must scale as $V$, and the number of occurrences of each transition during $\delta t$ starting from a typical state must scale as $V\delta t$, so that the continuity equation $\dot{\rho} = -\nabla \cdot \lambda$ is finite.

We are still considering population models with Markovian dynamics (hence with Poisson-distributed waiting times, called Poisson clocks, for each transition), and we want to describe the statistics of whole trajectories of the system. Since a trajectory is a sequence of small steps, a good starting point is to look at the distribution of those steps, i.e. the distribution of the rate of displacement of $\rho$, from any starting point, at any time, over a small time step $\delta t$. We have

$$P_{\delta t}(\dot{\rho}|\rho) = P[\rho(t + \delta t) = \rho + \delta t\dot{\rho} | \rho(t) = \rho] \quad (7.11)$$

with

$$P_{\delta t}(\dot{\rho}|\rho) = \sum_{\lambda; \dot{\rho}=-\nabla \cdot \lambda} P_{\delta t}(\lambda|\rho), \quad (7.12)$$

i.e. the probability of a displacement in $\rho$ is the sum of probabilities of all currents that produce that specific displacement. Those are, in principle, joint probabilities of currents, which are complicated. However, given the scaling discussed above, and given the set of possible transitions $(x \to y)$, we have in fact a quasi-independence

$$P_{\delta t}(\lambda|\rho) \approx \prod_{x,y} P_{\delta t}(\lambda_{yx}|\rho). \quad (7.13)$$

For each transition $x \to y$, observing a flux $\lambda_{yx}$ means that the Poisson clock of that transition, with rate $Vw_{yx}(\rho)$, must ring $\delta tV\lambda_{yx}$ times during a time $\delta t$, with probability

$$P_{\delta t}(\lambda_{yx}|\rho) = \frac{(\delta t V w_{yx}(\rho))^{\delta t V \lambda_{yx}}}{[\delta t V \lambda_{yx}]!} e^{-\delta t V w_{yx}(\rho)}. \quad (7.14)$$

In the limit $V \to \infty$, $\delta t \to 0$, $V\delta t \to \infty$, i.e. large volume, short time, large number of jumps, we can perform a Stirling approximation:

$$P_{\delta t}(\lambda_{yx}|\rho) \approx e^{-\delta t V L_{yx}(\lambda_{yx};\rho)} \quad \text{with} \quad L_{yx}(\lambda_{yx};\rho) = \lambda_{yx} \ln \left( \frac{\lambda_{yx}}{w_{yx}(\rho)} \right) - \lambda_{yx} + w_{yx}(\rho). \quad (7.15)$$

At this stage, we have entered large deviations territory: the function $L_{yx}$ is a rate function for one of the currents, conditioned on the state from which it is produced. From here on, we will combine LDFs to construct other LDFs and SCGFs, and mainly play around inside the exponentials. Most of the calculations and formulas will be made on those objects, which are all logarithms of probabilities, and translated back into probabilities only at the end.
Multiplying the probabilities for each transition, we get the joint LDF for all the currents, which we call detailed Lagrangian for reasons which will become clear later:

$$P_{\delta t}(\lambda|\rho) \approx e^{-N\delta tL(\lambda;\rho)} \quad \text{with} \quad L(\lambda;\rho) = \sum_{x,y} \lambda_{yx} \ln\left(\frac{\lambda_{yx}}{w_{yx}(\rho)}\right) - \lambda_{yx} + w_{yx}(\rho). \quad (7.16)$$

Since a displacement is a contraction of currents (i.e. a linear combination), we can use the contraction principle to obtain the LDF of $\dot{\rho}$, which we call the standard Lagrangian

$$P_{\delta t}(\dot{\rho}|\rho) \approx e^{-N\delta t\mathcal{L}(\dot{\rho};\rho)} \quad \text{with} \quad \mathcal{L}(\dot{\rho};\rho) = \inf_{\lambda: \dot{\rho}=-\nabla\cdot\lambda} [L(\lambda;\rho)]. \quad (7.17)$$

In most cases, we are not able to perform the minimisation explicitly, so that this object remains implicitly defined, unlike the detailed Lagrangian which is explicit.

### Example

For the replicator, we have three currents $\lambda_{21}$, $\lambda_{12}$, $\lambda_{01}$, related respectively to the transitions with rate constants a, b and c, and we can compute

$$L(\lambda; \rho) = \lambda_{21} \ln \left( \frac{\lambda_{21}}{a\rho} \right) + \lambda_{12} \ln \left( \frac{\lambda_{12}}{b\rho^2} \right) + \lambda_{01} \ln \left( \frac{\lambda_{01}}{c\rho} \right) - \lambda_{21} - \lambda_{12} - \lambda_{01} + a\rho + b\rho^2 + c\rho. \quad (7.18)$$

The way those currents contribute to the change in concentration, encoded in the divergence operator defined above, is $\dot{\rho} = \lambda_{21} - \lambda_{12} - \lambda_{01}$. In this simple one-dimensional case, the minimisation of L on the currents conditioned on a fixed value of $\dot{\rho}$ can be performed explicitly, by substituting $\lambda_{21}$ by $\dot{\rho} + \lambda_{12} + \lambda_{01}$ and looking for the values of $\lambda_{12}$ and $\lambda_{01}$ where L is minimal. We leave it to the overzealous reader to check that the result is

$$\mathcal{L}(\dot{\rho}; \rho) = \dot{\rho} \ln \left( \frac{\dot{\rho} + \sqrt{\dot{\rho}^2 + 4a\rho(b\rho^2 + c\rho)}}{2a\rho} \right) - \frac{\dot{\rho}^2 - (2a\rho - \sqrt{\dot{\rho}^2 + 4a\rho(b\rho^2 + c\rho)})^2}{4}. \quad (7.19)$$

The same result can be obtained much more easily using Hamiltonians rather than Lagrangians, which we explain next.

We can alternatively do computations at the level of SCGFs, via the generating function over all possible numbers of jumps $k = \delta tV\lambda_{yx}$ for a given transition $(x \to y)$, with a conjugate variable $f_{yx}$:

$$G_{yx}(f_{yx}) = \sum_{k} e^{kf_{yx}} P_{\delta t}(k/\delta tV|\rho) = \sum_{k} e^{kf_{yx}} \frac{(\delta tV w_{yx}(\rho))^k}{k!} e^{-\delta tV w_{yx}(\rho)} = e^{\delta tV w_{yx}(\rho)(e^{f_{yx}} - 1)}. \quad (7.20)$$

We can express it in terms of a SCGF:

$$G_{yx}(f_{yx}) \approx e^{\delta tV H_{yx}(f_{yx};\rho)} \quad \text{with} \quad H_{yx}(f_{yx};\rho) = w_{yx}(\rho)(e^{f_{yx}} - 1). \quad (7.21)$$

Multiplying all the independent generating functions, we get the joint SCGF, which we call detailed Hamiltonian:

$$G(f) \approx e^{\delta t yx H(f;\rho)} \quad \text{with} \quad H(f;\rho) = \sum_{x,y} w_{yx}(\rho)(e^{f_{yx}} - 1). \quad (7.22)$$
Contracting at the level of the SCGF means specifying the variable to a more restricted value, which is much easier to perform than a minimisation: for instance, it might involve specifying a vector of independent variables to a value depending on only one parameter, rather than finding the optimal value of L under a single linear constraint involving all variables. In the present case, contracting from currents $\lambda$ to displacements $\dot{\rho} = -\nabla \cdot \lambda$ translates to specifying the dynamical variable $f$ to the gradient of a state variable $u$, i.e. $f = \nabla u$. We obtain the SCGF of displacements, which we call the standard Hamiltonian:

$$ \tilde{G}(u) \asymp e^{\delta t V \mathcal{H}(u; \rho)} \quad \text{with} \quad \mathcal{H}(u; \rho) = H(\nabla u; \rho), \qquad (7.23) $$

where $u$ is the variable conjugate to the displacement rate $\dot{\rho}$, as can be determined by writing the consistent Legendre scalar product

$$ u\dot{\rho} = -u(\nabla \cdot \lambda) = (\nabla u) \cdot \lambda = f \cdot \lambda. \qquad (7.24) $$

### Example:

For the replicator, we have

$$ H(f; \rho) = a\rho(e^{f_{21}} - 1) + b\rho^2(e^{f_{12}} - 1) + c\rho(e^{f_{01}} - 1). \qquad (7.25) $$

Specifying the entries of $f$ as

$$ \begin{bmatrix} f_{21} \\ f_{12} \\ f_{01} \end{bmatrix} = \nabla u = \begin{bmatrix} u \\ -u \\ -u \end{bmatrix}, \qquad (7.26) $$

where we remember that the gradient is minus the transpose of the divergence. Injecting this into $H$ yields the standard Hamiltonian

$$ \mathcal{H}(u; \rho) = a\rho(e^u - 1) + b\rho^2(e^{-u} - 1) + c\rho(e^{-u} - 1) = (a\rho - (b\rho^2 + c\rho)e^{-u})(e^u - 1), \qquad (7.27) $$

whose Legendre transform can be taken relatively easily, to obtain the result found above in (7.19).

### A few common functions you may encounter

This section is a short list of Lagrangians and Hamiltonians of common classes of models, for reference. They are all written at the most detailed level, where everything is explicit.

First, independent walkers on a network. Sites $x$, densities $\rho_x$, transition rates $r_{yx}$, currents $\lambda_{yx}$. The detailed Lagrangian and Hamiltonian are:

$$ L(\lambda; \rho) = \sum_{x,y} \lambda_{yx} \ln \left( \frac{\lambda_{yx}}{r_{yx}\rho_x} \right) - \lambda_{yx} + r_{yx}\rho_x, \quad H(f; \rho) = \sum_{x,y} r_{yx}\rho_x(e^{f_{yx}} - 1). \qquad (7.28) $$

Note that, due to the independence of the particles, the Hamiltonian is linear in $\rho$. We sometimes refer to such cases as linear processes.
Second, a chemical network. Species $x$, complexes $\gamma$, concentrations $\rho_x$, reaction constants $k_{\gamma'\gamma}$, currents $\lambda_{\gamma'\gamma}$, stoichiometries $\nu_x^\gamma$. The detailed Lagrangian and Hamiltonian are:

$$L(\lambda; \rho) = \sum_{\gamma, \gamma'} \lambda_{\gamma'\gamma} \ln \left( \frac{\lambda_{\gamma'\gamma}}{\rho^{\nu^\gamma} k_{\gamma'\gamma}} \right) - \lambda_{\gamma'\gamma} + \rho^{\nu^\gamma} k_{\gamma'\gamma} \quad , \quad H(f; \rho) = \sum_{\gamma, \gamma'} \rho^{\nu^\gamma} k_{\gamma'\gamma} (e^{f_{\gamma'\gamma}} - 1) \qquad (7.29)$$

The Hamiltonian is no longer linear, unless the stoichiometries are such that the particles are independent. The structure is however very much the same as above, because in both cases the distribution of individual currents is Poissonian.

Third, independent diffusions. Particles in continuous space $\mathbb{R}^d$, with positions $x$, local density $\rho(x)$, drift field $V$, diffusion matrix $D$, local currents $j(x)$. The detailed Lagrangian and Hamiltonian are:

$$L(j; \rho) = \iint (j - V\rho + D\nabla\rho) \cdot \frac{D^{-1}}{4\rho} (j - V\rho + D\nabla\rho) \, dx dy, \\ H(f; \rho) = \iint f \cdot (Df\rho + V\rho - D\nabla\rho) \, dx dy \qquad (7.30)$$

Note that, once again, the Hamiltonian is linear in $\rho$.

Fourth, interacting diffusions, such as the famous WASEP. Same parameters as above, but with an extra function: the mobility $\sigma(\rho)$, which is no longer proportional to $\rho$. The detailed Lagrangian and Hamiltonian are:

$$L(j; \rho) = \iint (j - V\sigma(\rho) + D\nabla\rho) \cdot \frac{D^{-1}}{4\sigma(\rho)} (j - V\sigma(\rho) + D\nabla\rho) \, dx dy, \\ H(f; \rho) = \iint f \cdot (Df\sigma(\rho) + V\sigma(\rho) - D\nabla\rho) \, dx dy \qquad (7.31)$$

The Hamiltonian is no longer linear, unless the mobility $\sigma(\rho)$ is proportional to $\rho$, i.e. unless the particles are independent. The structure is however the same as above, with $L$ being quadratic in $j$ and $H$ being quadratic in $f$, because in both cases the distribution of individual currents is Gaussian.

Finally, the Gaussian approximation of a chemical network, related to the Chemical Langevin Equation (CLE). Same parameters as the second example, but here the currents are bidirectional and written as $j_{\gamma'\gamma}$. The detailed Lagrangian and Hamiltonian are:

$$L(j; \rho) \sim \sum_{\gamma, \gamma'} \frac{(j_{\gamma'\gamma} - (\rho^{\nu^\gamma} k_{\gamma'\gamma} - \rho^{\nu^\gamma'} k_{\gamma\gamma'}))^2}{2 (\rho^{\nu^\gamma} k_{\gamma'\gamma} + \rho^{\nu^\gamma'} k_{\gamma\gamma'})}, \\ H(f; \rho) \sim \sum_{\gamma, \gamma'} f_{\gamma'\gamma} \left( \frac{f_{\gamma'\gamma}}{2} (\rho^{\nu^\gamma} k_{\gamma'\gamma} + \rho^{\nu^\gamma'} k_{\gamma\gamma'}) + (\rho^{\nu^\gamma} k_{\gamma'\gamma} - \rho^{\nu^\gamma'} k_{\gamma\gamma'}) \right) \qquad (7.32)$$

Those functions are simply the second order approximation of the functions above, around typicality ($\lambda_{\gamma'\gamma} = \rho^{\nu^\gamma} k_{\gamma'\gamma}$ and $f = 0$). This translates the fact that this is a Gaussian approximation. This approximation is quite bad for any observable which cares about rare events, such as entropy production, or any marker of irreversibility.
Path integrals

Having described the little steps, we can combine them to get whole paths. For a finite time step $\delta t$, the probability of a path is a simple product

$$P_t[\{\rho_k\}] = \prod_{k=1}^{K} P_{\delta t}(\dot{\rho}_k | \rho_k). \qquad (7.33)$$

Keeping only the initial and final conditions fixed, we can write a finite transition probability in terms of these paths, and also express the path probabilities in terms of currents rather than displacements, by introducing the continuity equation as a constraint. We get

$$\begin{aligned} P_t(\rho_K | \rho_0) &= \sum_{\{\rho_k\}} P_t[\{\rho_k\}] = \sum_{\{\rho_k\}} \prod_{k=0}^{K-1} P_{\delta t}(\rho_{k+1} | \rho_k) \\ &= \sum_{\{\lambda_k\}} \prod_{k=0}^{K-1} P_{\delta t}(\lambda_k | \rho_k) \delta(\rho_{k+1} - \rho_k - \nabla \cdot \lambda \delta t), \end{aligned} \qquad (7.34)$$

A $\delta t \to 0$ limit will turn this expression into a path integral, where the dynamical variable is either the displacement rate $\dot{\rho}$, or the currents $\lambda$, with an extra continuity constraint. This constraint, as well as the initial and final condition, can be implemented by delta functions.

$$\begin{aligned} P_t(\rho_t | \rho_0) &\asymp \int e^{-V} \int_0^t \mathcal{L}(\dot{\rho}; \rho) d\tau \delta(\rho(t) - \rho_t) \delta(\rho(0) - \rho_0) \, d[\rho(\tau)] \\ &\asymp \int e^{-V} \int_0^t L(\lambda; \rho) d\tau \delta[\dot{\rho} + \nabla \cdot \lambda] \delta(\rho(t) - \rho_t) \delta(\rho(0) - \rho_0) \, d[\lambda(\tau)]. \end{aligned} \qquad (7.35)$$

More generally, we can write the expectation of an observable $O_t$ at time $t$ starting from a distribution $P_0$, which we both write at the large deviations scale in terms of two rate functions $U$ and $\theta$

$$P_0(\rho) \asymp e^{-VU(\rho)} \quad \text{and} \quad O_t(\rho) \asymp e^{-V\theta(\rho)}. \qquad (7.36)$$

We get the following path integral as our generic object of study (many objects of interest can be written in that form):

$$\begin{aligned} \langle \mathcal{O}_t \rangle_{P_0} &\asymp \int e^{-V(\int \mathcal{L}(\dot{\rho}(\tau); \rho(\tau)) d\tau + U(\rho_0) + \theta(\rho_t))} d[\rho(\tau)] \\ &\asymp \int e^{-V(\int L(\lambda(\tau); \rho(\tau)) d\tau + U(\rho_0) + \theta(\rho_t))} \delta[\dot{\rho} + \nabla \cdot \lambda] d[\lambda(\tau)]. \end{aligned} \qquad (7.37)$$

As stated earlier, we will mostly work inside the exponentials. The arguments of those two last exponentials are called actions (respectively the standard action and the detailed action), and they have a lot of structure, which we will look into next.

7.1.2 Properties of the action

Looking at the expressions above, we must remember that $V$ is large, so that the path with least action exponentially dominates the path integral. We will describe this dominating path with the usual tools of analytical mechanics, and we will see why we called our LDFs Lagrangians and our SCGFs Hamiltonians.
Equations of motion

The least action path can be characterised through a functional differentiation. For the standard action

$$S[\dot{\rho}(\tau), \rho_0] = \int_{\tau=0}^{t} \mathcal{L}(\dot{\rho}(\tau); \rho(\tau))d\tau + U(\rho_0) + \theta(\rho_t), \quad (7.38)$$

the differentiation of the action with respect to $\rho(\tau)$ at every time $\tau$, with an infinitesimal displacement $\delta\rho(\tau)$, gives:

$$\int_{\tau=0}^{t} (\partial_{\rho}\mathcal{L} \delta\rho(\tau) + \partial_{\dot{\rho}}\mathcal{L} \delta\dot{\rho}(\tau)) d\tau + \partial_{\rho}U(\rho_0)\delta\rho_0 + \partial_{\rho}\theta(\rho_t)\delta\rho_t = 0. \quad (7.39)$$

It is preferable to get rid of the term $\delta\dot{\rho}(\tau)$, which we can do through an integration by parts, which yields

$$\int_{\tau=0}^{t} \left( \partial_{\rho}\mathcal{L} - \frac{d}{dt}\partial_{\dot{\rho}}\mathcal{L} \right) \delta\rho(\tau) \, d\tau + \partial_{\rho}\mathcal{L} \delta\rho(t) - \partial_{\dot{\rho}}\mathcal{L} \delta\rho(0) + \partial_{\rho}U(\rho_0)\delta\rho_0 + \partial_{\rho}\theta(\rho_t)\delta\rho_t = 0. \quad (7.40)$$

Since the expression on the left-hand side needs to vanish regardless of the choice of displacement, the coefficients of $\delta\rho(\tau)$ must vanish separately at every time. Cancelling the first integrated term will yield the standard Euler-Lagrange equation, which justifies our calling $\mathcal{L}$ a standard Lagrangian:

$$\partial_{\rho}\mathcal{L} - \frac{d}{dt}\partial_{\dot{\rho}}\mathcal{L} = 0. \quad (7.41)$$

Canceling the boundary terms fix the boundary conditions:

$$\partial_{\dot{\rho}}\mathcal{L}(\dot{\rho}_0; \rho_0) = \partial_{\rho}U(\rho_0) \quad \text{and} \quad \partial_{\dot{\rho}}\mathcal{L}(\dot{\rho}_t; \rho_t) = -\partial_{\rho}\theta(\rho_t). \quad (7.42)$$

Note that, in the case where the initial or final condition is fixed to a single state $\tilde{\rho}$, the corresponding equation becomes irrelevant and the boundary condition should be replaced by $\rho_0 = \tilde{\rho}$ or $\rho_t = \tilde{\rho}$.

For the detailed action, the continuity constraint is enforced through a Lagrange multiplier $\mu$ (e.g. from a Fourier transform of the delta function):

$$S_{\lambda, \mu}[\dot{\rho}(\tau), \rho_0] = \int_{\tau=0}^{t} L(\lambda(\tau); \rho(\tau))d\tau + U(\rho_0) + \theta(\rho_t) + \int_{\tau=0}^{t} \mu(\tau)(\dot{\rho}(\tau) + \nabla \cdot \lambda(\tau))d\tau. \quad (7.43)$$

The minimisation of the integral terms with respect to the three independent variables $\rho$, $\lambda$ and $\mu$ yields three independent parts:

$$\delta[L + \mu(\dot{\rho} + \nabla \cdot \lambda)] = [\partial_{\rho}L \delta\rho + \mu \delta\dot{\rho}] + [\partial_{\lambda}L \cdot \delta\lambda + \mu \nabla \cdot \delta\lambda] + [\dot{\rho} + \nabla \cdot \lambda] \delta\mu = 0. \quad (7.44)$$

Each part has to vanish separately, yielding three equations. As above, we can handle the first $\partial_{\rho}L \delta\rho + \mu \delta\dot{\rho} = 0$ through an integration by parts on time, whose boundary terms will end up in the boundary conditions. The second $\partial_{\lambda}L \cdot \delta\lambda + \mu \nabla \cdot \delta\lambda = 0$ can be handled by transferring $\nabla$ to the left, as a gradient, i.e. applying $\mu \nabla \cdot \delta\lambda = (-\nabla\mu) \cdot \delta\lambda$. These two manipulations yield two equations, in addition to the continuity equation from the third term $\dot{\rho} + \nabla \cdot \lambda = 0$:

$$\partial_{\rho}L = \dot{\mu} \quad \text{and} \quad \partial_{\lambda}L = \nabla\mu. \quad (7.45)$$
The second equation allows us to get rid of $\mu$ in any equation, at the cost of applying $-\nabla$ to the rest of the equation. For instance, combining the two equations above yields the detailed Euler-Lagrange equation

$$ \nabla \partial_{\rho} L + \frac{d}{dt} \partial_{\lambda} L = 0. \qquad (7.46) $$

Moreover, the boundary terms from the partial integration of $\mu \delta\dot{\rho}$, combined with the differentials of $U_0$ and $\theta_t$, give the boundary conditions

$$ \mu(0) = \partial_{\rho} U(\rho_0) \quad \text{and} \quad \mu(t) = -\partial_{\rho} \theta(\rho_t) \qquad (7.47) $$

so that

$$ \partial_{\lambda} L(\lambda_0; \rho_0) = \nabla \partial_{\rho} U(\rho_0) \quad \text{and} \quad \partial_{\lambda} L(\lambda_t; \rho_t) = -\nabla \partial_{\rho} \theta(\rho_t). \qquad (7.48) $$

Note that, being LDFs, both Lagrangians vanish with zero derivative at their most probable values:

$$ \mathcal{L}(\dot{\rho}^*(\rho); \rho) = L(\lambda^*(\rho); \rho) = 0 \quad \text{or} \quad \partial_{\dot{\rho}} \mathcal{L}(\dot{\rho}^*(\rho); \rho) = \partial_{\lambda} L(\lambda^*(\rho); \rho) = 0 \qquad (7.49) $$

**Example** No example here, calculations are too messy. Everything is nicer with Hamiltonians.

The Hamiltonians are the SCGFs of the Lagrangians, i.e. their Legendre transforms: in the detailed case,

$$ H(f; \rho) = \sup_{\lambda} [f \cdot \lambda - L(\lambda; \rho)] \quad \text{and} \quad L(\lambda; \rho) = \sup_{f} [f \cdot \lambda - H(f; \rho)] \qquad (7.50) $$

so that, if $L$ is differentiable and convex with respect to $\lambda$,

$$ H(f; \rho) + L(\lambda; \rho) = f \cdot \lambda \quad \text{with} \quad f = \partial_{\lambda} L \quad \text{or} \quad \lambda = \partial_{f} H, \qquad (7.51) $$

and in the standard case,

$$ \mathcal{H}(u; \rho) = \sup_{u} [u \dot{\rho} - \mathcal{L}(\dot{\rho}; \rho)] \quad \text{and} \quad \mathcal{L}(\dot{\rho}; \rho) = \sup_{\dot{\rho}} [u \dot{\rho} - \mathcal{H}(u; \rho)], \qquad (7.52) $$

so that, if $\mathcal{L}$ is differentiable and convex with respect to $\dot{\rho}$,

$$ \mathcal{H}(u; \rho) + \mathcal{L}(\dot{\rho}; \rho) = u \dot{\rho} \quad \text{with} \quad u = \partial_{\dot{\rho}} \mathcal{L} \quad \text{or} \quad \dot{\rho} = \partial_{u} \mathcal{H}. \qquad (7.53) $$

This automatically justifies their name, but let us still check that they verify nice equations.

The equations for optimal paths in terms of Hamiltonians can be found thus: first, from the Legendre transform, we have the force-flux relations

$$ \lambda = \partial_{f} H \quad \text{or} \quad \dot{\rho} = \partial_{u} \mathcal{H}. \qquad (7.54) $$

Second, we rewrite the Euler-Lagrange equation in terms of the Hamiltonian. For this, we first need to differentiate the Legendre relation with respect to $\rho$ to get

$$ \partial_{\rho} \mathcal{H} + \partial_{\rho} \mathcal{L} + \partial_{\rho} \dot{\rho} \partial_{\dot{\rho}} \mathcal{L} = u \partial_{\rho} \dot{\rho} = \partial_{\dot{\rho}} \mathcal{L} \partial_{\rho} \dot{\rho} \qquad (7.55) $$
so that we immediately get the Hamilton equations (where we recall the force-flux relation as one of the two equations):

$$ \dot{u} = -\partial_{\rho}\mathcal{H} \quad \text{with} \quad \dot{\rho} = \partial_{u}\mathcal{H}, \qquad (7.56) $$

and in the same way

$$ \dot{f} = \nabla\partial_{\rho}H \quad \text{with} \quad \dot{\rho} = -\nabla \cdot \partial_{f}H. \qquad (7.57) $$

The boundary conditions simply translate to

$$ u(0) = \partial_{\rho}U(\rho_{0}) \quad \text{and} \quad u(t) = -\partial_{\rho}\theta(\rho_{t}). \qquad (7.58) $$

and

$$ f(0) = \nabla\partial_{\rho}U(\rho_{0}) \quad \text{and} \quad f(t) = -\nabla\partial_{\rho}\theta(\rho_{t}). \qquad (7.59) $$

Both Hamiltonians have a few useful properties:

* $H$ and $\mathcal{H}$ are conserved by the dynamics.

* $H(0; \rho) = 0$ and $\mathcal{H}(0; \rho) = 0$.

* $H$ is convex in $u$ and $\mathcal{H}$ is convex in $f$.

### Example

for the replicator, we get

$$ \dot{u} = -a(e^{u} - 1) - 2b\rho(e^{-u} - 1) - c(e^{-u} - 1) \qquad (7.60) $$

with $\dot{\rho} = a\rho e^{u} - b\rho^{2}e^{-u} - c\rho e^{-u}$,

$$ [\dot{f}_{21}, \dot{f}_{12}, \dot{f}_{01}] = (a(e^{f_{21}} - 1) + 2b\rho(e^{f_{12}} - 1) + c(e^{f_{01}} - 1)) [-1, 1, 1] \qquad (7.61) $$

with $\dot{\rho} = a\rho e^{f_{21}} - b\rho^{2}e^{f_{12}} - c\rho e^{f_{01}}$.

One last thing we can do to be complete in our parallel with analytical mechanics is to examine the evolution of the static LDF at time $t$, i.e. the function $g_t(\rho_t)$ which is the logarithm of the probability of states at a fixed time $t$, starting, for instance, from an initial condition $\rho_0$ at $t = 0$. This LDF is precisely the action of the optimal path between $(\rho_0, t = 0)$ and $(\rho_t, t)$, i.e.

$$ g_t(\rho_t) = \int_{0,\rho_0}^{t} \mathcal{L}(\dot{\rho}^*; \rho^*) d\tau = \int_{0,\rho_0}^{t} L(\lambda^*; \rho^*) d\tau \equiv S_t^*(\rho_t, \rho_0). \qquad (7.62) $$

Differentiating with respect to time $t$, and being careful to propagate the action over the infinitesimal time difference, we get the expected Hamilton-Jacobi equation

$$ \partial_t g_t(\rho) = -\mathcal{H}(\partial_\rho g_t(\rho); \rho). \qquad (7.63) $$

Note: varying the initial time instead yields a second Hamilton-Jacobi equation, with opposite signs. The computations are left as an exercise for the enthusiastic reader.
### What it means

At this point, it is worth taking a step back from the mathematics to look at the physics of the systems we typically want to describe. For a system of particles diffusing in a fluid, the randomness in the current of particles has a very specific source: the collisions between the particles and the fluid molecules. At every point in time, these collisions exert a net force per unit time on the particles, which is random, and whose distribution is related to the momentum distribution of the fluid molecules (as characterised by their temperature) but is independent of the properties of the particles. Those properties only come into play when the random force produced by the fluid is translated into a random current of particles, depending in particular on their mobility (i.e. the ease with which they move in a medium).

The point here is the following: the dynamical variable $f$ of the detailed Hamiltonian can be interpreted as the random force produced by the medium. The relation between that force and the resulting current, called the force-flux relation, is simply

$$ \lambda = \partial_f H(f; \rho) \qquad (7.64) $$

up to a possible shift of the reference $f = 0$. In particular, when there are fixed external forces acting on the particles as well, they will appear additively with $f$ in this relation. This is particularly clear for systems with detailed balance, where the force deriving from the free energy of a system turns up as a special value of $f$. In more general settings, this can instead be taken as the definition of a force, being the quantity conjugate to currents through Legendre transforms.

Finally, the variable $u$ of the standard Hamiltonian can be interpreted as (minus) a random potential, in situations where all random forces derive from potentials.

### 7.1.3 Long time phenomenology

The framework we have presented is quite versatile: all kinds of observables can be analysed, or at least expressed in terms of solutions to known equations. It is however often complicated to get explicit results. One thing that usually makes things simpler is to look at long times, where systems converge to characteristic behaviours that are more clearly interpretable. For simplicity, we focus solely on the transition probability from $\rho_0$ to $\rho_t$, as an example to get an intuition on what features of the model are relevant for its long time behaviour.

#### Hamiltonian flow

Let us first recall the optimal value of the action over a finite time $t$, which tells us which path is the most likely, and what its probability is. It is useful to write it in terms of the Hamiltonian:

$$ \begin{aligned} S_t^*(\rho_t, \rho_0) & \equiv \lim_{V \to \infty} -\frac{1}{V} \ln(P_t(\rho_t|\rho_0)) = \inf_{\rho} \left[ \int_{0;\rho_0}^{t;\rho_t} \mathcal{L}(\dot{\rho}; \rho) d\tau \right] = \int_{0;\rho_0}^{t;\rho_t} \mathcal{L}(\dot{\rho}^*; \rho^*) d\tau \\ & = -t\mathcal{H}^* + \int_{0;\rho_0}^{t;\rho_t} u^* \dot{\rho}^* d\tau \end{aligned} \qquad (7.65) $$

where the asterisks indicate optimal values. The path we need to plug into the integral is a solution to the Hamilton equations with the correct duration $t$.

General calculations are difficult in this case, but a simple example can teach us a lot.
### Example

Let us look once again at the replicator. Since this is a 1D model, the solutions to the Hamilton equations are simply the level lines of the Hamiltonian, with appropriate boundary conditions:

$$ \mathcal{H}(u; \rho) = (a\rho - (b\rho^2 + c\rho)e^{-u})(e^u - 1) = \text{cst.} \qquad (7.66) $$

We can easily study those trajectories by drawing the phase portrait of this Hamiltonian. The boundary conditions will correspond to vertical lines $\rho = \text{constant}$, and we still need to select the trajectory (or trajectories) which join those two lines in the prescribed duration $t$. All of this is represented on the figure below, for three different durations, along with the shape of the trajectory $\rho(t)$ as a function of time.

<figure><img src="image_59.png" /><figcaption>Figure 7.1: Left: sketch of Hamiltonian trajectories close to a fixed point, for increasing times (blue to green to red). Right: evolution of $\rho(t)$ for the same trajectories.</figcaption></figure>

### Convergence of trajectories

The time that each trajectory takes is not a simple function of its length in phase-space, but the only way to produce a long trajectory is to get closer to the fixed point, where the velocities vanish. As a result, the bulk of the trajectory converges to the fixed point, with only finite boundary layers connecting it to each boundary condition. This convergence is guaranteed purely by the boundary conditions, and the critical point can be rightly called an attractor. The right plot of 7.1 shows the evolution of $\rho(t)$ for those same trajectories, with the grey one representing the asymptotic limit for infinite time (where the middle horizontal part should in fact be infinitely long). Each trajectory shows a plateau at a value that converges to $\rho^*$, of increasing length, and two boundary layers whose shapes also converge and are bounded between $\rho_0$ and $\rho_t$.

Although this is just a fairly trivial example with a single scalar state variable, it is somewhat representative of typical stochastic processes, due to a few generic properties of their Hamiltonians, or certain reasonable assumptions. Here is an overview of some relevant considerations:

* The manifold $u = 0$ is stable under the dynamics, and contains the deterministic evolution of the system, including all the attractors of that evolution. This restricted dynamics $\dot{\rho} = \partial_u \mathcal{H}(0; \rho)$ is not Hamiltonian on its own, and is allowed to converge even though the appropriate boundary conditions are just $\rho(0) = \rho_0$.
* For the large deviations principle to hold, we need to ensure that the trajectories we consider stay away from the boundaries of phase-space (e.g. from $n = 0$ in the example above, since we need $n$ to be of order $V$). A reasonable assumption to guarantee that is to consider only systems that are *globally stable*, meaning that their deterministic dynamics are such that they always go back towards a given compact region of phase-space. In one dimension, this means that $\dot{\rho}$ must be positive for $\rho$ small enough, and negative for $\rho$ large enough. In general, this implies that all heteroclines (i.e. trajectories going from one attractor to another) outside that region are stable manifolds, so the system cannot diverge towards infinity in any direction, and the Poincaré-Bendixon theorem then guarantees that said region contains at least one stable attractor. With density-type boundary conditions for the Hamilton equations, this guarantees convergence in the long time limit.

* In many simple and physically motivated models, all attractors belong to the $u = 0$ manifold. This means that the value of $H^*$ in the principal function (7.65) will necessarily be 0. There can also be closed trajectories with $H \neq 0$, which are not attractors of the deterministic dynamics, but all known examples are such that $H < 0$ so that they are not selected by the infimum.

* Those attractors are typically connected in two ways: through the $u = 0$ manifold, which contains the deterministic trajectories between them, and through another manifold $u^*(\rho)$ such that $H(u^*(\rho); \rho) = 0$, which contains less probable trajectories called *instantons*, of which we give an example below.

To summarise, the long-time behaviour of the Hamilton equations for stochastic processes is mostly characterised by the attractors of the deterministic dynamics along with their stable and unstable manifolds, which all verify $H = 0$. In practice, the only type of attractor that is manageable with current methods is fixed points (and some simple limit cycles).

## Computations of the extinction time in the replicator

To illustrate those last few points, let us answer a specific long time question: what is the probability for the replicator to go extinct in the long run, i.e. to reach density $\rho = 0$. To make things easier, we can take $c$ very small: this parameter was here to avoid having an absorbing state, but here we don't mind.

The special manifolds of $H$ verify

$$H(u; \rho) = 0 \quad \Rightarrow \quad u = 0 \quad \text{or} \quad u = \ln\left(\frac{b\rho^2 + c\rho}{a\rho}\right) \sim \ln\left(\frac{b\rho}{a}\right). \quad (7.67)$$

The probability rate to reach $\rho = 0$ along the escape trajectory, starting from the steady state $\rho = a/b$, is then given in terms of the corresponding action

$$S = \int u(\rho) \dot{\rho} d\tau = \int_{a/b}^{0} u(\rho) d\rho = \frac{a}{b} \quad (7.68)$$

so that

$$P_t(\rho_t = 0 | \rho_0 = a/b) \asymp e^{-Va/b} \quad (7.69)$$

so that the extinction timescale is its inverse, $\tau = e^{Va/b}$, and is exponentially large in $V$.
We can also compare this result to what we would get in the Chemical Langevin Equation approximation, to see how good of an approximation it is. The approximate Hamiltonian is

$$ \mathcal{H}(u; \rho) \sim u \left( a\rho - b\rho^2 - c\rho + \frac{u}{2} (a\rho + b\rho^2 + c\rho) \right) \qquad (7.70) $$

so that

$$ \mathcal{H}(u; \rho) = 0 \quad \Rightarrow \quad u = 0 \quad \text{or} \quad u = 2 \frac{b\rho + c - a}{b\rho + c + a} \sim 2 \frac{b\rho - a}{b\rho + a}. \qquad (7.71) $$

The probability rate to reach $\rho = 0$ along the escape trajectory, starting from the steady state $\rho = a/b$, is given in terms of the corresponding action

$$ S = \int_{a/b}^{0} u(\rho)d\rho = \frac{a}{b}(\ln(4) - 1) \approx 0.386 \frac{a}{b}, \qquad (7.72) $$

so that

$$ P_t(\rho_t = 0 | \rho_0 = a/b) \asymp (4/e)^{-Va/b} \qquad (7.73) $$

so that the extinction timescale is its inverse, $\tau = (4/e)^{Va/b}$. This is exponentially longer than the real value.

## 7.2 Path integral formalisms in physics: a comparison

To conclude, let us compare the path integral formalism we have described here with its equivalents within classical and quantum mechanics. There are many similarities in structure, but the specificities of each version make all the difference.

Let us first recap the relevant characteristics of the formalism in statistical mechanics, which is built starting from the Lagrangian:

* The Lagrangian is the logarithm of a classical probability distribution: it is positive and vanishes in the most likely state. The Lagrangian is physical, whereas the Hamiltonian is effective (it has no physical interpretation, as it is merely a generating function).

* The Hamiltonian vanishes along the most likely trajectory, and is constant along other local extrema, but other trajectories that do not verify the Hamilton equations are possible due to statistical fluctuations.

* The dynamical variable of the Hamiltonian is a random force or a potential.

* The boundary conditions of the Hamilton equations are mixed, allowing for convergence even with conservative dynamics.

In contrast, the formalism for quantum mechanics is built starting from the Hamiltonian:

* The Hamiltonian is the logarithm of a quantum probability distribution: it is positive and vanishes in the most likely state. It is physically interpretable as an energy, whereas the Lagrangian is effective.

* The Lagrangian vanishes along the typical trajectory, but other trajectories that do not verify the Euler-Lagrange equations are possible due to quantum fluctuations.

* The dynamical variable of the Hamiltonian is a physical momentum.
* The boundary conditions of the Hamilton equations are full initial conditions: trajectories cannot converge in time due to unitary dynamics.

In analytical mechanics, the starting point is the conservation in time of the total energy of a system:

* The Hamiltonian is an energy, and is usually positive. Its conservation in time leads to the Hamilton equation. Only the trajectories that verify the Hamilton equations are physical.

* The Lagrangian is effective, and is merely a convenient object with which to write the equations of motion of the system. Only the trajectories verifying the Euler-Lagrange equations are physical, as there are no fluctuations.

* The dynamical variable of the Hamiltonian is a physical momentum.

* The boundary conditions of the Hamilton equations are full initial conditions: trajectories cannot converge in time due to conservative dynamics.

For more details on large deviation theory and its applications in physics, along with illustrations, exercises, and moderate ranting about the definitional nightmare that is thermodynamics, the reader may refer to *Large deviations in statistical mechanics: from free energies to path integrals*, a set of lecture notes by myself (Alexandre Lazarescu), which will appear on arXiv some day.
# 8 Metastability In and Out of Equilibrium

Gianmaria Falasco and Massimo Bilancioni.<fnref n="21" /> Metastability indicates the existence of long-lived, yet quasi-stationary states, possibly maintained by continuous energy dissipation. Examples are the logic state of an electronic circuit and the homeostatic state of a biochemical reaction network. We review the definition of metastability based on the spectral theory of the generator of a Markov process. We present two prototypical examples representing jump processes in the large-size limit and diffusive dynamics in the small-temperature limit. For the latter, we introduce the WKB expansion of the Fokker-Planck equation, the quasi-potential and its relation with the life-time of metastable states. For detailed balance dynamics, we retrieve the Arrhenius law from the standpoint of large deviations theory. For dynamics subjected to nonconservative forces, we show the relation between the life-time of metastable states and dissipation.

## 8.1 Introduction

### 8.1.1 Intuitive idea

The existence of various dynamical regimes, widely separated by distinct timescales is called metastability [105]. In the simplest case, metastability manifests itself when a process seems stationary for a very long time before jumping to a new seemingly stationary state. Such states are local minima of the free energy for nondissipative systems [106], i.e. whose dynamics are detailed balanced. What are they, when detailed balance is not satisfied? How to characterize metastability in terms that are intrinsically dynamic and make no a priori reference to geometric concepts such as free energy landscapes? These questions are crucial far from equilibrium where the free energy and the stationary probability are no longer related to each other in a simple way.

### 8.1.2 Examples

To set the stage we list some notable examples of metastability the reader might be familiar with. Detailed balanced systems include:

* the Ising model (in contact with a single thermal bath) below its critical temperature. Upwards and downwards magnetizations are the two metastable states.

* Carbon under standard conditions in the diamond state. It will eventually relax to a lower free energy, i.e graphite.

* A supercooled liquid, i.e. a liquid cooled fast below the melting point in the absence of crystal nucleators, e.g. impurities. It will eventually solidify [107].

Metastable dynamics that are not detailed balanced comprise

* Climate. It has (at least) two metastable states, warm and snowball—paleoclimatology identified 2 full glaciations and tens of ice-ages. Transitions are induced by (weather) noise and (astronomical) time-periodic factors [108].

* Chemical reaction networks: gene regulatory networks have multiple phenotypic (or epigenetic) states which are metastable—each state is characterized by the activation of different gene patterns, resulting in different protein concentrations. Intrinsic (i.e. finite copy numbers) and extrinsic (i.e. variability) noise cause transitions [109].

<footnote n="21">²¹GF was the lecturer and wrote this chapter; MB was the angel.</footnote>
* Electronic circuits: in a bit-storage element (CMOS inverters) 0 and 1 states are metastable. Thermal noise induces transitions, resulting in computation errors [110].

## 8.2 Metastability for Markov processes

The phenomenon is understood within the spectral theory of the Markovian generator [111–114], which is summarized as follows. The dynamical equation for the probability density vector $|P_t\rangle$ at time t is

$$ \partial_t |P_t\rangle = R |P_t\rangle \qquad (8.1) $$

where the operator R is the generator of the stochastic process. Think of R as the irreducible rate matrix<sup>22</sup> for a Markov jump process on a finite state space of dimension N. States can be labelled by n such that $\langle n|P_t\rangle = P(n, t)$. The scalar product is standard matrix multiplication.

### 8.2.1 Spectral decomposition

The time evolution of the probability distribution can be expanded in the right eigenvectors $|\psi_{\alpha}^{(r)}\rangle$ of the operator R, with eigenfunctions $\psi_{\alpha}^{(r)}(n) = \langle n|\psi_{\alpha}^{(r)}\rangle$, defined by the equations

$$ R |\psi_{\alpha}^{(r)}\rangle = \lambda_{\alpha} |\psi_{\alpha}^{(r)}\rangle \qquad \langle \psi_{\alpha}^{(l)} | R = \lambda_{\alpha} \langle \psi_{\alpha}^{(l)} | . \qquad (8.2) $$

Therefore, we can write

$$ |P_t\rangle = \sum_{\alpha} b_{\alpha} |\psi_{\alpha}^{(r)}\rangle e^{\lambda_{\alpha} t}, \qquad (8.3) $$

where $b_{\alpha} = \langle \psi_{\alpha}^{(l)} | P_0 \rangle$ is the overlap of the initial condition with the left eigenvector $\langle \psi_{\alpha}^{(l)} |$. Since R is not Hermitian in general<sup>23</sup>, left and right eigenvectors are different but still form an orthonormal basis

$$ \langle \psi_{\alpha}^{(l)} | \psi_{\alpha'}^{(r)} \rangle = \delta_{\alpha\alpha'} \qquad I = \sum_{\alpha} |\psi_{\alpha}^{(r)}\rangle \langle \psi_{\alpha}^{(l)}|, \qquad (8.4) $$

and the eigenvalues $\lambda_{\alpha}$ are not necessarily real. However, since R is an irreducible stochastic matrix<sup>24</sup>, by the Perron-Frobenius theorem<sup>25</sup> its eigenvalue $\lambda_0 = 0$ with largest real part

* is nondegenerate,

* is associated to the left eigenfunction $\langle \psi_0^{(l)} | = \langle - |$, i.e.

$$ 0 = \langle - | R, \qquad (8.5) $$

with $\langle - | n \rangle$ constant,

* and to a right eigenfunction with positive components, which is the stationary probability distribution, i.e. $|\psi_0^{(r)}\rangle = |P_{\infty}\rangle$.

Hence, R has N nonpositive eigenvalues $\lambda_{\alpha}$ which can be ordered by their real part $\text{Re}(\lambda_{\alpha}) \ge \text{Re}(\lambda_{\alpha+1})$.

<footnote n="22"><sup>22</sup>The graph is connected, all states can be reached and so the dynamics is ergodic.</footnote>

<footnote n="23"><sup>23</sup>This means that in general it cannot even be diagonalized, it can only be put in a Jordan normal form. This happens if, for some values of the system's parameters (called exceptional points), there are degenerate eigenvalues. We assume that the smaller eigenvalues are not degenerate.</footnote>

<footnote n="24"><sup>24</sup>Actually, $e^R$ is a stochastic matrix, since it has only positive entries and $\langle - | e^R = \langle - |$ because of probability normalization.</footnote>

<footnote n="25"><sup>25</sup>Or its infinite-dimensional generalization, the Krein-Rutman theorem.</footnote>
<figure><img src="image_60.png" /><figcaption>Figure 8.1: Sketch of the spectrum of R (inverse of the eigenvalues) for a system displaying metastability.</figcaption></figure>

### 8.2.2 Dynamical identification of metastable states

We say that the generator depends on a parameter $\epsilon$ which we can imagine as the temperature, the inverse size of the system, etc. Metastability appears if there exist eigenvalues $\lambda_1, \dots, \lambda_M$ whose real part goes to 0 = $\lambda_0$ as $\epsilon \to 0$, while all others $\lambda_\alpha$ with $\alpha > M$ stay finite. We can then identify two well separated timescales:

* $t_{esc} = \min(\lambda_1^{-1}, \dots, \lambda_M^{-1})$ is the shortest escape time out of metastable states

* $t_{rel} = \max(\lambda_{M+1}^{-1}, \dots, \lambda_N^{-1})$ the longest time to relax within a metastable state.

The gap in the spectrum, whose inverse is related to $t_{esc} - t_{rel}$, corresponds to a diverging time scale separating the (fast relaxation) dynamics within metastable states from the slow (jump dynamics) between them.

This intuitive picture can be made precise. It can be shown [111] that on the intermediate timescale $t_{rel} \ll t \ll t_{esc}$ one can find

* a basis of $M+1$ right eigenvectors $|\rho_i\rangle$ which are positive normalized, stationary, i.e $R|\rho_i\rangle \simeq 0$, nonzero only on non-overlapping regions of the phase space;

* a basis of $M+1$ left eigenvectors $\langle q_i|$ such that $\langle q_i| R \simeq 0$ and satisfying the orthogonality condition $\langle q_i|\rho_j\rangle \simeq \delta_{ij}$.

Basically, $\rho_i(n)$ are the probability distributions of the metastable states and $q_i(n)$ are nonzero only on their basin of attraction (decaying rapidly to zero outside). Moreover, $q_i(n)$ gives the probability of reaching the metastable state $i$ starting from $n$ on the intermediate timescale. It can be used to define the transition state out of the metastable state $i$, which is located at $q_i(n) = 1/2$. They are called committor functions [115]. Note that they both are linear combinations of the eigenvectors $\langle \psi_\alpha^{(l)}|$ and $|\psi_\alpha^{(r)}\rangle$, which allow us to write the evolution operator as a $e^{tR} \simeq \sum_i |\rho_i\rangle \langle q_i|$ on the intermediate timescale $^{26}$.

Ultimately, if $\epsilon \to 0$ at finite $t$, the system probability can only converge to those $\rho_i$ that are selected by the initial condition, i.e. with $\langle q_i|P_0\rangle \neq 0$. Therefore, the long-time limit and the small-$\epsilon$ limit do not commute, because if $t \to \infty$ at any finite $\epsilon$ the probability distribution converges surely to the unique $|\psi_0^{(r)}\rangle$. Historically, this fact is known as Keizer's paradox [80, 116].

<footnote>$^{26}$Everywhere in this subsection the symbol $\simeq$ means that we ignore terms of order $e^{t\lambda_\alpha}$ with $\alpha > M$.</footnote>
### 8.2.3 Example: the Schlögl model

We consider the Schlögl model [80,117], an autocatalytic chemical system displaying bistability far from equilibrium. It consists in the following set of chemical reactions,

$$2X + A \underset{k_{-1}}{\overset{k_{+1}}{\rightleftharpoons}} 3X; \quad B \underset{k_{-2}}{\overset{k_{+2}}{\rightleftharpoons}} X. \qquad (8.6)$$

Here, species A and B will be considered to be chemostated reservoirs (i.e., with constant concentration), and the only degree of freedom for the problem is the number $n$ of molecules of species X. The reactions are assumed to take place in a well-mixed container of volume $V$. If $V$ is large in comparison to the molecular scale, its inverse can be taken to be the small parameter $\epsilon = 1/V$. Thus, in the macroscopic limit we deal with the concentration $x = n/V = n\epsilon$ of species X, and the scaled (by $V$) transition rates follow the mass-action law (valid for dilute ideal gasses):

$$r_1(x) = k_{+1}ax^2, \quad r_{-1}(x) = k_{-1}x^3, \quad r_2(x) = k_{+2}b, \quad r_{-2}(x) = k_{-2}x. \qquad (8.7)$$

Expanding the chemical master equation, multiplying by $x$ and integrating, in the limit $\epsilon \to 0$ (in which the probability peaks around the most likely value $\langle x \rangle$) we obtain the rate equation of chemistry

$$\dot{x} = r_1(x) - r_{-1}(x) + r_2(x) - r_{-2}(x). \qquad (8.8)$$

In figure (8.2) we show for the inverse of the volume $\epsilon = 1/30 \simeq 3.3 \times 10^{-2}$ the first two right and left eigenfunctions obtained numerically with Mathematica. Clearly, from their linear combination we can obtain the (positive, normalized) metastable states and the committors (1 on the respective attractor). Partial overlapping of the metastable states and the discrepancies with the $\epsilon = 0$ predictions are due to the finite value of $\epsilon$.

In figure (8.3) we show the absolute value of the difference between consecutive leading eigenvalues. Even at the moderately small value $\epsilon = 1/30$ it is apparent that $\lambda_1 \to 0$ while the other eigenvalues stay finite.

### 8.2.4 Detailed balance vs. nondetailed balance dynamics

The first important point is that a system with a convex energy function $U$ can develop metastable states when it is subjected to nonpotential forces. Such states will require continuous energy dissipation to exist (see examples above). A second point is that in the absence of detailed balance $R$ is not Hermitian, so the eigenvalue $\lambda_\alpha$ can have a nonzero imaginary part. Thus, a metastable state can be periodic, e.g. a stable limit cycle is reached as $\epsilon \to 0$ [118, 119]. Nevertheless, we do not explicitly consider such periodic metastable states hereafter, or more general quasi-periodic and chaotic ones. In general, $R$ is symmetrizable only if detailed balance holds. In this case the leading right eigenfunction is the equilibrium Gibbs distribution, e.g. $\psi_0^{(r)}(n) \propto e^{-\beta U(n)}$ in the canonical ensemble at inverse temperature $\beta$, such that a change of basis defines the symmetric operator $R_{sym} = e^{\beta U/2} Re^{-\beta U/2} = R_{sym}^\dagger$.

## 8.3 Overdamped Langevin dynamics: low temperature limit

The general framework for the description of metastability can be applied to systems with a continuous state space [120], which we analyze more explicitly in the following by using large deviations theory [94]. We consider the Langevin equation for $x \in \mathbb{R}^d$,

$$\dot{x} = F(x) + \sqrt{2\epsilon\xi} \qquad (8.9)$$
<figure><img src="image_61.png" /><figcaption>Figure 8.2: For $\epsilon^{-1} = 30$, first two right and left eigenvectors (in (a) and (b), respectively), and their linear combination giving the metastable states (c) and the committors (d). Vertical grey lines in the panels indicate the stable (lower left) and unstable (lower right) fixed points of the deterministic dynamics (8.8), corresponding to $\epsilon = 0$. The horizontal grey line in (d) is at $q_i = 1/2$ and singles out the transition state. Numerical evaluation is obtained by approximating $R$ with a $150 \times 150$ matrix (finite size effects are also visible at $n\epsilon \simeq 5$). Parameters are $k_1a = 4.1, k_{-1} = 1, k_2b = 1, k_{-2} = 4$.</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>Figure 8.3: For $\epsilon = 1/30$, absolute value of the difference between the first 5 consecutive eigenvalues. System's parameter as before.</figcaption></figure>
<figure><img src="image_63.png" /><figcaption>Figure 8.4: Example of a drift field composed by an underlying double well free energy potential $U(x)$ and a superimposed nonequilibrium forcing $f(x)$. The arrows represent schematically the drift’s direction and intensity.</figcaption></figure>

with (weak) noise, Gaussian with zero mean and correlation matrix $\langle \xi(t)\xi(0) \rangle = D(x)\delta(t)$ $^{27}$. The drift field has stable and unstable zeros, denoted $x_i^*$ and $x_\nu$, respectively. Namely, $F(x_i^*) = 0 = F(x_\nu)$ and with (assume nonsingular) Jacobian matrix $\nabla F(x)$ having only negative eigenvalues in $x_i^*$ (resp. at least a positive eigenvalue in $x_\nu$). Namely, $x_i^*$ are stabled fixed points of (8.9), while $x_\nu$ are saddle points.

We specify the drift field as a conservative part coming from an underlying free energy landscape $U(x)$ and a nonequilibrium forcing part $f$ (Fig. 8.4):

$$F(x) = -D(x) \cdot \nabla U(x) + f(x). \qquad (8.10)$$

We take an autonomous dynamics leading to a stationary probability density $P_\infty$ and current $J_\infty$, i.e. the stationary solution of the Fokker-Planck

$$\partial_t P(x,t) = -\nabla \cdot J(x,t) = -\nabla \cdot [F(x)P(x,t) - D(x)\epsilon \nabla P(x,t)]. \qquad (8.11)$$

Detailed balance dynamics corresponds to $f = 0$, in which case $x_i^*$ are the local minima of the energy $U$ and $x_\nu$ its saddles. This setup can describe, for example, interacting colloidal particles under a shear flow at low temperature$^{28}$.

### 8.3.1 Deterministic limit and typical fluctuations

For $\epsilon = 0$ the dynamics are deterministic, namely, $P(x,t) = \delta(x - x(t))$ where $x(t)$ is the solution of

$$\dot{x}(t) = F(x(t)). \qquad (8.12)$$

Equation (8.12) describes the relaxation to the stable fixed point $x_i^*$, where $i$ is selected by the initial condition $x(0) = x_0$. The dynamics are therefore nonergodic, with the state space partitioned into the various $i = 0, \dots, M$ basins of attractions.

To describe the low-$\epsilon$ regime, it is therefore natural to use the WKB ansatz

$$P(x,t) = W_{\epsilon}(x,t)e^{-\frac{1}{\epsilon}I(x,t)}, \qquad (8.13)$$

<footnote>$^{27}$All choices of stochastic calculus are equivalent to leading order in $\epsilon$. We will use Ito later on when we do calculations at leading order.</footnote>

<footnote>$^{28}$in the sense that $k_B T$ is small in comparison to the typical interaction energy and forcing.</footnote>
with $W_\epsilon(x, t)$ subexponential in $\epsilon$. In the language of large deviations theory,

$$I(x, t) = - \lim_{\epsilon \to 0} \epsilon \ln P(x, t), \qquad (8.14)$$

is called rate function and concentration of probability happens at speed $1/\epsilon$. The logarithmic equality (8.14) is usually written $P(x, t) \asymp e^{-\frac{1}{\epsilon}I(x,t)}$. Plugging (8.13) into (8.11) we obtain at the leading order

$$-\partial_t I(x, t) = H(x, \nabla I(x, t)) \qquad (8.15)$$

with Hamiltonian

$$H(x, \pi) = F(x) \cdot \pi + \pi \cdot D(x) \cdot \pi. \qquad (8.16)$$

Note that the Hamiltonian is the cumulant generating function of a Gaussian noise with mean $F(x)$ and covariance matrix $D(x)$. Equation (8.15) is a Hamilton-Jacobi equation for the action function $I(x, t)$ and conjugate momentum $\pi = \nabla I$. We can recast the problem in terms of the equivalent Hamiltonian equations

$$\dot{x} = F(x) + 2D(x) \cdot \pi, \quad \dot{\pi} = -\nabla F(x) \cdot \pi + \pi \cdot \nabla D(x) \cdot \pi. \qquad (8.17)$$

These equations describe the most likely trajectories followed by the system as $\epsilon \to 0$. They also follow from the minimization of the action functional$^{29}$

$$\mathcal{A}[\{x(t), \pi(t)\}_0^\tau] = \int_0^\tau dt(\pi(t) \cdot \dot{x}(t) - H(x(t), \pi(t))). \qquad (8.18)$$

From lecture 6 we know that such action determines the leading order probability of paths $\{x(t)\}_0^\tau$ as$^{30}$

$$\mathrm{Prob}[\{x(t)\}_0^\tau | x_0] \asymp \int \mathcal{D}\pi e^{-\frac{1}{\epsilon}\mathcal{A}[\{x(t), \pi(t)\}_0^\tau]}. \qquad (8.19)$$

The stationary rate function $I_\infty(x) = -\lim_{\epsilon \to 0} \epsilon \ln P_\infty(x)$ is obtained by solving

$$0 = H(x, \nabla I_\infty(x)), \qquad (8.20)$$

which means that we look for solution of (8.17) on the submanifold $H = 0$. If $f \neq 0$, $I_\infty$ generally has nondifferentiable points, so that (8.20) should be solved within each basin on attraction, giving a local rate function [124].

### 8.3.2 Instantons and escape rates

From the probability current introduced in (8.11), we can also derive another decomposition of the drift field

$$F(x) = -D(x)\nabla I_\infty(x) + v(x) \qquad (8.21)$$

<footnote>$^{29}$The boundary conditions are determined by the initial distribution $P(x, 0)$ and possibly from a constrain at the time of interest $\tau$.</footnote>

<footnote>$^{30}$This can be derived in a number of more or less rigorous ways. For example, starting from the Langevin equation and the Gaussian path probability of the noise $\{\xi(t)\}_{t=0}^\tau$ as in the MSRDJ formalism [121–123], or à la Feynman by time-slicing a formal solution of the Fokker-Planck equation [18].</footnote>
that does not make reference to the underlying microscopic details (energy and forcing), but rather uses dynamically emergent quantities, i.e the stationary rate function and the probability velocity

$$v(x) = \lim_{\epsilon \to 0} \frac{J_{\infty}(x)}{P_{\infty}(x)} \qquad v \cdot \nabla I_{\infty} = 0. \qquad (8.22)$$

The last property follows immediately from the substitution of (8.21) into (8.20). Note that $v(x)$ is zero for detailed balance dynamics (since $J_{\infty}$ is zero). Therefore, we check that

$$I_{\infty} = \mathcal{U} \quad \text{if} \quad f = 0 \qquad (8.23)$$

namely, the equilibrium canonical (Gibbs) distribution is the stationary probability distribution for detailed balance dynamics. It is an important fact that with detailed balance the subexponential prefactor $W$ is independent of the state $x$. It is a mere normalization constant, as can be checked by direct substitution into (8.11). This remains true for $f \neq 0$ for $\nabla \cdot f = 0$. In this case, the stationary probability remains Gibbsian even if $J_{\infty} \neq 0$. This kind of nondetailed balance dynamics is known to accelerate relaxation, i.e. the absolute value of the first nonzero eigenvalue $|\lambda_1|$ increases in the presence of a divergenceless $f \neq 0$ [125, 126].

This decomposition is useful to analyze the two main sets of solutions of (8.17), determined by the boundary conditions. We see that $\pi = 0$ is a solution giving the noiseless dynamics (8.12) which we have already seen as a relaxation to the stable fixed point $x_i^*$. Typical fluctuations (called optimal trajectories or instantons) are characterized by $\pi \neq 0$. One particular family of such trajectories is that happening on the manifold of Hamiltonian $H = 0$: from (8.20) we see that they correspond to $\pi(t) = \nabla I_{\infty}(x(t))$. Plugging that into the first of (8.17) we obtain

$$\dot{x} = F(x) + 2D(x) \cdot \nabla I_{\infty}(x) = D(x) \cdot \nabla I_{\infty}(x) + v(x). \qquad (8.24)$$

These dynamics take place in a 'reversed' rate function but with the same currents. They are trajectory, called instantons, that start arbitrarily close to the stable fixed point $x_i^*$ and reach $x_\nu$ in an arbitrarily long time $^{31}$. Notably, if $f = 0$ we find that deterministic dynamics and typical trajectories are one the time-reversed of the other, which is the content of detailed balance dynamics specialized to this set of trajectories.

If we use (8.20) in (8.19), and the insight that such instanton takes infinite time to escape a fixed point, we obtain the long-time transition probability

$$\lim_{t \to \infty} \mathrm{Prob}(x_i^* \to x) \asymp e^{-\frac{1}{\epsilon} \int_{x(0)}^{x(\infty)} dt \dot{x}(t) \cdot \nabla I_{\infty}(x(t))} = e^{-\frac{1}{\epsilon} [I_{\infty}(x) - I_{\infty}(x_i^*)]} \qquad (8.25)$$

Note that the trajectory could be continued from $x_\nu$ to the other stable fixed point $x_{j \neq i}^*$ without changing the result: the action along any relaxation trajectory is zero $^{32}$.

It can be proved [127] that for $\epsilon \to 0$ the escape time from a basin of attraction through the saddle $x_\nu$ is exponentially distributed with rate [128]

$$k_{\nu} \asymp \lim_{t \to \infty} \mathrm{Prob}(x_i^* \to x_{\nu}) \asymp e^{-\frac{1}{\epsilon} [I_{\infty}(x_{\nu}) - I_{\infty}(x_i^*)]} \qquad (8.26)$$

<footnote>$^{31}$This is because both $\nabla I_{\infty}(x)$ and $v(x)$ tend to zero close to the zeros of $F(I_{\infty}$ by definition, and $v$ as a consequence of (8.21)).</footnote>

<footnote>$^{32}$We could also start from any other point in the basin of attraction of $x_i^*$ without changing the result: For very long times we would first relax to $x_i^*$ (with a zero-action relaxation trajectory) and the follow the instanton.</footnote>
As a consequence, the mean first passage time to reach the unstable fixed point $x_ν$ from the basin of attraction of $x_i^*$ equals the inverse of the escape rate (8.26). This boils down to the Arrhenius law when $f = 0$; from eq. (8.23), it is

$$k_{\nu}^{eq} \approx e^{-\frac{1}{\epsilon}[\mathcal{U}(x_{\nu}) - \mathcal{U}(x_{i}^{*})]}. \qquad (8.27)$$

### 8.3.3 Thermodynamic bounds on the rate function and the escape rates

For $f$ arbitrary, we make use of the following relations

* The decomposition of the scaled entropy flow rate (in the reservoir) into the adiabatic and nonadiabatic component 3: $\dot{\sigma} = \dot{\sigma}_{ad} + \dot{\sigma}_{na}$ [129,130]. The entropy flow rate can be obtained by the log ratio of probabilities for the forward and backward trajectories conditioned to start on a given state, together with the drift field decomposition (8.21),

$$\sigma = \epsilon \ln \frac{\mathrm{Prob}[\{x(t)\}_0^\tau | x_0]}{\mathrm{Prob}[\{\tilde{x}(t)\}_0^\tau | x_\tau]} = \int_0^\tau dt \dot{x} \cdot D^{-1} \cdot F = \int_0^\tau dt \dot{x} \cdot D^{-1} \cdot (-D \cdot \nabla I_\infty + v). \qquad (8.28)$$

* The positivity of adiabatic entropy production rate along relaxation and instanton:

$$\dot{\sigma}_{ad} = \dot{x} \cdot D^{-1} \cdot v = (\pm D \nabla I_{\infty} + v) \cdot D^{-1} \cdot v = v \cdot D^{-1} \cdot v \geq 0. \qquad (8.29)$$

* the fact that the nonadiabatic entropy flow rate equals $\dot{\sigma}_{na} = -\dot{x} \cdot \nabla I_{\infty} = -\frac{d}{dt} I_{\infty}(x(t))$ along any trajectory.

Integrating the entropy production rate along the relaxation ($\sigma_{\nu \to i}$) and the instanton ($\sigma_{i \to \nu}$) and using the relation (8.26) between the rate function difference and the transition rate, we obtain

$$-\sigma_{\nu \to i} \leq \lim_{\epsilon \to 0} \epsilon \ln \kappa_{\nu} \leq \sigma_{i \to \nu}. \qquad (8.30)$$

With detailed balance, $\sigma_{\nu \to i} = -\sigma_{i \to \nu} = \mathcal{U}(x_\nu) - \mathcal{U}(x_i^*)$ because relaxation and instanton are the time-reversed of each other [131].

Moreover, the above equality holds even at first order in $f \to 0$, that is, around detailed balance [26]. Indeed, setting $I_\infty = \mathcal{U} + g$ (with $O(g) = O(f)$) in (8.20) and using the decomposition of the drift field we can obtain, neglecting terms $O(f)$ and higher

$$-\nabla \mathcal{U} \cdot D \cdot \nabla g = f \cdot \nabla \mathcal{U}. \qquad (8.31)$$

Here, the left-hand side is the time derivative of the correction $g$ along a relaxation trajectory with $f = 0$ and the right-hand side is minus the entropy flow rate along the same trajectory. Therefore, integrating over an infinite time from $x(0) = x_\nu$ to $x(\infty) = x_i^*$ we see that the correction $g$ to the energy is given by $\sigma_{\nu \to i} = -\sigma_{i \to \nu}$, and (8.30) holds as an equality in view of (8.26).

## 8.4 Jump processes in occupation-number space: macroscopic limit

The large deviations formalism can be applied to processes described by master equations that admit a deterministic limit. For example, in a state space coordinatized by the occupation number $n$ (individuals of a species in population dynamics, molecule number in chemical reactions, charge number in electronic circuits, photon number per frequency
<figure><img src="image_64.png" /><figcaption>Figure 8.5: Sketch of the typical dynamics on the three different timescales: short-time relaxation towards a local minimum of the stationary rate function, small (Gaussian) fluctuations around a stable fixed point, long-time transition to another attractor with rate $k_ν$.</figcaption></figure>

in scattering) there might exists a scale parameter (typically an inverse system size) $\epsilon$ such that $n\epsilon = x$ defines a finite continuous variable as $\epsilon \to 0$. If the transition matrix $R$ scales like $1/\epsilon$, we incur in a weak-noise limit similar to the one presented before. The Schlögl model given before is a particular example of this class of systems.

Here, we highlight only the main differences with respect to overdamped diffusion. From the WKB ansatz inserted into the Master equation we arrive at (8.15) with the Hamiltonian that is now the cumulant generating function of Poisson noise <fnref n="33" />

$$H(x, \pi) = \sum_{\rho} r_{\rho}(x)(e^{\Delta_{\rho}\pi} - 1) \qquad (8.32)$$

where $r_ρ(x)$ is the macroscopic (i.e. scaled by $1/\epsilon$) transition rate from state $x$ and $\Delta_ρ$ is the jump size of transition $ρ$.

* The resulting Hamiltonian equations are [132]

$$\begin{aligned} \dot{x} &= \sum_{\rho} \Delta_{\rho} r_{\rho}(x) e^{\Delta_{\rho}\pi} \\ \dot{\pi} &= -\sum_{\rho} \nabla r_{\rho}(x) (e^{\Delta_{\rho}\pi} - 1), \end{aligned} \qquad (8.33)$$

Only close to a stable fixed point (or time dependent-attractor) these equations can be linearized in $\pi$ and $x - x_i^*$ to obtain (the linear version of) (8.17). Clearly, this approximation is valid only for times $t \ll t_{esc} \approx k_ν^{-1}$ and correspond to van Kampen's system size expansion [133].

* The deterministic dynamics has a nonlinear gradient part plus the 'orthogonal' probability velocity

$$\dot{x} = -A(x) \cdot \nabla I_{\infty} + v \qquad (8.34)$$

<footnote n="33">³³Recall that a Markov jump process is characterized by a (conditional) exponential waiting time distribution, or equivalently, a (conditional) Poisson distribution of number of jumps per site.</footnote>
<figure><img src="image_65.png" /><figcaption>Figure 8.6: On very long times the dynamics can be coarse-grained into a discrete jump process from one attractor to the other</figcaption></figure>

with $A(x) = \int_0^1 d\theta (1-\theta) \partial_\pi^2 H(c, \pi)|_{\pi=\theta \nabla I_\infty}$ a nonnegative diffusion matrix depending on the rate function itself.

* The instanton is related to the relaxation by considering the adjoint dynamics which acts on single transition rates as $r_\rho(x) \mapsto r_{-\rho}(x)e^{-\Delta_\rho \nabla I_\infty(x)}$. This in general does not boil down to flipping the sign of the rate function and leaving $v$ untouched.

* The bounds (8.30) remain valid.

## 8.5 Markov jump process on the space of attractors

On very long times of the order $t \ge t_{esc}$ we can describe the dynamics as a jump process from one attractor to the other with transition rates $k_\nu$ and $p^{(i)}(t)$ the probability of being within the attractor $i$ (Fig. 8.6). Here $\nu$ labels all the transitions in both directions, counted as positive and negative, respectively. Because of the previous result (8.27), local detailed balance on the rates $k_\nu$ holds trivially when $f=0$ but also when $f$ is small, in the form [26]

$$ \frac{k_\nu}{k_{-\nu}} = e^{-[\mathcal{U}(o(-\nu)) - \mathcal{U}(o(\nu)) + \sigma_{o(\nu) \to o(-\nu)}]} \quad (8.35) $$

where $o(\pm\nu)$ is the origin of the transition trough $\pm\nu$ and $\sigma_{x_i^* \to x_j^*}$ is the entropy flow the most probable path between $x_i^*$ and $x_j^*$. As mentioned above, solutions of (8.20) only give the rate function within each basin of attraction $I_\infty^{(i)}(x)$, which is determined up to a constant $\alpha_i$. Intuitively, this is because to obtain $I^{(i)}$ we first took the limit $\epsilon \to 0$, before $t \to \infty$, thus the system is stuck within a basin of attraction. Such constants have to be fixed relative to each other. They give the relative weight of each attractor, i.e. $p_i$, obtained by solving the master equation over the space of attractors [124, 128]. The rate function is finally obtained by taking at each $x$ the maximum $I_\infty^{(i)}(x) + \alpha_i$ over $i$.
### 8.5.1 Example: the Schlögl model (cont.)

The system is bistable, so we can write the long time Markov jump dynamics in terms of a 2X2 stochastic matrix

$$ \begin{pmatrix} -k_1 & k_{-1} \\ k_1 & -k_{-1} \end{pmatrix} \qquad (8.36) $$

whose eigenvalues are 0 and $-(k_1 + k_{-1})$. One can check [80] that they correspond to the first two eigenvalues of the generator of the Schlögl model.

## 8.6 Conclusions

We conclude by listing the differences that emerged above between detailed balance and nondetailed balance dynamics.

* **Geometric vs dynamic view:** in detailed balance systems, it is possible to understand metastability in terms of the underlying free energy landscape. This geometrical perspective is lost in nondetailed balance systems where metastability may be entirely due to nonconservative forces. In this case its understanding is based on the dynamics: one looks at the behavior of the smallest system's eigenvalues.

* **Generator R:** in the detailed balance case, the generator of the stochastic process is symmetrizable. This means that all of its eigenvalues are real. In the nondetailed balance case, this is no longer true and one may have complex eigenvalues leading to relaxation with oscillations, limit cycles or more complex time-dependent attractors.

* **Rate function:** $I_\infty(x) = U(x)$ for detailed balance dynamics. When nonconservative forces are introduced, the relation is modified. Close to detailed balance, the correction to $U$ is given by the entropy dissipated by the nonconservative force on the typical trajectory.

* **Instanton:** the instanton is the time-reverse of the relaxation for detailed balance systems and their associated entropy productions are equal and opposite $\sigma_{\nu \to i} = -\sigma_{i \to \nu}$. When nonconservative forces are present, the instanton and the relaxation may follow different trajectories, and an inequality holds for the entropy productions $-\sigma_{\nu \to i} \le \sigma_{i \to \nu}$ (at first order in $f \to 0$ it becomes the above equality).

* **Coarse-grained dynamics:** on longer time-scales, $t \ge t_{esc}$, the dynamics resembles a jump process among mesostates. For detailed balance systems a local detailed balance condition emerges on the mesoscopic rates. The dissipation is entirely due to jumps among attractors, whereas in the nondetailed balance case, there is a contribution to the dissipation coming from the microscopic currents inside each basin of attraction.
# 9 Martingale Approach for First-Passage Problems

**Izaak Neri and Adarsh Raghu.**<sup>34</sup> We explore nonequilibrium thermodynamics for first-passage processes with martingale theory. We use the gambler's ruin problem as a unifying theme.

## 9.1 Introduction

First-passage problems describe processes that have a finite, random termination time. Examples are escape problems — such as, a particle that escapes from a bounded domain or a system that escapes from a metastable state (i.e., Kramers' escape problem) — and diffusion controlled reactions. For processes in equilibrium, the time it takes for a particle to escape a bounded domain decreases exponentially as a function of the energy barrier that separates it from the outside world, as described by the van't Hoff-Arrhenius law [2, 134, 135].

The focus of these notes is on first-passage problems in nonequilibrium scenarios, where the underlying stochastic process does not adhere to detailed balance; see Fig. 9.1 for illustrations of two examples of escape problems in nonequilibrium, statistical physics. In nonequilibrium scenarios, we are often interested in the escape problem of a macroscopic current out of a bounded interval. An instance of such a scenario is observed in molecular motors, like kinesin-1, which convert free energy into mechanical work while attached to a biofilament. Due to the polarity of biofilaments, the exerted force on the filament is biased towards one of its ends. Notably, as the filament has a finite length, the motor's work required to transport a cargo to one of the filament's ends is evaluated at a random time.

As thermodynamicists, our primary goal is to establish thermodynamic connections between quantities that characterize first-passage processes, such as splitting probabilities and mean first-passage times, and thermodynamic quantities, such as the entropy production rate. Our focus is on obtaining universal relations that are applicable to a wide range of stochastic processes far from thermal equilibrium, similar to the van't Hoff-Arrhenius law that holds for equilibrium processes [135].

To analyze stochastic processes that are far from thermal equilibrium, we employ stochastic thermodynamics [35, 136, 137]. However, instead of analyzing stochastic thermodynamics over ensembles of trajectories that terminate at a fixed time, we investigate ensembles of trajectories that terminate at a trajectory-dependent termination time (i.e., a random termination time). Recent works [138–140] have demonstrated that martingale theory significantly aids in studying thermodynamics at random termination times leading to the discovery of several universal relations in nonequilibrium thermodynamics; see also the earlier work Ref. [141] on the use of martingale in stochastic thermodynamics.

In particular, we review two results in this context. Firstly, we review the second law of thermodynamics at stopping times [140, 142, 143] that extends the classical second law of thermodynamics, which is valid for processes terminating at a fixed time, to processes that terminate at a random time. Secondly, we review a recently discovered relation that connects the splitting probability of a current with its mean first-passage time and the rate of dissipation in a nonequilibrium process [144–147]. This relation can be interpreted as a trade-off between uncertainty (expressed as the ruin probability $p_-$), speed (quantified by the mean first passage time $\langle \hat{T} \rangle$), and thermodynamic cost (quantified by the rate of dissipation $\sigma$). We explain these results by making an analogy with the gambler's

<footnote n="34"><sup>34</sup>IN was the lecturer and AR was the angel. Both authors wrote the chapter.</footnote>
ruin problem, as studied originally by Blaise Pascal, and then extending this problem to thermodynamic setups.

The organization of the Lecture Notes is as follows: Section 9.2 provides a review of the gambler's ruin problem, which serves as a basic illustration of a first-passage problem with two boundaries. In this section, we solve the gambler's ruin problem through the classical approach that relies on difference equations. Section 9.3 revises some essentials from martingale theory, and in Sec. 9.4, we utilise martingales to solve the gambler's ruin problem. In Sec. 9.5, we introduce various thermodynamic versions of the gambler's ruin problem and explain how to utilise martingale theory to solve some of them. These lecture notes conclude with a discussion of some of the questions that arose during the lectures at the workshop.

## 9.2 Gambler's ruin problem

The gambler's ruin problem was developed in 1656 by Blaise Pascal [148] and is often used in textbooks to introduce the first-passage problem and the mathematical tools that come with it, see e.g. Chapter XIV of Feller's introductory textbook book on probability theory [149]. Furthermore, from a physics perspective it can be seen as an elementary example of a first-passage problem in a nonequilibrium process.

First, we briefly discuss the history of this problem, which is interesting in itself, as it is part of the Scientific Revolution that took place in the 17th century. Second, we frame the problem and provide its standard solution using difference equations.

### 9.2.1 A brief history of the problem

Although probability theory is commonly used in our everyday reasoning, its mathematical foundations were established relatively late. Nevertheless, during the 17th century's Scientific Revolution the development of probability theory was bound to happen due to various circumstances.

In this period there was a surge of interest in gambling amongst the wealthy and noble classes in Europe leading to the establishment of gambling houses, such as, the Ridotto in Venice (1638). Several puzzles in probability theory had appeared in the literature and were discussed by prominent thinkers of the time, notably, the problem of points that was known from Luca Paccioli's work "Summa de Arithmetica, Geometria, Proportioni et Proportionalità" (Summary of Arithmetic, Geometry, Proportion and Proportionality), published in 1494. Blaise Pascal got interested in probability problems through his interactions with Antoine Gombaud, who was a French writer (also known as Chevalier de Méré), but also an amateur mathematician and a gambler. Against this backdrop, Pascal developed the gambler's ruin problem, two years after his famous correspondence with Pierre de Fermat on the problem of points in 1654 [150].

The first in print version of the gambler's ruin problem appeared soon after in Cristiaan Huygens' 1660 paper on probability entitled "Van Rekeningh in Spelen van Geluck" (On calculations for games on chance) [151], and a year later it was translated by his mathematics teacher Schooten into Latin as "De Ratiociniis in Ludo Aleae", and was published in Schooten's "Exercitationum Mathematicarum"; note that Huygens visited Paris for the first time in 1655 — a few months after he identified the rings of Saturn and discovered Saturn's moon Titan, — where he met with French scholars and became aware of their interest in probability theory. At the end of Huygens' paper he presents five problems to the reader, the last of which became better known as the gambler's ruin problem. Huygens' formulation goes as follows:
<figure><img src="image_66.png" /><figcaption>Figure 9.1: Examples of (nonequilibrium) processes with a finite, random termination time. *Gambler's ruin problem* [Left Panel]: a particle represents the wealth of gambler A and moves on a finite segment until it reaches either the “-”-end (indicating that gambler A has lost their stake) or the “+”-end (indicating that gambler B has lost their stake). The ratio of the odds of the two gamblers is denoted by $\omega_+/\omega_-$. *Molecular motor bound to a filament* [Right Panel]: A two-headed molecular motor, such as kinesin-1, attaches to a filament, after which it converts the free energy $\Delta\mu$, obtained from the hydrolysis of adenosine triphosphate (ATP) into adenosine diphosphate (ADP) and an inorganic phosphate (P), into mechanical work ($f_{\text{mech}}$), yielding a stochastic motion biased towards the filament's plus end (+). When the motor reaches one of its end points, then it detaches from the filament, and its enzymatic activity stops.</figcaption></figure>

"A and B have both 12 coins and play a game of chance with three dice under these conditions: if a throw shows 11 pips, then A must give a coin to B; however, if a throw shows 14 pips, then B must give a coin to A; the game is won by the player that has first all the coins.". Huygens also reveals the solution: "The odds of A with respect to B are 244140625 to 282429536481."

Huygens' publication generated significant interest in probability theory, notably from Jacob Bernoulli. In 1713, eight years after Bernoulli's death, his famous work "Ars Conjectand" was published. In this book, he reformulates the gamblers' ruin problem in a way that closely resembles its modern version, viz., gambler A has m coins, gambler B has n coins, and the odds of A winning to B or a to b. Jacob Bernoulli shows that the chance for A to win the game, and thus B to lose the game, is [152]

$$p_{+}^{A} = p_{-}^{B} = \frac{a^{n}(a^{m} - b^{m})}{a^{m+n} - b^{m+n}} \qquad (9.1)$$

Equation (9.1) implies that for $b > a$ and large values of $m$, the ruin probability $p_{-}^{B}$ decays exponentially as a function of $n$, i.e.,

$$p_{-}^{B} = \exp\left(n \ln \frac{a}{b}\right) \left(1 + O\left(\exp\left(m \ln \frac{a}{b}\right)\right)\right) \qquad (9.2)$$

where O represents the big O notation. Formulas similar to Eqs. (9.1) and (9.2) reoccur in the present Lecture Notes.

### 9.2.2 Problem statement and solution

For either gambler A or B, the gambler's ruin problem can be seen as a first-passage problem in which a particle hops on a segment until it reaches one of two endpoints, after which the process terminates, as depicted in the left Panel of Fig. 9.1. If the position of
the random walker, denoted by $\hat{X}_n \in [\ell] = \{0, 1, \dots, \ell\}$ with $n \in \mathbb{N}$ a discrete time index, represents the stakes of gambler A, then the absorbing boundary at the origin corresponds to the event where gambler A loses all of its stakes, while the other endpoint represents the event where gambler B loses all their stakes in the game. Note that the main quantities of interest, such as the profit of gambler A, are ensemble averages of $\hat{X}$ evaluated at the random termination time $\hat{T}$.

Let us formalise the above. Consider a Markov process $\hat{X}_n \in [\ell]$ governed by the transition matrix

$$P(\hat{X}_n = x | \hat{X}_{n-1} = y) = \omega_+ \delta_{x,y-1} + \omega_- \delta_{x,y+1} \quad (9.3)$$

for $y \in [\ell] \setminus \{0, \ell\}$, with

$$\omega_+ = 1 - \omega_-, \quad (9.4)$$

with absorbing boundary conditions

$$p(\hat{X}_n = x | \hat{X}_{n-1} = 0) = \delta_{x,0} \quad \text{and} \quad p(\hat{X}_n = x | \hat{X}_{n-1} = \ell) = \delta_{x,\ell}, \quad (9.5)$$

and with initial state $\hat{X}_0 = x_0 \in [\ell]$. When the particle reaches the boundaries $x = 0$ or $x = \ell$, then the process terminates, in the sense that the particle stops moving, and hence we have a process of finite, random termination time given by

$$\hat{T} = \min \{ n \in \mathbb{N} \cup \{0\} : \hat{X}_n = \ell \text{ or } \hat{X}_n = 0 \}. \quad (9.6)$$

Hence, each realisation $\hat{X}_0^{\hat{T}} = \{\hat{X}_0, \hat{X}_1, \dots, \hat{X}_{\hat{T}}\}$ has a different duration $\hat{T}$ and $\hat{X}_{\hat{T}} \in \{0, \ell\}$.

Classically, the gambler's ruin problem is solved with difference equations [149]. In what follows, we illustrate the solution for $\omega_+ \neq 1/2$, and we refer to Ref. [149] for more details, including the derivation for $\omega_+ = 1/2$. Specifically, we determine the main observables of interest, viz., the probability of ruin $p_x^-$, the average profit $\langle \hat{X}_{\hat{T}} \rangle$, the average termination time $\langle \hat{T} \rangle$, and the generating functions of the termination time $\hat{T}$ at the negative boundary, $g_x^-$, and at the positive boundary, $g_x^+$.

**Probability of ruin.** The probability of ruin,

$$p_x^- := P(\hat{X}_{\hat{T}} = 0 | \hat{X}_0 = x) \quad (9.7)$$

solves the difference equation

$$p_x^- = \omega_+ p_{x+1}^- + \omega_- p_{x-1}^- \quad (9.8)$$

with boundary conditions $p_\ell^- = 0$ and $p_0^- = 1$. The solution to these equations are unique, as follows from the Maximum Principle for harmonic functions, see Ref. [153] or Theorem 2.1 in [154].

First we consider the Eq. (9.8) without boundary conditions. In this case, $p_x^-$ should be seen as a function of $x$ rather than the splitting probability in the gambler's ruin problem. For functions of the form $p_x^- = \alpha^x$, where $\alpha \in \mathbb{R}^+$, we obtain from Eq. (9.8) that $\alpha = \omega_+\alpha^2 + \omega_-$. For $\omega_+ \neq 1/2$, this quadratic equation contains two solutions, namely, $\alpha = 1$ and $\alpha = \omega_- / \omega_+$. Using the linearity of the difference Eq. (9.8), we obtain

$$p_x^- = c_1 + c_2 \left( \frac{\omega_-}{\omega_+} \right)^x, \quad (9.9)$$
<figure><img src="image_67.png" /><figcaption>Figure 9.2: Top: Ruin probability $p_x^-$ as a function of $x/\ell$ for parameters $\ell = 10$, $\omega_- = 0.55$. Bottom: The expected profit $\langle \hat{X}_{\hat{T}} | \hat{X}_0 = x \rangle$ in the gambler's ruin problem for the same parameters $\ell = 10$, $\omega_- = 0.55$.</figcaption></figure>

where $c_1$ and $c_2$ are arbitrary constants fixed by the boundary conditions. Second, we impose the boundary conditions on $p_x^-$ to obtain the splitting probability

$$p_x^- = \frac{(\omega_- / \omega_+)^\ell - (\omega_- / \omega_+)^x}{(\omega_- / \omega_+)^\ell - 1} \quad (9.10)$$

for the gambler's ruin problem. We have plotted $p_x^-$ as a function of $x/\ell$ in the top Panel of Fig. 9.2. Notice that Eq. (9.10) is equivalent to Bernoulli's formula Eq. (9.1) when identifying $n = x$, $m = \ell - x$, $b = \omega_+$, and $a = \omega_-$.

**Profit at the termination time.** The average profit $\langle \hat{X}_{\hat{T}} \rangle - \langle \hat{X}_0 \rangle$ at the termination time is given by

$$\langle \hat{X}_{\hat{T}} \rangle - \langle \hat{X}_0 \rangle = p_x^+(\ell - x) - xp_x^- = \frac{\ell - x - (\omega_- / \omega_+)^x (\ell - x (\omega_- / \omega_+)^{\ell - x})}{1 - (\omega_- / \omega_+)^\ell} \quad (9.11)$$

In the bottom panel of Fig. 9.2, we plot $\langle \hat{X}_{\hat{T}} \rangle$ as a function of $x/\ell$. Note that the profit is nonnegative for $\omega_+ \ge \omega_-$, a fact that we demonstrate in Sec. 9.4 with martingale theory.

**Mean first-passage time.** The mean duration of the process,

$$t_x := \langle \hat{T} | \hat{X}_0 = x \rangle, \quad (9.12)$$

solves

$$t_x = \omega_+ t_{x+1} + \omega_- t_{x-1} + 1 \quad (9.13)$$

with boundary conditions $t_0 = 0$ and $t_\ell = 0$. Solving these difference equations for $\omega_+ \ne 1/2$ towards $t_x$, we obtain [149]

$$t_x = -\frac{x}{\omega_+ - \omega_-} + \frac{\ell}{\omega_+ - \omega_-} \frac{1 - (\omega_- / \omega_+)^x}{1 - (\omega_- / \omega_+)^\ell} \quad (9.14)$$
**Generating function.** The fluctuations of the termination time at the negative boundary are determined by the generating function

$$g_x^-(s) := \sum_{n=0}^{\infty} \mathbb{P}(\hat{T} = n, \hat{X}_{\hat{T}} = 0 | \hat{X}_0 = x) s^n, \qquad (9.15)$$

which solves the difference equation

$$g_x^-(s) = \omega_+ s g_{x+1}^-(s) + \omega_- s g_{x-1}^-(s) \qquad (9.16)$$

with boundary conditions $g_0^-(s) = 1$ and $g_\ell^-(s) = 0$. Solving towards $g_x(s)$, we obtain for $\omega_+ \neq 1/2$ [149]

$$g_x^-(s) = \left( \frac{\omega_-}{\omega_+} \right)^x \frac{\lambda_+^{\ell-x}(s) - \lambda_-^{\ell-x}(s)}{\lambda_+^\ell(s) - \lambda_-^\ell(s)} \qquad (9.17)$$

where

$$\lambda_{\pm}(s) := \frac{1 \pm \sqrt{1 - 4\omega_{-}\omega_{+}s^{2}}}{2\omega_{+}s}. \qquad (9.18)$$

Analogously, for

$$g_x^+(s) := \sum_{n=0}^{\infty} \mathbb{P}(\hat{T} = n, \hat{X}_{\hat{T}} = \ell | \hat{X}_0 = x) s^n, \qquad (9.19)$$

we find

$$g_x^+(s) = \frac{\lambda_+^x(s) - \lambda_-^x(s)}{\lambda_+^\ell(s) - \lambda_-^\ell(s)}. \qquad (9.20)$$

Note that setting $s = 0$, we recover $g_x^+(0) = p_x^+$ and $g_x^-(0) = p_x^-$.

## 9.3 Martingales: Definition, Represenations, and Properties

As noted by Jean-André Ville in his PhD thesis [155], martingales provide an alternative framework to solve the gambler's ruin problem. We provide here a brief overview of martingale theory, and in the next section we discuss how to use martingales to solve the gambler's ruin problem.

### 9.3.1 Definition

Let $\hat{X}_m \in \mathcal{X}$ be a stochastic process taking values in a finite set $\mathcal{X}$, and let $m \in \mathbb{N}$ be a discrete time index. We denote a sequence of (random) variables by

$$\hat{X}_0^n := (\hat{X}_0, \hat{X}_1, \dots, \hat{X}_n) \qquad (9.21)$$

Let $\mathbb{P}$ be the probability measure describing the statistics of $X$, and we denote its path probability by

$$p(x_0^n) := \mathbb{P}(\hat{X}_0^n = x_0^n). \qquad (9.22)$$

For simplicity, we assume in these lecture notes that both the phase space $\mathcal{X}$ and time $n$ are discrete, but all the concepts can be extended into a continuous setting, see Ref. [143].

A martingale $\hat{M}_n(\hat{X}_0^n)$, relative to $\hat{X}$ and $\mathbb{P}$, is a real-valued process defined on $\hat{X}_0^n$ for which it holds that [156]

$$\langle |\hat{M}_m| \rangle < \infty \qquad (9.23)$$

and

$$\langle \hat{M}_n | \hat{X}_0^m \rangle = \hat{M}_m \qquad (9.24)$$
for all $m \in [n] := \{0, 1, \dots, n\}$. Here the averages $\langle \cdot \rangle$ are taken with respect to the probability measure $\mathbb{P}$. Using the definition of conditional probabilities, we get

$$ \langle \hat{M}_n | \hat{X}_0^m \rangle = \sum_{x_{m+1}^n \in \mathcal{X}^{m-n}} \frac{p(x_0^n)}{p(x_0^m)} M_n(x_0^n). \qquad (9.25) $$

We require the first condition, Eq. (9.23), as it guarantees that the conditional expectation $\langle \hat{M}_n | \hat{X}_0^m \rangle$ exists, see e.g., Theorem 10.1.1 in [157]. The second condition (9.24) is the defining property of the martingale. If we replace the condition (9.24) by the inequality

$$ \langle \hat{M}_n | \hat{X}_0^m \rangle \geq \hat{M}_m \qquad (9.26) $$

then we speak of a submartingale process.

We give a brief historical overview on martingales. The martingale condition Eq. (9.24), appears in Paul Lévy's book [158] published in 1937 on the sums of random variables as a technical condition required to extend the classical central limit theorem to sums of dependent random variables. However, Paul Lévy does not define martingales and does not study their properties. It is, around the same time, Jean-André Ville who introduces stochastic processes satisfying the condition Eq. (9.24) in his PhD thesis [155], and he coined them martingales. In chapter V of his thesis, Ville proves several properties of martingales and uses them to study, amongst others, the gambler's ruin problem. Developments of martingale theory accelerated after the second world war, mainly after the publication of Joseph Doob's book on stochastic processes [159]; Doob, aware of Jean-Ville's work, borrows its name. Several important results in martingale theory are named after Doob, such as, Doob's regularity theorem [160], which plays an important role in defining stochastic processes on spaces of right-continuous trajectories, the martingale convergence theorems [159], and Doob's stopping theorems [159], on which we will rely in this short review.

Nowadays, martingales are models for stocks prices under the efficient market hypothesis [161], and in physics martingales describe important functionals in nonequilibrium, thermodynamics, such as, the exponentiated negative entropy production [138, 140–142] and the exponential of the housekeeping heat [162].

### 9.3.2 Martingale representations

We discuss two important examples of processes that are martinagles.

*   **Nonnegative martingales:** Let $\hat{M}_n$ be a nonnegative martingale relative to $\hat{X}$ and $\mathbb{P}$, and additionally we assume that $\langle \hat{M}_0 \rangle = 1$. Then $\hat{M}_n$ is the ratio of two path probabilities. Indeed, if we define

$$ \mathbb{Q}(\hat{X}_0^n = x_0^n) := \hat{M}_n(x_0^n) \mathbb{P}(\hat{X}_0^n = x_0^n), \qquad (9.27) $$

then $\mathbb{Q}(\hat{X}_0^n = x_0^n) \geq 0$ and

$$ \sum_{x_0^n \in \mathcal{X}^n} \mathbb{Q}(\hat{X}_0^n = x_0^n) = \langle \hat{M}_n \rangle = 1. \qquad (9.28) $$

The converse is also true. If $q(x_0^n) := \mathbb{Q}(\hat{X}_0^n = x_0^n)$ is a path probability that is absolutely continuous with respect to $p(x_0^n) := \mathbb{P}(\hat{X}_0^n = x_0^n)$, i.e., $p(x_0^n) = 0 \Rightarrow q(x_0^n) = 0$, then

$$ \hat{M}_n := q(\hat{X}_0^n)/p(\hat{X}_0^n) \qquad (9.29) $$
exists and is a martingale. Indeed, $\langle |\hat{M}_n| \rangle = \langle \hat{M}_n \rangle = 1$, and for all $0 \le m \le n$ it
holds that

$$
\begin{align}
\langle \hat{M}_n | \hat{X}_0, \hat{X}_1, \dots, \hat{X}_m \rangle
&= \sum_{x_{m+1} \in \mathcal{X}^{m-n}} \mathbb{P}(\hat{X}_{m+1}^n = x_{m+1}^n | \hat{X}_0^m) \frac{q(\hat{X}_0^m, x_{m+1}^n)}{p(\hat{X}_0^m, x_{m+1}^n)} \nonumber \\
&= \sum_{x_{m+1}^n \in \mathcal{X}^{m-n}} \frac{p(\hat{X}_0^m, x_{m+1}^m)}{p(\hat{X}_0^m)} \frac{q(\hat{X}_0^m, x_{m+1}^n)}{p(\hat{X}_0^m, x_{m+1}^n)} \nonumber \\
&= \frac{q(\hat{X}_0^m)}{p(\hat{X}_0^m)} = M_m, \tag{9.30}
\end{align}
$$

where we have used $p(\hat{X}_0^m, x_{m+1}^n)$ to denote the path probability of the sequence $(\hat{X}_0, \hat{X}_1, \dots, \hat{X}_m, x_{m+1}, \dots, x_n)$. Notice that in the last line we have used the marginalisation properties

$$
\sum_{x_{m+1}^n \in \mathcal{X}^{m-n}} q(\hat{X}_0^m, x_{m+1}^m) = q(\hat{X}_0^m) \quad (9.31)
$$

of path probabilities.

• Harmonic functions: Let $\hat{X} \in \mathcal{X}$ be a Markov process with transition matrix

$$
\mathbf{R}_{xy} := \mathbb{P} \left( \hat{X}_{n+1} = x | \hat{X}_n = y \right), \quad (9.32)
$$

and let $\mathbf{h}(x)$ be a left eigenvector of $\mathbf{R}$ associated with the unit eigenvalue, i.e.,

$$
\mathbf{h}(y) = \sum_{x \in \mathcal{X}} \mathbf{h}(x) \mathbf{R}_{xy}. \qquad (9.33)
$$

It then holds that

$$
\hat{h}(\hat{X}_n) := \mathbf{h}(\hat{X}_n) \qquad (9.34)
$$

is martingale. Indeed,

$$
\langle \hat{h}(\hat{X}_n) | \hat{X}_{n-1} \rangle = \sum_{x \in \mathcal{X}} \mathbf{h}(x) \mathbf{R}_{x\hat{X}_{n-1}} = \hat{h}(\hat{X}_{n-1}). \quad (9.35)
$$

Analogously,

$$
\langle \hat{h}(\hat{X}_n) | \hat{X}_{n-2} \rangle = \sum_{x_1, x_2 \in \mathcal{X}} \mathbf{h}(x_2) \mathbf{R}_{x_2 x_1} \mathbf{R}_{x_1, \hat{X}_{n-1}} = \hat{h}(\hat{X}_{n-2}) \quad (9.36)
$$

and so forth.

The functions $h$ are also called *harmonic functions*. Note that these should not be confused with the right eigenvectors

$$
\mathbf{p}_{\infty}(x) = \sum_{y \in \mathcal{X}} \mathbf{R}_{xy} \mathbf{p}_{\infty}(y), \qquad (9.37)
$$

which represent stationary distributions. For ergodic processes, there exists exactly one stationary distribution, and correspondingly, also exactly one harmonic function, viz., $\mathbf{h}(x) = c$ for all $x \in \mathcal{X}$, where $c$ is an arbitrary constant independent of
x. Hence, to have a nontrivial harmonic functions we need to make the process nonergodic.

A possibility to make the process nonergodic is to introduce absorbing states. Specifically, if we set

$$h(x) = 1 \quad \text{if} \quad x \in \mathcal{X}_+ \qquad (9.38)$$

and

$$h(x) = 0 \quad \text{if} \quad x \in \mathcal{X}_- \qquad (9.39)$$

then $h(x)$ is the probability to reach $\mathcal{X}_+$ before reaching $\mathcal{X}_-$ when the initial state $X_0 = x$. Hence, splitting probabilities are examples of harmonic functions, and the difference Eq. (9.8) in Sec. 9.2.2 is equivalent to the Eq. (9.35).

### 9.3.3 Doob's optional stopping theorem

We end the first part of this lectures with Doob's optional stopping theorem, which is the main theorem we will be using.

Doob's optional stopping theorem states the following, see Theorem 4.1.1. in Ref. [163].

**Theorem 1** Let $\hat{M}_n$ be a martingale with respect to $\hat{X}$ and P, and let $\hat{T}$ be a stopping time with respect to $\hat{X}$ and P. If one of the following two conditions holds, viz.,

* $\hat{T} \le n$ for some $n$;

* or, $\hat{T} < \infty$ and there exists some $c \in \mathbb{R}^+$ such that $|\hat{M}_n| \le c$ for all $n \le \hat{T}$;

then

$$\langle \hat{M}_{\hat{T}} \rangle = \langle \hat{M}_0 \rangle. \qquad (9.40)$$

If the abovementioned conditions do not apply, then Doob's optional stopping theorem is in general not valid. For example, for an unbiased random walk $X_n$ with $\omega_+ = \omega_- = 1/2$, which is a martingale, and for $\hat{T} = \min\{n \ge 0 : \hat{X}_n = x_+\}$, it holds that $\langle \hat{X}_{\hat{T}} \rangle = x_+ \ne \langle \hat{X}_0 \rangle$, and thus Doob's optional stopping theorem does not apply. One way to bound $\hat{M}_n$ for $n \le \hat{T}$ is by defining $\hat{T}$ as the escape time from a bounded set.

For submartingales $M_n$, it holds analogously that $\langle \hat{M}_{\hat{T}} \rangle \ge \langle \hat{M}_0 \rangle$.

## 9.4 Martingale solution to the gambler's ruin problem

We now solve the gambler's ruin problem with martingales.

### 9.4.1 The average profit

To illustrate the use of martingale theory, let us first consider the average profit $\langle \hat{X}_{\hat{T}} \rangle$, which is nonnegative as we have seen from explicitly deriving the formula Eq. (9.11). Using martingale theory, this result is seen as a direct consequence of the submartingale property of $\hat{X}_n$. Indeed, since for $\omega_+ \ge \omega_-$,

$$\langle \hat{X}_n | \hat{X}_0^m \rangle \ge \hat{X}_m, \qquad (9.41)$$

the process $\hat{X}_n$ is a submartingale, and hence Doob's optional stopping theorem for submartingales applies, yielding

$$\langle \hat{X}_{\hat{T}} \rangle \ge \langle \hat{X}_0 \rangle. \qquad (9.42)$$
Thus, the profit $\hat{X}_{\hat{T}} - \hat{X}_0$ is on average nonnegative. Notably, this result holds for any stopping strategy $\hat{T}$ as long as $\hat{X}_n$ is bounded for all times $n < \hat{T}$ (i.e., the total amount of money is finite). Conversely, when $\omega_+ \le \omega_-$, then

$$ \langle \hat{X}_{\hat{T}} \rangle \le \langle \hat{X}_0 \rangle, \qquad (9.43) $$

and the gambler is losing money, no matter which betting strategy, determined by $\hat{T}$, is used.

## 9.4.2 Family of martingale processes

To obtain the splitting probabilities and the generating function of $\hat{T}$, we introduce a family of martingales.

First, we define the process $\hat{X}_n$ as in Sec. 9.2.2 with the transition matrix Eq. (9.3) but without absorbing boundary conditions, i.e., $\hat{X} \in \mathbb{Z}$. The process

$$ \hat{Z}_n := \exp\left(z\hat{X}_n + nf(z)\right), \qquad (9.44) $$

with

$$ f(z) := -\ln[\omega_+ \exp(z) + \omega_- \exp(-z)] \qquad (9.45) $$

is then a martingale with respect to $\hat{X}$ and $\mathbb{P}$ for all values $z \in \mathbb{R}$.

Indeed, condition (9.24) follows from

$$ \begin{aligned} \langle \hat{Z}_n | \hat{X}_0, \hat{X}_1, \dots, \hat{X}_m \rangle &= \sum_{x \in \mathbb{Z}} \exp(zx + nf(z)) \mathbb{P}(\hat{X}_n = x | \hat{X}_0^m) \\ &= \sum_{n_+ = 0}^{n-m} \exp\left(z(\hat{X}_m + n_+ - (n-m-n_+)) + nf(z)\right) \binom{n-m}{n_+} \omega_+^{n_+} \omega_-^{n-m-n_+} \\ &= \exp\left(z\hat{X}_m + mf(z)\right) \exp\left\{\left((n-m)f(z)\right)\right\} \sum_{n_+ = 0}^{n-m} \binom{n-m}{n_+} (e^z \omega_+)^{n_+} (e^{-z} \omega_-)^{n-m-n_+} \\ &= \hat{Z}_m \left( \frac{\omega_+ e^z + \omega_- e^{-z}}{\exp(-f(z))} \right)^{n-m} = \hat{Z}_m, \end{aligned} \qquad (9.46) $$

where $n_+$ represents the number of forward steps taken by $\hat{X}$ when moving from $\hat{X}_m$ towards $\hat{X}_n$.

## 9.4.3 Splitting probabilities

We use the martingales $\hat{Z}_n$ to determine the splitting probabilities in the gambler's ruin problem. Using Doob's optional stopping theorem (Theorem 1) for $\hat{Z}_n$ and the stopping time $\hat{T} = \min\{n \ge 0 : \hat{X}_n \notin (0, \ell)\}$, we get

$$ p_x^+ \langle \exp(z\ell + \hat{T}f(z)) | \hat{X}_{\hat{T}} = \ell \rangle + p_x^- \langle \exp(\hat{T}f(z)) | \hat{X}_{\hat{T}} = 0 \rangle = \exp(zx_0). \qquad (9.47) $$

Setting

$$ z = \ln\left(\frac{1-q}{q}\right), \qquad (9.48) $$

the nontrivial root of $f(z)$ (i.e., $f(z) = 0$ and $z \ne 0$), Eq. (9.47) can be solved together with

$$ p_x^+ + p_x^- = 1, \qquad (9.49) $$

yielding Eq. (9.10) for the gambler's ruin probability.
### 9.4.4 Generating functions

To obtain the generating functions $g_x^+(s)$ and $g_x^-(s)$, we consider Eq. (9.47) at the two roots $z = z_\pm(s)$ that solve the equation

$$f(z_\pm) = \ln(s). \qquad (9.50)$$

This yields the two equations

$$\exp(z_\pm \ell) g_x^+(s) + g_x^-(s) = \exp(z_\pm x) \qquad (9.51)$$

where $g_x^\pm$ are the generating functions of $\hat{T}$ defined in Sec. 9.2.2. Solving Eq. (9.50) towards $z_\pm$ gives

$$z_\pm = \ln \lambda_\pm(s) \qquad (9.52)$$

where $\lambda_\pm(s)$ are as defined in Eq. (9.18). Solving the Eqs. (9.51) towards $g_x^\pm$ yields Eqs. (9.17) and (9.20).

## 9.5 Thermodynamic versions of the gambler's ruin problem

We define thermodynamic versions of the gambler's ruin problem, which account for the attributes of nonequilibrium processes. As nonequilibrium processes are distinguished by currents $\hat{J}_n$ with a non-zero average value, $\langle \hat{J}_n \rangle > 0$, we focus on first-passage problems of currents $\hat{J}_n$ in Markov processes $\hat{X}_n$. Note that currents $\hat{J}_n$ are functionals defined on the trajectories $\hat{X}_0^n$, and hence they are in general nonMarkovian processes, which complicates the analysis of the corresponding first-passage problem. Nevertheless, using martingale theory we will be able to solve several first-passage problems of currents in nonequilibrium processes $X$.

### 9.5.1 General setup

Let $\hat{X}_n \in \mathcal{X}$, with $n \in \mathbb{N}$ and $\mathcal{X}$ a discrete set, be a stochastic process whose statistics described by $\mathbb{P}$ are Markovian with transition matrix $\mathbf{R}$, see Eq. (9.32). We assume that $\hat{X}$ is irreducible and that there exists a unique stationary, probability mass function $p_\infty(x)$ that solves Eq. (9.37) for all $x \in \mathcal{X}$. In other words, the process is ergodic, see Theorem 4.1 in Ref. [154]. In addition, we assume that $\mathbf{R}_{xy} > 0$ whenever that $\mathbf{R}_{yx} > 0$ so that the time-reversed Markov process exists.

An empirical integrated current $\hat{J}$ is a stochastic process of the form

$$\hat{J}_n := \sum_{x,y \in \mathcal{E}} c_{xy} \hat{J}_n(y|x) \qquad (9.53)$$

where the empirical edge current

$$\hat{J}_n(y|x) := \hat{N}_n(y|x) - \hat{N}_n(x|y) \qquad (9.54)$$

is the difference between the number of jumps $\hat{N}_n(y|x)$ in $\hat{X}_0^n$ from $x$ to $y$ and the number of jumps $\hat{N}_n(x|y)$ from $y$ to $x$ in the same trajectory, $c_{xy} \in \mathbb{R}$ are constant coefficients, and $\mathcal{E}$ is the set of pairs $(x, y)$ for which $\mathbf{R}_{xy} \neq 0$. In a stationary process, we denote the average current by

$$j(y|x) = \mathbf{R}_{yx}\mathbf{p}_\infty(x) - \mathbf{R}_{xy}\mathbf{p}_\infty(y). \qquad (9.55)$$

The gambler's ruin problem in $\hat{J}$ reads

$$\hat{T}_J := \min \left\{ n \ge 0 : \hat{J}_n \notin (-J_-, J_+) \right\}. \qquad (9.56)$$
where the thresholds $J_-$, $J_+ \ge 0$, and we assume, without loss of generality, that $\langle \hat{J}_n \rangle > 0$. The probability of ruin is then defined by

$$p_J^- := \mathbb{P} \left( \hat{J}_{\hat{T}_J} \le -J_- \right). \qquad (9.57)$$

The thermodynamic cost of a process can be quantified in terms of the stochastic entropy production [35, 136, 137]

$$\hat{S}_n := \hat{S}^{\text{sys}}(\hat{X}_n) - \hat{S}^{\text{sys}}(\hat{X}_0) + \hat{S}_n^{\text{env}} \qquad (9.58)$$

with

$$\hat{S}^{\text{sys}}(\hat{X}_n) = -\ln \mathbf{p}_{\infty}(\hat{X}_n) \qquad (9.59)$$

the Shannon entropy capturing the information theoretic content of the state $\hat{X}_n$, and

$$\hat{S}_n^{\text{env}} := \frac{1}{2} \sum_{(x,y) \in \mathcal{E}} \ln \frac{\mathbf{R}_{xy}}{\mathbf{R}_{yx}} \hat{J}_n(y|x) \qquad (9.60)$$

is the environment entropy change. The environment entropy is based on the principle of local detailed balance [25], which can be seen as a stochastic version of the local equilibrium concept [164]. At stationarity, the average rate of dissipation is

$$\sigma := \langle \hat{S}_n \rangle / n = \frac{1}{2} \sum_{(x,y) \in \mathcal{E}} \mathbf{R}_{xy} \mathbf{p}_{\infty}(y) \ln \frac{\mathbf{R}_{xy}}{\mathbf{R}_{yx}}. \qquad (9.61)$$

In what follows, we discuss the gambler's ruin problem for three examples of currents $\hat{J}$.

### 9.5.2 Gambler's ruin problem for entropy production

The first case we consider is for $\hat{J} = \hat{S}$, the stochastic entropy production. The stopping problem for entropy production reads

$$\hat{T}_S = \min \left\{ n \ge 0 : \hat{S}_n \notin (-S_-, S_+) \right\}. \qquad (9.62)$$

This case is relevant for two reasons: (i) in the uni-cyclic case, the macroscopic current is asymptotically proportional to the entropy production; (ii) using the stopping problem $\hat{T}_{\hat{S}}$ will establish relationships that resemble fluctuation relations for entropy production, albeit with regard to characteristics of entropy production at random times.

The calculation of the ruin probability $p_S^-$ and the average thermodynamic cost $\langle \hat{S}_{\hat{T}_S} \rangle$ is facilitated by the observation that the process $\exp(-\hat{S}_n)$ is a martingale [138, 141, 143]. Indeed,

$$e^{-\hat{S}_n} = \frac{\tilde{p}(\hat{X}_0^n)}{p(\hat{X}_0^n)} \qquad (9.63)$$

is the ratio of the two path probabilities [35, 136, 137]

$$p(\hat{X}_0^n) = \mathbf{p}_{\infty}(\hat{X}_0) \prod_{i=1}^{n} \mathbf{R}_{\hat{X}_i \hat{X}_{i-1}} \qquad (9.64)$$

and

$$\tilde{p}(\hat{X}_0^n) = \mathbf{p}_{\infty}(\hat{X}_0) \prod_{i=1}^{n} \tilde{\mathbf{R}}_{\hat{X}_i \hat{X}_{i-1}} \qquad (9.65)$$
with

$$ \tilde{\mathbf{R}}_{xy} = \mathbf{R}_{yx} \frac{\mathbf{P}_{\infty}(x)}{\mathbf{P}_{\infty}(y)}. \qquad (9.66) $$

Hence, according to Eq. (9.30), $e^{-\hat{S}(t)}$ is a martingale. Note that

$$ \tilde{p}(\Theta_n(\hat{X}_0^n)) = p(\hat{X}_0^n), \qquad (9.67) $$

where $\Theta_n(\hat{X}_0^n) = (\hat{X}_n, \hat{X}_{n-1}, \dots, \hat{X}_0)$ is the time-reversed trajectory, and hence $\tilde{p}$ describes the statistics of an observer that goes backwards in time.

Using Doob's optional stopping theorem, see Sec. 9.3.3, we find that

$$ \langle e^{-\hat{S}_{\hat{T}_S}} \rangle = 1, \qquad (9.68) $$

which is a version of the Integral Fluctuation Theorem at Stopping Times [140]. Further, using that

$$ \mathbb{P}(\hat{T}_S < \infty) = p_S^+ + p_S^- = 1, \qquad (9.69) $$

we obtain

$$ p_S^+ \langle e^{-\hat{S}_{\hat{T}_S}} \rangle_+ + p_S^- \langle e^{-\hat{S}_{\hat{T}_S}} \rangle_- = 1. \qquad (9.70) $$

If $\hat{S}_{\hat{T}_S} \in \{-S_-, S_+\}$, which includes the continuum limit, then Eqs. (9.69) and (9.70) yield

$$ p_S^- = \frac{e^{S_+} - 1}{e^{S_+ + S_-} - 1} \qquad (9.71) $$

for the ruin probability. Notice that Eq. (9.71) can be identified with Bernoulli's formula Eq. (9.1) when we set $S_+ = m \ln(b/a)$ and $S_- = n \ln(b/a)$, with $a, b \in (0, 1)$ arbitrary numbers. Analogously, setting $S_+ = (\ell - x) \ln(\omega_+/\omega_-)$ and $S_- = x \ln(\omega_+/\omega_-)$ in Eq. (9.71), we obtain the Eq. (9.10) for the ruin probability in the gambler's ruin problem. For processes that are non continuous, we get the inequality [140]

$$ p_S^- \le \frac{1}{e^{S_-} - e^{-S_+}} \qquad (9.72) $$

for the probability of ruin.

Interestingly, the results Eqs. (9.71) and (9.72) for $p_-$ allow us to precisely quantify the negative fluctuations of the entropy production in terms of its infimum value

$$ \hat{S}_{\inf} := \inf_{n \ge 0} \hat{S}_n. \qquad (9.73) $$

Indeed, observe that

$$ \lim_{S_+ \to \infty} p_S^- = \mathbb{P}(\hat{S}_{\inf} \le -S_-). \qquad (9.74) $$

Therefore, Eq. (9.72) yields in the limit of $S^+ \gg 1$ the inequality

$$ \mathbb{P}(\hat{S}_{\inf} \le -s) \le \exp(-s) \qquad (9.75) $$

such that on average

$$ \langle \hat{S}_{\inf} \rangle \ge -1, \qquad (9.76) $$

with equalities attained in the continuum limit. The inequality (9.75) readily implies the bound

$$ \mathbb{P}(\hat{S}_n \le -s) \le \exp(-s) \qquad (9.77) $$
<figure><img src="image_68.png" /><figcaption>Figure 9.3: Two examples of stochastic processes that model the dynamics of two-headed molecular motors, such as, kinesin-1. The motor states (D:T), (T:φ), and so forth, denote the chemical states of the rear and front motor head, viz., T stands for ATP bound, D for ADP bound, D' for ADP-P bound, and φ for nucleotide free. Left Panel: Model taken from Ref. [165] with three cycle currents, corresponding to forward motion (j<sub>f</sub>), backward motion (j<sub>b</sub>), and a current with no motion yielding a futile cycle [165]. Notably, the motor position along the biofilament is given by the edge current J<sub>n</sub><sup>2→5</sup>. Right Panel: Molecular motor model of Ref. [166]. In this case, there are two cycles corresponding with forward motion (with currents j<sub>f,1</sub> and j<sub>f,2</sub>, respectively), a cycle yielding backward motion (with current j<sub>b</sub>), and several futile cycles.</figcaption></figure>

that follows from the integral fluctuation relation $\langle \exp(-\hat{S}_n) \rangle = 1$ [35]. This illustrates how martingale methods can be used to obtain bounds on negative fluctuations of entropy production that are tighter than those obtained from fixed-time fluctuation relations in stochastic thermodynamics.

Furthermore, if $\hat{S}_{\hat{T}_S} \in \{-S_-, S_+\}$, then we obtain

$$ \langle \hat{S}_{\hat{T}_S} \rangle = p_S^+ S_+ - p_S^- S_- = \frac{S_+(e^{S_-} - 1) - S_-(1 - e^{-S_+})}{e^{S_-} - e^{-S_+}} \quad (9.78) $$

for the average thermodynamic cost at the termination time. Setting $S_+ = (\ell - x) \ln(\omega_+/\omega_-)$ and $S_- = x \ln(\omega_+/\omega_-)$ in Eq. (9.78), yields Eq. (9.11) for the average profit at the termination time. More generally, we have the inequality

$$ \langle \hat{S}(\hat{T}_S) \rangle \geq S_+ - \frac{S_+ + S_-}{e^{S_-} - e^{-S_+}} \geq 0, \quad (9.79) $$

which is a specific instance of the second law of thermodynamics at stopping times [140, 142]. This law states that on average the entropy production increases, even when we stop the process at a random stopping time. Key for the second law of thermodynamics at stopping times to hold is that the random time $\hat{T}$ obeys causality, i.e., it is a functional of the trajectory $X_0^{\hat{T}}$ up to the stopping time.

### 9.5.3 Gambler's ruin problem for edge currents

In general, currents of interest in nonequilibrium processes do not match the entropy production. Indeed, processes typically exhibit multiple currents, and the stochastic entropy production is a specific linear combination of those currents, see Eqs. (9.58-9.60). Therefore, in the present section, we consider the simplest example of a current that is not the stochastic entropy production, viz., an edge current $\hat{J} = \hat{J}(y|x)$. The corresponding
stopping problem is given by

$$ \hat{T}(y|x) = \min \left\{ n \ge 0 : \hat{J}_n(y|x) \notin (-J_-, J_+) \right\}. \qquad (9.80) $$

Note that in general macroscopic currents are not edge currents, but nevertheless, there exist interesting examples of processes for which this is the case. An example is the position of two-headed molecular motors bound to a biofilament as described by the model developed in Ref. [165], whose transition graph is shown in the Left Panel of Fig. 9.3. In this model, the position of the motor along the filament is given by the edge current $\hat{J}_n(5|2)$.

From the solution of the gambler's ruin problem for $\hat{S}$, we have learned that it is sufficient to find a martingale associated with $\hat{J}_n(y|x)$ to solve its gambler's ruin problem. However, $\exp(-a\hat{J}_n(y|x))$ is not a martingale, except for the trivial case of $a=0$. Nevertheless, as shown in Ref. [147], we can find a martingale that is asymptotically equivalent to $\exp(-a\hat{J}_n(y|x))$, i.e., up to terms of the order $O_n(1)$, for a certain value of $a$. Indeed, the process

$$ \hat{M}_n(y|x) := \frac{p_{\infty}^{(x,y)}(\hat{X}_0)q_{\infty}(\hat{X}_n)}{\underset{\infty}{\stackrel{(x,y)}{p}}(\hat{X}_n)p_{\infty}(\hat{X}_0)} \exp \left( -a^*(y|x)\hat{J}_n(y|x) \right) \qquad (9.81) $$

is a martingale, where

$$ a^*(y|x) := \ln \frac{\mathbf{R}_{yx} \mathbf{p}_{\infty}^{(x,y)}(x)}{\mathbf{R}_{xy} \mathbf{p}_{\infty}^{(x,y)}(y)} \qquad (9.82) $$

is the so-called effective affinity [86, 167, 168], where $p_{\infty}^{(x,y)}$ is the stationary state of the Markov chain $(\hat{X}, \mathbb{P}^{(x,y)})$ obtained from $(\hat{X}, \mathbb{P})$ by removing the edge $(x, y)$ from the Markov transition graph, and where $q_{\infty}$ is the stationary state of an auxiliary process $(\hat{X}, \mathbb{Q})$ with the transition matrix $\mathbf{V}$ with entries

$$ \mathbf{V}_{uv} = \begin{cases} \frac{\mathbf{p}_{\infty}^{(x,y)}(u)}{\mathbf{p}_{\infty}^{(x,y)}(v)} \mathbf{R}_{vu}, & \text{if } (u,v) \in \mathcal{E} \setminus \{(x,y), (y,x)\}, \\ \mathbf{R}_{uv}, & \text{if } (u,v) \in \{(x,y), (y,x)\}. \end{cases} \qquad (9.83) $$

In the Appendix 9.7, we perform a calculation analogous to the one in Ref. [147] for continuous time Markov chains showing that $\hat{M}_n(y|x)$ is a Radon-Nikodym derivative process of the form Eq. (9.29), and hence a martingale with respect to $\hat{X}$ and $\mathbb{P}$.

Consequently, Doob's optional stopping theorem in Sec. 9.3.3 applies to the martingale $\hat{M}_n(y|x)$ and the stopping time $\hat{T}(y|x)$, yielding

$$ \langle \hat{M}_{\hat{T}(y|x)}(y|x) | \hat{X}_0 = x_0 \rangle = \langle \hat{M}_0(y|x) | \hat{X}_0 = x_0 \rangle = \frac{\mathbf{q}_{\infty}(x_0)}{\mathbf{p}_{\infty}(x_0)}. \qquad (9.84) $$

Defining the splitting probabilities

$$ p_{(y|x)}^{\pm}(x_0) = \mathbb{P} \left( \hat{J}_{\hat{T}(y|x)}(y|x) = J_{\pm} | X_0 = x_0 \right), \qquad (9.85) $$

which obey

$$ p_{(y|x)}^{+}(x_0) + p_{(y|x)}^{-}(x_0) = 1, \qquad (9.86) $$

we find from Eqs. (9.84) and (9.86) that

$$ \bar{p}_{(y|x)}(x_0) = \frac{\exp(a^*(y|x)J_+)\zeta(x_0) - 1}{\exp(a^*(y|x)(J_- + J_+))\zeta(x) - 1}, \qquad (9.87) $$
where

$$ \zeta(u) = \frac{\mathbf{p}_{\infty}^{(x,y)}(y)\mathbf{q}_{\infty}(u)}{\mathbf{p}_{\infty}^{(x,y)}(u)\mathbf{q}_{\infty}(y)}, \quad \text{for } u \in \mathcal{X}. \qquad (9.88) $$

Analogously to the study of the infimum statistics of entropy production, Ref. [147] uses Eq. (9.87) to determine the probability mass function of the infimum $\hat{J}_{\text{inf}}(y|x)$ of $\hat{J}(y|x)$ and finds that it is geometrically distributed; this follows readily from taking the limit $J_+ \to \infty$, see also Ref. [169]. Remarkably, this property is unique to edge currents can thus be used to identify whether a measured current $\hat{J}$ is an edge current; note that we assume here that the observer can measure $\hat{J}$ but not $\hat{X}$.

Since the thermodynamic interpretation of the stationary distributions $p_{\infty}$, $p_{\infty}^{(x,y)}$, and $q_{\infty}$ is not entirely clear, we take the limits $J_+ \gg 1$ and $J_- \gg 1$ in Eq. (9.87), with the ratio $J_+/J_-$ fixed. In this limit, the formula (9.87) simplifies into

$$ p_{(y|x)}^{-}(x_0) = \exp\left(-a_{(y|x)}^{*} J_{-}(1 + o_{J_{\min}}(1))\right) \qquad (9.89) $$

with $J_{\min} = \min\{J_{-}, J_{+}\}$. Additionally, using the Wald-like equality [170,171],

$$ \langle \hat{T} \rangle = \frac{J_{+}}{j(y|x)}(1 + o_{J_{\min}}(1)), \qquad (9.90) $$

we get that

$$ |\log p_{(y|x)}^{-}(x_0)| = \frac{J_{-}}{J_{+}} a^{*}(y|x) \langle \hat{T} \rangle j(y|x) (1 + o_{J_{\min}}(1)). \qquad (9.91) $$

The Eq. (9.91) has an appealing thermodynamic meaning when realising that [86,168]

$$ a^{*}(y|x)j(y|x) \le \sigma, \qquad (9.92) $$

yielding [145]

$$ |\log p_{(y|x)}^{-}(x_0)| \le \frac{J_{-}}{J_{+}} \langle \hat{T} \rangle \sigma (1 + o_{J_{\min}}(1)). \qquad (9.93) $$

The Eq. (9.93) expresses a trade-off between speed ($\langle \hat{T} \rangle$), dissipation ($\sigma$), and uncertainty ($|\log p_{-}|$) in a nonequilibrium process $\hat{X}$. It is a universal inequality in the sense that it applies to arbitrary Markov chains $\hat{X}$ and, importantly, it remains valid for generic currents $\hat{J}$ as we discuss in the next subsection.

### 9.5.4 Gambler's ruin problem for generic currents

Although edge currents play a significant role in certain specific models — such as for the position of a molecular motor described by the model illustrated in the Left Panel of Fig. 9.3 — it is vital to study the first-passage problem of general currents in nonequilibrium processes, as macroscopic currents are typically not edge currents. This can be seen in the model illustrated in the Right Panel of Fig.9.3, which is an alternative model for a two-headed molecular motor that takes into consideration multiple pathways by which the motor can move forwards or backwards.

Therefore, we consider in this section the gambler's ruin problem $\hat{T}_J$, described by Eq. (9.56), for a general current $\hat{J}$ of the form given by Eq. (9.53).

It should be noted that the derivations we present in this section are based on continuous time Markov chains, as this is the setup considered in Ref. [145]. Therefore, we use the time index $t \in \mathbb{R}^+$ instead of $n \in \mathbb{N}$. Nevertheless, several results for the gambler's ruin problem reviewed here also apply to discrete time Markov chain; we will come back on this matter in more detail.