###### Contents

  1. I Introduction and Problem Formulation
    1. I-A Properties of Graph Laplacians
    2. I-B Stability of the Mathieu Equation
  2. II Full-network parametric forcing
  3. III Subnetwork Parametric Excitation
    1. III-A Perturbation Analysis
    2. III-B Comparison with Numerical Calculations
    3. III-C Controllability Interpretation
    4. III-D Generalization to Subnetwork Excitation
  4. IV Stability Analysis for Tori Graphs
    1. IV-A Ring Networks with more General Forcing
    2. IV-B Higher-Dimensional Tori
    3. IV-C The Effect of Connectivity
  5. V Conclusion and Further Questions
  6. A Appendix
    1. A-A Vector-valued multiple scale perturbation analysis
      1. A-A1 In the stable region
      2. A-A2 At critical frequencies $\kappa_{0}$, $\sqrt{\kappa_{0}}\left|{\scriptstyle\Omega}_{l}\pm{\scriptstyle\Omega}_{m}\right|=1$
      3. A-A3 Interpreting the multiple-scale analysis stability criteria
    2. A-B The multiple scale method with eigenvalue multiplicity
    3. A-C Multiple scale perturbation analysis: unstable modes
    4. A-D Eigenvector test for controllability

# Parametric Resonance in Networked Oscillators

Karthik Chikmagalur and Bassam Bamieh

###### Abstract

We investigate parametric resonance in oscillator networks subjected to
periodically time-varying oscillations in the edge strengths. Such models are
inspired by the well-known parametric resonance phenomena for single
oscillators, as well as the potential rich phenomenology when such parametric
excitations are present in a variety of applications like deep brain
stimulation, AC power transmission networks, as well as vehicular flocking
formations. We consider cases where a single edge, a subgraph, or the entire
network is subjected to forcing, and in each case, we characterize an
interesting interplay between the parametric resonance modes and the
eigenvalues/vectors of the graph Laplacian. Our analysis is based on a novel
treatment of multiple-scale perturbation analysis that we develop for the
underlying high-dimensional dynamic equations.

##  I Introduction and Problem Formulation

Parametric resonance is a ubiquitous phenomenon occurring at a wide range of
scales, from undesirable roll-motion of large ships [1], to actuated micro-
cantilever systems [2] where parametric resonance is used for high sensitivity
mass measurements. The distinction between this phenomenon and the more
familiar harmonic resonance is that the periodically oscillating signal is a
parameter in the system’s dynamics, rather than an additive input as is the
case with harmonic resonance.

The simplest mathematical model of parametric resonance is the Mathieu or Hill
equation111When $f$ is a pure sinusoid, e.g. $f(t)=a\cos(\omega t)$, it is
referred to as the Mathieu equation, while it is termed Hill’s equation when
$f$ has higher harmonics representing a more general periodic signal. The
qualitative behavior is similar in either case. [3, 4]

$\displaystyle\ddot{x}(t)+\left(k+f(t)\rule{0.0pt}{8.99994pt}\right)x(t)$
$\displaystyle=0$ (1)
$\displaystyle\Leftrightarrow~{}~{}~{}\ddot{x}(t)+k\,x(t)$
$\displaystyle=-f(t)\,x(t)$ (2)

where $f$ is a periodic signal which we refer to as the “parametric
excitation”. When $k>0$, Eq. 1 can be thought of as a Mass-Spring system with
a periodically time-varying spring constant $k+f(t)$. Alternatively, the
parametric excitation term $f(t)\,x(t)$ in Eq. 2 can be thought of as an
applied periodic force on the Mass-Spring system which is modulated (in time)
by the system’s state $x$. The unexcited system (i.e. with $f(t)=0$) has a
natural frequency of $\sqrt{k}$, and when the frequency of the periodic,
parametric excitation $f$ is tuned to certain fractions of this natural
frequency, i.e. at $\sqrt{k}/n,~{}n=\tfrac{1}{2},1,2,3,\ldots$, the system
above exhibits an exponentially growing instability for arbitrarily small
amplitudes of $f$. Thus in contrast to harmonic resonance, parametric
resonance occurs at several frequencies, and the resulting trajectories grow
exponentially rather than linearly. In many models, there are nonlinearities
in Eq. 1 (e.g. stiffening springs) that serve mainly to saturate those
exponentially growing oscillations into stable limit cycles.

When $k<0$ in Eq. 1, this system can be considered as the linearization of the
Kapitza pendulum [5] undergoing vibrational stabilization. In this case,
without parametric excitation (i.e. $f(t)=0$), the system is unstable. For
certain ranges of amplitudes and frequencies of $f$, the system can be
stabilized. This is stabilization without sensor feedback, commonly known as
vibrational stabilization. A compelling argument can be made [5] that
vibrational stabilization and parametric resonance are two different
manifestations of one underlying phenomenon, namely that with parametric,
periodic excitation, the stability of certain systems can be “flipped”, i.e.
from stable to unstable in the parametric resonance case, and from unstable to
stable in the vibrational stabilization case. We note here that since
vibrational stabilization has historically been studied using averaging
methods [6, 7], it is sometimes incorrectly assumed that it requires the use
of high excitation frequencies. In fact, it can occur at any frequency, but
the range of stabilizing amplitudes (of $f$) becoming smaller as the frequency
is reduced. One of the compelling aspects of vibrational control (as well as
parametric resonance) is that it is a form of sensorless control which does
not require feedback in order to change a system’s stability properties, and
may therefore be useful in setting where sensing is either unavailable or
difficult to implement.

The starting point of the present work is the fascinating observation made in
[8, 9] regarding possible mechanisms behind the clinical technique of Deep
Brain Stimulation (DBS). Roughly speaking, the observation is that DBS
essentially modulates the connection strengths between axons and dendrites,
and that DBS is plausibly a “network effect” in an neural network whose edges
have periodically varying strengths. This is then formalized mathematically
[8, 10] as a vibrational control (through periodic edge-strength modulation)
problem of an oscillator network.

In this paper, we take an alternative approach based on studying parametric
resonance, rather than vibrational stabilization in oscillator networks. The
basic model we adapt is for a network of undamped, second-order linear
oscillators of the following form

$\displaystyle\ddot{\phi}(t)+L\,\phi(t)=0,$ (3)

where $L$ is an $n\times n$ graph Laplacian, and $\phi\in\mathbb{R}^{n}$ is a
vector of phases where the $i^{\text{th}}$ component is the phase of the
$i^{\text{th}}$ oscillator.

Examples of such models are the linearized “swing dynamics” of AC power
networks [11], 2nd order models of flocking in vehicular formations [12], and
classical models of molecular vibrations [13]. A useful physical analogy [11]
for the model (3) is that the $ij^{\text{th}}$ entry of the Laplacian $L$ is
the coefficient of a spring-force coupling between the phases $\phi_{i}$ and
$\phi_{j}$.

We study parametric resonance in (3) by assuming that edge strengths (namely
the entries of the Laplacian $L$) are perturbed periodically at a given
frequency. More precisely, the phase dynamics are given by

$\displaystyle\ddot{\phi}+\left(L+\epsilon\,P\,\cos(\omega
t)\rule{0.0pt}{8.99994pt}\right)\phi=0,$ (4)

where $P$ is the graph Laplacian of a subnetwork, and $\epsilon$ and $\omega$
are the strength and frequency of the edge weight forcing respectively. This
subnetwork $P$ can be a single edge, a proper subnetwork, or the full network
itself. Eq. 4 is the higher-order analogue of Mathieu’s equation Eq. 1. By
generalizing the periodic function from a one-harmonic sinusoid to a general
${\scriptstyle T}$-periodic function, we can also obtain the higher-order
analogue of Hill’s

$\displaystyle\ddot{\phi}+\left(L+\epsilon\,P\,f(t)\rule{0.0pt}{8.99994pt}\right)\phi=0,\quad
f(t+{\scriptstyle T})=f(t).$ (5)

In either case, an important assumption we make in this paper is that the the
subnetwork $P$ is perturbed at a _single_ frequency, or by a single scalar
function $f(t)$.

The following two extreme cases serve as a motivation for our analysis.

  * •

$P=L$, the network itself. This is the “global” case, where the entire network
is acted on by the same parametric excitation. This is perhaps the case in DBS
where all neural connection are subject to the same external electromagnetic
field.

  * •

$P=E_{ik}:=e_{ik}e_{ik}^{\star}$, where $e_{ik}$ is the vector with $1$ and
$-1$ in the $i^{\text{th}}$ and $k^{\text{th}}$ positions and $0$ elsewhere.
It is the Laplacian of an $n$-node graph with a single edge from node $i\to
k$. This is the case of a “single line” perturbation as what might happen in
an AC power transmission network where a single line is perturbed by say a
thunderstorm.

We believe this problem formulation is novel and maybe of interest in the
application areas mentioned. In the present paper, we are interested in the
following basic questions about Eq. 4.

  * •

Does the system experience parametric resonance (instability) in the limit
$\epsilon\to 0$?

  * •

How are the destabilizing (“critical”) frequencies of the external excitation
related to the graph structure, $L$?

  * •

Which parametric resonance modes are the most robust or the easiest to excite?

  * •

How does the choice of $P$ affect the stability of the system? Is perturbing
certain edge weights more likely to cause parametric resonance than others?

This paper is organized as follows. We cover the special, simpler case of
"global" or full-network forcing ($P=L$) in Section II. The more challenging
single-line case ($P=E_{ik}$) is treated in Section III using a vector-valued
generalization of multiple-scale perturbation analysis developed in the
appendices. In Section III-D we extend the stability results obtained for
single line perturbations to the more general case of parametric forcing of an
arbitrary subset of network edges. Section III-C provides an alternate
interpretation of the stability results from the perspective of system
controllability. Section IV applies the main stability result to the case of
periodic edge perturbations of ring networks. The simple topology of rings and
periodic lattices (tori) allows us to investigate the effect of network size
and connectivity on its susceptibility to parametric resonance. Finally,
Section V summarizes our findings and concludes with open questions and
possible generalizations of the method.

The next two subsections gather preliminaries on graph Laplacians in Section
I-A used in subsequent analysis, as well as the basic stability properties of
the Mathieu equation in Section I-B.

###  I-A Properties of Graph Laplacians

We briefly review here basic properties of graph Laplacians that we use in our
analysis. We consider only undirected graphs, and thus $L$ or the subgraph
Laplacian $P$ are real symmetric, positive semi-definite with at least one
eigenvalue at $0$ with eigenvector $\mathbbm{1}$, i.e.
$L\mathbbm{1}=0\mathbbm{1}$.

The other eigenvalues give the natural frequencies of oscillation of this
system as the square roots of the eigenvalues

$\displaystyle L\,\mathbf{v}_{k}$
$\displaystyle={\scriptstyle\Omega}_{k}^{2}\,\mathbf{v}_{k},\quad{\scriptstyle\Omega}_{k}\in\mathbb{R}$
$\displaystyle{\scriptstyle\Omega}_{k}$
$\displaystyle\neq{\scriptstyle\Omega}_{m}\implies\mathbf{v}_{k}^{*}\,\mathbf{v}_{m}=0,$

and the mutual orthogonality of the eigenvectors is due to the symmetry of the
Laplacians.

Also of interest to us is the Hamiltonian matrix $\mathcal{L}$ that acts as
the generator of system (3)

$\displaystyle\mathcal{L}:=\begin{bmatrix}0&I\\\ \text{-}L&0\end{bmatrix}$ (6)

The non-zero eigenvalues and eigenvectors of $\mathcal{L}$ are related to
those of $L$ by

$\displaystyle
L\,\mathbf{v}={\scriptstyle\Omega}^{2}\,\mathbf{v}\implies\mathcal{L}\begin{bmatrix}\mathbf{v}\\\
\pm j{\scriptstyle\Omega}\mathbf{v}\end{bmatrix}=\pm
j{\scriptstyle\Omega}\begin{bmatrix}\mathbf{v}\\\ \pm
j{\scriptstyle\Omega}\mathbf{v}\end{bmatrix}$ (7)

Since $L$ has non-negative real eigenvalues, all eigenvalues of $\mathcal{L}$
are on the imaginary axis. Further, each zero eigenvalue of $L$ corresponds to
a Jordan block of $\mathcal{L}$ of size $2$ with eigenvalue zero. For a
network with one connected component, the corresponding eigenvector is
$\mathbbm{1}$, the vector of $1$s.

As a result, the modes of Eq. 4 corresponding to the zero eigenvalue are
always unstable, even at $\epsilon=0$. However, these modes are physically
unimportant: they correspond to drift in the average value of the state $\phi$
and do not otherwise influence the stability of the system.

###  I-B Stability of the Mathieu Equation

Consider the Mathieu equation

$\ddot{x}(t)\,+\,\left(\omega_{\rm n}^{2}+\epsilon\,\cos(\omega
t)\rule{0.0pt}{8.99994pt}\right)x(t)~{}=~{}0,$ (8)

where $\omega_{\rm n}=\sqrt{k}$ is the unexcited (i.e. $\epsilon=0$) natural
frequency of the system. For any amplitude $\epsilon$ and frequency $\omega$
of the parametric excitation, Equation Eq. 8 is a linear periodic system (with
period $2\pi/\omega$), and therefore its stability properties can be
determined using the eigenvalues of its monodromy matrix, which in general
must be computed numerically. Figure 1 shows the regions of
stability/instability in the $(\epsilon,\omega)$ plane. The grey-shaded areas
are the unstable regions, while the remainder are the (neutrally) stable
regions.

The unstable regions “touch” the $\epsilon=0$ axis at discrete multiples of
$\omega$, $\tfrac{\omega}{\omega_{n}}=2,1,\tfrac{1}{2},\tfrac{1}{3},\ldots$.
These represent excitation frequencies for which the system (8) can be
destabilized with arbitrarily small excitation amplitudes. They are referred
to as “Arnold tongues”, and can also be found using perturbation methods. The
behavior of each Arnold tongue around $\epsilon\approx 0$ depends on its
index. The first tongue at $\omega=2\omega_{\rm n}$ expands linearly for small
epsilon, i.e. the stability boundary behaves like
$\omega(\epsilon)=2\omega_{\rm n}\pm a_{1}\epsilon+O(\epsilon^{2})$, while the
stability boundary at the second tongue behaves like
$\omega(\epsilon)=\omega_{\rm n}\pm a_{2}\epsilon^{2}+O(\epsilon^{3})$. More
generally, the “width” of the $i^{\text{th}}$ tongue behaves like
$\epsilon^{n}$, which makes each successive tongue “thinner”. This has
important implications for observing parametric resonance phenomena in systems
with small damping, as damping will effectively “lift” the tongues away from
the $\epsilon=0$ axis as shown in Fig. 1 (right). For the higher order tongues
(i.e. $i\geq 2$), this effect makes the threshold of $\epsilon$ where
instability is observed much higher than $\epsilon=0$, and consequently, only
the first tongue is observed in applications where there is any amount of
damping.

