# Minimum Action Method for Nonequilibrium Phase Transitions

Ruben Zakine  Courant Institute, New York University, 251 Mercer Street, New
York, New York 10012, USA  Eric Vanden-Eijnden  Courant Institute, New York
University, 251 Mercer Street, New York, New York 10012, USA

(August 28, 2024)

###### Abstract

First-order nonequilibrium phase transitions observed in active matter, fluid
dynamics, biology, climate science, and other systems with irreversible
dynamics are challenging to analyze because they cannot be inferred from a
simple free energy minimization principle. Rather the mechanism of these
transitions depends crucially on the system’s dynamics, which requires us to
analyze them in trajectory space rather than in phase space. Here we consider
situations where the path of these transitions can be characterized as the
minimizer of an action, whose minimum value can be used in a nonequilibrium
generalization of the Arrhenius law to calculate the system’s phase diagram.
We also develop efficient numerical tools for the minimization of this action.
These tools are general enough to be transportable to many situations of
interest, in particular when the fluctuations present in the microscopic
system are non-Gaussian and its dynamics is not governed by the standard
Langevin equation. As an illustration, first-order phase transitions in two
spatially-extended nonequilibrium systems are analyzed: a modified Ginzburg-
Landau equation with a chemical potential which is non-gradient, and a
reaction-diffusion network based on the Schlögl model. The phase diagrams of
both systems are calculated as a function of their control parameters, and the
paths of the transitions, including their critical nuclei, are identified.
These results clearly demonstrate the nonequilibrium nature of the
transitions, with differing forward and backward paths.

##  I Introduction

Materials with identical microscopic constitutions can be found in very
different macroscopic states when external conditions, such as temperature or
pressure, vary. A major achievement of equilibrium statistical mechanics is to
give a first-principle explanation of these phase transitions. The theory
posits the existence of a distribution, for example of Boltzmann-Gibbs type,
that gives the probability of finding the microscopic system in any of its
possible configurations. Macroscopic properties like the system’s density,
magnetization [1], population number in the ground state [2], etc., can then
be deduced by enumerating all the microscopic configurations consistent with a
given value of the chosen macroscopic observable, and identifying which of
these values is most likely. More formally, if we denote by $\phi$ the
variable used to characterize the macroscopic state of the system, its
statistical weight is obtained by summing the system’s probability
distribution over all microscopic states consistent with a given realization
of $\phi$. Performing this calculation typically require sophisticated tools
such as renormalization group theory [3], the replica method or the cavity
method [4, 5], along with tools from large deviation theory [6]. It
generically shows that the statistical weight of $\phi$ is asymptotically
given by $\exp(-V(\phi)/\epsilon)$, where $V(\phi)$ is some free energy to be
calculated, and $\epsilon$ is a small parameter that tends to zero in the
thermodynamic limit when the number of microscopic constituents tends to
infinity: as a result, the theory predicts that the system will be found in
the macroscopic state $\phi$ of minimum free energy with probability one in
this limit. This also explains phase transitions: they take place when the
topology of the free energy $V(\phi)$ changes as a control parameter, like the
temperature or some applied external field, is varied. For example, if
$V(\phi)$ has two wells whose relative depths change with the control
parameter, a first-order phase transition occurs when the deepest well becomes
more shallow than the other well, and as a result, the macroscopic state of
the system changes from the first to the second.

The statistical mechanics approach to phase transitions rests on the
assumption that the probability distribution of the microscopic system is
known. This information is available for equilibrium systems, whose
microscopic dynamics is time-reversible. In this work, however, we are
primarily interested in nonequilibrium systems, whose dynamics is
irreversible. Except for some special situations where it can be computed
exactly [7, 8] or asymptotically e.g. via some thermodynamic mapping [9, 10,
11, 12], the invariant distribution of these systems is not known in general.
Yet, these systems too can undergo phase transitions. Examples include driven
systems arising from active matter [13, 14, 15, 16, 17, 18], fluid dynamics
[19, 20], biology [21], neuroscience [22, 23, 24], climate science [25, 26],
etc. The description of such nonequilibrium phase transitions requires a
generalization of the equilibrium statistical mechanics approach, in which we
must consider the probability of trajectories rather than configurations.

###  I.1 Minimum Action Principle

Table 1: Correspondence of formalisms |  Equilibrium approach in phase space |  Nonequilibrium approach in trajectory space
---|---|---
Microscopic weights |  Probability distribution (e.g. Boltzmann-Gibbs distribution) of the microscopic variables in state space. |  Probability distribution (e.g. path integral) of the microscopic trajectories.
Macroscopic variable |  Map from microscopic state space to coarse-grained macroscopic space (e.g. spin values to magnetization). |  Map from microscopic trajectories in phase space to macroscopic trajectories.
Macroscopic weights |  Marginal distribution of the macroscopic variables, and associated free energy |  Marginal distribution of the macroscopic trajectories, and associated action.
Macroscopic predictions |  Free energy minimization |  Action minimization

Even if the stationary distribution of nonequilibrium systems is not known in
general, we can often write down the probability distribution of their
trajectories, using e.g. path integral approaches such as the Martin-Siggia-
Rose-Janssen-De Dominicis [27, *janssen1976, *de_dominicis1976] or the Doi-
Peliti formalism [30, *peliti1985], or Girsanov theorem. This offers the
possibility to generalize the micro-to-macro mapping to trajectory space
rather than phase space: that is, enumerate all the microscopic trajectories
leading to the same evolution of a macroscopic variable, and thereby deduce
the probability weight of these macroscopic trajectories. While these
calculations are again to be performed on a case-by-case basis, by analogy
with the equilibrium setup we can deduce some of the generic features of the
result. Let us discuss those next.

Assuming again that the macroscopic state of the system can be described by a
variable or field $\phi$ in some differentiable manifold $\mathcal{M}$ (for
example $\mathbb{R}^{d}$, $\mathbb{S}^{d}$, or $L_{2}(\mathbb{R}^{d})$),
working in trajectory space amounts to calculating the probability weight of a
macroscopic path $\\{\phi(t)\\}_{t\in[0,T]}$ by enumerating the microscopic
trajectories consistent with $\\{\phi(t)\\}_{t\in[0,T]}$ and summing over
their probability distribution. Generically we expect the result of this sum
to indicate that the weight of the macroscopic path
$\\{\phi(t)\\}_{t\in[0,T]}$ is asymptotically given by the factor
$\exp(-S_{T}[\phi]/\epsilon)$, where $\epsilon$ is again a small parameter
that goes to zero as the number of microscopic constituents in the system goes
to infinity and $S_{T}[\phi]$ is an action, that is, a functional of
$\\{\phi(t)\\}_{t\in[0,T]}$. The specific form of this action depends on the
problem under consideration (examples will be given below), but it typically
takes the form of an integral over a Lagrangian

$S_{T}[\phi]=\int_{0}^{T}L(\phi,\dot{\phi})dt.$ (1)

where $\dot{\phi}=d\phi/dt\in T_{\phi}\mathcal{M}$. The action $S_{T}[\phi]$
is the non-equilibrium generalization of the free energy $V(\phi)$ and
minimization of $S_{T}[\phi]$ allows us to quantify the probability and
mechanism of various macroscopic events in the limit as $\epsilon\to 0$. In
particular:

1\. The probability that the system started in state $\phi_{a}$ at time $t=0$
ends up in state $\phi_{b}$ at time $T$, is obtained by summing
$\exp(-S_{T}[\phi]/\epsilon)$ over all paths with these end points. When
$\epsilon\ll 1$, the path with minimum action dominates this sum, which means
that the aforementioned probability is asymptotically given by

$\mathbb{P}(\phi(T)=\phi_{b}|\phi(0)=\phi_{a})\asymp\exp\left(-\inf
S_{T}[\phi]/\epsilon\right)$ (2)

where the minimization is taken over all paths $\\{\phi(t)\\}_{t\in[0,T]}$
such that $\phi(0)=\phi_{a}$ and $\phi(T)=\phi_{b}$ and $\asymp$ means
exponential asymptotics, i.e the ratio of the logarithm of both sides in (2)
tends to 1 as $\epsilon\to 0$. The minimizer of the action also gives the
pathway by which the macroscopic transition event occurs with probability one
in this limit.

2\. The non-equilibrium invariant distribution of the system can be
characterized similarly via the quasipotential defined as

$V_{\phi_{a}}(\phi_{b})=\inf_{T>0}\inf S_{T}[\phi]$ (3)

where the inner minimization is again taken over all paths
$\\{\phi(t)\\}_{t\in[0,T]}$ such that $\phi(0)=\phi_{a}$ and
$\phi(T)=\phi_{b}$. The quasipotential $V_{\phi_{a}}(\phi_{b})$ plays a role
analogous to the free energy barrier from state $\phi_{a}$ to $\phi_{b}$, and
it can be used to identify the possible phases and formulate an equivalent of
Arrhenius law. More precisely: $\phi_{a}$ is a metastable phase if
$V_{\phi_{a}}(\phi)\leq V_{\phi}(\phi_{a})$ for all $\phi$ in a vicinity of
$\phi_{a}$ (i.e. $\phi_{a}$ is the non-equilibrium equivalent of a local
minimum on the free energy); and if $\phi_{a}$ and $\phi_{b}$ are the only two
metastable phases in the system, the asymptotic rates of transition from
$\phi_{a}$ to $\phi_{b}$ and $\phi_{b}$ to $\phi_{a}$ are respectively given
by

$k_{a,b}\asymp e^{-V_{\phi_{a}}(\phi_{b})/\epsilon},\quad k_{b,a}\asymp
e^{-V_{\phi_{b}}(\phi_{a})/\epsilon}$ (4)

3\. Eq. (4) is a non-equilibrium generalization of Arrhenius law. It implies
that the relative probability to find the system in states $\phi_{a}$ or
$\phi_{b}$ on its non-equilibrium invariant distribution is asymptotically
given by

$\mathbb{P}(\phi_{b})/\mathbb{P}(\phi_{a})\asymp
e^{-[V_{\phi_{a}}(\phi_{b})-V_{\phi_{b}}(\phi_{a})]/\epsilon}.$ (5)

As a result, with probability 1 as $\epsilon\to 0$, the system is in state
$\phi_{a}$ if $V_{\phi_{a}}(\phi_{b})>V_{\phi_{b}}(\phi_{a})$ and state
$\phi_{b}$ if $V_{\phi_{a}}(\phi_{b})<V_{\phi_{b}}(\phi_{a})$. By analyzing
how the quasipotential varies in terms of the system’s control parameters we
can thereby identify non-equilibrium phase transitions that arise when
$V_{\phi_{b}}(\phi_{a})=V_{\phi_{a}}(\phi_{b})$ and characterize their
mechanism—the details of these calculations will be given below. We also refer
the reader to Fig. 1 for a graphical illustration in a toy system of non-
equilibrium phase transition whose detection requires the formalism above.

Figure 1: Example of nonequilibrium phase transition in a toy two-dimensional
system with two metastable states (black disks) and action
$S_{T}[x,y]=\frac{1}{2}\int_{0}^{T}\left(|\dot{x}-f(x,y)|^{2}+|\dot{y}-f(y,x)|^{2}\right)dt$,
where $f(x,y)=D(y-x)+x-x^{3}+h+\nu(x-y)^{2}$. The parameters $D$ and $h$ are
fixed to $0.5$ and $-0.1$, respectively, while $\nu$ is used as a control
parameter. The flow of the vector field $f$ for $\nu=0.5$ is shown in the top
panel and $\nu=1.5$ in the bottom panel. This flow is non-gradient (i.e. the
deterministic dynamics is not steepest descent over an energy) has two stable
fixed points, $\phi_{a}$ (bottom left) and $\phi_{b}$ (top right), which solve
$f(x,y)=f(y,x)=0$, which are the possible phases in this toy example. The blue
line represents the most probable transition path (i.e. the minimizer of the
action) from $\phi_{a}$ to $\phi_{b}$ and the orange line is the most probable
path from $\phi_{b}$ to $\phi_{a}$. The thickest line indicates the path with
larger rate: that is, $\phi_{a}$ is the stable phase when $\nu=0.5$ (top
panel), while $\phi_{b}$ is the stable phase when $\nu=1.5$ (bottom panel).
The transition paths were calculated with the method developed in this paper,
and cross-checked using GMAM [32].

###  I.2 Hamiltonian formalism

The minimum action principle described in the last section offers a way to
study transition events and phase transitions in non-equilibrium systems.
Concrete predictions however rest on our ability to: (i) derive the Lagrangian
used in the action (1) and (ii) minimize this action as needed in (2) and (3).

Like in the equilibrium case, resolving the first issue is again complicated
in general and requires to be handled on a case-by-case basis. When these
calculations can be done (see Sec. II for a list of examples), one often
deduces that $L(\phi,\dot{\phi})$ is given as the Legendre-Fenchel transform
of a Hamiltonian $H(\phi,\theta)$:

$L(\phi,\dot{\phi})=\sup_{\theta\in
T_{\phi}\mathcal{M}}\left(\langle\theta,\dot{\phi}\rangle-H(\phi,\theta)\right),$
(6)

where $\langle\cdot,\cdot\rangle$ denotes the scalar product in
$T_{\phi}\mathcal{M}$ and $\theta$ is a field conjugate to $\phi$ whose
physical meaning will be explained below. The form of the Hamiltonian
$H(\phi,\theta)$ is also problem-dependent but it is known in some instances,
see Sec. II.1.

Using (6), the minimization of the action can then be formulated as a min-max
problem:

$\inf_{\phi}S_{T}(\phi)=\inf_{\phi}\sup_{\theta}\int_{0}^{T}\left(\langle\dot{\phi},\theta\rangle-H(\phi,\theta)\right)dt.$
(7)

where the supremum is taken over all $\\{\theta(t)\\}_{t\in[0,T]}$ and the
infimum over all $\\{\phi(t)\\}_{t\in[0,T]}$ such that $\phi(0)=\phi_{a}$ and
$\phi(T)=\phi_{b}$. To get the quasipotential, we must also consider an extra
minimization over all $T>0$ to (7), while other applications may require
adding terms to (7) or modifying the boundary conditions for this min-max
problem. In most cases, these calculations must be performed numerically. One
of the main goal of this paper is to develop robust numerical methods for
these computations. These methods aim to be general enough to be applicable to
a wide variety of systems that fit the framework above; here we will also use
them to solve some non-trivial examples involving spatially-extended systems
undergoing non-equilibrium phase transitions.

###  I.3 Related Works

The Euler-Lagrange equations associated with the min-max problem (7) are
Hamilton’s classical equations:

$\dot{\phi}=\partial_{\theta}H,\qquad\dot{\theta}=-\partial_{\phi}H.$ (8)

What makes the problem non-standard, however, are the boundary conditions
imposed on $\phi(t)$ at $t=0$ and $t=T$. The nature of these boundary
conditions suggests to use shooting methods [33], as was proposed e.g. in
[34], but such methods scale badly with dimension, or can even be ill-posed
for the problems we are interested in, for which the equation for $\theta$ in
(8) cannot be integrated forward in time. Shooting methods are also hard to
use when $T=\infty$, which typically arises when we consider $\inf_{T>0}\inf
S_{T}[\phi]$.

To get around this difficulty, the minimum action method (MAM) proposed in
[35] evolves the whole trajectory $\\{\phi(t)\\}_{t\in[0,T]}$ while keeping
$\phi(0)=\phi_{a}$ and $\phi(T)=\phi_{b}$ fixed. This amounts to performing
gradient descent (GD) on the action in the landscape of all authorized paths
satisfying these boundary conditions. Introducing the artificial optimization
time $\tau$, GD results in the following evolution equation for
$\\{\phi(\tau,t)\\}_{\tau\geq 0,t\in[0,T]}$

$\displaystyle\partial_{\tau}\phi=-\frac{\delta
S_{T}[\phi]}{\delta\phi(t)},\quad\phi(\tau,0)=\phi_{a},\quad\phi(\tau,T)=\phi_{b},$
(9)

or using the Lagrangian formulation of the action

$\displaystyle\partial_{\tau}\phi=-\left[\frac{\partial
L}{\partial\phi}-\frac{d}{dt}\left(\frac{\partial
L}{\partial\dot{\phi}}\right)\right],$ (10)

with the same boundary conditions at $t=0,T$. The main drawback of MAM is that
it involves the Lagrangian $L(\phi,\dot{\phi})$ rather than the Hamiltonian
$H(\phi,\theta)$, and there are many problems of interest where the latter is
explicitly available but the former is not. Solving (10) then requires one to
perform $\max_{\theta}[\langle\dot{\phi},\theta\rangle-H(\phi,\theta)]$ for
all $t\in[0,T]$ at each iteration step in $\tau$ to numerically get an
estimate of the function $\vartheta(\phi,\dot{\phi})$ such that

$L(\phi,\dot{\phi})=\langle\dot{\phi},\vartheta(\phi,\dot{\phi})\rangle-H(\phi,\vartheta(\phi,\dot{\phi}))$
(11)

Proceeding similarly, we can also obtain numerical estimates for the
derivative $\partial L/\partial\phi$ and $(d/dt)\partial L/\partial\dot{\phi}$
appearing on the right hand-side of (10). This approach was used in [36]. The
downside is that (10) is a partial differential equation in physical time $t$,
optimization time $\tau$, and possibly space as well when $\phi$ and $\theta$
are fields; writing efficient numerical solvers for such equations typically
requires one to use implicit schemes for numerical stability and/or
efficiency, and such schemes are hard to design without explicit knowledge of
$\vartheta(\phi,\dot{\phi})$. This is why here we want to bypass the
computation of $\vartheta(\phi,\dot{\phi})$ and solve the min-max problem in
(7) concurrently by treating the minimization over $\phi$ and the maximization
over $\theta$ on equal footings.

If we now turn our attention to the quasipotential (3), an extra minimization
over all $T>0$ must be added to the min-max problem (7). The MAM can be
generalized to handle this problem by using a reparametrization of the path
$\\{\phi(t)\\}_{t\in[0,T]}$ by arclength rather than physical time. This
formulation leads to the geometric minimum action method (GMAM), in which the
minimization over $T$ is performed explicitly beforehand [32]. GMAM was
further developed and used in [37, 36] and a variant of it was also recently
proposed in [38] to compute the quasipotential by a spectral decomposition of
paths and an optimization of basis coefficients. However, these methods are
Lagrangian based, with the issues discussed before. We will show below how to
use the ideas behind GMAM in a Hamiltonian approach and thereby get an
efficient method to solve $\inf_{T>0}\inf_{\phi}S_{T}[\phi]$.

