Most probable transition paths in piecewise-smooth stochastic differential equations

During rare event transitions from one metastable state to another, a stochastic system will follow paths according to some distribution, which in the limit of vanishing noise strength is generally singly-peaked along a most probable transition path. Following [23], for a smooth vector field $\mathbfcal{F}$, the most probable transition path between $\mathbf{x}_{0}$ and $\mathbf{x}_{f}$ can be defined as follows. We first define an admissible set of transition paths ${\mathcal{A}^{(t_{0},t_{f})}}$ by

$${\mathcal{A}^{(t_{0},t_{f})}}=\{\bm{\alpha}\in H^{1}([t_{0},t_{f}];{\mathbb{R}}^{n}):\bm{\alpha}(t_{0})=\mathbf{x}_{0}\text{ and }\bm{\alpha}(t_{f})=\mathbf{x}_{f}\},$$

where $H^{1}([t_{0},t_{f}];{\mathbb{R}}^{n})$ is the Hilbert space of weakly differentiable curves $\bm{\alpha}$ satisfying $\bm{\alpha}\in H^{1}$ if and only if $\int_{t_{0}}^{t_{f}}|\bm{\alpha}|^{2}+|\dot{\bm{\alpha}}|^{2}\,dt<\infty$. The most probable transition path $\bm{\alpha}^{*}\in{\mathcal{A}^{(t_{0},t_{f})}}$ is then defined to be the global minimizer of the Freidlin-Wentzell rate functional ${I^{(t_{0},t_{f})}}:{\mathcal{A}^{(t_{0},t_{f})}}\mapsto\overline{{\mathbb{R}}},$ given by

$${I^{(t_{0},t_{f})}}[\bm{\alpha}]=\int_{t_{0}}^{t_{f}}\left\|\dot{\bm{\alpha}}(t)-\mathbfcal{F}(\bm{\alpha}(t),t)\right\|^{2}\,dt,$$

(1)

so that $\bm{\alpha}^{*}=\operatorname*{arg\,min}_{\bm{\alpha}\in{\mathcal{A}^{(t_{0},t_{f})}}}{I^{(t_{0},t_{f})}}[\bm{\alpha}]$. For this functional, the necessary condition satisfied by second differentiable minimizers is the following Euler-Lagrange equations [25]:

$\displaystyle\ddot{\bm{\alpha}}=\mathbfcal{F}_{t}+\dot{\bm{\alpha}}\left(\nabla\mathbfcal{F}-\nabla\mathbfcal{F}^{\intercal}\right)+\nabla\mathbfcal{F}^{\intercal}\mathbfcal{F}.$

Here, the subscript $t$ refers to the partial derivative with respect to time.

We determine the most probable transition path of a noise-induced tipping event when the drift is piecewise-smooth. In general, for systems with a piecewise-smooth drift, the appropriate rate functional to minimize is not known. Minimizers of the Freidlin-Wentzell rate functional may not be well-defined when a region of the vector field is not continuous, or even Lipschitz continuous. It is not clear how the most probable path may traverse across a switching manifold, which may have an attracting or repelling sliding region that introduces more complex dynamics. Specifically, letting $\mathbf{F}^{\pm}:{\mathbb{R}}^{n}\mapsto{\mathbb{R}}^{n}$ be smooth vector fields, we consider a system of stochastic differential equations of the form

$$d\mathbf{x}_{t}=\mathbf{F}(\mathbf{x}_{t})dt+\sigma d\mathbf{W}_{t},$$

(2)

where $\mathbf{x}=(x,\mathbf{y})\in{\mathbb{R}}^{n}$ such that $x\in{\mathbb{R}}$ and $\mathbf{y}\in{\mathbb{R}}^{n-1}$, $\sigma\in{\mathbb{R}}$, $\mathbf{W}=(W_{1},\ldots W_{n})$ is an $n$-dimensional Wiener process, and $\mathbf{F}:{\mathbb{R}}^{n}\setminus\{x=0\}\mapsto{\mathbb{R}}^{n}$ is defined by