Figure 1:  The stability diagram for the undamped Mathieu equation (left), and
for the Mathieu equation with slight damping (right), where the gray-shaded
regions are the unstable regions in the $(\epsilon,\omega)$ plane of the
parametric excitation amplitude $\epsilon$ and frequency $\omega$. In the
undamped case, the system is unstable for arbitrarily small excitation
amplitude $\epsilon$ at excitation frequencies that are $2,1,1/2,1/3,...$
times the natural frequency $\omega_{\rm n}$. The shaded areas near
$\omega\approx\omega_{\rm n}/k,~{}k=1/2,1,2,3,4,..$ are called the “Arnold
tongues”, which get progressively thinner with increasing $k$. Thus with
slight damping, only the first tongue at $\omega\approx 2\omega_{\rm n}$
represents an instability phenomenon that can be triggered with small
excitation amplitude $\epsilon$. Only the first three tongues are shown in the
diagrams above.

The stability diagram for the more general Hill equation is qualitatively
similar to that of Fig. 1. The shapes of the instability regions in that case
are generally different, and depend on the details of the periodic function
$f$. However, the asymptotic behavior (around $\epsilon\approx 0$) of the
Arnold tongues is the same as shown in the figure.

##  II Full-network parametric forcing

We now consider the case where $P=L$ in model (4). This represents a scenario
where all the edges of the network are subjected to the same oscillatory
perturbation simultaneously. This case is the simplest to analyze since Eq. 4
becomes

$\displaystyle\ddot{\phi}+\left(\rule{0.0pt}{8.99994pt}L+\epsilon\,L\,\cos(\omega
t)\right)\phi$ $\displaystyle=0,$ $\displaystyle\Leftrightarrow\hskip
20.00003pt\ddot{\phi}+\left(\rule{0.0pt}{8.99994pt}1+\epsilon\cos(\omega
t)\right)L\,\phi$ $\displaystyle=0.$ (9)

Since $L$ has a full set of independent eigenvectors, this system is
diagonalizable using $L\,V=V\,\Omega^{2}$ where $V$ is the matrix whose
columns are the eigenvectors, and $\Omega^{2}$ is a diagonal matrix made up of
the corresponding eigenvalues $\left\\{{\scriptstyle\Omega}^{2}_{k}\right\\}$.
The result is a collection of uncoupled Mathieu’s equations:

$\displaystyle L\,V=V\,\Omega^{2}\implies$
$\displaystyle\left(V^{\text{-}1}\ddot{\phi}\right)+\left(\rule{0.0pt}{8.99994pt}1+\epsilon\cos(\omega
t)\right)\Omega^{2}\left(V^{\text{-}1}\phi\right)=0.$ (10)

Defining $p(t):=V^{\text{-}1}\ddot{\phi}(t)$, the $j^{\text{th}}$ such system
is

$\displaystyle\ddot{p}_{j}(t)+{\scriptstyle\Omega}_{j}^{2}\left(1+\epsilon\cos(\omega
t)\rule{0.0pt}{8.99994pt}\right)p_{j}(t)=0.$ (11)

This is the scalar Mathieu equation for an oscillator with unperturbed
“natural frequency” ${\scriptstyle\Omega}_{j}$ and parametric forcing of
amplitude ${\scriptstyle\Omega}_{j}^{2}\epsilon$ at frequency $\omega$. For
each $j$, the stability properties of the system (11) are the same as the
corresponding Mathieu equation (8) as shown in Fig. 1.

Figure 2:  Stability diagrams for networked second order oscillators
corresponding connected in a triangle (top) and square (bottom) layout with
specified nominal edge weights. The stability diagrams are overlapping copies
of scaled versions of Figure 1, one per mode of the system. The region of
instability is the union of the unstable regions for each mode.

The stable region in $(\omega,\epsilon)$ space for (9) is thus the
intersection of the stable regions for all the modes $p_{j},\ j=1,\dots n$.
Fig. 2 shows such stability diagrams for two examples with 3 and 4 node
networks. The basic pattern is as follows. For an $n$-node connected network,
the Laplacian has $n-1$ non-zero eigenvalues, and therefore by (11) there are
$n-1$ tongues of each type. More precisely, there are $n-1$ first order
tongues at excitation frequencies of

$\omega~{}=~{}2\,{\scriptstyle\Omega}_{j},\hskip 20.00003ptj=1,\ldots,n-1,$

with $\left\\{{\scriptstyle\Omega}_{j}^{2}\right\\}$ being the non-zero
eigenvalues of the Laplacian $L$. Similarly, there are $n-1$ second-order
tongues at excitation frequencies $\omega={\scriptstyle\Omega}_{j}$, and so
on.

Fig. 2 shows the first, second, and third-order tongues for two networks of 3
and 4 nodes respectively. For example, the 3-node network has two tongues of
each type whose locations are determined by the eigenvalues of the Laplacian.
As mentioned in Section I-B, only the first-order tongues are
phenomenologically relevant in systems with any amount of damping.

##  III Subnetwork Parametric Excitation

In this section we consider periodic perturbation of the weight of a portion
of the network. We begin with the case of single-edge forcing as the analysis
is very similar to the more general case of a subnetwork as shown in Section
III-D.

In the single-edge forcing case, $P=e_{ik}e_{ik}^{\star}=:E_{ik}$, the
Laplacian of the $n$-node graph with a single edge from node $i\to k$. Then
system (4) becomes

$\displaystyle\ddot{\phi}+\left(L+\epsilon\,E_{ik}\cos(\omega
t)\rule{0.0pt}{8.99994pt}\right)\phi=0$ (12)

Unlike the case of the full network forcing, we cannot treat this system as a
set of uncoupled second-order oscillators since the pair $L,E_{ik}$ will
typically not be simultaneously diagonalizable. Instead, we will find its
critical parametric forcing frequencies and stability behavior using a
perturbation expansion around $\epsilon=0$. This analysis is quite delicate,
and is presented in Section III-A. It relies on a vector-valued multiple scale
perturbation analysis that we develop in the appendices. This is a
generalizations of the standard scalar-valued multiple scale analysis [3, 14]
and we believe is one contribution of the present paper.

We now summarize the main stability result of this section before presenting
the perturbation technique. For the reasons outlined earlier, we only consider
the first-order instability tongues since they are the most phenomenologically
important.

Our main result is summarized next.

  * •

All first-order tongues of the stability diagram of (12) occur at the
following excitation frequencies

$\omega_{lm}={\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m},\hskip
20.00003pt{\scriptstyle\Omega}_{l},{\scriptstyle\Omega}_{m}\in\sqrt{{\rm
eigs}(L)-\\{0\\}},$ (13)

where ${\scriptstyle\Omega}_{l},{\scriptstyle\Omega}_{m}$ are any two non-zero
natural frequencies.

  * •

At each $\omega_{lm}$, the behavior of the stability boundary of that tongue
is

