# Most probable escape paths in periodically driven nonlinear oscillators

Lautaro Cilenti<EMAIL_ADDRESS>Department of Mechanical Engineering,
University of Maryland, College Park, MD 20742, USA Maria K. Cameron
<EMAIL_ADDRESS>Department of Mathematics, University of Maryland, College
Park, MD 20742, USA  Balakumar Balachandran<EMAIL_ADDRESS>Department of
Mechanical Engineering, University of Maryland, College Park, MD 20742, USA

###### Abstract

The dynamics of mechanical systems such as turbomachinery with multiple blades
are often modeled by arrays of periodically driven coupled nonlinear
oscillators. It is known that such systems may have multiple stable
vibrational modes, and transitions between them may occur under the influence
of random factors. A methodology for finding most probable escape paths and
estimating the transition rates in the small noise limit is developed and
applied to a collection of arrays of coupled monostable oscillators with cubic
nonlinearity, small damping, and harmonic external forcing. The methodology is
built upon the action plot method (Beri et al. 2005) and relies on the large
deviation theory, optimal control theory, and the Floquet theory. The action
plot method is promoted to non-autonomous high-dimensional systems, and a
method for solving the arising optimization problem with discontinuous
objective function restricted to a certain manifold is proposed. The most
probable escape paths between stable vibrational modes in arrays of up to five
oscillators and the corresponding quasipotential barriers are computed and
visualized. The dependence of the quasipotential barrier on the parameters of
the system is discussed.

> Response transitions from one dynamic state to another can occur due to
> random perturbations in a variety of systems, including mechanical and
> structural systems. Some examples are sensor arrays, energy harvesters, and
> rotating machinery. The aim of the present work is to elucidate these
> transitions, in particular, when the random perturbations are weak. To that
> end, a methodology based on Action Plot Method from Large Deviation Theory
> and Optimal Control Theory is developed and illustrated for a range of
> periodically forced systems.

##  1 Introduction

Systems consisting of multiple beam-like and plate-like structures including
turbomachinery, micro-electromechanical systems (MEMS), and vibration energy
harvesters (VEH) can exhibit stable periodic mechanical oscillations that can
be unfavorable to structural integrity and/or performance of these systems. In
this work, the authors study the effect of random pertubations (noise) in
physics-based models of these systems, as noise can induce transitions into
and out of undesired stable modes. Many of these systems can be modeled as
coupled Duffing oscillator [1] arrays with external periodic forcing,
$\bm{f}(t)$, and a weak white noise, $\sqrt{\epsilon}\bm{\eta(t)}$, accounting
for random influences, as shown below:

$\ddot{\bm{x}}+\delta_{c}\dot{\bm{x}}+\alpha\bm{x}+\beta\bm{x}^{3}+\nu\bm{D}_{N}\bm{x}=\bm{f}(t)+\sqrt{\epsilon}\bm{\eta(t)},\quad\bm{x}\in\mathbb{R}^{N}.$
(1)

Here, $N$ is the number of oscillators, $\bm{x}$ is the vector of beam tip
displacements from their respective equilibrium positions, $\delta_{c}$ is a
damping related coefficient, $\alpha$ and $\beta$ are scalar constants that
render each oscillator monostable, $\bm{x}^{3}$ denotes the elementwise cube
of the vector, $\nu$ is the inter-oscillator coupling coefficient, and
$\bm{D}_{N}$ is the coupling matrix that takes the following form:

$\displaystyle\bm{D}_{1}$
$\displaystyle=0,\quad\bm{D}_{2}=\left[\begin{array}[]{rr}1&-1\\\
-1&1\end{array}\right],$ (4) $\displaystyle\bm{D}_{N}$
$\displaystyle=\left[\begin{array}[]{rrrr}2&-1&&-1\\\ -1&2&-1&\\\
&\ddots&\ddots&\ddots\\\ -1&&-1&2\end{array}\right],\quad N>2.$ (9)

Figure 1: Schematic of linearly coupled nonlinear spring– mass-damper systems
arranged in a cyclical array.

It is noted that if the number of oscillators is greater than two, they are
coupled in a circular array configuration, as depicted in Fig. 1.

Setting the parameter $\epsilon$ in the noise term to zero, the authors obtain
the corresponding deterministic ordinary differential equations (ODE) system:

$\ddot{\bm{x}}+\delta_{c}\dot{\bm{x}}+\alpha\bm{x}+\beta\bm{x}^{3}+\nu\bm{D}_{N}\bm{x}=\bm{f}(t),\quad\bm{x}\in\mathbb{R}^{N}.$
(10)

It is known that the ODE system (10) may admit multiple attractors [2, 3, 4,
5, 6]; that is, stable periodic solutions that attract other solutions from
appropriate open subsets of the initial condition space
$(\bm{x}(0),\dot{\bm{x}}(0))$. In Fig. 2, the authors show the response of the
system (10) to the forcing $\bm{f}$ with frequency $\omega$ for a circular
array of $N=5$ coupled Duffing oscillators. Each periodic solution
$\bm{q}(t)\equiv(\bm{x}(t),\dot{\bm{x}}(t))$ to (10) is represented by a
$L_{2}$-like norm

$||\bm{q}(t)||_{L_{2}}=\sqrt{\int_{0}^{1}\sum_{i}^{2N}\left(q_{i}\left(sT\right)\right)^{2}ds},$
(11)

where $T=\nicefrac{{2\pi}}{{\omega}}$ is the period.

Figure 2: Frequency response of a circular array of $N=5$ coupled forced
Duffing oscillators as a function of excitation frequency, $\omega$. At the
frequency $\omega=1.4$, unique attractors of the array are highlighted. This
figure was generated by using the continuation package AUTO2007 [7]. The
parameters in the ODE system (10) were set to $\alpha=1$, $\beta=0.3$,
$\delta=0.1$, $\nu=0.01$, and the periodic excitation was of the form
$f(t)=0.4\cos(\omega t)$. Solid (dotted) lines in the figure correspond to
stable (unstable) periodic solutions of the system (10). Here, it is of
interest to examine how random perturbations influence the transition from one
stable solution to another.

Each attractor is encoded by a sequence of symbols
$X_{1}X_{2}X_{3}X_{4}X_{5}$, $X_{j}\in\\{\rm{L},\rm{H}\\}$, $j=1,\ldots,5$,
where L and H stand, respectively, for “low amplitude” and “high amplitude”.
For example, LHLLL refers to the stable periodic solution where oscillator 2
has high amplitude while the rest of the oscillators have low amplitudes. From
Fig. 2, it is discernible that if the frequency $\omega$ of the external
forcing is far from the primary resonance, there exists at least one stable
periodic solution. Furthermore, at a certain range of
$\omega_{1}<\omega<\omega_{2}$, there are a total of at least 32 periodic
attractors corresponding to which each oscillator is in either a low or high
amplitude state. This region with multiple stable solutions is referred to as
the hysteresis region of the frequency response. In this region, the state
(stable response) of the system is sensitive to initial (prior) conditions.
Note that some solution branches in Fig. 2 coincide due to circular symmetry
or closeness of their representations by $||\bm{q}(t)||_{L_{2}}$. The
phenomenon wherein a smaller subset of the system (one or more oscillators but
not all) consists of high amplitude responses is called _localization_. In
this work, the authors refer to _response localization_ exclusively as the
stable response mode in which only _one_ oscillator oscillates with a high
amplitude while all of other oscillators oscillate with low amplitudes.

Response localizations have been demonstrated in experiments of macro-
cantilever beam arrays, wherein the beams represent the blades or components
of turbomachinery and VEHs [8, 9, 10, 11, 12]. The localized high amplitude
oscillations can, on one hand, be destructive in turbomachinery [2, 13, 14,
15, 16] and, on the other hand, be beneficial for energy harvesting systems
[17, 18].

It is remarked that localization in the system governed by (10) is a nonlinear
phenomenon: it would not occur if $\beta=0$. However, localization may occur
also in linear systems due to inhomogeneity [16, 19, 13, 20, 17]; that is, if
$\beta=0$ and $\alpha$ is replaced with the vector
$(\alpha_{1},\ldots,\alpha_{N})$. In the systems with $N=3$ and $N=5$, here,
all stable response localizations have symmetric, anti-phase, neighboring
oscillators.

Adding a stochastic perturbation (noise) to an array of coupled Duffing
oscillators; that is, replacing ODE (10) with (1) can facilitate transitions
between various attractors (modes) of these systems [11]. Therefore, in an
application, noise may be used to transition a system out of an undesirable
localized mode in a bladed rotor or may help a VEH to get into the desirable
high amplitude regime.

##  2 Objectives and Structure

Figure 3:  A flowchart overview of the main concepts in this work. The authors
begin by presenting background on a noise influenced system model, large
deviation theory, and optimal control theory. A methodology built on the
action plot method combined with the Floquet theory and numerical optimization
is used to compute most probable escape paths and the quasipotential barriers
of escapes from attractors of the system model.

The authors achieve the following two objectives here.

  1. 1.

Conduct a case study on the noise-induced most probable escape paths out of
undesired oscillations modes such as response localizations that exist only
under the influence of deterministic periodic excitations in circular arrays
of nonlinear oscillators.

  2. 2.

Develop a computational tool for accomplishing Objective 1; that is, for
finding the quasipotential barriers and most probable escape paths of
nonlinear systems with multiple attractors subjected to both deterministic
periodic excitations and weak random perturbations.

By achieving Objective 1, the authors contribute to a highly active area of
research and motivate Objective 2, since there is no alternative methodology
to compute these escape paths and quasipotential barriers for these high
dimensional systems with low damping and deterministic time-dependent external
excitations to the best of the authors’ knowledge.

The focus of this work is on the weak white noise limit; that is the limit
$\epsilon\to 0$. In a nonlinear system with multiple attractors, the
probability of a transition between its attractors is nonzero no matter how
small the white noise intensity is [21]. The possibility of these rare events
is interesting when the timescale of deterministic dynamics is different from
another timescale of interest. With regard to turbomachinery, one can imagine
blades that may vibrate with sub-second periods, but the timescale of
continuous operation of the system can range from hours to years. Given this
duration of operation, rare events become relevant. Furthermore, the
asymptotic estimates for the weak noise limit often remain relevant for a
broad range of noise intensities.

The current work relies on four building blocks:

  1. 1.

The authors use Freidlin’s and Wentzell’s large deviations theory (LDT) [21]
to define the action functional (FW action) whose minimizers are the most
probable escape paths from basins of attractors of (10). The corresponding
minimal values of the action, which are the values of a certain function
called the _quasipotential_ , serve as a quantitative measure for the
difficulty of escape from an attractor and the consequent transition.

  2. 2.

The optimal control theory and Pontryagin’s maximum principle [22] are used to
identify the set of paths called _Hamiltonian paths_ , which are the
candidates for being the minimizers of the FW action. As a result, the problem
of minimization of the FW action reduces to an optimization problem on the set
of initial conditions for the Hamiltonian paths.

  3. 3.

The action plot method [23] provides a nice recipe for cutting the
dimensionality of this optimization problem by a factor of two through
linearization of the system and approximation of the _Lagrangian manifold_
containing the desired minimizing initial conditions.

  4. 4.

Following Ref. [24], the authors use the Floquet theory for linear ODEs with
periodic coefficients [25] to adapt this recipe for periodic attractors of
non-autonomous stochastic differential equations (SDEs) (1).

The use of these building blocks is illustrated via the flow chart in Fig. 3.

The resulting optimization problem for the optimal set of initial conditions
for the Hamiltonian paths is still very challenging because the objective
function is discontinuous and this function’s minima lie next to its “cliffs”.
In the case of a $2$-dimensional system, i.e., only one oscillator, this
problem can be solved graphically as it has been done in references [23, 26].
For the higher dimensional problems considered in this work, the authors
propose an optimization algorithm _stochastic unit vector (SUV)_ for finding
the minimizers of the FW action.

First, the developed methodology is applied to a single Duffing oscillator to
quantify transitions between this system’s low amplitude and high amplitude
attractors. This example is highly instructive as it is visual and illuminates
the technical difficulties for $N>1$. Next, the authors compute the most
probable escape paths between attractors of arrays consisting of two, three,
and five coupled Duffing oscillators. Finally, the dependencies of the
quasipotential values characterizing the difficulty of escapes from the basins
of attractors on the external forcing frequency $\omega$ and inter-oscillator
coupling coefficient $\nu$ are visualized. These plots are helpful in
validating the obtained results.

The rest of the paper is organized as follows. In Section 3, the authors
discuss related work in the cases of relatively large noise intensity and the
weak noise limit. In Section 4, they lay out fundamental background from large
deviations theory and optimal control theory. The proposed methodology is
detailed in 5. Results obtained for a single forced Duffing oscillator, two
coupled forced Duffing oscillators, and two circular arrays of $N=3$ and $N=5$
Duffing oscillators are described in Sections 6, 7, and 8, respectively. The
paper is concluded with a discussion in Section 9. Details and parameter
values for the computations conducted in this work are found in Appendix A.
Table LABEL:tab:NotationTable in Appendix A is provided as a nomenclature
reference for a reader.

##  3 Related work

Nonlinear oscillators with external periodic forcing and added noise have been
a subject of active research for at least the last two decades.

###  3.1 Relatively large noise

Experimental and numerical studies conducted by using methods limited to
_sufficiently large noise intensity_ have yielded a qualitative understanding
of the dynamics of these systems. Approximations to the stationary probability
density response of mono and coupled second-order nonlinear systems were
constructed in references [27, 28, 29, 30]. The method based on the
construction of the _prehistory probability distributions_ [31] was applied to
a single chaotic oscillator to find the most probable escape paths from basins
of its attractors [32]. The average dynamics of the nonlinear systems was
approximated by deriving a series of nonlinear moment equations from the
Fokker-Planck equation and integrating them [33, 34, 35], and by means of the
cell mapping method [36, 37]. Transient probability density functions for
forced mono and coupled Duffing oscillators were generated by utilizing new
developments in the path integral method [38, 39, 40, 41]. A comparison
between simulations of circular arrays of six and more oscillators revealed no
fundamental differences with regards to the transitions from a response
localization to another stable mode [42].