On the analytical side, the quasipotential is related to the solution
$V(\phi)$ of the Hamilton-Jacobi equation [39, 40]

$0=H(\phi,\partial V/\partial\phi).$ (12)

This equation can be used to deduce some structural properties of the
quasipotential. It is a central object in problems that can be tackled through
macroscopic fluctuation theory (MFT), and the review [41] provides numerous
examples and useful insights. Ref. [42] also discusses how to perturbatively
solve this equation in systems that are close to equilibrium. In general,
however, (12) needs to be solved numerically, which is nontrivial since it is
a complicated partial differential equation (or even a functional equation
when $\phi$ is a field). In dimension 2 or 3 this can be done globally using
fast marching methods like the one discussed in [43, *gan_cameron_2017,
*dahiya_cameron_2018]. In dimension higher than 3, these methods become
inapplicable, and (12) must be solved locally by the method of characteristics
using the variational formulation of this equation: this brings us back to
solving $\inf_{T>0}\inf_{\phi}S_{T}[\phi]$.

###  I.4 Organization and Main Contributions

The remainder of this paper is organized as follows: In Sec. II, we give more
details about the minimum action principle that is at the core of our
approach. To this end, in Sec. II.1 we first present typical classes of
stochastic dynamical systems that display metastability and non-equilibrium
phase transitions, and are amenable to analysis via action minimization. For
illustration, we also use our method to compute the transition paths between
metastable states in two low-dimensional benchmark models, namely the Maier-
Stein model [34] and the Schlögl model [46]. In Sec. II.2 we discuss the
features of the minimum action framework in a general setup, and summarize the
main outputs of the approach.

In Sec. III we present our method for solving the min-max problem (7). The
scheme is based on performing gradient descent-ascent (GDA) on the objective,
which has the advantage that it can be formulated directly in the Hamiltonian
setup. The main issue we need to address is that, in our context, the GDA
equations are partial differential equations of hyperbolic type in some
optimization time $\tau$ and physical time $t$, to be solved with boundary
conditions at $t=0$ and $t=T$. In Sec. III.1 we show that a simple linear
change of variables in these equations allows us to reformulate them in a way
that is convenient for numerical solution via Strang splitting. In Sec. III.2
we then generalize this scheme to compute the quasipotential when an extra
minimization over $T$ is added to (7). This is done by using a geometric
formulation in which the paths $\\{\phi,\theta\\}_{t\in[0,T]}$ are
reparametrized using normalized arclength. This allows us to handle the
minimization of $T$ efficiently and calculate paths whose duration in physical
time is infinite but whose length remains finite.

In Secs. IV and V we then use the methods we propose to analyze two spatially
extended nonequilibrium systems that display first-order phase transitions.
More specifically, in Sec. IV we study a modified Ginzburg-Landau (GL)
dynamics subject to an additive Gaussian white noise. The noiseless evolution
of the field is non-gradient with two stable fixed points; the noise makes
these points metastable and we must resort to minimum action algorithms to
compute the nonequilibrium transition pathways between them. The action along
these paths allows us to estimate the relative probability of the metastable
states. We use this procedure to compute the phase diagram of the system in
function of two control parameters. In Sec. V, we study a spatially extended
version of the Schlögl model, in which the fluctuations are driven both by
diffusion of the microscopic molecules and reactions between them. This system
displays a first-order non-equilibrium phase transition in terms of the
diffusivity of the molecules, which we characterize. We also show that the
predictions of the minimum action approach explain the transition events
observed in the microscopic system.

Concluding remarks are given in Sec. VI and some technical developments are
deferred to several Appendices.

##  II Problem Setup and Interpretation

The aim of this section is to provide a better motivation of the minimum
action principle introduced in the introduction and pinpoint some of its key
predictive features. For the reader’s convenience we begin by listing a
collection of motivating problems where the formalism applies. In Sec. II.2 we
will then put the approach in a broader context and explain how to use it to
analyze metastability.

###  II.1 A collection of motivating problems

The first two examples a and b involve no coarse-graining from micro-to-macro
and are included because they are simple and transparent; the last three
examples c, d, and e requires one to define proper macroscopic variables to
derive the minimum action principle and its Hamiltonian.

##### Diffusion in detailed balance:

Consider the motion of a particle $x(t)\in\mathbb{R}^{d}$ whose evolution is
governed by the overdamped Langevin equation

$\dot{x}=-\nabla U(x)+\sqrt{2kT}\,\eta(t)$ (13)

where $U(x)$ is some potential, $kT$ is the product of Boltzmann constant $k$
and the temperature $T$, and $\eta(t)$ is a white-noise. The dynamics (13) is
in detailed balance with respect to the Boltzmann-Gibbs probability density
function $\rho(x)=Z^{-1}e^{-U(x)/kT}$ where
$Z=\int_{\mathbb{R}^{d}}e^{-U(x)/kT}dx$. If the potential $U(x)$ has multiple
local minima, and the temperature $kT$ is much smaller than the barriers
between them, (13) displays metastability: the system stays confined for a
long time in the well around a minimum of $U(x)$ before finally hopping to
another well where the process repeats. In this example, we can use a WKB
expansion to analyze the Fokker-Planck equation associated with (13). The
eikonal equation obtained at leading order in $kT$ is a Hamilton-Jacobi
equation whose Hamiltonian is given by

$H(x,\theta)=-\langle\theta,\nabla U(x)\rangle+|\theta|^{2}.$ (14)

This Hamiltonian is the one to be used in the min-max problem (7), as can also
be proven rigorously using Freidlin-Wentzell LDT [39]. In this example, the
quasipotential $V_{x_{a}}(x_{b})$ can be calculated explicitly. If $x_{a}$ and
$x_{b}$ are the locations of two local minima of $U(x)$ with adjacent wells,
the path minimizing the action (3) is the minimum energy path between these
two points, and $V_{x_{a}}(x_{b})$ is given by

$V_{x_{a}}(x_{b})=U(x_{s})-U(x_{a})$ (15)

where $U(x_{s})$ is the energy of the saddle point of minimum height (aka
mountain pass) between $x_{a}$ ad $x_{b}$. Thus we recover the Arrhenius law
for the rate of transition from $x_{a}$ to $x_{b}$,

$k_{a,b}\asymp e^{-V_{x_{a}}(x_{b})/kT}=e^{-[U(x_{s})-U(x_{a})]/kT}.$ (16)

a)

b)

Figure 2: Optimal reaction path in the Maier-Stein model [34]: the evolution
of $(x,y)\in\mathbb{R}^{2}$ is governed by the SDE $\dot{x}=x-x^{3}-\beta
xy^{2}+\sqrt{\epsilon}\eta_{1}$,
$\dot{y}=-(1+x^{2})y+\sqrt{\epsilon}\eta_{2}$, where $\eta_{1,2}$ are
independent Gaussian white noises, and $\beta$ some parameter (here
$\beta=10$). The flow lines of the noiseless model ($\epsilon=0$) are shown as
grey lines in panel a), and the background color indicates the magnitude of
the drift (darker means larger): there are two stable fixed points
$x_{a}=(-1,0)$ and $x_{b}=(1,0)$ (black disks), and one unstable fixed point
$x_{c}=(0,0)$ (black cross). When $\epsilon$ is small but finite, these two
fixed points become metastable, and the noise induces transitions between
them: the most likely path from $x_{a}$ to $x_{b}$ is the minimum action path,
which is shown as the blue line. Also shown as an red dashed line is the
heteroclinic orbit between $x_{a}$ and $x_{b}$: this orbit is different from
the minimum action path, indicative of a nonequilibrium transition. Panel b)
shows the increment of the action along the minimum action path and the
heteroclinic orbit, confirming that the former is more likely.

a)

b)

Figure 3: Optimal reactions paths in the bistable Schlögl model presented in
Sec. V, with two reactive compartments and where the particles are allowed to
jump from one compartment to another at rate $\gamma$. This system belongs to
the class of reaction networks. Here the Hamiltonian is given by
$H=H^{R}+H^{D}$ with
$H^{R}=\sum_{i=1}^{2}w^{+}(x_{i})(e^{\theta_{i}}-1)+w^{-}(x_{i})(e^{-\theta_{i}}-1)$
and $H^{D}=\sum_{i=1}^{2}\gamma
x_{i}(e^{\theta_{i-1}-\theta_{i}}+e^{\theta_{i+1}-\theta_{i}}-2)$. The flow
lines of the mass-action law are shown in panel a), along with its stable
fixed points (black dots) and unstable fixed points (black crosses). When the
number of agents is large but finite, these stable fixed points become
metastable states, and the most likely paths between them are shown as full
lines in blue and orange. Also shown in dashed red line is the heteroclinic
orbit. All three paths differ, indicative of a a nonequilibrium transition.
Panel b) shows the increment along the action of the two minimum action paths,
indicating that the bottom left state is most stable under random
fluctuations.

##### Diffusion out-of-equilibrium:

The picture above can be generalized to systems whose evolution is described
by the stochastic differential equation

$\dot{x}=b(x)+\sqrt{\epsilon}\sigma(x)\eta(t),$ (17)

even if this equation is not in detailed balance, i.e. it is not possible to
write the drift $b(x)$ as $-D(x)\nabla U(x)+kT\nabla\cdot D(x)$ for
$D(x)=(\sigma\sigma^{T})(x)$ and some potential $U(x)$. Metastability is
observed with (17) in situations where the noiseless deterministic system
$\dot{x}=b(x)$ has multiple stable fixed points and the noise amplitude is
small but finite: the system hovers for long times in the basin around one of
these fixed points, but a noise-driven transition to another basin eventually
occurs. These transitions can be described by the minimum action principle
using the Hamiltonian which can again be obtained via WKB analysis of the
Fokker Planck equation associated with (17). It is given by

$H(x,\theta)=\langle b(x),\theta\rangle+\tfrac{1}{2}|\sigma(x)\theta|^{2}.$
(18)

Given two stable fixed point $x_{a}$ and $x_{b}$ of $\dot{x}=b(x)$ with
adjacent basins of attraction, it is no longer possible in general to solve
(3) analytically and calculate $V_{x_{a}}(x_{b})$—to do so requires numerical
tools of the type developed below. Still, we know that the rates of transition
from $x_{a}$ to $x_{b}$ and $x_{b}$ to $x_{a}$ satisfy respectively

$k_{a,b}\asymp e^{-V_{x_{a}}(x_{b})/\epsilon},\qquad k_{b,a}\asymp
e^{-V_{x_{b}}(x_{a})/\epsilon},$ (19)

and, with probability $1$ as $\epsilon\to 0$, the system performs the
transition by following the optimal path minimizing (3)–for an illustration in
the context of Maier-Stein model [34], see Fig. 2. The results in SDE with
small noise of this type can be made rigorous using Freidlin-Wentzell theory
of large deviations (LDT) [39].

##### Reaction networks:

Consider a well-stirred chemical network between $M$ chemical species, where
the quantity of species $i$ is denoted $X_{i}$. Define the population vector
$X=(X_{1},\ldots,X_{M})^{T}$, and assume that there are $R$ reaction channels
with rates $w_{j}(x)$ and change (stoichiometric) vectors
$\nu_{j}\in\mathbb{Z}^{M}$, i.e.

$X\xrightarrow[]{w_{j}(X)}X+\nu_{j},\qquad j=1,\ldots,R.$ (20)

When the typical number of agents, $\Omega$, tends to infinity, the dynamics
of $x=X/\Omega$ is captured by the mass-action law

$\dot{x}=\sum_{j=1}^{R}w_{j}(x)\nu_{j},$ (21)

and the minimum action principle is useful to quantify the effects of
fluctuations when the number of agents is large but finite. In particular,
metastability arises if (21) has multiple stable fixed points (see Ref. [47,
48]), and it can be analyzed using the Hamiltonian

$H(x,\theta)=\sum_{j=1}^{R}w_{j}(x)\left(e^{\langle\nu_{j},\theta\rangle}-1\right).$
(22)

Here too this Hamiltonian can be obtained rigorously via Freidlin-Wentzell
LDT, or formally via WKB analysis of the system’s master equation [47] (see
also the appendix of [16] for a pedagogical derivation), or via a Doi-Peliti
field theory computation [30, *peliti1985, 49, 50, 51]. If $x_{a}$ and $x_{b}$
denote two stable fixed points of (21) with adjacent basins of attraction, the
rates of transition from $x_{a}$ to $x_{b}$ and $x_{b}$ to $x_{a}$ are
respectively given by

$k_{a,b}\asymp e^{-\Omega V_{x_{a}}(x_{b})},\qquad k_{b,a}\asymp e^{-\Omega
V_{x_{b}}(x_{a})},$ (23)

and with probability $1$ as $\Omega\to\infty$, when the networks performs the
transition, $X/\Omega$ follows the optimal path minimizing (3). This
minimization needs again to be performed numerically in general. As an
illustrating example, Fig. 3 displays reaction paths in the bistable Schlögl
model with two reactive compartments: this model belongs to the class of
reaction-diffusion networks that will be properly introduced in Sec. V.1. Note
that if we reduce the number of compartments to one, the quasipotential of
this model can be explicitly obtained from the Hamiltonian, see Fig. 4.

##### Interacting particle systems:

Consider $N$ particles $x_{i}\in\Lambda\subset\mathbb{R}^{d}$, $i=1,\ldots,N$,
that evolve according to

$\dot{x}_{i}=b(x_{i})+\frac{1}{N}\sum_{j=1}^{N}k(x_{i},x_{j})+\sigma(x_{i})\,\eta_{i}(t),$
(24)

where $b(x)$ is a drift as in (17), $k(x,y)$ is some interaction kernel, and
$\eta_{i}(t)$ are independent white-noises. To analyze such interacting
particle systems, it is convenient to introduce the empirical density of the
particles, $\rho_{N}(t,x)=N^{-1}\sum_{i=1}^{N}\delta(x-x_{i}(t))$. As
$N\to\infty$, this empirical density converges towards the density $\rho(t,x)$
that satisfies McKean-Vlasov equation

$\partial_{t}\rho=-\nabla\cdot\left(B(x,[\rho])\rho\right)+\tfrac{1}{2}\nabla\nabla:(D(x)\rho)$
(25)

where $B(x,[\rho])=b(x)+\int_{\mathbb{R}^{d}}k(x,y)\rho(y)dy$ and
$D(x)=(\sigma\sigma^{T})(x)$. Large fluctuations away from the mean-field
dynamics (25) can be captured by the minimum action principle, by using the
Hamiltonian

$\displaystyle H(\rho,\theta)$
$\displaystyle=\int_{\mathbb{R}^{d}}\nabla\theta(x)\cdot b(x)\rho(x)dx$ (26)
$\displaystyle+\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\nabla\theta(x)\cdot
k(x,y)\rho(x)\rho(y)dxdy$
$\displaystyle+\frac{1}{2}\int_{\mathbb{R}^{d}}|\sigma(x)\nabla\theta(x)|^{2}\rho(x)dx.$

This Hamiltonian can again be derived rigorously using LDT [52], and it
formally follows from a WKB analysis, now performed on the functional master
equation for the empirical particle density.

a)

b)

Figure 4: a) Contour plot of the Hamiltonian $H(\phi,\theta)$ of the bistable
Schlögl model introduced in Sec. V, Eq. (71), with a single reactive
compartment containing a large number of particles $\Omega$. The two black
disks show the stable fixed points $\phi_{a}$ (left disk) and $\phi_{b}$
(right disk) of the law of mass action valid when $\Omega\to\infty$, and the
orange square shows its unstable fixed point at $\phi_{s}=1$. The solid and
dashed white line shows the level set $H(\phi,\theta)=0$ on which the solution
to Hamilton’s equations (8) evolve in the limit when $T\to\infty$. The
horizontal dashed line at $\theta=0$ corresponds to the noiseless dynamics of
the mass-action law by which the system relaxes to one of the two stable fixed
points $\phi_{a}$ or $\phi_{b}$ by moving in the direction of the white
arrows. The solid white line, to follow in the direction of the black arrows,
shows the solutions of Hamilton’s equations on $H(\phi,\theta)=0$ at
$\theta\neq 0$: these trajectories indicate the most probable path by which
the fluctuations can drive the system away from $\phi_{a}$ and $\phi_{b}$, and
induce transitions when reaching $\phi_{s}$. In this example, the action
required to drive the system away from $\phi_{a}$ to some $\phi\leq\phi_{s}$
in its basin of attraction is the area between the white and the dashed solid
lines, going from $\phi_{a}$ to $\phi$: this area is also the quasipotential
$V_{\phi_{a}}(\phi)$ for $\phi\leq\phi_{s}$, and $V_{\phi_{b}}(\phi)$ for
$\phi\geq\phi_{s}$ can be defined similarly. b) Global quasipotential
$V(\phi)$ obtained by gluing the two quasipotentials $V_{\phi_{a}}(\phi)$ and
$V_{\phi_{b}}(\phi)$ at $\phi=\phi_{s}$ after vertical shifting; $V(\phi)$ is
also the viscosity solution of the Hamilton-Jacobi equation (12).

##### Fast-slow systems

Consider a system made of a pair of variables $(x,y)\in\mathbb{R}^{d\times D}$
whose evolution is governed by