$\displaystyle\omega(\epsilon)$ $\displaystyle=\omega_{lm}\pm
a_{lm}\epsilon+O(\epsilon^{2})$ $\displaystyle\mbox{where}\hskip
20.00003pta_{lm}$
$\displaystyle=\frac{\left|\mathbf{v}_{m}^{*}\,e_{ik}e_{ik}^{\star}\,\mathbf{v_{l}}\right|}{2\left({\scriptstyle\Omega}_{l}\,{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(14)

The first statement (13) says that all first order tongues occur at sums of
the natural frequencies of the unperturbed system. Those in turn are given by
the square roots of the non-zero eigenvalues of the network Laplacian. Since
there are $n-1$ of those, then the number of first-order instability tongues
is $\left((n-1)^{2}+(n-1)\right)/2$ with possible repetitions if some natural
frequencies are repeated.

The quantity $a_{lm}$ (14) can be visualized as the “width” of the
$lm^{\text{th}}$ tongue as illustrated in Figure 3. This quantity can be used
as a proxy for the susceptibility of the system to instability when excited
near that frequency since it captures (up to first order) the width of the
excitation-frequency range that causes a parametric instability. We note that
this number depends not only on the eigenvalues
${\scriptstyle\Omega}_{l}^{2},{\scriptstyle\Omega}_{m}^{2}$ of $L$, but on the
corresponding eigenvectors $\mathbf{v}_{l},\mathbf{v}_{m}$ as well. Thus
first-order tongues in this case will typically come with a large variety of
widths, some of which may actually be zero.

Figure 3: Depiction of the first-order tongues of the single-edge excitation
model (12). All possible locations of such tongues are given by sums
$\omega_{lm}={\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}$ of natural
frequencies ${\scriptstyle\Omega}_{k}\in\sqrt{{\rm eigs}(L)}$ of the
unperturbed system, which in turn are determined by the eigenvalues of the
network Laplacian $L$. The “widths” $a_{lm}$ of each tongue are determined by
both eigenvalues and eigenvectors of $L$ by (14).

###  III-A Perturbation Analysis

We assume that $\phi(t;\epsilon)$ is analytic in $\epsilon$, so we may write

$\displaystyle\phi(t;\epsilon)$
$\displaystyle=:\phi_{0}(t)+\epsilon\phi_{1}(t)+\epsilon^{2}\phi_{2}(t)+\mathcal{O}(\epsilon^{3}),$
(15)

where $\phi_{k}(t)$ are differentiable functions of time.

Applying this expansion to system (12) and collecting powers of $\epsilon$, we
get the triangular system

$\displaystyle\mathcal{O}(\epsilon^{0}):\ $
$\displaystyle\ddot{\phi}_{0}+L\phi_{0}=0$ (16)
$\displaystyle\mathcal{O}(\epsilon^{1}):\ $
$\displaystyle\ddot{\phi}_{1}+L\phi_{1}=\text{-}\cos(\omega t)E_{ik}\phi_{0}$
(17) $\displaystyle\vdots$ $\displaystyle\mathcal{O}(\epsilon^{n}):\ $
$\displaystyle\ddot{\phi}_{n}+L\phi_{n}=\text{-}\cos(\omega
t)E_{ik}\phi_{n\text{-}1}.$ (18)

$\phi_{k}$ at each order of $\epsilon$ behaves as a coupled harmonic
oscillator with additive input $\phi_{k\text{-}1}$. $\phi_{0}(t)$ is the
solution to the unperturbed, linearized oscillator dynamics (3):

$\displaystyle\phi_{0}(t)=e^{j\sqrt{L}t}(V\mathbf{c})+e^{\text{-}j\sqrt{L}t}(V\mathbf{c})^{\star}$
(19)

for some $\mathbf{c}\in\mathbb{C}^{n}$.

To analyze the $\mathcal{O}(\epsilon^{1})$ term, it will be useful to work in
the basis of the eigenvectors of $L$. With $LV=V\Omega^{2}$ (Section I-A), we
can write

$\displaystyle V^{\text{-}1}\phi_{0}(t)=e^{j\Omega
t}\mathbf{c}+e^{\text{-}j\Omega t}\mathbf{c}^{\star}$ (20)

for appropriate initial conditions $\mathbf{c}$. The
$\mathcal{O}(\epsilon^{1})$ system is

$\displaystyle\left(V^{\text{-}1}\ddot{\phi}_{1}\right)+\Omega^{2}\left(V^{\text{-}1}\phi_{1}\right)=-V^{\text{-}1}E_{ik}V\cos(\omega
t)\left(V^{\text{-}1}\phi_{0}\right)$
$\displaystyle=-\tfrac{1}{2}V^{\text{-}1}E_{ik}V\cos(\omega t)\left(e^{j\Omega
t}\mathbf{c}+e^{\text{-}j\Omega t}\mathbf{c}^{\star}\right)$
$\displaystyle=-\tfrac{1}{2}V^{\text{-}1}E_{ik}V\left(e^{j(\Omega\pm\omega
I)t}\mathbf{c}+e^{\text{-}j(\Omega\pm\omega I)t}\mathbf{c}^{\star}\right)$
(21)

Each forcing term in system (21) is thus of the form

$\displaystyle-\tfrac{1}{2}\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\left(e^{j(\Omega_{l}\pm\omega)t}c_{l}+e^{j(\Omega_{l}\pm\omega)t}c^{\star}_{l}\right)$
(22)

where $l,m=1,2,\dots,n$. System (21) is neutrally stable when no natural
frequency $\pm j{\scriptstyle\Omega}$ equals any frequency of the forcing $\pm
j{\scriptstyle\Omega}\pm j\omega$. Indeed, if any of the forcing frequencies
equals one of these modal frequencies, so that
$j{\scriptstyle\Omega}_{m}=j{\scriptstyle\Omega}_{l}+j\omega$, say, the
corresponding mode of the solution experiences harmonic resonance and is
unstable. In the parlance of perturbation theory, the corresponding forcing
term is _secular_. Additional care has to be taken in this line of reasoning
when the eigenvalues of $L$ have multiplicity $>1$. This case is handled in
Section A-B.

Consequently the $\mathcal{O}(\epsilon)$ term of the solution (15) is unstable
when

$\displaystyle\pm{\scriptstyle\Omega}_{m}=\pm{\scriptstyle\Omega}_{l}\pm\omega$
(23)

for any $l,m=1,2,\dots,n$.

By assumption, $\omega>0$, and all the eigenvalues ${\scriptstyle\Omega}$ are
non-negative. So this list of conditions reduces to the following two:

$\displaystyle\omega$
$\displaystyle=\left|{\scriptstyle\Omega}_{m}+{\scriptstyle\Omega}_{l}\right|$
(24) $\displaystyle\omega$
$\displaystyle=\left|{\scriptstyle\Omega}_{m}-{\scriptstyle\Omega}_{l}\right|$
(25) $\displaystyle l,m$ $\displaystyle=1,2,\dots,n$

The critical parametric forcing frequency, if it exists, must be the sum or
difference of the square roots of the eigenvalues of the graph Laplacian. This
condition is subject to the validity of the regular perturbation analysis, and
also requires that the coefficient of the relevant additive forcing term
($\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}$) corresponding to the
eigenvalues ${\scriptstyle\Omega}_{l}$, ${\scriptstyle\Omega}_{m}$ of $L$ be
non-zero (Eq. 22). This requirement is investigated in more detail in Section
III-C.

The boundedness of the $\mathcal{O}(\epsilon^{2})$ and higher order terms in
the ansatz (15) can be investigated similarly. Assuming $\phi_{1}(t)$ is
bounded, it has the form

$\displaystyle V^{\text{-}1}\phi_{1}(t)=$ $\displaystyle e^{j\Omega
t}\mathbf{c}_{0}+e^{j(\Omega\pm\omega
I)t}\mathbf{c}_{1}+e^{\text{-}j(\Omega\pm\omega I)t}\mathbf{c}_{2}$
$\displaystyle+\text{complex conjugates,}$ (26)

where $\mathbf{c}_{*}$ are constants that depend on the initial conditions.
Then $\phi_{2}(t)$ is the solution of the $\mathcal{O}(\epsilon^{2})$ system

$\displaystyle
V^{\text{-}1}\ddot{\phi}_{2}+\Omega^{2}\,V^{\text{-}1}\phi_{2}=-\tfrac{V^{\text{-}1}E_{ik}V}{2}\left(e^{j\omega
t}+e^{-j\omega t}\right)V^{\text{-}1}\phi_{1}$ (27)

which contains additive forcing at frequencies
${\scriptstyle\Omega}\pm\omega$, as before, as well as at
${\scriptstyle\Omega}\pm 2\omega$. Thus the $\phi_{2}(t)$ system experiences
harmonic forcing whenever any of these frequencies equals a natural frequency
of the system, or

$\displaystyle\pm j{\scriptstyle\Omega}_{m}$ $\displaystyle=\pm
j{\scriptstyle\Omega}_{l}\pm j\omega,\text{ or}$ $\displaystyle\pm
j{\scriptstyle\Omega}_{m}$ $\displaystyle=\pm j{\scriptstyle\Omega}_{l}\pm
2j\omega$

for any $l,m=1,2,\dots,n$. The first of these conditions does not hold since,
by assumption, $\phi_{1}(t)$ is stable. As before (Eq. 24 and Eq. 25), this
list of conditions reduces to
$\omega=\left|{\scriptstyle\Omega}_{m}\pm{\scriptstyle\Omega}_{l}\right|/2$.
By a similar argument, we can show that the $\mathcal{O}(\epsilon^{n})$ term
$\theta_{n}(t)$ is unstable when

$\displaystyle\omega$
$\displaystyle=\tfrac{1}{n}\left|{\scriptstyle\Omega}_{m}+{\scriptstyle\Omega}_{l}\right|$
(28) $\displaystyle\omega$
$\displaystyle=\tfrac{1}{n}\left|{\scriptstyle\Omega}_{m}-{\scriptstyle\Omega}_{l}\right|$
(29)

Fig. 4 compares (to $\mathcal{O}(\epsilon)$) the numerically computed
stability diagram for networked oscillators with our predictions for the
critical parametric forcing frequencies. We observe that:

  * •

While all the critical frequencies (as $\epsilon\to 0$) to first order are
accounted for, criteria (24) and (25) generate “false positives” as well.

  * •

The regular perturbation analysis gives us no information on the slopes or
widths of the tongues. Thus we cannot order by severity (or excitability of
parametric resonance) either the critical frequencies or choices of the edge
to be perturbed.

To refine the instability conditions, we apply a multiple time-scale
perturbation analysis to our original system (12). This allows us to
investigate the behavior of the system in the vicinity of the critical
frequencies in $(\omega,\epsilon)$ parameter space, resolve the “false
positives” and obtain a measure of the instability at each critical frequency.

Figure 4:  The numerically computed stability diagram for a triangular graph
with unequal edge weights. The overlaid red lines are the critical parametric
forcing frequencies predicted (to first order) by the regular perturbation
analysis of Section III-A. The narrow tongues at $\omega\approx 2.9$ and
$\omega\approx 4.7$ correspond to the predictions from higher order
perturbations (Eq. 28 and Eq. 29). We observe that while the Arnold tongue
positions are predicted correctly by the perturbation analysis, the prediction
of $2.85$ is a false positive. The perturbation analysis is refined using a
multiple time-scale method in Section A-A.

We do this in the following sequence:

  1. 1.

Investigate the stability of the system along straight lines in the
$(\omega,\epsilon)$ parameter space through the critical frequencies. Each
line is determined by the $\omega$-intercept $\omega_{0}$ and the slope
$\left.\tfrac{\mathrm{d}\omega}{\mathrm{d}\epsilon}\right\rvert_{\epsilon=0}=:a$.

  2. 2.

Expand $\phi(t)$ along the fast and slow time scales $t$ and $\epsilon t$.
This decomposition is non-unique.

  3. 3.

System stability requires both the fast and slow time dynamics to be stable.
Enforce this condition to give us both (i) a suitable decomposition of the
dynamics and (ii) stability criteria for the $\omega$ intercept and $a$.

The details of the multiple scale perturbation analysis are presented in
Section A-A. We summarize the results here:

######  Theorem 1 (Instability criteria)

Consider the curve in the $(\omega,\epsilon)$ parameter plane given by
$\omega(\epsilon)=\omega_{0}+a\epsilon+\mathcal{O}(\epsilon^{2})$. Let
${\scriptstyle\Omega}_{m}^{2}$ and $\mathbf{v}_{m}$ be the eigenvalues and
eigenvectors of the graph Laplacian $L$, for $m=1,2,\dots,n$. In the limit
$\epsilon\to 0$, system (12) is unstable at $(\omega(\epsilon),\epsilon)$ if

$\displaystyle\omega_{0}={\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m},\
l,m\in 1,2,\dots,n$
$\displaystyle\left|a\right|<\frac{\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}.$
(30)

We interpret the result as follows:

  * •

The stability criteria obtained from the regular perturbation analysis suggest
that both $\pm\omega_{0}={\scriptstyle\Omega}_{l}-{\scriptstyle\Omega}_{m}$
and $\omega_{0}={\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}$ cause
instability. However the method of multiple time-scales rules out the former.
This is the first source of false positives in the results suggested by the
regular perturbation analysis.

  * •

$a=0$ corresponds to a vertical line in the $(\omega,\epsilon)$ parametric
plot. This implies that to $\mathcal{O}(\epsilon)$, the stability boundary has
a symmetric V-shape, where the critical angle from the vertical to either side
is
$\pm\tan^{-1}\left(\tfrac{\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}\right)$.

  * •

If $\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}=0$, then the system is stable
around $(\omega_{0},0)$. This is the second source of false positives in the
regular perturbation analysis. It is also consistent with the result of the
controllability analysis (Section III-C), which posits the same criterion for
instability.

###  III-B Comparison with Numerical Calculations

Fig. 5 compares the approximation (30) to the results of numerical computation
for the stability of a test graph. We observe that the largest Arnold tongue
structures correspond to the critical values of $\omega$ predicted by the
$\mathcal{O}(\epsilon)$ term in the perturbation analysis. The slopes of the
tongues as $\epsilon\to 0$ are captured accurately by the multiple time-scale
analysis above.

Figure 5: A comparison of the numerically computed stability diagram with the
stability estimates from perturbation analysis for swing dynamics on
triangular graph with one time-periodically perturbed edge. The wide Arnold
tongue positions in $(\omega,\epsilon)$ space (black curves) correspond to
first order perturbation analysis and are predicted correctly (dashed red
lines). The multiple time-scale analysis from section A-A eliminates the false
positives observed in figure 4. Additionally, the estimate of the slope of the
Arnold tongues as $\epsilon\to 0$ computed using (81) (solid red lines) match
the numerical results closely, and thus provide a simple way of predicting the
degree of instability of each critical parametric forcing frequency $\omega$
from the properties of the graph Laplacian. The narrower, partially resolved
Arnold tongues correspond to higher order terms in the perturbation expansion
(15).

###  III-C Controllability Interpretation

We can understand the stability criterion (30) using a controllability
argument. The oscillator dynamics with single edge forcing are given by Eq.
12, or

$\displaystyle\ddot{\phi}+\left(L+\epsilon E_{ik}\cos(\omega t)\right)\phi=0.$

Following the perturbation expansion (15) and solution (60), we have, up to
$\mathcal{O}(\epsilon)$

$\displaystyle\phi_{0}(t)$ $\displaystyle=V\left(e^{j\Omega
t}\mathbf{c}+e^{\text{-}j\Omega t}\mathbf{c}^{\star}\right)$

where $\phi_{0}(t)$ is the solution to the unperturbed, linearized swing
equations, and appears as an additive input in the equation for the
$\mathcal{O}(\epsilon)$ term $\phi_{1}(t)$:

$\displaystyle\ddot{\phi}_{1}(t)$
$\displaystyle=-L\phi_{1}(t)-E_{ik}\cos(\omega t)\phi_{0}(t)$
$\displaystyle=-L\phi_{1}(t)-\underbrace{E_{ik}V}_{\mathcal{B}}\underbrace{e^{j(\Omega+\omega
I)t}\mathbf{c}}_{u_{1}(t)}-\underbrace{E_{ik}V}_{\mathcal{B}}\underbrace{e^{\text{-}j(\Omega\text{-}\omega
I)t}\mathbf{c}}_{u_{2}(t)}$
$\displaystyle\quad-\mathcal{B}u_{1}^{\star}(t)-\mathcal{B}u_{2}^{\star}(t).$

Here, $\mathcal{B}:=E_{ik}V$. The inputs $u_{1}(t)$ and $u_{2}(t)$ are
composed of sinusoids at the frequencies $(\pm\Omega_{l}\pm\omega)$, which may
destabilize the system if (following conditions (24) and (25)) this equals
$\pm{\scriptstyle\Omega}_{m}$, one of the modes of the unforced system.
However, this also requires the system to be controllable from those inputs.
For instance, when
$j{\scriptstyle\Omega}_{m}=j\left({\scriptstyle\Omega}_{l}+\omega\right)$ for
some $l,m$, we have

$\displaystyle\ddot{\phi}_{1}=$
$\displaystyle-L\phi_{1}(t)-\tfrac{1}{2}\underbrace{E_{ik}\mathbf{v}_{l}}_{\mathcal{B}_{l}}\underbrace{e^{j(\Omega_{l}+\omega)t}c_{l}}_{u_{l}(t)}$
$\displaystyle+\text{other inputs}$ (31)

where the other inputs include the complex conjugate of $u_{l}(t)$, as well as
inputs that do not destabilize the system. Applying the eigenvector
controllability test to system (31) rewritten in first order form (Section
A-D) gives us the following criterion for uncontrollability:

$\displaystyle\begin{bmatrix}\mathbf{v}_{m}^{\star}&\text{-}j{\scriptstyle\Omega}_{m}^{\text{-}1}\mathbf{v}_{m}^{\star}\end{bmatrix}\begin{bmatrix}0&0\\\
E_{ik}&0\end{bmatrix}\begin{bmatrix}\mathbf{v}_{l}\\\
j{\scriptstyle\Omega}_{l}\mathbf{v}_{l}\end{bmatrix}$
$\displaystyle\stackrel{{\scriptstyle?}}{{=}}0$
$\displaystyle\implies\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}$
$\displaystyle\stackrel{{\scriptstyle?}}{{=}}0.$ (32)

When this is the case, the mode with frequency $\pm j{\scriptstyle\Omega}_{m}$
is not controllable from $u_{l}(t)$, and thus the harmonic resonance predicted
by the instability criteria (24) and (25) cannot occur. The controllability-
based criterion (32) thus provides an alternative interpretation of the
criterion obtained by the multiple-scale perturbation analysis (30).

###  III-D Generalization to Subnetwork Excitation

The methods described in sections II, III-A and Section A-A can be extended to
the parametric forcing of any subnetwork $P$ of the graph $L$, where $P$ and
$L$ are as in Eq. 4. The critical frequencies are still subject to the same
conditions (28) and (29), so they are unchanged. The perturbation analysis is
unchanged but for replacing the single edge Laplacian $E_{ik}$ with $P$ in
(30). The slopes of the first-order Arnold tongues are thus given by

$\displaystyle\left|a_{lm}\right|=$
$\displaystyle\frac{\left|\mathbf{v}_{m}^{\star}P\mathbf{v}_{l}\right|}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(33)