Experimental studies of the effect of noise on response localizations [11, 43,
44] have demonstrated that likelihood of transitions from one stable mode to
another depends on the deterministic excitation frequency of the forced
system. Although all of these studies provide relevant insights into the
probabilistic behavior of the considered stochastic nonlinear mechanical
systems, the methodologies used in the earlier studies are not amenable for
studying the dynamics in the weak noise limit.

###  3.2 Weak noise limit

The present work is devoted to the study of quasipotential barriers and the
most probable escape paths (MPEPs) from basins of various attractors of arrays
of periodically forced Duffing oscillators in the limit of the white noise
intensity tending to zero. Central to the developed methodology is the _action
plot method_ introduced in the work of Beri et al. [23] for finding MPEPs in
two-dimensional nongradient SDEs. An important advantage of the action plot
method is that it is suitable for systems with inertia and external forcing.
For example, Chen et al. [26] used this method to find a MPEP for a chaotic
oscillator [32] and validated it by matching it with the computed prehistory
probability distribution. The action plot method implies a graphical solution
to the optimization problem that arises – this is the reason for its name –
and hence this approach in its original form [23] is limited to low
dimensional systems (single oscillator). The main methodological contribution
of the present work is the extension of the action plot method to high
dimensional systems with periodic forcing and the development of a
minimization technique for solving the optimization problem for initial
conditions for the MPEP.

In addition, the authors would like to highlight a recent work [45] that
offers an asymptotic technique for computing the quasipotential barriers in a
bistable multidimensional oscillator with small noise and small periodic
forcing based on an asymptotic expansion of the FW action in powers of the
scaling parameter used for the forcing. While this approach is feasible for
high dimensions, it is not suitable for the oscillator arrays studied in the
present work as they are monostable in the absence of the periodic forcing.
Moreover, the periodic forcing is not assumed to be small in this work.

##  4 Background

###  4.1 Large Deviations Theory

In this section, the authors provide a brief background from large deviations
theory, which is used in this work to define the quasipotential.

The Freidlin-Wentzell large deviations theory [21] primarily deals with
_autonomous_ SDEs with _non-degenerate_ noise

$d\bm{q}=\bm{b}(\bm{q})dt+\sqrt{\varepsilon}\,\bm{G}(\bm{q})\,d\bm{W}(t),$
(12)

where $d{\bf W}(t)$ is the standard Brownian motion. The drift field $\bm{b}$
and the diffusion matrix $\bm{Q}(\bm{q})=\bm{G}\bm{G}^{\top}$ are assumed to
be continuously differentiable.

The work performed by random perturbations can be quantified with the action
functional $S_{\tau}(\bm{\varphi})$. When $\bm{Q}$ is invertible; that is, the
noise in equation (12) is non-degenerate, the Freidlin-Wentzell action is
defined on the set of absolutely continuous paths for $\bm{\varphi}$ by [21]

$S_{\tau}(\bm{\varphi})=\frac{1}{2}\int_{0}^{\tau}\left\|\dot{\bm{\varphi}}-\bm{b}(\bm{\varphi})\right\|^{2}_{\bm{Q}^{-1}(\bm{\varphi})}dt,$
(13)

where $\|\cdot\|^{2}_{\bm{Q}^{-1}}$ denotes the inner product
$\langle\cdot,\bm{Q}^{-1}\cdot\rangle$. One takeaway from (13) is that the
action is zero if the path $\bm{\varphi}$ follows a deterministic trajectory
of the corresponding ODE system

$\dot{\bm{q}}=\bm{b}(\bm{q}).$ (14)

The probability for a stochastic trajectory of SDEs (12) to follow a small
$\delta$-tube surrounding a given absolutely continuous path $\bm{\varphi}$
scales with $\epsilon$ as $\exp\\{-S_{\tau}(\bm{\varphi})/\epsilon\\}$ (see
Refs. [21, 46]).

The _quasipotential $U_{\mathcal{A}}(\bm{q})$ at a point $\bm{q}$ with respect
to an attractor_ $\mathcal{A}$ of (14) is defined as the infimum of the action
taken over all possible absolutely continuous paths starting at $\mathcal{A}$
and ending at $\bm{q}$ and all travel times $\tau$:

$U_{\mathcal{A}}(\bm{q})=\inf_{\bm{\varphi},\tau}\\{S_{\tau}(\bm{\varphi})\,|\,\bm{\varphi}(0)\in\mathcal{A},\leavevmode\nobreak\
\bm{\varphi}(\tau)=\bm{q}\\}.$ (15)

Note that this infimum is achieved at $\tau=\infty$ as it takes infinitely
long time to escape from the attractor in the weak noise limit. Furthermore,
if the attractor is not a point attractor but a limit cycle or a more complex
object, the infimum is achieved via an escape path of infinite length, as it
makes infinitely many revolutions in neighborhood of the attractor while
escaping from it. This minimizing path is called the _minimum action path_ or
_MAP_.

Let $D$ be a domain with a smooth boundary $\partial D$ lying in the basin of
$\mathcal{A}$. The quasipotential allows us to estimate the expected exit time
$\tau_{exit}$ from $D$ up to the exponential order [21, 46]:

$\lim_{\varepsilon\to 0}\varepsilon\textrm{ log
}E\left[\tau_{exit}\right]=\inf_{\bm{q}\in\partial D}U_{\mathcal{A}}(\bm{q})$
(16)

Therefore, in the limit $\epsilon\rightarrow 0$, the system escapes the domain
$D$ after a waiting time logarithmically equivalent to
$\exp\left\\{-\inf_{\bm{q}\in\partial
D}U_{\mathcal{A}}(\bm{q})/\epsilon\right\\}$. Moreover, the escape path lies
in a small tube surrounding the minimizer of the FW action (13), i.e. the MAP,
with probability tending to one as $\epsilon\rightarrow 0$. Consequently, this
MAP is often referred to as the _most probable escape path (MPEP)_. In the
present work, the authors study the escape from the basin of attraction of
$\mathcal{A}$. Therefore $D$ is chosen to be the basin of attraction and
$\partial D$ is the basin boundary.

The generalization of these results of LDT for systems with inertia has been
addressed only partially in the mathematical literature, to the best of the
authors’ knowledge. This has been done for the dynamics governed by
$\ddot{\bm{x}}+\delta_{c}\dot{\bm{x}}+\bm{b}(\bm{x})=\sqrt{\varepsilon}\bm{\dot{W}}(t)$
for the case where $\bm{b}$ is either uniformly Lipschitz [47] or a gradient
of a smooth potential function [48] $\bm{b}=-\nabla\bm{V}$. The conjecture
that similar results for the MPEPs and the expected exit times hold for
oscillators of the form (1) has been in physical literature since at least
1979 – see reference [49]. This conjecture is rooted in Feynman’s path
integral theory [50]. A rigorous proof of an estimate for the exit time
similar to (16) for the dynamics governed by SDEs (1) is yet to be given to
the best of the authors’ knowledge.

Finally, it worth mentioning that there is strong numerical evidence
suggesting that the MPEP for a nonlinear oscillator with external periodic
forcing is the MAP obtained by matching the MAP with the plot of the
prehistory probability distribution [26].

###  4.2 Definitions for periodically driven nonlinear oscillators

Now consider SDE (1) with periodic forcing $\bm{f}(t)$ of period $T$ and the
corresponding ODE (10). There is an important difference between a periodic
solution of an autonomous ODE and a periodic solution of the time-dependent or
non-autonomous ODE (10). Any time shift of a periodic solution of an
autonomous ODE gives another periodic solution of the same ODE. In contrast,
only time shifts of periodic solutions to (10) that are multiples of the
period result in other periodic solutions to (10). Therefore, every point
$(\bm{x},\bm{v})$ of a periodic solution to (10) is associated with a
particular phase $\theta:=t\thinspace{\rm mod}\thinspace T$. Hence, it is
natural to make ODE (10) autonomous by introducing the phase variable

$\theta:=t\thinspace{\rm mod}\thinspace
T\in\mathbb{S}^{1}_{T}\equiv\mathbb{S}^{1}(T/{2\pi}),$

where $\mathbb{S}^{1}_{T}$ is a circle of radius ${T}/{2\pi}$. The resulting
autonomous ODE on the manifold $\mathbb{R}^{2N}\times\mathbb{S}^{1}_{T}$ is:

$\displaystyle\dot{\bm{x}}$ $\displaystyle=\bm{v},$
$\displaystyle\dot{\bm{v}}$
$\displaystyle=-\delta_{c}\bm{v}-\bm{K}(\bm{x})+\bm{f}(\theta),$ (17)
$\displaystyle\dot{\theta}$ $\displaystyle=1,$ $\displaystyle\bm{x}$
$\displaystyle\in\mathbb{R}^{N},\quad\bm{v}\in\mathbb{R}^{N},\quad\theta\in\mathbb{S}^{1}_{T}.$

Here, the notation
$\bm{K}(x):=\alpha\bm{x}+\beta\bm{x}^{3}+\nu\bm{D}_{N}\bm{x}$ is introduced
for brevity and convenience of generalization. The attracting periodic
trajectories of the ODE (10) are the attractors of (17).

The corresponding SDE (1) is made autonomous in the same way:

$\begin{split}&d\bm{x}=\bm{v}\,dt\\\
&d\bm{v}=\left(-\delta_{c}\bm{v}-\bm{K}(\bm{x})+\bm{f}(\theta)\right)dt+\sqrt{\varepsilon}\,d{\bm{W}}(t)\\\
&d\theta=dt\\\
&\bm{x}\in\mathbb{R}^{N},\quad\bm{v}\in\mathbb{R}^{N},\quad\theta\in\mathbb{S}^{1}_{T}.\\\
\end{split}$ (18)

The action functional for SDE (18) on the set of absolutely continuous paths
$\\{[\bm{\varphi}(t),\theta(t)]\leavevmode\nobreak\ |\leavevmode\nobreak\
t_{0}\leq t\leq\tau\\}\subset\mathbb{R}^{N}\times\mathbb{S}^{1}_{T}$ is
defined as

$\begin{split}S_{[t_{0},\tau]}([\bm{\varphi},\theta])&=\frac{1}{2}\int_{t_{0}}^{\tau}||\ddot{\bm{\varphi}}+\delta_{c}\dot{\bm{\varphi}}+\bm{K}(\bm{\varphi})-\bm{f}(\theta)||^{2}dt\\\
\end{split}$ (19)

According to Dykman’s and Krivoglaz’s conjecture [49], which is consistent
with Feyman’s path integral theory [50], and supported by a strong numerical
evidence [26], the minimizers of (19) have the following significance. Let
$[\bm{\varphi}^{\ast}(t),\theta^{\ast}(t)]$, $t_{0}\leq t\leq\tau$, be a
minimizer of (19) amongst all paths satisfying
$[\bm{\varphi}(t_{0}),\theta(t_{0})]=[\bm{\varphi}^{\ast}(t_{0}),\theta^{\ast}(t_{0})]$
and
$[\bm{\varphi}(\tau),\theta(\tau)]=[\bm{\varphi}^{\ast}(\tau),\theta^{\ast}(\tau)]$.
Then the probability for a trajectory of (18) starting at
$[\bm{\varphi}^{\ast}(t_{0}),\theta^{\ast}(t_{0})]$ and ending at
$[\bm{\varphi}^{\ast}(\tau),\theta^{\ast}(\tau)]$ to follow a small tube
around a path connecting these points is maximized if this path a minimizer of
(19). Moreover, this probability scales with the parameter $\epsilon$ in (18)
as $\exp\\{-S_{[t_{0},\tau]}([\bm{\phi}^{\ast},\theta^{\ast}])/\epsilon\\}$.
Therefore, in order to find the most probable escape path from a basin of
attractor $\hat{\bm{A}}_{0}$ of (17), one should minimize (19) with respect to
all paths and all times $t_{0}\in[0,T)$ and $\tau\in[0,\infty]$ satisfying the
boundary conditions $[\bm{\varphi}(t_{0}),\theta(t_{0})]\in\hat{\bm{A}}_{0}$,
$[\bm{\varphi}(\tau),\theta(\tau)]\in\partial B_{0}$ where $B_{0}$ denotes the
basin of $\hat{\bm{A}}_{0}$.

###  4.3 Pontryagin’s Maximum Principle

Once the action functional (the cost functional) is known, one can use optimal
control theory [22, 51] to simplify the problem of minimizing it. The relevant
ODE with an $N-dimensional$ control function $\bm{u}(t)$ is given by

$\displaystyle\dot{\bm{x}}$ $\displaystyle=\bm{v},$
$\displaystyle\dot{\bm{v}}$
$\displaystyle=-\delta_{c}\bm{v}-\bm{K}(\bm{x})+\bm{f}(\theta)+\bm{u}(t),$
(20) $\displaystyle\dot{\theta}$ $\displaystyle=1,$ $\displaystyle\bm{x}$
$\displaystyle\in\mathbb{R}^{N},\quad\bm{v}\in\mathbb{R}^{N},\quad\theta\in\mathbb{S}^{1}_{T}.$

Expressing $\bm{u}$ from the equation for $\dot{\bm{v}}$ and comparing it to
the action functional (19) we obtain the cost functional

$S(\bm{u})=\frac{1}{2}\int_{t_{0}}^{\tau}\|\bm{u}\|^{2}\,dt.$ (21)

The optimal control problem consists in finding a function $\bm{u}$ that
minimizes the cost functional subject to some user-chosen constraints.

Suppose that ODEs (17) admit $n_{a}>1$ attractors, $\hat{\bm{A}}_{k}$, $0\leq
k\leq n_{a}-1$. The authors denote the basin of $\hat{\bm{A}}_{k}$ by $B_{k}$
and the basin boundary by $\partial B_{k}$, $0\leq k\leq n_{a}-1$. Index 0 is
reserved for the initial attractor, from which the most probable escape path
emerges. Thus, the constraints are set as

$(\bm{x}(t_{0}),\bm{v}(t_{0}),\theta(t_{0}))\in\hat{\bm{A}}_{0},\quad(\bm{x}(\tau),\bm{v}(\tau),\theta(\tau))\in\partial
B_{0}.$ (22)

The optimal control problem

$\displaystyle
S(\bm{u})=\frac{1}{2}\int_{t_{0}}^{\tau}\|\bm{u}\|^{2}\,dt\rightarrow\min$
(23) subject to (20) and (22)

will be solved with the aid of _Pontryagin’s Maximum Principle_ [22] whose
basic form is the following. Suppose a controlled system is governed by a
system of ODEs

