# Stability of a Parametrically Driven, Coupled Oscillator System: An
Auxillary Function Method Approach

Andrew N. McMillan
Department of Mathematics
University of Michigan, Ann Arbor
Ann Arbor, MI 48108
<EMAIL_ADDRESS>
& Yin Lu Young
Department of Mechanical Engineering and
Department of Naval Architecture and Marine Engineering
University of Michigan, Ann Arbor
Ann Arbor, MI 48108
<EMAIL_ADDRESS>
& Mary Robinson
Department of Physics
University of Michigan, Ann Arbor
Ann Arbor, MI 48108
<EMAIL_ADDRESS>

###### Abstract

Coupled, nonlinear oscillators are often studied in applied biology, physics,
fluids, and many other disciplines. In this paper, we study a parametrically
driven, coupled oscillator system where the individual oscillators are
subjected to varying frequency and phase with a focus on the influence of the
damping and coupling parameters away from parametric resonance frequencies. In
particular, we study the key long-term statistics of the oscillator system’s
trajectories and stability. We present a novel, robust and computationally
efficient method come to be known as an auxillary function method for long-
time averages, and we pair this method with classical, perturbative-asymptotic
analysis to corroborate the results of this auxillary function method. These
paired methods are then used to compute the regions of stability for a coupled
oscillator system. The objective is to explore the influence of higher order,
coupling effects on the stability boundary across a broad range of modulation
frequencies, including frequencies away from parametric resonances. We show
that both simplified and more general asymptotic methods can be dangerously
un-conservative in predicting the true regions of stability due to high order
effects caused by coupling parameters. The differences between the true
stability boundary and the approximate stability boundary can occur at
physically relevant parameter values in regions away from parametric
resonance. The differences between the solutions depends on the specific
parameters of the system, as explained in the results section. As an
alternative to asymptotic methods, we show that the auxillary function method
for long-time averages is an efficient and robust means of computing true
regions of stability across all possible initial conditions.

##  1 Introduction

The study of coupled oscillators is ubiquitous in both applied and theoretical
fields for their wide range of modeling capabilities [35, 18, 23], as well as
their rich mathematical structure [16, 17]. In applied fields, qualitative
analysis of an oscillator system’s trajectories or long term statistics in the
presence of varying parameters or initial conditions is of primary importance
[19, 36]. The stability of coupled oscillators are also of great interest in
many fields including engineering, quantum optics, bio-chemistry, etc [37, 38,
39]. Moreover, in theoretical work, the presence of parametric terms prohibits
classical, linear analysis, and novel analysis tools are an interest of
development. In this paper, we present a novel, theoretical method to study
the long term statistics of a parametrically driven, coupled oscillator
system.

To motivate this study, we provide several, key applications of coupled
oscillators. In the biological sciences, coupled oscillators have extensively
been employed in the study of circadian rhythms where individual circadian
oscillations are generated in suprachiasmatic nuclei by an internal regulatory
network [20, 1]. A general interest is to establish both when and how these
individual oscillations become synchronized or out of synchrony [21], and key
quantities, such as measures of synchrony, are realized as time averages of
the underlying oscillator’s correlations [1, 21]. Hence, computing parameter
and initial condition dependent time averages is of pivotal importance and
frequently challenging. Coupled oscillators have also been explored in the
neuroscience and neural network communities, respectively. In neuroscience,
for example, the Kuromoto model has received widespread attention. Past and
on-going work focuses on the study of long time average frequencies between
individual oscillators and investigates how these individual oscillators
become globally phase locked [22, 24]. Furthermore, in the study of neural
networks, data-driven models attempt to construct potentially nonlinear
dynamical models via data sampling [34, 40, 41]. This requires a potentially
large data set for the model to be trained on and one’s choice of data and
assumptions about the form of the model have far reaching consequences on the
stability of the underlying dynamics and the accuracy of model predictions.
Additionally, coupled oscillators have received widespread attention in the
work on fluids. For instance, extensive work has been done on the phenomena of
phase lock-ins for vortex induced vibrations [27]. The general interest is to
establish which range of frequencies induce unstable and growing amplitudes
for the underlying dynamics. Finally, there are a variety of applications in
physics. Networks of optical parametric oscillators are being investigated as
a way to simulate classical Ising (i.e. spin-$\frac{1}{2}$) systems. Such
simulators have exciting potential applications in quantum computation [6]. In
addition, the behavior of coupled parametric oscillators have been used to
study the response of microelectromechanical systems (MEMs) [33]. These
examples provide only a small glimpse of the wide range of applications that
these coupled oscillator systems have across a variety of scientific
disciplines.

Frequently, a primary interest is in studying the effects of the oscillators
coupling strength, the type of coupling, the phase differences between the
oscillators, and the magnitude, frequency and phase of the driving force on
the long term statistics of the system’s solutions. Although most real world
applications involve oscillators with differing phases or differing
frequencies, most of the analysis that appears in the literature has ignored
both of these essential features while also ignoring the effects of coupling
strength, especially at frequencies away from parametric resonance [6].
Ignoring frequencies away from parametric resonance frequencies is natural if
one is interested in studying phase-locking or if one has energy harvesting
applications in mind. However, real-world systems typically operate over a
broad range of parametric, modulation frequencies. Moreover, many engineering
and biological applications including coupled oscillator systems are designed
to operate away from parametric resonance to avoid excessive vibration, noise,
and accelerate fatigue. We demonstrate and explain in Section IV the
importance of considering higher order effects across parameter space for the
stability boundaries and long-time statistics particularly if the interest is
in regions away from parametric resonance frequencies.

The reason these essential features have primarily been ignored is almost
exclusively due to the problem’s difficulty and the limitations of
perturbative asymptotic methods [26]. From a computational perspective, direct
numerical simulation (DNS) is incredibly expensive. Moreover, one is limited
to the study of trajectories that are numerically stable, which pivotally
depends on one’s choice of initial conditions for nonlinear oscillators.
Frequently, one is interested in unstable solutions in particular for the
“control of chaos” [25]. These concerns naturally pose the problems of which
initial conditions to choose and a robust method to determine whether one’s
findings are generically sensitive to initial data. From a theoretical
perspective, analysis is notoriously difficult due to the presence of
parametric, nonlinear terms, particularly for coupled systems. Research is
frequently concerned with knowing the boundaries of parametric instabilities,
especially when the parametric term has a frequency that is different than a
multiple of the system’s natural frequency. Therefore, approximate solutions
are frequently derived via perturbative, asymptotic methods [28, 29]. However,
these asymptotic methods both fail to capture the full range of the parameters
in question and the effects of the oscillator’s coupling term, which can make
the higher order terms non-negligible. As high order terms are frequently
neglected, asymptotics foregoes when or how these higher order terms become
relevant [30].

In this paper, we demonstrate that our proposed method, come to be known as an
auxillary function method for long-time averages, is incredibly robust in
addressing the difficulties previously explicated. On computational grounds,
our method is computationally efficient and foregoes the need to manually
check various initial conditions. In fact, our method determines extremal time
averages across all initial conditions. On theoretical grounds, our method is
robust even in the face of highly non-linear dynamics and does not require
approximate expansions or restrictions to only small parameter values.
Therefore, we can recover the stability boundaries and long term statistics
without needing to ignore the higher order, non-linear terms that are
disregarded in perturbative asymptotic analysis. Hence, the use of the
auxillary function method for stability analysis of nonlinear coupled
oscillators is the primary and novel contribution of our work.

Our objective is to study a parametrically driven, coupled oscillator system
with the auxillary function method for long-time averages. We demonstrate the
robustness of our auxillary function method in computing essential long time
averages of the system’s trajectories and determine stability boundaries for
varying system parameters. We will also compare our results against the
results predicted by perturbative asymptotic analysis and investigate the
effects of higher order, coupling terms on the system’s stability and long-
term statistics across a broad range of physically relevant modulation
frequencies at and away from parametric resonance.

##  2 Auxillary Function Method for Long-time Averages

To introduce the auxiliary function method for long-time averages, we focus on
determining upper bounds for time averages of functions of the dynamical
variable for autonomous ordinary differential equations (ODEs).

Consider $\textbf{x}(t)\in\mathbb{R}^{d}$ satisfying

$\dot{\textbf{x}}=\textbf{f}(\textbf{x})$ (1)

for $\textbf{f}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ a continuously
differentiable vector field. When there is no confusion, we denote
$\textbf{x}_{i}(t)$ and $\textbf{f}_{i}(\textbf{x})$ as the i-th component of
the vectors x and $\textbf{f}(\textbf{x})$, respectively.