When the subnetwork $P$ in (4) is the full network $L$, we expect the
perturbation analysis criterion (30) to reduce to the results of Section II.
When $P=L$, all cross terms of the form
$\mathbf{v}_{m}^{\star}P\mathbf{v}_{l}=\mathbf{v}_{m}^{\star}L\mathbf{v}_{l}$
with $m\neq l$ vanish, since
$\mathbf{v}_{m}^{\star}L\mathbf{v}_{l}={\scriptstyle\Omega}_{l}^{2}\
\mathbf{v}_{m}^{\star}\mathbf{v}_{l}=0$. Thus the only critical frequencies
are of the form
$\omega_{0}=\left|{\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{l}\right|/n=2{\scriptstyle\Omega}_{l}/n$,
for $n=1,2,\dots,{\scriptstyle N}$. These are the parametric resonance
frequencies of ${\scriptstyle N}$ decoupled Mathieu oscillators, which agrees
with the diagonalization argument of Section II.

##  IV Stability Analysis for Tori Graphs

As an example, we apply the stability analysis results of Section III to ring
and tori graphs. We can obtain analytical expressions or bounds for all the
frequencies and tongue widths for parametric resonance on such networks, which
gives us some insight into how the stability characteristics vary with
increasing graph size and connectivity.

We begin with the one-dimensional torus, i.e. a ring graph. With
$\theta:=\frac{2\pi}{N}$, the eigenvalues and eigenvectors of the
${\scriptstyle N}$-node ring graph are [15]

$\displaystyle{\scriptstyle\Omega}_{m}^{2}$
$\displaystyle=2(1-\cos(m\theta)),\quad m=0,\dots,{\scriptstyle N}\text{-}1$
(34) $\displaystyle\mathbf{v}_{m}$
$\displaystyle=\begin{bmatrix}1&e^{jm\theta}&e^{j2m\theta}&\dots&e^{j(N\text{-}1)m\theta}\end{bmatrix}^{\mathrm{T}}$
(35)

For the ring graph, we can thus find the critical frequencies (24) and (25)
and the slopes of the corresponding Arnold tongues (30) analytically. The
frequency corresponding to the eigenvalue pair $(l,m)$ is

$\displaystyle\omega_{lm}=$
$\displaystyle{\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m},\quad(l,m=1,\dots,{\scriptstyle
N}\text{-}1)$ $\displaystyle=$ $\displaystyle
2^{1/2}\left((1-\cos(l\theta))^{1/2}+(1-\cos(m\theta))^{1/2}\right)$ (36)

Since ${\scriptstyle\Omega}_{m}^{2}={\scriptstyle\Omega}_{N\text{-}m}^{2}$,
all eigenvalues of the Laplacian (except possibly one) have multiplicity $2$.
The corresponding eigenvectors are the complex conjugates $\mathbf{v}_{m}$ and
$\mathbf{v}_{N\text{-}m}=\mathbf{\bar{v}}_{m}$.

The first-order Arnold tongues in $\omega,\epsilon$ space are then given by

$\displaystyle\omega(\epsilon)=\omega_{lm}\pm
a_{lm}\epsilon+\mathcal{O}(\epsilon^{2})$ (37)

where $\pm a_{lm}$ is the slope of the tongue at critical frequency
$\omega_{lm}$. $a_{lm}$ is the maximum attainable value of
${\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}}/{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(Section A-B), where $E_{ik}=e_{ik}e_{ik}^{\star}$ is the single edge ($i\to
k$) Laplacian corresponding to the parametrically forced edge, and
$\mathbf{v}_{m}$ and $\mathbf{v}_{l}$ range over the corresponding
$2$-dimensional eigenspaces. If $m\neq{\scriptstyle N}/2$ and
$l\neq{\scriptstyle N}/2$,

$\displaystyle\left|a_{lm}\right|=$
$\displaystyle\frac{\left\|\begin{pmatrix}\mathbf{v}_{m}^{\star}\\\
\mathbf{v}_{N\text{-}m}^{\star}\end{pmatrix}E_{ik}\begin{pmatrix}\mathbf{v}_{l}&\mathbf{v}_{N\text{-}l}\end{pmatrix}\right\|_{2}}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
$\displaystyle=$
$\displaystyle\frac{\left\|\begin{matrix}\mathbf{v}_{m}^{\star}e_{ik}\\\
\mathbf{\bar{v}}_{m}^{\star}e_{ik}\end{matrix}\right\|_{2}\left\|\begin{matrix}e_{ik}^{\star}\mathbf{v}_{l}&e_{ik}^{\star}\mathbf{\bar{v}}_{l}\end{matrix}\right\|_{2}}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
$\displaystyle=$
$\displaystyle\frac{\left|\mathbf{v}_{m}^{\star}e_{ik}\right|}{{\scriptstyle\Omega}_{l}^{1/2}}\frac{\left|e_{ik}^{\star}\mathbf{v}_{l}\right|}{{\scriptstyle\Omega}_{m}^{1/2}}=\frac{\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|}{\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(38) $\displaystyle=$ $\displaystyle\frac{1}{{\scriptstyle
N}}\frac{\left|1-e^{\text{-}jm\theta}\right|}{(2\text{-}2\cos(\text{-}m\theta))^{1/4}}\frac{\left|1-e^{jl\theta}\right|}{(2\text{-}2\cos(l\theta))^{1/4}}$
(39)

If $m={\scriptstyle N}/2$ or $l={\scriptstyle N}/2$, then the corresponding
eigenspace is one-dimensional, and the slope $\left|a_{lm}\right|$ is lower by
a factor of $2^{1/2}$:

$\displaystyle\left|a_{lm}\right|$
$\displaystyle=\frac{\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|}{\nu\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}},\text{
with}$ (40) $\displaystyle\nu$
$\displaystyle=\begin{cases}1&l\neq{\scriptstyle N}/2,m\neq{\scriptstyle
N}/2\\\ 2^{1/2}&l={\scriptstyle N}/2\text{ or }m={\scriptstyle N}/2\\\
2&l=m={\scriptstyle N}/2\end{cases}$

Expression (39) is separable and symmetric in $m$ and $l$, so we can find the
most critical frequencies by maximizing the function
$\delta(t):=\left|1-e^{jt}\right|(2-2\cos(t))^{\text{-}1/4}$, where
$\left|a_{lm}\right|=\frac{1}{\nu{\scriptstyle
N}}\delta(\text{-}m\theta)\delta(l\theta)$. This function attains a maximum of
$2^{1/2}$ at $t=\pi$, corresponding to $l=m=\lfloor{\scriptstyle N}/2\rfloor$.
We make the following observations about the stability of the ring graph under
parametric excitation of one of its edges:

  * •

All critical frequencies for parametric instability at small forcing
amplitudes are contained in $\omega\in(0,4]$ (Eq. 34). The critical
frequencies are independent of the choice of parametrically forced edge(s).

  * •

The size of the unstable region in $(\omega,\epsilon)$ space around any
critical frequency is inversely proportional to the number of nodes
${\scriptstyle N}$ in the graph.

  * •

As measured from the vertical line, the widest Arnold tongue in the parameter
space $(\omega,\epsilon)$ has slope between ${\scriptstyle 1/N}$ and
$2/{\scriptstyle N}$.

Fig. 6 shows the $\mathcal{O}(\epsilon)$ approximation to the Arnold tongues
(Eq. 37) for ring networks of varying sizes.

Figure 6: Visualization of the stability diagram for ring networks with single
edge perturbation. The four networks considered have (top-left) ${\scriptstyle
N}=4$, (top-right) $8$, (bottom-left) $16$ and $24$ nodes. All parametric
forcing frequencies causing parametric resonance are contained in the interval
$(0,4]$. The most critical frequency corresponds to the widest Arnold tongue
and is colored purple. For sufficiently large ${\scriptstyle N}$, this tongue
is at $\omega=4$ or close to it, and the slopes of all tongues as measured
from the vertical are inversely proportional to ${\scriptstyle N}$, the number
of nodes.

###  IV-A Ring Networks with more General Forcing

These analytical results for ring networks can be extended to the case of more
general forcing, when more than one edge is forced parametrically at a single
frequency $\omega$. Let $P$ be the graph Laplacian corresponding to the
parametrically forced edges, so that the full graph Laplacian is $L+\epsilon
P\cos(\omega t)$ (Section I). With some abuse of notation, we refer to the
graph corresponding to $P$ as $P$ as well. The perturbation analysis is
unchanged but for replacing the single edge Laplacian $E_{ik}$ with $P$. As
before, with ${\scriptstyle\Omega}_{m}^{2}={\scriptstyle\Omega}_{l}^{2}$ and
$\mathbf{v}_{m}=\mathbf{\bar{v}}_{N\text{-}m}$the slopes of the Arnold tongues
are given by (Section III-D)

$\displaystyle\left|a_{lm}\right|=$
$\displaystyle\frac{\left\|\begin{pmatrix}\mathbf{v}_{m}^{\star}\\\
\mathbf{v}_{{\scriptscriptstyle
N}\text{-}m}^{\star}\end{pmatrix}P\begin{pmatrix}\mathbf{v}_{l}&\mathbf{v}_{N\text{-}l}\end{pmatrix}\right\|_{2}}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
$\displaystyle=$
$\displaystyle\frac{1}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}\left\|\begin{matrix}\mathbf{v}_{m}^{\star}P\mathbf{v}_{l}&\mathbf{v}_{m}^{\star}P\mathbf{\bar{v}}_{l}\\\
\mathbf{\bar{v}}_{m}^{\star}P\mathbf{v}_{l}&\mathbf{\bar{v}}_{m}^{\star}P\mathbf{\bar{v}}_{l}\end{matrix}\right\|_{2}$
$\displaystyle=$
$\displaystyle\frac{2^{1/2}}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}\left\|\begin{matrix}\mathbf{v}_{m}^{\star}P\mathbf{v}_{l}\\\
\mathbf{\bar{v}}_{m}^{\star}P\mathbf{v}_{l}\end{matrix}\right\|_{2}$
$\displaystyle=$
$\displaystyle\frac{\left(\left|\mathbf{v}_{m}^{\star}P\mathbf{v}_{l}\right|^{2}+\left|\mathbf{\bar{v}}_{m}^{\star}P\mathbf{v}_{l}\right|^{2}\right)^{1/2}}{\left(2{\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(41)

Accounting for one-dimensional eigenspaces corresponding to $m={\scriptstyle
N}/2$ or $l={\scriptstyle N}/2$, we get

$\displaystyle\left|a_{lm}\right|$
$\displaystyle=\frac{\left(\left|\mathbf{v}_{m}^{\star}P\mathbf{v}_{l}\right|^{2}+\left|\mathbf{\bar{v}}_{m}^{\star}P\mathbf{v}_{l}\right|^{2}\right)^{1/2}}{\nu\left(2{\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(42)

where $\nu$ is the same as in Eq. 40.

Since $P$ can be written as the sum of single edge graph Laplacians
($P=\sum\limits_{i\to k}E_{ik}$, where $i\to k$ are the edges of $P$ with
$\left|i-k\right|=1$), we can simplify each term in this expression further:

$\displaystyle\left|\mathbf{v}_{m}^{\star}P\mathbf{v}_{l}\right|=$
$\displaystyle\left|\sum\limits_{\begin{subarray}{c}i\to k\\\
\left|i\text{-}k\right|=1\end{subarray}}\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|$
$\displaystyle=$ $\displaystyle\frac{1}{{\scriptstyle
N}}\left|\left(1-e^{j\theta l}\right)\left(1-e^{\text{-}j\theta
m}\right)\sum\limits_{i\to k}e^{j\theta i(l\text{-}m)}\right|,$

with a similar expression for
$\left|\mathbf{\bar{v}}_{m}^{\star}P\mathbf{v}_{l}\right|$. Both expressions
attain a maximum of $4{\scriptstyle M/N}$ at $l=m=\left\lfloor
N/2\right\rfloor$, allowing us to draw the following conclusions.

  * •

As in the single edge case, all critical frequencies for parametric
instability at small forcing amplitudes are contained in $\omega\in(0,4]$.

  * •

The size of the unstable region in $(\omega,\epsilon)$ space near any critical
frequency is inversely proportional to the number of nodes ${\scriptstyle N}$,
and the the widest tongue has slope (as measured from the vertical) between
${\scriptstyle M/N}$ and $2{\scriptstyle M/N}$.

When the whole network is perturbed (${\scriptstyle M=N}$), the analysis of
Section II applies, and the slope of the widest Arnold tongue is $\pm 1$, a
well known result for Mathieu’s equation.

###  IV-B Higher-Dimensional Tori

We can study the effect of the connectivity of the ring network on its
stability under parametric edge excitation by extending the results (39) and
(42) to tori in higher dimensions. We consider a ring network that is a
$d$-fold periodic lattice with ${\scriptstyle N}^{d}$ nodes, where each node
has degree $d$. The Laplacian of this graph is an ${\scriptstyle
N}^{d}\times{\scriptstyle N}^{d}$ matrix, whose spectral properties govern the
slopes of the Arnold tongues as in (30). Since this matrix is $d$-fold
circulant [15], its action on a vector $\mathbf{v}$ is equivalent to a
$d$-fold circular convolution:

$\displaystyle L_{a}\mathbf{v}$ $\displaystyle=\mathbf{a}\star\mathbf{v}$
$\displaystyle\left(L_{a}\mathbf{v}\right)_{(p_{1},\dots,p_{d})}$
$\displaystyle=\sum\limits_{\begin{subarray}{c}q_{1},\dots,q_{d}\in\mathbb{Z}_{N}\end{subarray}}\mathbf{a}_{(p_{1}-q_{1}),\dots,(p_{d}-q_{d})}\mathbf{v}_{(q_{1},\dots,q_{d})}$
$\displaystyle\implies\left(L_{a}\mathbf{v}\right)_{\mathbf{p}}$
$\displaystyle=\sum\limits_{\mathbf{q}\in\mathbb{Z}_{N}^{d}}\mathbf{a}_{\mathbf{p}-\mathbf{q}}\mathbf{v}_{\mathbf{q}}$

where $L_{a}$ is the Laplacian and $\mathbf{a}$ is one of its columns that
fully determines $L_{a}$. We index $\mathbf{a}$ and $\mathbf{v}$ as
multidimensional arrays with $d$ indices
$\mathbf{p}:=\left(p_{1},\dots,p_{d}\right)$ and
$\mathbf{q}:=\left(q_{1},\dots,q_{n}\right)$ in $\mathbb{Z}_{N}^{d}$, which
makes the convolution (and thus the spectrum of $L_{a}$) convenient to
represent. The eigenvectors of $L_{a}$ are the eigenvectors of the $d$-fold
circular shift operator, and the eigenvalues are the entries of the
$d$-dimensional DFT of $\mathbf{a}$. Indexing the eigenvalues by
$(m_{1},\dots,m_{d})=:\mathbf{m}\in\mathbb{Z}_{N}^{d}$, we have

$\displaystyle{\scriptstyle\Omega}^{2}_{\mathbf{m}}$
$\displaystyle=2\sum\limits_{i=1}^{d}\left(1\text{-}\cos(m_{i}\theta)\right))$
(43) $\displaystyle\mathbf{v}^{(\mathbf{m})}_{(q_{1},q_{2},\dots,q_{d})}$
$\displaystyle=e^{j\theta\left(m_{1}q_{1}+\dots+m_{d}q_{d}\right)}$
$\displaystyle\implies\mathbf{v}^{(\mathbf{m})}_{\mathbf{q}}$
$\displaystyle=e^{j\theta\mathbf{m}\cdot\mathbf{q}}$ (44)

where $\mathbf{v}^{(\mathbf{m})}$ is the eigenvector corresponding to
${\scriptstyle\Omega}_{\mathbf{m}}^{2}$ and $\theta=2\pi/{\scriptstyle N}$, as
before. The critical frequencies (24) and (25) are then

$\displaystyle\omega_{\mathbf{l}\mathbf{m}}={\scriptstyle\Omega}_{\mathbf{l}}+{\scriptstyle\Omega}_{\mathbf{m}}$
$\displaystyle=\sqrt{2}\left[\left(\sum\limits_{i=1}^{d}1\text{-}\cos(l_{i}\theta)\right)^{\frac{1}{2}}+\left(\sum\limits_{i=1}^{d}1\text{-}\cos(m_{i}\theta)\right)^{\frac{1}{2}}\right]$
(45)

Because of the $d$-fold symmetry in the expression for
${\scriptstyle\Omega}_{\mathbf{m}}^{2}$
($\cos(m_{i}\theta)=\cos(({\scriptstyle N}\text{-}m_{i})\theta)$), the
eigenvalue ${\scriptstyle\Omega}_{\mathbf{m}}^{2}$ has multiplicity $2^{d}$,
and the slope of the corresponding Arnold tongue involves the $2$-norm of a
matrix (Section A-B). Let $V_{\mathbf{m}}$ and $V_{\mathbf{l}}$ be matrices
whose columns span the eigenspace corresponding to
${\scriptstyle\Omega}_{\mathbf{m}}^{2}$ and
${\scriptstyle\Omega}_{\mathbf{l}}^{2}$ respectively, and
$E_{\mathbf{ik}}\equiv e_{\mathbf{ik}}e_{\mathbf{ik}}^{\star}$ be the single
edge Laplacian corresponding to the edge
$\mathbf{i}\leftrightarrow\mathbf{k}$, where
$\mathbf{i},\mathbf{k}\in\mathbb{Z}_{N}^{d}$. Note that by the topology of the
ring lattice, $\mathbf{i}$ and $\mathbf{k}$ are identical in all indices
except for one, which differs by one. Without loss of generality, we may
assume this is the first index. The slopes of the Arnold tongues (Eq. 30)
corresponding to perturbations in edge $\mathbf{i}\to\mathbf{k}$ are then

$\displaystyle\left|a_{\mathbf{l}\mathbf{m}}\right|=$
$\displaystyle\max_{\left|\mathbf{z}^{(\mathbf{m})}\right|=\left|\mathbf{z}^{(\mathbf{l})}\right|=1}\frac{\mathbf{z}^{(\mathbf{m})\star}V_{\mathbf{m}}^{\star}E_{ik}V_{\mathbf{l}}\mathbf{z}^{(\mathbf{l})}}{2\left({\scriptstyle\Omega}_{\mathbf{l}}{\scriptstyle\Omega}_{\mathbf{m}}\right)^{1/2}}$
$\displaystyle=$
$\displaystyle\frac{\|V_{\mathbf{m}}^{\star}E_{ik}V_{\mathbf{l}}\|_{2}}{2\left({\scriptstyle\Omega}_{\mathbf{l}}{\scriptstyle\Omega}_{\mathbf{m}}\right)^{1/2}}$
(46) $\displaystyle\geq$
$\displaystyle\frac{\left|\mathbf{v}^{(\mathbf{m})\star}E_{\mathbf{ik}}\mathbf{v}^{(\mathbf{l})}\right|}{2\left({\scriptstyle\Omega}_{\mathbf{l}}{\scriptstyle\Omega}_{\mathbf{m}}\right)^{1/2}}=\frac{\left|(\mathbf{v}^{(\mathbf{m})\star}e_{\mathbf{ik}})(e_{\mathbf{ik}}^{\star}\mathbf{v}^{(\mathbf{l})})\right|}{2\left({\scriptstyle\Omega}_{\mathbf{l}}{\scriptstyle\Omega}_{\mathbf{m}}\right)^{1/2}}$
(47) $\displaystyle=$ $\displaystyle\frac{1}{2{\scriptstyle
N}^{d}}\frac{\left|e^{j\theta\mathbf{m}\cdot\mathbf{i}}\right|\left|1-e^{\text{-}jm_{1}\theta}\right|}{{\scriptstyle\Omega}_{\mathbf{m}}^{1/2}}\frac{\left|e^{j\theta\mathbf{l}\cdot\mathbf{i}}\right|\left|1-e^{jl_{1}\theta}\right|}{{\scriptstyle\Omega}_{\mathbf{l}}^{1/2}}$
$\displaystyle=$ $\displaystyle\frac{1}{2{\scriptstyle
N}^{d}}\frac{\left|1-e^{\text{-}jm_{1}\theta}\right|}{{\scriptstyle\Omega}_{\mathbf{m}}^{1/2}}\frac{\left|1-e^{jl_{1}\theta}\right|}{{\scriptstyle\Omega}_{\mathbf{l}}^{1/2}}$
(48)

This is analogous to expression (39) for the one-dimensional ring, except in
that we derive a lower bound on the slope of the widest Arnold tongue instead
of an exact estimate. This lower bound is inversely proportional to the total
number of nodes (${\scriptstyle N}^{d}$ here), and attains its maximum of
$1/{\scriptstyle N}^{d}$ when $m_{1}=l_{1}=\left\lfloor{\scriptstyle
N}/2\right\rfloor$, $m_{q}=l_{q}=0,\ q=2,\dots,d$. With the number of nodes
fixed, this behavior is identical to the one-dimensional case and independent
of the dimension of the lattice. All critical frequencies are contained in
$\omega\in(0,4d^{1/2}]$, and the widest Arnold tongue in the parameter space
$(\omega,\epsilon)$ thus has slope larger than $1/{\scriptstyle N}^{d}$.

Fig. 7 shows the $\mathcal{O}(\epsilon)$ approximation to the Arnold tongues
for ring lattices of various sizes with $d=2$. As a result of the eigenvalue
multiplicity, there are fewer critical frequencies, and the corresponding
Arnold tongues are wider than for the one-dimensional ring graph with the same
number of nodes.

Figure 7: Comparison of the stability diagram for ring lattices of dimension
$1$ (left) and $2$ (right). The top row is for graphs with 25 nodes, the
bottom row for graphs with 64. The widest Arnold tongue in each case is
colored purple. The Laplacian for a two-dimensional lattice has much higher
eigenvalue multiplicity, the effect of which is to produce fewer critical
frequencies and wider Arnold tongues. (Compare with Fig. 8.)

###  IV-C The Effect of Connectivity

Narrower Arnold tongues are harder to observe in physical systems, especially
in the presence of damping. Thus Eq. 46 suggests that the increasing the
connectivity of the graph has a destabilizing effect, as the tongues are wider
(Fig. 7). However, this is a consequence of the effect of the eigenvalue
multiplicity, a non-generic factor, and not the increased connectivity of the
graph. We investigate this by perturbing the edge weights of the ring graph
slightly, which makes the eigenvalues distinct in the generic case. This turns
the lower bound calculations for the slope of the widest tongues (Eq. 38 and
Eq. 47) into equalities, and ring graph connectivity has little to no effect
on the distribution of tongue slopes. Fig. 8 compares the stability diagrams
for one and two-dimensional ring graphs with the same number of nodes but with
the effect of eigenvalue multiplicity removed.

Figure 8: Comparison of the stability diagram (to $\mathcal{O}(\epsilon)$) for
ring lattices of dimension $1$ (left) and $2$ (right) without eigenvalue
multiplicity. The top row is for graphs with 25 nodes, the bottom row for
graphs with 64, and the nominal edge weights of $1$ are perturbed by a normal
variable with standard deviation $0.01$, which makes the Laplacian eigenspaces
one-dimensional in the generic case. The widest Arnold tongue in each case is
colored purple. With the effects of eigenvalue multiplicity removed, the
slopes of the widest tongues for the two cases are very close to each other,
implying a negligible effect of graph connectivity on the stability of the
dynamics. The widest tongues are no longer guaranteed to be near $\omega=4$,
since the graph Laplacians are not circulant. Compare with Fig. 7.

Fig. 9 compares all the slopes and the number of first-order Arnold tongues
for one and two-dimensional ring graphs with perturbed edge weights and
various graph sizes. The increased connectivity of the $2$-dimensional lattice
does not affect the stability of the dynamics.

Figure 9: (Above) Slopes of all Arnold tongues for one and two-dimensional
ring graphs of the same size with perturbed edge weights, as a function of the
size of the network. The two are offset slightly for legibility. (Below) The
number of Arnold tongues for both networks as a function of the network size.
There is no significant difference in the distribution of slopes or the number
of distinct Arnold tongues between the two kinds of graphs.

##  V Conclusion and Further Questions

We draw the following tentative conclusions about parametric resonance via
periodic edge perturbations in networked oscillators:

  * •

The stability characteristics of parametric resonance via edge perturbations
of linear oscillator dynamics on undirected graphs can be inferred from the
spectrum of the graph Laplacian. The critical parametric forcing frequencies
are sums of pairs of the natural frequencies of the graph, which are
themselves the square roots of the Laplacian eigenvalues.

  * •

Each critical frequency is thus associated with two eigenvalues. The slope or
width of the Arnold tongue structure in parameter space originating at that
frequency is a function of the corresponding eigenvectors. It is large when a
corresponding eigenvector varies significantly between the nodes connecting
the perturbed edges, or when the eigenvalues have high multiplicity. If either
eigenvector has the same value on these nodes, there is no Arnold tongue at
the corresponding frequency and the dynamics remains stable.

  * •

The Arnold tongue slopes are generally proportional to the number of edges
being perturbed, and inversely proportional to the size of the network. They
do not appear to depend strongly on the degree of network connectivity.

While we choose parametric forcing at a single frequency to illustrate the
method, these results extend to more general periodic forcing. The
perturbation analysis treats the forcing as additive input to a linear system,
so each harmonic of the forcing period acts as a separate input. The system
thus experiences parametric resonance if any harmonic of the forcing frequency
equals the sum of two oscillator natural frequencies.

Parametric resonance in real-world oscillator networks differs from the models
considered in this work because of damping. In parametric oscillators the
effect of damping on the stability diagrams is to “raise” the Arnold tongues
(Fig. 1), an effect more pronounced in narrower or higher-order tongues. Thus
only the resonances corresponding to the wider and lower-order tongues are
observable in physical systems. Approximating these realizable tongue
positions and widths for the purposes of design, to avoid or excite parametric
resonance, can be a computationally intensive undertaking for large networks.
The multiple time-scale perturbation analysis provides a simple criterion for
the same that depends only on knowledge of the graph Laplacian.

A natural extension of the analysis presented here is thus to consider the
effect of damping. In its current form, the perturbation analysis depends on a
diagonalization argument and can be readily extended to cover only very
specific cases of damping: the damping coefficient matrix must commute with
the graph Laplacian. This is the case, for example, when each oscillator
experiences “self-damping” with the same damping coefficient. A new approach
is required for a more general treatment.

It would be useful to know if the multiple scale perturbation method can
predict the unstable modes in addition to the positions and slopes of the
Arnold tongues. Preliminary analysis suggests that this is the case. To first
order in the forcing strength, the position and slope of each Arnold tongue
depends on two modes (eigenvalues and eigenvectors) of the graph Laplacian,
and the unstable slow-time dynamics in these tongues is a combination of these
modes as well (see Section A-C). Fig. 10 compares, for a ring graph, the
actual unstable mode obtained via numerical computation of the Monodromy map
with the corresponding modes of the Laplacian used in the perturbation
analysis. We observe a qualitative similarity between the actual and predicted
mode shapes.

Figure 10: Comparison of the unstable mode shapes as predicted by the
numerically computed Monodromy map (left, large circle) and the multiscale
perturbation analysis (right, small circles) for a $64$ node ring graph with
one time-periodically varying edge weight, at two critical parametric forcing
frequencies (Eq. 30). The effect of eigenvalue multiplicity has been removed
by perturbing the nominal edge weights slightly, and the value of the mode at
each node is plotted radially. Above and below: this comparison is shown at
two critical parametric forcing frequencies. At both frequencies, the mode
shape as predicted by the perturbation analysis appears to be in the span of
the two eigenvectors of the Laplacian corresponding to the critical frequency
(Section A-C). We observe that the true unstable modes are visually similar to
one or a combination of the two Laplacian eigenvectors.

## References

  * [1] T. Fossen and H. Nijmeijer, _Parametric resonance in dynamical systems_. Springer Science & Business Media, 2011.
  * [2] K. Moran, C. Burgner, S. Shaw, and K. Turner, “A review of parametric resonance in microelectromechanical systems,” _Nonlinear theory and its applications, IEICE_ , vol. 4, no. 3, pp. 198–224, 2013.
  * [3] A. H. Nayfeh and D. T. Mook, _Nonlinear Oscillations_. Wiley, May 1995.
  * [4] W. Magnus and S. Winkler, _Hill’s Equation_. Courier Corporation, Jan. 2004.
  * [5] J. M. Berg and I. M. Wickramasinghe, “Vibrational control without averaging,” _Automatica_ , vol. 58, pp. 72–81, Aug. 2015.
  * [6] F. Bullo, “Averaging and vibrational control of mechanical systems,” _SIAM Journal on Control and Optimization_ , vol. 41, no. 2, pp. 542–562, 2002.
  * [7] S. Meerkov, “Principle of vibrational control: Theory and applications,” _IEEE Transactions on Automatic Control_ , vol. 25, no. 4, pp. 755–762, 1980\.
  * [8] Y. Qin, D. S. Bassett, and F. Pasqualetti, “Vibrational Control of Cluster Synchronization: Connections with Deep Brain Stimulation,” Apr. 2022.
  * [9] M. H. Histed, V. Bonin, and R. C. Reid, “Direct activation of sparse, distributed populations of cortical neurons by electrical microstimulation,” _Neuron_ , vol. 63, no. 4, pp. 508–522, 2009.
  * [10] A. M. Nobili, Y. Qin, C. A. Avizzano, D. S. Bassett, and F. Pasqualetti, “Vibrational stabilization of complex network systems,” in _2023 American Control Conference (ACC)_. IEEE, 2023, pp. 1980–1985.
  * [11] F. Dörfler and F. Bullo, “Synchronization in complex networks of phase oscillators: A survey,” _Automatica_ , vol. 50, no. 6, pp. 1539–1564, 2014\.
  * [12] B. Bamieh, M. Jovanovic, P. Mitra, and S. Patterson, “Effect of topological dimension on rigidity of vehicle formations: Fundamental limitations of local feedback,” in _2008 47th IEEE Conference on Decision and Control_. IEEE, 2008, pp. 369–374.
  * [13] Y. Chalopin, F. Piazza, S. Mayboroda, C. Weisbuch, and M. Filoche, “Universality of fold-encoded localized vibrations in enzymes,” _Scientific reports_ , vol. 9, no. 1, p. 12835, 2019.
  * [14] C. M. Bender and S. A. Orszag, _Advanced Mathematical Methods for Scientists and Engineers I_. New York, NY: Springer New York, 1999.
  * [15] B. Bamieh, M. R. Jovanovic, P. Mitra, and S. Patterson, “Coherence in Large-Scale Networks: Dimension-Dependent Limitations of Local Feedback,” _IEEE Transactions on Automatic Control_ , vol. 57, no. 9, pp. 2235–2249, Sep. 2012.
  * [16] V. Yakubovich and V. Starzhinskii, _Linear Differential Equations with Periodic Coefficients_ , 1st ed. Wiley, 1975, vol. 1.
  * [17] J. P. Hespanha, _Linear Systems Theory_. Princeton University Press, Feb. 2018.

##  Appendix A Appendix

###  A-A Vector-valued multiple scale perturbation analysis

When one of the conditions (24) or (25) is met, we investigate the behavior of
the system near a critical $\omega$.

To facilitate the multiple time-scale analysis, we rescale time so that the
parametric forcing frequency $\omega$ is uniformly $1$. Under the
transformation $\eta\leftarrow\omega t$, we get

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}$
$\displaystyle=\omega\frac{\mathrm{d}}{\mathrm{d}\eta}$
$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}$
$\displaystyle=\omega^{2}\frac{\mathrm{d}^{2}}{\mathrm{d}\eta^{2}}$

so that system (12) becomes

$\displaystyle\ddot{\phi}(t)$ $\displaystyle+\left(L+\epsilon
E_{ik}\cos(\omega t)\right)\phi(t)=0$ $\displaystyle\updownarrow$
$\displaystyle\phi^{\prime\prime}(\eta)$
$\displaystyle+\frac{1}{\omega^{2}}\left(L+\epsilon
E_{ik}\cos(t)\right)\phi(\eta)=0$ (49)

Defining $\kappa:=\frac{1}{\omega^{2}}$, we get

$\displaystyle\phi^{\prime\prime}(\eta)+\left(\kappa L+\kappa\epsilon
E_{ik}\cos(t)\right)\phi(\eta)=0$ (50)

For certain values of $\kappa$, this system is unstable as $\epsilon\to 0$.
From (23), this is when

$\displaystyle\pm j\sqrt{\kappa}{\scriptstyle\Omega}_{l}=\pm
j\sqrt{\kappa}{\scriptstyle\Omega}_{k}\pm j$
$\displaystyle\implies\sqrt{\kappa}\left|{\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{k}\right|=1$
(51)
$\displaystyle\implies\sqrt{\kappa}\left|{\scriptstyle\Omega}_{l}-{\scriptstyle\Omega}_{k}\right|=1$
(52)

for some $l,k=0,1,\dots,n$. Before proceeding further, we switch our notation
back to $t$ from $\eta$. Our rescaled system is thus

$\displaystyle\ddot{\phi}(t)+\kappa\left(L+\epsilon
E_{ik}\cos(t)\right)\phi(t)=0$ (53)

This system exhibits the instabilities studied in the previous section for
specific critical values of $\kappa$ (since $\kappa=\frac{1}{\omega^{2}}$). We
investigate the stability of the system (for small $\epsilon$) in the vicinity
of these $\kappa$. Specifically, we are interested in the behavior of the
system (53) along curves in $\kappa,\epsilon$ space:

$\kappa(\epsilon)=\kappa_{0}+\kappa_{1}\epsilon+\kappa_{2}\epsilon^{2}+\mathcal{O}(\epsilon^{3})$
(54)

This formulation allows us to estimate not only the critical values of
$\kappa_{0}$ for instability, but also the slopes (or widths) of the
corresponding Arnold tongues $\kappa_{1}$. Figure 11 shows visualizations of
$\kappa(\epsilon)$.

Figure 11: Examples of test curves $\kappa(\epsilon)$ used for multiple time-
scale perturbation analysis of system (53), overlaid on top of a stability
diagram for a linear time-periodic oscillator in $(\kappa,\epsilon)$ space. In
gray are the Arnold tongues for the system. The test curves are represented in
the diagram as the colored streaks, and algebraically using Equation (54). We
expect the stability of the system at points on the test curves for small
$\epsilon$ to depend on both the $\kappa$-intercept $\kappa_{0}$ and the slope
$\kappa_{1}$ measured from the vertical. Along the green curve (with
$\kappa_{0}=1.3$), the system is stable for all possible slopes $\kappa_{1}$.
For the curves with $\kappa_{0}=0.5$, the stability of the system depends on
the slope $\kappa_{1}$ – the system is unstable at all points on the red curve
and stable at all points on the cyan curve. The objective of the multiple
time-scale analysis is to estimate both the critical points $\kappa_{0}$ and
the critical slopes $\kappa_{1}$, or their equivalents in the original scaling
in the original $(\omega,\epsilon)$ space (Equation (12)).

We treat the configuration of the system $\phi(t)$ as a function of two time
scales, a fast time $t$ and a slow time scale $\tau:=\epsilon t$. This allows
us to study the effect of the slow time behavior of the system on the fast
time dynamics. Specifically, we set

$\Phi(t,\tau;\epsilon):=\phi(t;\epsilon)$ (55)

Note that it is always possible to find a $\Phi(t,\tau;\epsilon)$ that
satisfies this relation – the trivial case is when $\Phi$ does not depend on
its second argument $\tau$. Since this decomposition is non-unique,
$\Phi(t,\tau;\epsilon)$ is not well-defined yet, and will be uniquely
determined in the following analysis.

We can now expand $\Phi(t,\tau;\epsilon)$ in powers of $\epsilon$:

$\displaystyle\Phi(t,\tau;\epsilon)=\sum\limits_{n=0}^{N}\epsilon^{n}\Phi_{n}(t,\tau)+\mathcal{O}(\epsilon^{N+1})$
(56)

We relate the derivatives as follows:

$\displaystyle\frac{\mathrm{d}\phi}{\mathrm{d}t}(t;\epsilon)$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}t}\Phi(t,\tau;\epsilon)=\frac{\partial\Phi}{\partial
t}+\frac{\partial\Phi}{\partial\tau}\frac{\mathrm{d}\tau}{\mathrm{d}t}$
$\displaystyle=\frac{\partial\Phi}{\partial
t}+\epsilon\frac{\partial\Phi}{\partial\tau}$
$\displaystyle\frac{\mathrm{d}^{2}\phi}{\mathrm{d}t^{2}}(t;\epsilon)$
$\displaystyle=\frac{\partial^{2}\Phi}{\partial
t^{2}}+2\epsilon\frac{\partial\Phi}{\partial
t}\frac{\partial\Phi}{\partial\tau}+\frac{\partial^{2}\Phi}{\partial\tau^{2}}$

or, in more compact notation

$\displaystyle\partial_{1}(\cdot)$
$\displaystyle:=\frac{\partial(\cdot)}{\partial t},\
\partial_{2}(\cdot):=\frac{\partial(\cdot)}{\partial\tau}$
$\displaystyle\implies\frac{\mathrm{d}}{\mathrm{d}t}\phi$
$\displaystyle=\partial_{1}\Phi+\epsilon\partial_{2}\Phi$
$\displaystyle\frac{\mathrm{d}^{2}\phi}{\mathrm{d}t^{2}}$
$\displaystyle=\partial_{1}^{2}\Phi+2\epsilon\partial_{1}\partial_{2}\Phi+\partial_{2}^{2}\Phi$

Together, (53) is transformed to

$\displaystyle\partial_{1}^{2}\Phi+(\kappa_{0}+\kappa_{1}\epsilon+\dots)$
$\displaystyle(L+E_{ik}\cos(t))\Phi=$
$\displaystyle-2\epsilon\partial_{1}\partial_{2}\Phi-\epsilon^{2}\partial_{2}^{2}\Phi$
(57)

Substituting for $\kappa$ from equation (54), for $\Phi(t,\tau;\epsilon)$ from
the series expansion (56) and collecting powers of $\epsilon$, we get to
$\mathcal{O}(\epsilon)$:

$\displaystyle\mathcal{O}(\epsilon^{0}):\quad$
$\displaystyle\partial_{1}^{2}\Phi_{0}(t,\tau)+\kappa_{0}L\Phi_{0}(t,\tau)=0$
(58) $\displaystyle\mathcal{O}(\epsilon^{1}):\quad$
$\displaystyle\partial_{1}^{2}\Phi_{1}(t,\tau)+\kappa_{0}L\Phi_{1}(t,\tau)=$
$\displaystyle\quad\text{-}\left(\kappa_{1}L+\kappa_{0}E_{ik}\cos(t)+2\partial_{1}\partial_{2}\right)\Phi_{0}(t,\tau)$
(59)

As in the regular perturbation analysis (19), the $\mathcal{O}(\epsilon^{0})$
system is a linear, unforced harmonic oscillator whose solution is a matrix
exponential:

$\Phi_{0}(t,\tau)=V\left(e^{j\sqrt{\kappa_{0}}\Omega
t}\mathbf{c}(\tau)+e^{\text{-}j\sqrt{\kappa_{0}}\Omega
t}\mathbf{d}(\tau)\right)$ (60)

However, the “initial conditions” $\mathbf{c}(\tau)$ and $\mathbf{d}(\tau)$
are functions of the slow time variable $\tau$ and yet to be determined.

To analyze the $\mathcal{O}(\epsilon^{1})$ system, it will be useful to work
in the basis of the eigenvectors of $L$. With $LV=V\Omega^{2}$ and
$V^{\star}=V^{\text{-}1}$,

$\displaystyle V^{\star}\Phi_{0}(t,\tau)$
$\displaystyle=e^{j\sqrt{\kappa_{0}}\Omega
t}\mathbf{c}(\tau)+e^{\text{-}j\sqrt{k_{0}}\Omega t}\mathbf{d}(\tau)$ (61)

and the $\mathcal{O}(\epsilon^{1})$ system is

$\displaystyle\partial_{1}^{2}(V^{\star}\Phi_{1}(t,\tau))+$
$\displaystyle\kappa_{0}V^{\star}LV(V^{\star}\Phi_{1}(t,\tau))=$
$\displaystyle\text{-}V^{\star}(\kappa_{1}L+\kappa_{0}E_{ik}\cos(t)+2\partial_{1}\partial_{2})V(V^{\star}\Phi_{0}).$
(62)

Substituting for $V^{\star}\Phi_{0}$ from (61), the additive forcing term
becomes

$\displaystyle-V^{\star}\left(\kappa_{1}L+\kappa_{0}E_{ik}\cos(t)+2\partial_{1}\partial_{2}\right)V\left(V^{\star}\Phi_{0}\right)$
$\displaystyle=$
$\displaystyle-\left(\kappa_{1}\Omega^{2}+\kappa_{0}V^{\star}E_{ik}V\cos(t)+2\partial_{1}\partial_{2}I\right)\times$
$\displaystyle\quad(e^{j\sqrt{\kappa_{0}}\Omega
t}\mathbf{c}(\tau)+e^{\text{-}j\sqrt{\kappa_{0}}\Omega t}\mathbf{d}(\tau))$
$\displaystyle=$ $\displaystyle-e^{j\sqrt{\kappa_{0}}\Omega
t}\left(\kappa_{1}\Omega^{2}\,\mathbf{c}(\tau)+2j\sqrt{\kappa_{0}}\Omega\,\mathbf{c}^{\prime}(\tau)\right)$
$\displaystyle-e^{\text{-}j\sqrt{\kappa_{0}}\Omega
t}\left(\kappa_{1}\Omega^{2}\,\mathbf{d}(\tau)-2j\sqrt{\kappa_{0}}\Omega\,\mathbf{d}^{\prime}(\tau)\right)$
$\displaystyle-\frac{\kappa_{0}}{2}V^{\star}E_{ik}V\left(e^{j(\sqrt{\kappa_{0}}\Omega+I)t}+e^{(\sqrt{\kappa_{0}}\Omega-I)t}\right)\mathbf{c}(\tau)$
$\displaystyle+\frac{\kappa_{0}}{2}V^{\star}E_{ik}V\left(e^{\text{-}j(\sqrt{\kappa_{0}}\Omega+I)t}+e^{\text{-}j(\sqrt{\kappa_{0}}\Omega-I)t}\right)\mathbf{d}(\tau)$
(63)

where we have used the facts that

  * •

$e^{j\sqrt{\kappa_{0}}\Omega t}$ commutes with $\kappa_{1}\Omega^{2}$ by the
properties of the matrix exponential, and

  * •

$\cos(t)=\frac{1}{2}\left(e^{jt}+e^{\text{-}jt}\right)$.

For stability, we require both the fast and slow time dynamics to be stable.
The former requires choices of $\mathbf{c}(\tau)$ and $\mathbf{d}(\tau)$ for
which the forcing terms that cause harmonic resonance in system (59) (and are
_secular_) vanish. This uniquely determines $\Phi_{0}(0,\tau)$. The latter
imposes conditions on $\mathbf{c}(\tau)$ and $\mathbf{d}(\tau)$, and in turn
on the system parameters $\omega$ and $\epsilon$.

We make the following observations:

  1. 1.

Every term involving $e^{j\sqrt{\kappa_{0}}\Omega t}$ or
$e^{\text{-}j\sqrt{\kappa_{0}}\Omega t}$ is secular, since these are modes of
the unforced system.

  2. 2.

If $\kappa_{0}$ corresponds to $\omega$ at the base of an Arnold tongue
(_i.e._ is “critical”), then $j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}=\pm
j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}\pm 1$ for some $l,m$. Then the
corresponding term involving
$\exp\left(j\left(\pm\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}\pm
1\right)t\right)$ is secular.