$\dot{\bm{q}}=\bm{b}(\bm{q},\bm{u}).$ (24)

Let the cost functional be

$S(\bm{u})=\int_{t_{0}}^{t_{1}}L(\bm{q},\bm{u})\,dt.$ (25)

Then, the _control theory Hamiltonian_ is defined by

$H(\bm{q},\bm{u},\bm{p}):=\bm{p}^{\top}\bm{b}(\bm{q},\bm{u})-L(\bm{q},\bm{u}),$
(26)

where $\bm{p}$ is a costate variable. _Pontryagin’s maximum principle_ states
that if $\bm{u}^{\ast}(t)$ is the minimizer of (25) and $\bm{q}^{\ast}(t)$ is
the corresponding optimal trajectory of (24) then there exists a function
$\bm{p}^{\ast}(t)$, that, together with $\bm{q}^{\ast}(t)$ satisfies the
Hamilton canonical equations [22]

$\left(\begin{array}[]{c}\dot{\bm{q}}\\\
\dot{\bm{p}}\end{array}\right)=\left(\begin{array}[]{r}\nabla_{\bm{p}}H\\\
-\nabla_{\bm{q}}H\end{array}\right).$ (27)

The costate variables $\bm{p}$ are the momenta in the context of classical
mechanics [52] but they do not have that physical meaning in the current
problem.

Furthermore, according to Pontryagin’s maximum principle, the Hamiltonian is
identically equal to zero and achieves its maximum along the optimal
trajectory $(\bm{q}^{*},\bm{p}^{*})$ of (27):

$0=H(\bm{q}^{*},\bm{u}^{*},\bm{p}^{*})=\max_{\bm{u}\in\mathbb{U}}H(\bm{q},\bm{u},\bm{p}),$
(28)

where $\mathbb{U}$ is the set of admissible controls.

For the control problem (23), the state and costate variables are

$\hat{\bm{q}}\equiv\left(\begin{array}[]{c}\bm{q}\\\
\theta\end{array}\right)=\left(\begin{array}[]{c}\bm{x}\\\ \bm{v}\\\
\theta\end{array}\right)\quad{\rm
and}\quad\hat{\bm{p}}\equiv\left(\begin{array}[]{c}\bm{p}\\\
p_{\theta}\end{array}\right)=\left(\begin{array}[]{c}\bm{p}_{\bm{x}}\\\
\bm{p}_{\bm{v}}\\\ p_{\theta}\end{array}\right).$

The optimal control $\bm{u}^{\ast}$ can be found using the fact that the
Hamiltonian achieves its maximum at the optimal control $\bm{u}^{*}$.
Therefore

$\displaystyle 0=$
$\displaystyle\nabla_{\bm{u}}H(\hat{\bm{q}},\bm{u}^{\ast},\hat{\bm{p}})$
$\displaystyle=$
$\displaystyle\nabla_{\bm{u}}\left(\bm{p}_{\bm{x}}^{\top}\bm{v}+\bm{p}_{\bm{v}}^{\top}\left[-\delta_{c}\bm{v}-\bm{K}(\bm{x})+\bm{f}(\theta)+\bm{u}^{\ast}\right]\right.$
$\displaystyle\left.+p_{\theta}-\frac{1}{2}\|\bm{u}^{\ast}\|^{2}\right)$
$\displaystyle=$ $\displaystyle\bm{p}_{\bm{v}}-\bm{u}^{\ast}.$ (29)

Hence $\bm{u}^{\ast}=\bm{p}_{\bm{v}}^{\ast}$ and the optimal path satisfies
the Hamiltonian dynamical system

$\left(\begin{array}[]{c}\dot{\bm{x}}\\\ \dot{\bm{v}}\\\ \dot{\theta}\\\
\dot{\bm{p_{x}}}\\\ \dot{\bm{p_{v}}}\\\
\dot{p}_{\theta}\end{array}\right)=\left(\begin{array}[]{c}\bm{v}\\\
-\delta_{c}\bm{v}-\bm{K}(\bm{x})+\bm{f}(\theta)+\bm{p_{v}}\\\ 1\\\
\nabla_{\bm{x}}\left(\bm{p_{v}}^{\top}\bm{K}(\bm{x})\right)\\\
-\bm{p_{x}}+\delta_{c}\bm{p_{v}}\\\
-\bm{p}_{\bm{v}}^{\top}\bm{f}^{\prime}(\theta)\end{array}\right).$ (30)

Any attractor $\hat{\bm{A}}_{k}$ of the ODEs (17) can be mapped onto a saddle-
type invariant set $\bar{\bm{A}}_{k}$ of the corresponding Hamiltonian ODEs
(30) by including the costate variables $\hat{\bm{p}}$ and setting them to
zero. Since the first boundary condition in (22) is set at the attractor
$\hat{\bm{A}}_{0}$, the corresponding trajectory of (30) passing through
$(\hat{\bm{q}}(t_{0}),\bm{0})$ will stay forever at the corresponding saddle-
type invariant set $\bar{\bm{A}}_{0}$. On the other hand, a whole family of
trajectories called the _Hamiltonian paths_ emanates from $\bar{\bm{A}}_{0}$
at time $t_{0}=-\infty$ and constitutes the $(2N+1)$-dimensional unstable
manifold of $\bar{\bm{A}}_{0}$ called the _Lagrangian manifold_. This will be
clear from the structure of the Jacobian matrix obtained in Section 5.1.

Therefore, the optimal control problem (23) is reduced to finding the initial
condition for the optimal Hamiltonian path at some finite moment in time. This
is done by surrounding the attractor $\hat{\bm{A}}_{0}$ with a small
neighborhood and searching for the initial conditions for the optimal
Hamiltonian path in the intersection of the boundary of this neighborhood and
the Largangian manifold. This construction is elaborated in Section 5.1.

##  5 Methodology

The methodology is developed for SDEs of the form (18) where the external
forcing $\bm{f}(t)$ is periodic. It is assumed that the corresponding ODEs
(17) admit $n_{a}\geq 2$ attractors $\hat{\bm{A}}_{0}$, $\ldots$,
$\hat{\bm{A}}_{n_{a}-1}$ with basins $B_{0}$, $\ldots$, $B_{n_{a}-1}$, and
these _attractors are periodic orbits_. The intersections of basin boundaries
$\partial B_{k}$ and $\partial B_{l}$ are denoted by $\partial B_{kl}$, $0\leq
k,l\leq n_{a}-1$.

Let the system governed by SDEs (18) be initially located near attractor
$\hat{\bm{A}}_{0}$. The implication of Pontryagin’s maximum principle is that
the minimizer of the action functional (19) is a Hamiltonian path given by
(30). As soon as the Hamiltonian path reaches the basin boundary $\partial
B_{0}$, the authors reset and fix the costates $\bm{p_{x}}$ and $\bm{p_{v}}$
to zero. From that point, the path continues a trajectory of (17) restricted
to basin boundary $\partial B_{0}$. This trajectory approaches an attractor
$S$ for the dynamics restricted to $\partial B_{0}$. For the systems studied
in this work, $S$ is always a saddle cycle. Let $S\equiv S_{0k}$ belongs to
$\partial B_{0k}$ separating basins of $\hat{\bm{A}}_{0}$ and
$\hat{\bm{A}}_{k}$, $k\neq 0$. Then, one can complete the transition path from
$\hat{\bm{A}}_{0}$ to $\hat{\bm{A}}_{k}$ by adding a trajectory going from
$S_{0k}$ to $\hat{\bm{A}}_{k}$. Since the contribution to action (13) along
any trajectory is zero, the total action along the transition path from
$\hat{\bm{A}}_{0}$ to $\hat{\bm{A}}_{k}$ is equal to the action along the
Hamiltonian path. The authors’ numerical experiments suggest (see Figs. 14,
15, 20, 22, 24, 27) that the maximum likelihood transition path from
$\hat{\bm{A}}_{0}$ to $\hat{\bm{A}}_{k}$ is a concatenation of $(i)$ a
Hamiltonian path that approaches $S_{0k}$ that is infinitely long, $(ii)$ the
saddle periodic orbit $S_{0k}$, and $(iii)$ a trajectory of (17) starting
infinitely close to $S_{0k}$ and approaching $\hat{\bm{A}}_{k}$. The key
challenge is that the initial condition for the optimal Hamiltonian path is
unknown.

###  5.1 The Lagrangian Manifold

The action plot method [23] utilizes Pontryagin’s maximum principle to reduce
the problem of finding the most probable escape path (MPEP) to the problem of
finding the initial condition for the Hamiltonian path and offers an elegant
way to simplify this problem. This is done by linearizing the Hamiltonian
system near an invariant set corresponding to an attractor of the original ODE
and approximating the unstable Lagrangian manifold by a linear relationship
between $\hat{\bm{q}}(0)$ and $\hat{\bm{p}}(0)$. The construction of the
Lagrangian manifold adapted for the case of periodic forcing is detailed
below. The Floquet theory [25] for linear ODEs with periodic coefficients is
used to facilitate finding the unstable manifold. The latter is inspired by
Ref. [24].

The equations for $\theta$ and $p_{\theta}$ in the Hamiltonian ODE (30) can be
omitted, as they do not affect the other equations. The resulting non-
autonomous Hamiltonian ODE system is given by

$\left(\begin{array}[]{c}\dot{\bm{x}}\\\ \dot{\bm{v}}\\\ \dot{\bm{p_{x}}}\\\
\dot{\bm{p_{v}}}\\\ \end{array}\right)=\left(\begin{array}[]{c}\bm{v}\\\
-\delta_{c}\bm{v}-\bm{K}(\bm{x})+\bm{f}(t)+\bm{p_{v}}\\\
\nabla_{\bm{x}}\left(\bm{p_{v}}^{\top}\bm{K}(\bm{x})\right)\\\
-\bm{p_{x}}+\delta_{c}\bm{p_{v}}\end{array}\right).$ (31)

The periodic trajectories of (30) are naturally mapped onto periodic
trajectories of (31) by dropping of the components $\theta$ and $p_{\theta}$.
Let ${\hat{\bm{A}}}_{0}=:(\bm{x}(t),\bm{v}(t),\theta(t))$ be an attractor of
the ODE system (17). It corresponds to a saddle-type periodic trajectory
$\bar{\bm{A}}_{0}$ of the Hamiltonian dynamical system (30). In turn,
$\bar{\bm{A}}_{0}$ corresponds to a periodic solution
$\tilde{\bm{A}}_{0}:=(\bm{x}(t),\bm{v}(t),\bm{0},\bm{0}))$ of the non-
autonomous Hamiltonian dynamical system (31). The authors’ goal is to find the
set of the initial conditions for the Hamiltonian paths forming the Lagrangian
manifold. To do so, one needs to linearize ODEs (31) near the periodic
solution $\tilde{\bm{A}}_{0}$. The Jacobian of the right-hand side of (31) is

$\begin{split}&J(\bm{x},\bm{v},\bm{p}_{\bm{x}},\bm{p}_{\bm{v}},t_{0})\\\
=&\left(\begin{array}[]{cccc}\bm{0}&\bm{I}&\bm{0}&\bm{0}\\\
-\bm{J}_{\bm{K}}(\bm{x})&-\delta_{c}\bm{I}&\bm{0}&\bm{I}\\\
\nabla_{\bm{x}}\nabla_{\bm{x}}\left[\bm{p}_{\bm{v}}^{\top}\bm{K}(\bm{x})\right]&\bm{0}&\bm{0}&\bm{J}_{\bm{K}}(\bm{x})^{\top}\\\
\bm{0}&\bm{0}&-\bm{I}&\delta_{c}\bm{I}\\\ \end{array}\right),\\\ \end{split}$
(32)

and its value at $(\hat{\bm{x}}(t),\hat{\bm{v}}(t),\bm{0},\bm{0}))$ is

$\begin{split}&J(\hat{\bm{x}},\hat{\bm{v}},\bm{0},\bm{0},t_{0})\\\
=&\left(\begin{array}[]{cccc}\bm{0}&\bm{I}&\bm{0}&\bm{0}\\\
-\bm{J}_{\bm{K}}(\hat{\bm{x}})&-\delta_{c}\bm{I}&\bm{0}&\bm{I}\\\
\bm{0}&\bm{0}&\bm{0}&\bm{J}_{\bm{K}}(\hat{\bm{x}})^{\top}\\\
\bm{0}&\bm{0}&-\bm{I}&\delta_{c}\bm{I}\\\ \end{array}\right)\equiv
J_{\bm{\tilde{A}}_{0}}(t_{0}).\end{split}$ (33)

In equations (32) and (33), $\bm{J}_{\bm{K}}$ denotes the Jacobian matrix of
the vector field $\bm{K}$. The linear ODE with periodic coefficients of period
$T$ approximating (31) near $\tilde{\bm{A}}_{0}$ is of the form

$\dot{\bm{z}}=J_{\bm{\tilde{A}}_{0}}(t_{0})\bm{z}.$ (34)

Equation (34) has a fundamental solution matrix $\bm{\Psi}(t)$. According to
Floquet’s theorem [25], the monodromy matrix

$\bm{\Psi}(t_{0}+T)\bm{\Psi}(t_{0})^{-1}$ (35)

provides a coordinate change that results in the Poincaré mapping for the
corresponding linearized autonomous system with
$\theta_{0}=t_{0}\thinspace{\rm mod}\thinspace T\in\mathbb{S}^{1}_{T}$. The
eigenvalues of the monodromy matrix, called _characteristic multipliers_ , are
intrinsic to the ODE, as they remain invariant for all $t_{0}\in\mathbb{R}$.
The corresponding matrix of eigenvectors is a periodic function of $t_{0}$.
The structure of the matrix $J_{\bm{\tilde{A}}_{0}}(t_{0})$ implies that the
monodromy matrix has $2N$ eigenvalues less than 1 (stable) and $2N$
eigenvalues greater than 1 (unstable). The $4N\times 2N$ matrix functions
$\bm{\Phi}^{s}_{\theta_{0}}$ and $\bm{\Phi}^{u}_{\theta_{0}}$ are formed,
respectively, by the stable and unstable eigenvectors. At each fixed
$\theta_{0}\in[0,T]$, the span of $\bm{\Phi}^{u}_{\theta_{0}}$ is the linear
unstable manifold. Since ODE (34) describes the time evolution of the first-
order perturbation to the periodic solution $\tilde{\bm{A}}_{0}$, the
Lagrangian manifold near $\tilde{\bm{A}}_{0}$ is approximated by the sum of
trajectories of the Hamiltonian ODE (31) passing through the set of initial
conditions

