# Complete Synchronization of Coupled Oscillators Based on the Contraction
Theory

Brian Y Zhang

###### Abstract

In this paper, the contraction theory is introduced and applied to the
complete synchronization of complex networks of coupled oscillators. Based on
the contraction theory, the conditions that make the virtual systems of the
coupled-oscillator networks in three different constructions, including all-
to-all topology, star topology, and ring topology, contracting are given,
under which the coupled-oscillator networks are completely synchronized
regardless of the initial conditions. As concrete examples, the conditions
that the networks consisting of linearly coupled van der Pol and Duffing van
der Pol oscillators are completely synchronized in any initial conditions are
obtained. It is found that the virtual system of the network of linearly
coupled Rayleigh van der Pol oscillators cannot be contracting for any
coefficients, but its contraction region can be found. Using the trapping
region method, the conditions for the networks of linearly coupled Rayleigh
van der Pol oscillators with different coefficients that start in some
specific regions to be completely synchronized are obtained.

Department of Physics, University of Colorado Denver

Campus Box 157, PO Box 173364

Denver, CO 80217-3364, USA

Keywords: nonlinear coupled oscillators; synchronization; contraction theory;
Rayleigh van der Pol oscillators

##  1 Introduction

The study of coupled oscillators and complex systems appears in many fields of
research, such as mathematics, engineering, robotics, and biology [1-4]. A
large number of complex systems that exist in nature and in the real world can
be described by networks. The nodes in a network are used to represent the
individual oscillators in a complex system, and the links in a network are
used to represent the couplings between oscillators in a complex system. The
development of network science has greatly promoted the in-depth study and
development of complex systems. The synchronization phenomenon is one of the
hot topics of interest in the study of coupled oscillators and network
science. The synchronization is crucial to understand the physical mechanisms
of many seemingly disparate or unrelated collective behaviors of complex
systems [5-6]. Early synchronization phenomena were discovered and observed in
a variety of different man-made devices. With the development of scientific
research, there have been many researchers studying and discussing
synchronization phenomena and their related problems. It has been found that
synchronization is widely present in various naturally occurring collective
phenomena in nature, such as Josephson junctions, nanomechanics, synchronous
flickering of fireflies, collective chirping of crickets, recoil atomic
lasers, etc [7-11]. Therefore, it is important to explore the emergence of
these self-organized cooperative states and reveal their intrinsic dynamical
mechanisms not only for understanding the group dynamics of complex systems,
but also for conducting related experiments and potential applications
[12-13].

The phenomenon that both phase and amplitude are synchronized is called
complete synchronization. Contraction theory, which originated form fluid
dynamics and differential geometry, has been established as an effective tool
for analyzing the complete synchronization behaviors of nonlinear networks
with given state spaces [14-15]. By constructing a virtual system of the
network and obtaining its contraction region, it is possible to theoretically
predict whether the network will be completely synchronized. Sometimes we
cannot construct a virtual system of a coupled oscillator network, but if we
can build its “virtual network” and thus construct a higher dimensional
virtual system, then we can still find the conditions to synchronize the
original network. A key problem in applying contraction theory is how to
confirm that the trajectories of all oscillators in a network are confined to
the contraction region of the corresponding virtual system. For some
particular coupled-osillator networks, we can adjust their coefficients so
that the contraction regions cover the entire phase plane, which means that
every trajectory remains within the contraction regions no matter where it
starts. However, for most networks, the contraction regions are finite in the
phase plane regardless of coefficients. A feasible solution is to find the
trapping region of the virtual system and make sure that the trapping region
is contained in the contraction region, then as long as the trajectories of
all oscillators start inside the trapping region, they will remain there and
thus remain in the contraction region. Nevertheless, the virtual system of a
coupled-oscillator network usually depends on the trajectory of every
oscillator in the original network, so finding the trapping region requires
some analysis of the properties of the original coupled-oscillator system.
This will be discussed in detail in Section 4 using linearly coupled Rayleigh
van der Pol oscillators.

The paper is organized as follows. In Section 2, I introduce basic contraction
theory. In Section 3, I describe the application of contraction theory to the
complete synchronization of the coupled-oscillator networks in different
network constructions, including all-to-all topology, star topology, and ring
topology. In Section 4, I study the dynamics of linearly coupled van der Pol,
Duffing van der Pol, and Rayleigh van der Pol oscillators and give the
conditions for each of them to achieve complete synchronization.

##  2 Methods

Definition 2.1 (contracting region). Given a system equation of the form

$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x}(t)),$ (1)

where $\boldsymbol{x}\in\mathbb{R}^{n}$ is a set of state variables and
$\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a smooth nonlinear vector
function. A contraction region $C$ of this system is a region of the state
space where the symmetric part of the Jacobian $\left(\frac{\partial
f}{\partial\boldsymbol{x}}\right)_{s}$ is uniformly negative definite, i.e.,
$C\text{:=}\left\\{\boldsymbol{x}\in\mathbb{R}^{n}:\lambda_{\max}(\boldsymbol{x})<0\right\\}$,
where $\lambda_{\max}(\boldsymbol{x})$ is the largest eigenvalue of
$\frac{1}{2}\left(\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}+\left(\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}\right)^{T}\right)$.

######  Theorem 2.1.

If any two trajectories, $\boldsymbol{x}_{1}(t)$ and $\boldsymbol{x}_{2}(t)$,
of a system of the form (1), starting from different initial conditions,
remain within $C$, then they converge exponentially to each other.

###### Proof.

See Appendix A.

∎

In order to apply this theorem to synchronization of coupled oscillators, the
concept of a virtual system needs to be introduced [16]. Such a system is
expressed in the old state variables of the original network of the coupled
oscillators and in a set of new virtual state variables. The virtual system
can recover the trajectory of each node in the original network by
substituting the state variables of each node for the virtual state variables.
Thus, by constructing a virtual system, the multiple coupled systems are
“condensed” into a single system as described in Eq. (1) to which the
contraction theory can be applied.

##  3 Problem Statement

Consider a network of nonlinear systems

$\dot{\boldsymbol{x}}_{i}=\boldsymbol{f}_{i}\left(\boldsymbol{x}_{i}\right)+\boldsymbol{u}_{i},i=1,2,\cdots,N,$
(2)

where
$\boldsymbol{x}_{i}=\left({x_{i}}_{1}\;{x_{i}}_{2}\;\cdots\;{x_{i}}_{n}\right){}^{T}\in\mathbb{R}^{n}$
is the state variables of node $i$,
$\boldsymbol{f}_{i}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ are nonlinear smooth
vector functions, $\boldsymbol{u}_{i}$ are coupling functions which depend on
the difference between the state variables of $i^{\text{th}}$ node and those
of others that are linked with it.

Definition 3.1 (complete synchronization). The network described by (2) is
said to be completely synchronized if $\forall j=1,\cdots,n$,
$\exists\boldsymbol{x}_{i}(0)\in\mathbb{R}^{n}$ for $i=1,\cdots,N$ such that

$\lim_{t\to\infty}x_{1j}(t)=\lim_{t\to\infty}x_{2j}(t)=\cdots=\lim_{t\to\infty}x_{Nj}(t)\neq
0.$

######  Theorem 3.1.

Consider a network of the form (2). Assume that there exists a virtual system
of this network:

$\dot{\boldsymbol{y}}=\boldsymbol{\Phi}\left(\boldsymbol{y},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right),$
(3)

where $\boldsymbol{y}\in\mathbb{R}^{n}$ is the virtual state variable, such
that

$\begin{split}\Phi\left(\boldsymbol{x}_{1},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right)&=\boldsymbol{f}_{1}\left(\boldsymbol{x}_{1}\right)+\boldsymbol{u}_{1},\\\
\Phi\left(\boldsymbol{x}_{2},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right)&=\boldsymbol{f}_{2}\left(\boldsymbol{x}_{2}\right)+\boldsymbol{u}_{2},\\\
&\vdots\\\
\Phi\left(\boldsymbol{x}_{N},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right)&=\boldsymbol{f}_{N}\left(\boldsymbol{x}_{N}\right)+\boldsymbol{u}_{N}.\end{split}$

Then if the trajectories of the oscillators,
$\boldsymbol{x}_{1}(t),\boldsymbol{x}_{2}(t),\text{...},\boldsymbol{x}_{N}(t)$,
remain within the contraction region of the virtual system (3) for some
initial conditions, they are completely synchronized.

###### Proof.

Notice that all trajectories of network (2) are just particular solutions to
the virtual $\boldsymbol{y}$-system as described in Eq. (3). Since the
trajectories of (2) all remain within the contraction region of (3) as the
system evolves, Theorem 2.1 implies that all oscillators in (2) converge into
the same synchronous state exponentially. ∎

Theorem 3.1 shows that constructing virtual systems is the key to studying
comlete synchronization problems based on contraction theory. The virtual
system does not have to be unique, we can choose any one of them that is easy
to obtain the contraction region. However, not all networks have corresponding
virtual systems of the form (3). To apply contraction theory, we need to make
some assumptions about the network (2). The two most basic assumptions are
that all coupled oscillators in the network are identical, i.e.,
$\boldsymbol{f}_{1}=\boldsymbol{f}_{2}=\cdots=\boldsymbol{f}_{N}$, and that
the couplings are linear, i.e., $\boldsymbol{u}_{i}$ are linear
transformations from $\mathbb{R}^{n}$ to $\mathbb{R}^{n}$ for all $i$ from 1
to $N$. Note that a linear transform
$\boldsymbol{u}_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ must be left-
multiplication by some $n\times n$ matrix $A_{i}$. For a two-coupled-
oscillator system, these two assumptions are sufficient to guarantee the
existence of its virtual system.

Below, I will use the following symbols and notations:

“$M\prec\ 0$ ($\succ 0$)” means that the matrix $M$ is negative (positive)
definite;

“$M\preceq 0$ ($\succeq 0$)” means that $M$ is negative (positive) semi-
definite;

$J_{\boldsymbol{f}}$ refers to the Jacobian of $\boldsymbol{f}$, i.e.,
$J=\frac{\partial\boldsymbol{f}(\boldsymbol{y})}{\partial{\boldsymbol{y}}}$;

and subscript $s$ indicates the symmetric part of a matrix.

######  Proposition 3.2.

For two identical coupled oscillators with linear couplings described by

$\begin{split}&\dot{\boldsymbol{x}}_{1}=\boldsymbol{f}\left(\boldsymbol{x}_{1}\right)+A_{1}(\boldsymbol{x}_{2}-\boldsymbol{x}_{1})\\\
&\dot{\boldsymbol{x}}_{2}=\boldsymbol{f}\left(\boldsymbol{x}_{2}\right)+A_{2}(\boldsymbol{x}_{1}-\boldsymbol{x}_{2}),\end{split}$
(4)

the complete synchronization arises if $\boldsymbol{x}_{1}(t)$ and
$\boldsymbol{x}_{2}(t)$ are always inside the region where
${J_{\boldsymbol{f}}}_{s}-(A_{1}+A_{2})_{s}\prec 0$.

###### Proof.

Notice that the virtual system of (12) can be constrcted as follows

$\dot{\boldsymbol{y}}=f\left(\boldsymbol{y}\right)+A_{1}\boldsymbol{x}_{2}+A_{2}\boldsymbol{x}_{1}-(A_{1}+A_{2})\boldsymbol{y}.$
(5)

This is because substituting $\boldsymbol{x}_{1}$ for $\boldsymbol{y}$ in Eq.
(4) yields the first equation of Eq. (4) and, similarly, inserting
$\boldsymbol{x}_{2}$ into Eq. (4) gives the second equation of Eq. (4).

Write the virtual system (4) as
$\dot{\boldsymbol{y}}=\boldsymbol{\Phi}\left(\boldsymbol{y},\boldsymbol{x}_{1},\boldsymbol{x}_{2}\right)$.
We have

$\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{y}}=J_{\boldsymbol{f}}-(A_{1}+A_{2}).$