We first consider the simpler case, _i.e_ of non-critical $\kappa_{0}$.

####  A-A1 In the stable region

When $\kappa_{0}$ is not critical, all forcing terms causing harmonic
resonance in system (59) are of the form

$\displaystyle{\scriptstyle\Omega}_{l}\left[c_{l}^{\prime}(\tau)-\frac{j}{2}\frac{\kappa_{1}}{\kappa_{0}}\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}_{l}c_{l}(\tau)\right]\exp\left(j\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}_{l}t\right)$
$\displaystyle{\scriptstyle\Omega}_{l}\left[d_{l}^{\prime}(\tau)+\frac{j}{2}\frac{\kappa_{1}}{\kappa_{0}}\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}_{l}d_{l}(\tau)\right]\exp\left(-j\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}_{l}t\right)$

where $c_{l}$ is the $l^{\mathrm{th}}$ component of $\mathbf{c}(\tau)$. There
are $n$ such components of $\mathbf{c}$ and $\mathbf{d}$, corresponding to the
exponentials of $\pm j\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}_{l}t$,
$l=1,2,\dots,n$. By assumption, none of the terms involving
$\exp\left(j\left(\pm\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}\pm
1\right)t\right)$ are secular.

If we pick $\mathbf{c}(\tau)$ and $\mathbf{d}(\tau)$ such that these
(independently) vanish, then $\Phi_{1}(t,\tau)$ is stable in $t$, and we fully
determine $\Phi_{0}(t,\tau)$ in the process. In the parlance of perturbation
theory, we eliminate the secular terms in equation (59).