$\left\\{\tilde{\bm{A}}_{0}(\theta_{0})+\hat{\gamma}\bm{\Phi}^{u}_{\theta_{0}}\bm{c}\,\leavevmode\nobreak\
\vline\leavevmode\nobreak\
\,\bm{c}\in\mathbb{R}^{2N},\theta_{0}\in\mathbb{S}^{1}_{T}\right\\},$ (36)

where $\hat{\gamma}>0$ is a scaling parameter. It is worth noting that the
choice of $\hat{\gamma}$ would be arbitrary if equation (31) were linear, and
$\theta_{0}$ could be set to zero if equation (31) were autonomous. However,
since (31) is nonlinear and nonautonomous, $\hat{\gamma}|\bm{c}|$ needs to be
small and all $\theta_{0}\in[0,T)$ should be considered.

Any escape path in $\mathbb{R}^{2N}\times\mathbb{S}^{1}$ leaving
$\hat{\bm{A}}_{0}$ must cross through a torus that surrounds
$\hat{\bm{A}}_{0}$. For each fixed $\theta_{0}\in\mathbb{S}^{1}_{T}$, the
point $\hat{\bm{A}}_{0}(\theta_{0})$ can be completely surrounded by a
$2N-1$-dimensional sphere embedded in $\mathbb{R}^{2N}$. This allows one to
further reduce the optimization space as follows. Consider the sphere of
radius $\gamma$:

$\begin{split}\mathbb{S}_{\gamma}^{2N-1}=\\{\bm{q}_{\gamma}\in\mathbb{R}^{2N}\textrm{
such that }|\bm{q}_{\gamma}|=\gamma\\}.\end{split}$ (37)

The values of $\bm{q}_{\gamma}$ that are in the Lagrangian manifold must
satisfy

$\begin{split}\bm{A}_{0}(\theta_{0})+\hat{\gamma}\bm{Z}_{\bm{q},\theta_{0}}\bm{c}&=\bm{A}_{0}(\theta_{0})+\bm{q}_{\gamma},\\\
\end{split}$ (38)

where $\bm{A}_{0}(\theta_{0})$ consists of $\bm{x}$\- and $\bm{v}$-components
of $\tilde{\bm{A}}_{0}(\theta_{0})$ and the matrix,
$\bm{Z}_{\bm{q},\theta_{0}}$ is a sub-matrix of dimension $2N\times 2N$
extracted from $\bm{\Phi}^{u}_{\theta_{0}}$ such that

$\left[\begin{matrix}\bm{Z}_{\bm{q},\theta_{0}}\\\
\bm{Z}_{\bm{p},\theta_{0}}\end{matrix}\right]\equiv\bm{\Phi}^{u}_{\theta_{0}}.$
(39)

It follows from (38) that

$\hat{\gamma}\bm{c}=\bm{Z}_{\bm{q},\theta_{0}}^{-1}\bm{q}_{\gamma},$ (40)

which can be used to solve for the initial costate perturbations:

$\begin{split}\bm{p}_{\gamma,\theta_{0}}&=\bm{Z}_{\bm{p},\theta_{0}}\hat{\gamma}\bm{c}\\\
&=\bm{Z}_{\bm{p},\theta_{0}}\bm{Z}_{\bm{q},\theta_{0}}^{-1}\bm{q}_{\gamma}.\end{split}$
(41)

With the aid of (37) and (41), one can define the space of initial conditions
as

$W\equiv\left\\{\tilde{\bm{A}}_{0}(\theta_{0})+\left[\begin{matrix}\bm{q}_{\gamma}\\\
\bm{p}_{\gamma,\theta_{0}}\end{matrix}\right]\,\leavevmode\nobreak\
\vline\leavevmode\nobreak\
\,\bm{q}_{\gamma}\in\mathbb{S}^{2N-1}_{\gamma},\theta_{0}\in\mathbb{S}^{1}_{T}\right\\},$
(42)

which is now a $2N$-dimensional manifold.

In this work, the authors have used $\gamma=10^{-20}$ for $N=1$, a single
oscillator, and $\gamma=10^{-15}$ for the $N>1$ oscillators. Note that it is
not worth it to pick $\gamma$ too small as there is a trade-off between the
error in the approximation of the Lagrangian manifold and the error
accumulating due to numerical integration of Hamiltonian paths.

###  5.2 The Escape Time

The escape cost from the basin of $\hat{\bm{A}}_{0}$ for a Hamiltonian path
$(\bm{q},\bm{p})$ defined by (31) and the initial condition
$(\bm{q}_{0},\bm{p}_{0})$ of the form (42) is found by numerical integration
of the action (19) until the moment of time $\tau(\bm{q}_{\gamma},\theta_{0})$
when the path reaches the basin boundary $\partial B_{0}$:

$\tau(\bm{q}_{\gamma},\theta_{0})=\inf_{t\in\mathbf{R}}\\{\bm{q}(t)\in\partial
B_{0},\leavevmode\nobreak\ \|\bm{q}(\theta_{0})\|=\gamma\\}.$ (43)

The escape time for most Hamiltonian paths is finite when radius $\gamma$ in
(42) is positive, except for when the path is a heteroclinic orbit. In fact,
the results presented in this paper suggest that it is exactly the case that
the most probable escape paths approach attractors as $t\rightarrow-\infty$
and approach saddle cycles on basin boundaries as $t\rightarrow+\infty$.

###  5.3 Optimization and Implementation

Figure 4: An illustration for the origin of discontinuities in the cost
function (LABEL:AP_Cost) on the example of a single bistable Duffing
oscillator $\ddot{{x}}+\dot{{x}}-{{x}}+{x}^{3}=\sqrt{\epsilon}{\eta}(t)$. This
oscillator has stable equilibria at $(\pm 1,0)$ and a saddle at the origin.
The basin of $(-1,0)$ is shaded cyan. The optimal Hamiltonian escape path and
a collection of other Hamiltonian paths with close initial conditions are
shown in red and black, respectively. It is apparent that overshooting initial
conditions results in a much longer escape route and a sharp increase of
escape cost.

Initialization _Select an initial point $\bm{x}_{0}$ on a $d$-dimensional
manifold $\mathcal{M}$ and an initial step size $\sigma$. Set $i=0$. Select
the minimal step size $\sigma_{\min}$. _

Main Body

while  _$\sigma\geq\sigma_{\min}$ _ do

1: Sample a random unit vector $\bm{e}_{1}\in T_{\bm{x}_{i}}\mathcal{M}$ (the
tangent space of $\mathcal{M}$ at $\bm{x}_{i}$) from a uniform distribution.

2: Complete $\bm{e}_{1}$ to an orthonormal basis
$[\bm{e}_{1},\bm{e}_{2},...,\bm{e}_{d}]$ in the tangent space
$T_{\bm{x}_{i}}\mathcal{M}$.

3: Consider a collection of points
$Y=\\{\,\bm{y}_{j}=\bm{x}_{i}+\sigma\,\bm{e}_{j},\,j=[1,2,...,d],$

$\quad\quad\quad\bm{y}_{j}=\bm{x}_{i}-\sigma\,\bm{e}_{j},\quad
j=[d+1,d+2,...,2d]\\}$.

4: Calculate the cost $C(\bm{y}_{j})$ at each $\bm{y}_{j}$, $1\leq j\leq 2d$.

5: if  _$\min_{1\leq j\leq 2d}C(\bm{y}_{j})\geq C(\bm{x}_{i})$ _ then

Remain at $\bm{x}_{i}$ and reduce the step size:
$\bm{x}_{i+1}=\bm{x}_{i},\leavevmode\nobreak\ \sigma=\sigma/2.$

else

Set $\tilde{\bm{x}}=\arg\min_{1\leq j\leq 2d}C(\bm{y}_{j})$ and project
$\tilde{\bm{x}}$ onto the manifold $\mathcal{M}$: $\bm{x}_{i+1}={\sf
Proj}_{\mathcal{M}}(\tilde{\bm{x}})$.

end if

6: $i=i+1$.

end while

Algorithm 1 Stochastic Unit Vectors

Based on the argument presented in Sections 5.1 and 5.2 the authors reduce the
optimization problem (23) for the most probable escape path (MPEP) to

$\begin{split}&C(\bm{q}_{\gamma},\theta_{0})=\frac{1}{2}\int_{\theta_{0}}^{\tau(\bm{q}_{\gamma},\theta_{0})}||\bm{p}_{\bm{v}}(t)||^{2}dt\rightarrow\min\\\
\end{split}$ (44)

where $\bm{p_{v}}\equiv\bm{u}$ follows the dynamics of (31) with initial
conditions defined by $\eqref{AP_ICReduced}$. The integral time bounds are
defined by $\eqref{AP_ICReduced}$ and (43). The authors minimize
(LABEL:AP_Cost) with respect to both $\bm{q}_{\gamma}$ and $\theta_{0}$ in the
manifold $\mathcal{M}:=\mathbb{S}^{2N-1}_{\gamma}\times\mathbb{S}^{1}_{T}$.
The contribution of $\theta_{0}$ to the minimizers is discussed in Section 6.

A set of $n_{\rm init}$ initial points $\bm{q}_{\gamma}$ and $\theta_{0}$ is
sampled uniformly on the manifold
$\mathcal{M}:=\mathbb{S}^{2N-1}_{\gamma}\times\mathbb{S}^{1}_{T}$. The
corresponding initial conditions for the Hamiltonian paths are defined
according to (42). Next, each Hamiltonian path (31) together with its cost
functional (19) is integrated until a crossing of the basin boundary is
detected. The corresponding moment of time is
$\tau(\bm{q}_{\gamma},\theta_{0})$. The found values of the cost functional
are first sorted in the ascending order, and then $n_{s}<n_{\rm init}$ initial
points $\bm{q}_{\gamma}$ and $\theta_{0}$ corresponding to $n_{s}$ smallest
values of the cost functional are selected as the candidates for local
minimizers. The values of $n_{\rm init}$ and $n_{s}$ are found in Table
LABEL:tab:ComputationalCostTable in the appendix.

The cost function (LABEL:AP_Cost) is discontinuous and its minima lie at the
bottoms of its “cliffs” – see Figure 10 ahead depicting the escape cost for a
single oscillator. Such behavior of the cost function can be understood as
follows. Consider a sequence of Hamiltonian paths with initial conditions
$(\bm{q}_{\gamma}^{j},\theta_{0}^{j})$ escaping from the basin of the
attractor $\hat{\bm{A}}_{0}$ whose escape points approach a saddle cycle lying
at the boundary of $B_{0}$ and whose escape costs
$C(\bm{q}_{\gamma}^{j},\theta_{0}^{j})$ decrease. The goal of numerical
minimization of (LABEL:AP_Cost) is to generate a sequence of Hamiltonian paths
converging to a limiting path that has minimal escape cost and approaches the
saddle cycle infinitely long. Suppose the optimizer “overshoots”. Then the
path will fail to exit the basin of $\hat{\bm{A}}_{0}$ near the saddle and
proceed to an extra loop causing a jump in the escape cost. This phenomenon is
illustrated on a simple example of a bistable Duffing oscillator
$\ddot{{x}}+\dot{{x}}-{{x}}+{x}^{3}=\sqrt{\epsilon}{\eta}(t)$ in Figure 4.

Figure 5: A flow chart of the process of evaluating the cost function given an
input $\bm{q}_{\gamma}$ and $\theta_{0}$. A detailed version of the process
shown here is described in the pseudocode of Algorithm 2.

Input: Initial conditions $\bf{q}_{\gamma},\theta_{0}$

Initialization

i: Define $t_{i}\leftarrow\theta_{0}$, $i\leftarrow 0$

ii: Define initial condition $(\bm{q}(t_{i}),\bm{p}(t_{i}))\in W$ using
$\bf{q}_{\gamma},\theta_{0}$

iii: Let index0 be the index of the initial attractor.

iv: Define a function, ${\tt Aindex}(\bm{q}(t))$, that returns $m$ if
$\|{\bm{q}(t)-\bm{A}_{m}(t)}\|\leq R\quad\forall m=[0,1,...,n_{a}-1]$ or
noindex otherwise. Here, $R$ is a small radius and $\bm{A}_{m}(t)$ is the
position of attractor $m$ at time $t$.

v: Define the system’s differential equations as $ODE_{S}$

vi: Define the system’s corresponding Hamiltonian differential equations as
$ODE_{H}$

Escape via Integration

while  __ ${\sf index}==0$__ do

1:  $i\leftarrow i+1,\,t_{i}\leftarrow t_{i-1}+dT$

2: Evaluate $(\bm{q}(t_{i}),\bm{p}(t_{i}))$ by integrating $ODE_{H}$ forward
in time from $t_{i-1}$ to $t_{i}$ with initial condition
$(\bm{q}(t_{i-1}),\bm{p}(t_{i-1}))$.

3: $j\leftarrow 0,\,\check{t}_{j}\leftarrow
t_{i},\,\check{\bm{q}}(\check{t}_{j})\leftarrow\bm{q}(t_{i})$

4: while  __${\tt Aindex}(\check{\bm{q}}(\check{t}_{j}))=={\tt noindex}$__ do

a:  $j\leftarrow j+1,\,\check{t}_{j}\leftarrow\check{t}_{j-1}+dT$

b: Evaluate $\check{\bm{q}}(\check{t}_{j})$ by integrating $ODE_{S}$ forward
in time from $\check{t}_{j-1}$ to $\check{t}_{j}$ with initial condition
$\check{\bm{\bm{q}}}(\check{t}_{j-1})$.

end while

5: ${\tt index}\leftarrow{\tt Aindex}(\check{\bm{q}}(\check{t}_{j}))$

end while

Refine the Escape Time

i: Define upper and lower time bounds for the escape time:
$\tau_{H}=t_{i},\tau_{L}=t_{i-1}$

ii: $d\tau\leftarrow\tau_{H}-\tau_{L}$

while  __$d\tau >\epsilon$__ do

1: $\tau\leftarrow(\tau_{L}+\tau_{H})/2$

2: $j\leftarrow
0,\,\check{t}_{j}\leftarrow\tau,\,\check{\bm{q}}(\check{t}_{j})\leftarrow\bm{q}(\tau)$