$\left\\{\begin{aligned} \dot{x}&=f(x,y)\\\
\dot{y}&=\alpha^{-1}b(x,y)+\alpha^{-1/2}\sigma(x,y)\eta(t)\end{aligned}\right.$
(27)

where $\alpha>0$ measure the separation of time-scale between $x$ and $y$. For
small $\alpha$, this separation is large and $y$ evolves much faster than $x$.
In particular, when $\alpha\to 0$, the dynamics of $x$ is effectively captured
by the deterministic limiting equation

$\dot{x}=F(x)$ (28)

where $F(x)$ is obtained by averaging $f(x,y)$ over the stationary
distribution of the equation for $y$ at $x$ fixed, assuming that it exists.
This equation defines the so-called virtual fast process, which on the fast
time scale $\tau=\alpha t$ reads

$dy_{x}/d\tau=b(x,y_{x})+\sigma(x,y_{x})\eta(\tau)\quad\text{({\hbox{x}}
frozen)}$ (29)

Denoting by $\mathbb{E}^{x}$ the expectation over the stationary distribution
of $y_{x}(\tau)$, the function $F$ entering (28) is given by

$F(x)=\mathbb{E}^{x}f(x,y_{x}).$ (30)

If we want to analyze the effect of fluctuations on the dynamics of $x$ when
$\alpha$ is small but finite, we can use the minimum action principle with a
Hamiltonian that can be derived from LDT [53, 54, 55, 56]

$H(x,\theta)=\log\mathbb{E}^{x}\exp\left(\langle\theta,f(x,y_{x})\rangle\right)$
(31)

In general, this Hamiltonian will need to be calculated numerically, which
makes slow-fast systems of the type above more difficult to treat than the
models previously discussed in this section. Still, provided that we can
design some numerical routine to estimate $H$ as well as its derivatives
$\partial_{x}H$ and $\partial_{\theta}H$, the numerical methods presented
below are applicable in the context of slow-fast systems too. In models of
this type $\epsilon=\alpha$.

###  II.2 General setup

The problems listed in Sec. II.1 all share a Hamiltonian with the following
features:

  * A1.

$H(\phi,0)=0$ for all $\phi\in\mathcal{M}$;

  * A2.

$H(\phi,\theta)$ is strictly convex in $\theta$ for all $\phi\in\mathcal{M}$.

  * A3.

$H(\phi,\theta)$ is twice differentiable in both its arguments;

Assumptions A1 and A2 follow from the fact, generically, the Hamiltonian can
be expressed as a cumulant generating function, i.e. in the form of an
expectation generalizing (31):

$H(\phi,\theta)=\log\mathbb{E}^{\phi}\exp\left(\langle\theta,F(\phi,y_{\phi})\rangle\right)$
(32)

where $F$ is problem-dependent and the expectation is taken over the
statistics of some underlying process $y_{\phi}$ conditional on $\phi$ being
fixed. Assumption A3 is added for simplicity, as it guarantees that Hamilton’s
equations (8) are well posed.

The aim of this section is to discuss at generic level the meaning we can give
to the min-max problem (7) assuming that the Hamiltonian satisfies these
assumptions. In particular, we will show that the conjugate field $\theta$
appearing in the min-max problem (7) can be generically interpreted as
measuring the mean effects of the fluctuations in the system dynamics needed
to achieve some rare event, and the minimum of the action as the total cost of
these fluctuations from which the probability of the event can be estimated as
well as its mechanism and rate (if we add the minimization over $T>0$) .

####  II.2.1 Mean behavior

To begin, notice that for the interpretation that $\theta$ measures the
effects of the fluctuations to be consistent, the stochastic system under
consideration should, in some appropriate limit in which the fluctuations
disappear, satisfy the deterministic evolution equation obtained by setting
$\theta=0$ in Hamilton’s equations (8):

$\dot{\phi}=\partial_{\theta}H(\phi,\theta=0).$ (33)

For example, returning to the problems mentioned in Sec. II.1, (33) reduces to
the ODE $\dot{x}=b(x)$ as $\epsilon\to 0$ in the SDE (17) with Hamiltonian
(18); to the law of mass action (21) as $\Omega\to\infty$ for the reaction
network (20) with Hamiltonian (22); to the McKean-Vlasov equation as
$N\to\infty$ for the interacting particle system (24) with Hamiltonian (26);
and to (28) as $\alpha\to 0$ in the slow-fast system (27) with Hamiltonian
(31). More generally, under our assumptions, the solution of (33) is indeed a
special solution with $\theta(t)=0$ of Hamilton’s equations (8) since
Assumptions A1 and A3 imply that $\partial_{\phi}H(\phi,0)=0$ for all $\phi$.

####  II.2.2 Impact of the fluctuations

At the same time, the solution to the deterministic evolution (33) with
$\phi(0)=\phi_{a}$ does not satisfies $\phi(T)=\phi_{b}$ in general.
Therefore, for general boundary conditions $\phi(0)=\phi_{a}$ and
$\phi(T)=\phi_{b}$, the solution of (8) must have $\theta(t)\not=0$—in the
minimum action framework this is a reflection that in the original system
fluctuations are needed to drive the system’s trajectory away from the
solution of (33), and the value of $\theta(t)\not=0$ allows us to quantify the
cost/probability of observing the event $\phi(T)=\phi_{b}$ given that
$\phi(0)=\phi_{a}$. Specifically, under the strict convexity Assumption A2, we
have

$\langle\theta,\partial_{\theta}H(\theta,\phi)\rangle-H(\phi,\theta)\geq-H(\phi,0)=0,$
(34)

with equality if and only if $\theta=0$. Since
$\dot{\phi}=\partial_{\theta}H(\phi,\theta)$ along the solution to Hamilton’s
equations (8), we deduce that along this solution

$\inf_{\phi}S_{T}[\phi]=\inf_{\phi}\int_{0}^{T}\left(\langle\theta,\partial_{\theta}H(\phi,\theta)\rangle-H(\phi,\theta)\right)dt\geq
0$ (35)

Thus the action $S_{T}[\phi]$ can be indeed interpreted as a cost, which is
zero only if $\theta(t)=0$ (i.e. when the event can occur without
fluctuations), and is strictly positive otherwise (i.e. when the event
requires fluctuations). For the problems listed in Sec. II.1 the form of the
integrand $\langle\theta,\partial_{\theta}H\rangle-H$ is
$|\sigma(x)\theta|^{2}\geq 0$ for the SDE (17);
$\sum_{j=1}^{R}a_{j}(x)\left(\langle\nu_{j},\theta\rangle
e^{\langle\nu_{j},\theta\rangle}-e^{\langle\nu_{j},\theta\rangle}+1\right)\geq
0$ for the reaction network (20);
$\frac{1}{2}\int_{\mathbb{R}^{d}}|\sigma(x)\nabla\theta(x)|^{2}\rho(x)dx\geq
0$ for the interacting particle system (24); and
$\mathbb{E}^{x}f(x,y_{x})e^{\langle\theta,f(x,y_{x})\rangle}/\mathbb{E}^{x}e^{\langle\theta,f(x,y_{x})\rangle}-H$
for the slow-fast system (27).

In this interpretation, minimizing $S_{T}[\phi]$ amounts to minimizing the
cost of the fluctuations or, equivalently, finding the most likely fluctuation
that drives the event $\phi(T)=\phi_{b}$ given that $\phi(0)=\phi_{a}$,
leading to the asymptotic estimate (2) for the probability of the event.

####  II.2.3 Long-time limit.

Turning now our attention to the problem $\inf_{T\geq
0}\inf_{\phi}S_{T}[\phi]$, its interpretation is easiest if we assume that the
noiseless equation (33) has $N$ stable fixed points $\phi_{1}$, $\phi_{2}$, …,
$\phi_{N}$, and the basins of attraction of these fixed points under (33),
denoted respectively as $B_{1},\ldots B_{N}$, partition
$\mathcal{M}=\cup_{i=1}^{N}\bar{B}_{i}$ (note that $B_{i}\cap B_{j}=\emptyset$
if $i\not=j$ by definition). In this case, we can calculate the
quasipotentials $V_{\phi_{i}}(\phi_{j})$ of every $\phi_{i}$ and $\phi_{j}$
with $i\not=j$ such that these points have adjacent basins: these are defined
as $V_{\phi_{i}}(\phi_{j})=\inf_{T>0}\inf_{\phi}S_{T}[\phi]$ for all paths
$\\{\phi(t)\\}_{T\in[0,T]}$ such that $\phi(0)=\phi_{i}$, $\phi(T)=\phi_{j}$,
and $\phi(t)\in\overline{B_{i}\cup B_{j}}$ for all $t\in[0,T]$. Consistently
we also set $V_{\phi_{i}}(\phi_{j})=+\infty$ if $B_{i}$ and $B_{j}$ have no
common boundary. These quasipotentials quantify the cost of the fluctuations
needed to escape $\phi_{i}$ conditional on entering $\phi_{j}$ next, and these
costs can be used to deduce the asymptotic rate of these transition events.
More precisely, as the parameter $\epsilon$ measuring the amplitude of the
macroscopic fluctuations tends to zero, the system dynamics can be
approximated by a Markov jump process (MJP) between the metastable states
$\phi_{1}$, …, $\phi_{N}$, with rates given asymptotically by

$k_{i,j}=e^{-V_{\phi_{i}}(\phi_{j})/\epsilon}\qquad
i,j=1,\ldots,N\quad(i\not=j).$ (36)

Questions about the asymptotic behavior of the system’s dynamics can be
answered by analyzing this MJP. Since the rates are all vanishing
exponentially at different rates as $\epsilon\to 0$, this process is quite
singular and can be analyzed by the method of decomposition into cycles
developed by Freidlin and Wentzell [39] (see also [40]). This method is
however rather intricate, and in practice it is often simpler to solve
specific questions in the MJP directly (e.g. what is its invariant
distribution or what is the mean first passage time from state $i$ to state
$j$), and take the limit as $\epsilon\to 0$ afterward. Note also that the
quasipotentials $V_{\phi_{i}}(\phi)$ for $i=1,\ldots,N$ can be used to
construct a global nonequilibrium potential $\phi$, solution of the Hamilton-
Jacobi equation (12) [39, 40]; since this construction is not doable in
practice for the high-dimensional examples we are interested in, we will not
dwell upon it here. We will however look at the minimizers of $\inf_{T\geq
0}\inf S_{T}[\phi]$, which can exist if the path $\\{\phi(t)\\}_{t\in[0,T]}$
is reparametrized using arclength rather than physical on any minimizing
sequence: these geometric paths give the mechanism of the transition between
$\phi_{i}$ and $\phi_{j}$.

The statements made in this section summarize what can be deduced from the
minimum action principle when it applies. As stated early, these results can
be proven rigorously in specific cases using e.g. tools for large deviation
theory (LDT), in the small noise [39] and weak interaction limit [52], or in
more complicated setups involving hydrodynamic limits [57, 58]. We can however
envision situations beyond the realm of LDT where min-max problem like (7)
with a Hamiltonian $H$ satisfying Assumptions A1, A2, and A3 arise (e.g. from
the Martin-Siggia-Rose-Jensen-De Dominicis [27, *janssen1976,
*de_dominicis1976] or the Doi-Peliti formalism [30, *peliti1985], or from
macroscopic fluctuation theory [41, 59]), and can be interpreted as above. The
methods that we propose below are of general purpose and were designed with
such general situations in mind.

##  III Computational Aspects

###  III.1 Min-max on the action as a wave propagation problem

Algorithm 1 : Action Minimization by Gradient Descent Ascent

1:Inputs: $M\in\mathbb{N}$; a path $\\{\phi^{0}_{i}\\}_{i\in I}$ with
$\phi_{0}^{0}=\phi_{a}$ and $\phi_{M}^{0}=\phi_{b}$; the functions $f(u,v)$
and $g(u,v)$; $T>0$, $\Delta\tau>0$, $\alpha>0$.

2:Initialization: For every $i\in I$, take $\theta^{0}_{i}=0$, and set
$u_{i}^{0}=\phi_{i}^{0}+\alpha\theta_{i}^{0}$ and
$v_{i}^{0}=\phi_{i}^{0}-\alpha\theta_{i}^{0}$; set $\Delta t=T/M$.

3:for $n\geq 0$ do

4: Update $u$ with an implicit upwind scheme, namely, solve
$\\{u^{n+1}_{i}\\}_{i\in I}$ sequentially from $i=M$ to $i=0$ using:
$\displaystyle\begin{cases}u^{n+1}_{M}=-v_{M}^{n}+2\phi_{b}\\\\[4.0pt]
\dfrac{u_{i}^{n+1}-u_{i}^{n}}{\Delta\tau}=\dfrac{u_{i+1}^{n+1}-u_{i}^{n+1}}{\Delta
t}+f(u^{n}_{i+1},v^{n}_{i+1}),\qquad i=M-1,\dots,0\end{cases}$

5: Update $v$ with an implicit upwind scheme, namely, solve
$\\{v_{i}^{n+1}\\}_{i\in I}$ sequentially from $i=0$ to $i=M$ using:
$\displaystyle\begin{cases}v^{n+1}_{0}=-u_{0}^{n+1}+2\phi_{a}\\\\[4.0pt]
\dfrac{v_{i}^{n+1}-v_{i}^{n}}{\Delta\tau}=-\dfrac{v_{i}^{n+1}-v_{i-1}^{n+1}}{\Delta
t}+g(u_{i-1}^{n+1},v_{i-1}^{n}),\qquad i=1,\dots,M\end{cases}$

6: Compute $\\{\phi_{i}^{n+1}=\frac{1}{2}(u_{i}^{n+1}+v_{i}^{n+1})\\}_{i\in
I}$ and
$\\{\theta_{i}^{n+1}=\frac{1}{2}\alpha^{-1}(u_{i}^{n+1}-v_{i}^{n+1})\\}_{i\in
I}$ (if needed).

In this section, we discuss how to solve the min-max problem stated in (7),
i.e.

$\inf_{\phi}S_{T}[\phi]=\min_{\phi}\max_{\theta}I_{T}(\phi,\theta)$ (37)

where we defined the functional

$I_{T}(\phi,\theta)=\int_{0}^{T}\left(\langle\dot{\phi},\theta\rangle-H(\phi,\theta)\right)dt$
(38)

and the optimization is to be performed over trajectories
$\\{\phi(t),\theta(t)\\}_{t\in[0,T]}$ subject to the boundary conditions
$\phi(0)=\phi_{a}$, $\phi(T)=\phi_{b}$. If the functional $I_{T}$ is convex
with respect to $\phi$ and concave with respect to $\theta$, then it is well-
known [60, 61] that this min-max problem can be solved by amortizing the
minimization and maximization over small alternating steps of steepest descent
in $\phi$ and steepest ascent in $\theta$. If these steps are infinitesimal in
some artificial optimization time $\tau$, this gradient descent ascent (GDA)
method leads to the evolution equation

$\displaystyle\partial_{\tau}\phi=-\alpha\delta
I_{T}/\delta\phi,\qquad\alpha\partial_{\tau}\theta=\delta I_{T}/\delta\theta,$
(39)

where for convenience we have introduced a parameter $\alpha>0$ that sets the
relative time scales over which $\phi$ and $\theta$ evolve – in the jargon of
GDA this is referred to two time scale GDA [62]. Calculating the functional
derivatives, the system (39) is explicitly given by

$\displaystyle\displaystyle\partial_{\tau}\phi=\alpha\partial_{t}\theta+\alpha\partial_{\phi}H,\qquad\displaystyle\alpha\partial_{\tau}\theta=\partial_{t}\phi-\partial_{\theta}H.$
(40)

These equations for $\\{\phi(\tau,t),\theta(\tau,t)\\}$ are to be solved with
the boundary conditions (in physical time $t$)

$\displaystyle\phi(\tau,t=0)=\phi_{a},\qquad\phi(\tau,t=T)=\phi_{b}.$ (41)

for some initial conditions (in optimization time $\tau$)

$\displaystyle\phi(\tau=0,t)=\phi^{0}(t),\qquad\theta(\tau=0,t)=\theta^{0}(t),$
(42)

with $\phi^{0}(t)$ such that $\phi^{0}(0)=\phi_{a}$ and
$\phi^{0}(T)=\phi_{b}$.

It is easy to see that the fixed points (in $\tau$) of (40) are solution to
Hamilton’s equations (8) that satisfy $\phi(0)=\phi_{a}$ and
$\phi(T)=\phi_{b}$. In Appendix A we show that: (i) there is a one-to-one
correspondence between the fixed points of (40) and the critical points of the
action $S_{T}[\phi]$, and (ii) if $\alpha$ is small enough these fixed points
are stable if and only if they are local minimizers of the action. Thus
solving (40) is indeed a way to perform $\inf_{\phi}S_{T}[\phi]$. For
illustrative purposes, we derive in Appendix B how the GDA converges to the
instanton for an Ornstein-Uhlenbeck process.

Let us now show how to put (40) in a form that is convenient for numerical
integration. Since (40) is an hyperbolic system of partial differential
equations (PDEs), it is useful to diagonalize the problem and introduce the
fields $u=\phi+\alpha\theta$ and $v=\phi-\alpha\theta$ that propagate along
characteristics and verify

$\displaystyle\partial_{\tau}u$ $\displaystyle=\partial_{t}u+f(u,v)$ (43)
$\displaystyle\partial_{\tau}v$ $\displaystyle=-\partial_{t}v+g(u,v),$ (44)

where we have defined

$\displaystyle f(u,v)$
$\displaystyle=\alpha\partial_{\phi}H-\partial_{\theta}H$ (45) $\displaystyle
g(u,v)$ $\displaystyle=\alpha\partial_{\phi}H+\partial_{\theta}H.$ (46)

The boundary conditions now only involve the propagating fields

$\displaystyle v(\tau,t=0)=-u(\tau,t=0)+2\phi_{a}$ (47) $\displaystyle
u(\tau,t=T)=-v(\tau,t=T)+2\phi_{b}.$ (48)

This formulation shows that the system made of (43) and (44) is well-posed
under these boundary conditions (see [63]) since the fields $v$ and $u$
propagate respectively forward and backward in physical time $t$ as the
optimization time $\tau$ increases. It also immediately suggests an algorithm
to solve Eqs. (43) and (44) based on Strang splitting [64]: to update the
fields at every iteration step in $\tau$, first update $u$ at $v$ fixed by
propagating the final condition at $t=T$ for $u$ in (48) towards $t=0$ using
(43) with forward differentiation in $t$, then update $v$ at $u$ fixed by
propagating the initial condition at $t=0$ for $v$ in (47) towards $t=T$ using
(44) with backward differentiation in $t$. Of course, the convergence does not
change if the algorithm starts by updating $v$ before updating $u$.

In practice, the continuous paths $\phi(\tau,t)$ and $\theta(\tau,t)$ are
discretized in physical time on $M+1$ points with index $i\in
I=\\{0,\cdots,M\\}$ such that $T=M\Delta t$, and we use index
$n\in\mathbb{N}_{0}$ to encode the evolution of the paths in optimization time
$\tau$ using steps of size $\Delta\tau$, so that for any field $\psi(\tau,t)$,
$\psi^{n}_{i}\equiv\psi^{n}(i\Delta t)$. The details are given in Algorithm 1.
The stability of the code relies on two important features. First, advection
of the fields is treated with an implicit upwind scheme for both $v$ and $u$.
Second, reaction terms $g$ and $f$ are also evaluated on an upwind grid point
with respect to the direction of advection. Reaction terms could also be
evaluated on site $i$ but we empirically found that the upwind implementation
strongly stabilizes the code when dealing with spatially extended diffusive
fields. The stability analysis of the numerical scheme is detailed in Appendix
C.

We should also emphasize that our interest resides in the fixed point of the
dynamics that solves Hamilton’s equation in physical time $t$, rather than in
the details of the dynamics in algorithmic time $\tau$, which has no physical
relevance. This consideration enjoins us to look for the largest time-step
$\Delta\tau$ that still provides a converging algorithm. The time-step $\Delta
t$, however, crucially needs to remains smaller than some characteristic time
$t_{c}$ needed to correctly resolve the dynamics of the instanton. This issue
is discussed in more details in Appendix C.

Finally, it is important to mention that a higher-order finite difference
stencil for the advection of the fields can be implemented while keeping the
same algorithmic complexity. Such a scheme can significantly improve
computation time since we need a smaller number of grid points to get the same
accuracy as the first-order scheme, introduced in the text for purpose of
simplicity. The second-order scheme is presented in Appendix D.

###  III.2 Geometric formulation on unbounded time intervals

Algorithm 2 : Geometric Action Minimization by Gradient Descent Ascent

1:Inputs: $M\in\mathbb{N}$; two stable fixed points $\phi_{a}$ and $\phi_{b}$
of the noiseless dynamics where
$H(\phi_{a,b},0)=\partial_{\theta}H(\phi_{a,b},0)=0$; a path
$\\{\hat{\phi}^{0}_{i}\\}_{i\in I}$ with $\hat{\phi}_{0}^{0}=\phi_{a}$ and
$\hat{\phi}_{M}^{0}=\phi_{b}$, such that
$|\hat{\phi}^{0}_{i+1}-\hat{\phi}^{0}_{i}|$ is constant in $i$; the functions
$f(u,v)$, $g(u,v)$, and $\lambda(u,v)$; $\Delta\tau>0$, $\alpha>0$.

2:Initialization: For every $i\in I$, take $\hat{\theta}^{0}_{i}=0$, and set
$u_{i}^{0}=\hat{\phi}_{i}^{0}+\alpha\hat{\theta}_{i}^{0}$ and
$v_{i}^{0}=\hat{\phi}_{i}^{0}-\alpha\hat{\theta}_{i}^{0}$; set $\Delta s=1/M$.

3:for $n\geq 0$ do

4: Update $u$ with an implicit upwind scheme, namely, solve
$\\{u^{n+1}_{i}\\}_{i\in I}$ sequentially from $i=M$ to $i=0$ using:
$\displaystyle\begin{cases}u^{n+1}_{M}=-v_{M}^{n}+2\phi_{b}\\\\[4.0pt]
\dfrac{u_{i}^{n+1}-u_{i}^{n}}{\Delta\tau}=\lambda_{i}(u^{n},v^{n})\dfrac{u_{i+1}^{n+1}-u_{i}^{n+1}}{\Delta
s}+f(u_{i}^{n},v_{i}^{n}),\qquad i=M-1,\dots,0\end{cases}$

5: Update $v$ with an implicit upwind scheme, namely, solve
$\\{v^{n+1}_{i}\\}_{i\in I}$ sequentially from $i=0$ to $i=M$ using:
$\displaystyle\begin{cases}v^{n+1}_{0}=-u_{0}^{n+1}+2\phi_{a}\\\\[4.0pt]
\dfrac{v_{i}^{n+1}-v_{i}^{n}}{\Delta\tau}=-\lambda_{i}(u^{n+1},v^{n})\dfrac{v_{i}^{n+1}-v_{i-1}^{n+1}}{\Delta
s}+g(u_{i}^{n+1},v_{i}^{n}),\qquad i=1,\ldots,M\\\ \end{cases}$

6: Compute $\\{\bar{\phi}^{n+1}=\frac{1}{2}(u^{n+1}+v^{n+1})\\}_{i\in I}$ and
$\\{\bar{\theta}^{n+1}=\frac{1}{2}\alpha^{-1}(u^{n+1}-v^{n+1})\\}_{i\in I}$.

7: Interpolate $\\{(\bar{\phi}_{i}^{n+1},\bar{\theta}_{i}^{n+1})\\}_{i\in I}$
onto a path $\\{(\hat{\phi}_{i}^{n+1},\hat{\theta}_{i}^{n+1})\\}_{i\in I}$
such that $|\hat{\phi}^{n+1}_{i+1}-\hat{\phi}^{n+1}_{i}|$ is constant in $i$,
as in the string method.

8: Set $\\{u^{n+1}=\hat{\phi}^{n+1}+\alpha\hat{\theta}^{n+1}\\}_{i\in I}$ and
$\\{v_{i}^{n+1}=\hat{\phi}_{i}^{n+1}-\alpha\hat{\theta}_{i}^{n+1}\\}_{i\in
I}$.

Let us now turn to the min-max problem

$V_{\phi_{a}}(\phi_{b})=\inf_{T>0}\inf_{\phi}S_{T}[\phi]=\inf_{T>0}\inf_{\phi}\sup_{\theta}I_{T}(\phi,\theta)$
(49)

where $I_{T}(\phi,\theta)$ is the functional defined in (38) and the inner
min-max is again to be performed over trajectories
$\\{\phi(t),\theta(t)\\}_{t\in[0,T]}$ subject to the boundary conditions
$\phi(0)=\phi_{a}$, $\phi(T)=\phi_{b}$.

In general a trajectory $\\{\phi(t),\theta(t)\\}_{t\in[0,T]}$ of finite time
duration $T$ cannot solve (49), i.e. we can always reduce the value of the
inner min-max by increasing $T$. This lack of optimizer complicates the
solution of (49). To proceed, it is useful to follow the strategy of the
geometric minimum action method (GMAM) in [32], and parametrize the physical
time as $t(s)$ for $s\in[0,1]$. Writing $dt/ds=\lambda^{-1}(s)$, and
introducing $(\hat{\phi}(s),\hat{\theta}(s))=(\phi(t(s)),\theta(t(s))$, the
minimization over $T$ in (49) can be turned into a minimization over
$\lambda$:

$V_{\phi_{a}}(\phi_{b})=\min_{\lambda\geq
0}\min_{\hat{\phi}}\max_{\hat{\theta}}\int_{0}^{1}\left(\langle\hat{\phi}^{\prime},\hat{\theta}\rangle-\lambda^{-1}H(\hat{\phi},\hat{\theta})\right)ds.$
(50)

where $\hat{\phi}^{\prime}=d\hat{\phi}/ds$ and the min-max is to be performed
over the triple $\\{\hat{\phi}(s),\hat{\theta}(s),\lambda(s)\\}_{s\in[0,1]}$
subject to the boundary conditions $\hat{\phi}(0)=\phi_{a}$,
$\hat{\phi}(1)=\phi_{b}$. Since we have added degrees of freedom by
representing $T$ by the function $\lambda(s)$, we can add a constraint on the
parametrization of $\hat{\phi}(s)$, e.g. by imposing that
$|\hat{\phi}^{\prime}|=cst$ (in $s$), in which case $s$ is normalized
arclength along the path $\hat{\phi}$. This choice is convenient numerically
as it guarantees that discretization points in $s$ will be uniformly
distributed along the path.

The min-max problem (50) is now in a form that can be solved using GDA,
similar to what we did for (37). As shown in Appendix E, it is however
convenient to treat $\lambda$ separately, and this leads to the following GDA
equations (compare (40))

$\displaystyle\partial_{\tau}\hat{\phi}=\alpha\lambda\partial_{s}\hat{\theta}+\alpha\partial_{\hat{\phi}}H,\qquad\alpha\partial_{\tau}\hat{\theta}=\lambda\partial_{s}\hat{\phi}-\partial_{\hat{\theta}}H,$
(51)

with $\lambda$ given by

$\displaystyle\lambda=\frac{|\langle\partial_{\theta}H,\partial_{s}\hat{\phi}\rangle|+\sqrt{\max(0,\langle\partial_{\theta}H,\partial_{s}\hat{\phi}\rangle^{2}-4H|\partial_{s}\hat{\phi}|^{2})}}{2|\partial_{s}\hat{\phi}|^{2}}.$
(52)

Eqs. (51) are to be solved under the constraint that
$|\partial_{s}\hat{\phi}|=cst$ (in $s$, not $\tau$), with the boundary
conditions (in $s$)

$\displaystyle\hat{\phi}(\tau,s=0)=\phi_{a},\qquad\hat{\phi}(\tau,s=1)=\phi_{b},$
(53)

for some initial conditions (in optimization time $\tau$)

$\displaystyle\hat{\phi}(\tau=0,s)=\hat{\phi}^{0}(s),\qquad\hat{\theta}(\tau=0,s)=\hat{\theta}^{0}(s),$
(54)

with $\hat{\phi}^{0}(s)$ such that $\hat{\phi}^{0}(0)=\phi_{a}$ and
$\hat{\phi}^{0}(1)=\phi_{b}$.

Similar to what we did with (40), the system of equations (51) as well as the
expression (52) for $\lambda$ can be put in a form suitable for numerical
solution by rewriting them in terms of the fields
$u=\hat{\phi}+\alpha\hat{\theta}$ and $v=\hat{\phi}-\alpha\hat{\theta}$. For
brevity we will not write these equations explicitly and refer the reader to
Algorithm 2 for their discretized version. At convergence, $\hat{\phi}$ and
$\hat{\theta}$ satisfy

$\displaystyle\lambda\hat{\phi}^{\prime}=\partial_{\theta}H,\qquad\lambda\hat{\theta}^{\prime}=-\partial_{\phi}H.$
(55)

If we insert the first of these equation in (52), we can reorganize this
equation into

$\displaystyle\lambda=\frac{\lambda|\hat{\phi}^{\prime}|^{2}+\sqrt{\lambda^{2}|\hat{\phi}^{\prime}|^{4}-4H|\hat{\phi}^{\prime}|^{2}}}{2|\hat{\phi}^{\prime}|^{2}}.$
(56)

This equation shows that at convergence we must also have

$\forall s\in[0,1]\ :\ H(\hat{\phi}(s),\hat{\theta}(s))=0.$ (57)

This is consistent with the fact that the term involving $\lambda^{-1}H$ in
(50) can also be interpreted as a Lagrangian multiplier term added to the
objective function to enforce the constraint that $H=0$, see Ref. [32] for
details.

In practice the continuous paths $\hat{\phi}(\tau,s)$ and
$\hat{\theta}(\tau,s)$ are discretized in $s$ on a grid of $M+1$ points with
index $i\in I=\\{0,\cdots,M\\}$ such that $1=M\Delta s$, and we use the index
$n\in\mathbb{N}_{0}$ to encode their evolution with step size $\Delta\tau$ in
artificial time $\tau$ so that for any field $\psi(\tau,s)$,
$\psi^{n}_{i}\equiv\psi^{n}(i\Delta s)$. The details are given in Algorithm 2,
which is a modified version of Algorithm 1 that includes the step of
reparametrization of the path. Note that we have kept the implicit upwind
discretization for the advection term. Also, the reaction terms $f$ and $g$,
and the coefficient $\lambda$ should be evaluated at the same grid point, as
originating from the same term $\lambda^{-1}H(\hat{\phi},\hat{\theta})$ in
(50). We found that evaluating these terms at the target grid point $i$ is
better for stability in general.

Algorithms 1 and 2 are simple to implement and only require evaluating the
first derivative of the Hamiltonian $H(\phi,\theta)$ with respect to its
arguments, just as we would have to do to solve Hamilton’s equations (8). It
was used to calculate the paths shown in Figs. 1, 2, and 3. In the more
complicated examples we treat below, because the state variable is a field
that depends on space as well as time, and as a result the PDEs (43) and (44)
involve spatial derivative too, some additional consideration must be given to
the way we discretize space and evaluate these derivatives to ensure that the
resulting scheme is numerically stable. As usual with PDEs, this issue needs
to be addressed on a case by case basis, depending of the nature of the PDE.

##  IV First-Order Phase Transitions in a Nonequilibrium Ginzburg-Landau (GL)
system

A possible starting point of the method borrows from Landau’s theory of phase
transitions where, relying on the symmetries of the system, we postulate the
free energy of the macroscopic variable of interest $\phi$ rather than
deriving it from a microscopic distribution [65, 66, 67, 68]. When the system
is in equilibrium, this procedure is well understood [69, 70]. For
nonequilibrium systems, macroscopic fluctuation theory (MFT) [71, 41, 59]
offers a generalization of Landau’s approach where we directly starts from the
action. In what follows, we apply a similar approach to study phase
transitions in a modified Ginzburg-Landau system.

###  IV.1 Nonequilibrium GL dynamics

Following the general framework introduced in [72, 73], the stochastic
evolution of a non-conservative field $\rho(t,x)$ can formally be described by
the Langevin equation

$\displaystyle\partial_{t}\rho=-\mu([\rho],x)+\sqrt{2\epsilon}\,\eta,$ (58)

where the drift $\mu([\rho],x)$ can be interpreted as a chemical potential,
$\eta(t,x)$ is a standard Gaussian white noise in space and time whose
amplitude is measured by $\epsilon>0$, and where for simplicity we have set
the mobility to 1. For a field $\rho$ defined on a one-dimensional domain, say
$[0,1]$, the action corresponding to this equation is

$\displaystyle
S_{T}[\rho]=\frac{1}{2}\int_{0}^{T}\int_{0}^{1}|\partial_{t}\rho+\mu([\rho],x)|^{2}dxdt,$
(59)

and it will be the subject of our investigations.

The dynamics (58) is in detailed balance when $\mu=\mu_{E}$, with $\mu_{E}$
such that it can be written as a derivative of a free energy functional
$\mathcal{F}[\rho]$,

$\mu_{E}([\rho],x)=\frac{\delta\mathcal{F}[\rho]}{\delta\rho(x)}.$ (60)

In this case, at equilibrium, the field configurations are distributed
according to the Gibbs measure associated with the free energy $\mathcal{F}$,
at temperature $\epsilon$.

Here we will be mostly interested in the active version of dynamics (58), when
$\mu=\mu_{E}+\mu_{A}$ with $\mu_{A}$ that cannot be cast into the form (60),
or, equivalently, does not satisfy the functional Schwarz relation [16, 12]:

$\displaystyle\frac{\delta\mu_{A}([\rho],x)}{\delta\rho(y)}-\frac{\delta\mu_{A}([\rho],y)}{\delta\rho(x)}\not=0.$
(61)

When $\mu_{A}\not=0$, the stationary distribution of the field configurations,
if it exists, is a nonequilibrium distribution which is not available in
closed form.

For concreteness, we will focus on the following example that displays a
nonequilibrium first-order phase transition: we assume that the field
$\rho(t,x)$ is one-dimensional, with $x\in[0,1]$ and periodic boundary
conditions, and we use chemical potentials given by

$\displaystyle\mu_{E}([\rho],x)=-(D\partial_{x}^{2}\rho+\rho-\rho^{3}+h)$ (62)
$\displaystyle\mu_{A}([\rho])=-\kappa\int_{0}^{1}\rho^{2}(y)dy,$ (63)

such that (58) becomes

$\displaystyle\partial_{t}\rho=D\partial^{2}_{x}\rho+\rho-\rho^{3}+h+\kappa\int_{0}^{1}\rho^{2}(y)dy+\sqrt{2\epsilon}\,\eta.$
(64)

Here $D>0$ is a diffusion constant which effectively depends on the system
size in the dimensionless variables we use (reducing $D$ is equivalent to
enlarging the domain size), $h$ the strength of an externally applied field,
and $\kappa$ the strength of the nonequilibrium coupling. The chemical
potential $\mu_{E}$ is borrowed from the Ginzburg-Landau $\phi^{4}$ theory,
often referred to as Model A [69]. A direct functional integration shows that
$\mu_{E}=\delta\mathcal{F}/\delta\rho$ with

$\mathcal{F}[\rho]=\int_{0}^{1}\left(\tfrac{1}{2}D|\partial_{x}\rho|^{2}+\tfrac{1}{4}(1-\rho(x))^{2}-h\rho(x)\right)dx.$
(65)

The chemical potential $\mu_{A}$ was chosen because the phenomenology that it
brings remains straightforward without being trivial. Specifically, $\mu_{A}$
acts as an additional uniform applied field, that is active and depends
nonlinearly on the value of $\rho$, rather than being externally imposed.
While $\mu_{E}$ drives the system towards the minimizers of the energy (65),
which are homogeneous state solutions of $\rho-\rho^{3}+h$, $\mu_{A}$
homogeneously pushes the field upward when $\kappa>0$ and downward when
$\kappa<0$. Therefore, the applied field $h$ and $\mu_{A}$ have competing
effects when $h$ and $\kappa$ have opposite signs. Since there is a region in
the $(\kappa,h)$ space where the noiseless dynamics has two stable fixed
points (see Sec. IV.2), this means that the system can undergo a
nonequilibrium first-order phase transition for critical values of $h$ and
$\kappa$ which we will determine in Sec. IV.3.

Figure 5: Phase diagram of the modified Ginzburg-Landau system. The phase
diagram is divided in 4 regions. Region (i): $\rho_{-}$ is the only phase.
Region (ii): coexistence region where the noiseless dynamics has two stable
homogeneous fixed points $\rho_{\pm}$, and one unstable homogeneous fixed
point $\rho_{c}$. Region (ii)a: $\rho_{-}$ is the stable phase. Region (ii)b:
$\rho_{+}$ is the stable phase. Region (iii): $\rho_{+}$ is the only phase.
Yellow solid line: $h_{c}^{+}(\kappa)$. Purple solid line:
$h_{c}^{-}(\kappa)$. Purple dashed line: phase transition curve between region
(ii)a and (ii)b obtained by the minimum action method. Grey dashed line: phase
transition curve obtained by treating the transition as if the system were in
equilibrium (wrong prediction).

a)

b)