By Definition 2.1, the contractoin region of (5 is the region in which
$(\frac{\partial\boldsymbol{\Phi}}{\partial\boldsymbol{y}})_{s}={J_{\boldsymbol{f}}}_{s}-(A_{1}+A_{2})_{s}\prec
0$. Therefore, by Theorem 3.1, if $\boldsymbol{x}_{1}(t)$ and
$\boldsymbol{x}_{2}(t)$ remain within this region, (4) is completely
synchronized. ∎

For a non-identical two-coupled-oscillator system with linear couplings, one
way to construct its virtual system is

$\dot{\boldsymbol{y}}=(\boldsymbol{y}-\boldsymbol{x}_{2})\frac{\boldsymbol{f}_{1}\left(\boldsymbol{y}\right)}{\boldsymbol{x}_{1}-\boldsymbol{x}_{2}}+(\boldsymbol{y}-\boldsymbol{x}_{1})\frac{\boldsymbol{f}_{2}\left(\boldsymbol{y}\right)}{(\boldsymbol{x}_{2}-\boldsymbol{x}_{1})}+A_{1}\boldsymbol{x}_{2}+A_{2}\boldsymbol{x}_{1}-(A_{1}+A_{2})\boldsymbol{y}.$

However, this system forms singularities when the trajectories of the two
oscillators of the original system coincide, and the right-hand-side function
is not smooth in this case. As a result, the virtual system so constructed
does not make sense. Therefore, we cannot apply the contraction theory to
coupled non-identical oscillators. This is a drawback of applying contraction
theory to synchronization problems, because in reality many coupled
oscillators are not exactly identical. In fact, the coupled non-identical
oscillators cannot generally completely synchronized at an exponential rate.
Another fact is that in the case of nonlinear couplings, a virtual system of
the two-coupled-oscillator system can be easily constructed if the coupling
between oscillators is the difference of some nonlinear function of their
state variables instead of a nonlinear function of the difference of their
state variables [17]. We will not discuss nonlinear coupling in the following
paper.

To construct the virtual system of a network containing more than two
oscillators, in addition to having identical oscillators and linear couplings,
the network must possess some particular structure and the coupling must have
some sort of symmetry. In the following we discuss three common types of
network topology - all-to-all topology, star topology, and ring topology.

######  Proposition 3.3.

Consider $N$ all-to-all coupled identical oscillators with symmetric linear
couplings described by

$\dot{\boldsymbol{x}}_{i}=\boldsymbol{f}\left(\boldsymbol{x}_{i}\right)+\sum_{j=1}^{N}A_{j}(\boldsymbol{x}_{j}-\boldsymbol{x}_{i}),\forall
i=1,\cdots,N,$ (6)

Here, the “symmetric coupling” means that each oscillator has the same effect
on all other oscillators. This network is completely synchronized if
$\boldsymbol{x}_{i}(t)$ are always inside the region where
${J_{\boldsymbol{f}}}_{s}-\sum_{j=1}^{N}{A_{i}}_{s}\prec 0$ for all $i$.

###### Proof.

Notice that this coupled-oscillator system has a virtual system

$\dot{\boldsymbol{y}}=\boldsymbol{f}\left(\boldsymbol{y}\right)+\sum_{j=1}^{N}A_{j}\boldsymbol{x}_{j}-(\sum_{j=1}^{N}A_{j})\boldsymbol{y}.$
(7)

The rest of the proof follows by similar arguments as in Proposition 3.2. ∎

######  Proposition 3.4.

Consider $N$ identical oscillators coupled by a star network with symmetric
linear couplings described by

$\displaystyle\dot{\boldsymbol{x}}_{1}$
$\displaystyle=\boldsymbol{f}\left(\boldsymbol{x}_{1}\right)+\sum_{j=2}^{N}A_{j}(\boldsymbol{x}_{j}-\boldsymbol{x}_{1}),$
(8a) $\displaystyle\dot{\boldsymbol{x}}_{i}$
$\displaystyle=\boldsymbol{f}\left(\boldsymbol{x}_{i}\right)+A_{1}(\boldsymbol{x}_{1}-\boldsymbol{x}_{i}),\text{
}\forall i=2,\cdots,N,$ (8b)

This network is completely synchronized if $\boldsymbol{x}_{1}(t)$ and
$\boldsymbol{x}_{2}(t)$ are always inside the region where
${J_{\boldsymbol{f}}}_{s}-\sum_{j=1}^{N}{A_{j}}_{s}\prec 0$ and
$\boldsymbol{x}_{i}(t)$ are always inside the region where
${J_{\boldsymbol{f}}}_{s}-{A_{1}}_{s}\prec 0$ for all $i$ from 2 to $N$.

###### Proof.

Notice that the network (8) does not have a single virtual system, but it is
easy to see that the virtual system of (8b) is

$\dot{\boldsymbol{y}}=\boldsymbol{f}\left(\boldsymbol{y}\right)+A_{1}(\boldsymbol{x}_{1}-\boldsymbol{y})$
(9)

Since $\boldsymbol{x}_{i}(t)$ remain in the contraction region of (9) where
${J_{\boldsymbol{f}}}_{s}-{A_{1}}_{s}\prec 0$ at all times for all $i$ from 2
to $N$, by Theorem 3.1, all the oscillators except the first one converge
toward the same synchronous state exponentially. Hence after a while Eq. (8a)
reduces to

$\dot{\boldsymbol{x}}_{1}=\boldsymbol{f}\left(\boldsymbol{x}_{1}\right)+\sum_{j=2}^{N}A_{j}(\boldsymbol{x}_{2}-\boldsymbol{x}_{1}),$

which describes the state of the first oscillator in the network. Then the
virtual system of the subnetwork consisting of only the first two oscillators
is

$\dot{\boldsymbol{y}}=\boldsymbol{f}\left(\boldsymbol{y}\right)+A_{1}\boldsymbol{x}_{1}+(\sum_{j=2}^{N}A_{j})\boldsymbol{x}_{2}-(\sum_{j=1}^{N}A_{j})\boldsymbol{y}.$
(10)

Therefore, the first two oscillators are completely synchronized as
$\boldsymbol{x}_{1}(t)$ and $\boldsymbol{x}_{2}(t)$ are always inside the
contraction region of (10) where ${J_{f}}_{s}-\sum_{j=1}^{N}{A_{j}}_{s}\prec
0$, and thus the whole newtork (8) is completely synchronized by transitivity.
∎

The limitation of applying the above propositions is that when the symmetric
part of the Jacobian of the virtual system of a network is not uniformly
negative definite on the entire state space $D$, it is difficult to determine
whether the trajectories of the oscillators in this network are confined in
the contraction region of the virtual system. Note that we cannot say that a
trajectory will stay in the contraction region at all future times simply
because it starts there. One way to determine whether the trajectory will
remain within the contraction region by the position of its starting point is
to find a trapping region, i.e, a compact subset $N$ of the state space such
that the flow of the system is inward everywhere on the boundary of $N$, that
is inside the contraction region. If we can find such a trapping region $N$,
then every trajectory that starts within $N$ will remain there for all future
time. These are discussed in detail in Section 4.

In general, we cannot find virtual systems of bidirectional ring networks in
which each oscillator affects its two neighboring oscillators and of
unidirectional ring networks in which each oscillator affects only its next
oscillator. We now introduce a multivariate “virtual network” as follows. For
a network (2), we try to construct

$\dot{\boldsymbol{y}_{i}}=\boldsymbol{\Phi}\left(\boldsymbol{y}_{i},\boldsymbol{y}_{i+1},\text{...},\boldsymbol{y}_{i+N-1},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right),i=1,\cdots,N,$

for virtual variables $\boldsymbol{y}_{1},\text{...},\boldsymbol{y}_{N}$,
which satisfies

$\boldsymbol{\Phi}\left(\boldsymbol{x}_{i},\boldsymbol{x}_{i+1},\text{...},\boldsymbol{x}_{i+N-1},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right)=\boldsymbol{f}\left(\boldsymbol{x}_{i}\right)+\boldsymbol{u}_{i},\forall
i=1,\cdots,N.$ (11)

where all subscripts are calculated modulo $N$. We can turn this virtual
network into a single $nN$-dimensinal virtual system of $\boldsymbol{Y}$ by
concatenating all original $n$-dimensional virtual state vectors
$\boldsymbol{y}_{i}$ to form an $nN$-dimensional virtual state variable
$\boldsymbol{Y}$, that is

$\boldsymbol{Y}=\left(\boldsymbol{y}_{1}\text{ }\boldsymbol{y}_{2}\text{
}\cdots\text{ }\boldsymbol{y}_{N}\right){}^{T}\in\mathbb{R}^{nN}.$

Then we have

$\dot{\boldsymbol{Y}}=\boldsymbol{\Psi}\left(\boldsymbol{Y},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right)\equiv\begin{pmatrix}\boldsymbol{\Phi}\left(\boldsymbol{y}_{1},\boldsymbol{y}_{2},\text{...},\boldsymbol{y}_{N},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right)\\\
\boldsymbol{\Phi}\left(\boldsymbol{y}_{2},\text{...},\boldsymbol{y}_{N},\boldsymbol{y}_{1},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right)\\\
\vdots\\\
\boldsymbol{\Phi}\left(\boldsymbol{y}_{N},\boldsymbol{y}_{1},\text{...},\boldsymbol{y}_{N-1},\boldsymbol{x}_{1},\boldsymbol{x}_{2},\text{...},\boldsymbol{x}_{N}\right)\end{pmatrix}$
(12)

######  Theorem 3.5.

If network (2) has a virtual system (12) which is contracting with respect to
$\boldsymbol{Y}$, then the network is completely synchronized regardless of
the initial conditions.

###### Proof.

Since the function $\boldsymbol{\Phi}$ satisfies (11), the trajectory
$\left(\boldsymbol{x}_{1}\text{ }\boldsymbol{x}_{2}\text{ }\cdots\text{
}\boldsymbol{x}_{N}\right){}^{T}$ formed by concatenating the trajectories of
all the oscillators in the original network (10) and all its cyclic
permutations are particular solutions of the system (12). In particular,
$\left(\boldsymbol{x}_{1}\text{ }\boldsymbol{x}_{2}\text{ }\cdots\text{
}\boldsymbol{x}_{N}\right){}^{T}$ and $\left(\boldsymbol{x}_{2}\text{
}\cdots\text{ }\boldsymbol{x}_{N}\text{ }\boldsymbol{x}_{1}\right){}^{T}$ are
two particular solutions of (12). Since the system is contracting, they
converge exponentially to each other, which implies that
$\boldsymbol{x}_{1}\sim\boldsymbol{x}_{2},\text{ }\cdots,\text{
}\boldsymbol{x}_{N-1}\sim\boldsymbol{x}_{N},\text{ and
}\boldsymbol{x}_{N}\sim\boldsymbol{x}_{1}$, where
$\boldsymbol{a}\sim\boldsymbol{b}$ denotes that $\boldsymbol{a}$ and
$\boldsymbol{b}$ are converged. Therefore, all oscillators converge into the
same synchronous state vector. ∎

We now consider ring networks consisting of $N$ identical coupled oscillators
whose links are either unidirectional with identical linear coupling or
bidirectional with different linear couplings. Such unidirectional and
bidirectional ring networks can be expressed as the following two equations
respectively.

$\dot{\boldsymbol{x}}_{i}=\boldsymbol{f}\left(\boldsymbol{x}_{i}\right)+A(\boldsymbol{x}_{i-1}-\boldsymbol{x}_{i}),\text{
}\forall i=1,\cdots,N.$ (13a)
$\dot{\boldsymbol{x}}_{i}=\boldsymbol{f}\left(\boldsymbol{x}_{i}\right)+A_{i+1,i}(\boldsymbol{x}_{i+1}-\boldsymbol{x}_{i})+A_{i-1,i}(\boldsymbol{x}_{i-1}-\boldsymbol{x}_{i}),\text{
}\forall i=1,\cdots,N,$ (13b)

where $A_{i,j}$ denotes coupling from the $i^{th}$ to $j^{th}$ node, and the
subscripts $i,j$ are calculated modulo $N$.

Wang and Slotine [16] discussed the conditions for synchronization of
unidirectional ring network when $A$ is a square symmetric matrix, i.e,
$A^{T}=A$, and of bidirectional ring network when the couplings are
interactional, i.e., $A_{i,j}=A_{j,i}$. In the following, I relax both of
these conditions.

######  Theorem 3.6.

The unidirectional ring network (13a) is completely synhcronized if
${J_{f}}_{s}+A_{s}\prec 0$ on $D$, $A_{s}\succ 0$, and
$4A_{s}-AA_{s}^{-1}A^{T}\succeq 0$.

###### Proof.

The virtual network of (13a) can be constructed as

$\dot{\boldsymbol{y}}_{i}=\boldsymbol{f}\left(\boldsymbol{y}_{i}\right)+A(\boldsymbol{y}_{i-1}-\boldsymbol{y}_{i}),\text{
}i=1,\cdots,N.$ (14)

It is easy to check Eq. (14) by replacing $\boldsymbol{y}_{j}$ with
$\boldsymbol{x}_{j}$ for all $j$. Write (14) in the form of (12) by
concatenating $\boldsymbol{y}_{i}$. Then the Jacobian of $\boldsymbol{\Psi}$
is

$\begin{split}\frac{\partial\boldsymbol{\Psi}}{\partial{\boldsymbol{Y}}}=&\left(\begin{array}[]{cccc}J_{\boldsymbol{f}1}&0&\cdots&0\\\
0&J_{\boldsymbol{f}2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&J_{\boldsymbol{f}N}\end{array}\right)\\\
-&\left(\begin{array}[]{ccccc}A&0&\cdots&0&-A\\\ -A&A&0&\cdots&0\\\
\vdots&\vdots&\ddots&&\vdots\\\ \vdots&\vdots&&\ddots&\vdots\\\
0&\cdots&0&-A&A\end{array}\right),\end{split}$

where
$J_{\boldsymbol{f}i}=\frac{\partial\boldsymbol{f}(\boldsymbol{y}_{i})}{\partial{\boldsymbol{y}_{i}}}$.
The symmetric part of the Jacobian is

$\begin{split}(\frac{\partial\boldsymbol{\Psi}}{\partial{\boldsymbol{Y}}})_{s}=&\left(\begin{array}[]{cccc}{J_{\boldsymbol{f}1}}_{s}&&&\\\
&{J_{\boldsymbol{f}2}}_{s}&&\\\ &&\ddots&\\\
&&&{J_{\boldsymbol{f}N}}_{s}\end{array}\right)-\left(\begin{array}[]{cccc}A_{s}&-\frac{A^{T}}{2}&&-\frac{A}{2}\\\
-\frac{A}{2}&A_{s}&\ddots&\\\ &\ddots&\ddots&-\frac{A^{T}}{2}\\\
-\frac{A^{T}}{2}&&-\frac{A}{2}&A_{s}\end{array}\right)\\\
=&\left(\begin{array}[]{cccc}{J_{\boldsymbol{f}1}}_{s}+A_{s}&&&\\\
&{J_{\boldsymbol{f}2}}_{s}+A_{s}&&\\\ &&\ddots&\\\
&&&{J_{\boldsymbol{f}N}}_{s}+A_{s}\end{array}\right)\\\
&-\left(\begin{array}[]{cccc}2A_{s}&-\frac{A^{T}}{2}&&-\frac{A}{2}\\\
-\frac{A}{2}&2A_{s}&\ddots&\\\ &\ddots&\ddots&-\frac{A^{T}}{2}\\\
-\frac{A^{T}}{2}&&-\frac{A}{2}&2A_{s}\end{array}\right),\end{split}$

where the empty elements are 0 blocks. By Theorem 3.5, we need to prove that
$(\frac{\partial\boldsymbol{\Psi}}{\partial{\boldsymbol{Y}}})_{s}$ is
uniformly negative definite under the given conditions. For an arbitrary
$\boldsymbol{Y}$, since ${J_{\boldsymbol{f}}}_{s}+A_{s}\prec 0$ on $D$, we
have ${J_{\boldsymbol{f}i}}_{s}+A_{s}\prec 0$ for all $i$, and thus the first
part in $(\frac{\partial\boldsymbol{\Psi}}{\partial{\boldsymbol{Y}}})_{s}$ is
uniformly negative definite. According to the addition rule for definite
matrices, all that remains to be proved is that the second part is a positive
semi-definite matrix given $A_{s}\succ 0$ and $4A_{s}-AA_{s}^{-1}A^{T}\succeq
0$.

Define block matrices:

$\begin{split}\tilde{A}_{1}=&\left(\begin{array}[]{cc}A_{s}&-\frac{A^{T}}{2}\\\
-\frac{A}{2}&A_{s}\\\ \end{array}\right),\\\
\tilde{A}_{2}=&\left(\begin{array}[]{cc}A_{s}&-\frac{A}{2}\\\
-\frac{A^{T}}{2}&A_{s}\\\ \end{array}\right).\end{split}$

We recall the following lemma

######  Lemma 3.7.

Let $M$ be a $2\times 2$ symmetric block matrix of the form

$M=\left(\begin{array}[]{cc}A&B\\\ B^{T}&D\\\ \end{array}\right).$

Suppose $D$ ($A$) is invertible. If $D\succ 0$ ($A\succ 0$), then $M\succeq 0$
iff the Schur complement of the block $D$ ($A$) of the matrix $M$, $M/D$
($M/A$), is positive semi-definite, i.e., $M/D=A-BD^{-1}B^{T}\succeq 0$
($M/A=D-B^{T}A^{-1}B\succeq 0$).

###### Proof.

See [18]. ∎

######  Lemma 3.8.

Let $M=\left(\begin{array}[]{cc}A&B\\\ B^{T}&D\\\ \end{array}\right)$ be a
positive semi-definite $2\times 2$ block matrix. Define $M^{m,n}$ to be a
block matrix of the form

$M^{m,n}=\left(\begin{array}[]{ccccc}&\vdots&&\vdots&\\\
\cdots&(A)_{mm}&\cdots&(B)_{mn}&\cdots\\\ &\vdots&&\vdots&\\\
\cdots&(B^{T})_{nm}&\cdots&(D)_{nn}&\cdots\\\ &\vdots&&\vdots&\\\
\end{array}\right).$

where $(H)_{ij}$ denotes that the block $H$ is in the $i^{th}$ row and
$j^{th}$ column of the matrix $M^{m,n}$. Here, all the blocks in $M^{m,n}$ are
zero matrices, except for the four blocks written out. Then for any unequal
$m$ and $n$, $M\succeq 0$ iff $M^{m,n}\succeq 0$.

This lemma can be easily proved by noticing that
$\boldsymbol{x}^{T}M\boldsymbol{x}=\boldsymbol{y}^{T}M^{m,n}\boldsymbol{y}$,
where $\boldsymbol{x}$ is an arbitrary vector which has dimensions that are
compatible with matrix $M$, and $\boldsymbol{x}$ is exported to vector
$\boldsymbol{y}$ by arbitrarily choosing the missing elements to have
appropriate dimensions according to $M^{m,n}$.

Observe that the second part in
$(\frac{\partial\boldsymbol{\Psi}}{\partial{\boldsymbol{Y}}})_{s}$ can be
written as $\sum_{i=1}^{N-1}{\tilde{A}_{1}}^{i,i+1}+{\tilde{A}_{2}}^{1,N}$.
Since we have assumed that $A_{s}\succ 0$ and $4A_{s}-AA_{s}^{-1}A^{T}\succeq
0$, Lemma 3.7 implies that both $\tilde{A}_{1}$ and $\tilde{A}_{2}$ are
positive semi-definite. Then, using Lemma 3.8, we have
${\tilde{A}_{1}}^{i,i+1}\succeq 0$ for $i$ from 1 to $N-1$ and
${\tilde{A}_{2}}^{1,N}\succeq 0$, and thus the second part in
$(\frac{\partial\boldsymbol{\Psi}}{\partial{\boldsymbol{Y}}})_{s}$ is positive
semi-definite. Therefore,
$(\frac{\partial\boldsymbol{\Psi}}{\partial{\boldsymbol{Y}}})_{s}$ is negative
definite, and it follows that (13a) is completely synhcronized. ∎

######  Corollary 3.9.

The unidirectional ring network (13a) with identical symmetric coupling matrix
$A$ is completely synhcronized if ${J_{\boldsymbol{f}}}_{s}+A\prec 0$ on $D$,
$A\succ 0$.

###### Proof.

Since $A$ is symmetric, $A=A^{T}=A_{s}$. Since $A\succ 0$, $A$ is invertible.
Then $4A_{s}-AA_{s}^{-1}A^{T}=3A\succeq 0$. By Theorem 3.6, the network is
completely synhcronized. ∎

######  Theorem 3.10.

The bidirectional ring network (13b) is completely synhcronized if
${J_{\boldsymbol{f}}}_{s}\prec 0$ and for each $i$ from 1 to $N$ at least one
of the following two conditions is satisfied

1\. ${A_{i+1,i}}_{s}\succ 0$ and
${A_{i,i+1}}_{s}-\frac{1}{4}(A_{i,i+1}+A_{i+1,i}^{T}){A^{-1}_{i+1,i}}_{s}(A_{i+1,i}+A_{i,i+1}^{T})\succeq
0$;

2\. ${A_{i,i+1}}_{s}\succ 0$ and
${A_{i+1,i}}_{s}-\frac{1}{4}(A_{i+1,i}+A_{i,i+1}^{T}){A^{-1}_{i,i+1}}_{s}(A_{i,i+1}+A_{i+1,i}^{T})\succeq
0$.

###### Proof.

The virtual network of (13b) can be constructed as

$\dot{\boldsymbol{y}}_{i}=\boldsymbol{f}\left(\boldsymbol{y}_{i}\right)+A_{i+1,i}(\boldsymbol{y}_{i+1}-\boldsymbol{y}_{i})+A_{i-1,i}(\boldsymbol{y}_{i-1}-\boldsymbol{y}_{i}),\text{
}i=1,\cdots,N.$ (15)

Define block matrices:

$\tilde{A}_{j,i}=\left(\begin{array}[]{cc}{A_{j,i}}_{s}&-\frac{(A_{j,i}+A_{i,j}^{T})}{2}\\\
-\frac{(A_{i,j}+A_{j,i}^{T})}{2}&{A_{i,j}}_{s}\\\ \end{array}\right),$

Writing the virtual network (15) in the form of virtual system (12), the
symmetric part of $\frac{\partial\boldsymbol{\Psi}}{\partial{\boldsymbol{Y}}}$
is
$D_{{J_{\boldsymbol{f}}}_{s}}-\sum_{i=1}^{N-1}{\tilde{A}_{i+1,i}}^{i,i+1}-{\tilde{A}_{N,1}}^{1,N}$,
where $D_{{J_{\boldsymbol{f}}}_{s}}$ is a block diagonal matrix whoes the
main-diagonal blocks are ${J_{\boldsymbol{f}i}}_{s}$. The fact that
${J_{\boldsymbol{f}}}_{s}$ is negative definite guarantees that
$D_{{J_{\boldsymbol{f}}}_{s}}\prec 0$. By Lemma 3.7, for each $i$ from $1$ to
$N$, when ${A_{i+1,i}}_{s}$ is invertable, $\tilde{A}_{i+1,i}\succeq 0$ and
$\tilde{A}_{i,i+1}\succeq 0$ iff ${A_{i+1,i}}_{s}\succ 0$ and
${A_{i,i+1}}_{s}-\frac{1}{4}(A_{i,i+1}+A_{i+1,i}^{T}){A^{-1}_{i+1,i}}_{s}(A_{i+1,i}+A_{i,i+1}^{T})\succeq
0$; while when ${A_{i,i+1}}_{s}$ is invertable, $\tilde{A}_{i+1,i}\succeq 0$
and $\tilde{A}_{i,i+1}\succeq 0$ iff ${A_{i,i+1}}_{s}\succ 0$ and
${A_{i+1,i}}_{s}-\frac{1}{4}(A_{i+1,i}+A_{i,i+1}^{T}){A^{-1}_{i,i+1}}_{s}(A_{i,i+1}+A_{i+1,i}^{T})\succeq
0$. Thus, we have $\tilde{A}_{i+1,i}\succeq 0$ for $i$ from $1$ to $N-1$ and
$\tilde{A}_{N,1}\succeq 0$. As a consequence of Lemma 3.8,
${\tilde{A}_{i+1,i}}^{i,i+1}\succeq 0$ for $i$ from $1$ to $N-1$ and
${\tilde{A}_{N,1}}^{1,N}\succeq 0$. Then
$(\frac{\partial\boldsymbol{\Psi}}{\partial{\boldsymbol{Y}}})_{s}\prec 0$.
Hence (13b) is completely synhcronized. ∎

######  Corollary 3.11.

The bidirectional ring network (13b) with interactional couplings,
$A_{i,j}=A_{j,i}$, is completely synhcronized if
${J_{\boldsymbol{f}}}_{s}\prec 0$ and ${A_{i,i+1}}_{s}\succ 0$ for all $i$
from 1 to $N$.

###### Proof.

Since ${A_{i,i+1}}_{s}\succ 0$, ${A_{i,i+1}}_{s}$ is invertible. Notice that
when $A_{i,j}=A_{j,i}$, two conditions in Theorem 3.10 become the same and
that we have

$\begin{split}&{A_{i+1,i}}_{s}-\frac{1}{4}(A_{i+1,i}+A_{i,i+1}^{T}){A^{-1}_{i,i+1}}_{s}(A_{i,i+1}+A_{i+1,i}^{T})\\\
=&{A_{i,i+1}}_{s}-{A_{i,i+1}}_{s}{A^{-1}_{i,i+1}}_{s}{A_{i,i+1}}_{s}\\\
=&0\succeq 0.\end{split}$

Therefore, consition 2 in Theorem 3.10 is satisfied, and thus the network is
completely synchronized. ∎

##  4 Networks of the (Duffing/Rayleigh) van der Pol Systems and Numerical
Simulations

We first consider a network of $N$ all-to-all coupled identical van der Pol
oscillators with symmetric linear couplings. In the absence of coupling a
single van der Pol oscillator has following equation of motion

$\ddot{x}+\left(\alpha x^{2}-\gamma\right)\dot{x}+\omega^{2}x=0.$ (16)

It is usually written in two different two-dimensional forms, one by setting
$y=\dot{x}$ and the other by setting
$y=\frac{1}{\omega}(\dot{x}+\frac{\alpha}{3}x^{3}-\gamma x)$. I use the latter
because in this case the Jacobian $J_{\boldsymbol{f}}$ has a globally negative
semi-definite symmetric part if $\alpha\geqslant 0$ and $\gamma<0$. To see
this, we first note that Eq. (16) can be written as

$\begin{split}\ddot{x}+\alpha x^{2}\dot{x}-\gamma\dot{x}&=-\omega^{2}x\\\
\frac{d}{dt}(\dot{x}+\frac{\alpha}{3}x^{3}-\gamma x)&=\omega(-\omega
x).\end{split}$

Hence the two-dimensional form of Eq. (16) can be written as

$\left(\begin{array}[]{l}\dot{x}\\\
\dot{y}\end{array}\right)=\boldsymbol{f}(x,y)=\left(\begin{array}[]{c}\omega
y-\frac{\alpha}{3}x^{3}+\gamma x\\\ -\omega x\end{array}\right).$ (17)

Then the symmetric part of the Jacobian of $\boldsymbol{f}$ is

${J_{\boldsymbol{f}}}_{s}=\left(\begin{array}[]{cc}\gamma-\alpha x^{2}&0\\\
0&0\\\ \end{array}\right).$

It is easy to see that ${J_{f}}_{s}\preceq 0$ for all $x$ and $y$ if
$\alpha\geqslant 0$ and $\gamma\leqslant 0$. Consider a network of $N$ all-to-
all coupled identical van der Pol oscillators with symmetric linear couplings
described by Eq. (14). Write

$\sum_{j=1}^{N}{A_{j}}_{s}=\left(\begin{array}[]{cc}A^{11}_{ts}&A^{12}_{ts}\\\
A^{21}_{ts}&A^{22}_{ts}\\\ \end{array}\right)$ (18)

with $A^{12}_{ts}=A^{21}_{ts}$. We can then establish the following theorem.

######  Theorem 4.1.

For a network of $N$ all-to-all coupled identical van der Pol oscillators with
symmetric linear couplings described by Eq. (14) in which $\boldsymbol{f}$ is
given by Eq. (17) , if $\alpha\geqslant 0$,
$\gamma<A^{11}_{ts}-\frac{(A^{12}_{ts})^{2}}{A^{22}_{ts}}$, and
$A^{22}_{ts}>0$, then the network is completely synchronized.

###### Proof.

According to Proposition 3.3, the network is completely synchronized if

${J_{\boldsymbol{f}}}_{s}-\sum_{j=1}^{N}{A_{j}}_{s}=\left(\begin{array}[]{cc}\gamma-\alpha
x^{2}-A^{11}_{ts}&-A^{12}_{ts}\\\ -A^{21}_{ts}&-A^{22}_{ts}\\\
\end{array}\right)$

is uniformly negative definite on $D$, or, equivalently,
$-({J_{\boldsymbol{f}}}_{s}-\sum_{j=1}^{N}{A_{j}}_{s})\succ 0$ on $D$. Since
$\gamma<A^{11}_{ts}-\frac{(A^{12}_{ts})^{2}}{A^{22}_{ts}}$, $\alpha\geqslant
0$ and $A^{22}_{ts}>0$, we have

$-(\gamma-\alpha
x^{2}-A^{11}_{ts})\geqslant-(\gamma-A^{11}_{ts})\geqslant-(\gamma-A^{11}_{ts}+\frac{(A^{12}_{ts})^{2}}{A^{22}_{ts}})>0,$

and

$\displaystyle-(\gamma-A^{11}_{ts})>\frac{(A^{12}_{ts})^{2}}{A^{22}_{ts}}$
$\displaystyle-(\gamma-A^{11}_{ts})A^{22}_{ts}>(A^{12}_{ts})^{2}$
$\displaystyle-(\gamma-A^{11}_{ts})A^{22}_{ts}-(A^{12}_{ts})^{2}>0.$

Consequently,

$-(-\gamma+\alpha
x^{2}+A^{11}_{ts})^{-1}\geqslant-(-\gamma+A^{11}_{ts})^{-1},$

and

$\begin{split}0<&A^{22}_{ts}-(-\gamma+A^{11}_{ts})^{-1}(A^{12}_{ts})^{2}\\\
\leqslant&A^{22}_{ts}-(-\gamma+\alpha
x^{2}+A^{11}_{ts})^{-1}(A^{12}_{ts})^{2}\\\
=&A^{22}_{ts}-A^{21}_{ts}(-\gamma+\alpha
x^{2}+A^{11}_{ts})^{-1}A^{12}_{ts}.\end{split}$

Lemma 3.7 then implies
$-({J_{\boldsymbol{f}}}_{s}-\sum_{j=1}^{N}{A_{j}}_{s})\succ 0$. ∎

Next we consider $N$ identical van der Pol oscillators coupled by a star
network with symmetric linear couplings described by Eq. (8) where the sum of
coupling matrices is given by Eq. (18) and ${A_{1}}_{s}$ is given by

${A_{1}}_{s}=\left(\begin{array}[]{cc}A^{11}_{1s}&A^{12}_{1s}\\\
A^{21}_{1s}&A^{22}_{1s}\\\ \end{array}\right)$

with $A^{12}_{1s}=A^{21}_{1s}$. Then using Proposition 3.4 and following
similar arguments as in the proof of Theorem 4.1, we have

######  Theorem 4.2.

For $N$ identical van der Pol oscillators coupled by a star network with
symmetric linear couplings described by Eq. (8) in which $\boldsymbol{f}$ is
given by Eq. (17) , if $\alpha\geqslant 0$,
$\gamma<\text{min}\\{A^{11}_{1s}-\frac{(A^{12}_{1s})^{2}}{A^{22}_{1s}},A^{11}_{ts}-\frac{(A^{12}_{ts})^{2}}{A^{22}_{ts}}\\}$,
and $\text{min}\\{A^{22}_{1s},A^{22}_{ts}\\}>0$, then the network is
completely synchronized.

######  Example.

Consider a star network of six coupled identical van der Pol oscillators
described by Eqs. (8) and (17) with $\alpha=1$, $\omega=1$, and $\gamma=5$,
where the coupling matrices are given by

$\displaystyle A_{1}=\left(\begin{array}[]{cc}8&1\\\ 3&4\\\
\end{array}\right),\;A_{2}=\left(\begin{array}[]{cc}3&3\\\ 4&-5\\\
\end{array}\right),\;A_{3}=\left(\begin{array}[]{cc}7&-2\\\ -5&-2\\\
\end{array}\right),$ $\displaystyle A_{4}=\left(\begin{array}[]{cc}1&2\\\
4&2\\\ \end{array}\right),\;A_{5}=\left(\begin{array}[]{cc}8&-3\\\ 1&10\\\
\end{array}\right),\;A_{6}=\left(\begin{array}[]{cc}-3&-4\\\ 2&-6\\\
\end{array}\right).$

In this case,
$\text{min}\\{A^{22}_{1s},A^{22}_{ts}\\}=\text{min}\\{4,3\\}=3>0$, and

$5=\gamma<\text{min}\\{A^{11}_{1s}-\frac{(A^{12}_{1s})^{2}}{A^{22}_{1s}},A^{11}_{ts}-\frac{(A^{12}_{ts})^{2}}{A^{22}_{ts}}\\}=\text{min}\\{7,21\\}=7.$

By Theorem 4.2, the network is completely synchronized.

Figure 1 shows the complete synchronization phenomenon of this network. It can
be seen that all the six oscillators converge into synchronous $x$ and $y$
states very quickly as expected.

(a) The time series of $x_{i}$ for $i$ from 1 to 6

(b) The time series of $y_{i}$ for $i$ from 1 to 6

Figure 1: The time series of self-sustained six identical van der Pol
oscillators coupled by a star network with symmetric linear couplings

We now introduce a Duffing term, $x^{3}$, into the van der Pol equation (16)
so that the euqation becomes

$\ddot{x}+\left(\alpha x^{2}-\gamma\right)\dot{x}+(\omega^{2}+\varepsilon
x^{2})x=0.$ (19)

This is the equation of motion of a Duffing van der Pol oscillator. Suppose
$\alpha\geqslant 0$, $-(\omega^{2}+1)<\gamma<-1$, $\omega\neq 0$, and
$\varepsilon=\frac{\alpha}{3}$. Let $\gamma^{\prime}=\gamma+1$. Since
$-(\omega^{2}+1)<\gamma<-1$, $-\omega^{2}<\gamma^{\prime}<0$. We can therefore
let ${\omega^{\prime}}^{2}=\omega^{2}+\gamma^{\prime}>0$. Equation (19) can
thus be written as

$\displaystyle\ddot{x}+\left(\alpha
x^{2}-\gamma^{\prime}+1\right)\dot{x}+({\omega^{\prime}}^{2}-\gamma^{\prime}+\frac{\alpha}{3}x^{2})x=0$
$\displaystyle\ddot{x}+\left(\alpha
x^{2}-\gamma^{\prime}\right)\dot{x}=-{\omega^{\prime}}^{2}x-\dot{x}-\frac{\alpha}{3}x^{3}+\gamma^{\prime}x$
$\displaystyle\frac{d}{dt}(\dot{x}+\frac{\alpha}{3}x^{3}-\gamma^{\prime}x)=\omega^{\prime}[-\omega^{\prime}x-\frac{1}{\omega^{\prime}}(\dot{x}+\frac{\alpha}{3}x^{3}-\gamma^{\prime}x)]$

Making a transformation
$y=\frac{1}{\omega^{\prime}}(\dot{x}+\frac{\alpha}{3}x^{3}-\gamma^{\prime}x)$,
we get a two dimensional system

$\left(\begin{array}[]{l}\dot{x}\\\
\dot{y}\end{array}\right)=\boldsymbol{f}(x,y)=\left(\begin{array}[]{c}\omega^{\prime}y-\frac{\alpha}{3}x^{3}+\gamma^{\prime}x\\\
-\omega^{\prime}x-y\end{array}\right).$ (20)

The symmetric part of the Jacobian of $\boldsymbol{f}$ is

${J_{\boldsymbol{f}}}_{s}=\left(\begin{array}[]{cc}\gamma^{\prime}-\alpha
x^{2}&0\\\ 0&-1\\\ \end{array}\right).$

${J_{\boldsymbol{f}}}_{s}\prec 0$ as $\alpha\geqslant 0$ and
$\gamma^{\prime}<0$. Using Corollary 3.11, we conclude the following theorem

######  Theorem 4.3.

For a bidirectional ring network of $N$ identical coupled Duffing van der Pol
oscillators with balanced linear couplings described by Eq. (13b) in which
$\boldsymbol{f}$ is given by Eq. (20) , if $\alpha\geqslant 0$,
$\gamma^{\prime}<0$, and ${A_{i,i+1}}_{s}\succ 0$ for all $i$, then the
network is completely synchronized.

For a unidirectional ring network of identical coupled Duffing van der Pol
oscillators, we have the following example

######  Example.

Consider a unidirectional ring network of six coupled identical Duffing van
der Pol oscillators described by Eqs. (13a) and (20) with $\alpha=1$,
$\omega=1$, and $\gamma=-1$, where the identical coupling matrix is given by

$A=\left(\begin{array}[]{cc}3/4&1/4\\\ 1/8&1/2\\\ \end{array}\right),$

so that $A_{s}\succ 0$, $4A_{s}-AA_{s}^{-1}A^{T}\succ 0$, and

${J_{\boldsymbol{f}}}_{s}+A_{s}=\left(\begin{array}[]{cc}-1/4-x^{2}&3/16\\\
3/16&-1/2\\\ \end{array}\right).$

Notice that

$\displaystyle-(-\frac{1}{2})>0$
$\displaystyle-(-\frac{1}{4}-x^{2})-(\frac{1}{2})^{-1}(-\frac{3}{16})^{2}=\frac{23}{128}+x^{2}>0.$

By Lemma 3.7, $-({J_{\boldsymbol{f}}}_{s}+A_{s})\succ 0$, or
${J_{\boldsymbol{f}}}_{s}+A_{s}\prec 0$. Theorem 3.6 then implies the network
is completely synchronized.

Figure 2 shows the time series of $x$ and $y$ states of this network. It can
be seen that the network achieves complete synchronization before it decays to
the equilibrium position.

(a) The time series of $x_{i}$ for $i$ from 1 to 6

(b) The time series of $y_{i}$ for $i$ from 1 to 6

Figure 2: The time series of overdamped six identical Duffing van der Pol
oscillators coupled by a unidirectional ring network with an identical linear
coupling

Next, we introduce a Rayleigh term, $\dot{x}^{3}$, into the van der Pol
equation (16) to get a Rayleigh van der Pol oscillator

$\ddot{x}+\left(\alpha
x^{2}+\beta\dot{x}^{2}-\gamma\right)\dot{x}+\omega^{2}x=0.$ (21)

Letting $y=\dot{x}$, its two dimensional form can be written as

$\left(\begin{array}[]{l}\dot{x}\\\
\dot{y}\end{array}\right)=\boldsymbol{f}(x,y)=\left(\begin{array}[]{c}y\\\
-\left(\alpha x^{2}+\beta y^{2}-\gamma\right)y-\omega^{2}x\end{array}\right).$
(22)

We assume all-to-all symmetric coupling as described in Eq. (14). We also
assume full-state linear coupling, i.e., the coupling matrices $A_{i}$ are
diagonal for all $i$. If $x_{i}$ and $y_{i}$ represent the “position” and
“velocity” of $i^{\text{th}}$ node, respectively, then full-state coupling
means that each node adjusts its own position and velocity according to the
difference between its position and velocity with other nodes, respectively
[19]. We can write

${J_{\boldsymbol{f}}}_{s}-\sum_{j=1}^{N}{A_{j}}_{s}=\left(\begin{array}[]{cc}-\frac{c}{N}\sum_{j=1}^{N}a_{j}&\frac{1-2\alpha
xy-\omega^{2}}{2}\\\ \frac{1-2\alpha xy-\omega^{2}}{2}&-(\alpha x^{2}+3\beta
y^{2}-\gamma)-\frac{c}{N}\sum_{j=1}^{N}b_{j}\\\ \end{array}\right).$ (23)

where $c$ is a constant coupling strength, and $a_{j}$ and $b_{j}$ are
coefficients that represent the coupling effect from $j^{\text{th}}$ node to
others. In this case, ${J_{\boldsymbol{f}}}_{s}$ is neither negative definite
nor negative semi-definite. Furthermore, the definiteness of
${J_{\boldsymbol{f}}}_{s}-\sum_{j=1}^{N}{A_{j}}_{s}$ cannot be determined
globally since it depends on the position $(x,y)$ in the virtual phase plane
$D$. We can nevertheless find the contraction region of the virtual system of
this network.

######  Lemma 4.4.

For a network of $N$ all-to-all coupled identical Rayleigh van der Pol
oscillators with symmetric linear couplings described by Eq. (14) in which
$\boldsymbol{f}$ is given by Eq. (22), if $\frac{c}{N}\sum_{j=1}^{N}a_{j}>0$,
then the network is completely synchronized if all trajectories,
$(x_{k}(t)\text{ }y_{k}(t))^{T}$, start within the region
$C_{RvdP}:=\left\\{(x\text{
}y)^{T}\in\mathbb{R}^{2}:\frac{c^{2}}{N^{2}}\sum_{i,j=1}^{N}a_{i}b_{j}+\mu(x,y)\frac{c}{N}\sum_{j=1}^{N}a_{j}-\nu^{2}(x,y)>0\right\\}$
and ramain there for all $k$ from 1 to $N$, where $\mu(x,y)=\alpha
x^{2}+3\beta y^{2}-\gamma$ and $\nu(x,y)=\frac{1-2\alpha xy-\omega^{2}}{2}$.

###### Proof.

Since the coupling is all-to-all and symmetric, Proposition 3.3 implies that
the network is completely synchronized if $(x_{i}(t)\text{ }y_{i}(t))^{T}$ are
always in the region where
${J_{\boldsymbol{f}}}_{s}-\sum_{j=1}^{N}{A_{j}}_{s}\prec 0$ for all $i$. Using
Eq. (23) and Lemma 3.7,
${J_{\boldsymbol{f}}}_{s}-\sum_{j=1}^{N}{A_{j}}_{s}\prec 0$ iff
$-{J_{\boldsymbol{f}}}_{s}+\sum_{j=1}^{N}{A_{j}}_{s}\succ 0$ iff
$\frac{c}{N}\sum_{j=1}^{N}a_{j}>0$ and $(\alpha x^{2}+3\beta
y^{2}-\gamma)+\frac{c}{N}\sum_{j=1}^{N}b_{j}-\frac{(1-2\alpha
xy-\omega^{2})^{2}}{4}(\frac{c}{N}\sum_{j=1}^{N}a_{j})^{-1}>0$. ∎

To determine whether all trajectories will remain within $C_{RvdP}$, we need
to find a trapping region of the virtual system of this network that is
contained in $C_{RvdP}$. We make a further assumption that the coupling
matrices are identical and symmetric, i.e., $a_{i}=b_{i}=a$ for all $i$. The
dynamics of this network can be expressed as

$\begin{split}\dot{x}_{i}&=y_{i}+\frac{c}{N}a\sum_{j=1}^{N}(x_{j}-x_{i})\\\
\dot{y}_{i}&=-\left(\alpha{x_{i}}^{2}+\beta{y_{i}}^{2}-\gamma\right)y_{i}-\omega^{2}x_{i}+\frac{c}{N}a\sum_{j=1}^{N}(y_{j}-y_{i}).\end{split}$
(24)

Using the LaSalle’s invariance principle, we can obtain the following lemma.
The proof is given in Appendix B.

######  Lemma 4.5.

For a network of $N$ all-to-all coupled identical Rayleigh van der Pol
oscillators with identical full-state coupling described by Eq. (24), the
origin is globally asymptotically stable if $\alpha\geqslant 0$,
$\beta\geqslant 0$, $ca\geqslant 0$, $\omega\neq 0$, and $\gamma<0$.

As a result, the solution of Eq. (24) with $\alpha\geqslant 0$,
$\beta\geqslant 0$, $ca\geqslant 0$, $\omega\neq 0$, and $\gamma<0$ will
enventually converge to the origion regardless of the initial conditions.
However, as stated in Definition 3.1, usually 0 is not treated as a
synchronous state, so we want to show that all overdamped coupled Rayleigh van
der Pol oscillators in the network (24) can reach complete synchronization
before they stop oscillating.

Using Eq. (7), the virtual system of (24) can be constrcuted as follows:

$\left(\begin{array}[]{l}\dot{x}\\\
\dot{y}\end{array}\right)=\left(\begin{array}[]{c}y+\frac{ca}{N}(\sum_{j=1}^{N}x_{j}-Nx)\\\
-\left(\alpha x^{2}+\beta
y^{2}-\gamma\right)y-\omega^{2}x+\frac{ca}{N}(\sum_{j=1}^{N}y_{j}-Ny)\end{array}\right).$
(25)

######  Lemma 4.6.

The disk of radius $r$

${R_{RvdP}}_{\omega^{2}=1}:=\left\\{(x,y)\in\mathbb{R}^{2}:H(x,y)=x^{2}+y^{2}\leqslant
r^{2}\right\\}$

is a trapping region of the virtual system (25) with $\alpha\geqslant 0$,
$\beta\geqslant 0$, $ac\geqslant 0$, $\omega^{2}=1$, and $\gamma<0$.

###### Proof.

Recall that the gradient of $H$ at a point on a level curve of $H$ is
perpendicular to the level curve and pointing in the direction of fastest
increase. To verify that the vector field of (25) points inward on the
boundary of the disk ${R_{RvdP}}_{\omega^{2}=1}$ at all times, it suffices to
show that $\nabla
H\cdot\langle\dot{x},\dot{y}\rangle=2x\dot{x}+2y\dot{x}\leqslant 0$ on the
circle $x^{2}+y^{2}=r^{2}$ at all times. At any point $(x,y)$ on
$x^{2}+y^{2}=r^{2}$, we have

$\begin{split}\nabla
H\cdot\langle\dot{x}(t),\dot{y}(t)\rangle=&2x\dot{x}(t)+2y\dot{y}(t)\\\
=&2xy+2\frac{ca}{N}(\sum_{j=1}^{N}xx_{j}(t)-Nx^{2})\\\ &-2\left(\alpha
x^{2}+\beta
y^{2}-\gamma\right)y^{2}-2xy+2\frac{ca}{N}(\sum_{j=1}^{N}yy_{j}(t)-Ny^{2})\\\
\leqslant&2\frac{ca}{N}(\sum_{j=1}^{N}(xx_{j}(t)+yy_{j}(t))-N(x^{2}+y^{2}))\\\
=&2\frac{ca}{N}(\sum_{j=1}^{N}rr_{j}(t)\cos\theta_{i}(t)-Nr^{2})\\\
\leqslant&2\frac{ca}{N}(\sum_{j=1}^{N}rr_{j}(t)-Nr^{2}),\end{split}$

where $r_{j}(t)$ represent the distances between the origin and the point
$(x_{j}(t),y_{j}(t))$ at time $t$ in the phase plane of (25) and
$\theta_{i}(t)$ represent the angles between the position vectors $\langle
x,y\rangle$ and $\langle x_{j}(t),y_{j}(t)\rangle$. By Lemma 4.5, the origin
is a globally asymptotically stable point of the original network, which means

$\sum_{j=1}^{N}r_{j}^{2}(t)\leqslant\sum_{j=1}^{N}r_{j}^{2}(0),$

for all $t>0$. Suppose that the starting points $({x_{j}}_{0},{y_{j}}_{0})$ of
all oscillators lie in ${R_{RvdP}}_{\omega^{2}=1}$, $r_{j}(0)\leqslant r$ for
all $j$ from 1 to $N$. Then

$\sum_{j=1}^{N}r_{j}^{2}(t)\leqslant\sum_{j=1}^{N}{x^{2}_{j}}_{0}+{y^{2}_{j}}_{0}\leqslant
Nr^{2}.$

Using the Cauchy–Schwarz inequality, we have

$\frac{(\sum_{j=1}^{N}r_{j}(t))^{2}}{N}\leqslant\sum_{j=1}^{N}r_{j}^{2}(t)\leqslant
Nr^{2}.$

It follows that

$\sum_{j=1}^{N}r_{j}(t)\leqslant Nr,$

or,

$\sum_{j=1}^{N}rr_{j}(t)\leqslant Nr^{2},$

Therefore,

$\nabla H\cdot\langle\dot{x}(t),\dot{y}(t)\rangle\leqslant
2\frac{ca}{N}(\sum_{j=1}^{N}rr_{j}(t)-Nr^{2})\leqslant 0$

∎

By Lemma 4.4, when $\omega^{2}=1$, the contraction region $C_{RvdP}$ of the
virtual system (25) becomes

$\begin{split}{C_{RvdP}}_{\omega^{2}=1}=&\left\\{(x,y)\in\mathbb{R}^{2}:c^{2}a^{2}+ca\mu(x,y)-\nu^{2}(x,y)>0\right\\}\\\
=&\left\\{(x,y)\in\mathbb{R}^{2}:g(x,y)>0\right\\}.\end{split}$

where $g(x,y)=c^{2}a^{2}-ca\gamma+ca\alpha x^{2}+3ca\beta
y^{2}-\alpha^{2}x^{2}y^{2}$.

######  Lemma 4.7.

Provided $ca>\max\\{\gamma,0\\}$, ${R_{RvdP}}_{\omega^{2}=1}$ is contained in
${C_{RvdP}}_{\omega^{2}=1}$ in the following four cases

1\. $3\beta r^{2}>\gamma-ca$ and $\alpha=0$;

2\. $r^{2}\leqslant\frac{ca\alpha-3ca\beta}{\alpha^{2}}$ and $3\beta
r^{2}>\gamma-ca$ with $0\neq\alpha>3\beta$;

3\. $r^{2}\leqslant\frac{3ca\beta-ca\alpha}{\alpha^{2}}$ and $\alpha
r^{2}>\gamma-ca$ with $0\neq\alpha<3\beta$;

4\. $r^{2}\geqslant\max\\{\frac{ca\alpha-3ca\beta}{\alpha^{2}},\frac{3ca\beta-
ca\alpha}{\alpha^{2}}\\}$ and
$\frac{ca\alpha+3ca\beta-2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}<r^{2}<\frac{ca\alpha+3ca\beta+2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}$
with $\alpha$ and $\beta$ satisfying $3ca\alpha\beta+\alpha^{2}(ca-\gamma)>0$.

###### Proof.

Since $ca>0$ and $ca>\gamma$, $c^{2}a^{2}-ca\gamma=ca(ca-\gamma)>0$. Take an
arbitrary point $(x^{*},y^{*})$ in ${R_{RvdP}}_{\omega^{2}=1}$. We want to
show $g(x^{*},y^{*})>0$ so that $(x^{*},y^{*})\in{C_{RvdP}}_{\omega^{2}=1}$.

Case 1: $3\beta r^{2}>\gamma-ca$ and $\alpha=0$. In this case,

$g(x^{*},y^{*})=G(y^{*})=ca(ca-\gamma)+3ca\beta{y^{*}}^{2}=ca(ca-\gamma+3\beta{y^{*}}^{2}),$

where $y^{*}\in[-r,r]$. For $\beta\geqslant 0$, $G(y)$ has a minimum at $y=0$,
and we have $G(0)=ca(ca-\gamma)>0$. For $\beta<0$ and $y\in[-r,r]$, $G(y)$ has
minima at $y=\pm r$. Since $3\beta r^{2}>\gamma-ca$, $G(\pm
r)=ca(ca-\gamma+3\beta r^{2})>0$. Consequently, $G(y)>0$ for all $y\in[-r,r]$.
Therefore, $g(x^{*},y^{*})=G(y^{*})>0$.

Case 2: $r^{2}\leqslant\frac{ca\alpha-3ca\beta}{\alpha^{2}}$ and $3\beta
r^{2}>\gamma-ca$ with $0\neq\alpha>3\beta$.

The partial derivatives of $g(x,y)$ are

$\begin{split}g_{x}(x,y)&=2ca\alpha x-2\alpha^{2}y^{2}x\\\
g_{y}(x,y)&=6ca\beta y-2\alpha^{2}x^{2}y.\end{split}$

Notice that $g_{x}(x,y)=0$ when $x=0$ or $y^{2}=\frac{ca}{\alpha}$ (if
$\alpha>0$) and that $g_{y}(x,y)=0$ when $y=0$ or
$x^{2}=\frac{3ca\beta}{\alpha^{2}}$ (if $\alpha\neq 0$ and $\beta\geqslant
0$). Thus, for $\beta<0$, $g(x,y)$ only has one critical point, $(0,0)$, at
which $g(0,0)=ca(ca-\gamma)>0$. For $\beta\geqslant 0$, $g(x,y)$ has five
critical points, $(0,0)$,
$(\pm\sqrt{\frac{3ca\beta}{\alpha^{2}}},\pm\sqrt{\frac{ca}{\alpha}})$, if
$\alpha>0$. However,

$r^{2}\leqslant\frac{ca\alpha-3ca\beta}{\alpha^{2}}\leqslant\frac{ca\alpha+3ca\beta}{\alpha^{2}},$

which indicates that there is only one critical point, $(0,0)$, in the
interior of ${R_{RvdP}}_{\omega^{2}=1}$ at which $g(0,0)>0$.

On the boundary, $x^{2}+y^{2}=r^{2}$, of ${R_{RvdP}}_{\omega^{2}=1}$, we have

$\begin{split}g(x,y)=h(x)\equiv&ca(ca-\gamma)+ca\alpha
x^{2}+3ca\beta(r^{2}-x^{2})-\alpha^{2}x^{2}(r^{2}-x^{2})\\\
=&\alpha^{2}x^{4}+(ca\alpha-3ca\beta-\alpha^{2}r^{2})x^{2}+ca(ca-\gamma)+3ca\beta
r^{2}>0\end{split}$

$h^{\prime}(x)=4\alpha^{2}x^{3}+2(ca\alpha-3ca\beta-\alpha^{2}r^{2})x.$

Since $r^{2}\leqslant\frac{ca\alpha-3ca\beta}{\alpha^{2}}$,
$ca\alpha-3ca\beta-\alpha^{2}r^{2}\geqslant 0$, so $h^{\prime}(x)=0$ only at
$x=0$. Furthermore, $h^{\prime}(x)>0$ when $x>0$ and $h^{\prime}(x)<0$ when
$x<0$. It follows that $h(x)$ achieves a minimum at $x=0$.

$h(0)=ca(ca-\gamma)+3ca\beta r^{2}=ca(ca-\gamma+3\beta).$

Since $3\beta r^{2}>\gamma-ca$, $h(0)>0$. Then $h(x)>0$ for all $x$. Thus,
$g(x^{*},y^{*})>0$.

Case 3: $r^{2}\leqslant\frac{3ca\beta-ca\alpha}{\alpha^{2}}$ and $\alpha
r^{2}>\gamma-ca$ with $0\neq\alpha<3\beta$. For $\alpha<0$, $g(x,y)$ has only
one critical point $(0,0)$ at which $g(0,0)=ca(ca-\gamma)>0$. For $\alpha>0$,

$r^{2}\leqslant\frac{3ca\beta-
ca\alpha}{\alpha^{2}}<\frac{ca\alpha+3ca\beta}{\alpha^{2}}.$

Only one critical point in ${R_{RvdP}}_{\omega^{2}=1}$ is still $(0,0)$, and
$g(0,0)>0$.

On the boundary of ${R_{RvdP}}_{\omega^{2}=1}$, $g(x,y)=h(x)$ for
$x\in[-r,r]$. Since $\alpha<3\beta$, $h^{\prime}(x)=0$ at $x=0$ and at

$x=x_{\pm}=\pm\frac{\sqrt{-2(ca\alpha-3ca\beta-\alpha^{2}r^{2})}}{2\alpha}.$

We have

$h^{\prime\prime}(x)=12\alpha^{2}x^{2}+2(ca\alpha-3ca\beta-\alpha^{2}r^{2}),$

so that

$\displaystyle h^{\prime\prime}(0)=$ $\displaystyle
2(ca\alpha-3ca\beta-\alpha^{2}r^{2})<0,$ $\displaystyle
h^{\prime\prime}(x_{1\pm})=$
$\displaystyle-4(ca\alpha-3ca\beta-\alpha^{2}r^{2})>0.$

It follows that $h(x)$ reaches its global minimum value at $x=x_{+}$ or at
$x=x_{-}$. But since $r^{2}\leqslant\frac{3ca\beta-ca\alpha}{\alpha^{2}}$,

$\displaystyle 2\alpha^{2}r^{2}\leqslant-2(ca\alpha-3ca\beta)$ $\displaystyle
4\alpha^{2}r^{2}\leqslant-2(ca\alpha-3ca\beta-\alpha^{2}r^{2})$ $\displaystyle
r^{2}\leqslant\frac{-2(ca\alpha-3ca\beta-\alpha^{2}r^{2})}{4\alpha^{2}}=x^{2}_{\pm}$

Therefore, for $x\in[-r,r]$, $h(x)$ achieves its minimum at $x=r$ or at
$x=-r$. Since $\alpha r^{2}>\gamma-ca$,

$h(\pm r)=ca\alpha r^{2}+ca(ca-\gamma)=ca(\alpha r^{2}+ca-\gamma)>0.$

So $h(x)>0$ for all $x\in[-r,r]$. Therefore, $g(x^{*},y^{*})>0$.

Case 4:
$r^{2}\geqslant\max\\{\frac{ca\alpha-3ca\beta}{\alpha^{2}},\frac{3ca\beta-
ca\alpha}{\alpha^{2}}\\}$ and
$\frac{ca\alpha+3ca\beta-2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}<r^{2}<\frac{ca\alpha+3ca\beta+2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}$
with $3ca\alpha\beta+\alpha^{2}(ca-\gamma)>0$.

Notice from $3ca\alpha\beta+\alpha^{2}(ca-\gamma)>0$ that $\alpha\neq 0$.
$g(x,y)$ has one critical point, $(0,0)$, if $\alpha<0$ or $\beta<0$,
otherwise $g(x,y)$ has five critical points, $(0,0)$ and
$(\pm\sqrt{\frac{3ca\beta}{\alpha^{2}}},\pm\sqrt{\frac{ca}{\alpha}})$, that
could lie in the interior of ${R_{RvdP}}_{\omega^{2}=1}$ since

$\frac{ca\alpha+3ca\beta}{\alpha^{2}}\leqslant\frac{ca\alpha+3ca\beta+2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}$

for $\alpha>0$ and $\beta>0$. $g(0,0)>0$ as before. At points
$(\pm\sqrt{\frac{3ca\beta}{\alpha^{2}}},\pm\sqrt{\frac{ca}{\alpha}})$,

$\displaystyle g(x,y)=$ $\displaystyle
ca(ca-\gamma)+ca\alpha\frac{3ca\beta}{\alpha^{2}}+3ca\beta\frac{ca}{\alpha}-\alpha^{2}\frac{3ca\beta}{\alpha^{2}}\frac{ca}{\alpha}$
$\displaystyle=$ $\displaystyle
ca(ca-\gamma)+\frac{3c^{2}a^{2}\beta}{\alpha}+\frac{3c^{2}a^{2}\beta}{\alpha}-\frac{3c^{2}a^{2}\beta}{\alpha}$
$\displaystyle=$ $\displaystyle ca(ca-\gamma)+\frac{3c^{2}a^{2}\beta}{\alpha}$
$\displaystyle=$
$\displaystyle\frac{ca}{\alpha^{2}}[\alpha^{2}(ca-\gamma)+3ca\alpha\beta]>0.$

On the boundary of ${R_{RvdP}}_{\omega^{2}=1}$, we have seen in Case 3 that
$h(x)$ reaches its global minimum value at $x=x_{+}$ or at $x=x_{-}$ because
$r^{2}\geqslant\frac{ca\alpha-3ca\beta}{\alpha^{2}}$. In this case,
$x_{\pm}\in[-r,r]$ since $r^{2}\geqslant\frac{3ca\beta-ca\alpha}{\alpha^{2}}$.
Thus, for $x\in[-r,r]$, $h(x)$ has a minimum at $x=x_{+}$ or at $x=x_{-}$.
Note that

$\displaystyle h(x_{\pm})=l(r^{2})\equiv$
$\displaystyle\frac{-(ca\alpha-3ca\beta-\alpha^{2}r^{2})^{2}}{4\alpha^{2}}+ca(ca-\gamma)+3ca\beta
r^{2}$ $\displaystyle=$
$\displaystyle-\frac{\alpha^{2}}{4}r^{4}+\frac{ca\alpha+3ca\beta}{2}r^{2}-\frac{(ca\alpha-3ca\beta)^{2}}{4\alpha^{2}}+ca(ca-\gamma).$

Since $3ca\alpha\beta+\alpha^{2}(ca-\gamma)>0$,
$3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)>0$. Observe that at
$r^{2}=\frac{ca\alpha+3ca\beta\pm
2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}$,
$l(r^{2})=0$. Hence $l(r^{2})>0$ given

$\displaystyle\frac{ca\alpha+3ca\beta-2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}<r^{2}$
$\displaystyle<\frac{ca\alpha+3ca\beta+2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}.$

Thus, $h(x)>0$ for all $x\in[-r,r]$, and therefore $g(x^{*},y^{*})>0$.

Therefore, ${R_{RvdP}}_{\omega^{2}=1}\subset{C_{RvdP}}_{\omega^{2}=1}$ in all
four cases. ∎

Due to Lemma 4.6, any trajectory starting within ${R_{RvdP}}_{\omega^{2}=1}$
of the network (24) with appropriate coefficients will remain there and thence
in ${C_{RvdP}}_{\omega^{2}=1}$ due to Lemma 4.7. Consequently, Lemma 4.4
implies the following theorem.

######  Theorem 4.8.

Consider a network of $N$ all-to-all coupled identical Rayleigh van der Pol
oscillators with an identical full-state linear coupling described by Eq. (24)
with $\alpha\geqslant 0$, $\beta\geqslant 0$, $ca>0$, $\gamma<0$, and
$\omega^{2}=1$. In the cases stated in Lemma 4.7, the network is completely
synchronized if its starting points, $({x_{i}}_{0},{y_{i}}_{0})$, are taken in
the region ${R_{RvdP}}_{\omega^{2}=1}$ for all $i$ from 1 to $N$.

Note that this theorem is based on the contraction theory, so the rate of
synchronization is exponential due to Theorem 2.1. However, as can be seen in
Lemma 4.9 below, the rate of convergence of the network to the origin is also
exponential if the initial conditions are in a neighborhood of the origin.
Therefore, for initial conditions in ${R_{RvdP}}_{\omega^{2}=1}$ and close to
the origin, we usually cannot determine whether the network can reach complete
synchronization before it decays to the equilibrium position without numerical
methods. In general, for large coupling constants $c$ and $a$, the complete
synzhronization can be achieved first. For initial conditions in
${R_{RvdP}}_{\omega^{2}=1}$ and far from the origin, in most cases the network
will get completely synchronied before the oscillations stop.

######  Example.

Consider a network of six all-to-all coupled identical Rayleigh van der Pol
oscillators with an identical full-state linear coupling described by Eq. Eq.
(24) with $\alpha=1$, $\beta=1$, $\gamma=-0.1$, $\omega=1$, $a=1$, $c=1$, and
$N=6$. Suppose the starting points of this system are

$\displaystyle(x_{1}(0),y_{1}(0))=(1,1.5),\;(x_{2}(0),y_{2}(0))=(0.3,-1.2),$
$\displaystyle(x_{3}(0),y_{3}(0))=(-0.5,-1.4),\;(x_{4}(0),y_{4}(0))=(0.5,-0.5),$
$\displaystyle(x_{5}(0),y_{5}(0))=(1.7,0.9),\;(x_{6}(0),y_{6}(0))=(0.3,-1).$

Choose the trapping region, ${R_{RvdP}}_{\omega^{2}=1}$, to be a disk of
radius $r=2$. Then all the starting points are in ${R_{RvdP}}_{\omega^{2}=1}$.
By Lemma 4.7 (case 4 of the lemma), ${R_{RvdP}}_{\omega^{2}=1}$ is included in
the contraction region, ${C_{RvdP}}_{\omega^{2}=1}$, of the virtual system of
this network since

$\displaystyle
4=r^{2}\geqslant\max\\{\frac{ca\alpha-3ca\beta}{\alpha^{2}},\frac{3ca\beta-
ca\alpha}{\alpha^{2}}\\}=2,$
$\displaystyle-0.05\approx\frac{ca\alpha+3ca\beta-2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}<r^{2}$
$\displaystyle<\frac{ca\alpha+3ca\beta+2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}\approx
8.05,$ $\displaystyle 3ca\alpha\beta+\alpha^{2}(ca-\gamma)=4.1>0.$

By Theorem 4.8, the network is completely synchronized.

It can be seen from Fig. 3 that the trajectories, $(x_{i}(t),y_{i}(t))$, of
all oscillators in this network converge to the origin as predicted by Lemma
4.5, but all oscillators quickly reach complete synchronization before the
amplitudes of the oscillations decrease to zero.

(a) The time series of $x_{i}$ for $i$ from 1 to 6

(b) The time series of $y_{i}$ for $i$ from 1 to 6

(c) Phase portraits

Figure 3: The time series and phase portraits of underdamped six all-to-all
coupled identical Rayleigh van der Pol oscillators with an identical full-
state linear coupling

Next we consider the case when $\gamma>0$. In this case, the origin may not be
an asymptotically stable equilibrium point and there may be a periodic
behavior of the system. We want to find the conditions that make the system
(24) self-sustained so that a stable limit cycle exists in its state space.

######  Lemma 4.9.

Suppose $ac>0$ and $\omega\neq 0$. When $|\gamma|\geqslant 2|\omega|$ and
$\gamma<0$, the origin is a locally exponentially stable stationary point for
the system (24) (overdamped oscillators for $|\gamma|>2|\omega|$). When
$|\gamma|<2|\omega|$ and $\gamma<2ca$, the system (24) experiences a Hopf
bifurcation at the origin when $\gamma$ crosses 0, and the origin is locally
exponentially stable for $\gamma<0$ (underdamped oscillators) and unstable for
$\gamma>0$, and the bifurcation is supercritical if $\alpha\geqslant 0$ and
$\beta\geqslant 0$.

###### Proof.

The Jacobian of the linearized system (24) about the origin is given by

$J=\left(\begin{array}[]{cc}-\frac{ca}{N}L_{N}&I_{N}\\\
-\omega^{2}I_{N}&\gamma I_{N}-\frac{ca}{N}L_{N}\\\ \end{array}\right),$

where $I_{N}$ is an $N\times N$ identical matrix and $L_{N}$ is an $N\times N$
Laplacian matrix of complete graph whose elements are given by

${L_{N}}_{ij}=\left\\{\begin{array}[]{lll}N-1&\text{ if { }}i=j\\\ -1&\text{
otherwise }\end{array}\right..$

Hence an eigenvalue $\lambda$ of $J$ satisfies

$\text{det}\left(\begin{array}[]{cc}-\frac{ca}{N}L_{N}-\lambda I_{N}&I_{N}\\\
-\omega^{2}I_{N}&(\gamma-\lambda)I_{N}-\frac{ca}{N}L_{N}\\\
\end{array}\right)=0.$

Since $I_{N}$ is an identical matrix, the lower two blocks in matrix $J$
commute. Using a property of block matrix determinant, we have

$\begin{split}0=&\text{det}\left(\begin{array}[]{cc}-\frac{ca}{N}L_{N}-\lambda
I_{N}&I_{N\times N}\\\
-\omega^{2}I_{N}&(\gamma-\lambda)I_{N}-\frac{ca}{N}L_{N}\\\
\end{array}\right)\\\ =&\text{det}((-\frac{ca}{N}L_{N}-\lambda
I_{N})((\gamma-\lambda)I_{N}-\frac{ca}{N}L_{N})+\omega^{2}I^{2}_{N})\\\
=&\text{det}(\frac{c^{2}a^{2}}{N^{2}}L^{2}_{N}+\frac{ca}{N}(2\lambda-\gamma)L_{N}+[\omega^{2}-\lambda(\gamma-\lambda)]I_{N}).\end{split}$

It can be calculated that

${L^{2}_{N}}_{ij}=\left\\{\begin{array}[]{lll}N^{2}-N&\text{ if { }}i=j\\\
-N&\text{ otherwise }\end{array}\right..$

We therefore have

$\begin{split}0=\text{det}\left(\begin{array}[]{cccc}p(\lambda)&q(\lambda)&\cdots&q(\lambda)\\\
q(\lambda)&p(\lambda)&\cdots&q(\lambda)\\\ \vdots&\vdots&\ddots&\vdots\\\
q(\lambda)&q(\lambda)&\cdots&p(\lambda)\end{array}\right)_{N\times
N},\end{split}$

where

$\displaystyle p(\lambda)=$
$\displaystyle\lambda^{2}+(2ca\frac{N-1}{N}-\gamma)\lambda+c^{2}a^{2}-ca\gamma+\frac{ca\gamma-c^{2}a^{2}}{N}+\omega^{2},$
$\displaystyle q(\lambda)=$
$\displaystyle-\frac{2ca}{N}\lambda+\frac{ca\gamma-c^{2}a^{2}}{N}.$

######  Lemma 4.10.

For $N\times N$ matries $W_{N}$ and $S_{N}$ whose elements are defined as

${W_{N}}_{ij}=\left\\{\begin{array}[]{ll}a&\text{ if { }}i=j\\\
b&\text{otherwise}\end{array}\right.$

and

${S_{N}}_{ij}=\left\\{\begin{array}[]{ll}a&\text{ if { }}i=j\neq 1\\\
b&\text{otherwise}\end{array}\right.$

respectively, we have

$\begin{split}\text{det}[W_{N}]=&(a-b)^{N-1}[a+(N-1)b]\\\
\text{det}[S_{N}]=&b(a-b)^{N-1}.\end{split}$ (26)

###### Proof.

The case when $N=1$ is trivial. We assume Eq. (26) is true for $N\geqslant 1$
and choose any row or column to calculate the determinants of $W_{N+1}$ and
$S_{N+1}$. By switching rows, one can find

$\begin{split}\text{det}[W_{N+1}]=&a(\text{det}[W_{N}])-Nb(\text{det}[S_{N}])\\\
=&a(a-b)^{N-1}[a+(N-1)b]-Nb^{2}(a-b)^{N-1}\\\ =&(a-b)^{N-1}[a^{2}+(N-1)ab-
Nb^{2}]\\\ =&(a-b)^{N-1}[(a-b)(a+Nb)]\\\ =&(a-b)^{N}(a+Nb),\end{split}$

and

$\begin{split}\text{det}[S_{N+1}]=&b(\text{det}[W_{N}])-Nb(\text{det}[S_{N}])\\\
=&b(a-b)^{N-1}[a+(N-1)b]-Nb^{2}(a-b)^{N-1}\\\ =&(a-b)^{N-1}[ab-b^{2}]\\\
=&b(a-b)^{N}.\end{split}$

By mathematical induction, Lemma 4.10 has been proved. ∎

Using Lemma 4.10, $\lambda$ satisfies

$(p-q)^{N-1}[p+(N-1)q](\lambda)=0.$ (27)

One can find

$(p-q)(\lambda)=\lambda^{2}+(2ca-\gamma)\lambda+c^{2}a^{2}-ca\gamma+\omega^{2},\\\
$ (28)

$[p+(N-1)q](\lambda)=\lambda^{2}-\gamma\lambda+\omega^{2}.$ (29)

It follows that two roots of the function (28) are a pair of solutions of Eq.
(27) with multiplicity $N-1$ and two roots of the function (29) are a pair of
solutions of Eq. (27) with multiplicity $1$. One can find the roots of (28)
and (29) are

$\displaystyle\lambda_{1\pm}=$
$\displaystyle\frac{1}{2}(-2ac+\gamma\pm\sqrt{\gamma^{2}-4\omega^{2}}),$
$\displaystyle\lambda_{2\pm}=$
$\displaystyle\frac{1}{2}(\gamma\pm\sqrt{\gamma^{2}-4\omega^{2}})$

respectively. Assume that $ac\geqslant 0$ and $\omega\neq 0$. Then for
$\gamma<0$ and $|\gamma|\geqslant 2|\omega|$, $\lambda_{1\pm}$ and
$\lambda_{2\pm}$ are all negative reals since

$\gamma+\sqrt{\gamma^{2}-4\omega^{2}}<0.$

Therefore, in this case, $J$ has $2N$ negative real eigenvalues, which implies
the origin is locally exponentially stable. For $|\gamma|>2|\omega|$, each
pair of roots is distinct from each other, so the oscillations are overdamped.
For $|\gamma|<2|\omega|$, if $\gamma<2ca$, $J$ has the same $N-1$ pairs of
complex conjugate eigenvalues with negative real parts ,
{$\lambda_{1+}$,$\lambda_{1-}$}, and a pair of complex conjugate eigenvalues,
{$\lambda_{2+}$,$\lambda_{2-}$}, which cross the imaginary axis due to a
variation of $\gamma$.

Note that

$\displaystyle\left.\text{Re}(\lambda_{2\pm})(\gamma)\right|_{\gamma=0}=$
$\displaystyle 0$
$\displaystyle\left.\frac{d\text{Re}(\lambda_{2\pm})(\gamma)}{d\gamma}\right|_{\gamma=0}=$
$\displaystyle\frac{1}{2}>0,$

where $\text{Re}(\lambda_{2\pm})$ denote the real parts of $\lambda_{2\pm}$.
According to the Hopf bifurcation theorem, a Hopf bifurcation arises from the
origin at $\gamma=0$, and the origin is locally exponentially stable for
$\gamma<0$ and unstable for $\gamma>0$. For $\gamma<0$, the oscillations are
underadmped since $J$ has $2N$ complex conjugate eigenvalues with negative
real parts in this case. If we further have $\alpha\geqslant 0$ and
$\beta\geqslant 0$, then the origin is globally asymptotically stable
according to Lemma 4.5, which means that there is no limit cycle around the
origin for $\gamma<0$. Thus, the Hopf bifurcation must be supercritical. ∎

By the Hopf bifurcation theorem for multi-dimensional systems, for
sufficiently small values of $\gamma\neq 0$ a continuous family of stable
limit cycles $\\{\xi(\gamma)\\}$, with the parameter $\gamma$, bifurcates from
the origin into the region $\gamma>0$ in a multi-dimensional state space when
the supercritical Hopf bifurcation occurs in the coupled-oscillator system
(24). As can be seen in Fig. 5 in Appendix C, for fixed values of $\alpha$ and
$\beta$ and a sufficiently small variable $\gamma>0$, the “projection” of the
limit cycle $\xi(\gamma)$ on each $x_{i}y_{i}$-plane is roughly a circle,
$\xi_{i}(\gamma)$, surrounding the origin, and all $\xi_{i}(\gamma)$ have the
same radius, $r(\gamma)$, that depends only on $\gamma$. We can put all
$\xi_{i}(\gamma)$ in the phase plane of the virtual system (25) of the network
(24). In fact, any two circles, $\xi_{i}(\gamma)$ and $\xi_{j}(\gamma)$, that
lie inside $C_{RvdP}$ must coincide, otherwise the projections,
$(x_{i}(t),y_{i}(t))$ and $(x_{j}(t),y_{j}(t))$, of the trajectory of (24)
that starts near $\xi(\gamma)$ would never converge. Therefore, coincidence of
$\xi_{i}(\gamma)$ for all $i$ is a necessary condition for the network (24)
that has a stable limit cycle $\xi(\gamma)$ and starts near it to be
completely synchronized. For a sufficiently small value of $\gamma$, the
origin is in fact the only equilibrium point of the system (24) with
$\alpha\geqslant 0$, $\beta\geqslant 0$, and $ca>0$, and the bifurcated limit
cycle $\xi(\gamma)$ is also unique. So $\xi(\gamma)$ is a globally stable
limit cycle, which means that every trajectory of (24) with $\alpha\geqslant
0$, $\beta\geqslant 0$, $ca>0$, and a sufficiently small $\gamma>0$ that does
not start at the origin will spiral into $\xi(\gamma)$ as time approaches
infinity.

When a globally stable limit cycle $\xi(\gamma)$ that evolves around the
origin exists in a coupled-oscillator system (24), it is plausible to assume
that if the projected trajectories $(x_{j}(t),y_{j}(t))$ of (24) start in a
disk of radius $R$ that contains $\xi_{j}(\gamma)$ for all $j$ from 1 to $N$,
then we have

$\sum_{j=1}^{N}x_{j}^{2}(t)+y_{j}^{2}(t)\leqslant NR^{2}.$ (30)

Then we can see from the proof of Lemma 4.6 that at any point, $(x,y)$, on the
boundary of the disk ${R_{RvdP}}_{\omega^{2}=1}$ of radius $r$ we have

$\begin{split}\nabla
H\cdot\langle\dot{x}(t),\dot{y}(t)\rangle=&2x\dot{x}(t)+2y\dot{y}(t)\\\
=&-2\left(\alpha x^{2}+\beta y^{2}-\gamma\right)y^{2}\\\
&+2\frac{ca}{N}(\sum_{j=1}^{N}(xx_{j}(t)+yy_{j}(t))-N(x^{2}+y^{2}))\\\
\leqslant&2\gamma(x^{2}+y^{2})+2\frac{ca}{N}(\sum_{j=1}^{N}(xx_{j}(t)+yy_{j}(t))-N(x^{2}+y^{2}))\\\
=&2\frac{ca}{N}(\sum_{j=1}^{N}rr_{j}(t)\cos\theta(t)-N(1-\frac{\gamma}{ca})r^{2})\\\
\leqslant&2\frac{ca}{N}(\sum_{j=1}^{N}rr_{j}(t)-N(1-\frac{\gamma}{ca})r^{2}).\end{split}$

Suppose that
${x^{2}_{j}}_{0}+{y^{2}_{j}}_{0}\leqslant(1-\frac{\gamma}{ca})^{2}r^{2}$ for
all $j$ and that $(1-\frac{\gamma}{ca})r>r(\gamma)$. By the Cauchy–Schwarz
inequality and Eq. (30), we have

$\frac{(\sum_{j=1}^{N}r_{j}(t))^{2}}{N}\leqslant\sum_{j=1}^{N}r_{j}^{2}(t)\leqslant
N(1-\frac{\gamma}{ca})^{2}r^{2}.$

Hence,

$\sum_{j=1}^{N}rr_{j}(t)\leqslant N(1-\frac{\gamma}{ca})r^{2}.$

Therefore, $\nabla H\cdot\langle\dot{x}(t),\dot{y}(t)\rangle\leqslant 0$ on
the boundary of ${R_{RvdP}}_{\omega^{2}=1}$. It follows that the disk,
${R_{RvdP}}_{\omega^{2}=1}$, having radius
$r>r(\gamma)(1-\frac{\gamma}{ca})^{-1}$ is a “vague” trapping region of the
virtual system (25) with $\alpha\geqslant 0$, $\beta\geqslant 0$, $ca>0$,
$\omega^{2}=1$, and with a sufficiently small value of $\gamma$ satisfying
$0<\gamma<2\min\\{ca,1\\}$ in the sence that every trajectory that starts
within a smaller disk of radius $(1-\frac{\gamma}{ca})r$ which is included in
${R_{RvdP}}_{\omega^{2}=1}$ will remain in ${R_{RvdP}}_{\omega^{2}=1}$ as the
system (25) evolves, and that the corresponding network (24) is completely
synchronized in the cases given in Lemma 4.7 if
${x^{2}_{i}}_{0}+{y^{2}_{i}}_{0}\leqslant(1-\frac{\gamma}{ca})^{2}r^{2}$ for
all $i$ from 1 to $N$ and $\sum_{i=1}^{N}{x^{2}_{i}}_{0}+{y^{2}_{i}}_{0}\neq
0$.

I do not have a proof for the inequality (30). However, this is not a crucial
problem if we do not limited to exponential synchronization because any
trajectory of (24) that starts everywhere except the origin would eventually
get infinitely close to $\xi(\gamma)$ on which (30) holds for $R\geqslant
r(\gamma)$.

Another problem is that the value of $r(\gamma)$ needs to be obatined
numerically. By an empirical rule of thumb for supercritical Hopf
bifurcations, the size of $\xi(\gamma)$ grows continuously from zero and
increases proportionally to $\sqrt{\gamma}$, so $r(\gamma)\approx
k\sqrt{\gamma}$. As you can see in Fig. 5, it can be found numerically that
the coefficient of proportionality $k$ depends on the nonlinear coefficients
$\alpha$ and $\beta$ in the way that $k$ increases with decrease in $\alpha$
or $\beta$. It is difficult to obtain a further relation of $k$ with $\alpha$
and $\beta$, but we can see from Fig. 5 (c) that $k\approx 1$ when
$\alpha=\beta=1$.

######  Example.

Consider a network of six all-to-all coupled identical Rayleigh van der Pol
oscillators with an identical full-state linear coupling described by Eq. (24)
with $\alpha=1$, $\beta=1$, $\gamma=0.6$, $\omega=1$, $a=1$, $c=1$, and $N=6$.
It can be seen from Fig. 6 (c) that for $\alpha=\beta=1$,
$r(\gamma)\approx\sqrt{\gamma}=\sqrt{0.6}$.

Choose the “vague” trapping region, ${R_{RvdP}}_{\omega^{2}=1}$, to be a disk
of radius $r=2>r(\gamma)(1-\frac{\gamma}{ca})^{-1}\approx 1.94$. The starting
points of this system are

$\displaystyle({x_{1}}_{0},{y_{1}}_{0})=(0.6,0.5),\;({x_{2}}_{0},{y_{2}}_{0})=(0.3,-0.2),$
$\displaystyle({x_{3}}_{0},{y_{3}}_{0})=(-0.7,-0.3),\;({x_{4}}_{0},{y_{4}}_{0})=(0.5,-0.5),$
$\displaystyle({x_{5}}_{0},{y_{5}}_{0})=(0.1,0.4),\;({x_{6}}_{0},{y_{6}}_{0})=(0,0).$

We have
${x^{2}_{i}}_{0}+{y^{2}_{i}}_{0}\leqslant(1-\frac{\gamma}{ca})^{2}r^{2}=0.64$
for all $i$ from 1 to 6, and it is obvious that
$\sum_{i=1}^{6}{x^{2}_{i}}_{0}+{y^{2}_{i}}_{0}\neq 0$.

By Lemma 4.7 (case 4 of the lemma), ${R_{RvdP}}_{\omega^{2}=1}$ is included in
${C_{RvdP}}_{\omega^{2}=1}$ since

$\displaystyle
4=r^{2}\geqslant\max\\{\frac{ca\alpha-3ca\beta}{\alpha^{2}},\frac{3ca\beta-
ca\alpha}{\alpha^{2}}\\}=2,$ $\displaystyle
4-2\sqrt{3.4}=\frac{ca\alpha+3ca\beta-2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}<r^{2}$
$\displaystyle<\frac{ca\alpha+3ca\beta+2\sqrt{3c^{2}a^{2}\alpha\beta+\alpha^{2}ca(ca-\gamma)}}{\alpha^{2}}=4+2\sqrt{3.4},$
$\displaystyle 3ca\alpha\beta+\alpha^{2}(ca-\gamma)=3.4>0.$

Then by the above analysis, we can predict that this network is completely
synchronized, which can be seen in Fig. 4.

(a) The time series of $x_{i}$ for $i$ from 1 to 6

(b) The time series of $y_{i}$ for $i$ from 1 to 6

(c) Phase portraits

Figure 4: The time series and phase portraits of six self-sustained all-to-all
coupled identical Rayleigh van der Pol oscillators with an identical full-
state linear coupling

Determining the direction of the Hopf bifurcation in the system (24) at
$\gamma=0$ is nontrivial when $\alpha$ or $\beta$ is negative. Moreover, the
large number of oscillators in the network makes it even more difficult. For
simplicity, we only consider the case of two coupled oscillators, i.e., $N=2$.
We use the techniques described in [20].

For a system of the form (1) with the bifurcation parameter $\nu$, let
$J_{\boldsymbol{f}}(\nu)$ be the Jacobian at the equilibrium point
$\boldsymbol{x}_{*}(\nu)$, i.e.,
$J_{\boldsymbol{f}}(\nu)=\frac{\partial\boldsymbol{f}(\boldsymbol{x};\nu)}{\partial{\boldsymbol{x}}}|_{\boldsymbol{x}=\boldsymbol{x}_{*}}$.
Suppose $J_{\boldsymbol{f}}(\nu)$ has $m$ pairs of complex conjugate
eigenvalues, $\\{\lambda_{i+},\lambda_{i-}\\}$, where $\lambda_{i+}$ have
positive imaginary parts and $\lambda_{i-}$ have negative ones, and $n-2m$
real eigenvalues, ${\lambda_{r}}_{j}$. Let
$\text{Re}(\lambda_{1+})\geqslant\text{Re}(\lambda_{2+})\geqslant\cdots\geqslant\text{Re}(\lambda_{m+})$
and let
${\lambda_{r}}_{1}\geqslant{\lambda_{r}}_{2}\geqslant\cdots\geqslant{\lambda_{r}}_{n-2m}$.
The system can be rewritten as
$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x};\nu)=J_{\boldsymbol{f}}(\nu)\boldsymbol{x}+\boldsymbol{h}(\boldsymbol{x})$,
where $\boldsymbol{h}(\boldsymbol{x})$ is the nonlinear term.

The Hopf bifurcation theorem gives that the Hopf bifurcation occurs at
$\nu_{c}$ where $\nu_{c}$ satiesfies $\text{Re}(\lambda_{1+}(\nu_{c}))=0$,
$\frac{d\text{Re}(\lambda_{1+}(\nu_{c}))}{d\nu}\neq 0$,
$\text{Im}(\lambda_{1+}(\nu_{c}))\neq 0$, $\text{Re}(\lambda_{i+}(\nu_{c}))<0$
for $i$ from 2 to $m$, and ${\lambda_{r}}_{j}(\nu_{c})<0$ for $j$ from 1 to
$n-2m$. Form a matrix
$P=(\text{Re}(\boldsymbol{v}_{1+})\;-\text{Im}(\boldsymbol{v}_{1+})\;\cdots\;\text{Re}(\boldsymbol{v}_{m+})\;-\text{Im}(\boldsymbol{v}_{m+})\;{\boldsymbol{v}_{r}}_{1}\;\cdots\;{\boldsymbol{v}_{r}}_{n-2m})$,
where $\boldsymbol{v}_{i+}$ and ${\boldsymbol{v}_{r}}_{j}$ are (generalized)
eigenvectors of $J_{\boldsymbol{f}}(\nu_{c})$ corresponding to
$\lambda_{i+}(\nu_{c})$ and ${\lambda_{r}}_{j}(\nu_{c})$ respectively. We have
the following lemma.

######  Lemma 4.11.

$P$ is invertable and $P^{-1}J_{\boldsymbol{f}}(\nu_{c})P$ is in the Jordan
canonical form, i.e.,

$P^{-1}J_{\boldsymbol{f}}(\nu_{c})P=\left(\begin{array}[]{r@{}c|c@{}l}&\begin{smallmatrix}0&-\text{Im}(\lambda_{1}(\nu_{c}))\\\
\text{Im}(\lambda_{1}(\nu_{c}))&0\rule[-3.01389pt]{0.0pt}{6.02777pt}\end{smallmatrix}&\mbox{\large
0}\\\ \hline\cr&\mbox{\large 0}&\mbox{D}\end{array}\right),$ (31)

where $D$ is a block diagonal matrix of original dimension $n-2$, whose
diagonals are Jordan blocks $D_{i}$ that are either of the form

$\left(\begin{array}[]{cccc}B_{i}&I_{2\times 2}&&\\\ &B_{i}&\ddots&\\\
&&\ddots&I_{2\times 2}\\\ &&&B_{i}\\\ \end{array}\right)$

for $i$ from 2 to $m$, where the empty elements are 0 blocks, $I_{2\times 2}$
is the two dimensional identical matrix and

$B_{i}=\left(\begin{array}[]{cc}\text{Re}(\lambda_{i+}(\nu_{c}))&-\text{Im}(\lambda_{i+}(\nu_{c}))\\\
\text{Im}(\lambda_{i+}(\nu_{c}))&\text{Re}(\lambda_{i+}(\nu_{c}))\\\
\end{array}\right)$

for complex eigenvalues $\lambda_{i+}(\nu_{c})$, or of the form

$\left(\begin{array}[]{cccc}{\lambda_{r}}_{j}(\nu_{c})&1&&\\\
&{\lambda_{r}}_{j}(\nu_{c})&\ddots&\\\ &&\ddots&1\\\
&&&{\lambda_{r}}_{j}(\nu_{c})\\\ \end{array}\right),$

for real eigenvalue ${\lambda_{r}}_{j}(\nu_{c})$ with $i$ from 1 to $n-2m$,
where the empty elements are 0.

Perforem coordinate transformations from $x$ to $y$ frames by
$\boldsymbol{x}=P\boldsymbol{y}+\boldsymbol{x}_{*}(\nu_{c})$. Then the system
at $\nu=\nu_{c}$ is converted into

$\begin{split}\dot{\boldsymbol{y}}=P^{-1}\dot{\boldsymbol{x}}=&P^{-1}\boldsymbol{f}(P\boldsymbol{y}+\boldsymbol{x}_{*}(\nu_{c}))\\\
=&P^{-1}J_{\boldsymbol{f}}(\nu_{c})P\boldsymbol{y}+P^{-1}J_{\boldsymbol{f}}(\nu_{c})\boldsymbol{x}_{*}(\nu_{c})+\boldsymbol{h}(P\boldsymbol{y}+\boldsymbol{x}_{*}(\nu_{c}))\\\
\equiv&\boldsymbol{F}(\boldsymbol{y}),\end{split}$

where $\boldsymbol{F}=\langle
F_{1},\cdots,F_{n}\rangle:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a vector-valued
function. Then
$\frac{\partial\boldsymbol{F}(\boldsymbol{y})}{\partial{\boldsymbol{y}}}|_{\boldsymbol{y}=\boldsymbol{0}}=P^{-1}J_{\boldsymbol{f}}(\nu_{c})P$.
Thus, by Lemma 4.10, the Jacobian for $\boldsymbol{y}$-system at
$\boldsymbol{y}=\boldsymbol{0}$ is in the Jordan canonical form. We calculate
the following quantities at $\boldsymbol{y}=\langle
y_{1},\cdots,y_{n}\rangle=\boldsymbol{0}$.

$\begin{split}g_{11}=&\frac{1}{4}[\frac{\partial^{2}F_{1}}{\partial{y_{1}}^{2}}+\frac{\partial^{2}F_{1}}{\partial{y_{2}}^{2}}+i(\frac{\partial^{2}F_{2}}{\partial{y_{1}}^{2}}+\frac{\partial^{2}F_{2}}{\partial{y_{2}}^{2}})],\\\
g_{20}=&\frac{1}{4}[\frac{\partial^{2}F_{1}}{\partial{y_{1}}^{2}}-\frac{\partial^{2}F_{1}}{\partial{y_{2}}^{2}}+2\frac{\partial^{2}F_{2}}{\partial
y_{1}\partial
y_{2}}+i(\frac{\partial^{2}F_{2}}{\partial{y_{1}}^{2}}-\frac{\partial^{2}F_{2}}{\partial{y_{2}}^{2}}-2\frac{\partial^{2}F_{1}}{\partial
y_{1}\partial y_{2}})],\\\
g_{02}=&\frac{1}{4}[\frac{\partial^{2}F_{1}}{\partial{y_{1}}^{2}}-\frac{\partial^{2}F_{1}}{\partial{y_{2}}^{2}}-2\frac{\partial^{2}F_{2}}{\partial
y_{1}\partial
y_{2}}+i(\frac{\partial^{2}F_{2}}{\partial{y_{1}}^{2}}-\frac{\partial^{2}F_{2}}{\partial{y_{2}}^{2}}+2\frac{\partial^{2}F_{1}}{\partial
y_{1}\partial y_{2}})],\\\
G_{21}=&\frac{1}{8}[\frac{\partial^{3}F_{1}}{\partial{y_{1}}^{3}}+\frac{\partial^{3}F_{1}}{\partial{y_{1}}\partial{y_{2}}^{2}}+\frac{\partial^{3}F_{2}}{\partial{y_{1}}^{2}\partial{y_{2}}}+\frac{\partial^{3}F_{2}}{\partial{y_{2}}^{3}}+i(\frac{\partial^{3}F_{2}}{\partial{y_{1}}^{3}}+\frac{\partial^{3}F_{2}}{\partial{y_{1}}\partial{y_{2}}^{2}}-\frac{\partial^{3}F_{1}}{\partial{y_{1}}^{2}\partial{y_{2}}}-\frac{\partial^{3}F_{1}}{\partial{y_{2}}^{3}})],\\\
g_{21}=&G_{21}+\sum_{i=1}^{n-2}(2G^{i}_{110}w^{i}_{11}+G^{i}_{101}w^{i}_{20}),\end{split}$

where

$\begin{split}G^{j-2}_{110}=&\frac{1}{2}[\frac{\partial^{2}F_{1}}{\partial
y_{1}\partial y_{j}}+\frac{\partial^{2}F_{2}}{\partial y_{2}\partial
y_{j}}+i(\frac{\partial^{2}F_{2}}{\partial y_{1}\partial
y_{j}}-\frac{\partial^{2}F_{1}}{\partial y_{2}\partial y_{j}})],\\\
G^{j-2}_{101}=&\frac{1}{2}[\frac{\partial^{2}F_{1}}{\partial y_{1}\partial
y_{j}}-\frac{\partial^{2}F_{2}}{\partial y_{2}\partial
y_{j}}+i(\frac{\partial^{2}F_{2}}{\partial y_{1}\partial
y_{j}}+\frac{\partial^{2}F_{1}}{\partial y_{2}\partial y_{j}})]\end{split}$

for $j=3,\cdots,n$, and where $\boldsymbol{w}_{11}=\langle
w^{1}_{11},\cdots,w^{n-2}_{11}\rangle$ and $\boldsymbol{w}_{20}=\langle
w^{1}_{20},\cdots,w^{n-2}_{20}\rangle$ are the solutions of

$\begin{split}&D\boldsymbol{w}_{11}=-\boldsymbol{h}_{11},\\\
(&D-2\text{Im}(\lambda_{1}(\nu_{c}))iI_{n-2})\boldsymbol{w}_{20}=-\boldsymbol{h}_{20},\end{split}$

respectively, with $\boldsymbol{h}_{11}=\langle
h^{1}_{11},\cdots,h^{n-2}_{11}\rangle$ and $\boldsymbol{h}_{20}=\langle
h^{1}_{20},\cdots,h^{n-2}_{20}\rangle$ whose elements are defined by

$\begin{split}h^{j-2}_{11}=&\frac{1}{4}(\frac{\partial^{2}F_{j}}{\partial{y_{1}}^{2}}+\frac{\partial^{2}F_{j}}{\partial{y_{2}}^{2}}),\\\
h^{j-2}_{20}=&\frac{1}{4}(\frac{\partial^{2}F_{j}}{\partial{y_{1}}^{2}}-\frac{\partial^{2}F_{j}}{\partial{y_{2}}^{2}}-2i\frac{\partial^{2}F_{j}}{\partial
y_{1}\partial y_{2}})\end{split}$

for $j=3,\cdots,n$, and with matrix $D$ defined in Eq. (31). We then calculate

$\beta_{2}=\text{Re}[\frac{i}{\text{Im}(\lambda_{1}(\nu_{c}))}(g_{20}g_{11}-2|g_{11}|^{2}-\frac{|g_{02}|^{2}}{3})+g_{21}].$

The limit cycle is stable (unstable) if $\beta_{2}<0$ ($\beta_{2}>0$) on the
side of $\nu=\nu_{c}$ where the limit cycle appears, which indicates the Hoft
bifurcation is supercritical (subcritical) [20].

For system (24) with $N=2$, we have

$\displaystyle\boldsymbol{f}(\boldsymbol{x};\gamma)=\left(\begin{array}[]{c}x_{3}+\frac{ca(x_{2}-x_{1})}{2}\\\
x_{4}+\frac{ca(x_{1}-x_{2})}{2}\\\
-(\alpha{x_{1}}^{2}+\beta{x_{3}}^{2}-\gamma)x_{3}-\omega^{2}x_{1}+\frac{ca(x_{4}-x_{3})}{2}\\\
-(\alpha{x_{2}}^{2}+\beta{x_{4}}^{2}-\gamma)x_{4}-\omega^{2}x_{2}+\frac{ca(x_{3}-x_{4})}{2}\\\
\end{array}\right),$ $\displaystyle
P=\left(\begin{array}[]{cccc}0&{|\omega|}^{-1}&0&-{|\omega|}^{-1}\\\
0&{|\omega|}^{-1}&0&{|\omega|}^{-1}\\\ 1&0&-1&0\\\ 1&0&1&0\\\
\end{array}\right),$ $\displaystyle g_{11}=g_{20}=g_{02}=0,$ $\displaystyle
G_{21}=-\frac{\alpha+3\beta\omega^{2}}{4\omega^{2}},$
$\displaystyle\boldsymbol{h}_{11}=\boldsymbol{h}_{20}=\boldsymbol{0},$
$\displaystyle\boldsymbol{w}_{11}=\boldsymbol{w}_{20}=\boldsymbol{0},$
$\displaystyle g_{21}=G_{21},$

and

$\beta_{2}=-\frac{\alpha+3\beta\omega^{2}}{4\omega^{2}}.$

If $\omega\neq 0$ and $\alpha+3\beta\omega^{2}>0$ ($<0$), then $\beta_{2}<0$
($>0$), and thus a supercritical (subcritical) Hopf bifurcation occurs. For
systems with more oscillators, we can analyze the direction of the Hopf
bifurcation in the same way, which will not be discussed in this paper.
Combined with Lemma 4.9, we can conclude

######  Theorem 4.12.

For a two-oscillator system described by Eq. (24) with $ca>0$ and $N=2$, when
$|\gamma|<2\omega$ and $\gamma<2ca$, the origin is locally exponentially
stable for $\gamma<0$ and unstable for $\gamma>0$, and a continuous one-
parameter family of stable (unstable) limit cycles bifurcate from the origin
into the region $\gamma>0$ ($<0$) if $\alpha+3\beta\omega^{2}>0$ ($<0$).

In the case of a supercritical Hoft bifurcation, the conclusion about the
conditions for synchronization of the network (24) with $\alpha,\beta\geqslant
0$ also applies here for the case $\alpha+3\beta\omega^{2}>0$ with negative
value of $\alpha$ or $\beta$. However, there are two issues to be noted, one
is that the shape of the limit cycle becomes distorted more rapidly from a
circle with increasing $\gamma$ compared to the case for
$\alpha,\beta\geqslant 0$, and the other is that the stability of the limit
cycle becomes local since the origin may not be a unique equilibrium point
even for a very small value of $\gamma$ when $\alpha$ or $\beta$ is negative.
Therefore, in the case where either $\alpha<0$ or $\beta<0$, the radius of the
“vague” trapping region ${R_{RvdP}}_{\omega^{2}=1}$ of the virtual system (25)
with $ca>0$, $\omega^{2}=1$, and $\alpha+3\beta>0$, can be taken to be only
slightly greater than $r(\gamma)(1-\frac{\gamma}{ca})^{-1}$.

##  5 Conclusion

In this paper, the contraction theory is introduced to analyze the complete
synchronization of nonlinear coupled oscillators. It can be seen that on the
assumptions of identical coupled oacillators and linear couplings, we can find
the contraction regions of the virtual systems of the coupled-oscillatot
networks in different constructions, including all-to-all topology, star
topology, and ring topology. For linearly coupled van der Pol and Duffing van
der Pol oscillators, the contraction regions cover the entire phase plane of
the corresponding virtual systems when the coefficients of the networks meet
certain conditions, thus allowing all oscillatos in the networks to be
completely synchronized under arbitrary initial conditions. In contrast, the
contraction region of the virtual system of a network of the linearly coupled
Rayleigh van der Pol oscillators is always finite, but we can find a trapping
region and the conditions that make it included in the contraction region so
that the network is completely synchronized if the trajectories of all
oscillators in the network start within such a trapping region. Contraction
theory has so far been applied to the study of complete synchronization and
anti-synchronization of coupled oscillators. To show the anti-synchronization
of two oscillators with states $\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{2}$,
respectively, we only need to prove that the oscillator with state
$\boldsymbol{x}_{1}$ and an oscillator with state $-\boldsymbol{x}_{2}$ are
completely synchronized. The future research direction could be to apply the
contraction theory to phase synchronization and phase-locking phenomena and to
combine it with other traditional methods such as the Kuramoto model and the
Hilbert transform method.

## Acknowledgement

I remain indebted to my professors Randy Tagg and Masoud Asadi-Zeydabadi for
giving me the golden opportunity to do this ressearch on the synchronization
of coupled oscillators and for their valuable inputs in the early stages of
this study.

## References

  * [1] Mirollo, R.E. and Strogatz, S.H., 1990. Synchronization of pulse-coupled biological oscillators. SIAM Journal on Applied Mathematics, 50(6), pp.1645-1662.
  * [2] Roy, P.K., Hens, C., Grosu, I. and Dana, S.K., 2011. Engineering generalized synchronization in chaotic oscillators. Chaos: An Interdisciplinary Journal of Nonlinear Science, 21(1), p.013106.
  * [3] Golubitsky, M., Stewart, I., Buono, P.L. and Collins, J.J., 1998. A modular network for legged locomotion. Physica D: Nonlinear Phenomena, 115(1-2), pp.56-72.
  * [4] Winfree, A.T., 1967. Biological rhythms and the behavior of populations of coupled oscillators. SIAM Journal on Applied Mathematics, 50(6), pp.1645-1662.
  * [5] Pikovsky, A., Kurths, J., Rosenblum, M. and Kurths, J., 2003. Synchronization: a universal concept in nonlinear sciences (No. 12). Cambridge university press.
  * [6] Collins, J., 2003. Sync: The Emerging Science of Spontaneous Order. Nature, 422(6928), pp.117-118.
  * [7] Mirollo, R.E. and Strogatz, S.H., 1990. Synchronization transitions in a disordered Josephson series array. Physical review letters, 76(3), p.404.
  * [8] Cross, M.C., Zumdieck, A., Lifshitz, R. and Rogers, J.L., 2004. Synchronization by nonlinear frequency pulling. Physical review letters, 93(22), p.224101.
  * [9] Ermentrout, B., 1991. An adaptive model for synchrony in the firefly Pteroptyx malaccae. Journal of Mathematical Biology, 29(6), pp.571-585.
  * [10] Jiang, Y., 2012. Modeling Accoustic Synchronization of Crickets.
  * [11] Javaloyes, J., Perrin, M. and Politi, A., 2008. Collective atomic recoil laser as a synchronization transition. Physical Review E, 78(1), p.011108.
  * [12] Montbrió, E. and Pazó, D., 2018. Kuramoto model for excitation-inhibition-based oscillations. Physical review letters, 120(24), p.244101.
  * [13] Xu, C., Gao, J., Sun, Y., Huang, X. and Zheng, Z., 2015. Explosive or Continuous: Incoherent state determines the route to synchronization. Scientific reports, 5(1), pp.1-11.
  * [14] Simpson-Porco, J.W. and Bullo, F., 2014. Contraction theory on Riemannian manifolds. Systems $\&$ Control Letters, 65, pp.74-80.
  * [15] Lohmiller, W. and Slotine, J.J.E., 1998. On contraction analysis for non-linear systems. Automatica, 34(6), pp.683-696.
  * [16] Wang, W. and Slotine, J.J.E., 2005. On partial contraction analysis for coupled nonlinear oscillators. Biological cybernetics, 92(1), pp.38-53.
  * [17] Ruzitalab, A., Farahi, M.H. and Erjaee, G., 2018. Partial contraction analysis of coupled fractional order systems. Journal of Applied Mathematics, 2018.
  * [18] Gallier, J., 2020. The Schur complement and symmetric positive semidefinite (and definite) matrices (2019). URL https://www. cis. upenn. edu/jean/schur-comp. pdf.
  * [19] Alderisio, F., Bardy, B.G. and Di Bernardo, M., 2016. Entrainment and synchronization in networks of Rayleigh$--$van der Pol oscillators with diffusive and Haken$--$Kelso$--$Bunz couplings. Biological cybernetics, 110(2), pp.151-169.
  * [20] Hassard, B.D., Hassard, B.D., Kazarinoff, N.D., Wan, Y.H. and Wan, Y.W., 1981. Theory and applications of Hopf bifurcation (Vol. 41). CUP Archive.

## Appendix

##  Appendix A Proof of Theorem 2.1

Consider a nonlinear system of the form (1). Let $\boldsymbol{x}_{0}$ and
$\boldsymbol{x}^{\prime}_{0}$ be two different but close initial points of the
system on the state space. Assume that the state space of the system is a flat
space equipped with Cartesian coordinates.

Let
$\delta\boldsymbol{x}=\boldsymbol{x}\left(t;\boldsymbol{x}_{0}\right)-\boldsymbol{x}^{\prime}\left(t;\boldsymbol{x}^{\prime}_{0}\right)$
be an infinitesimal displacement between two different trajectories of the
system at the same time. By chain rule, it can be seen from Eq. (1) that

$\delta\dot{\boldsymbol{x}}=\delta\boldsymbol{f}(\boldsymbol{x}(t))=\frac{\partial\boldsymbol{f}(\boldsymbol{x}(t))}{\partial\boldsymbol{x}}\delta\boldsymbol{x}.$

It implies that the rate of change of the squared “distance” between two
trajectories $\delta s^{2}=\delta\boldsymbol{x}^{T}\delta\boldsymbol{x}$ in
the flat state space is:

$\displaystyle\frac{d\left(\delta\boldsymbol{x}^{T}\delta\boldsymbol{x}\right)}{dt}=2\delta\boldsymbol{x}^{T}\delta\dot{\boldsymbol{x}}=2\delta\boldsymbol{x}^{T}\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}\delta\boldsymbol{x}=\delta\boldsymbol{x}^{T}\left(\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}+\left(\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}\right)^{T}\right)\delta\boldsymbol{x}$
$\displaystyle=2\delta\boldsymbol{x}^{T}\left(\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}\right)_{s}\delta\boldsymbol{x}\leq
2\lambda_{\max}(\boldsymbol{x})\delta\boldsymbol{x}^{T}\delta\boldsymbol{x}.$