For the purposes of the computation of long term statistics, it becomes
natural to ask which trajectories of Equation (1) realize the optimal time
average and their corresponding averaged values. For a general quantity of
interest $\Phi(\textbf{x})$, we define its long-time average along the
trajectories $\textbf{x}(t)$ with $\textbf{x}(0)=\textbf{x}_{0}$ by:

$\overline{\Phi}(\textbf{x}_{0}):=\limsup_{T\rightarrow+\infty}\frac{1}{T}\int_{0}^{T}\Phi(\textbf{x}(t))dt.$
(2)

As the underlying dynamics are of physical relevance, it is natural to assume
that all trajectories $\textbf{x}(t)$ lie in a compact subset B of the phase
space $\mathbb{R}^{d}$. We are interested in the maximal long-time average
among all trajectories (eventually) remaining in B, i.e.,

$\overline{\Phi}^{*}=\max_{\textbf{x}_{0}\in\text{B}}\overline{\Phi}(\textbf{x}_{0}).$
(3)

Upper bounds on averages can be deduced using the fact that time derivatives
of bounded functions average to zero. This elementary observation follows from
the fact that for every $\text{V}(\textbf{x})\in C^{1}(\text{B})$—the set of
continuously differentiable functions on B—we have

$\begin{split}0&=\limsup\limits_{T\rightarrow+\infty}\frac{\text{V}(\textbf{x}(T))-\text{V}(\textbf{x}(0))}{T}\\\
&=\overline{\frac{d}{dt}\text{V}(\textbf{x}(\cdot))}\\\
&=\overline{\textbf{f}(\textbf{x}(\cdot))\cdot\nabla\text{V}(\textbf{x}(\cdot))}.\end{split}$
(4)

We hereafter refer to any such $\text{V}(\textbf{x})\in C^{1}(\text{B})$ as
“auxiliary” functions. Note that Equation (4) holds for any auxiliary function
so there is an infinite family of functions with the same time average as
$\Phi(\textbf{x})$. In particular,

$\overline{\Phi}(\textbf{x}_{0})=\overline{\Phi(\textbf{x}_{0})+\textbf{f}(\textbf{x}_{0})\cdot\nabla\text{V}(\textbf{x}_{0})}.$
(5)

For any auxiliary function one obtains an upper-bound on
$\overline{{\Phi}}^{*}$ by bounding the right hand-side point-wise on B and
subsequently maximizing the left hand side over initial data $\textbf{x}_{0}$:

$\overline{{\Phi}}^{*}\leq\max_{\textbf{x}\in\text{B}}[\Phi(\textbf{x})+\textbf{f}(\textbf{x})\cdot\nabla\text{V}(\textbf{x})].$
(6)

The sharpest upper bound on $\overline{{\Phi}}^{*}$ is then

$\overline{{\Phi}}^{*}\leq\inf_{\text{V}\in
C^{1}(\text{B})}\max_{\textbf{x}\in\text{B}}\
[\Phi(\textbf{x})+\textbf{f}(\textbf{x})\cdot\nabla\text{V}(\textbf{x})],$ (7)

and we note that the minimization on the right hand side is convex in the
auxillary function $V$. The remarkable fact is that the inequality in Equation
(7) can actually be made an equality. The details of the proof can be found in
[12], but the key ideas are that maximizing time averages can be realized as
maximizing phase space averages against invariant measures via Birkhoff’s
Ergodic Theorem; maximizing phase space averages over invariant probability
measures can be realized as a Lagrange multiplier problem, where the auxillary
function $V$ can be interpreted as a Lagrange multiplier; finally, the
resulting max-min problem can be written as a min-max problem using a variety
of abstract min-max theorems over infinite dimensional spaces.

Note that the bound $\overline{\Phi}^{*}\leq U$ follows from Equation (6) if

$U-\Phi(\textbf{x})-\textbf{f}(\textbf{x})\cdot\nabla V\geq
0\,\forall\textbf{x}\in\text{B},$ (8)

and if one restricts to polynomial dynamics, then the problem reduces to a
convex optimization problem subject to a local, non-negative polynomial
constraint.

When B is a semi-algebraic set 111A subset B in $\mathbb{R}^{d}$ is called
semi-algebraic if it is defined by a finite collection of polynomial equations
of the form $P(\textbf{x})=0$ and $Q(\textbf{x})>0$ for
$\textbf{x}\in\mathbb{R}^{d}$., one can augment the constraint in Equation (8)
with polynomials that define B; this is called the S-procedure. The non-
negative polynomial constraint in Equation (8) can then be replaced with a sum
of squares222For $\textbf{x}\in\mathbb{R}^{d}$, a polynomial $p(\textbf{x})$
is said to be a sum of squares if $p(\textbf{x})=\Sigma p_{i}(\textbf{x})^{2}$
for some finite collection of polynomials. The set of all sum of squares
polynomials over $\textbf{x}\in\mathbb{R}^{d}$ will be denoted
$\mathcal{S}[\textbf{x}]$. constraint, which has an equivalent semi-definite
program formulation. In the resulting semi-definite program, the polynomial
degree of the auxillary function $V$ is a degree of freedom, where increasing
the degree of $V$ improves the sharpness of the computed bounds. When
increasing the degree of $V$ no longer improves the bound, we say the bounds
are sharp, at least to numerical precision; see [12] for discussions and
reviews on sum of squares programs via semi-definite programming.

One should note that the above formulation is completely dependent on the
original dynamical system in Equation (1) being autonomous. For our purposes,
the model in question will have non-linear, non-autonomous trigonometric
dependence. Therefore, the dynamics generically take on the form:

$\dot{\textbf{x}}=\textbf{f}(\textbf{x},\cos(\omega t),\sin(\omega t)).$ (9)

The canonical way of making systems of the form seen in Equation (9)
autonomous is to introduce an additional coordinate $x_{d+1}=t$ and extend the
system dimension from $d$ to $d+1$. However, this method unfortunately
introduces an unbounded dependent variable while retaining non-polynomial
dependence on it. This problem can be circumvented by introducing two new
dynamical variables satisfying the polynomial sub-system

$\begin{split}\dot{x}_{d+1}=&-\omega x_{d+2}\\\ \dot{x}_{d+2}=&\omega
x_{d+1}\\\ \text{Subject to:}\quad&x^{2}_{d+1}+x^{2}_{d+2}=1.\end{split}$ (10)

One should note that equations of the form in Equation (10) are frequently
referred to as differential-algebraic systems- that is, a differential
equation subject to an algebraic constraint. In order to write the system in
Equation (10) as a standard dynamical system, the form seen in Equation (9),
one needs to omit the origin as a potential solution either via assumption or
other techniques.

Additionally, one can also formulate equivalent autonomous polynomial dynamics
for both quasiperiodic and substantially more complex
$\frac{2\pi}{\omega}$-periodic time dependences in the vector field. Employing
a new pair of dynamical variables like those in Equation (10) for each
independent frequency allows for quasiperiodic time dependence, at least for
quasiperiodicity involving only a finite number of independent frequencies.
Other $\frac{2\pi}{\omega}$-periodic time functions can be expressed as finite
linear combinations of $\cos(n\omega t)$ and $\sin(n\omega t)$, each of which
in turn is a finite polynomial combination of $\cos(\omega t)$ and
$\sin(\omega t)$. The overall order of the dynamical system necessarily
increases but autonomous polynomial dynamics are still sufficient to capture
the systems’ dynamics; see [13] for further details.

There are a few requisite remarks that need to be made in view of this
auxillary function method. The first is that this auxillary function method
computes the maximal long time average of trajectories across all initial
conditions, but it does not give the knowledge of which initial conditions or
trajectories attain the maximal long time average. If one were interested in
the analysis of the particular solution that achieves the maximal long-time
average, other methods would need to be developed. This is perhaps one
potential drawback of the method. However, in many applications, the knowledge
of these extremal averages is sufficient for the study at hand. This method is
particularly well-suited for systems subject to random initial conditions and
when long-time statistics are of the utmost importance. Finally, we note for
the interested reader that there are extensions of this auxillary function
method to partial differential equations applications; see the work of [31, 2,
3] for explication.

##  3 Model Setup

Our general model of interest is a parametrically driven, coupled oscillator
system of the form

$\begin{split}&\ddot{x}_{1}+{\omega_{0}}^{2}[1+h\sin(\gamma
t+\frac{\phi}{2})]x_{1}+\omega_{0}g\dot{x}_{1}-\omega_{0}r\dot{x}_{2}=0\\\
&\ddot{x}_{2}+{\omega_{0}}^{2}[1+h\sin(\gamma
t-\frac{\phi}{2})]x_{2}+\omega_{0}g\dot{x}_{2}+\omega_{0}r\dot{x}_{1}=0.\end{split}$
(11)