$$\mathbf{F}(\mathbf{x})=\begin{cases}\mathbf{F}^{+}(\mathbf{x}),&x>0,\\ \mathbf{F}^{-}(\mathbf{x}),&x<0.\end{cases}$$

(3)

We let $S_{+}=\{\mathbf{x}\in{\mathbb{R}}^{n}:x>0\}$, $S_{-}=\{\mathbf{x}\in{\mathbb{R}}^{n}:x<0\}$, and $\Sigma=\{\mathbf{x}\in{\mathbb{R}}^{n}:x=0\}$. The set $\Sigma$ is often called the switching manifold or discontinuity boundary [51]. In general, $\Sigma$ may be defined as the zero level set of a smooth function $H:{\mathbb{R}}^{n}\mapsto{\mathbb{R}},$ but we assume for simplicity that $\Sigma$ is a hyperplane in ${\mathbb{R}}^{n-1}$.

We define the deterministic skeleton of Equation (2) as the dynamical system

$$\dot{\mathbf{x}}=\mathbf{F}(\mathbf{x}).$$

(4)

We assume that there exist asymptotically stable fixed points $\mathbf{x}_{0}\in S_{+}$ and $\mathbf{x}_{f}\in S_{-}$ of Equation (4), such that $\mathbf{F}^{+}(\mathbf{x}_{0})=0$ and $\mathbf{F}^{-}(\mathbf{x}_{f})=0$. By noise-induced transitions, we mean realizations of Equation (2) that transition from the basin of attraction of $\mathbf{x}_{0}$ to that of $\mathbf{x}_{f}$. More precisely, we define a noise-induced transition on the interval $[t_{0},t_{f}]$ as a solution to the stochastic boundary value problem given by realizations of Equation (2) that satisfy $\mathbf{x}(t_{0})=\mathbf{x}_{0}$ and $\mathbf{x}(t_{f})=\mathbf{x}_{f}$. Here the boundary conditions mean that the process is conditioned to transition from $\mathbf{x}_{0}$ to $\mathbf{x}_{f}$.

Dynamics that occur entirely within the smooth regions $S_{\pm}$ can be fully described by regular (smooth) dynamical systems theory. On the other hand, dynamics on a switching manifold $\Sigma$ may not be defined and are typically imposed a posteriori, e.g. using Filippov’s convex combination [51, 52, 53]. In general, the lack of smoothness requirements across the switching manifold allows for more diverse phenomena than in smooth systems of similar dimension, since the vector field across $\Sigma$ may be discontinuous or even point in opposite directions. Regions of $\Sigma$ where either occurs are called sliding regions. A solution that reaches $\Sigma$ at a sliding region may “slide” along it. On the other hand, crossing regions occur on $\Sigma$ where sliding is not possible. We differentiate these regions using the following notation:

1. 1.

   Attracting sliding regions: $\Sigma_{A}=\{\mathbf{x}\in\Sigma:F_{1}^{+}(0,\mathbf{y})\leq 0,\,F_{1}^{-}(0,\mathbf{y})\geq 0\}$,

2. 2.

   Repelling sliding regions: $\Sigma_{R}=\{\mathbf{x}\in\Sigma:F_{1}^{+}(0,\mathbf{y})\geq 0,\,F_{1}^{-}(0,\mathbf{y})\leq 0\},$

3. 3.

   Positive crossing regions: $\Sigma_{+}=\{\mathbf{x}\in\Sigma:F^{\pm}_{1}(0,\mathbf{y})>0\}$,

4. 4.

   Negative crossing regions: $\Sigma_{-}=\{\mathbf{x}\in\Sigma:F^{\pm}_{1}(0,\mathbf{y})<0\}$,

where ${F}^{\pm}_{1}$ is the first component of $\mathbf{F}^{\pm}$.

When solutions may slide along the switching manifold $\Sigma,$ we impose a flow using the Filippov convex method [53, 51]. That is, we define the sliding flow as a convex combination of $\mathbf{F}^{\pm},$