where $\lambda_{\max}(\boldsymbol{x})$ is the largest eigenvalue of
$\left(\frac{\partial\boldsymbol{f}(\boldsymbol{x})}{\partial\boldsymbol{x}}\right)_{s}$.
We therefore have

$\frac{d\left(\delta\boldsymbol{x}^{T}\delta\boldsymbol{x}\right)}{\delta\boldsymbol{x}^{T}\delta\boldsymbol{x}}\leqslant
2\lambda_{\max}(\boldsymbol{x}(t))dt.$

Integrating both sides from $0$ to $t$ gives

$\ln\left[\delta s^{2}\right]-\ln\left[\delta s^{2}_{t=0}\right]\leqslant
2\int_{0}^{t}\lambda_{\max}(\boldsymbol{x}(\tau))\,d\tau,$

or,

$\displaystyle\delta s^{2}\leqslant$ $\displaystyle\delta
s^{2}_{t=0}\exp\left[2\int_{0}^{t}\lambda_{\max}(\boldsymbol{x}(\tau))\right]$
$\displaystyle\delta s\leqslant$ $\displaystyle\delta
s_{t=0}\exp\left[\int_{0}^{t}\lambda_{\max}(\boldsymbol{x}(\tau))\right].$

Since $\lambda_{\max}(\boldsymbol{x}(t))$ is uniformly negative in the
contraction region $C$, the infinitesimal distance between any two
trajectories that are confined in $C$ exponentially converges to zero. Using
path integration, any two trajectories of the system that remain in $C$
converge exponentially to each other.