In Equation (11), $\omega_{0}$ denotes the proper frequency of the oscillators
and $g$ is the intrinsic loss term, which is taken to be equal for both
oscillators for the sake of simplicity. The $h$ and $\gamma$ terms are the
intensity and frequency of the parametric stiffness terms, respectively. We
focus on $h\in[0,1]$ and $\gamma\in[-3,3]$, as this is the regime which
encapsulates most physically relevant phenomena. The coupling strength $r$
describes an energy preserving coupling between the oscillators, which
corresponds to rotations in the $(x_{1},x_{2})$ plane and preserves the
system’s total energy. We note that in the limit as $r\rightarrow 0$ that the
system in Equation (11) becomes equivalent to two, decoupled parametric
oscillators described by two, damped Mathieu equations. Finally, $\phi$
denotes the phase difference between the oscillators. Just as in [6], we will
focus on the cases $\phi=0,\frac{\pi}{2},\text{and}\,\pi$, and we remark that
the system in Equation (11) for the aforementioned $\phi$ values exhibits
three potential resonance frequencies at $\gamma=2\Omega_{r}$ and
$\gamma=2\Omega_{r}\pm\omega_{0}r$, where
$\Omega_{r}=\omega_{0}\sqrt{1+(r^{2}-g^{2})/4}$.

For the purposes of the auxillary function method for long-time averages, the
system in Equation (11) can be written as a first order, coupled system with
exclusively polynomial terms:

$\begin{split}\dot{x}_{1}&=y_{1}\\\ \dot{x}_{2}&=y_{2}\\\
\dot{y}_{1}&=-\omega_{0}^{2}[1+hx_{3}]x_{1}-\omega_{0}gy_{1}+\omega_{0}ry_{2}\\\
\dot{y}_{2}&=-\omega_{0}^{2}[1+hx_{3}]x_{2}-\omega_{0}gy_{2}-\omega_{0}ry_{1}\\\
\dot{x_{3}}&=\gamma x_{4}\\\ \dot{x}_{4}&=-\gamma x_{3}\\\ &\text{Subject to:
$x_{3}^{2}+x_{4}^{2}=1$},\end{split}$ (12)

where $\phi=0$. The system in Equation (12) can be written similarly for
additional values of $\phi$ via elementary trigonometric identities. The
constraint $x_{3}^{2}+x_{4}^{2}=1$ is enforced via the S-procedure and forces
the variables $x_{3}$ and $x_{4}$ to be uniquely determined. However, we
remark that the S-procedure is not enforced when $\gamma=\phi=0$, as Equation
(12) becomes independent of trigonometric terms. We employ the auxillary
function method for long-time averages on the quantity of interest

$\Phi=x_{1}^{2}+x_{2}^{2}.$ (13)

Hence, computing $\overline{\Phi}^{*}$ gives the maximal long-time average of
the summed, mean squared amplitudes, which allows us to distinguish between
small and large amplitude, parameter dependent solutions. The intuition being
that small or large amplitude solutions should generically correlate with
asymptotically stable or unstable solutions, respectively. Since the maximal
long-time averages will be either zero or infinite in the regions of stability
and instability, respectively, the auxillary function method will allow us to
trace out the threshold between stability and instability. The fineness of the
threshold, moreover, only depends on one’s mesh size for the system parameters
333The mesh size is 80 by 80 for every figure in this paper.. With $\Phi$ as
defined in Equation (12), the sum of squares program to solve becomes:

$\begin{gathered}\min\,\,U\\\ \text{s.t.}\quad
U-\Phi-\textbf{f}(x_{1},x_{2},y_{1},y_{2},x_{3},x_{4})\cdot\nabla V+\\\
S(x_{1},x_{2},y_{1},y_{2},x_{3},x_{4})(1-x_{3}^{2}-x_{4}^{2})\in\mathcal{S}[x_{1},x_{2},y_{1},y_{2},x_{3},x_{4}]\\\
S(x_{1},x_{2},y_{1},y_{2},x_{3},x_{4})(1-x_{3}^{2}-x_{4}^{2})\in\mathcal{S}[x_{1},x_{2},y_{1},y_{2},x_{3},x_{4}],\end{gathered}$
(14)

where $\textbf{f}(x_{1},x_{2},y_{1},y_{2},x_{3},x_{4})$ is the polynomial
vector field as defined in Equation (11), $V$ is an auxillary function
polynomial whose degree may vary, and we have enforced the S-procedure at all
$\gamma$ values other than zero. The convex optimization problems in this
paper were solved using Mosek[4] paired with Yalmip[5].

We remark that similar computations using the auxillary function method have
already appeared in [13], where the model of interest was the Duffing
oscillator. In particular, it was shown that the auxillary function method
accurately reproduces the Duffing equation’s frequency response curve and
parameter dependent hysteresis phenomena as derived by harmonic balance down
to numerical precision. This further validates the auxillary function method
as a tool to study oscillator dynamics.

To the author’s knowledge, this is the first time that the auxillary function
method has been applied to studying nonlinear, coupled oscillators or has been
used to explicitly study dynamical stability. Moreover, as we suspect and as
we will show, the dynamics of Equation (11) are substantially more complicated
than that which can be described by asymptotic methods. The reason is that in
contrast to the Duffing equation, our model has two coupled, nonlinear
oscillators and the non-autonomous forcing appears parametrically instead of
externally. Hence, the synchronization of the two oscillators plays a crucial
role in their stability and that synchronization pivotally depends on one’s
choice of system parameters.

In Section V. and in a similar fashion to [13], the parameter dependent
computations of $\overline{\Phi}^{*}$ will be used to recover stability region
boundaries and the validity of these results will be corroborated via
asymptotic analysis.

##  4 Asymptotic Analysis

We will corroborate the findings of the auxillary function method for long-
time averages via comparison with the approximate solutions established by
Floquet theory. This is a standard method which is discussed by many sources
on asymptotic analysis; the specific reference used here is [32].

Consider $x(t)\in\mathbb{R}$ satisfying

$\ddot{x}(t)+F(t)x(t)=0,$ (15)

where $F(t)$ is a periodic function with period $T$. Equation (15) is a second
order differential equation and thus has two linearly independent solutions
$x_{1}(t)$ and $x_{2}(t)$.

Since $x_{1}(t)$ and $x_{2}(t)$ are linearly independent, they span the
solution space of Equation (15) In particular, there must exist
$\textbf{M}\in\mathbb{R}^{2\times 2}$ such that

$\begin{bmatrix}x_{1}(t+T)\\\
x_{2}(t+T)\end{bmatrix}=\textbf{M}\begin{bmatrix}x_{1}(t)\\\
x_{2}(t)\end{bmatrix},$ (16)

where M is called the _monodromy matrix_ of $x_{1}(t)$ and $x_{2}(t)$. A
monodromy matrix can be constructed for any pair of linearly independent
solutions of Equation (15). One can show that the determinant of the monodromy
matrix is $1$ and eigenvalues of the monodromy matrix are independent of the
pair of solutions that one chooses to start with.

Suppose $\lambda_{1}$, $\lambda_{2}$ are eigenvalues of a monodromy matrix M
corresponding to solutions $x_{1}(t)$, $x_{2}(t)$. For each $\lambda_{j}$
there exists a linear combination $x_{\lambda_{j}}(t)$ of $x_{1}(t)$ and
$x_{2}(t)$ such that

$x_{\lambda_{j}}(t+T)=\lambda_{j}x_{\lambda_{j}}(t).$ (17)

_Floquet’s theorem_ states that we may write

$x_{\lambda_{j}}=e^{-i\mu_{j}t}u_{\lambda_{j}}(t),$ (18)

where $u_{\lambda_{j}}(t)$ is periodic with period $T$ and
$\mu_{j}\in\mathbb{C}$ is such that $\lambda_{j}=e^{-i\mu_{j}T}$. Since
$u_{\lambda_{j}}(t)$ is periodic, one may expand $u_{\lambda_{j}}(t)$ as a
Fourier series and recover the analytic solution $x_{\lambda_{j}}(t+T)$ to
within arbitrary accuracy by recursively solving for the Fourier amplitudes.

For Equation (11) with $\phi=0$, we generally follow the method outlined in
[6], where the study focused on frequencies at or near parametric resonance
due to an interest in phase-locking, but our novel contribution is that we do
not limit our attention to resonant frequencies.

The equations in Equation (11) can be decoupled by performing a change of
basis and defining

$x_{\pm}(t)=x_{1}\pm ix_{2},$ (19)

Adding the first line of Equation (11) to $\pm i$ times the second line yields

$\ddot{x}_{\pm}+\omega_{0}^{2}\left(1+h\sin(\gamma
t)\right)x_{\pm}+\omega_{0}\left(g\pm ir\right)\dot{x}_{\pm}=0,$ (20)

where the solutions in the new basis $x_{\pm}(t)$ exhibit real and imaginary
loss terms. We then let