$$\mathbf{F}^{s}(0,\mathbf{y})=\lambda\mathbf{F}^{+}(0,\mathbf{y})+(1-\lambda)\mathbf{F}^{-}(0,\mathbf{y})$$

(5)

with

$$\lambda(\mathbf{y})=\frac{F_{1}^{-}(0,\mathbf{y})}{F_{1}^{-}(0,\mathbf{y})-F_{1}^{+}(0,\mathbf{y})}\in[0,1].$$

(6)

Notice that Equation (6) fixes $\lambda(\mathbf{y})$. This flow naturally arises in the context of our main result without a priori imposing it; see Theorem 3.8. Also, note that $\mathbf{F}^{s}(0,\mathbf{y})$ may vanish for some point $(0,\mathbf{y})\in\Sigma$. Since it is not an equilibrium of the smooth flow, we say that such a point $\mathbf{x}_{s}=(0,\mathbf{y}_{s})\in\Sigma$ is a pseudoequilibrium if it is an equilibrium of the sliding flow; i.e., for some $\lambda\in(0,1),$ $\mathbf{F}^{s}(0,\mathbf{y}_{s})=\mathbf{0}$.

Filippov’s convex combination interprets solutions of System (3),(4) as a continuous curve $\mathbf{x}(t)$ satisfying the following conditions:

$$\begin{cases}\dot{\mathbf{x}}(t)=\mathbf{F}(\mathbf{x}),\\ \mathbf{x}(0)=\mathbf{x}_{0},\end{cases}$$

(7)

where $\mathbf{F}:{\mathbb{R}}^{n}\mapsto{\mathbb{R}}^{n}$ is defined as

$$\mathbf{F}(\mathbf{x})=\begin{cases}\mathbf{F}^{+}(\mathbf{x}),&\mathbf{x}\in S_{+}\cup\Sigma_{+},\\ \mathbf{F}^{-}(\mathbf{x}),&\mathbf{x}\in S_{-}\cup\Sigma_{-},\\ \mathbf{F}^{s}(\mathbf{x}),&\mathbf{x}\in\Sigma_{A}\cup\Sigma_{R}.\end{cases}$$

(8)

That is, for points in time in which a solution curve $\mathbf{x}(t)$ intersects $\Sigma$ it will either cross $\Sigma,$ tracking the flow of either $\mathbf{F}^{+}$ or $\mathbf{F}^{-}$ depending on whether the curve is crossing into $S_{+}$ or $S_{-}$, or it will track $\Sigma$ in the sliding regions until entering a crossing region. Note, implicit in this definition is that solution curves are differentiable everywhere except on $\Sigma$.

To provide a brief illustration, Figure 1(a)-(c) shows some possible dynamics across $\Sigma$ in a two-dimensional system, with (a) $\Sigma_{+}$, (b) $\Sigma_{\pm}$ and $\Sigma_{A}$, and (c) $\Sigma_{\pm}$ and $\Sigma_{R}$. In (a) the flow traverses across $\Sigma$ continuously but not smoothly; in (b) and (c) the flow may also remain on $\Sigma,$ in intervals indicated by the blue line (for $\Sigma_{A}$) and the green line (for $\Sigma_{R}$). Note that the flow leaving a repelling sliding region is non-unique in forward time, while the flow leaving an attracting sliding region is non-unique in backward time. In addition, complex dynamics may appear on $\Sigma$ itself, depending on the manner in which the imposed vector field on $\Sigma$ transitions between $S_{+}$ and $S_{-}$. Figure 1(d) shows one example of a path traversing from a stable equilibrium in $S_{+}$, through the switching manifold in an attracting sliding region using dynamics imposed using Filippov’s convex combination [53], then traversing to a stable equilibrium in $S_{-}$. The system in this illustration is a piecewise-smooth version of the Lorenz 63 model,