The above terms vanish when $\mathbf{c}(\tau)$ and $\mathbf{d}(\tau)$ satisfy
the ODEs

$\displaystyle\mathbf{c}^{\prime}(\tau)-j\frac{\kappa_{1}}{2\kappa_{0}}\kappa_{0}^{\frac{1}{2}}\Omega\mathbf{c}(\tau)=0$
$\displaystyle\mathbf{d}^{\prime}(\tau)+j\frac{\kappa_{1}}{2\kappa_{0}}\kappa_{0}^{\frac{1}{2}}\Omega\mathbf{d}(\tau)=0$
$\displaystyle\implies$
$\displaystyle\mathbf{c}^{\prime}(\tau)=j\frac{\kappa_{1}}{2\kappa_{0}^{\frac{1}{2}}}\Omega\mathbf{c}(\tau)$
(64)
$\displaystyle\mathbf{d}^{\prime}(\tau)=-j\frac{\kappa_{1}}{2\kappa_{0}^{\frac{1}{2}}}\Omega\mathbf{d}(\tau)$
(65)

The solution to this is

$\displaystyle\mathbf{c}(\tau)=\exp{\left(j\tfrac{\kappa_{1}}{2\kappa_{0}^{\frac{1}{2}}}\Omega\tau\right)}\mathbf{c}(0),$
$\displaystyle\quad\mathbf{d}(\tau)=\exp{\left(\text{-}j\tfrac{\kappa_{1}}{2\kappa_{0}^{\frac{1}{2}}}\Omega\tau\right)}\mathbf{d}(0)$

This choice of $\mathbf{c}(\tau)$ and $\mathbf{d}(\tau)$ makes
$\Phi_{1}(t,\tau)$ stable. Since $\Omega$ is diagonal with non-negative
entries, $\mathbf{c}(\tau)$ and $\mathbf{d}(\tau)$ are bounded, and
$\Phi_{0}(t,\tau)$ is stable as well. This is the case for all $\kappa_{1}$,
_i.e._ irrespective of our choice of test curve in $(\kappa,\epsilon)$ space.
This is in line with our expectations from the analysis of section III-A,
since every critical $\kappa_{0}$ should satisfy one of the conditions (51) or
(52). Our test curve $\kappa(\epsilon)$ corresponds to the green curve in
figure 11. When starting from a non-critical $\kappa_{0}$, the system remains
stable along all curves in the $(\epsilon,\kappa)$ parameter space (and thus
in the $(\epsilon,\omega)$ domain) as $\epsilon\to 0$.

Further, $\Phi_{0}(t,\tau)$ is given by

$\displaystyle
V^{\star}\left(e^{j\sqrt{\kappa_{0}}\Omega\left(1+\epsilon\frac{\kappa_{1}}{2\kappa_{0}}\right)t}\mathbf{c}(0)+e^{-j\sqrt{\kappa_{0}}\Omega\left(1+\epsilon\frac{\kappa_{1}}{2\kappa_{0}}\right)t}\mathbf{d}(0)\right)$
(66)

which is a sinusoidal response at a perturbed frequency. The change in the
frequency is proportional to $\epsilon$ and the slope $\kappa_{1}$.

####  A-A2 At critical frequencies $\kappa_{0}$,
$\sqrt{\kappa_{0}}\left|{\scriptstyle\Omega}_{l}\pm{\scriptstyle\Omega}_{m}\right|=1$