c)

Figure 6: Nonequilibrium GL system for $D=5\times 10^{-3}$, $\kappa=1$,
$h=-0.5$, hence $\rho_{+}=1.45161$ and $\rho_{-}=-0.854638$. For these values,
we have $V_{\rho_{-}}(\rho_{+})<V_{\rho_{+}}(\rho_{-})$, indicating that
$\rho_{+}$ is the stable phase. Upper panels: contourplots of the paths
$\hat{\rho}(s,x)$; middle panels: contourplots of the conjugate momentum
$\hat{\theta}(s,x)$; lower panels: Lagrangian along the paths. Panels a):
minimum action path from $\rho_{-}$ to $\rho_{+}$; panels b): minimum action
path from $\rho_{+}$ to $\rho_{-}$. Notice the strong difference between the
forward and backward paths. The black dashed line at $s=0.5$ marks the
critical nucleus, past which $\theta$ and $L$ are both zero as they should.
Panel c): heteroclinic orbit joining $\rho_{-}$ and $\rho_{+}$ calculated by
the string method. Since the ascent is set as the reverse descent here, we
plot on one graph the conjugate momentum and the Lagrangian in both the
forward and reversed directions.

###  IV.2 Phase boundaries for coexisting homogeneous fixed points

Numerical evidence indicates that the stable fixed points of the noiseless
dynamics (i.e. Eq. (64) with $\epsilon=0$) are homogeneous states. As a
result, they are solutions to $\rho-\rho^{3}+h+\kappa\rho^{2}=0$. In the
domain where this equation has three real roots, $\rho_{-}$, $\rho_{c}$ and
$\rho_{+}$ with $\rho_{-}<\rho_{c}<\rho_{+}$, $\rho_{-}$ and $\rho_{+}$ are
stable fixed points of the noiseless dynamics whereas $\rho_{c}$ is an
unstable point. The coexistence region in the parameter space $(\kappa,h)$
where both $\rho_{-}$ and $\rho_{+}$ are present is marked as region (ii) in
Fig. 5; it is where $h\in[h_{c}^{-},h_{c}^{+}]$, with $h_{c}^{-}$ and
$h_{c}^{+}$ given by