$$\begin{pmatrix}\dot{x}\\ \dot{y}\\ \dot{z}\end{pmatrix}=\begin{pmatrix}\sigma_{-}\left(y-y_{0}-x+x_{0}\right)\\ (x-x_{0})(\rho_{-}-z+z_{0})-y+y_{0}\\ (x-x_{0})(y-y_{0})-\beta_{-}(z-z_{0})\end{pmatrix},\;x<0,\mbox{ and}\;\begin{pmatrix}\dot{x}\\ \dot{y}\\ \dot{z}\end{pmatrix}=\begin{pmatrix}\sigma_{+}\left(y-x\right)\\ x(\rho_{+}-z)-y\\ xy-\beta_{+}z\end{pmatrix},\;x>0,$$

(9)

originally studied in [56]. Here, $\sigma_{\pm},\rho_{\pm},\beta_{\pm}>0$ are parameters and $x_{0},y_{0},z_{0}\in{\mathbb{R}}$. For this example, $S_{\pm}$ and $\Sigma$ are defined as previously, with $n=3$.

For our problem, noise-induced transitions must cross $\Sigma$ and the Freidlin-Wentzell rate functional, Equation (1), must be modified to account for the discontinuity in $\mathbf{F}$ and possible discontinuity-induced dynamics. For example, a most probable path that reaches a repelling sliding region $\Sigma_{R}$ leaves it non-uniquely due to the non-uniqueness of the flow, as illustrated in the case studies in Sections 4 and 5. Additionally, if one were to naively piece together the most probable paths of the smooth regions and neglect possible contribution from the manner in which the path crosses the switching manifold, this may lead to paths that do not reflect a global minimum of the rate functional.

We resolve these issues by smoothing out $\mathbf{F}$ via mollification with a compactly supported, smooth, radially symmetric kernel of characteristic width $\varepsilon$, and consider the sequence of minimizers to the rate functional Equation (1) as $\varepsilon\rightarrow 0$. We show for any piecewise-smooth system of the form given by System (2),(3) that the most probable path minimizes Equation (1) in the smooth regions $S_{\pm}$, plus an additional functional whose contribution represents time spent sliding in a crossing region of the switching manifold. For $\mathbf{F}=(F_{1},\mathbf{G})$ such that $F_{1}:{\mathbb{R}}^{n}\setminus\Sigma\mapsto{\mathbb{R}}$ and $\mathbf{G}:{\mathbb{R}}^{n}\setminus\Sigma\mapsto{\mathbb{R}}^{n-1},$ and $\bm{\alpha}=(\alpha,\bm{\beta})\in{\mathcal{A}^{(t_{0},t_{f})}}$ and $0<\varepsilon\ll 1$, the appropriate rate functional with both contributions is

	$\displaystyle{\overline{I}^{(t_{0},t_{f})}}[\bm{\alpha}]=$	$\displaystyle\int_{\{t\in S_{+}\cup S_{-}\}}\|\dot{\bm{\alpha}}(t)-\mathbf{F}(\bm{\alpha}(t))\|^{2}\,dt$
		$\displaystyle+\int_{\{t\in\Sigma_{\pm}\}}\min_{\lambda\in[0,1]}\left\{\left[\lambda F_{1}^{+}(0,\bm{\beta})+(1-\lambda)F_{1}^{-}(0,\bm{\beta})\right]^{2}\right.$
		$\displaystyle\left.\mathmakebox[\widthof{$ +\int_{\mathcal{I}_{\Sigma}[\bm{\alpha}]}\min\qquad$}][l]{}+\left|\dot{\bm{\beta}}-\lambda\mathbf{G}^{+}(0,\bm{\beta})-(1-\lambda)\mathbf{G}^{-}(0,\bm{\beta})\right|^{2}\right\}\,dt.$

We make the following observations about this result:

1. 1.

   The form of the above rate functional is independent of the chosen mollifier.

2. 2.

   When there is a sliding region, the most probable path may track the Filippov dynamics; that is, no additional contribution is made during times where the most probable path slides via the imposed flow on the switching manifold.