$x_{\pm}(t)=e^{-\frac{(g\pm ir)\omega_{0}t}{2}}y_{\pm}(t)$ (21)

for $y_{\pm}(t)$ to be determined. Upon substituting Equation (21) into the
decoupled system of Equation (20), we arrive at:

$\ddot{y}_{\pm}+\omega_{0}^{2}\left[1-\frac{(g\pm ir)^{2}}{4}+h\sin(\gamma
t)\right]y_{\pm}(t)=0.$ (22)

Since the coefficient of $y_{\pm}$ in Equation (22) is periodic with period
$T=\frac{2\pi}{\gamma}$, then according to Floquet Theory, we may write

$y_{\pm}(t)=e^{-i\mu t}f_{\pm}(t),$ (23)

where $f_{\pm}(t)$ are some periodic functions with period
$T=\frac{2\pi}{\gamma}$. Wherever $\text{Im}(\mu)<0$, $y_{+}(t)$ will be
unstable as $t\rightarrow\infty$. Since $f_{\pm}(t)$ are periodic, we may
express them as Fourier series

$f_{\pm}(t)=\sum_{n\in\mathbb{N}}A_{n}^{\pm}e^{in\gamma t},$ (24)

where $A_{n}^{\pm}$ denotes the amplitude of the nth Fourier coefficient for
$f_{\pm}(t)$, respectively. If we substitute Equation (24) into Equation (22)
and collect the coefficients of $e^{in\gamma t}$, we get the recursion
relation for the Fourier coefficients $A_{n}^{\pm}$ as derived by [6]:

$\begin{gathered}D_{\pm,n}(\mu)A_{n}^{\pm}+i\frac{\omega_{0}h^{2}}{2}(A_{n+1}^{\pm}-A_{n-1}^{\pm})=0,\end{gathered}$
(25)

where
$D_{n}(\mu)=\omega_{0}^{2}-(n\gamma-\mu^{2})+i\omega_{0}(g+ir)(n\gamma-\mu)$.
We note that if the driving intensity $h$ is less than one, $A_{n}^{\pm}$ is
coupled most strongly to $A_{n\pm 1}^{\pm}$, with coupling proportional to
$h^{2}$. Without loss of generality, we will consider $n=0$. Therefore, one
recovers a matrix equation of the form:

$\text{Det}\begin{pmatrix}\textbf{D}_{-1}(\mu)&\textbf{M}(h)&0\\\
\textbf{M}^{*}(h)&\textbf{D}_{0}(\mu)&\textbf{M}(h)\\\
0&-\textbf{M}^{*}(h)&\textbf{D}_{+1}(\mu)\end{pmatrix}=0,$ (26)

where one defines

$\begin{split}\textbf{D}_{n}(\mu)&=\begin{pmatrix}D_{+,n}(\mu)&0\\\
0&\textbf{D}_{-,n}(\mu)\end{pmatrix}\\\
\textbf{M}(h)&=\begin{pmatrix}\frac{ih\omega_{0}^{2}}{2}&0\\\
0&\frac{ih\omega_{0}^{2}}{2}\end{pmatrix}.\end{split}$ (27)

Note that in Equation (26) that enforcing the determinant to vanish ensures
that there are non-trivial solutions to Equation (22). In contrast, it was
also shown in [6] that the simplifying assumptions of $r=g=0$ and ignoring the
coupling effects of $D_{+1}$ yield the matrix equation:

$\text{Det}\begin{pmatrix}\omega_{0}^{2}-(\gamma-\mu)^{2}&\frac{ih\omega_{0}^{2}}{2}\\\
-\frac{ih\omega_{0}^{2}}{2}&\omega_{0}^{2}-\mu^{2}\\\ \end{pmatrix}=0,$ (28)

The matrix equation in Equation (28) is a consequence of assuming the coupling
effects as well as the higher order effects are negligible in the asymptotic
analysis. Hence, we hereafter refer to the solutions of Equation (28), which
is the same as presented in [6], as the simplified, uncoupled solution.

Similar to [6], we focus on the solutions of Equation (26) and Equation (28)
that display parametric resonance at
$\gamma=2\Omega_{r},2\Omega_{r}\pm\omega_{0}r$; we present the solutions in
the Results section. We show that the simplifying assumptions have highly non-
trivial and significant effects on the system’s regions of stability at and
away from parametric resonance frequencies.

##  5 Results

We first compare the stability boundary of the results established via
Equation (26) with the stability boundary of the simplified, uncoupled
solution of Equation (28). We then compare the stability boundaries as
predicted by asymptotic analysis with the stability boundary as predicted by
the auxillary function method for long time averages.

###  5.1 Influence of Higher Order Effects

Indeed, we employ MATLAB’s Solve function to find the roots in $\mu$ of
Equation (26), whose solutions are sufficiently complicated and long to not be
included in the text 444The expressions for the roots can be found on the lead
author’s web-page: sites.google.com/andrewmcmillan/Research.; This yields six
solutions. In Figure 1 555The authors remark that implicitly solving for the
roots of $\mu$ is computationally sensitive to very small numerical error.
Hence, the stability boundary seen in Figure 1 had miniscule variations
depending on the software used. However, the general shape of the stability
boundary is consistent across software., we plot the contours of the imaginary
part of one of the combined boundaries of the six solutions of our more
generalized result in comparison to the simplified, uncoupled solution to
Equation (28) for a wide range of parameter values:

Figure 1: The simplified, uncoupled asymptotic solution to Equation (28)
(shown with colored, solid contour lines) and the more general asymptotic
solution to Equation (26) (shown with dashed, gray contour lines). The markers
$\text{I}_{\text{S}}$, $\text{I}_{\text{U}}$, $\text{II}_{\text{S}}$, and
$\text{II}_{\text{U}}$ denote the stable (S) or unstable (U) regions of the
simplified (I) and general (II) solutions, respectively. The parameter values
are $r=.2$, $g=.01$, and $\phi=0$.

In Figure 1, we see that the more general, higher order asymptotic analysis
solution is a dramatic and surprising improvement on the simplified, uncoupled
solution, as the region of instability for the simplified, uncoupled
asymptotic result is a subset of the region of instability for the more
general asymptotic result. While there is good agreement between the two
solutions for $|\gamma|>1$, there is a large region within the parameter
space, $|\gamma|\leq 1$, for which the simplified, uncoupled asymptotic
analysis solution predicts stability, while the more general, higher order
solution reveals instability.

The fact that the simplified, uncoupled asymptotic expansion stability
solution fails so drastically for $|\gamma|\leq 1$ is quite a surprising
finding. This means the simplified solutions are dangerously un-conservative
in predicting regions of stability. However, it then becomes natural to ask if
the more general, higher order solution also fails to fully describe the
stability boundary of the solutions to Equation (11). That is, we may ask if
including successively higher order terms in the asymptotic expansion would
lead to such drastic changes in the stability boundary as seen between the two
solutions in Figure 1. In order to address this question, we compare the
stability regions as predicted by the general asymptotic analysis with the
stability regions predicted by the auxillary function method for long time
averages.

###  5.2 Asymptotic Analysis vs. The auxillary function method

Upon computing the stability boundary with the auxillary function method for
long time averages using the expression in Equation (13), we arrive at the
solid black line seen Figure 2:

Figure 2: There is strong agreement between our auxillary function method
(shown with a thick, solid black line) with the theoretical predictions of the
simplified, uncoupled asymptotic analysis (shown with colored, solid contour
lines) and the more general, higher order asymptotic analysis (shown with
dashed, gray contour lines) for $r=.2$, $g=.01$, and $\phi=0$. The open blue
circles and red crosses indicate where direct numerical simulation results are
shown in Figure 3 for stable and unstable points, respectively.

In Figure 2, we see that the more general asymptotic analysis result agrees
quite well with the results of the auxillary function method. In particular
and most noticeably, both methods capture a pair of narrow, protruding tongues
that occur at $\gamma\in\\{-1,1\\}$. However, there is still a quite large
range of $\gamma$ values for which the auxillary function method predicts
potential instability, but the asymptotic analysis solution to Equation (26)
predicts stability.

In order to validate the extended region of instability as predicted by the
auxillary function method for long time averages results, we choose four
points within the parameter space to perform direct numerical simulation using
ode45 in MATLAB.The results show that the auxillary function method for long
time averages is corroborated via direct numerical simulation.

Figure 3:  The time histories predicted by DNS for the four points marked in
Figure 2. The left two plots correspond to the two blue, open circles just
outside the region of instability for Figure 2, and the two right plots
correspond to the two red crosses just inside the region of instability for
Figure 2. In the above plot, an initial condition vector IC: [a,b,c,d]
corresponds to $x_{1}(0)=a$, $x_{2}(0)=b$, $y_{1}(0)=c$, and $y_{2}(0)=d$ for
Equation (17) with $r=0.2$, $g=0.01$, and $\phi=0$.