##  Appendix B Proof of Lemma 4.5

Note that the origin is an equilibrium point. Since $\omega\neq 0$, we can
construct a continuously differentiable function
$V(x_{i},y_{i})=\sum_{i=1}^{N}(\frac{{x_{i}}^{2}}{2}+\frac{{y_{i}}^{2}}{2\omega^{2}})$.
This is a Lyapunov function for system (24). To see this, we take the
derivative of $V$ with respect to time along an arbitrary trajectory of Eq.
(24).

$\begin{split}\dot{V}(x_{i},y_{i})=&\sum_{i=1}^{N}(x_{i}\dot{x}_{i}+\frac{y_{i}\dot{y}_{i}}{\omega^{2}})\\\
=&\sum_{i=1}^{N}(x_{i}y_{i}-\frac{ca}{N}x_{i}\sum_{j=1}^{N}(x_{i}-x_{j})\\\
&-\frac{y_{i}}{\omega^{2}}\left(\alpha{x_{i}}^{2}+\beta{y_{i}}^{2}-\gamma\right)y_{i}-\frac{y_{i}}{\omega^{2}}\omega^{2}x_{i}-\frac{ca}{N\omega^{2}}y_{i}\sum_{j=1}^{N}(y_{i}-y_{j}))\\\
\leqslant&\frac{\gamma}{\omega^{2}}\sum_{i=1}^{N}{y_{i}}^{2}-\frac{ca}{N}\sum_{i,j=1}^{N}x_{i}(x_{i}-x_{j})-\frac{ca}{N\omega^{2}}\sum_{i,j=1}^{N}y_{i}(y_{i}-y_{j})\\\
=&\frac{\gamma}{\omega^{2}}\sum_{i=1}^{N}{y_{i}}^{2}-\frac{ca}{2N}\sum_{i,j=1}^{N}({x_{i}}^{2}-2x_{i}x_{j}+{x_{i}}^{2})-\frac{ca}{2N\omega^{2}}\sum_{i,j=1}^{N}({y_{i}}^{2}-2y_{i}y_{j}+{y_{i}}^{2}).\end{split}$