3: while  _${\tt Aindex}(\check{\bm{q}}(\check{t}_{j}))=={\tt noindex}$ _ do

a:  $j\leftarrow j+1,\,\check{t}_{j}\leftarrow\check{t}_{j-1}+dT$

b: Evaluate $\check{\bm{q}}(\check{t}_{j})$ by integrating $ODE_{S}$ forward
in time from $\check{t}_{j-1}$ to $\check{t}_{j}$ with initial condition
$\check{\bm{q}}(\check{t}_{j-1})$.

end while

4: if  __${\tt Aindex}(\check{\bm{q}}(\check{t}_{j}))=={\tt index0}$_ _ then

$\tau_{L}\leftarrow\tau$

else

$\tau_{H}\leftarrow\tau$

end if

5: $d\tau\leftarrow\tau_{H}-\tau_{L}$

end while

Evaluate the Cost

Evaluate the cost by integrating (LABEL:AP_Cost) along the path
$(\bm{q}(t),\bm{p}(t))$ on the time interval $\theta_{0}\leq t\leq\tau$.

Algorithm 2 The Cost Function: Step 4 in Algorithm 1

A variety of gradient-based and gradient-free methods and their combinations
have been explored as optimizers for (LABEL:AP_Cost). The gradient of the cost
function has been approximated by finite differences.

The gradient-based methods include the gradient descent and the Fletcher-
Reeves nonlinear conjugate gradient methods projected on the manifold
$\mathcal{M}=\mathbb{S}^{2N-1}_{\gamma}\times\mathbb{S}^{1}_{T}$. The
gradient-free techniques that the authors have used are the Nelder-Mead
simplex method and an originally developed _stochastic unit vectors_ method
(SUV) summarized in the pseudocode in Algorithm 1. While gradient-based
methods are more efficient while stepping within smooth regions of the cost
function than Nelder-Mead or SUV, they are unable to crawl along the edges of
the bottoms of the cliffs. Hence, as a proximity of a cliff is detected, a
gradient-free method continues the search. The authors have found that there
is no significant advantage in starting minimization with a gradient-based
technique. Furthermore, the brute-force SUV algorithm turned out to work the
most robustly and reliably on the problems addressed in this work.

As with all gradient free methods that cannot guarantee convergence to a
global minimum, it is best to validate the minimizer with a secondary
approach. The authors increase the probability of finding a global minimum by
first sampling the cost function space with a large number of initial points
prior to beginning the optimization as mentioned earlier in this section. The
authors validate minimizers in this work by checking that they approach saddle
cycles, and that the resulting quasipotential values follow a trend with
respect to parameters such as frequency and coupling. Since the optimization
and the initial sampling are stochastic, the authors also validate the
minimizers by running the optimization multiple times with different randomly
picked initial conditions.

The other challenge of numerical minimization of the cost function
(LABEL:AP_Cost) for two or more oscillators, i.e., when the dimension $d$ of
the manifold $\mathcal{M}$ is four or higher, is that the basins of attraction
must be sampled online during optimization as it is impossible to store high-
resolution data of the basins of attraction or the basin boundaries for $d\geq
4$. Only the attractors are determined offline and stored. Therefore, Step 4
of Algorithm 1; that is, the calculation of the exit cost that involves the
integration of Hamiltonian paths and the detection of crossing basin boundary
requires elaboration. Each Hamiltonian path $(\bm{q}(t),\bm{p}(t))$ traced by
the minimization algorithm is checked at equispaced moments of time $t_{i}$,
$i=1,2,\ldots$ for having left or not the basin of the initial attractor. This
is done by integrating the original unperturbed ODE (10) starting from the
point $\bm{q}(t_{i})=(\bm{x}(t_{i}),\dot{\bm{x}}(t_{i}))$ at time $t_{i}$
forward in time. If it is found that that $(\bm{q}(t_{i}),t_{i}{\rm\,mod\,}T)$
does not belong to the basin of the initial attractor, tracing of the
Hamiltonian path is stopped and a refinement procedure is implemented to
determine the escape time $\tau(\bm{q}_{\gamma},\theta_{0})\in(t_{i-1},t_{i})$
more precisely. The evaluation of exit cost is detailed in Algorithm 2 and
illustrated in Fig. 5. An important parameter in this Algorithm is $dT$, the
time interval between any two consecutive check times $t_{i}$ and $t_{i+1}$.
If $dT$ is too small, the cost function can take a long time to evaluate. If
$dT$ is too large, the unstable dynamics of the Hamiltonian system blow up
during the time span, which increases the integration time and can cause the
algorithm to miss the initial escape from the basin. It was found that a good
default value for $dT$ is half-period, $T/2$. However, when the basins are
shallow or small, a smaller value of $dT$ should be used. The code is set up
to decrease $dT$ and restart automatically if a blow-up of a Hamiltonian path
is detected. The code is written in MATLAB. The ODE solver ode45 with absolute
and relative error tolerances of $10^{-8}$ is used for numerical integration.

The location of the attractors and the frequency response diagrams in sections
6,7, and 8 were computed by using a combination of the shooting method and a
continuation package. The shooting method was used to generate a stable
solution of the system at one excitation frequency. The continuation package
AUTO2007 [7] was used to continue the solution branch for different frequency
values. The use of the continuation procedure allows for determining the
stable solutions branches, the unstable solution branches, the Floquet
multipliers of each solution, and the bifurcation points with respect to the
excitation frequency. This additional information is useful for visualization
and understanding of the behavior of the system, but not necessary to solve
the optimization problem for the most probable escape path.

##  6 Forced Duffing Oscillator

Figure 6: Frequency response of a forced Duffing oscillator as a function of
excitation frequency, $\omega$. Attractor and saddle locations are highlighted
at the frequency of $\omega=1.4$. This figure was generated by using the
continuation package AUTO2007 [7].

Figure 7: Poincaré Map of the Basins of Attraction at the frequency of
$\omega=1.4$.

This section is devoted to a detailed analysis of noise-driven transitions
between the high amplitude and low amplitude attractors of a single monostable
Duffing oscillator with external harmonic forcing. The methodology introduced
in the previous section is applied to find the MPEPs. The case of a single
oscillator is instructive and visual, as its study illuminates the
difficulties of the problem and explains the authors’ technical decisions.

The dynamics of the Duffing oscillator studied in this section is governed by
SDE (1) with $N=1$. The values of the parameters chosen are
$\alpha=1,\,\beta=0.3,\,\delta=0.1,$ and $\,F(t)=0.4\text{ cos}(\omega t)$.
The frequency response of this oscillator is shown in Fig. 6. As the frequency
$\omega$ of the external force increases, the dynamics of the Duffing
oscillator undergoes bifurcations at $\omega_{1}\approx 1.27$ and at
$\omega_{2}\approx 1.57$. At $0<\omega<\omega_{1}$, at least one periodic
solution exists and its $L_{2}$-norm is an increasing function of $\omega$. At
$\omega=\omega_{1}$, a pair of a stable and unstable periodic solutions are
born. As a result, there exist two stable periodic solutions for
$\omega_{1}<\omega<\omega_{2}$, one of which with the smaller $L_{2}$-norm is
born at $\omega=\omega_{1}$, while the other with higher $L_{2}$-norm one
continues to exist from lower values of $\omega$. They will be referred to as
the low amplitude attractor and the high amplitude attractor, respectively. At
$\omega=\omega_{2}$, the high amplitude periodic solution annihilates with the
unstable periodic solution born at $\omega=\omega_{1}$. Only the low amplitude
attractor exists for $\omega>\omega_{2}$, and its $L_{2}$-norm decreases as
$\omega$ increases.

For the frequency of $\omega=1.4$, the Poincaré map of the periodic solutions
and the basins of attractions of each attracting solution at $\theta=0$ is
shown in Fig. 7.

The autonomous Hamiltonian system corresponding to this oscillator is given by

$\left(\begin{array}[]{c}\dot{q}_{1}\\\ \dot{q}_{2}\\\ \dot{p}_{1}\\\
\dot{p}_{2}\\\ \dot{\theta}\end{array}\right)=\left(\begin{array}[]{c}q_{2}\\\
-\beta\,{q_{1}}^{3}-\alpha\,q_{1}+p_{2}-\delta_{c}\,q_{2}+F\cos(\omega\theta)\\\
p_{2}\,\left(3\,\beta\,{q_{1}}^{2}+\alpha\right)\\\ \delta_{c}\,p_{2}-p_{1}\\\
1\end{array}\right),$ (45)

where
$(q_{1},q_{2},p_{1},p_{2},\theta)\in\mathbb{R}^{4}\times\mathbb{S}^{1}_{T}$.
The action functional for this system is:

$\begin{split}S(\varphi,\theta)&=\frac{1}{2}\int_{\theta_{0}}^{\tau}\,p_{2}^{2}\,dt\\\
&=\frac{1}{2}\int_{\theta_{0}}^{\tau}\,\left(\beta\,\varphi^{3}+\alpha\,\varphi+\mathrm{\ddot{\varphi}}+\delta_{c}\,\mathrm{\dot{\varphi}}-F\cos(\omega\theta)\right)^{2}dt.\end{split}$
(46)

The Lagrangian manifold for this systems is three-dimensional, and the set $W$
of initial conditions for the Hamiltonian paths is, respectively, two-
dimensional:
$(\bm{q}_{\gamma},\theta_{0})\in\mathbb{S}^{1}_{\gamma}\times\mathbb{S}^{1}_{T}$.

Figure 8: The cost (46) to escape from the high amplitude attractor as a
function of $\zeta$ and the initial time $\theta_{0}=t_{0}\thinspace{\rm
mod}\thinspace T$. The variable $\zeta$ is the angle around a circle that
completely surrounds the initial attractor at time $\theta_{0}$. Three lowest
local minima at $\theta_{0}=0$ are marked with red stars. They correspond to
the minimum escape cost up to numerical errors.

Figure 9: The cost to escape from the high amplitude attractor as a function
of $\zeta$ for fixed phase $\theta_{0}=0$. The variable $\zeta$ is the angle
around a circle that completely surrounds the initial attractor at
$\theta_{0}$. Note that the minimizers of the cost function arise immediately
next to discontinuities.

Figure 10: The escape action from the high amplitude attractor is minimized
with respect to $\zeta$ at each initial phase $\theta_{0}$. Notice that with
some small error tolerance, the minimum action does not change with
$\theta_{0}$.

The escape cost from the high amplitude attractor to the low amplitude
attractor is displayed in Fig. 8 as a function of
$(\zeta,\theta_{0})\in\mathbb{S}^{1}_{\gamma}\times\mathbb{S}^{1}_{T}$. The
variable $\zeta$ is the angle around the circle $\mathbb{S}^{1}_{\gamma}$.
Fig. 8 is used to elucidate the complexity of the optimization problem.
Discontinuities of the cost function form “cliffs”. The cost function has
multiple local minima located at the bottom of the cliffs. Furthermore, the
cost function along the curves delineating the bottom of the cliffs is nearly
constant, as one can see in Figs. Fig. 8, 9, and 10 – its variation along
these curves is comparable to the numerical errors. Therefore, there are whole
families of escape paths with values of cost functional matching up to
numerical errors. Hence, MPEPs are non-unique at least up to the available
numerical accuracy.

The Poincaré section of one of the MPEPs corresponding to the global minimum
of the cost functional (46), specifically, the one corresponding to the
initial condition marked with the bottom red star in Fig. 8, is shown in Fig.
11. Here and throughout the rest of the paper, the Poincaré sections all
correspond to the phase angle $\theta=0$. The solid marks are the
intersections of the MPEP with the plane $\theta=0$. The MPEP goes through one
revolution between each pair of consecutive points of its Poincaré section.
The line segments connecting the marks are added only to indicate the ordering
of these intersections – they are _not_ the projections of the MPEP segments
onto the plane $\theta=0$.

The Poincaré section of the MPEPs defined by the initial conditions
corresponding to the three red stars in Fig. 8 are shown in Fig. 12. These
MPEPs escape from the high amplitude attractor and approach the saddle cycle
lying in the basin boundary of the high and low amplitude attractors. The
Poincaré section of a MPEP going from the low amplitude attractor to the
saddle cycle is shown in Fig. 13. Figures 12 and 13 suggest that _all MPEPs_
have infinitely long escape from the attractors and an infinitely long
approach to the saddle cycles lying in the basin boundaries.

Figure 11: The Poincaré section of one MPEP that escapes the high amplitude
attractor. This escape path begins at the attractor, rotates clockwise in the
Poincaré section, and ends at the saddle.

Figure 12: The Poincaré sections of three MPEPs, which escape the high
amplitude attractor with $\theta_{0}=0$, are plotted on top of each other.
Notice that these escape paths all overlap in phase space, but cross the
Poincaré section at different times. These MPEPs and all other MPEPs of this
system share the same quasipotential, share this same characteristic in phase
space, and have the same probability of occurrence.

Figure 13: The Poincaré section of a MPEP that escapes the low amplitude
attractor. This escape path begins at the attractor, rotates counter-clockwise
in the Poincaré section, and ends at the saddle.

Figure 14: Euclidean norm,
$|\bm{q}(t)|=\sqrt{\sum_{i}^{2N}q_{i}\left(t\right)^{2}}$, at all time steps
of a MPEP that escapes the high amplitude attractor (solid) and of an
unperturbed trajectory that descends from the saddle to the low amplitude
attractor (dotted).

Figure 15: Euclidean norm of a MPEP that escapes the low amplitude attractor
(solid) and of an unperturbed trajectory that ascends from the saddle to the
high amplitude attractor (dotted).

For comparison, the Poincaré section of a trajectory of (10) going from the
saddle cycle to the high amplitude attractor is visualized in Fig. 16. Note
that both Poincaré sections of the MPEPs from the high amplitude attractor to
the saddle cycle and of the trajectory going the other way around are
clockwise.

Figure 16: The Poincaré section of a trajectory of equation (3) for $N=1$ that
is drawn to the saddle cycle as $t\rightarrow-\infty$, is rotated clockwise in
the Poincaré section, and is drawn to the high amplitude attractor as
$t\rightarrow\infty$. The rotation of the trajectory in phase space matches
the rotation of the escape path of that same basin; this is common in lower
dimensional systems such as the bistable unforced Duffing oscillator shown in
Fig. 3, and is noteworthy because one could naively expect the most probable
escape paths to move along the space of a trajectory but in the opposite
direction, and this is not the case.