That is, for points marked by red crosses inside the region of instability,
there is a blow up of the solutions, with the time histories shown on the
right hand side of Figure 3. However, for points marked by open, blue circles
just outside the instability region experience decay, with the time histories
shown on the left hand side of Figure 3.

Additionally, Figure 3 displays three trajectories for three randomly chosen
initial conditions for each of the four chosen points in parameter space. We
also note that the above figure only plots the trajectory of $x_{1}(t)$ for
our model, but when we plotted the trajectories for $x_{2}(t)$ a similar
pattern of blow-up or stability still held except that $x_{2}(t)$ blew up at a
rate much slower in comparison to $x_{1}(t)$. This is to be expected as the
coupling term appears negatively in the equation for $x_{1}(t)$ and positively
in the equation for $x_{2}(t)$. Hence, the $x_{2}(t)$ oscillator feeds energy
into the $x_{1}(t)$ oscillator, so we would expect a faster blow-up for
$x_{1}(t)$.

This is perhaps a case study example of how the auxillary function method can
lend itself to finding regions of instability where perturbative methods fail.
Strikingly, the perturbative method for this system seems to fail quite
drastically with a large portion of the stability diagram not being captured
by the perturbative method but instead by the auxillary function method.
Therefore, the auxillary function method has the advantage of being able to
recover the true stability boundary both at and away from parametric
resonance, while doing so in a computationally efficient way. We remark that
the auxillary function method is quite computationally efficient despite the
polynomial representation, seen in Equation (12) , containing five degrees of
freedom. In particular, the CPU times corresponding to the computation of
$\overline{\phi^{*}}$ at the points
$(\gamma,h)=(-1.13,.650),(-1.00,.650),(-.660,.390),$ and $(-.660,.490)$ are
3.4844, 2.6744, 2.8574, and 2.5737 seconds for a degree 8 auxillary function
and a degree 6 S-procedure enforcement on nothing more than a standard laptop
using a single core, 2.2 GHZ processor666The source code for the various
points is also available on the lead author’s web-page..

In light of the above results, it then becomes natural to analyze how varying
the system parameters changes the above regions of stability. In particular,
we study the effects of varying the coupling strength in the next section.

###  5.3 Effects of Damping and Coupling Parameters

In order to study the affects of varying the coupling strength, we consider
Equation (21) for $r=.4$ and compare the stability boundary with the prior
results of the both auxillary function method as well as the asymptotic
analysis for $r=.2$ and $g=.01$. We do not consider negative values for $r$
because the coupling term appears negatively in the first equation and
positively in the second equation of Equation (21); hence, a sign change of
$r$ just swaps the role of $x_{1}(t)$ and $x_{2}(t)$, respectively. Also, we
focus on $\gamma\in[-1.5,1.5]$, as this is where the discrepancies between the
auxillary function method and asymptotic analysis were established in the
previous section.

Performing the same procedure with the auxillary function method as in the
previous section as well as incorporating the predictions of the asymptotic
analysis, we find in Figure 4 that the region of stability for the equations
with $r=.2$ is a subset for the region of the stability for the equations with
$r=.4$. Intuitively, one would expect this to be true because as the r-value
increases, the rate of energy transfer from one oscillator to another
increases and hence initial transients are far less likely to become unstable.

Figure 4: The stability boundary as predicted by the auxillary function method
is shown with a thick, solid black line and a thick, dotted gray line for
$r=.2$ and $r=.4$, respectively, while the asymptotic analysis boundary is
shown with a thin, purple dashed line and a thin, blue dotted line for $r=.2$
and $r=.4$, respectively, with $g=.01$ and $\phi=0$.

Figure 4 shows that the discrepancies between the auxillary function method
and the asymptotic analysis solution for $r=.4$ persist. Moreover, the
discrepancies seem to grow as the value of $r$ is increased. Therefore, as the
coupling strength between the oscillators increases, the asymptotic analysis
solution to Equation (26) increasingly fails to capture the true stability
boundary of the system, because the coupling terms amplify the relative
importance of the higher order terms, particularly away from parametric
resonance as the sine and cosine terms do not cancel. We also note that Figure
4 is for $g=.01$, but if one is interested in the affects of $g$, we can
similarly compute regions of stability for fixed $r$ and different g values.
The expected result would simply be a shift up of the region of stability
appearing in Figure 1. The reason is that $g$ appears on the damping term and
hence increasing the damping coefficient should increase the size of the
region of stability.

###  5.4 Effects of the Parametric Phase

Finally, we study the effects of varying the phase value $\phi$. In a fashion
similar to [6], we focus on $\phi=0,\frac{\pi}{2},\pi$. However, it is
important to remark that [6] was unable to determine a closed form, asymptotic
solution for varying $\phi$. We also could not find a closed form solution to
the analog of Equation 26 with $\phi\neq 0$. The ability to find the stability
boundary for general system parameters proves to be another benefit of the
auxillary function method. Indeed, performing the analogous computations via
the auxillary function method and choosing the same $\Phi$ that appears in
Equation (13), we arrive at Figure 5:

Figure 5: The stability boundary as predicted by the auxillary function method
for $\phi=0$ (shown with a thick, solid black line), $\phi=\frac{\pi}{2}$
(shown with a thin, solid blue line), and $\phi=\pi$ (shown with a thin,
dashed gray line).

Figure 5 shows that the auxillary function method reproduces the protruding
tongues around $\gamma=\pm 2\Omega_{r}$ and $\gamma=\pm
2\Omega_{r}\pm\omega_{0}r$ with $\phi=0$, $\phi=\pi$, and $\phi=\frac{\pi}{2}$
corresponding to one, two, and three tongues respectively. Moreover, note that
varying $\phi$ has non-trivial effects throughout the entirety of parameter
space, and not only about $\gamma=\pm 2\Omega_{r},\pm
2\Omega_{r}\pm\omega_{0}r$, as varying the phase of the parametric oscillator
term can lead to more localized tongue formation.

Since the phase of a parametric oscillator can often be random, it is
important to identify the instability boundary when varying phases are taken
into consideration, especially when there is non-trivial disagreement between
the true stability boundary and the boundary as predicted by asymptotic
analysis.

##  6 Conclusion and Discussion

In summary, we have investigated the higher order effects caused by coupling
parameters on the stability boundary of a non-linear, parametrically driven
coupled oscillator system across a broad range of modulation frequencies. We
have shown that the stability boundary of a non-linear, parametrically driven
oscillator system as predicted by simplified, second order asymptotic methods
ignoring the coupling terms between the oscillators can differ quite
substantially at modulation frequencies away from parametric resonance
frequencies. The simplified, asymptotic solution is un-conservative compared
to a full, second order asymptotic solution when coupling terms are not
ignored, which can in turn differ and still be non-trivially un-conservative
in comparison to the true stability boundary at modulation frequencies away
from the parametric resonance frequencies. The differences are caused by both
neglecting the coupling terms and higher order effects.

This is a primary drawback of critical importance for asymptotic methods. That
is, the validity of the asymptotic results depends pivotally on one’s choice
of the approximations. However, there is currently no way to know which order
will be sufficient to capture the true instability boundary, which depends on
the strength of the coupling terms, $g$ and $r$, as well as the phase of the
parametric oscillator term $\phi$. The solution is expected to be even more
sensitive to these system parameters if the resonance frequency of the
individual oscillations are different or if additional non-linearities are
present.

Hence, we have shown that the auxillary function method for long-time averages
is an efficient and robust means of computing the true long-time averages and
true regions of stability across all possible initial conditions without the
need of ad-hoc approximations. Moreover, this auxillary function method has
the advantage of being able to compute regions of stability both at and away
from parametric resonance.

The differences between the true stability boundary and the approximate
stability boundary are immaterial if one is operating withing a region of
parameter space where both boundaries agree. However, if one is operating
within a region of parameter space for which they disagree, this may be quite
problematic for both experimental and real-world implementation- especially
without the knowledge of operating within one of these regions of
disagreement. This point is exacerbated by the fact that, at least in studying
Equation (11), we discovered two, very narrow protruding tongues in parameter
space for which the system was potentially unstable. Moreover, these tongues
occur at very naturally occuring, non-pathological values of the driving
frequency $\gamma$.

In the context of applications, our results have profound implications on the
reliance of asymptotic methods across a variety of applied sciences. If one
was interested in establishing a dynamical system via data-driven methods, the
instability of these very narrow and protruding regions of parameter space
could be completely missed, especially if data sampling was on the order of
the size of the narrow tongues. If one was interested in studying synchrony
behavior or global phase locking, asymptotic methods may fail to capture
sensitivity to a system’s parameters.