$\displaystyle
h_{c}^{-}=-\frac{1}{27}\left(\sqrt{\kappa^{2}+3}+\kappa\right)^{2}\left(2\sqrt{\kappa^{2}+3}-\kappa\right),$
(66)

shown as a blue line in Fig. 5, and

$\displaystyle
h_{c}^{+}=\frac{1}{27}\left(\kappa-\sqrt{\kappa^{2}+3}\right)^{2}\left(2\sqrt{\kappa^{2}+3}+\kappa\right),$
(67)

shown as a yellow line in Fig. 5. Exactly on these boundaries, only two real
roots coexist, and one state is thus marginally stable. In regions (i) and
(iii) only one stable state exists.

Our next goal will be to analyze the relative stability of $\rho_{-}$ and
$\rho_{+}$ under the effect of the noise, i.e. derive the phase diagram of the
system. Even though we lack a free energy that yields the stationary measure,
we expect $\rho_{+}$ to be the stable phase for $h,\kappa>0$ and $\rho_{-}$
for $h,\kappa<0$. However, when $h$ and $\kappa$ have opposite signs (and thus
opposite effects on $\rho$), determining the most likely phase becomes
nontrivial.

###  IV.3 Extracting the minimum action paths

To assess whether $\rho_{+}$ or $\rho_{-}$ is the stable phase in the
coexistence region, we will compute the difference of the minimal actions
$\Delta S\equiv V_{\rho_{-}}(\rho_{+})-V_{\rho_{+}}(\rho_{-})$, obtained from
the minimum action paths from $\rho_{-}$ to $\rho_{+}$ and vice-versa. The
Hamiltonian entering the action is the one associated with Eq. (58):

$\displaystyle
H(\rho,\theta)=\langle-\mu[\rho],\theta\rangle_{L^{2}}+\langle\theta,\theta\rangle_{L^{2}},$
(68)

where the scalar product of two functions $f$ and $g$ is given by $\langle
f,g\rangle_{L^{2}}=\int_{0}^{1}f(x)g(x)dx$. To obtain the phase diagram, these
calculations must be repeated for a set of values $(\kappa,h)$ in the
coexistence region to compute the minimal action difference as a function of
these parameters, $\Delta S(\kappa,h)$: From Eq. (5), the line of phase
transition is then the curve where $\Delta S(\kappa,h)=0$.

In practice we use Algorithm 2 with $M=400$ copies along the path, $N_{x}=64$
points of space discretization, $\Delta\tau=10^{-3}$, and $\alpha=0.33$, and
we monitor convergence by looking at the decay of the action. The result of
these computations is shown in Fig. 5, where the purple dashed line frontier
in the phase diagram corresponding to $\Delta S(\kappa,h)=0$. The stable phase
is $\rho_{-}$ below the line (region (ii)a) , and $\rho_{+}$ above it
(region(ii)b).

For comparison, we also compute the phase diagram under the (wrong) assumption
that the escape paths were given by the heteroclinic orbits followed in a
time-reversed way. This would have to be the case in equilibrium by time-
reversal symmetry. While these escape paths are incorrect in general in
nonequilibrium systems, their respective cost in the action gives an upper
bound on the actual minima $V_{\rho_{-}}(\rho_{+})$ and
$V_{\rho_{+}}(\rho_{-})$. Denoting by $S^{\mathrm{het}}_{\rho_{-}}(\rho_{+})$
and $S^{\mathrm{het}}_{\rho_{+}}(\rho_{-})$ the actions along the heteroclinic
orbit, we compute for instance $S^{\mathrm{het}}_{\rho_{-}}(\rho_{+})$ as

$\displaystyle
S^{\mathrm{het}}_{\rho_{-}}(\rho_{+})=\int_{0}^{1}\left(\langle\partial_{s}\hat{\rho},\hat{\theta}\rangle-\lambda^{-1}H(\hat{\rho},\hat{\theta})\right)ds,$
(69)

where the path $\hat{\rho}(s)$ and the parametrization $\lambda(s)$
($s\in[0,1]$) have been obtained by the string method [74, 75] that identifies
the heteroclinic orbit between $\hat{\rho}(0)=\rho_{-}$ and
$\hat{\rho}(1)=\rho_{+}$, and $\theta$ solves

$\displaystyle\lambda\partial_{s}\hat{\rho}=\partial_{\theta}H(\hat{\rho},\hat{\theta}).$
(70)

As a sanity check, we did verify that one always has
$V_{\rho_{-}}(\rho_{+})\leq S^{\mathrm{het}}_{\rho_{-}}(\rho_{+})$ and
$V_{\rho_{+}}(\rho_{-})\leq S^{\mathrm{het}}_{\rho_{+}}(\rho_{-})$, namely,
that the minimizer of the action is always smaller than the action along the
heteroclinic orbit. It is also worth noticing that these bonds offer no
information about the location of the phase transition line: The line where
$\Delta S^{\mathrm{het}}(\kappa,h)=0$ is plotted as the grey dashed line in
Fig. 5, and it is different from the actual transition line $\Delta
S(\kappa,h)=0$

The minimum action paths also give physical insights about the mechanism of
the transition: the contour plot in $(s,x)$ space of these paths are shown in
Fig. 6 for the specific value $(\kappa,h)=(1,-0.5)$ (which is in region
(ii)a): panel (a) shows the forward path from $\rho_{-}$ to $\rho_{+}$, and
panel (b) the reversed path from $\rho_{+}$ to $\rho_{-}$. Also shown in panel
(c) of the figure is the heteroclinic orbit. The actual path
$\\{\hat{\rho}(s)\\}_{s\in[0,1]}$ are shown in the first row, while the second
row displays the conjugate momentum $\\{\hat{\theta}(s)\\}_{s\in[0,1]}$, and
the third the action increment (i.e. the Lagrangian)
$\langle\hat{\theta}(s),\hat{\rho}^{\prime}(s)\rangle$ along the path.

As can be seen in Fig. 6 the forward and the backward minimum action paths are
different, and they cross the separatrix (marked as a dashed black line in the
figure) at different places: that is, the critical nucleii for the forward and
backward transitions are different (for the backward path this ‘critical
nucleus’ is actually flat). This is a signature of time-symmetry breaking that
can be intuitively explained as follows: When $\kappa>0$, as in Fig. 6, the
nonequilibrium term $\kappa\int\rho^{2}dx$ favors the movement from $\rho_{-}$
to $\rho_{+}$ but opposes the one from $\rho_{-}$ to $\rho_{+}$. In the
forward path from $\rho_{-}$ to $\rho_{+}$, it is better to have
$\int\rho^{2}dx$ large, which favors nucleation; conversely, in the backward
path from $\rho_{+}$ to $\rho_{-}$, it is better to have $\int\rho^{2}dx$
small, and the way to minimize this quantity given the value of its changing
mean is to have $\rho$ spatially uniform. Notice however that this is a finite
size effect: If $D$ were decreased to even smaller values, we would observe
nucleation events in both directions (albeit different ones in each).

Note also that the interesting part of these minimum action paths is their
escape half, where the noise is needed and the action increment is therefore
positive: it is the first half for the forward path in panel (a) and the
second half for the backward path in panel (b); for the heteroclinic orbit
shown in panel (c), we display both halves at once since the forward and
reversed paths are symmetric. In all situations, past the critical nucleus,
the paths simply follow the noiseless dynamics, and this half of the path may
not be unique if the critical nucleus has more than one unstable direction.
This non-uniqueness has no impact on the action however, since the Lagrangian
is zero along the solution of the noiseless dynamics.

The effective nonequilibrium Ginzburg-Landau-like dynamics we have considered
in this section can be modified to include different nonequilibrium terms,
like for example $\mu_{A}=\kappa|\partial_{x}\rho|^{2}$ which naturally
appears in the coarse-grained field description of interface growth phenomena
[76] and active matter systems [15, 72, 77, 78]. We could also use our
approach to compute the phase diagram of modified Cahn-Hilliard systems, which
naturally emerge in active matter field theories [11, 72, 73, 36, 78, 79].

##  V Phase transitions in a bistable reaction-diffusion system

###  V.1 The Schlögl model

In 1972, Schlögl introduced the following chemical reaction network [46]

$\displaystyle A\xrightleftharpoons[k_{1}]{k_{0}}X,\qquad
2X+B\xrightleftharpoons[6k_{3}]{2k_{2}}3X,$ (71)

with microscopic rates $k_{i}>0$, and where the concentration of $A$ and $B$
are held constant. In a certain regime of the reaction rates, this system
displays metastability between a low density and a high density phase. Here we
will consider a spatially extended variant of this model introduced by Tănase-
Nicola and Lubensky [48], in which a one dimensional domain is split into
$L\in\mathbb{N}$ well-stirred compartments: the molecules only react within
their compartment, and randomly jump to neighboring ones with rate $\gamma>0$.
We also impose periodic boundary conditions.

When the number $n_{i}$ of molecules in compartment $i\in\\{1,\dots,L\\}$ is
large, it is convenient to introduce the rescaled $\rho_{i}=n_{i}/\Omega$,
where $\Omega\gg 1$ is the typical number of molecules per compartments. In
the limit as $\Omega\to\infty$, the law of mass action for $\rho_{i}$ is a
(discrete) reaction-diffusion equation

$\dot{\rho}_{i}=\gamma(\rho_{i+1}+\rho_{i-1}-2\rho_{i})+w_{+}(\rho_{i})-w_{-}(\rho_{i})$
(72)

where $w_{+}(\rho_{i})=\lambda_{0}+\lambda_{2}\rho_{i}^{2}$ and
$w_{-}(\rho_{i})=\lambda_{1}\rho_{i}+\lambda_{3}\rho_{i}^{3}$ are the rescaled
reaction rates with $\lambda_{i}=k_{i}\Omega^{i-1}$ (see App. G). (72) can be
written as a gradient flow:

$\dot{\rho}_{i}=-\partial_{\rho_{i}}E(\rho),$ (73)

where $\rho=(\rho_{1},\dots,\rho_{L})$ and we introduced

$E(\rho)=\sum_{i=1}^{L}\left(\tfrac{1}{2}\gamma(\rho_{i+1}-\rho_{i})^{2}+U(\rho_{i})\right)$
(74)

with

$U(\rho_{i})=-\lambda_{0}\rho_{i}-\tfrac{1}{3}\lambda_{2}\rho_{i}^{3}+\tfrac{1}{2}\lambda_{1}\rho_{i}^{2}+\tfrac{1}{4}\lambda_{3}\rho_{i}^{4}.$
(75)

We recognize a (discrete) Ginzburg-Landau free energy, from which we conclude
that the stable fixed points of (72) are homogeneous states as long as
$\gamma$ is not too small. For the value of
$(\lambda_{0},\lambda_{1},\lambda_{2},\lambda_{3})$ that we will consider
here, there are two such fixed points, $\rho_{i}=\rho_{\pm}$ for all
$i\in\\{1,\dots,L\\}$, where $\rho_{-}<\rho_{+}$ are the largest and smallest
roots of $w_{+}(z)=w_{-}(z)$.

The gradient structure of (73) may suggest that $\rho_{-}$ is the stable phase
if $E(\rho_{-})<E(\rho_{+})$ whereas $\rho_{+}$ is if
$E(\rho_{-})>E(\rho_{+})$. This conclusion is however incorrect, as already
observed by Tănase-Nicola and Lubensky [48] who showed that the system
undergoes a nonequilibrium first-order phase transition when $\gamma$ changes.
Since $E(\rho_{\pm})=U(\rho_{\pm})$ and, therefore, is independent of
$\gamma$, the phase transition cannot be predicted by analyzing $E(\rho)$
only: the system is not in detailed balance with respect to the Gibbs measure
associated to $E(\rho)$. What was not determined in [48] is the critical value
$\gamma_{c}$ at which the phase transition occurs. This is the question we
address next.

Figure 7: Phase diagram of the Schlögl model in the limit of a large number of
particles jumping between well-stirred reactive boxes. The parameters
$\lambda_{i}$ are fixed and we vary the jump rate $\gamma$. For
$\gamma<\gamma_{c}=5\pm 0.1$, we have
$V_{\rho_{+}}(\rho_{-})>V_{\rho_{-}}(\rho_{+})$, which indicates that
$\rho_{+}$ is stable state phase, whereas $\rho_{-}$ is for
$\gamma>\gamma_{c}$. We use the parameters from [36]: $\lambda_{0}=0.8$,
$\lambda_{1}=2.9$, $\lambda_{2}=3.1$, $\lambda_{3}=1$, and took $L=40$.

a)

b)

c)

Figure 8:  Spatially extended Schlögl model for $L=40$ reactive boxes,
$(\lambda_{0},\lambda_{1},\lambda_{2},\lambda_{3})=(0.8,2.9,3.1,1)$, and
$\gamma=12$. For these values we have $\rho_{s}=1$ and $\rho_{+}=1.6$. The
figure is organized as Fig. 6. As spatial coordinate, we use $x=i/L$ and plot
the transitions paths and associated momenta as if they were continuous in
space.

###  V.2 Change of relative stability with increasing jumping rate

The spatially extended Schlögl model is a reaction network of the type
considered in Sec. II.1, and its phase diagram can be analyzed by minimizing
the action associated with a Hamiltonian similar to (22). Using the structure
of the model, it is natural to decompose $H=H^{R}+H^{D}$, with $H^{R}$
accounting for the reaction and $H^{D}$ for the jumps:

$\displaystyle H^{R}(\rho,\theta)$
$\displaystyle=\sum_{i=1}^{L}w_{+}(\rho_{i})(e^{\theta_{i}}-1)+w_{-}(\rho_{i})(e^{-\theta_{i}}-1),$
(76) $\displaystyle H^{D}(\rho,\theta)$
$\displaystyle=\gamma\sum_{i=1}^{L}\rho_{i}(e^{\theta_{i-1}-\theta_{i}}+e^{\theta_{i+1}-\theta_{i}}-2),$
(77)

We use this Hamiltonian in the action that we minimize using Algorithm 2 to
calculate the quasipotentials $V_{\rho_{-}}(\rho_{+})$ and
$V_{\rho_{-}}(\rho_{+})$. We repeat these calculations for different values of
the jump rate $\gamma$ while keeping the rates fixed at
$(\lambda_{0},\lambda_{1},\lambda_{2},\lambda_{3})=(0.8,2.9,3.1,1)$, for which
$\rho_{-}=0.5$ and $\rho_{+}=1.6$. We use $L=40$ compartments and in Algorithm
2 we set $M=400$ ($\Delta s=2.5\times 10^{-3}$), $\alpha=1$, and
$\Delta\tau=0.01$. The graphs of $V_{\rho_{-}}(\rho_{+})$ and
$V_{\rho_{-}}(\rho_{+})$ versus $\gamma$ are shown in Fig. 7. These results
indicate that the nonequilibrium first-order phase transition occurs at
$\gamma_{c}\simeq 5$: $\rho_{+}$ is the stable phase for $\gamma>\gamma_{c}$,
while $\rho_{-}$ is the stable one for $\gamma<\gamma_{c}$. We stress again
that this result cannot be deduced by looking at $E(\rho)$ even though the law
of mass action can be written as the gradient flow (73). This is of course not
a contradiction: the Schlögl model is not in detailed balance and lacks time
reversal symmetry because this property also depends on the nature of the
noise.

The non-equilibrium nature of the phase transition can be confirmed by looking
at the transition paths from $\rho_{-}$ and $\rho_{+}$ and vice-versa. They
are shown in Fig. 8, in which we use the same plotting conventions as in Fig.
6. These results are for $\gamma=12$, when $\rho_{-}$ is the stable phase. We
can see that the forward (panel a) and the backward (panel b) paths are
different and go through different critical nucleii (marked as a dashed
vertical black line on the graphs). These path are also different from the
heteroclinic orbit (panel c). If we increase the value of $\gamma$, the
forward path eventually become homogeneous too (results not show): this is
consistent with the fact that at high $\gamma$, the system behaves essentially
as one single well-stirred compartment. In that limit we can calculate the
nonequilibrium steady distribution of the system and use it to calculate
$V_{\rho_{-}}(\rho_{+})$ and $V_{\rho_{-}}(\rho_{+})$: this is done in App. G
and gives the same values as the ones obtained by Algorithm 2 when $\gamma$ is
large and the transition paths are both homogeneous. Conversely, if we
decrease the value of $\gamma$, the backward path becomes inhomogeneous too
(results not shown). This transition from homogeneous to inhomogeneous
backward path occur around $\gamma=7$. How to gain intuition about these
changes of behavior is harder on this model than in the GL model of Sec. IV
because the noise is non-Gaussian, and ultimately the shape of the minimum
action paths depend on a complex interplay between many effects in the
dynamics.

###  V.3 Comparison with microscopic simulations