Notice that $\sum_{i,j=1}^{N}{x_{i}}^{2}=\sum_{j,i=1}^{N}{x_{j}}^{2}$ and that
$\sum_{i,j=1}^{N}{y_{i}}^{2}=\sum_{j,i=1}^{N}{y_{j}}^{2}$. We therefore have

$\begin{split}\dot{V}(x_{i},y_{i})\leqslant&\frac{\gamma}{\omega^{2}}\sum_{i=1}^{N}{y_{i}}^{2}-\frac{ca}{2N}\sum_{i,j=1}^{N}({x_{i}}^{2}-2x_{i}x_{j}+{x_{j}}^{2})-\frac{ca}{2N\omega^{2}}\sum_{i,j=1}^{N}({y_{i}}^{2}-2y_{i}y_{j}+{y_{j}}^{2})\\\
\leqslant&\frac{\gamma}{\omega^{2}}\sum_{i=1}^{N}{y_{i}}^{2}-\frac{ca}{2N}\sum_{i,j=1}^{N}(x_{i}-x_{j})^{2}-\frac{ca}{2N\omega^{2}}\sum_{i,j=1}^{N}(y_{i}-y_{j})^{2}\\\
\leqslant&\frac{\gamma}{\omega^{2}}\sum_{i=1}^{N}{y_{i}}^{2}\leqslant
0.\end{split}$