We rigorously prove the above result using a technique from the calculus of variations called $\Gamma$-convergence to show that ${\overline{I}^{(t_{0},t_{f})}}$ is the $\Gamma$-limit of the sequence of rate functionals for the mollified system [57, 58]. The $\Gamma$-limit is the natural notion of a limiting functional in this context in the sense that minimizers of the rate functional with the mollified drift converge (weakly) to a minimizer of ${\overline{I}^{(t_{0},t_{f})}}$ in the limit as $\varepsilon\rightarrow 0$. While $\Gamma$-convergence has been used to study the convergence of minimizers of the Onsager-Machlup functional in the small noise limit [59, 60, 61, 62, 63], to our knowledge we are the first to use $\Gamma$-convergence to compute an appropriate extension of the Freidlin-Wentzell rate functional for piecewise smooth systems.

The utility of the derived rate functional is demonstrated through two case studies. Each case study also presents distinct phenomena in their most probable transition paths that are not possible in systems with smooth drift. The first case study analyzes the most probable transition path in a two-dimensional piecewise-linear system, a simple demonstrative case in which the typical rate functional formulation of the most probable path breaks down. Notably, in this case the most probable path follows the switching manifold along an attracting sliding region, but it does not follow the switching manifold along a repelling sliding region.

The second case study analyzes the most probable path in a simple one-dimensional periodically-forced piecewise-smooth system, constructed similarly to some conceptual models for Arctic energy balance [5, 64]. A significant portion of the analysis of the deterministic dynamics and the Monte Carlo simulations were first derived in Zanetell’s Master’s thesis [65]. Here we again observe the breakdown of the typical formulation of the most probable path. Additionally, we observe an emergent dynamical phenomenon, which we describe as noise-induced sliding, in which the most probable path follows the switching manifold in a crossing region, despite the energy cost of the associated rate functional.

In both case studies, we note in some regimes the most probable path(s) predicted by our result do not match observed Monte Carlo simulations. Several factors may contribute to this discrepancy, including needing a smaller noise strength $\sigma,$ a higher-order functional in the smoothing parameter $\varepsilon,$ or including the correction term in the Onsager-Machlup functional. Such possibilities will be explored in future work.

The paper is organized as follows: in Section 2 we briefly highlight notation we use and some technical assumptions on $\mathbf{F}$. Then in Section 3 we derive the mollified vector field and use $\Gamma-$convergence to determine the appropriate limiting functional for calculating the most probable path for stochastic differential equations of the form given by System (2),(3). Next, in Section 4 we provide a simple case study of a planar linear piecewise-smooth system as an illustration of the most probable path calculated in Section 3. Finally, in Section 5 we provide a case study that extends the most probable path derivation to the case of a one-dimensional non-autonomous system.

We begin this section by recalling the definition of a (Friedrichs) mollifier, following the presentation in [66]. Let $\zeta:{\mathbb{R}}\mapsto{\mathbb{R}}$ be a smooth even bump function, defined by

$$\zeta(x)=\begin{cases}Ch(x),&|x|<1,\\ 0,&|x|\geq 1,\end{cases}$$

(14)