Figure 9: Spatially extended Schlögl model: Relaxation of the microscopic
system with a large number of particles per site starting from a step profile
at $t=0$. The microscopic dynamics (dashed lines) closely follows the solution
of (72) (solid lines) with same initial condition, and converges to the
stationary fixed point $\rho=\rho^{+}$ within the same physical time $t\simeq
18$. Parameters: $\lambda_{0}=0.8$, $\lambda_{1}=2.9$, $\lambda_{2}=3.1$,
$\lambda_{3}=1$, $\gamma=4$, $L=40$, $\Omega=16\times 10^{4}$.  Figure 10:
Logarithm of the average nucleation time $\tau_{+,-}$ needed to reach the
separatrix starting from a state metastable state $\rho_{+}$ close to linear
instability, as a function of the typical number of particles $\Omega$ per
site. In this regime, large deviation theory can be tested with Monte Carlo
simulations in reasonable time. LDT does not predict the prefactor in front of
the exponential scaling. Blue dots: results of MC simulations. Dashed line:
slope predicted by theory $V_{\rho_{+}}(\rho_{-})=9.931\times 10^{-4}$.
Parameters: $\lambda_{0}=0.728$, $\lambda_{1}=2.9$, $\lambda_{2}=3.1$,
$\lambda_{3}=1$, $\gamma=20$, $L=40$.

In this section we corroborate the conclusions of Sec. V.2 by performing
Markov Chain Monte Carlo (MCMC) simulations of the microscopic system.

For the parameters value used in Sec. V.2 it is not possible to calculate this
way the phase diagram shown in Fig. 7: this is because $\Omega$ needs to be
large in order for the phases $\rho_{-}$ and $\rho_{+}$ to be (meta)stable
under the noise, and the timescales of the forward or backward transition
between the phases $\rho_{-}$ and $\rho_{+}$ is of the order of the nucleation
times $e^{\Omega V_{\rho_{-}}(\rho_{+})}$ and $e^{\Omega
V_{\rho_{+}}(\rho_{-})}$. These timescales are too big to be accessible with
MCMC. We can however perform two types of experiments to test the results of
the minimum action principle:

First we can check that MCMC simulations of the microscopic system initiated
with some nonuniform profile in the compartments behave as predicted by Eq.
(72). The results of these simulations are shown in Fig. 9 where we compare
the evolution of a step profile: the microscopic system follows the
deterministic dynamics (72), as expected. To do these calculation, we fix
$\Omega=1.6\times 10^{5}$, set the number of molecules in each compartment
$i\in\\{1,\dots,L\\}$ to be $n_{i}=\lfloor\rho_{i}\Omega\rfloor$ for all $i$,
and simulate the microscopic dynamics exactly with Gillespie algorithm [80,
81] using the microscopic rates $k_{i}=\lambda_{i}\Omega^{1-i}$.

Second, to check the prediction of the minimum action principle in terms of
nucleation times, we can change the rates $\lambda_{i}$ to make one of two
phases (say $\rho_{+}$) only very weakly metastable, i.e. such that
$V_{\rho_{+}}(\rho_{-})\ll 1$. We can then calculate the mean escape time
$\tau_{+,-}$ from this state towards $\rho_{-}$, repeat this calculation for
different values of $\Omega$ and check that $\tau_{+,-}\asymp e^{\Omega
V_{\rho_{+}}(\rho_{-})}$. The result is shown in Fig. 10, which confirms that
the prediction from the minimum action framework explains the microscopic
simulations.

###  V.4 Continuous limit

Finally, let us consider the continuous space limit of the model by sending
the number of compartments $L\to\infty$. To this end, let us set
$\rho_{i}=\tilde{\rho}(x_{i})/L$ and $\theta_{i}=\tilde{\theta}(x_{i})$, with
$x_{i}=i/L$ and where $\tilde{\rho}(x)$ and $\tilde{\theta}(x)$ are fields on
$x\in[0,1]$. Let us also set $\gamma=DL^{2}$ for some diffusion coefficient
$D>0$, and $\lambda_{i}=\tilde{\lambda}_{i}L^{i-1}$ for some rescaled rates
$\tilde{\lambda}_{i}$. Assuming that $D$ and $\tilde{\lambda}_{i}$ are $O(1)$
in $L$, in the limit as $L\to\infty$, it is easy to verify that the
Hamiltonians $H^{R}$ and $H^{D}$ in (76) and (77) now become

$\displaystyle\begin{split}\tilde{H}^{R}(\tilde{\rho},\tilde{\theta})&=\int_{0}^{1}\left(\tilde{w}_{+}(\tilde{\rho})(e^{\tilde{\theta}}-1)+\tilde{w}_{-}(\tilde{\rho})(e^{-{\tilde{\theta}}}-1)\right)dx\\\
\tilde{H}^{D}(\tilde{\rho},\tilde{\theta})&=D\int_{0}^{1}\left(\tilde{\rho}\partial^{2}_{x}\tilde{\theta}+\tilde{\rho}(\partial_{x}\tilde{\theta})^{2}\right)dx,\end{split}$
(78)

where we defined
$\tilde{w}_{+}(\tilde{\rho})=\tilde{\lambda}_{0}+\tilde{\lambda}_{2}\tilde{\rho}^{2}$,
$\tilde{w}_{-}(\tilde{\rho})=\tilde{\lambda}_{1}\tilde{\rho}+\tilde{\lambda}_{3}\tilde{\rho}^{3}$.
Note that in this limit the discrete Poisson jumps of the molecules are
approximated as a Gaussian noise on the density: we have recovered the
multiplicative Gaussian noise that appears in the Dean-Kawasaki equation [82,
51]. The structure of the Poisson noise of the reaction is however left
unaffected by the limit.

We used Algorithm 2 to calculate transition pathways at continuous level by
minimizing the action associated with $\tilde{H}=\tilde{H}^{R}+\tilde{H}^{D}$.
The pathways (not shown) are not significantly different from those shown in
Fig. 8 when $L=1/40$: this indicates that the system was already close to its
continuous limit at that value of $L$.

##  VI Conclusion

In summary, the analysis of first-order phase transitions and other activated
processes in nonequilibrium systems can be reduced to the minimization of an
action in situations where these processes are rare and occur via reproducible
pathways. The approach can be justified rigorously within the framework of
large deviation theory (LDT), and minimum action principles can also be
derived formally in other instances using e.g. the Martin-Siggia-Rose-Janssen-
De Dominicis [27, *janssen1976, *de_dominicis1976] or the Doi-Peliti [30,
*peliti1985] formalisms. The minimum of the action can be used to generalize
Arrhenius law, and its minimizer to explain the mechanism of the transitions,
including the shape of the critical nucleus that serves as transition state.
Concrete predictions however rest on our ability to solve this minimization
problem, which often needs to be done numerically.

Here we developed algorithms to perform these calculations, both in finite and
infinite times. These algorithms are designed to be used directly within the
Hamiltonian formulation of the action, which leads to a min-max problem, and
do not require the user to calculate the Lagrangian beforehand. In particular
it can be used for systems where the fluctuations are non-Gaussian and the
Lagrangian is typically unavailable in closed form. The applicability of our
method was tested on two nontrivial examples involving spatially extended
systems undergoing nonequilibrium phase transitions: a modified Ginzburg-
Landau equation perturbed by noise and a reaction-diffusion system based on
the Schlögl model. In both cases the method allowed us to calculate the phase
diagram of the systems, compute the paths of the transitions, and identify the
critical nucleus.

We hope that our work will pave the way for a systematic approach to study
activated processes in other interesting examples, like active matter systems
undergoing a motility-induced phase separation or a flocking transition, or
other systems which display phase transitions but are not in detailed balance.
These applications will depend on the possibility to derive an appropriate
minimum action principle, potentially through coarse-graining, which is a non-
trivial question on its own.

###### Acknowledgements.

We thank Jasna Brujic, Tobias Grafke, Tobias Kühn, and Frédéric van Wijland
for useful comments. This work was supported by the Materials Research Science
and Engineering Center (MRSEC) program of the National Science Foundation
under Grants No. NSF DMR-1420073 and in part by Grant No. NSF DMR-1710163. RZ
would like to thank Laboratoire MSC Paris for hospitality. RZ and EVE would
also like to thank Center for Data Science ENS Paris for hospitality.

##  Appendix A Convergence of the min-max for $\alpha$ small but finite

We would like to show that if a path is stable with respect to the Lagrangian
minimization, then the path is also stable with respect to the Hamiltonian
min-max algorithm, for $\alpha$ small enough. We start from (39), but in this
section, we conveniently rescale artificial time $\tau=\tilde{\tau}/\alpha$,
we set $\tilde{\alpha}=\alpha^{2}$ and drop the tilde such that the evolution
equations read

$\displaystyle\begin{cases}\partial_{\tau}x=\partial_{t}\theta+\partial_{x}H\\\
\alpha\partial_{\tau}\theta=\partial_{t}x-\partial_{\theta}H.\end{cases}$ (79)

We focus on the dynamical system subjected to an additive Gaussian white
noise, as presented in Eqs. (17) and (18), i.e the Hamiltonian takes the form

$\displaystyle H(x,\theta)=\langle
b(x),\theta\rangle+\frac{1}{2}|\theta|^{2}.$ (80)

We assume that a minimum action path $\\{(x^{*},\theta^{*})\\}_{t\in[0,T]}$
has been obtained. We look at a perturbed path $(x^{*}+X,\theta^{*}+\Theta)$
with $X=(x_{1},\cdots,x_{P})^{T}$ and
$\Theta=(\theta_{1},\cdots,\theta_{P})^{T}$ the perturbations, and we assess
the conditions for path relaxation to the minimum action path. The evolution
of the perturbation reads

$\displaystyle\begin{cases}\partial_{\tau}x_{p}=\dot{\theta}_{p}+\theta_{j}^{*}\dfrac{\partial^{2}b^{*}_{j}}{\partial
x_{k}\partial x_{p}}x_{k}+\theta_{j}\dfrac{\partial b^{*}_{j}}{\partial
x_{p}}\\\ \alpha\partial_{\tau}\theta_{p}=\dot{x}_{p}-\dfrac{\partial
b^{*}_{p}}{\partial x_{k}}x_{k}-\theta_{p},\end{cases}$ (81)

where we have used the Einstein convention for the sum on repeated indices,
and the shorthand notations $\frac{\partial b^{*}_{p}}{\partial
x_{k}}=\frac{\partial b_{p}}{\partial x_{k}}|_{x^{*}}$ and
$\frac{\partial^{2}b^{*}_{j}}{\partial x_{k}\partial
x_{p}}=\frac{\partial^{2}b_{j}}{\partial x_{k}\partial x_{p}}|_{x^{*}}$, for
any indices $j$, $k$, $p$.

The evolution of the fields can now be cast into the following form:

$\displaystyle\begin{cases}\partial_{\tau}X=\mathcal{L}_{1}X+\mathcal{L}_{2}\Theta\\\
\partial_{\tau}\Theta=\alpha^{-1}(\mathcal{L}_{3}X-\Theta)\end{cases}$ (82)

where

$\displaystyle(\mathcal{L}_{1})^{(pk)}$
$\displaystyle=\theta_{j}^{*}\dfrac{\partial^{2}b^{*}_{j}}{\partial
x_{k}\partial x_{p}}$ (83) $\displaystyle(\mathcal{L}_{2})^{(pk)}$
$\displaystyle=(\delta_{pk}\partial_{t}+\frac{\partial b_{k}^{*}}{\partial
x_{p}})$ (84) $\displaystyle(\mathcal{L}_{3})^{(pk)}$
$\displaystyle=(\delta_{pk}\partial_{t}-\frac{\partial b_{p}^{*}}{\partial
x_{k}})$ (85)

In the Lagrangian algorithm, the equation $\Theta=\mathcal{L}_{3}X$ is always
verified, so the evolution of the perturbation $X$ is simply given by

$\displaystyle\partial_{\tau}X=(\mathcal{L}_{1}+\mathcal{L}_{2}\mathcal{L}_{3})X.$
(86)

Any perturbation $X$ vanishes if the eigenvalues of the operator
$\mathcal{L}=\mathcal{L}_{1}+\mathcal{L}_{2}\mathcal{L}_{3}$ are all of
negative real part. The operator $\mathcal{L}$ is self-adjoint (since
$\mathcal{L}$ is the Hessian of the action $S$), and one can thus extract a
basis of normalized orthogonal eigenvectors $X_{n}^{\mathcal{L}}$. Let us
denote by $\mu_{n}$ the n-th eigenvalue and $X_{n}^{\mathcal{L}}$ the
corresponding eigenvector. We have
$\mathcal{L}X_{n}^{\mathcal{L}}=\mu_{n}X_{n}^{\mathcal{L}}$, with $\mu_{n}<0$.

Now, in the Hamiltonian algorithm, we would like to find the conditions under
which any perturbation $(X,\Theta)$ close to a path of minimum action vanishes
when evolving in artificial time $\tau$. The perturbation will vanish if and
only if the eigenvalues $\lambda_{n}$ of the linear operator given in Eq. (82)
have a negative real part. The eigenvalue $\lambda_{n}$ associated to the
eigenvector $(\Theta_{n},X_{n})^{T}$ should verify

$\displaystyle\begin{cases}\mathcal{L}_{1}X_{n}+\mathcal{L}_{2}\Theta_{n}=\lambda_{n}X_{n}\\\
\alpha^{-1}(\mathcal{L}_{3}X_{n}-\Theta_{n})=\lambda_{n}\Theta_{n}.\end{cases}$
(87)

The second equation in (87) yields
$\Theta_{n}=(1+\lambda_{n}\alpha)^{-1}\mathcal{L}_{3}X_{n}$, assuming that
there exists $\alpha>0$ such that $(1+\lambda_{n}\alpha)\neq 0$ for every $n$.
Injecting this result into the first equation of (87) yields a closed equation
for $X_{n}$ and $\lambda_{n}$

$\displaystyle(1+\lambda_{n}\alpha)(\mathcal{L}_{1}X_{n}-\lambda_{n}X_{n})+\mathcal{L}_{2}\mathcal{L}_{3}X_{n}=0.$
(88)

For $\alpha\ll 1$, the system (82) can be seen as a perturbation of the
Lagrangian problem, where an additional degree of freedom $\Theta$ relaxes to
$\mathcal{L}_{3}X$ on a fast time scale $1/\alpha$. This suggests to look for
eigenvalues with a specific form: (i) a first set of eigenvalues
$\lambda_{n}^{(1)}$ should be the perturbed eigenvalues $\mu_{n}$ with
perturbed $X_{n}^{\mathcal{L}}$ as corresponding eigenvectors, (ii) a second
set of eigenvalues $\lambda_{n}^{(2)}$ is expected to scale as
$O(\alpha^{-1})$ and encodes the fast relaxation of the variable $\Theta$ to
$\mathcal{L}_{3}X$.

Therefore we look for eigenvalues of a general form

$\lambda_{n}=\frac{\eta_{n}}{\alpha}+\nu_{n}+O(\alpha)$ (89)

with $\eta_{n}$ and $\nu_{n}$ of $O(1)$. Expanding at leading order
$O(\alpha^{-1})$ in (88) yields

$\displaystyle\eta_{n}(1+\eta_{n})X_{n}=0,$ (90)

implying that $\eta_{n}=0$ or $\eta_{n}=-1$. The case $\eta_{n}=0$ leads us to
consider the first case (i) mentioned above, in which we look for eigenvalues
$\lambda_{n}^{(1)}=\mu_{n}+\alpha\xi_{n}$, with $\xi_{n}=O(1)$ and we expand
(88) to leading order. At order $0$ in $\alpha$, we find that $X_{n}$ must
solve

$\displaystyle(\mathcal{L}_{1}+\mathcal{L}_{2}\mathcal{L}_{3})X_{n}-\mu_{n}X_{n}=0,$
(91)

where we recognize the operator
$\mathcal{L}=\mathcal{L}_{1}+\mathcal{L}_{2}\mathcal{L}_{3}$, which confirms
that we expand around the eigenvalues $X_{n}^{\mathcal{L}}$ of the Lagrangian
system. We write $X_{n}=X_{n}^{\mathcal{L}}+\alpha Y_{n}$ with
$\|Y_{n}\|=O(1)$. Now expanding (88) at order 1 in $\alpha$, and using the
relation $\mathcal{L}X_{n}^{\mathcal{L}}=\mu_{n}X_{n}^{\mathcal{L}}$, we get

$\displaystyle\mathcal{L}Y_{n}-\mu_{n}Y_{n}=\xi_{n}X_{n}^{\mathcal{L}}+\mu_{n}\mathcal{L}_{2}\mathcal{L}_{3}X_{n}^{\mathcal{L}}.$
(92)

By taking the scalar product with the eigenvector $X_{n}^{\mathcal{L}}$ on
both sides of (92), we obtain

$\displaystyle\xi_{n}=-\mu_{n}\langle
X_{n}^{\mathcal{L}},\mathcal{L}_{2}\mathcal{L}_{3}X_{n}^{\mathcal{L}}\rangle.$
(93)

The brackets $\langle\cdot,\cdot\rangle$ stand for the scalar product in
$L^{2}([0,T];{\mathbb{R}}^{P})$ which is given by

$\displaystyle\langle\phi,\psi\rangle=\int_{0}^{T}\phi_{p}(t)\psi_{p}(t)dt.$
(94)

Using the fact that $\mathcal{L}_{2}^{*}=-\mathcal{L}_{3}$, we have

$\displaystyle\xi_{n}=\mu_{n}\|\mathcal{L}_{3}X_{n}^{\mathcal{L}}\|^{2},$ (95)

where $\|\cdot\|$ is the norm associated with $\langle\cdot,\cdot\rangle$.
Inserting this result as well as $\eta_{n}=0$ in (89), we deduce that

$\displaystyle\lambda^{(1)}_{n}=\mu_{n}(1+\alpha\|\mathcal{L}_{3}X_{n}^{\mathcal{L}}\|^{2}+O(\alpha^{2})).$
(96)

This shows that for $\alpha$ small enough, the sign of the real part of
$\lambda_{n}$ and $\mu_{n}$ are the same.

Now, in case (ii) where $\eta_{n}=-1$, at leading order 1 in $\alpha$ (88)
becomes

$\displaystyle-\nu_{n}X_{n}=\mathcal{L}_{2}\mathcal{L}_{3}X_{n}.$ (97)

Hence using again the fact that $\mathcal{L}_{2}^{*}=-\mathcal{L}_{3}$,
$\nu_{n}$ verifies

$\displaystyle\nu_{n}=\frac{\|\mathcal{L}_{3}X_{n}\|^{2}}{\|X_{n}\|^{2}},$
(98)

where $X_{n}$ is an eigenvector of $\mathcal{L}_{2}\mathcal{L}_{3}$. Therefore
we have another set of eigenvalues given by

$\displaystyle\lambda^{(2)}_{n}=-\alpha^{-1}+\|\mathcal{L}_{3}X_{n}^{\mathcal{L}_{2}\mathcal{L}_{3}}\|^{2}+O(\alpha),$
(99)