Another way to visualize a MPEP is to plot the Euclidean norm of its
$\mathbf{q}$-components versus time. The advantage of this technique is that
it allows for visualization of escape paths in a system consisting of an
arbitrary number of oscillators. The Euclidean norm of the MPEP (the solid
line) corresponding to the bottom star in Fig. 8 is displayed in Fig. 14. This
MPEP is concatenated with a heteroclinic trajectory (the dotted line) going
from the saddle cycle to the low amplitude attractor. A similar visualization
of a MPEP going from the low- to the high amplitude attractor is presented in
Fig. 15.

Figure 17: (top) $L_{2}$ norm of the frequency response of the forced Duffing
oscillator with attractor branches highlighted in color in the hysteresis
region. (bottom) The quasipotential as a function of excitation frequency;
action values corresponding to escapes from the low and high amplitude
attractors are shown in blue and red, respectively.

The escape costs from the low- to high amplitude attractor and vice versa;
that is, the _quasipotential barriers_ , are plotted as functions of the
excitation frequency $\omega$ in Fig. 17. The quasipotential barrier to escape
the low amplitude attractor approaches zero near the bifurcation point of the
deterministic dynamics at $\omega\approx 1.27$. This agrees with experimental
and numerical observations in literature [11, 44], wherein noise is more
likely to induce a change in stable mode from the low- to the high amplitude
attractor near the jump up point of the hysteresis region. Similarly, the
quasipotential barrier for the escape from the high-to low amplitude attractor
decays to zero as the frequency approaches the second bifurcation point
$\omega=\omega_{2}\approx 1.57$ at which the high amplitude attractor
annihilates. Therefore, the quasipotential barrier serves as a quantifier of
the escape difficulty from the basin of the corresponding attractor. At
$\omega\approx 1.3825$, the quasipotential barriers for both escapes are the
same.

##  7 Two Coupled Oscillators

The case of two coupled Duffing oscillators is interesting because this is the
simplest example that illuminates the effect of coupling and sheds the light
on transition mechanisms in larger oscillator arrays. Its dimensionality is
high enough to require online detection of basin crossing by Hamiltonian
paths, as outlined in Algorithm 2. The values of the parameters
$\alpha=1,\,\beta=0.3,\,\delta=0.1,$ and $\,F(t)=0.4\text{ cos}(\omega t)$ are
the same as for the single oscillator in the previous section. The coupling
coefficient is set to be $\nu=0.01$.

Figure 18: Frequency response of two coupled forced Duffing oscillators as a
function of excitation frequency, $\omega$. Attractor and saddle locations are
shown at the frequency $\omega=1.4$. This figure was generated by using the
continuation package AUTO2007 [7].

With these parameter values, the deterministic system exhibits the frequency
response depicted in Fig. 18. In the hysteresis region
$\omega_{1}<\omega<\omega_{2}$, where, $\omega_{1}\approx 1.27$ and
$\omega_{2}\approx 1.57$, there are $9$ solution branches. Solid and dashed
curves represent stable and unstable periodic solutions, respectively. As
described in Section 1, letters $H$ and $L$ denote states of individual
oscillators moving along a high amplitude orbit and a low amplitude orbit,
respectively. The letter $S$ encodes an intermediate state in which an
oscillator is close to the saddle orbit of the single oscillator studied in
Section 6. The highest (lowest) branch corresponds to the stable periodic
solution in which both oscillators vibrate at high (low) amplitudes, the $HH$
and $LL$ attractors, respectively. Two overlapping stable branches (the $HL$
and $LH$ attractors) in the middle represent two symmetric stable periodic
solutions where the first oscillator (oscillator 1) has high amplitude while
the second oscillator (oscillator 2) has low amplitude and vice versa. The
remaining branches correspond to unstable periodic solutions of the saddle
type.

Figure 19: Hyperplanes of a Poincaré map of a MAP that escapes the $HH$
attractor.

Figure 20: Euclidean norm of a MAP that escapes the $HH$ attractor (solid) and
of an unperturbed trajectory that from the $SH$ saddle to the $LH$ attractor
(dotted).

For $N=2$ coupled oscillators, the optimization problem for the initial
condition for the optimal Hamiltonian path is four-dimensional. That is, the
optimal initial condition $(\bm{q}_{\gamma}^{\ast},\theta_{0}^{\ast})$ is
sought on $\mathbb{S}^{3}_{\gamma}\times\mathbb{S}_{T}^{1}$, the direct
product of the three-dimensional sphere of radius $\gamma$ embedded into
$\mathbb{R}^{4}$ and the circle of circumference $T$. This optimization
problem has been solved to find optimal escape paths from the attractors $HH$,
$LL$, and $HL$. Due to symmetry, the solution for the escape from the $LH$
attractor is found for the one for the $HL$ attractor by switching the
oscillators. A collection of optimal paths visualized in Figs. 19–24
corresponds to the excitation frequency $\omega=1.4$.

The escape path from the $HH$ attractor is displayed in Fig. 19 by plotting
the Poincaré sections of its projections onto $(x_{1},v_{1})$\- and
$(x_{2},v_{2})$-planes. After escaping from the basin of the $HH$ attractor
along this path, the system arrives at the $SH$ saddle. From there, it can
descend either to the $LH$ attractor or back to the $HH$ attractor. Fig. 19
shows that the first oscillator follows a path similar to the escape path from
the high- to low amplitude attractor of a single oscillator in Fig. 12, while
the motion of the second oscillator experiences a relatively small change
throughout this transition. Note the closeness of the projections of the $HH$
and $LH$ attractors onto the $(x_{2},v_{2})$-plane. The $L_{2}$-norm of the
$\bm{q}(t)$-component of the transition path from the $HH$ to the $LH$
attractor via the $SH$ saddle is plotted in Fig. 20.

Figure 21: Hyperplanes of a Poincaré map of a MPEP that escapes the $LL$
attractor.

Figure 22: Euclidean norm of a MPEP that escapes the $LL$ attractor (solid)
and of an unperturbed trajectory that ascends from the $LS$ saddle to the $LH$
attractor (dotted).

The escape path from the $LL$ attractor to the $LS$ saddle lying on the
boundary separating the $LL$ and $LH$ basins is shown in Figs. 21 and 22. The
path from $LL$ to $SL$ can be found using symmetry.

The escape path from the $HL$ attractor shown in Figs. 23 and 24 leads to the
$SL$ saddle. It has been found that the quasipotential of the $SL$ saddle with
respect to the $HL$ attractor is lower than that of the $HS$ saddle. This
means that the escape from $HL$ via $HS$ is exponentially less likely than the
escape via $SL$. Furthermore, no direct transition between the localized modes
$HL$ and $LH$ has been found. This suggests that a transition between $HL$ and
$LH$ can only be accomplished in the weak noise limit by visiting the $LL$ or
the $HH$ mode in-between.

Figure 23: Hyperplanes of a Poincaré map of a MPEP that escapes the $HL$
attractor.

Figure 24: Euclidean norm of a MPEP that escapes the $HL$ attractor (solid)
and of an unperturbed trajectory that descends from the $SL$ saddle to the
$LL$ attractor (dotted).

The MPEPs and the corresponding quasipotential barriers have been computed for
a set of values of the excitation frequency $\omega$ spanning the hysteresis
region $[\omega_{1},\omega_{2}]$. The quasipotential barriers between each
attractor and the corresponding lowest saddle are plotted as functions of
$\omega$ in Fig. 25. It is instructive to compare this figure with the similar
figure for the single oscillator (Fig. 17). It is evident that in the case of
weak coupling considered here, the transitions to and from the localized modes
in the two-oscillator system are well approximated by the transitions in the
single oscillator. The values of the quasipotential barriers for the escapes
from the $LL$ , $HL$ , and $HH$ attractors for $\omega=1.4$ and the coupling
coefficient $\nu=0.01$ are given in Table 1.

Furthermore, the following observation can be made:

  * •

Escape from the $LL$ mode requires smaller minimum action near the jump up
frequency and larger minimum action near the jump down frequency.

  * •

Escape from the $HH$ mode requires larger minimum action near the jump up
frequency and goes to zero near the jump down frequency.

  * •

The $\Lambda$-shape of the orange curve in Fig. 25 is caused by the fact that
there are two possible escapes from each localized mode: to $HH$ and to $LL$.
For $\omega<1.3825$, the required minimum action for the escape from $HL$ is
close but slightly less than the quasipotential barrier for the escape from
$LL$. For $\omega>1.3825$, the required minimum action for the escape from
$HL$ is close but slightly less than the quasipotential barrier for the escape
from $HH$. This reduction in quasipotential is due to the coupling: the
optimal escape paths for $N=2$ are sought in a larger space than for $N=1$.
Hence the quasipotential barriers for $N=2$ can only decrease in comparison to
those for the single oscillator.

Figure 25: (top) $L_{2}$ norm of the frequency response of the $N=2$ system
with attractor branches highlighted in color in the hysteresis region.
(bottom) The quasipotential as a function of excitation frequency; action
values corresponding to escapes from the $LL$, $HL$, $HH$ attractors are shown
in blue, yellow, and red respectively.

##  8 Circular Oscillator Arrays

Figure 26: Frequency response of a circular array of $N=3$ coupled forced
Duffing oscillators as a function of excitation frequency, $\omega$. Unique
attractors are shown at the frequency $\omega=1.4$. This figure was generated
by using the continuation package AUTO2007 [7].

Figure 27: Euclidean norm of a MPEP that escapes the $LHL$ attractor (solid)
and of an unperturbed trajectory that descends from the $LSL$ saddle to the
$LLL$ attractor (dotted).

Figure 28: Euclidean norm of a MPEP that escapes the $LHLLL$ attractor (solid)
and of an unperturbed trajectory that descends from the $LSLLL$ saddle to the
$LLLLL$ attractor (dotted).

The application to circular three- and five-oscillator arrays presented in
this section demonstrates the viability of the proposed methodology for high-
dimensional systems and sheds light of noise-induced transitions in larger
circular arrays.

The systems under consideration are governed by SDE (1) with $N=3$ and $N=5$,
i.e., with three and five oscillators, respectively. They are models of
vibrations of blades of turbomachinery or cyclically arranged VEHs. The
parameter values are the same as for $N=2$. The frequency responses of $N=3$
and $N=5$ systems are shown, respectively, in Figs. 26 and 2. Provided that
the coupling is weak, the number of stable vibrating modes grows dramatically
with $N$. For $N=3$, this number is at least 8, while for $N=5$, it is at
least 32, as mentioned in the introduction. Two of the stable modes correspond
to all oscillators vibrating at low or high amplitudes, while the rest are
localized modes in which some oscillators have high amplitudes while others
have low ones. Each localized mode has a number of equivalent localized modes
obtained by circular shift and/or changing direction of rotation of oscillator
indexing. Moreover, nonequivalent localized modes with equal numbers of high
amplitude oscillators have very close time-averaged $L_{2}$-norms of their
$\bm{q}$-components (11) resulting in indistinguishable points in the
frequency-response diagram. The authors have chosen to focus here on the study
of escapes from localized modes for $N=3$ and $N=5$ in which one oscillator
has high amplitude while the other $N-1$ oscillators have low amplitudes. The
escapes from these modes are the most relevant to noise influenced
turbomachinery.

The most probable escape paths from the basins of the $LHL$ and $LHLLL$
attractors of the $N=3$ and $N=5$ systems at frequency $\omega=1.4$ and
coupling coefficient $\nu=0.01$ are shown in Figs. 27 and 28, respectively.
After these escapes, the systems land onto the attractors in which all
oscillators have low amplitudes. These escape paths are similar to the MPEP
leaving the localized mode of the $N=2$ system: the path of the oscillator
originally at high amplitude resembles the escape path of the single
oscillator, $(N=1)$, from its high attractor, while the motions of the other
oscillators are only slightly affected. The quasipotential barriers for these
escapes, $U_{LHL}$ and $U_{LHLLL}$ are slightly smaller than the
quasipotential barrier $U_{H}$ for the escape of the single oscillator from
its $H$ attractor basin. Furthermore, as expected, the following trend is
observed:

$U_{H}>U_{HL}>U_{LHL}\approx U_{LHLLL}.$

The values of these barriers are given in Table 1.

The dependence of the quasipotential barrier for the escape from the $LHL$
attractor on the excitation frequency $\omega$ in the hysteresis region for
three oscillators is very similar to the one for the escape from the $HL$
attractor for $N=2$ – compare Figs. 29 and 25.

Figure 29: (top) The $L_{2}$-norm (11) of the frequency response of the $N=3$
system with attractor branches highlighted in color in the hysteresis region.
(bottom) The quasipotential corresponding to escapes from the $LHL$ attractor
as a function of excitation frequency.

Figure 30: The quasipotential barriers as a function of coupling coefficient $\nu$ that correspond to the following: (a) escapes from the $H$ attractor of the single Duffing oscillator, (b) escapes from the localized mode of the single side coupling ($N=2$) system, and (c) and (d) escapes from the localized mode with a single oscillator being at high amplitude in circular arrays ($N=3$) and ($N=5$), respectively. Table 1: The quasipotential barriers for the escapes from localized modes at which only one oscillator has high amplitude for $N=1,2,3,5$ and from all modes for $N=2$ at the frequency of $\omega=1.4$ with $\nu=0.01$.  $N$ | Initial Attractor | Final Attractor | $S(\varphi)$
---|---|---|---
$1$ | $H$ | $L$ | $0.129$
$2$ | $HL$ | $LL$ | $0.125$
$3$ | $LHL$ | $LLL$ | $0.122$
$5$ | $LHLLL$ | $LLLLL$ | $0.122$
1-4 $1$ | $L$ | $H$ | $0.180$
$2$ | $LL$ | $HL$ | $0.179$
1-4 $2$ | $HH$ | $HL$ | $0.155$

The dependencies of the quasipotential barriers $U_{H}$, $U_{HL}$, $U_{LHL}$,
and $U_{LHLLL}$ on the coupling coefficient $\nu$ at $\omega=1.4$ are plotted
in Fig. 30. The barriers for the circular array are lower than those for the
single oscillator. In turn, the barriers for $N=3$ and $N=5$ are smaller than
for $N=2$. Moreover, the difference between these barriers grows with the
increase in the inter-oscillator coupling coefficient $\nu$. Note that the
basins of localized modes get shallower as $\nu$ increases, and as a
consequence, it is difficult to find most probable escape paths for
$\nu>0.07$. Furthermore, the localized modes become unstable as $\nu$ reaches
$0.08$; hence, localizations do not exist for strong inter-oscillator
couplings, which is also noted in prior studies [42].