When $\kappa_{0}$ is critical, _i.e._ when either (51) or (52) is true,
forcing terms in the $\mathcal{O}(\epsilon)$ dynamics (59) involving
$\exp\left(j\left(\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}\pm
1\right)t\right)$ also cause harmonic resonance and are secular. Suppose one
of the relations
$\sqrt{\kappa_{0}}\left|{\scriptstyle\Omega}_{l}\pm{\scriptstyle\Omega}_{m}\right|=1$
holds for a given pair of indices $l$ and $m$. Restricting our attention to
equations in (62) with forcing terms involving the frequencies
$\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}$ and
$\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}$, we get

$\displaystyle\partial_{1}^{2}\left(V^{\star}\Phi_{0}(t,\tau)\right)_{m}+\kappa_{0}{\scriptstyle\Omega}_{m}^{2}\left(V^{\star}\Phi_{0}(t,\tau)\right)_{m}=$
$\displaystyle-e^{j\sqrt{\kappa_{0}}\Omega_{m}t}\left[\kappa_{1}{\scriptstyle\Omega}^{2}_{m}+2j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}c_{m}^{\prime}(\tau)\right]$
$\displaystyle-\tfrac{1}{2}\kappa_{0}\mathbf{v}^{\star}_{m}E_{ik}\mathbf{v}_{l}c_{l}(\tau)$
$\displaystyle\times\left(e^{j(\sqrt{\kappa_{0}}\Omega_{l}+1)t}+e^{j(\sqrt{\kappa_{0}}\Omega_{l}\text{-}1)t}\right)$
(67)
$\displaystyle\partial_{1}^{2}\left(V^{\star}\Phi_{0}(t,\tau)\right)_{l}+\kappa_{0}{\scriptstyle\Omega}_{l}^{2}\left(V^{\star}\Phi_{0}(t,\tau)\right)_{l}=$
$\displaystyle-e^{j\sqrt{\kappa_{0}}\Omega_{l}t}\left[\kappa_{1}{\scriptstyle\Omega}^{2}_{l}+2j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}c_{l}^{\prime}(\tau)\right]$
$\displaystyle-\tfrac{1}{2}\kappa_{0}\mathbf{v}^{\star}_{l}E_{ik}\mathbf{v}_{m}c_{m}(\tau)$
$\displaystyle\times\left(e^{j(\sqrt{\kappa_{0}}\Omega_{m}+1)t}+e^{j(\sqrt{\kappa_{0}}\Omega_{m}\text{-}1)t}\right)$
(68)

with similar terms involving $\mathbf{d}(\tau)$ and the frequencies
$\text{-}\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}$ and
$\text{-}\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}$. Here we have assumed that
the eigenspaces corresponding to eigenvalues ${\scriptstyle\Omega}_{l}$ and
${\scriptstyle\Omega}_{m}$ are $1$-dimensional. The case of eigenvalue
multiplicity is handled in Section A-B.

To find conditions on $c_{l}(\tau)$ and $c_{m}(\tau)$ such that
$\Phi_{1}(t,\tau)$ is stable in $t$, we need to collect terms involving like
exponentials. This organization differs between the two cases
$\sqrt{\kappa_{0}}\left|{\scriptstyle\Omega}_{l}\pm{\scriptstyle\Omega}_{m}\right|=1$.

##### Case I:
$\sqrt{\kappa_{0}}\left|{\scriptstyle\Omega}_{l}-{\scriptstyle\Omega}_{m}\right|=1$

This condition applies when the difference between the square roots of
eigenvalues of the Laplacian $L$ equals the original parametric forcing
frequency $\omega$, since
$\left|{\scriptstyle\Omega}_{l}-{\scriptstyle\Omega}_{m}\right|=1/\sqrt{\kappa_{0}}=\omega$.

Without loss of generality, we can pick
${\scriptstyle\Omega}_{l}\geq{\scriptstyle\Omega}_{m}$, so that

  * •

$\exp\left(j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}t\right)=\exp\left(j\left(\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}+1\right)t\right)$
and

  * •

$\exp\left(j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}t\right)=\exp\left(j\left(\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}-1\right)t\right)$,

corresponding to forcing terms in (67) and (68):

$\displaystyle\left[c_{l}^{\prime}(\tau)-j\tfrac{\kappa_{1}}{2\kappa_{0}}\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}_{l}c_{l}(\tau)-j\,\tfrac{1}{4}\kappa_{0}^{\text{-}\frac{1}{2}}{\scriptstyle\Omega}_{l}^{\text{-}1}\mathbf{v}_{l}^{\star}E_{ik}\mathbf{v}_{m}c_{m}(\tau)\right]$
$\displaystyle\times\text{-}2j{\scriptstyle\Omega}_{l}e^{\left(j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}t\right)}$
(69)
$\displaystyle\left[c_{m}^{\prime}(\tau)-j\tfrac{\kappa_{1}}{2\kappa_{0}}\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}_{m}c_{m}(\tau)-j\tfrac{1}{4}\kappa_{0}^{\text{-}\frac{1}{2}}{\scriptstyle\Omega}_{m}^{\text{-}1}\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}c_{l}(\tau)\right]$
$\displaystyle\times\text{-}2j{\scriptstyle\Omega}_{m}e^{\left(j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}t\right)}$
(70)

The forcing terms corresponding to all other modes
(${\scriptstyle\Omega}\neq{\scriptstyle\Omega}_{l},{\scriptstyle\Omega}_{m}$)
are of the form
$\left[c^{\prime}(\tau)-j\tfrac{\kappa_{1}}{2\kappa_{0}}\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}\mathbf{c}(\tau)\right]e^{\left(j\sqrt{\kappa_{0}}{\scriptstyle\Omega}t\right)}$,
for which the analysis from Section A-A1 carries over.

We pick $c_{l}(\tau)$ and $c_{m}(\tau)$ such that both these terms vanish.

$\displaystyle
c_{l}^{\prime}(\tau)\,\text{-}\,j\tfrac{\kappa_{1}}{2\kappa_{0}}\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}_{l}c_{l}(\tau)\,\text{-}\,j\tfrac{1}{4}\kappa_{0}^{\text{-}\frac{1}{2}}{\scriptstyle\Omega}_{m}^{\text{-}1}\mathbf{v}_{l}^{\star}E_{ik}\mathbf{v}_{m}c_{m}(\tau)=0$
(71) $\displaystyle
c_{m}^{\prime}(\tau)\,\text{-}\,j\tfrac{\kappa_{1}}{2\kappa_{0}}\kappa_{0}^{\frac{1}{2}}{\scriptstyle\Omega}_{m}c_{m}(\tau)\,\text{-}\,j\tfrac{1}{4}\kappa_{0}^{\text{-}\frac{1}{2}}{\scriptstyle\Omega}_{l}^{\text{-}1}\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}c_{l}(\tau)=0$
(72)

The complex conjugate terms combine similarly:

  * •

$\exp\left(-j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}t\right)=\exp\left(j\left(-\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}-1\right)t\right)$
and

  * •

$\exp\left(-j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}t\right)=\exp\left(j\left(-\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}+1\right)t\right)$,

giving us complex conjugates of equations (71) and (72) but involving
$d_{l}(\tau)$ and $d_{m}(\tau)$.

We can rewrite this system of equations as:

$\displaystyle\begin{bmatrix}c_{l}^{\prime}(\tau)\\\
c_{m}^{\prime}(\tau)\end{bmatrix}$
$\displaystyle=\frac{j}{2\sqrt{\kappa_{0}}}\underbrace{\begin{bmatrix}\kappa_{1}{\scriptstyle\Omega}_{l}&\frac{\kappa_{0}(\mathbf{v}_{l}^{\star}E_{ik}\mathbf{v}_{m})}{2{\scriptstyle\Omega}_{l}}\\\
\frac{\kappa_{0}(\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l})}{2{\scriptstyle\Omega}_{m}}&\kappa_{1}{\scriptstyle\Omega}_{m}\end{bmatrix}}_{A}\begin{bmatrix}c_{l}(\tau)\\\
c_{m}(\tau)\end{bmatrix}$ (73) $\displaystyle\text{or
}\mathbf{c}^{\prime}(\tau)$ $\displaystyle=\frac{j}{2\sqrt{\kappa_{0}}}A\
\mathbf{c}(\tau)$

This choice of $\mathbf{c}(\tau)$ makes $\Phi_{1}(t,\tau)$ stable in $t$. For
bounded slow-time behavior of $\Phi_{0}(t,\tau)$, we require
$\mathbf{c}(\tau)$ to be stable, or the eigenvalues of $A$ to be real (because
of the $j$ prefactor). The characteristic polynomial of $A$ is

$\displaystyle
y^{2}-\underbrace{\left({\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}\right)\kappa_{1}}_{\text{trace}(A)}y$
$\displaystyle+\underbrace{\left(\kappa_{1}^{2}{\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}-(\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l})^{2}\frac{\kappa_{0}^{2}}{4{\scriptstyle\Omega}_{m}{\scriptstyle\Omega}_{l}}\right)}_{\text{determinant}(A)}=0.$

Requiring its discriminant to be real, we get

$\displaystyle\left({\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}\right)^{2}\kappa_{1}^{2}-4\kappa_{1}^{2}{\scriptstyle\Omega}_{1}{\scriptstyle\Omega}_{m}$
(74)
$\displaystyle+(\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l})^{2}\frac{\kappa_{0}^{2}}{{\scriptstyle\Omega}_{m}{\scriptstyle\Omega}_{l}}\geq
0$ $\displaystyle\implies$
$\displaystyle\left({\scriptstyle\Omega}_{l}-{\scriptstyle\Omega}_{m}\right)^{2}\kappa_{1}^{2}+(\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l})^{2}\frac{\kappa_{0}^{2}}{{\scriptstyle\Omega}_{m}{\scriptstyle\Omega}_{l}}\geq
0$ (75)

This discriminant is always non-negative, implying that there is no
instability at $(\kappa_{0},\epsilon)$, or along any curve
$\kappa(\epsilon)=\kappa_{0}+\kappa_{1}\epsilon$. Despite our initial
assessment, the condition
$\sqrt{\kappa_{0}}\left|{\scriptstyle\Omega}_{l}-{\scriptstyle\Omega}_{m}\right|=1$
does not cause instability. This is one of the sources of false positives in
the regular perturbation analysis of section III-A. An alternative derivation
of this result is provided by the Krein-Gel’fand-Lidskii theorem [16], in
terms of the stability properties of perturbed time-dependent Hamiltonian
systems.

##### Case II:
$\sqrt{\kappa_{0}}\left|{\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}\right|=1$

This condition applies when the sum of the square roots of eigenvalues of the
Laplacian $L$ is the original parametric forcing frequency $\omega$, since
$\left|{\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}\right|=\frac{1}{\sqrt{\kappa_{0}}}=\omega$.

By a similar analysis, the relevant secular terms are obtained by combining
the coefficients of

  * •

$\exp\left(j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}t\right)=\exp\left(j\left(-\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}+1\right)t\right)$

  * •

$\exp\left(-j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}t\right)=\exp\left(j\left(\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}-1\right)t\right)$,

which are

$\displaystyle\left[c_{l}^{\prime}(\tau)-\tfrac{j}{2}\tfrac{\kappa_{1}}{\kappa_{0}}\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}c_{l}(\tau)-\tfrac{j}{4}\tfrac{1}{\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}}\mathbf{v}_{l}^{\star}E_{ik}\mathbf{v}_{m}d_{m}(\tau)\right]$
$\displaystyle\times\text{-}2{\scriptstyle\Omega}_{l}je^{\left(j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}t\right)}$
(76)
$\displaystyle\left[d_{m}^{\prime}(\tau)+\tfrac{j}{2}\tfrac{\kappa_{1}}{\kappa_{0}}\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}d_{m}(\tau)+\tfrac{j}{4}\tfrac{1}{\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}}\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}c_{l}(\tau)\right]$
$\displaystyle\times\text{-}2{\scriptstyle\Omega}_{m}je^{\left(j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}t\right)}$
(77)

The complex conjugate terms combine similarly:

  * •

$\exp\left(-j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}t\right)=\exp\left(j\left(\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}-1\right)t\right)$

  * •

$\exp\left(j\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}t\right)=\exp\left(j\left(-\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}+1\right)t\right)$,

giving us complex conjugates of terms (76) and (77), but involving
$c_{m}(\tau)$ and $d_{l}(\tau)$.

We pick $c_{l}(\tau)$ and $d_{m}(\tau)$ so that these terms vanish, which
gives us the system of equations

$\displaystyle\begin{bmatrix}c_{l}^{\prime}(\tau)\\\
d_{m}^{\prime}(\tau)\end{bmatrix}$
$\displaystyle=\frac{j}{2\sqrt{\kappa_{0}}}\underbrace{\begin{bmatrix}\kappa_{1}{\scriptstyle\Omega}_{l}&\frac{\kappa_{0}(\mathbf{v}_{l}^{\star}E_{ik}\mathbf{v}_{m})}{2{\scriptstyle\Omega}_{l}}\\\
-\frac{\kappa_{0}(\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l})}{2{\scriptstyle\Omega}_{m}}&-\kappa_{1}{\scriptstyle\Omega}_{m}\end{bmatrix}}_{A}\begin{bmatrix}c_{l}(\tau)\\\
d_{m}(\tau)\end{bmatrix}$ (78) $\displaystyle\text{or
}\mathbf{c}^{\prime}(\tau)$ $\displaystyle=\frac{j}{2\sqrt{\kappa_{0}}}A\
\mathbf{c}(\tau).$

As before, this choice of $\mathbf{c}(\tau)$ makes $\Phi_{1}(t,\tau)$ stable
in $t$. For bounded slow-time behavior of $\Phi_{0}(t,\tau)$, we require
$\mathbf{c}(\tau)$ to be stable, or the eigenvalues of $A$ to be real (because
of the $j$ prefactor). For the discriminant of the characteristic polynomial
of $A$ to be non-negative, we require