with $X_{n}^{\mathcal{L}_{2}\mathcal{L}_{3}}$ an eigenvector of
$\mathcal{L}_{2}\mathcal{L}_{3}$. Since $\lambda^{(2)}_{n}<0$ if $\alpha$ is
small enough, the associated eigenvectors are always stable. Therefore the
stability of the fixed points of the Hamiltonian system is determined by the
sign of $\lambda^{(1)}_{n}$, which is the same as the sign of the eigenvalues
$\mu_{n}$ of the Lagrangian system for $\alpha$ small enough.

##  Appendix B Convergence of the GDA in an analytically solvable case

To gain insight about the convergence of the GDA algorithm, we consider an
Ornstein-Uhlenbeck (OU) process in one dimension for which the evolution
equations of the GDA are amenable to analytic solution. The example is also
relevant since it is the dynamics verified by each Fourier mode of a freely
diffusive field in a one-dimensional periodic box. The stability of the
numerical scheme associated with this example will be analyzed in Appendix C.

For an OU process, the Hamiltonian is given by

$\displaystyle H(x,\theta)=-\zeta x\theta+\frac{1}{2}\theta^{2},$ (100)

where $\zeta>0$ is the stiffness of the confining harmonic potential. We look
for the instanton joining $x_{a}$ to $x_{b}\not=x_{a}$, with $x_{a}$ possibly
nonzero (i.e. not the stable point of the deterministic dynamics). The GDA
equations read

$\displaystyle\begin{cases}\partial_{\tau}x=\alpha(\partial_{t}\theta-\zeta\theta)\\\
\alpha\partial_{\tau}\theta=\partial_{t}x+\zeta x-\theta,\end{cases}$ (101)

subject to the boundary conditions

$\displaystyle x(\tau,0)=x_{a},\quad x(\tau,T)=x_{b},$ (102)

and to the initial conditions

$\displaystyle x(0,t)=x_{0}(t),\quad\theta(0,t)=\theta_{0}(t).$ (103)

Again, as in the body of text, we have introduced a scale $\alpha>0$ that
controls the relative evolution of $x$ and $\theta$ in time $\tau$. By taking
the derivative with respect to $\tau$ of the first equation in (101) and
applying the operator $\partial_{t}-\zeta$ to the second equation, a closed
PDE for the function $x(\tau,t)$ can be obtained:

$\displaystyle\partial_{\tau}^{2}x=-\frac{1}{\alpha}\partial_{\tau}x+\left(\partial_{t}^{2}x-\zeta^{2}x\right).$
(104)

For $\alpha>0$, we set $\gamma(t)=x_{a}(1-t/T)+x_{b}t/T$ and define

$\displaystyle\varphi(\tau,t)\equiv x(\tau,t)-\gamma(t)$ (105)

such that $\varphi(\tau,0)=\varphi(\tau,T)=0$. To proceed, let us use the
Fourier decomposition of $\varphi$:

$\displaystyle\varphi(\tau,t)=\sum_{q=1}^{\infty}\varphi_{q}(\tau)\sin(q\pi
t/T),$ (106)

where we have taken the convention

$\displaystyle\varphi_{q}(\tau)=\frac{2}{T}\int_{0}^{T}\varphi(\tau,t)\sin(q\pi
t/T)dt.$ (107)

Inserting (105) and (106) into (104) and projecting, we arrive at the equation
verified by the modes $\varphi_{q}$ ($q\in\mathbb{N}$):

$\displaystyle\alpha\partial_{\tau}^{2}\varphi_{q}+\partial_{\tau}\varphi_{q}+\alpha\omega_{q}^{2}\varphi_{q}=\alpha\zeta^{2}\frac{2(x_{a}-x_{b}(-1)^{q})}{\pi
q},$ (108)

with $\omega_{q}^{2}=\zeta^{2}+(\pi q/T)^{2}$. We recognize the equation for
the damped harmonic oscillator, implying that the modes will eventually decay
exponentially in $\tau$ to their stationary value. Setting

$\Delta_{q}\equiv 1-4\alpha^{2}\omega_{q}^{2}$ (109)

and defining

$\displaystyle\bar{\varphi}_{q}=\frac{\zeta^{2}(x_{a}-x_{b}(-1)^{q})}{\pi
q\omega_{q}^{2}}.$ (110)

we have: for $\Delta_{q}>0$

$\displaystyle\varphi_{q}(\tau)=\bar{\varphi}_{q}$
$\displaystyle+C_{1q}e^{-\tau/(2\alpha)}e^{\sqrt{\Delta_{q}}\tau/(2\alpha)}$
(111)
$\displaystyle+C_{2q}e^{-\tau/(2\alpha)}e^{-\sqrt{\Delta_{q}}\tau/(2\alpha)};$

for $\Delta_{q}<0$,

$\displaystyle\varphi_{q}(\tau)=\bar{\varphi}_{q}$
$\displaystyle+C_{1q}e^{-\tau/(2\alpha)}\cos\big{(}\sqrt{|\Delta_{q}|}\tau/(2\alpha)\big{)}$
(112)
$\displaystyle+C_{2q}e^{-\tau/(2\alpha)}\sin\big{(}\sqrt{|\Delta_{q}|}\tau/(2\alpha)\big{)};$

and for $\Delta_{q}=0$,

$\displaystyle\varphi_{q}(\tau)=\bar{\varphi}_{q}+(C_{1q}+C_{2q}\tau)e^{-\tau/(2\alpha)}.$
(113)

In each cases the constants $C_{1q}$ and $C_{2q}$ are then determined by the
initial conditions $x_{0}(t)$ and $\theta_{0}(t)$.

This computation addresses the question of the convergence rate, which is
different for each mode $q$. For a given $q\geq 1$, the best decay rate
$\lambda_{q}$ is obtained when choosing $\alpha=1/(2\omega_{q})$, which
corresponds to the damped critical regime. In practice, we need to use the
same $\alpha$ for all modes $q$. The convergence to the final path is then
determined by the smallest rate, which must be chosen as large as possible.
Since the decay rates $\lambda_{q}$ are ordered according to
$\lambda_{1}\leq\lambda_{2}\leq...\leq\lambda_{q}$, this prescribes the choice
of $\alpha=1/(2\omega_{1})$ that maximises the convergence rate of the mode
$q=1$, for which we have $\lambda_{1}=\omega_{1}$. This also yields identical
rate $\lambda_{q}=\omega_{1}$ for all modes $q$ and the modes $q>1$ will
display damped oscillations.

##  Appendix C Stability of the numerical scheme

In this section, we analyze the stability of Algorithm 1 for the Ornstein-
Uhlenbeck system introduced in Appendix B. Starting from (101), where $x$ is
subject to the boundary conditions (102), we set $u=x+\alpha\theta$ and
$v=x-\alpha\theta$. The GDA equations in these variables becomes

$\displaystyle\begin{cases}\partial_{\tau}u=\partial_{t}u+f(u,v)\\\
\partial_{\tau}v=-\partial_{t}v+g(u,v)\end{cases}$ (114)

with

$\displaystyle f(u,v)$
$\displaystyle=\left(\zeta+\frac{1}{\alpha}\right)v-\frac{1}{\alpha}u$ (115)
$\displaystyle g(u,v)$
$\displaystyle=\left(-\zeta+\frac{1}{\alpha}\right)u-\frac{1}{\alpha}v,$ (116)

and with the boundary conditions

$\displaystyle v(\tau,t=0)$ $\displaystyle=-u(\tau,t=0)+2x_{a}$ (117)
$\displaystyle u(\tau,t=0)$ $\displaystyle=-v(\tau,t=0)+2x_{b}.$ (118)

Since (114) are linear equations, Algorithm 1 can be expressed in terms of
matrix multiplications. To this end, let us define $r=\Delta\tau/\Delta t$,
the ratio between the algorithm evolution time step and the physical time
step, and set $U^{n}=(u_{0}^{n},\cdots,u_{M}^{n})^{T}$ and
$V^{n}=(v_{0}^{n},\cdots,v_{M}^{n})^{T}$, such that
$(U^{n},V^{n})^{T}=(u_{0}^{n},\cdots,u_{M}^{n},v_{0}^{n},\cdots,v_{M}^{n})^{T}$.
The scheme we prescribe in Algorithm 1 thus writes in a block-matrix form,
where each block $K$, $L$, $\\{G_{i}\\}_{1\leq i\leq 4}$ is a
$(M+1)\times(M+1)$ matrix, and where $\mathbf{1}$ is the identity matrix: we
solve $(U^{n+1},V^{n+1})^{T}$ such that

$\displaystyle\begin{pmatrix}K&0\\\
0&\mathbf{1}\end{pmatrix}\begin{pmatrix}U^{n+1}\\\ V^{n}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}G_{1}&G_{2}\\\
0&\mathbf{1}\end{pmatrix}\begin{pmatrix}U^{n}\\\ V^{n}\end{pmatrix}+R_{1},$
(119) $\displaystyle\begin{pmatrix}\mathbf{1}&0\\\
0&L\end{pmatrix}\begin{pmatrix}U^{n+1}\\\ V^{n+1}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\mathbf{1}&0\\\
G_{3}&G_{4}\end{pmatrix}\begin{pmatrix}U^{n+1}\\\ V^{n}\end{pmatrix}+R_{2},$
(120)

with

$\displaystyle K$ $\displaystyle=\begin{pmatrix}1+r&-r&0&\dots&0\\\
0&1+r&-r&\ddots&\vdots\\\ \vdots&\ddots&\ddots&\ddots&0\\\ \vdots&&0&1+r&-r\\\
0&&&0&1\end{pmatrix};$ (121) $\displaystyle L$
$\displaystyle=\begin{pmatrix}1&0&0&\dots&0\\\ -r&1+r&0&\ddots&0\\\
0&\ddots&\ddots&\ddots&\vdots\\\ \vdots&&-r&1+r&0\\\ 0&&0&-r&1+r\end{pmatrix}$
(122)

and, setting $z\equiv\Delta\tau(\zeta+\frac{1}{\alpha})$, and
$\bar{z}\equiv\Delta\tau(-\zeta+\frac{1}{\alpha})$,

$\displaystyle G_{1}$
$\displaystyle=\begin{pmatrix}1&-\frac{\Delta\tau}{\alpha}&0&\dots&0\\\
0&1&-\frac{\Delta\tau}{\alpha}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&0\\\ \vdots&&0&1&-\frac{\Delta\tau}{\alpha}\\\
0&&0&0&0\end{pmatrix};$ (123) $\displaystyle G_{2}$
$\displaystyle=\begin{pmatrix}0&z&0&\dots&0\\\ 0&0&z&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&0\\\ \vdots&&0&0&z\\\ 0&&0&0&-1\end{pmatrix};$
(124) $\displaystyle G_{3}$ $\displaystyle=\begin{pmatrix}-1&&&&\\\
\bar{z}&0&&\mathbf{0}&\\\ 0&\ddots&\ddots&&\\\ \vdots&&\bar{z}&0&\\\
0&&0&\bar{z}&0\end{pmatrix}$ (125) $\displaystyle G_{4}$
$\displaystyle=\begin{pmatrix}0&&&&\\\
-\frac{\Delta\tau}{\alpha}&1&&\mathbf{0}&\\\ 0&\ddots&\ddots&&\\\
\vdots&&-\frac{\Delta\tau}{\alpha}&1&\\\
0&&0&-\frac{\Delta\tau}{\alpha}&1\end{pmatrix},$ (126)

and with $R_{1}=(\mathbf{0}_{M},2x_{b},\mathbf{0}_{M+1})^{T}$ and
$R_{2}=(\mathbf{0}_{M+1},2x_{a},\mathbf{0}_{M})^{T}$, and where
$\mathbf{0}_{M}$ indicates a list with $M$ zeros.

The scheme is stable if the eigenvalues of the matrix $Q$

$\displaystyle Q\equiv\begin{pmatrix}\mathbf{1}&0\\\
0&L^{-1}\end{pmatrix}\begin{pmatrix}\mathbf{1}&0\\\
G_{3}&G_{4}\end{pmatrix}\begin{pmatrix}K^{-1}&0\\\
0&\mathbf{1}\end{pmatrix}\begin{pmatrix}G_{1}&G_{2}\\\
0&\mathbf{1}\end{pmatrix}$ (127)

are of module $\leq 1$. To check the stability for different values of
$\Delta\tau$ and $\Delta t$, we prescribe a value of $\zeta$, which defines a
typical relaxation time scale $t_{\zeta}=1/\zeta$ of the noiseless dynamics.
We then either prescribe the final time $T$ and vary the number of points $M$
along the path, or we fix the number of points and vary $\Delta t$, which
changes the final time $T$. In particular, it is important to check the
stability for a final time $T\gg t_{\zeta}$, such that a transition from
$x_{a}$ or $x_{b}$ to the critical point $x=0$ can follow the deterministic
flow on most of the path (such transitions only need noise close to the
critical point that must be reached exactly at time $t=T$). To correctly
resolve the instanton, we should then typically take $\Delta t<t_{\zeta}$.

In Fig. 11 we display the stability region of the scheme in space $(\Delta
t,\Delta\tau)$ for some values of the parameters. The result shows that the
stability region depends on $\zeta$ and $\alpha$. Since the fixed point (the
instanton) is independent of the dynamical parameters $\Delta\tau$ and
$\alpha$, we should take the best values that bring stability and convergence.
Interestingly sending $\alpha\to 0$ is not the best choice for stability, even
if it looks appealing at first sight since it corresponds to solving the
Legendre-Fenchel transform
$\operatorname*{argmax}_{\theta}(\langle\dot{x},\theta\rangle-H(x,\theta))$
for a given $x$. It is instead preferable to take $\alpha=O(1)$ to update $x$
and $\theta$ on similar time scales. It is also worth noticing that the use of
an implicit scheme for advection has relaxed the Courant-Friedrichs-Lewy (CFL)
stability condition, i.e. one can take $\Delta\tau>\Delta t$ and still have a
stable scheme. Unfortunately, it is clear from the graphs that the limit
$\Delta\tau\to\infty$ is not stable. This limit is however interesting since
it corresponds to an infinitely fast update of $u^{n+1}$ at $v^{n}$ fixed,
i.e. that $u^{n+1}$ solves in the continuous limit

$\displaystyle 0=\partial_{t}u^{n+1}+f(u^{n},v^{n})$ (128)

at $v^{n}$ fixed, and, following, $v^{n+1}$ solves

$\displaystyle 0=-\partial_{t}v^{n+1}+g(u^{n+1},v^{n+1})$ (129)

at $u^{n+1}$ fixed. This greedy procedure would bring the convergence to the
fixed point in a few steps, when it converged. Here instead, $\Delta\tau$ must
remain finite to guarantee convergence.

a)$\zeta=1$, $\alpha=1$

b)$\zeta=5$, $\alpha=1$

c)$\zeta=5$, $\alpha=0.05$

d)$\zeta=5$, $\alpha=20$

e)$\zeta=1$, $\alpha=1$

f)$\zeta=5$, $\alpha=1$

Figure 11: Stability of the numerical scheme represented as the shaded zone in
the space $(\Delta t,\Delta\tau)$. The dashed line materializes the usual CFL
conditions for explicit advection scheme, $\Delta\tau=\Delta t$, which is
relaxed here since the scheme is semi-implicit. For (a), (b), (c) and (d) we
fix the number of discretization points $M=100$ and $T$ varies according to
$T=M\Delta t$. For (e) and (f) we fix $T=10t_{k}$, and $M$ is then given by
$M=T/\Delta t$.

Finally, let us discuss the issue of preconditioning the scheme for large
values of $\zeta$, which is relevant e.g. if $\zeta$ refers to the different
relaxation rates of the different Fourier modes of a diffusive field (we
typically have $\zeta(k)=Dk^{2}$ with $D$ the diffusion coefficient of the
field). Specifically, we would like to choose time step $\Delta\tau$ for which
the numerical scheme remains stable when evolving separately each Fourier
component in a semi-spectral method for solving PDEs. Preconditioning the
evolution (101) aims at keeping the same $\Delta\tau$ for each mode, and
allowing higher modes (large $\zeta(k)$) to relax slower than lower modes
(small $\zeta(k)$). This procedure does not change the final path, which still
solves Hamilton’s equations. Such preconditioning for the GDA reads

$\displaystyle\begin{cases}\partial_{\tau}u=(1+\frac{1}{\alpha}+\zeta)^{-1}[\partial_{t}u+f(u,v)],\\\
\partial_{\tau}v=(1+\frac{1}{\alpha}+\zeta)^{-1}[-\partial_{t}v+g(u,v)].\end{cases}$
(130)

The shape of the matrices $K$, $L$ and $\\{G_{i}\\}_{1\leq i\leq 4}$ is
unchanged but $\Delta\tau$ should be modified following the substitution rule:

$\displaystyle\Delta\tau\to\Delta\tau(1+\frac{1}{\alpha}+\zeta)^{-1},$ (131)

which also modifies $r$, $z$, and $\bar{z}$. We find that this procedure
ensures the stability of the scheme up to $\Delta\tau=1$, for all $\zeta$,
assuming $\Delta t$ is small enough.

The semi-analytical proof presented here in the Ornstein-Uhlenbeck setup
provides insights on the behavior of the system, but the problems we are
usually interested in display strong non-linearities, and their stability
cannot be analyzed through the spectrum of a linear operator. We keep in mind
however that a few ingredients should be reused and that they stabilize the
code in general: the advection should be treated with an implicit upwind
scheme, the reaction terms can be evaluated upwind, and preconditioning the
dynamics allows us to take the same time step $\Delta\tau$ for each Fourier
mode in a semi-spectral scheme.

##  Appendix D Higher-order scheme

In Algorithms 1 and 2 that we have presented in the text, the derivative of
$u$ and $v$ with respect to physical time $t$ (or w.r.t parametrization $s$)
is approximated with a first-order finite difference upwind derivative. This
finite difference scheme is straightforward to implement but it is only first-
order accurate, i.e. the error with respect to the analytical solution
decreases as $O(M^{-1})$ when $M$ increases, $M$ being the number of points
used to discretize the interval $[0,T]$ (or interval $[0,1]$ if working with
reparametrized time $s$). Yet, it is possible to implement a higher-order
stable finite difference scheme for the advection, while keeping large steps
$\Delta\tau$. This can be done by keeping the implicit and upwind features of
the scheme, while approximating the derivatives with the second-order
difference stencil. The second-order stencil involves two upwind grid points,
which means that it can only be used for points $i=2,...,M$ in forward
advection for $v_{i}$, and for $i=M-2,...,0$ for backward advection for
$u_{i}$. The values $v_{1}$ and $u_{M-1}$ will still be computed with a first-
order finite-difference stencil. We have checked that the higher-order scheme
indeed significantly improves the accuracy.