where $h:\mathbb{R}\rightarrow\mathbb{R}$ is a smooth function satisfying $h(x)\geq 0$, $h(-1)=h(1)=0$ and all of its derivatives vanish at $x=-1$ and $x=1$, e.g., $h(x)=\exp\left(1/\left(|x|^{2}-1\right)\right)$. Additionally, $C>0$ is determined so that $\int_{-\infty}^{\infty}\zeta(x)\,dx=1$; i.e., $C=1/\int_{-\infty}^{\infty}\,h(x)dx$. For each $\varepsilon>0,$ we define a sequence of functions $\zeta_{\varepsilon}(x)=\zeta\left(x/\varepsilon\right)/\varepsilon$ which satisfy $\int_{-\varepsilon}^{\varepsilon}\zeta_{\varepsilon}(x)\,dx=1$ and are compactly supported on $[-\varepsilon,\varepsilon]$. The mollification $\textbf{F}^{\varepsilon}$ of F in $x$ is then defined by convolution with $\zeta_{\varepsilon}$ in the $x-$direction:

	$\displaystyle\textbf{F}^{\varepsilon}(x,\mathbf{y})=\zeta_{\varepsilon}(x)*\textbf{F}(x,\mathbf{y})=$	$\displaystyle\int_{-\infty}^{0}\zeta_{\varepsilon}(x-s)\mathbf{F}^{-}(s,\mathbf{y})\,ds+\int_{0}^{\infty}\zeta_{\varepsilon}(x-s)\mathbf{F}^{+}(s,\mathbf{y})\,ds,$		(15)
	$\displaystyle=$	$\displaystyle\begin{cases}\displaystyle\int_{-\varepsilon}^{\varepsilon}\zeta_{\varepsilon}(u)\mathbf{F}^{-}(x-u,\mathbf{y})\,du,&x\leq-\varepsilon,\\ \displaystyle\int_{-\varepsilon}^{x}\zeta_{\varepsilon}(u)\mathbf{F}^{+}(x-u,\mathbf{y})\,du+\int_{x}^{\varepsilon}\zeta_{\varepsilon}(u)\mathbf{F}^{-}(x-u,\mathbf{y})\,du,&|x|<\varepsilon,\\ \displaystyle\int_{-\varepsilon}^{\varepsilon}\zeta_{\varepsilon}(u)\mathbf{F}^{+}(x-u,\mathbf{y})\,du,&x\geq\varepsilon,\\ \end{cases}$
	$\displaystyle=$	$\displaystyle\int_{-\varepsilon}^{\varepsilon}\zeta_{\varepsilon}(u)\left(\mathbf{F}^{-}(x-u,\mathbf{y})\mathbbm{1}_{\{u\geq x\}}(u)+\mathbf{F}^{+}(x-u,\mathbf{y})\mathbbm{1}_{\{u\leq x\}}(u)\right)du,$

where $\mathbbm{1}_{A}$ is the indicator function on the set $A$. The mollified Freidlin-Wentzell rate functional ${I_{\varepsilon}^{(t_{0},t_{f})}}:{\mathcal{A}^{(t_{0},t_{f})}}\mapsto\mathbb{R}$ is then defined by:

$${I_{\varepsilon}^{(t_{0},t_{f})}}[\bm{\alpha}]=\int_{t_{0}}^{t_{f}}\left|\dot{\bm{\alpha}}-{\mathbf{F}^{\varepsilon}}(\bm{\alpha})\right|^{2}\,dt.$$

(16)

Note, by construction, ${\mathbf{F}^{\varepsilon}}(\mathbf{x})$ is a smooth function that converges pointwise to $\mathbf{F}(\mathbf{x})$ as $\varepsilon\rightarrow 0$ except at $x=0$ [66]. However, since $\mathbf{F}(\mathbf{x})$ is not continuous this convergence is necessarily not uniform and thus it is not a priori clear which properties of $\mathbf{F}$ are preserved under mollification. The following lemmas ensure that ${\mathbf{F}^{\varepsilon}}$ also has $p$-growth and is asymptotically inward-flowing. These properties are crucial to establishing the existence of a minimum for ${I_{\varepsilon}^{(t_{0},t_{f})}},$ but to streamline our main results, the proofs of these lemmas are in Appendix A.

We now prove the existence of a minimizer of ${I_{\varepsilon}^{(t_{0},t_{f})}}$ using Tonelli’s direct method of the calculus of variations; see, e.g., [67, 58]. This procedure consists of first proving that a minimizing sequence of ${I_{\varepsilon}^{(t_{0},t_{f})}}$ is bounded with respect to the $H^{1}$ norm and thus has a weakly convergent subsequence in the $H^{1}$ topology. Note, it is necessary to work in the $H^{1}$ topology since closed and bounded sets are not necessarily compact in the strong topology. Nevertheless, upon passing to a subsequence, we can establish a limit point for the minimizing sequence. Finally, by proving that ${I_{\varepsilon}^{(t_{0},t_{f})}}$ is lower semi-continuous with respect to weak convergence in the $H^{1}$ topology, it follows that limit points of the minimizing sequence are indeed minimizers.

To begin the procedure for the direct method of the calculus of variations outlined in the previous paragraph, we first recall the definition of weak convergence in the $H^{1}$ topology in the context of our problem.