Since $\gamma<0$, the equality holds if and only if
$\boldsymbol{x}(t)=(x_{i}(t)\text{ }y_{i}(t))^{T}=\boldsymbol{0}$ for
$t\geqslant 0$. Note that

$\begin{split}V(\boldsymbol{0})=&0\\\ V(\boldsymbol{x})>&0\text{,
}\forall\boldsymbol{x}\neq\boldsymbol{0}\\\
V(\boldsymbol{x})\rightarrow&\infty\text{, as
}\|x\|\rightarrow\infty\end{split}$

Thus, by the LaSalle’s invariance principle, the origin is globally
asymptotically stable.

##  Appendix C The projections of limit cycles

(a) $\alpha=0.8$, $\beta=1$: $r(\gamma)\approx 1.03\sqrt{\gamma}$

(b) $\alpha=1$, $\beta=0.8$: $r(\gamma)\approx 1.08\sqrt{\gamma}$

(c) $\alpha=1$, $\beta=1$: $r(\gamma)\approx\sqrt{\gamma}$

(d) $\alpha=1.2$, $\beta=1$: $r(\gamma)\approx 0.973\sqrt{\gamma}$

(e) $\alpha=1$, $\beta=1.2$: $r(\gamma)\approx 0.935\sqrt{\gamma}$

Figure 5: The projections of stable limit cycles $\xi(\gamma)$ of six-coupled-
oscillator system (24) with $c=a=\omega=1$ and with different combinations of
$\alpha$, $\beta$, and $\gamma$ on $x_{i}y_{i}$-planes for all $i$ from 1 to
6. It can be seen that for fixed $\alpha$ and $\beta$, all projections are
roughly circles of the same radius $r(\gamma)$ that is directly proportional
to $\sqrt{\gamma}$ for small $\gamma>0$, and that the constant of
proportionality increases with decrease in $\alpha$ or $\beta$.