Hence, it is recommended that the instability region be determined and
considered when planning new experiments or full-scale implementation, as
disregarding the stability boundary or only relying on ad-hoc approximations
may have costly consequences.

## Dedication

We are greatly indebted to Dr. Charles Doering (1956- 2021)- a great mentor,
friend, and colleague. His insights, commentary and suggestions immensely
impacted this work, and the entire scientific community will continue to miss
his absence.

## Acknowledgments

Support for this research is provided by the U.S. Office of Naval Research,
contract No. N00014-18-C-1025 managed by Ms. Deborah Nalchajian.

## Author Declarations

The authors acknowledge no personal, professional, or financial conflicts of
interest.

## References

  * [1] Gonze D, Bernard S, Waltermann C, Kramer A, and Herzel H. (2005). Spontaneous synchronization of coupled circadian oscillators. J. Biophysics .
  * [2] Rosa, R., and Temam, R. (2020). Optimal minimax bounds for time and ensemble averages of dissipative infinite-dimensional systems with applications to the incompressible Navier-Stokes equations.
  * [3] Goulart, P., and Chernyshenko, S. (2021). Global stability analysis of fluid flows using sum-of-squares. Physica D: Nonlinear Phenomena, 241, 692-704.
  * [4] MOSEK ApS (2015). The MOSEK optimization toolbox for MATLAB manual.
  * [5] Lofberg, J. (2004). YALMIP : A toolbox for modeling and optimization in MATLAB. IEEE International Conference on Robotics and Automation.
  * [6] Marcello Calvanese Strinati, Leon Bello, Avi Pe’er, and Emanuele G. Dalla Torre (2019). Theory of coupled parametric oscillators beyond coupled Ising spins. Physical Review A 100.
  * [7] Cottone, F., Vocca, H., and Gammaitoni, L. (2009). Nonlinear Energy Harvesting. Phys. Rev. Lett., 102, 080601.
  * [8] A. Erturk, D. (2011). Broadband piezoelectric power generation on high-energy orbits of the bistable Duffing oscillator with electromechanical coupling. Journal of Sound and Vibration, 330, 2339-2353.
  * [9] Gammaitoni, L., Neri, I., and Vocca, H. (2009). Nonlinear oscillators for vibration energy harvesting. Applied Physics Letters, 94(16), 164102.
  * [10] Triplett, A. and D. Dane Quinn (2009). The Effect of Non-Linear Piezoelectric Coupling on Vibration-Based Energy Harvesting. Journal of Intelligent Material Systems and Structures, 20(16), 1959–67.
  * [11] E. Ott, C. and J.A. Yorke (1990). Controlling Chaos. Physical Review Letters, 64, 1196–1199.
  * [12] I. Tobasco, D. Goluskin, and C.R. Doering (2018). Optimal bounds and extremal trajectories for time averages in nonlinear dynamical systems. Physics Letters A, 208, 382–386.
  * [13] A.N. McMillan ad C.R. Doering (2018). Optimal time averages in non-autonomous nonlinear dynamical systems. Journal of Pure and Applied Functional Analysis, 208, 382–386.
  * [14] E. Esmailzadeh, and G. Nakhaie-jazar (1997). Periodic solution of a Mathieu-Duffing type equation. International Journal of Non-Linear Mechanics, 32, 905–912.
  * [15] Sethia, G., Sen, A., and Atay, F. (2008). Clustered Chimera States in Delay-Coupled Oscillator Systems. Phys. Rev. Lett., 100, 144102.
  * [16] Fujisaka, H., and Yamada, T. (1983). Stability Theory of Synchronized Motion in Coupled-Oscillator Systems. Progress of Theoretical Physics, 69.
  * [17] Heagy, J., Carroll, T., and Pecora, L. (1994). Synchronous chaos in coupled oscillator systems. Phys. Rev. E, 50, 1874–1885.
  * [18] Wilson, C., and Callaway, J. (2021). Coupled Oscillator Model of the Dopaminergic Neuron of the Substantia Nigra. Journal of Neurophysiology, 83(5), 3084–3100.
  * [19] D. Eliyahu, K. Sariri, A. Kamran, and M. Tokhmakhian (2002). Improving short and long term frequency stability of the opto-electronic oscillator. In Proceedings of the 2002 IEEE International Frequency Control Symposium and PDA Exhibition (Cat. No.02CH37234) (pp. 580-583).
  * [20] Reppert, S., and David R. Weaver (2002). Coordination of Circadian Timing in Mammals. Nature, 418(6901), 935-941.
  * [21] Komin, N., Murza A. C., Hernandez-Garcia E., and Toral R. (2011). Synchronization and Entrainment of Coupled Circadian Oscillators.. Interface Focus, 1(1), 167-176.
  * [22] Engelbrecht, J., and Mirollo, R. (2012). Structure of Long-Term Average Frequencies for Kuramoto Oscillator Systems. Phys. Rev. Lett., 109, 034103.
  * [23] Schwemmer, Michael A. and Lewis, Timothy J., The Theory of Weakly Coupled Oscillators, 2012,Springer New York, 3-31, 978-1-4614-0739-3, doi=10.1007/978-1-4614-0739-31.
  * [24] Strogatz, S. (2000). From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D Nonlinear Phenomena, 143(1-4), 1-20.
  * [25] E. Ott, C., and J.A. Yorke (1990). Controlling Chaos. Physical Review Letters, 64, 1196-1199.
  * [26] Beneke, M., and Zakharov, V. (1992). Improving large-order perturbative expansions in quantum chromodynamics. Phys. Rev. Lett., 69, 2472–2474.
  * [27] E. de Langre (2006). Frequency lock-in is caused by coupled-mode flutter. Journal of Fluids and Structures, 22(6), 783-791.
  * [28] Xiao-Jun Yang, and H.M. Srivastava (2015). An asymptotic perturbation solution for a linear oscillator of free damped vibrations in fractal medium described by local fractional derivatives. Communications in Nonlinear Science and Numerical Simulation, 29(1), 499-504.
  * [29] William E Caswell (1979). Accurate energy levels for the anharmonic oscillator and a summable series for the double-well potential in perturbation theory. Annals of Physics, 123(1), 153-184.
  * [30] Jayanta K. Ghosh (1994). Higher Order Asymptotics. NSF-CBMS Regional Conference Series in Probability and Statistics, 4, i–111.
  * [31] Goluskin, D., and Fantuzzi, G. (2019). Bounds on mean energy in the Kuramoto–Sivashinsky equation computed using semidefinite programming. Nonlinearity, 32, 1705-1730.
  * [32] Peter Miller (2006). Applied Asymptotic Analysis. American Mathematical Society.
  * [33] Lifschitz, R., and Cross, M. (2003). Response of parametrically driven nonlinear coupled oscillators with application to micromechanical and nanomechanical resonator arrays. Physical Review B, 67.
  * [34] Matteo Corbetta (2020). Application of sparse identification of nonlinear dynamics for physics-informed learning. Journal of Neurophysiology, 83(5), 3084–3100.
  * [35] Csaba, G., and Porod, W. (2020). Coupled oscillators for computing: A review and perspective. Applied Physics Reviews, 7(1), 011302.
  * [36] Zheng, Z., Hu, G., and Hu, B. (1998). Phase Slips and Phase Synchronization of Coupled Oscillators. Phys. Rev. Lett., 81, 5318–5321.
  * [37] Strogatz, S., and Mirollo, R. (1991). Stability of incoherence in a population of coupled oscillators. J Stat Phys, 63, 613-635.
  * [38] Fink, K., Johnson, G., Carroll, T., Mar, D., and Pecora, L. (2000). Three coupled oscillators as a universal probe of synchronization stability in coupled oscillator arrays. Phys. Rev. E, 61, 5080–5090.
  * [39] Ermentrout, G. (1985). The behavior of rings of coupled oscillators. J. Math. Biologu, 23, 55-74.
  * [40] Kyongmin Yeo (2019). Data-driven reconstruction of nonlinear dynamics from sparse observation. Journal of Computational Physics, 395, 671-689.
  * [41] Champion, K., Lusch, B., Kutz, J., and Brunton, S. (2019). Data-driven discovery of coordinates and governing equations. Proceedings of the National Academy of Sciences, 116(45), 22445–22451.


is combinatorially flat. Since the first map in the factorization above is
strict, the map

${\prod_{e}}_{\bigtimes_{v}\mathsf{GW}(X_{v})}\mathcal{D}_{e}^{\\{2\\}}\to\bigtimes_{v}\mathsf{Ev}^{\dagger}_{\gamma}(v)$

is also combinatorially flat.

Let $\Delta^{\dagger}$ denote the strict transform of the diagonal:

$\prod_{e}D_{e}\hookrightarrow\prod_{e}D_{e}\times
D_{e}=\prod_{v}\mathsf{Ev}_{\gamma}(v),$

under the blowup

$\bigtimes_{v}\mathsf{Ev}_{\gamma}(v)\to\prod_{v}\mathsf{Ev}_{\gamma}(v)$

######  Proposition 12.3.1.

The preimage of $\Delta^{\dagger}$ under the map

${\prod_{e}}_{\bigtimes_{v}\mathsf{GW}(X_{v})}\mathcal{D}_{e}^{\\{2\\}}\to\bigtimes_{v}\mathsf{Ev}_{\gamma}(v).$

is the relative diagonal locus, i.e. the fibered product of
$\mathcal{D}_{e}^{\Delta}$ over all edges $e$, with the fiber product taken
over $\bigtimes_{v}\mathsf{GW}(X_{v})$.

###### Proof.

Since the map is combinatorially flat, the preimage coincides with the fine
and saturated logarithmic pullback filling in:

${{\prod_{e}}_{\bigtimes_{v}\mathsf{GW}(X_{v})}\mathcal{D}_{e}^{\\{2\\}}}$${\prod_{e}D_{e}}$${\prod_{e}D_{e}\times
D_{e},}$

where the lower map is the factorwise diagonal. It therefore suffices to prove
the statement by interpreting the diagram in the category of stacks over
logarithmic schemes. The claim now follows by unwinding the functor of points
of the fine and saturated base change. ∎

As a consequence of the proposition, we obtain the first refinement of the
gluing formula. We continue with the notation $\mathcal{Y}$ for the
degeneration, $Y_{0}$ for its special fiber, and $\gamma$ for the rigid
tropical map.

In what follows, we fix subdivisions such that there exists a combinatorially
flat evaluation map:

$\mathsf{ev}^{\sim}:\bigtimes_{v}\mathsf{GW}(X_{v})\to\bigtimes_{v}\mathsf{Ev}^{\dagger}_{\gamma}(v).$

As we have used before, we let $\ell$ denote the product, taken over all edges
$e$ in $\overline{\gamma}$, of the least common multiplies of the dilation
factors $m_{e^{\prime}}$ of edges $e^{\prime}$ in $\gamma$ that map to $e$.

######  Corollary 12.3.2.

Given the morphism:

$q:\mathsf{GW}_{\gamma}(Y_{0})\to\bigtimes_{v}\mathsf{GW}(X_{v}),$

there is an equality of Chow homology classes

$q_{\star}\left[\mathsf{GW}_{\gamma}(Y_{0})\right]^{\mathsf{vir}}=\prod_{e}m_{e}\cdot\frac{1}{\ell(\gamma)}\cdot\left(\mathsf{ev}^{\sim}\right)^{\star}\Delta^{\dagger}\cap\left[\bigtimes_{v}\mathsf{GW}(X_{v})\right]^{\mathsf{vir}}.$

We now present one final simplification of the formula. The result is parallel
to Theorem 8.2.2 on the DT side. Choose any combinatorial flattening:

$\mathsf{ev}^{\dagger}:\prod_{v}\mathsf{GW}^{\dagger}_{\gamma}(X_{v})\to\prod_{v}\mathsf{Ev}_{\gamma}^{\dagger}(v).$

The number $\ell$ is as above.

######  Theorem 12.3.3.

Given the morphism

$\kappa:\mathsf{GW}^{\dagger}_{\gamma}(Y_{0})\to\prod_{v}\mathsf{GW}^{\dagger}_{\gamma}(X_{v}),$

there is an equality of Chow homology classes

$\kappa_{\star}\left[\mathsf{GW}^{\dagger}_{\gamma}(Y_{0})\right]^{\mathsf{vir}}=\prod_{e}m_{e}\cdot\frac{1}{\ell(\gamma)}\cdot\left(\mathsf{ev}^{\dagger}\right)^{\star}\Delta^{\dagger}\cap\left[\prod_{v}\mathsf{GW}^{\dagger}(X_{v})\right]^{\mathsf{vir}}.$

We stress that the key difference between the corollary that precedes this
result, and the theorem we have just stated, is that the second formula splits
into a vertex-by-vertex expression.

###### Proof.

The proof is exactly parallel to the DT side. Since it is important for our
main results, we nevertheless sketch the steps. First, we have proved in the
preceding corollary that there exists a combinatorially flat evaluation map,
namely $\mathsf{ev}^{\sim}$, such that the pullback of the strict transform of
the diagonal gives the required class associated to the rigid tropical curve
$\gamma$. Second, we move this from a statement about a particular
combinatorial flattening of the evaluation map to an arbitrary cone. We have
fixed such an arbitrary flattening from the outset:

$\mathsf{ev}^{\dagger}:\prod_{v}\mathsf{GW}^{\dagger}_{\gamma}(X_{v})\to\prod_{v}\mathsf{Ev}_{\gamma}^{\dagger}(v).$

We can also assume by performing further blowups, that there is a subdivision

$\bigtimes_{v}\mathsf{Ev}_{\gamma}^{\dagger}(v)\to\prod_{v}\mathsf{Ev}_{\gamma}^{\dagger}(v)$

We therefore end up with a commutative diagram:

${\bigtimes_{v}\mathsf{GW}_{\gamma}(X_{v})}$${\prod_{v}\mathsf{GW}^{\dagger}_{\gamma}(X_{v})}$${\bigtimes_{v}\mathsf{Ev}_{\gamma}^{\dagger}(v)}$${\prod_{v}\mathsf{Ev}_{\gamma}^{\dagger}(v).}$

Although this map is not commutative, we can replace the top left part of the
diagram with the fine and saturated fiber product: the result will still be a
modification of the product since it is a pullback of a subdivision. Formal
properties of the virtual class under pushforward guarantee that the
conclusion of Corollary 12.3.2 continues to hold with this replacement. Now,
since the right vertical is combinatorially flat, this fine and saturated
fibered diagram is also fibered in the category of stacks. We are now in
exactly the situation of the proof of Theorem 8.2.2 in the previous section.
After making the appropriate cosmetic changes, replacing the Hilbert schemes
of curves and points with the spaces of stable maps and their evaluation
spaces, the argument follows from a formally identical diagram chase. We
conclude. ∎

The discussion in this section concludes the simplified gluing formula that we
have input into the key corollaries of Section 8.