$\displaystyle\left({\scriptstyle\Omega}_{l}-{\scriptstyle\Omega}_{m}\right)^{2}\kappa_{1}^{2}+4{\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\kappa_{1}^{2}$
$\displaystyle-\left(\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right)^{2}\frac{\kappa_{0}^{2}}{{\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}}\geq
0$ $\displaystyle\implies$
$\displaystyle\left({\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}\right)^{2}\kappa_{1}^{2}-\left(\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right)^{2}\frac{\kappa_{0}^{2}}{{\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}}\geq
0$ $\displaystyle\implies$
$\displaystyle\kappa_{1}^{2}\geq\left(\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right)^{2}\frac{\kappa_{0}^{2}}{\left({\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}\right)^{2}{\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}}$

Since
${\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}=\frac{1}{\sqrt{\kappa_{0}}}$,
this stability condition can be simplified further, to give

$\left|\kappa_{1}\right|\geq\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|\frac{\kappa_{0}^{3/2}}{\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(79)

This suggests that there is a critical slope $\kappa_{1}$ below which the
solution at $(\kappa_{0}+\kappa_{1}\epsilon,\epsilon)$ is unstable. In figure
11, this corresponds to one of the two test curves (red and cyan) originating
at a point of instability, with the critical slope $\kappa_{1}$ determining
stability as $\epsilon\to 0$.

####  A-A3 Interpreting the multiple-scale analysis stability criteria

Our findings from the first-order multiple scale perturbation analysis are
presented in table I.

Condition on $\kappa_{0}$ | Condition on $\kappa_{1}$ | Stability
---|---|---
$\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}\pm\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}\neq\pm 1$ | - | Always stable
$\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}-\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}=\pm 1$ | - | Always stable
$\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{l}+\sqrt{\kappa_{0}}{\scriptstyle\Omega}_{m}=1$ | $\left|\kappa_{1}\right|\geq\frac{\kappa_{0}^{3/2}\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|}{\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$ | Stable
| $\left|\kappa_{1}\right|<\frac{\kappa_{0}^{3/2}\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|}{\left(\Omega_{l}\Omega_{m}\right)^{1/2}}$ | Unstable
TABLE I: Conditions for stability of the system (53), which represents
linearized swing dynamics on a graph where one edge weight is varied
periodically. Here, $\kappa$ lies on a test curve
$\kappa=\kappa_{0}+\epsilon\kappa_{1}+\mathcal{O}(\epsilon^{2})$, and we
impose conditions (derived in Section A-A) on $\kappa_{0}$ and $\kappa_{1}$
for stability at low forcing amplitudes $\epsilon$ (_i.e._ as $\epsilon\to
0$). In this limit, the system is unstable at a discrete set of parameters
$\kappa_{0}$ that depends on the eigenvalues and eigenvectors of the graph
Laplacian.

Since $\kappa_{1}$ represents the slope of the curve we are following in the
$(\kappa,\epsilon)$ plane at $(0,\kappa_{0})$, the critical value of
$\kappa_{1}$ corresponds to the tongue boundary (when the tongue exists). We
can reinterpret this in the $(\omega,\epsilon)$ plane via the scaling rule
$\kappa=1/\omega^{2}$. The curve
$\kappa(\epsilon)=\kappa_{0}+\epsilon\kappa_{1}+\mathcal{O}(\epsilon^{2})$
transforms to
$\omega(\epsilon)=\omega_{0}+a\epsilon+\mathcal{O}(\epsilon^{2})$, where
$\kappa_{0}=\frac{1}{\omega_{0}^{2}}$ and $a$ is
$\mathrm{d}\omega(\epsilon)/\mathrm{d}\epsilon|_{\epsilon=0}$. Then

$\displaystyle\kappa_{1}$
$\displaystyle:=\left.\frac{\mathrm{d}\kappa(\epsilon)}{\mathrm{d}\epsilon}\middle|_{\epsilon=0}\right.$
$\displaystyle=\left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}\frac{1}{\omega(\epsilon)^{2}}\middle|_{\epsilon=0}\right.$
$\displaystyle=-\left[\frac{2}{\omega(\epsilon)^{3}}\frac{\mathrm{d}\omega(\epsilon)}{\mathrm{d}\epsilon}\right]_{\epsilon=0}=-\frac{2}{\omega_{0}^{3}}\left.\frac{\mathrm{d}\omega(\epsilon)}{\mathrm{d}\epsilon}\middle|_{\epsilon=0}\right.$
$\displaystyle=-\frac{2}{\omega_{0}^{3}}a$ (80)

The stability condition (79) on $\kappa_{1}$ can then be expressed in terms of
$\omega$, since $\kappa_{0}=1/\omega_{0}^{2}$:

$\displaystyle\left|\kappa_{1}\right|$
$\displaystyle\geq\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|\frac{\kappa_{0}^{3/2}}{\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
$\displaystyle\implies\left|-\frac{2}{\omega_{0}^{3}}a\right|$
$\displaystyle\geq\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|\frac{1}{\omega_{0}^{3}\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
$\displaystyle\implies\left|a\right|$
$\displaystyle\geq\frac{\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(81)

From (80) and (81), we can state the _instability_ criteria for solutions of
the system (12) for $\epsilon\to 0$ at a point on the curve
$\omega(\epsilon)=\omega_{0}+a\epsilon+\mathcal{O}(\epsilon^{2})$ as

$\displaystyle\omega_{0}={\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m},\
l,m\in 1,2,\dots,n$
$\displaystyle\left|a\right|<\frac{\left|\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}\right|}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(82)

###  A-B The multiple scale method with eigenvalue multiplicity

The multiscale analysis of appendix A-A assumes that all eigenvalues of the
graph Laplacian have multiplicity $1$. This is not true for graphs with
symmetries, such as the ring networks considered in Section IV. The vector-
valued multi-scale analysis is still valid, but now Equations (67) and (68)
apply to any pair of eigenvectors $\mathbf{v}_{m}$ and $\mathbf{v}_{l}$ of the
Laplacian $L$ corresponding to the eigenvalues ${\scriptstyle\Omega}_{m}^{2}$
and ${\scriptstyle\Omega}_{l}^{2}$ respectively. Let $V_{m}$ and $V_{l}$ be
matrices whose columns form an orthonormal basis for the eigenspaces
corresponding to ${\scriptstyle\Omega}_{m}^{2}$ and
${\scriptstyle\Omega}_{l}^{2}$ respectively. Then we have

$\displaystyle\mathbf{v}_{m}=V_{m}\mathbf{z}_{m},\quad\mathbf{v}_{l}=V_{l}\mathbf{z}_{l}$

for appropriately sized vectors $\mathbf{z}_{m}$ and $\mathbf{z}_{l}$.
Further, $\left\|\mathbf{z}_{m}\right\|=1$ since
$\left\|\mathbf{v}_{m}\right\|=1$ and $V_{m}$ is unitary (similarly
$\mathbf{z}_{l}$). Consequently the instability criterion (30) is modified to

$\displaystyle\omega_{0}={\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m},\
l,m\in 1,2,\dots,n$
$\displaystyle\left|\omega^{\prime}(0)\right|<\max_{\left\|\mathbf{z}_{m}\right\|\text{=}\left\|\mathbf{z}_{l}\right\|\text{=}1}\frac{\left|\mathbf{z}_{m}^{\star}\left(V_{m}^{\star}E_{ik}V_{l}\right)\mathbf{z}_{l}\right|}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}.$
(83)

This maximum is attained when $\mathbf{z}_{m}$ and $\mathbf{z}_{l}$ are the
first left and right singular vectors of $V_{m}^{\star}E_{ik}V_{l}$
respectively, and the maximum is the $2$-induced norm of this matrix:

$\displaystyle\left|\omega^{\prime}(0)\right|<\frac{\left\|V_{m}^{\star}E_{ik}V_{l}\right\|_{2}}{2\left({\scriptstyle\Omega}_{l}{\scriptstyle\Omega}_{m}\right)^{1/2}}$
(84)

###  A-C Multiple scale perturbation analysis: unstable modes

In this section we verify that the modes of the unstable response at a
critical frequency
$\omega_{lm}=\left|{\scriptstyle\Omega}_{l}+{\scriptstyle\Omega}_{m}\right|$
(Eq. 82) are the corresponding eigenvectors $\mathbf{v}_{l}$ and
$\mathbf{v}_{m}$ of the graph Laplacian. From Eq. 61, the $\mathcal{O}(1)$
term in the expansion of the response is given by

$\displaystyle\Phi_{0}(t,\tau)=V\left(e^{j\sqrt{\kappa_{0}}\Omega
t}\mathbf{c}(\tau)+e^{\text{-}j\sqrt{\kappa_{0}}\Omega
t}\mathbf{d}(\tau)\right).$

At a critical frequency $\omega_{lm}$, $\mathbf{c}(\tau)$ and
$\mathbf{d}(\tau)$ satisfy (from Eq. 65 and Eq. 78)

$\displaystyle\begin{bmatrix}c_{k}(\tau)\\\ d_{k}(\tau)\end{bmatrix}$
$\displaystyle=\exp\left(j\frac{\kappa_{1}}{2\kappa_{0}^{1/2}}\begin{pmatrix}1&0\\\
0&\text{-}1\end{pmatrix}\tau\right)\begin{bmatrix}c_{k}(0)\\\
d_{k}(0)\end{bmatrix}$ $\displaystyle k\neq l,m$
$\displaystyle\begin{bmatrix}c_{l}(\tau)\\\ d_{m}(\tau)\end{bmatrix}$
$\displaystyle=\exp\left(j\frac{1}{2\sqrt{\kappa_{0}}}A\tau\right)\begin{bmatrix}c_{l}(0)\\\
d_{m}(0)\end{bmatrix}$ $\displaystyle\begin{bmatrix}c_{m}(\tau)\\\
d_{l}(\tau)\end{bmatrix}$
$\displaystyle=\exp\left(\text{-}j\frac{1}{2\sqrt{\kappa_{0}}}A\tau\right)\begin{bmatrix}c_{m}(0)\\\
d_{l}(0)\end{bmatrix}$

where $A$ is as in Eq. 78 and has eigenvalues with non-zero imaginary part.
$c_{l}$, $c_{m}$, $d_{l}$ and $d_{m}$ are the only components of $\mathbf{c}$
and $\mathbf{d}$ that grow unboundedly. $\Phi_{0}(t,\tau)$ is

$\displaystyle\Phi_{0}(t,\tau)$
$\displaystyle=\mathbf{v}_{l}(e^{j\sqrt{\kappa_{0}}\Omega_{l}t}c_{l}(\tau)+e^{\text{-}j\sqrt{\kappa_{0}}\Omega_{l}t}d_{l}(\tau))$
$\displaystyle+\mathbf{v}_{m}(e^{j\sqrt{\kappa_{0}}\Omega_{m}t}c_{m}(\tau)+e^{\text{-}j\sqrt{\kappa_{0}}\Omega_{m}t}d_{m}(\tau))$
$\displaystyle+\text{stable modes.}$

This choice of $\mathbf{c}$ and $\mathbf{d}$ makes the secular terms in the
$\mathcal{O}(\epsilon)$ dynamics (Eq. 62) vanish, so there are no other
unstable modes to $\mathcal{O}(\epsilon)$. The unstable response at
$\omega_{lm}$ as predicted by the multiple scale analysis thus involves only
the modes $\mathbf{v}_{l}$ and $\mathbf{v}_{m}$.

###  A-D Eigenvector test for controllability

We obtain a controllability criterion for the second order system (31) by
applying the eigenvector test [17] to it rewritten in first order form.

######  Lemma 1 (Eigenvector test for harmonic oscillator networks)

Let $L\in\mathbb{R}^{n\times n}$ be real-symmetric with
$L\mathbf{v}_{m}={\scriptstyle\Omega}_{m}^{2}\mathbf{v}_{m}$,
${\scriptstyle\Omega}_{m}\neq 0$. Let $\mathcal{B}\in\mathbb{C}^{n\times k}$.
The harmonic oscillator system

$\displaystyle\ddot{\phi}=-L\phi+\mathcal{B}u(t)$ (85)

is not controllable from the input $u(t)$ if
$\mathbf{v}_{m}^{\star}\mathcal{B}=0$ for some $m$.

###### Proof:

Rewritten in first order form, Eq. 85 is

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix}\phi\\\
\dot{\phi}\end{bmatrix}=\underbrace{\begin{bmatrix}0&I\\\
\text{-}L&0\end{bmatrix}}_{\mathcal{L}}\begin{bmatrix}\phi\\\
\dot{\phi}\end{bmatrix}+\begin{bmatrix}0\\\ \mathcal{B}\end{bmatrix}u(t)$ (86)

For ${\scriptstyle\Omega}_{m}\neq 0$, the eigenvalues and eigenvectors of the
generator $\mathcal{L}$ are given by the relation

$\displaystyle\mathcal{L}\begin{pmatrix}\mathbf{v}_{m}\\\ \pm
j{\scriptstyle\Omega}_{m}\mathbf{v}_{m}\end{pmatrix}=\pm
j{\scriptstyle\Omega}_{m}\begin{pmatrix}\mathbf{v}_{m}\\\ \pm
j{\scriptstyle\Omega}_{m}\mathbf{v}_{m}\end{pmatrix},$

which can be verified by direct computation.

Applying the eigenvector test to Eq. 86, a left eigenvector of $\mathcal{L}$
is in the left nullspace of
$\begin{bmatrix}0&\mathcal{B}\end{bmatrix}^{\mathrm{T}}$ iff

$\displaystyle\begin{bmatrix}\mathbf{v}_{m}^{\star}&\mp
j{\scriptstyle\Omega}_{m}^{\text{-}\star}\mathbf{v}_{m}^{\star}\end{bmatrix}\begin{bmatrix}0\\\
\mathcal{B}\end{bmatrix}=0$
$\displaystyle\implies\mathbf{v}_{m}^{\star}\mathcal{B}=0$

∎

Applying Lemma 1 to system (31), it is uncontrollable from the input
$u_{l}(t)=e^{j(\Omega_{l}+\omega)t}c_{l}$ iff

$\displaystyle\begin{bmatrix}\mathbf{v}_{m}^{\star}&\text{-}j{\scriptstyle\Omega}_{m}^{\text{-}1}\mathbf{v}_{m}^{\star}\end{bmatrix}\underbrace{\begin{bmatrix}0&0\\\
E_{ik}&0\end{bmatrix}\begin{bmatrix}\mathbf{v}_{l}\\\
j{\scriptstyle\Omega}_{l}\mathbf{v}_{l}\end{bmatrix}}_{\mathcal{B}}$
$\displaystyle=0$
$\displaystyle\implies\mathbf{v}_{m}^{\star}E_{ik}\mathbf{v}_{l}$
$\displaystyle=0.$ (87)

For ${\scriptstyle\Omega}_{m}=0$, $\mathcal{L}$ has a Jordan block of size
$2k$, where $L$ has $k$ zero eigenvalues, and no conclusion can be drawn for
general $\mathcal{B}$. But when $L$ and $E_{ik}$ are symmetric graph
Laplacians, as in Eq. 31, the mode corresponding to the zero eigenvalue of $L$
is $\mathbf{v}_{m}=\mathbbm{1}$. This mode is not controllable from $u_{l}(t)$
either, as evident from the PBH test [17]:

$\displaystyle\mathrm{rank}\left[\begin{array}[]{c:c}\lambda
I-\mathcal{L}&\begin{matrix}0\\\
E_{ik}\mathbf{v}_{l}\end{matrix}\end{array}\right]$ $\displaystyle=$
$\displaystyle\mathrm{rank}\left[\begin{array}[]{cc:c}\lambda I&\text{-}I&0\\\
L&\lambda
I&E_{ik}\mathbf{v}_{l}\end{array}\right]\stackrel{{\scriptstyle?}}{{<}}n$

When $\lambda=0$, this matrix drops rank, since $\mathbbm{1}^{\star}E_{ik}=0$,
and

$\displaystyle\begin{bmatrix}0&\mathbbm{1}^{\star}\end{bmatrix}\left[\begin{array}[]{cc:c}0&\text{-}I&0\\\
L&0&E_{ik}\mathbf{v}_{l}\end{array}\right]=0$