The following lemma establishes that boundedness of ${I_{\varepsilon}^{(t_{0},t_{f})}}[\bm{\alpha}_{m}]$ implies boundedness of $\bm{\alpha}_{m}$ with respect to the $H^{1}$ norm.

Finally, we use the direct method of the calculus of variations to prove the existence of a minimum.

Now that the existence of minimizers is established for the Freidlin-Wentzell rate ${I_{\varepsilon}^{(t_{0},t_{f})}}[\bm{\alpha}_{\varepsilon}]$ corresponding to the mollified system $\dot{\mathbf{x}}={\mathbf{F}^{\varepsilon}}$, we can establish a necessary condition for minimizers, which is that $\bm{\alpha}\in{\mathcal{A}_{\Sigma}^{(t_{0},t_{m})}}$ is a $C^{2}$ integral curve of the flow generated by the following system of Euler-Lagrange (EL) equations:

$\displaystyle\ddot{\mathbf{x}}=\mathbf{F}^{\varepsilon}_{t}+\left(\nabla{\mathbf{F}^{\varepsilon}}-\nabla{\mathbf{F}^{\varepsilon}}^{\intercal}\right)\dot{\mathbf{x}}+\nabla{\mathbf{F}^{\varepsilon}}^{\intercal}{\mathbf{F}^{\varepsilon}},$

(17)

where the subscript $t$ denotes the partial derivative with respect to time. It is convenient to introduce the scaling $z=x/\varepsilon,$ $z\in[-1,1],$ so the EL Equations (17) become the fast-slow system

	$\displaystyle\varepsilon\ddot{z}$	$\displaystyle=F^{\varepsilon}_{1,t}+\left(\nabla_{\mathbf{y}}F^{\varepsilon}_{1}-\mathbf{G}^{\varepsilon}_{z}\right)\mathbf{\dot{y}}+\mathbf{G}^{\varepsilon\,\intercal}_{z}\mathbf{G}^{\varepsilon},$
	$\displaystyle\ddot{\mathbf{y}}$	$\displaystyle=\mathbf{G}^{\varepsilon}_{t}+\left(\mathbf{G}^{\varepsilon}_{z}-\nabla_{\mathbf{y}}F^{\varepsilon}_{1}\right)\varepsilon\dot{z}+\left(\nabla_{\mathbf{y}}\mathbf{G}^{\varepsilon}-\nabla_{\mathbf{y}}\mathbf{G}^{\varepsilon\,\intercal}\right)\dot{\mathbf{y}}+\nabla_{\mathbf{y}}\mathbf{G}^{\varepsilon\,\intercal}\mathbf{G}^{\varepsilon}.$

Introducing the conjugate momenta $\varphi=\varepsilon\dot{z}-F^{\varepsilon}_{1}$ and $\bm{\psi}=\dot{\mathbf{y}}-\mathbf{G}^{\varepsilon}$ leads to the Hamiltonian form,

	$\displaystyle\varepsilon\dot{z}$	$\displaystyle=F^{\varepsilon}_{1}+\varphi,$		(18)
	$\displaystyle\dot{\mathbf{y}}$	$\displaystyle=\mathbf{G}^{\varepsilon}+\mathbf{\bm{\psi}},$
	$\displaystyle\dot{\varphi}$	$\displaystyle=-\varphi F^{\varepsilon}_{1,z}-\langle\mathbf{\bm{\psi}},\mathbf{G}^{\varepsilon}_{z}\rangle,$
	$\displaystyle\mathbf{\dot{\bm{\psi}}}$	$\displaystyle=-\varphi\nabla_{\mathbf{y}}F^{\varepsilon}_{1}-\langle\mathbf{\bm{\psi}},\nabla_{\mathbf{y}}\mathbf{G}^{\varepsilon}\rangle.$

The most probable path of System (2),(3) corresponds to solutions of the EL Equations (18) subject to boundary conditions for appropriate values of $t$. In the case studies of Sections 4 and 5, we use the gradient flow to numerically solve System (18) when an explicit solution is not available.