Fig. 30 indicates a very important fact that one can predict the qualitative
behavior of most probable escape paths and the values of the quasipotential
barriers for the escapes from the localized mode where only one oscillator has
high amplitude by comparing the results obtained for $N=3$ and $N=5$. It is
evident that the differences between the quasipotential barriers $U_{LHL}$ and
$U_{LHLLL}$ at the feasible range of $\nu$ are small and comparable with
numerical errors. Therefore, for weak inter-oscillator coupling facilitating
the existence of response localizations, one can consider the barrier
$U_{LHL}$ as an adequate approximation for $U_{LHLLL}$. This suggests that the
circular nature of coupling effectively limits the influence of the array on a
given oscillator to the influence of its nearest neighbors. As a consequence,
the quasipotential barriers for the escapes from the localized modes
$LHL\ldots L$ for $N>5$ are close to those for $N=3$ and $N=5$. Meanwhile, the
numerical solution is cheaper and more accurate for $N=3$ than for larger $N$.

##  9 Discussion

The contribution of this work is two-fold.

  1. 1.

A methodology for computing most probable escape paths and quasipotential
barriers from periodic attractors of arrays of nonlinear oscillators with
external periodic forcing and small white noise has been developed.

  2. 2.

This methodology has been applied for a case study of a single monostable
harmonically driven Duffing oscillator and circular arrays of $N=2$, 3, and 5
oscillators.

The methodology relies on the large deviations theory [21] and Feynman’s path
integral theory [50], exploits Pontryagin’s maximum principle for the optimal
control problems[22] and the Floquet theory for linear ODEs with periodic
coefficients [25], and uses the Action Plot Method [23] as a building block.
To the best of the authors’ knowledge, this is the first work in which the
Action Plot Method was promoted to high dimensions. The main challenges in
this promotion are:

  1. 1.

the arising optimization problem for the initial condition for the optimal
Hamiltonian path has a discontinuous objective function with multiple local
minima located at the bottoms of the “cliffs”;

  2. 2.

the final time for a Hamiltonian path is the time when it reaches the
appropriate basin boundary, and basin boundaries need to be identified during
the optimizations for higher dimensional systems.

The first challenge has been tackled by the introduction of the Stochastic
Unit Vectors (SUV) method, while the second one has been overcome via routine
checks in which basin the state of the system is in and refinement of the
escape time and location once a basin switch is detected. These techniques
have shown themselves to be very robust and suitable for high dimensions.

The maximal number of coupled oscillators $N=5$ considered in this work is not
the limit for the proposed methodology. However, as the number of oscillators
$N$ increases, the runtime of the optimization routine grows as $O(N^{2})$
(see Appendix A). Furthermore, the behavior of the green curve in Fig. 30
corresponding to $N=5$ suggests that the found MPEP and the corresponding
quasipotential barrier become less accurate. Therefore, in order to apply the
proposed methodology to larger oscillator arrays it is desirable to make the
optimization method more efficient.

The case study of circular oscillator arrays has shown an important trend: the
predictions for the most probable escape paths and quasipotential barriers for
response localizations in large arrays can be made based on the results for
small arrays. In particular, it has been shown that the results for $N=5$ are
well-approximated by those for $N=3$. This is consistent with the results
obtained by means of reduced-order modeling in Ref. [42].

The methodology described here has been generalized for different stiffness
models. In this case study, the authors chose
$K(\bm{x})=\alpha\bm{x}+\beta\bm{x}^{3}+\nu D_{N}\bm{x}$, but $K(\bm{x})$
could contain higher order nonlinear stiffness terms, higher order coupling
terms, or inhomogeneous stiffness models, for example. In particular,
quantifying quasipotentials for escapes from response localizations subjected
to noise in systems that have both cubic stiffness and cubic coupling could be
a direction of future work. It is known that such systems exhibit two types of
response localizations known as s-t mode and p-mode, and that the two types of
localized modes change stability depending on the cubic stiffness coefficient
and the cubic coupling coefficient. Quantifying the quasipotential for those
systems would reveal how robust the two types of localized modes are to noise.

The methodology developed in this work has an important limitation: the noise
term in SDE (1) is the standard Brownian motion. On each fixed interval
$\Delta t$, the increment of the Brownian motion is $N$-dimensional Gaussian
random variable with mean 0 and the covariance matrix $I_{N\times N}$, which
is unbounded. In physical experiments, only bounded amplitude and bounded
bandwidth noise can be realized. Bounded amplitude can be modeled by
truncating and renormalizing increments of Brownian motion. This will result
in rendering all Hamiltonian paths for which $\|\bm{p}_{\bm{v}}\|$ exceeds the
maximum allowed noise magnitude infeasible, i.e., the corresponding action
functional value is $+\infty$. Modifying the action functional for noise of
limited bandwidth is less straightforward. The extension for more realistic
noise models is left for the future.

##  10 Data Availability Statement

All data used in this work were generated within this work. The authors’ codes
are available upon request.

## Acknowledgments

The authors thank Professor Sandra Cerrai for a valuable discussion on the
large deviations theory. This work was partially supported by the Clark
Doctoral Fellowship (LC), NSF CAREER Grant DMS-1554907 (MC), AFOSR MURI grant
FA9550-20-1-0397 (MC), and NSF grant CMMI-1760366 and associated data science
supplements (BB). A part of this research project was conducted by using
computational resources available at the Maryland Advanced Research Computing
Center (MARCC).

##  11 Appendix A: Notes on Computational Cost and Nomenclature

Table 2: Parameters and computational cost Description |  Symbol |  Example 1 |  Example 2 |  Example 3
---|---|---|---|---
Number of oscillators |  $N$ |  1 |  2 |  5
Excitation Frequency |  $\omega$ |  1.4 |  1.4 |  1.4
Coupling Coefficient |  $\nu$ |  - |  0.01 |  0.01
Dimension of cost function |  $d$ |  2 |  4 |  10
Sphere Magnitude |  $\gamma$ |  $10^{-20}$ |  $10^{-15}$ |  $10^{-15}$
Basin check interval |  $dT$ |  $\pi/\omega$ |  $\pi/\omega$ |  $\pi/\omega$
Min. Escape Radius |  $r_{A_{0}}$ |  0.1 |  0.2 |  0.5
Min. Attracted Radius |  $r_{A_{k}}$ |  0.1 |  0.2 |  0.5
Number of initial samples of cost space |  $n_{init}$ |  300 |  1000 |  1000
Number of parallel optimizations |  $n_{s}$ |  10 |  15 |  15
Initial SUV step size |  $\sigma$ |  0.25 |  0.25 |  0.25
Maximum optimization iterations |  $i_{m}$ |  35 |  35 |  35
Min. number of optimizations to keep |  $n_{m}$ |  5 |  5 |  5
Number of available processors |  $c$ |  24 |  24 |  24
1-5 Average Processor Efficiency ($\%$) |  $e$ |  $54\%$ |  $69\%$ |  $81\%$
Real time cost (seconds) |  $t$ |  160 |  365 |  4942

This appendix includes an estimate of the computational cost of the
optimization algorithm for finding MPEPs and three examples of parameter
choices. The purpose of this appendix is to facilitate the reproduction of the
results of this work by a third party with reasonable expectations of
computational cost.

To save computational cost, a $2N-1$-dimensional sphere of radius $r_{A_{0}}$
embedded in $\mathbb{R}^{2N}$ that surrounds the initial attractor
$\bm{A}_{0}(\theta)$ at a fixed moment of time, $\theta$, is introduced. Here,
$\bm{A}_{0}(\theta)$ describes the $(\bm{x},\bm{v})$ coordinates of the
attractor at time $\theta$. The radius $r_{A_{0}}$ is chosen so that if the
path is still within the sphere, the path has not left the current basin yet.
Similarly, $r_{A_{k}}$, is the radius of a $2N-1$-dimensional sphere that
surrounds attractor $A_{k}$. When computing the domains of attraction online,
it is assumed that if a trajectory enters the sphere surrounding attractor,
$A_{k}$ than the trajectory is in the basin of attractor $k$.

Table LABEL:tab:ComputationalCostTable is a list of the definitions of many of
the parameters relevant to the oscillator arrays studied in this work and to
the proposed numerical algorithms. Omitted parameters are of lesser influence
on the computational cost. $t_{i}$ is the cost of the initial sample set of
the cost function space and $t_{o}$ is time cost of the parallel optimizations
using the SUV algorithm. The order of $t_{o}$ is loosely bounded by an
inequality that depends on the number of optimizations running because the
algorithm is allowed to throw out high cost, poor improvement descents while
preserving useful descents. This is a gradient free optimization scheme on a
cost function with many local minima and discontinuities; finding a global
minima is not guaranteed.

One could estimate the influence of the parameters on the computational costs
as

$\begin{split}O(t_{i})&\approx O\left(\frac{d\times n_{init}}{c\times e\times
dT}\right)\\\ O\left(\frac{2d^{2}\times n_{m}}{c\times e\times dT}\right)&\leq
O(t_{o})\leq O\left(\frac{2d^{2}\times n_{s}\times i_{m}}{c\times e\times
dT}\right)\\\ t&\approx t_{i}+t_{o}\end{split}$ (47)

The principle idea of (47) is that the computational cost of evaluating the
cost function scales with $d\times n/dT$ and that the cost of optimization
using SUV scales the number of cost function evaluations by a multiple of
$2d$. With efficient implementations at the excitation frequency of
$\omega=1.4$, finding a solution to the $N=5$ system was a 10 dimensional
optimization problem that took $80$ minutes to solve with 24 processors in
MATLAB. Additional examples are shown in Table
LABEL:tab:ComputationalCostTable. These results were produced on the Maryland
Advanced Research and Computing Center’s (MARCC) Intel Haswell dual socket,
12-core processors, 2.5GHz, 30MB cache, 128GB RAM node.

The important nomenclature in this work is summarized in Table
LABEL:tab:NotationTable. This nomenclature is also described in the work as it
is used.

Table 3: Nomenclature Symbol |  Type |  Description
---|---|---
$N$ |  $\mathbb{W}$ |  Number of oscillators
$t$ |  $\mathbb{R}$ |  Time
$\bm{x}$ |  $\mathbb{R}^{N}$ |  Beam tip displacements
$\bm{v}$ |  $\mathbb{R}^{N}$ |  Beam tip velocities
$\bm{f}(t)$ |  $\mathbb{R}^{N}$ |  Deterministic periodic external excitation
$\bm{\eta}(t)$ |  $\mathbb{R}^{N}$ |  Stochastic external excitation
$\delta_{s}$ |  $\mathbb{R}$ |  Damping coefficient
$\alpha$ |  $\mathbb{R}$ |  Linear stiffness coefficient
$\beta$ |  $\mathbb{R}$ |  Cubic stiffness coefficient
$\nu$ |  $\mathbb{R}$ |  Linear coupling coefficient
$\epsilon$ |  $\mathbb{R}$ |  Noise intensity/multiplier
$\omega$ |  $\mathbb{R}$ |  Periodic excitation frequency
$T$ |  $\mathbb{R}$ |  Periodic excitation period
$\bm{D}_{N}$ |  $\mathbb{R}^{N\times N}$ |  Linear Coupling Matrix
$\bm{q}$ |  $\mathbb{R}^{2N}$ |  Vector of displacement states and velocity states $(\bm{x},\bm{v})$
$\bm{p}$ |  $\mathbb{R}^{2N}$ |  Vector of costates $(\bm{p_{x}},\bm{p_{v}})$
$\bm{u}$ |  $\mathbb{R}^{N}$ |  Vector of deterministic control inputs
$||\bm{q}||_{L_{2}}$ |  $\mathbb{R}$ |  $L_{2}$-like norm
$\bm{b}(\bm{q})$ |  $\mathbb{R}^{2N}$ |  Drift vector
$\bm{Q}(\bm{q})$ |  $\mathbb{R}^{2N\times 2N}$ |  Diffusion matrix
${S}_{\tau}(\bm{\varphi})$ |  $\mathbb{R}$ |  Action functional
$\bm{\varphi}$ |  $\mathbb{R}^{2N}$ |  Candidate escape path
$\tau$ |  $\mathbb{R}$ |  Travel time, final time boundary of $\varphi$, escape time
$\mathcal{A}$ |  char |  Attractor
$B$ |  char |  Basin of attraction
$\partial B$ |  char |  Basin boundary
$D$ |  char |  Domain
$\partial D$ |  char |  Domain boundary
$\mathbb{S}_{T}$ |  $\mathbb{S}^{1}$ |  A circle with circumference $T$
$\bm{K}(\bm{x})$ |  $\mathbb{R}^{N}$ |  Generalized restoring force function of $\bm{x}$
$\theta$ |  $\mathbb{S}_{T}$ |  A periodic time coordinate
$\hat{\bm{q}}$ |  $\mathbb{R}^{2N}\times\mathbb{S}_{T}$ |  Augmented state vector for autonomous system, $(\bm{x},\bm{v},\theta)$
$\hat{\bm{p}}$ |  $\mathbb{R}^{2N+1}$ |  Augmented costate vector for autonomous system, $(\bm{p_{x}},\bm{p_{v}},p_{\theta})$
$n_{a}$ |  $\mathbb{W}$ |  Number of attractors
$k$ |  $\\{0,1,...,n_{a}\\}$ |  Index of $n_{a}$ attractors
$k=0$ |  |  Index of the initial attractor
$A_{k}$ |  char |  The attractor with index k
$B_{k}$ |  char |  Basin of attraction of $A_{k}$
$\partial B_{k}$ |  char |  Basin boundary of $A_{k}$
$\partial B_{kl}$ |  char |  Intersection of Basin boundary of $A_{k}$ and $A_{l}$
$\bm{A}_{k}(\theta)$ |  $\mathbb{R}^{2N}$ |  The location in coordinates $(\bm{x}(\theta),\bm{v}(\theta))$ of attractor with index $k$ at time $\theta$
$\bm{A}_{0}$ |  $\mathbb{R}^{2N}$ |  The location of the initial attractor
$\hat{\bm{A}}_{k}$ |  $\mathbb{R}^{2N}\times\mathbb{S}_{T}$ |  The location in coordinates $(\bm{x},\bm{v},\theta)$ of attractor with index $k$
$\bar{\bm{A}}_{k}$ |  $\mathbb{R}^{4N+1}\times\mathbb{S}_{T}$ |  Augmented $\bar{\bm{A}}_{k}$ in Hamiltonian coordinates, $(\bm{x},\bm{v},\theta,\bm{0},\bm{0},0)$
$\tilde{\bm{A}}_{k}(\theta)$ |  $\mathbb{R}^{4N}$ |  Augmented $\bm{A}_{k}(\theta)$ in Hamiltonian coords. w/o explicit $\theta$ and $p_{\theta}$, $(\bm{x}(\theta),\bm{v}(\theta),\bm{0},\bm{0})$
$J$ |  $\mathbb{R}^{4N\times 4N}$ |  Jacobian of Hamiltonian system
$\bm{J_{K}}$ |  $\mathbb{R}^{N\times N}$ |  Jacobian of vector field $\bm{K}(\bm{x})$
$\bm{z}$ |  $\mathbb{R}^{4N}$ |  Coordinates of linearized system
$\bm{\Psi}(t)$ |  $\mathbb{R}^{4N\times 4N}$ |  Fundamental solution matrix
$\bm{\Phi}(t)^{s}_{\theta_{0}}$ |  $\mathbb{R}^{4N\times 2N}$ |  Stable eigenvectors of Monodromy matrix
$\bm{\Phi}(t)^{u}_{\theta_{0}}$ |  $\mathbb{R}^{4N\times 2N}$ |  Untable eigenvectors of Monodromy matrix
$\gamma$ |  $\mathbb{R}$ |  Small radius of a sphere that surrounds the initial attractor at fixed time
$\mathbb{S}_{\gamma}^{2N-1}$ |  $\mathbb{S}^{2N-1}$ |  A 2N sphere with radius $\gamma$
$\bm{q}_{\gamma}$ |  $\mathbb{S}_{\gamma}^{2N-1}$ |  A value of $\bm{q}$ with euclidean norm equal to $\gamma$
$dT$ |  $\mathbb{R}$ |  The amount of time to numerically map the path forward in time prior to checking if the path has entered a new basin of attraction
$L$ |  char |  ”Low amplitude”
$H$ |  char |  ”High amplitude”