The higher-order implementation of Algorithms 1 and 2 are detailed in
Algorithms 3 and 4, respectively. The algorithmic complexity of the higher-
order schemes is the same as the algorithmic complexity of lower-order ones.
Indeed, in both cases the implicit fields $u^{n+1}$ and $v^{n+1}$ simply solve
a triangular system.

Algorithm 3 : Action Minimization by Higher-Order Gradient Descent Ascent

1:Follow Steps 1 and 2 of Algorithm 1.

2:for $n\geq 0$ do

3: Update $u$ with an implicit upwind scheme, namely, solve
$\\{u^{n+1}_{i}\\}_{i\in I}$ sequentially from $i=M$ to $i=0$ using:
$\displaystyle\begin{cases}u^{n+1}_{M}=-v_{M}^{n}+2\phi_{b}\\\\[4.0pt]
\dfrac{u_{M-1}^{n+1}-u_{M-1}^{n}}{\Delta\tau}=\dfrac{u_{M}^{n+1}-u_{M-1}^{n+1}}{\Delta
t}+f(u^{n}_{M},v^{n}_{M}),\\\\[4.0pt]
\dfrac{u_{i}^{n+1}-u_{i}^{n}}{\Delta\tau}=\dfrac{-u_{i+2}^{n+1}+4u_{i+1}^{n+1}-3u_{i}^{n+1}}{2\Delta
t}+f(u^{n}_{i+1},v^{n}_{i+1}),\qquad i=M-2,\dots,0\end{cases}$

4: Update $v$ with an implicit upwind scheme, namely, solve
$\\{v_{i}^{n+1}\\}_{i\in I}$ sequentially from $i=0$ to $i=M$ using:
$\displaystyle\begin{cases}v^{n+1}_{0}=-u_{0}^{n+1}+2\phi_{a}\\\\[4.0pt]
\dfrac{v_{1}^{n+1}-v_{1}^{n}}{\Delta\tau}=-\dfrac{v_{1}^{n+1}-v_{0}^{n+1}}{\Delta
t}+g(u_{0}^{n+1},v_{0}^{n})\\\\[4.0pt]
\dfrac{v_{i}^{n+1}-v_{i}^{n}}{\Delta\tau}=-\dfrac{3v_{i}^{n+1}-4v_{i-1}^{n+1}+v_{i-2}^{n+1}}{2\Delta
t}+g(u_{i-1}^{n+1},v_{i-1}^{n}),\qquad i=2,\dots,M\end{cases}$

5: Compute $\\{\phi_{i}^{n+1}=\frac{1}{2}(u_{i}^{n+1}+v_{i}^{n+1})\\}_{i\in
I}$ and
$\\{\theta_{i}^{n+1}=\frac{1}{2}\alpha^{-1}(u_{i}^{n+1}-v_{i}^{n+1})\\}_{i\in
I}$ (if needed).

Algorithm 4 : Geometric Action Minimization by Higher-Order Gradient Descent
Ascent

1:Follow Steps 1 and 2 of Algorithm 2.

2:for $n\geq 0$ do

3: Update $u$ with an implicit upwind scheme, namely, solve
$\\{u^{n+1}_{i}\\}_{i\in I}$ sequentially from $i=M$ to $i=0$ using:
$\displaystyle\begin{cases}u^{n+1}_{M}=-v_{M}^{n}+2\phi_{b}\\\\[4.0pt]
\dfrac{u_{M-1}^{n+1}-u_{M-1}^{n}}{\Delta\tau}=\lambda_{M-1}(u^{n},v^{n})\dfrac{u_{M}^{n+1}-u_{M-1}^{n+1}}{\Delta
s}+f(u_{M-1}^{n},v_{M-1}^{n}),\\\\[4.0pt]
\dfrac{u_{i}^{n+1}-u_{i}^{n}}{\Delta\tau}=\lambda_{i}(u^{n},v^{n})\dfrac{-u_{i+2}^{n+1}+4u_{i+1}^{n+1}-3u_{i}^{n+1}}{2\Delta
s}+f(u_{i}^{n},v_{i}^{n}),\qquad i=M-2,\dots,0\end{cases}$

4: Update $v$ with an implicit upwind scheme, namely, solve
$\\{v^{n+1}_{i}\\}_{i\in I}$ sequentially from $i=0$ to $i=M$ using:
$\displaystyle\begin{cases}v^{n+1}_{0}=-u_{0}^{n+1}+2\phi_{a}\\\\[4.0pt]
\dfrac{v_{1}^{n+1}-v_{1}^{n}}{\Delta\tau}=-\lambda_{1}(u^{n+1},v^{n})\dfrac{v_{1}^{n+1}-v_{0}^{n+1}}{\Delta
s}+g(u_{1}^{n+1},v_{1}^{n}),\\\\[4.0pt]
\dfrac{v_{i}^{n+1}-v_{i}^{n}}{\Delta\tau}=-\lambda_{i}(u^{n+1},v^{n})\dfrac{3v_{i}^{n+1}-4v_{i-1}^{n+1}+v_{i-2}^{n+1}}{2\Delta
s}+g(u_{i}^{n+1},v_{i}^{n}),\qquad i=2,\ldots,M\\\ \end{cases}$

5: Follow Steps 5, 6 and 7 of Algorithm 2.

##  Appendix E Details on the calculation of $\lambda$ in the geometric
formulation

We detail here the steps that led us to consider the geometric formulation of
the gradient descent ascent (GDA) scheme with a specific treatment to
$\lambda$. We start from Eq. (50), and we notice that the minimization of
$\lambda$ can be performed before the minimization over ${\hat{\phi}}$. This
reads

$\min_{{\hat{\phi}}}\min_{\lambda\geq
0}\max_{{\hat{\theta}}}\int_{0}^{1}\left(\langle{\hat{\phi}}^{\prime},\hat{\theta}\rangle-\lambda^{-1}H({\hat{\phi}},\hat{\theta})\right)ds.$
(132)

The $\max$ on ${\hat{\theta}}$ can be performed pointwise in $s$ and this
gives the following relation that the minimum action path should satisfy:

$\displaystyle\lambda{\hat{\phi}}^{\prime}=\partial_{\hat{\theta}}H,$ (133)

which implicitly defines
${\hat{\theta}}(s)\equiv\vartheta({\hat{\phi}}(s),{\hat{\phi}}^{\prime}(s),\lambda(s))$.
Inserting this expression in (132), the minimization in $\lambda$ gives the
equation

$\displaystyle\langle{\hat{\phi}}^{\prime},\partial_{\lambda}\vartheta\rangle+\lambda^{-2}H-\lambda^{-1}\langle\partial_{\hat{\theta}}H,\partial_{\lambda}\vartheta\rangle=0,$
(134)

which can be simplified to obtain $\lambda$ explicitly. Indeed, by taking the
derivative of (133) with respect to $\lambda$, we get
${\hat{\phi}}^{\prime}=(\partial_{\hat{\theta}}^{2}H)\partial_{\lambda}\vartheta$,
or
$\partial_{\lambda}\vartheta=(\partial_{\hat{\theta}}^{2}H)^{-1}{\hat{\phi}}^{\prime}$,
since $C^{-1}\equiv(\partial_{\hat{\theta}}^{2}H)$ is invertible under
Assumption A3 that $H$ strictly convex in $\theta$. Inserting this expression
for $\partial_{\lambda}\vartheta$ in Eq. (134), and extracting the root, we
get

$\displaystyle\lambda=\frac{\langle\partial_{\hat{\theta}}H,C{\hat{\phi}}^{\prime}\rangle+\sqrt{\langle\partial_{\hat{\theta}}H,C{\hat{\phi}}^{\prime}\rangle^{2}-4H\langle{\hat{\phi}}^{\prime},C{\hat{\phi}}^{\prime}\rangle}}{2\langle{\hat{\phi}}^{\prime},C{\hat{\phi}}^{\prime}\rangle}.$
(135)

Eq. (52) follows from this equation if we modify a few terms to guarantee that
$\lambda\geq 0$ (which may not always be satisfied since (133) only holds at
convergence and not during the optimization). Note that equations (133) and
(135) impose $H=0$ along the trajectory, for any definite positive $C$. This
enjoins us to consider replacing $C$ by the identity in our numerical
algorithms in order to avoid computing $(\partial_{\hat{\theta}}^{2}H)$ and
its inverse. Note also that the value $H=0$ is the only possible value since
$H(0)=H(1)=0$ (endpoints are critical points) and that $H(s)=H(0)$ for every
$s\in[0,1]$ (Hamiltonian system). In this sense, the coefficient
$\lambda^{-1}$ can also be seen as a Lagrange multiplier enforcing $H=0$.
Finally, the $\min$ on ${\hat{\phi}}$ brings the second Hamilton equation

$\displaystyle\lambda{\hat{\theta}}^{\prime}=-\partial_{\hat{\phi}}H.$ (136)

##  Appendix F Modified Ginzburg-Landau dynamics discretized on two sites

a)

b)

c)

Figure 12: Comparison between paths in the low dimensional modified Ginzburg-
Landau system for different values of the final time $T=2,5,10$ and $20$.
Minimum action paths (left), Lagrangian (middle) and Hamiltonian (right) along
the paths. Orange squares (left): unstable fixed points. Orange disks: stable
fixed points. Background color: intensity of the force field, darker encodes
stronger force. As $T$ increases, the action decreases. For large $T$, the
Lagrangian plateaus in the vicinity of critical points. We notice that the
final Hamiltonian is constant along the trajectory. Parameters: $N=1200$,
$D=0.03$, $\kappa=0.26$, $\Delta\tau=0.002$.

To work in low dimensions, we discretize Eq. (64) on two sites with periodic
boundary conditions. Defining $x(t)$ and $y(t)$ the values of the field on
site 1 and site 2, respectively, the system is now equivalent to studying a
particle at position $(x,y)^{T}$ subjected to a nonequilibrium force. The
dynamics read

$\displaystyle\begin{cases}\dot{x}=b_{x}(x,y)+\sqrt{2\epsilon}\eta_{x}\\\
\dot{y}=b_{y}(x,y)+\sqrt{2\epsilon}\eta_{y},\end{cases}$ (137)

where, following Eq. (64), we set

$\displaystyle b_{x}(x,y)=8D(y-x)+x-x^{3}+\frac{\kappa}{2}(x^{2}+y^{2}),$
(138) $\displaystyle
b_{y}(x,y)=8D(x-y)+y-y^{3}+\frac{\kappa}{2}(x^{2}+y^{2}),$ (139)

and where $\eta_{x}$ and $\eta_{y}$ are Gaussian white noise of variance
unity. The term $\kappa(x^{2}+y^{2})/2$ is the two-site analogue of
$\kappa\int\rho^{2}$. For simplicity, we have taken $h=0$, such that the
noiseless system has always two stable fixed points in
$(x_{1},y_{1})=(q_{1},q_{1})$ and $(x_{2},y_{2})=(q_{2},q_{2})$, with
$q_{1}=(-\kappa-\sqrt{\kappa^{2}+4})/2$ and
$q_{2}=(-\kappa+\sqrt{\kappa^{2}+4})/2$, and one unstable fixed point at
$(x_{u},y_{u})=(0,0)$. The system displays time-reversal symmetry for
$\kappa=0$, and we therefore know that the dynamics follows a stochastic
gradient descent on the landscape
$V(x,y)=-x^{2}/2+x^{4}/4-y^{2}/2+y^{4}/4+4D(x-y)^{2}$, and that the minimum
action path in this case corresponds to the minimum energy path, which simply
follows the gradient ascent to escape the basin of attraction of one or the
other stable fixed point. For $\kappa\neq 0$ however, we need to resort to
numerical algorithms to find the new minimum action paths between the stable
fixed points. In this diffusive system, as presented in Sec. II.1, the
Hamiltonian is explicitly given by

$\displaystyle
H=b_{x}(x,y)\theta_{x}+\theta_{x}^{2}+b_{y}(x,y)\theta_{y}+\theta_{y}^{2},$
(140)

where $\theta_{x}$ and $\theta_{y}$ are the conjugated fields of $x$ and $y$,
respectively. The solution to Hamilton’s equation of motion is found via
iteration of the gradient ascent-descent procedure introduced above.

a)

b)

Figure 13: Comparison between paths in the low dimensional modified Ginzburg-
Landau system. Left panel: flow lines and paths. Right panel: action increment
along the paths parametrized by normalized arclength. (a) Dashed line:
heteroclinic orbit obtained with the string method. Purple line: forward
minimizer obtained with the geometric algorithm. Orange line: backward
minimizer obtained with the geometric algorithm. Squares: unstable fixed
points. Disks: stable fixed points. Background color: intensity of the force
field, darker encodes stronger force. Parameters: $N=600$, $D=0.03$,
$\kappa=0.26$, $\Delta\tau=0.01$.

The results can be found in Fig. 12. We have checked that the scheme perfectly
recovers the prediction of previous methods [35, 36, 38]. As expected, the
action decreases when the final time $T$ increases, see Fig. 12b). Also, the
Hamiltonian $H$ is conserved along the path since the final path solves
Hamilton’s equation of motion. In addition we notice that $H\to 0$ when $T$
increases.

Now, we know that the minimal action is reached when we allow $T\to\infty$,
see Sec. III.2. However, since the number of points along the trajectory must
remain finite, increasing $T$ translates into an increase of the physical time
step $\Delta t$ and decreases the resolution of the path. To overcome this
issue, we use the geometric parametrization of the path presented in Sec.
III.2. In Fig. 13, we see that the geometric parametrization allows us to
reach the minimum of the action. We have also checked that $H=0$ along the
trajectory.

This example has illustrated the approach on low-dimensional systems subjected
to an additive Gaussian noise. The method yields the optimal paths of finite
duration but also the paths of infinite time length if the geometric
formulation is used. Note that the numerical scheme can be adapted to
determine transition paths for spatially extended fields: this is what we do
in Sections IV and V.

##  Appendix G Large $\Omega$ expansion and probability density

Figure 14: Rescaled and shifted logarithms of the probability density on a
unique site in the Schlögl model when the local number of particle $\Omega$
increases. The probability is explicitly known for this equilibrium case, and
we show the convergence to the function $V_{\mathrm{eq}}(\rho)$ as the number
of particles in the box increases. For the parameters chosen here, $\rho_{-}$
is the stable phase, while looking at the naive Ginzburg-Landau potential
$U(\rho)=-\lambda_{0}\rho-\lambda_{2}\rho^{3}/3+\lambda_{1}\rho^{2}/2+\lambda_{3}\rho^{4}/4$
would wrongly predict that $\rho_{+}$ is the stable phase. Parameters:
$\lambda_{0}=0.8$, $\lambda_{1}=2.9$, $\lambda_{2}=3.1$, $\lambda_{3}=1$.

In this section, we focus on the case of a unique well-stirred compartment of
fixed volume. We work in the limit of a large number of particles, and we
choose the microscopic rates $k_{i}$ such that the mean-field equation always
displays two stable fixed points $n_{-}$ and $n_{+}$ that solve
$0=k_{0}-k_{1}n+k_{2}n(n-1)-k_{3}n(n-1)(n-2)$, and where $n$ is the number of
particles in the compartment. When all the particles coexist in the same
compartment, the evolution of the probability $P(n)$ to find a number $n$ of
particles $X$ at time $t$ is given by

$\displaystyle\begin{aligned}
\partial_{t}P(n)=&W_{+}(n-1)P(n-1)+W_{-}(n+1)P(n+1)\\\
&-(W_{+}(n)+W_{-}(n))P(n),\end{aligned}$ (141)

where the rates are

$\displaystyle W_{+}(n)$ $\displaystyle=k_{0}+k_{2}n(n-1)$ (142)
$\displaystyle W_{-}(n)$ $\displaystyle=k_{1}n+k_{3}n(n-1)(n-2).$ (143)

The stationary probability $P_{\mathrm{eq}}$ can be obtained explicitly

$\displaystyle
P_{\mathrm{eq}}(n)=K\prod_{i=1}^{n}\frac{W_{+}(i-1)}{W_{-}(i)}=K\exp\left(\sum_{i=1}^{n}\ln\frac{W_{+}(i-1)}{W_{-}(i)}\right)$
(144)

where $K$ is a normalization constant. We are interested in the case of a
large number of particles per compartment, such that we can extract a large
deviation principle. We denote $\Omega\gg 1$ the typical number of particles
in the compartment, we define the rescaled number of particles
$\rho=n/\Omega$, and we write the rates as

$\displaystyle W_{+}(n)$ $\displaystyle=\Omega[w_{+}(\rho)+O(\Omega^{-1})]$
(145) $\displaystyle W_{-}(n)$
$\displaystyle=\Omega[w_{-}(\rho)+O(\Omega^{-1})],$ (146)

with $w^{+}(\rho)=\lambda_{0}+\lambda_{2}\rho^{2}$ and
$w^{-}(\rho)=\lambda_{1}\rho+\lambda_{3}\rho^{3}$, and where the $\lambda_{i}$
are now rescaled reaction rates verifying $\lambda_{i}=k_{i}\Omega^{i-1}$.
This rescaling ensures that the deterministic mean-field dynamics

$\displaystyle\partial_{t}\rho=\lambda_{0}+\lambda_{2}\rho^{2}-\lambda_{1}\rho-\lambda_{3}\rho^{3}$
(147)

keeps the same fixed points $\rho_{-}$, $\rho_{s}$, $\rho_{+}$ when
$\Omega\to\infty$, with $\rho_{-}$ and $\rho_{+}$ stable, while $\rho_{s}$ is
unstable. Following [83, 47, 48, 84], using the WKB (or eikonal) approximation
and the continuum limit, the probability now becomes a probability density and
can be cast into the following form

$\displaystyle P_{\mathrm{eq}}(\rho)=K(\rho,\Omega)e^{-\Omega
V_{\mathrm{eq}}(\rho)},$ (148)

where the equilibrium potential is given by

$\displaystyle
V_{\mathrm{eq}}(\rho)=\int^{\rho}dy\ln\frac{w_{+}(y)}{w_{-}(y)}$ (149)

and where $K(\rho,\Omega)$ is a function with the property that for any