### References

  * [1] D. Abramovich and Q. Chen, Stable logarithmic maps to Deligne-Faltings pairs II, Asian J. Math., 18 (2014), pp. 465–488.
  * [2] D. Abramovich, Q. Chen, M. Gross, and B. Siebert, Decomposition of degenerate Gromov-Witten invariants, arXiv:1709.09864, (2017).
  * [3] D. Abramovich, Q. Chen, M. Gross, and B. Siebert, Punctured logarithmic curves, In preparation, (2017).
  * [4] D. Abramovich, Q. Chen, S. Marcus, M. Ulirsch, and J. Wise, Skeletons and fans of logarithmic structures, in Nonarchimedean and Tropical Geometry, M. Baker and S. Payne, eds., Simons Symposia, Springer, 2016, pp. 287–336.
  * [5] D. Abramovich and B. Fantechi, Orbifold techniques in degeneration formulas, arXiv:1103.5132, (2011).
  * [6] D. Abramovich and K. Karu, Weak semistable reduction in characteristic 0, Invent. Math., 139 (2000), pp. 241–273.
  * [7] D. Abramovich, S. Marcus, and J. Wise, Comparison theorems for Gromov–Witten invariants of smooth pairs and of degenerations, in Ann. Inst. Fourier, vol. 64, 2014, pp. 1611–1667.
  * [8] D. Abramovich and J. Wise, Birational invariance in logarithmic Gromov–Witten theory, Comp. Math., 154 (2018), pp. 595–620.
  * [9] K. Adiprasito, G. Liu, and M. Temkin, Semistable reduction in characteristic $0$, arXiv:1810.03131, (2018).
  * [10] A. A. Beĭlinson, J. Bernstein, and P. Deligne, Faisceaux pervers, in Analysis and topology on singular spaces, I (Luminy, 1981), vol. 100 of Astérisque, Soc. Math. France, Paris, 1982, pp. 5–171.
  * [11] P. Bousseau, Tropical refined curve counting from higher genera and lambda classes, arXiv:1706.07762, (2017).
  * [12] , A proof of n. takahashi’s conjecture for $(\mathbb{P}^{2},e)$ and a refined sheaves/gromov-witten correspondence, arXiv:1909.02992, (2019).
  * [13] P. Bousseau, H. Fan, S. Guo, and L. Wu, Holomorphic anomaly equation for $(\mathbb{P}^{2},E)$ and the Nekrasov-Shatashvili limit of local $\mathbb{P}^{2}$, Forum Math. Pi, 9 (2021), pp. Paper No. e3, 57.
  * [14] J. Bryan and R. Pandharipande, The local Gromov-Witten theory of curves, J. Amer. Math. Soc., 21 (2008), pp. 101–136.
  * [15] F. Carocci and N. Nabijou, Rubber tori in the boundary of expanded stable maps, arXiv preprint arXiv:2109.07512, (2021).
  * [16] R. Cavalieri, M. Chan, M. Ulirsch, and J. Wise, A moduli stack of tropical curves, Forum Math. Sigma, 8 (2020), pp. 1–93.
  * [17] Q. Chen, The degeneration formula for logarithmic expanded degenerations, J. Algebr. Geom., 23 (2014), pp. 341–392.
  * [18] , Stable logarithmic maps to Deligne-Faltings pairs I, Ann. of Math., 180 (2014), pp. 341–392.
  * [19] M. A. A. de Cataldo and L. Migliorini, The Chow groups and the motive of the Hilbert scheme of points on a surface, J. Algebra, 251 (2002), pp. 824–848.
  * [20] , The Chow motive of semismall resolutions, Math. Res. Lett., 11 (2004), pp. 151–170.
  * [21] W. Fulton, Intersection theory, vol. 2 of Ergebnisse der Mathematik und ihrer Grenzgebiete. 3. Folge. A Series of Modern Surveys in Mathematics [Results in Mathematics and Related Areas. 3rd Series. A Series of Modern Surveys in Mathematics], Springer-Verlag, Berlin, second ed., 1998.
  * [22] A. Gathmann and H. Markwig, The numbers of tropical plane curves through points in general position, J. Reine Angew. Math. (Crelle’s Journal), 2007 (2007), pp. 155–177.
  * [23] L. Göttsche and W. Soergel, Perverse sheaves and the cohomology of Hilbert schemes of smooth algebraic surfaces, Math. Ann., 296 (1993), pp. 235–245.
  * [24] I. Grojnowski, Instantons and affine algebras. I. The Hilbert scheme and vertex operators, Math. Res. Lett., 3 (1996), pp. 275–291.
  * [25] M. Gross, Remarks on gluing punctured logarithmic maps, arXiv preprint arXiv:2306.02661, (2023).
  * [26] M. Gross and B. Siebert, Logarithmic Gromov-Witten invariants, J. Amer. Math. Soc., 26 (2013), pp. 451–510.
  * [27] , Intrinsic mirror symmetry, arXiv:1909.07649, (2019).
  * [28] D. Holmes, S. Molcho, R. Pandharipande, A. Pixton, and J. Schmitt, Logarithmic double ramification cycles, arXiv:2207.06778, (2022).
  * [29] M. Kapranov and E. Vasserot, The cohomological hall algebra of a surface and factorization cohomology, arXiv:1901.07641, (2019).
  * [30] K. Kato, Logarithmic structures of Fontaine-Illusie, Algebraic analysis, geometry, and number theory (Baltimore, MD, 1988), (1989), pp. 191–224.
  * [31] G. Kempf, F. Knudsen, D. Mumford, and B. Saint-Donat, Toroidal embeddings I, Lecture Notes in Mathematics, 339 (1973).
  * [32] P. Kennedy-Hunt, Logarithmic Pandharipande-Thomas spaces and the secondary polytope, arXiv:2112.00809, (2021).
  * [33] P. Kennedy-Hunt, Q. Shafi, and A. U. Kumaran, Tropical refined curve counting with descendants, arXiv:2307.09436, (2023).
  * [34] J. Li, Stable morphisms to singular schemes and relative stable morphisms, J. Diff. Geom., 57 (2001), pp. 509–578.
  * [35] , A degeneration formula of GW-invariants, J. Diff. Geom., 60 (2002), pp. 199–293.
  * [36] J. Li and B. Wu, Good degeneration of Quot-schemes and coherent systems, Communications in Analysis and Geometry, 23 (2015), pp. 841–921.
  * [37] C. Manolache, Virtual pull-backs, J. Algebr. Geom., 21 (2012), pp. 201–245.
  * [38] D. Maulik, N. Nekrasov, A. Okounkov, and R. Pandharipande, Gromov–Witten theory and Donaldson–Thomas theory, I, Compos. Math., 142 (2006), pp. 1263–1285.
  * [39] , Gromov–Witten theory and Donaldson–Thomas theory, II, Compos. Math., 142 (2006), pp. 1286–1304.
  * [40] D. Maulik, A. Oblomkov, A. Okounkov, and R. Pandharipande, Gromov-Witten/Donaldson-Thomas correspondence for toric 3-folds, Invent. Math., 186 (2011), pp. 435–479.
  * [41] D. Maulik, R. Pandharipande, and R. P. Thomas, Curves on K3 surfaces and modular forms, J. Topol., 3 (2010), pp. 937–996.
  * [42] D. Maulik and D. Ranganathan, Logarithmic Donaldson-Thomas theory, preprint, (2020).
  * [43] G. Mikhalkin, Enumerative tropical geometry in ${\mathbb{R}^{2}}$, J. Amer. Math. Soc, 18 (2005), pp. 313–377.
  * [44] S. Molcho, Universal stacky semistable reduction, Israel J. Math., 242 (2021), pp. 55–82.
  * [45] S. Molcho and D. Ranganathan, A case study of intersections on blowups of the moduli of curves, arXiv:2106.15194, (2021).
  * [46] H. Nakajima, Heisenberg algebra and Hilbert schemes of points on projective surfaces, Ann. of Math. (2), 145 (1997), pp. 379–388.
  * [47] G. Oberdieck, Marked relative invariants and GW/PT correspondences, arXiv:2112.11949, (2021).
  * [48] A. Okounkov, Lectures on K-theoretic computations in enumerative geometry, in Geometry of moduli spaces and representation theory, vol. 24 of IAS/Park City Math. Ser., Amer. Math. Soc., Providence, RI, 2017, pp. 251–380.
  * [49] A. Okounkov and R. Pandharipande, The local Donaldson-Thomas theory of curves, Geom. Topol., 14 (2010), pp. 1503–1567.
  * [50] M. C. Olsson, Logarithmic geometry and algebraic stacks, Ann. Sci. Éc. Norm. Supér., 36 (2003), pp. 747–791.
  * [51] R. Pandharipande and A. Pixton, Gromov–witten/pairs descendent correspondence for toric 3–folds, Geometry & Topology, 18 (2014), pp. 2747–2821.
  * [52] , Gromov-witten/pairs correspondence for the quintic 3-fold, J. Amer. Math. Soc., 30 (2017), pp. 389–449.
  * [53] R. Pandharipande and R. P. Thomas, Curve counting via stable pairs in the derived category, Invent. Math., 178 (2009), pp. 407–447.
  * [54] J. Pardon, Universally counting curves in Calabi-Yau threefolds, arXiv:2308.02948, (2023).
  * [55] B. Parker, Tropical enumeration of curves in blowups of the projective plane, arXiv preprint arXiv:1411.5722, (2014).
  * [56] , Gluing formula for gromov-witten invariants in a triple product, arXiv preprint arXiv:1511.00779, (2015).
  * [57] , De rham theory of exploded manifolds, Geometry & Topology, 22 (2017), pp. 1–54.
  * [58] , Three dimensional tropical correspondence formula, Commun. Math. Phys., 353 (2017), pp. 791–819.
  * [59] , Tropical gluing formulae for gromov-witten invariants, arXiv:1703.05433, (2017).
  * [60] , Holomorphic curves in exploded manifolds: virtual fundamental class, Geom. Topol., 23 (2019), pp. 1877–1960.
  * [61] D. Ranganathan, Logarithmic Gromov-Witten theory with expansions, arXiv:1903.09006, (2019).
  * [62] D. Ranganathan and A. U. Kumaran, Logarithmic gromov-witten theory and double ramification cycles, arXiv:2212.11171, (2022).
  * [63] Stacks Project Authors, Stacks project.  http://stacks.math.columbia.edu, 2019.
  * [64] M. van Garrel, D. P. Overholser, and H. Ruddat, Enumerative aspects of the Gross-Siebert program, in Calabi-Yau varieties: arithmetic, geometry and physics, vol. 34 of Fields Inst. Monogr., Fields Inst. Res. Math. Sci., Toronto, ON, 2015, pp. 337–420.
  * [65] J. Wise, Moduli of morphisms of logarithmic schemes, Algebra Number Theory, 10 (2016), pp. 695–735.
  * [66] Y. Wu, Splitting of Gromov-Witten invariants with toric gluing strata, J. Algebr. Geom. (to appear), arXiv:2103.14780, (2021).