## References

## References

  * [1] G. Duffing, “Erzwungene schwingungen bei veränderlicher eigenfrequenz,” Vieweg u. Sohn, Braunschweig, vol. 7, 1918.
  * [2] A. Papangelo, F. Fontanela, A. Grolet, M. Ciavarella, and N. Hoffmann, “Multistability and localization in forced cyclic symmetric structures modelled by weakly-coupled Duffing oscillators,” Journal of Sound and Vibration, vol. 440, pp. 202–211, Feb. 2019.
  * [3] A. Grolet and F. Thouverez, “Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model,” Journal of Sound and Vibration, vol. 331, pp. 2911–2928, June 2012.
  * [4] A. Dick, B. Balachandran, and c. Mote, “Intrinsic localized modes in microresonator arrays and their relationship to nonlinear vibration modes,” Nonlinear Dynamics, vol. 54, pp. 13–29, Jan. 2008.
  * [5] T. Ikeda, Y. Harata, and K. Nishimura, “Intrinsic localized modes of harmonic oscillations in nonlinear oscillator arrays,” Journal of Computational and Nonlinear Dynamics, vol. 8, no. 4, 2013.  Publisher: American Society of Mechanical Engineers Digital Collection.
  * [6] B. Balachandran, E. Perkins, and T. Fitzgerald, “Response localization in micro-scale oscillator arrays: influence of cubic coupling nonlinearities,” International Journal of Dynamics and Control, vol. 3, no. 2, pp. 183–188, 2015.  Publisher: Springer.
  * [7] E. J. Doedel, T. F. Fairgrieve, B. Sandstede, A. R. Champneys, Y. A. Kuznetsov, and X. Wang, “AUTO-07P: Continuation and bifurcation software for ordinary differential equations,” tech. rep., 2007.
  * [8] J. Emad, A. F. Vakakis, and N. Miller, “Experimental nonlinear localization in a periodically forced repetitive system of coupled magnetoelastic beams,” Physica D: Nonlinear Phenomena, vol. 137, pp. 192–201, Mar. 2000.
  * [9] M. Kimura and T. Hikihara, “Coupled cantilever array with tunable on-site nonlinearity and observation of localized oscillations,” Physics Letters A, vol. 373, pp. 1257–1260, Mar. 2009.
  * [10] M. Kimura and T. Hikihara, “Experimental manipulation of intrinsic localized modes in macro-mechanical system,” Nonlinear Theory and Its Applications, IEICE, vol. 3, no. 2, pp. 233–245, 2012.
  * [11] E. Perkins, M. Kimura, T. Hikihara, and B. Balachandran, “Effects of noise on symmetric intrinsic localized modes,” Nonlinear Dynamics, vol. 85, no. 1, pp. 333–341, 2016.  ISBN: 0924-090X Publisher: Springer.
  * [12] B. Niedergesäß, A. Papangelo, A. Grolet, A. Vizzaccaro, F. Fontanela, L. Salles, A. Sievers, and N. Hoffmann, “Experimental observations of nonlinear vibration localization in a cyclic chain of weakly coupled nonlinear oscillators,” Journal of Sound and Vibration, vol. 497, p. 115952, 2021.  Publisher: Elsevier.
  * [13] J. A. Kenyon and J. H. Griffin, “Forced Response of Turbine Engine Bladed Disks and Sensitivity to Harmonic Mistuning,” Journal of Engineering for Gas Turbines and Power, vol. 125, pp. 113–120, Dec. 2002\.
  * [14] A. V. Srinivasan, “Flutter and Resonant Vibration Characteristics of Engine Blades,” Journal of Engineering for Gas Turbines and Power, vol. 119, pp. 742–775, Oct. 1997.
  * [15] M. P. Castanier and C. Pierre, “Modeling and Analysis of Mistuned Bladed Disk Vibration: Current Status and Emerging Directions,” Journal of Propulsion and Power, vol. 22, pp. 384–396, Mar. 2006.
  * [16] I. A. Sever, Experimental validation of turbomachinery blade vibration predictions.  PhD thesis, Imperial College London (University of London), 2004.
  * [17] P. Gardonio and M. Zilletti, “Vibration energy harvesting based on mode localization,” ISMA 2014)(Katholieke Universiteit Leuven, Belgium, 15–17 September 2014), 2014.
  * [18] E. Lefeuvre, A. Badel, C. Richard, L. Petit, and D. Guyomar, “A comparison between several vibration-powered piezoelectric generators for standalone systems,” Sensors and Actuators A: Physical, vol. 126, pp. 405–416, Feb. 2006.
  * [19] D. J. Ewins, “The effects of detuning upon the forced vibrations of bladed disks,” Journal of Sound and Vibration, vol. 9, no. 1, pp. 65–79, 1969\.  Publisher: Elsevier.
  * [20] Y. Yan, P. Cui, and H. Hao, “Vibration mechanism of a mistuned bladed-disk,” Journal of sound and vibration, vol. 317, no. 1-2, pp. 294–307, 2008.  Publisher: Elsevier.
  * [21] M. I. Freidlin and A. D. Wentzell, “Random perturbations,” in Random perturbations of dynamical systems, pp. 15–43, Springer, 1998.
  * [22] M. Bardi and I. Capuzzo-Dolcetta, Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations.  Springer Science & Business Media, 2008.
  * [23] S. Beri, R. Mannella, D. G. Luchinsky, A. Silchenko, and P. V. McClintock, “Solution of the boundary value problem for optimal escape in continuous stochastic systems and maps,” Physical Review E, vol. 72, no. 3, p. 036131, 2005.  Publisher: APS.
  * [24] L. Lin, H. Yu, and X. Zhou, “Quasi-potential calculation and minimum action method for limit cycle,” Journal of Nonlinear Science, vol. 29, no. 3, pp. 961–991, 2019.  Publisher: Springer.
  * [25] C. Chicone, Ordinary Differential Equations with Applications.  Springer-Verlag, New York, 1999.
  * [26] Z. Chen, Y. Li, and X. Liu, “Noise induced escape from a nonhyperbolic chaotic attractor of a periodically driven nonlinear oscillator,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 26, no. 6, p. 063112, 2016\.  Publisher: AIP Publishing LLC.
  * [27] U. von Wagner and W. V. Wedig, “On the Calculation of Stationary Solutions of Multi-Dimensional Fokker–Planck Equations by Orthogonal Functions,” Nonlinear Dynamics, vol. 21, pp. 289–306, Mar. 2000.
  * [28] U. Von Wagner, “On double crater-like probability density functions of a duffing oscillator subjected to harmonic and stochastic excitation,” Nonlinear Dynamics, vol. 28, no. 3-4, pp. 343–355, 2002.
  * [29] W. Martens, U. von Wagner, and V. Mehrmann, “Calculation of high-dimensional probability density functions of stochastically excited nonlinear mechanical systems,” Nonlinear Dynamics, vol. 67, pp. 2089–2099, Feb. 2012.
  * [30] A. Förster, L. Panning-von Scheidt, and J. Wallaschek, “Approximate Solution of the Fokker–Planck Equation for a Multidegree of Freedom Frictionally Damped Bladed Disk Under Random Excitation,” Journal of Engineering for Gas Turbines and Power, vol. 141, Sept. 2018.
  * [31] M. I. Dykman, P. V. E. McClintock, V. N. Smelyanski, N. D. Stein, and N. G. Stocks, “Optimal paths and the prehistory problem for large fluctuations in noise-driven systems,” Phys. Rev. Lett., vol. 68, pp. 2718–2721, May 1992\.  Publisher: American Physical Society.
  * [32] D. G. Luchinsky, S. Beri, R. Mannella, P. V. McClintock, and I. Khovanov, “Optimal fluctuations and the control of chaos,” International Journal of Bifurcation and Chaos, vol. 12, no. 03, pp. 583–604, 2002.  Publisher: World Scientific.
  * [33] S. Ramakrishnan and B. Balachandran, “Energy localization and white noise-induced enhancement of response in a micro-scale oscillator array,” Nonlinear dynamics, vol. 62, no. 1, pp. 1–16, 2010.  Publisher: Springer.
  * [34] E. Perkins and B. Balachandran, “Noise-enhanced response of nonlinear oscillators,” Procedia Iutam, vol. 5, pp. 59–68, 2012.  Publisher: Elsevier.
  * [35] E. Perkins, C. Chabalko, and B. Balachandran, “Noise-influenced transient energy localization in an oscillator array,” Nonlinear Theory and Its Applications, IEICE, vol. 4, no. 3, pp. 232–243, 2013.  ISBN: 2185-4106 Publisher: The Institute of Electronics, Information and Communication Engineers.
  * [36] W. Xu, Q. He, T. Fang, and H. Rong, “Global analysis of stochastic bifurcation in Duffing system,” International Journal of Bifurcation and Chaos, vol. 13, no. 10, pp. 3115–3123, 2003.  Publisher: World Scientific.
  * [37] W. Xu, Q. He, T. Fang, and H. Rong, “Stochastic bifurcation in Duffing system subject to harmonic excitation and in presence of random noise,” International Journal of Non-Linear Mechanics, vol. 39, no. 9, pp. 1473–1479, 2004.  ISBN: 0020-7462 Publisher: Elsevier.
  * [38] J. S. Yu and Y. K. Lin, “Numerical path integration of a non-homogeneous Markov process,” International Journal of Non-Linear Mechanics, vol. 39, no. 9, pp. 1493–1500, 2004.
  * [39] P. Kumar and S. Narayanan, “Modified path integral solution of Fokker–Planck equation: response and bifurcation of nonlinear systems,” Journal of computational and nonlinear dynamics, vol. 5, no. 1, 2010.
  * [40] S. Narayanan and P. Kumar, “Numerical solutions of Fokker–Planck equation of nonlinear systems subjected to random and harmonic excitations,” Probabilistic Engineering Mechanics, vol. 27, no. 1, pp. 35–46, 2012.  ISBN: 0266-8920 Publisher: Elsevier.
  * [41] L. Cilenti and B. Balachandran, “Transient probability in basins of noise influenced responses of mono and coupled Duffing oscillators,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 31, p. 063117, June 2021\.  Publisher: American Institute of Physics.
  * [42] B. Balachandran, T. Breunung, G. D. Acar, A. Alofi, and J. A. Yorke, “Dynamics of circular oscillator arrays subjected to noise,” Nonlinear Dynamics, pp. 1–14, 2022.  Publisher: Springer.
  * [43] V. Agarwal, X. Zheng, and B. Balachandran, “Influence of noise on frequency responses of softening Duffing oscillators,” Physics Letters A, vol. 382, no. 46, pp. 3355–3364, 2018.  ISBN: 0375-9601 Publisher: Elsevier.
  * [44] A. M. Alofi, Coupled Oscillator Arrays: Dynamics and Influence of Noise.  PhD Thesis, 2021.
  * [45] Y. Chao and M. Tao, “Parametric resonance for enhancing the rate of metastable transition.”.
  * [46] A. Dembo and O. Zeitouni, Large Deviations Techniques and Applications.  Springer, New York, 1998.
  * [47] Z. Chen and M. Freidlin, “Smoluchowski–Kramers approximation and exit problems,” Stochastics and Dynamics, vol. 5, no. 04, pp. 569–585, 2005\.  Publisher: World Scientific.
  * [48] M. Freidlin, “Some remarks on the smoluchowski–kramers approximation,” Journal of Statistical Physics, vol. 117, no. 3/4, pp. 617–634, 2004.
  * [49] M. I. Dykman and M. A. Krivoglaz, “Theory of fluctuational transitions between stable states of a nonlinear oscillator,” Sov. Phys. JETP, vol. 50, no. 1, pp. 30–37, 1979.
  * [50] R. P. Feynman and A. R. Hibbs, Quantum mechanics and path integrals.  McGraw-Hill, New York, 1965.
  * [51] D. E. Kirk, Optimal control theory: an introduction.  Courier Corporation, 2004.
  * [52] V. I. Arnold, Mathematical Methods of Classical Mechanics.  Springer-Verlag, New York, 1978.