# Matrix coupling and generalized frustration in Kuramoto oscillators

Guilhermo L. Buzanello  Ana Elisa D. Barioni  Marcus A. M. de Aguiar
Instituto de Física ‘Gleb Wataghin’, Universidade Estadual de Campinas,
Unicamp 13083-970, Campinas, SP, Brazil

###### Abstract

The Kuramoto model describes the synchronization of coupled oscillators that
have different natural frequencies. Among the many generalizations of the
original model, Kuramoto and Sakaguchi (KS) proposed a frustrated version that
resulted in dynamic behavior of the order parameter, even when the average
natural frequency of the oscillators is zero. Here we consider a
generalization of the frustrated KS model that exhibits new transitions to
synchronization. The model is identical in form to the original Kuramoto
model, but written in terms of unit vectors. Replacing the coupling constant
by a coupling matrix breaks the rotational symmetry and forces the order
parameter to point in the direction of the eigenvector with highest
eigenvalue, when the eigenvalues are real. For complex eigenvalues the module
of order parameter oscillates while it rotates around the unit circle,
creating active states. We derive the complete phase diagram for Lorentzian
distribution of frequencies using the Ott-Antonsen ansatz. We also show that
changing the average value of the natural frequencies leads to further phase
transitions where the module of the order parameter goes from oscillatory to
static.

##  I Introduction

In 1975 Kuramoto proposed a simple model of $N$ coupled oscillators that could
be solved analytically in the limit where $N$ goes to infinity Kuramoto [1975,
1984]. The oscillators are described by their phases $\theta_{i}$ and are
coupled according to the equations

$\dot{\theta}_{i}=\omega_{i}+\frac{K}{N}\sum_{j=1}^{N}\sin{(\theta_{j}-\theta_{i})}$
(1)

where $\omega_{i}$ are their natural frequencies, selected from a symmetric
distribution $g(\omega)$, $K$ is the coupling strength and $i=1,...,N$. The
complex order parameter

$z=pe^{i\psi}\equiv\frac{1}{N}\sum_{i=1}^{N}e^{i\theta_{i}}$ (2)

measures the degree of phase synchronization of the particles: disordered
motion implies $p\approx 0$ and coherent motion $p\approx 1$. The phase of $z$
has been largely ignored, as it converges to a constant that depends on the
initial conditions. Random initial conditions lead, therefore, to random
values of $\psi$ for each simulation.

Kuramoto showed that the onset of synchronization could be described, in
equilibrium, as a continuous phase transition, where $p$ remains very small
for $0<K<K_{c}$ and increases as $p=\sqrt{1-K_{c}/K}$ for $K>K_{c}$. Since
then, a large number of modifications and generalizations of his model have
been proposed, including different types of coupling functions Hong and
Strogatz [2011], Yeung and Strogatz [1999], Breakspear et al. [2010],
introduction of networks of connections (so that not all oscillators are
connected to each other) Rodrigues et al. [2016], Climaco and Saa [2019],
different distributions of the oscillator’s natural frequencies (including
frequencies proportional to the number of connections, leading to explosive
synchronization) Gomez-Gardenes et al. [2011], Ji et al. [2013], inertial
terms Acebrón et al. [2005], Dörfler and Bullo [2011], Olmi et al. [2014],
external periodic driving forces Childs and Strogatz [2008], Moreira and de
Aguiar [2019a, b] and coupling with particle swarms O’Keeffe et al. [2017,
2022].

More recently, interest has been shifted to understand oscillations in larger
dimensions. Chandra et al Chandra et al. [2019a] have shown that Kuramoto
oscillators could also be described by unit vectors
$\vec{\sigma_{i}}=(\cos{\theta_{i}},\sin{\theta_{i}})$ that rotate on the unit
circle. It is easy to show that if $\theta_{i}$ satisfies Kuramoto’s equation
(1) then

$\frac{d\vec{\sigma_{i}}}{dt}=\mathbf{W}_{i}\vec{\sigma_{i}}+\frac{K}{N}\sum_{j}[\vec{\sigma_{j}}-(\vec{\sigma_{i}}\cdot\vec{\sigma_{j}})\vec{\sigma_{i}}]$
(3)

where $\mathbf{W}_{i}$ is an anti-symmetric matrix containing the natural
frequency $\omega_{i}$:

$\mathbf{W}_{i}=\left(\begin{array}[]{cc}0&-\omega_{i}\\\
\omega_{i}&0\end{array}\right).$ (4)

The complex order parameter $z$, Eq.(2), is replaced by the vector

$\vec{p}=\frac{1}{N}\sum_{i}\vec{\sigma_{i}}=(p\cos\psi,p\sin\psi)$ (5)

describing the center of mass of the system. Eq.(3) can be naturally extended
to higher dimensions by simply considering unit vectors $\vec{\sigma}_{i}$ in
D-dimensions rotating on the surface of the corresponding (D-1)-sphere. It has
been shown, in particular, that the system exhibits discontinuous phase
transitions in odd dimensions and continuous transitions in even dimensions
Chandra et al. [2019a]. Also, complexity reduction similar to that proposed by
Ott and Antonsen Ott and Antonsen [2008] can be applied in any dimension
Chandra et al. [2019b], Barioni and de Aguiar [2021a, b].

In this paper we explore the original two-dimensional case (Section II)
generalizing the coupling constant $K$ to a $2\times 2$ matrix K that can be
separated into a rotation plus a symmetric matrix. The corresponding Kuramoto
equations are then written as a generalization of the Kuramoto-Sakaguchi
model. This allows us to construct a complete phase diagram (Section III)
displaying three different regions: no synchrony, phase tuned states and
active states. In addition, varying the average frequency of the oscillators
we find conditions for the stability of these state. Finally, we compare our
analytical model to the numerical simulations (Section IV).

##  II Matrix coupling

Kuramoto’s model in vector form can be naturally extended if the coupling
constant $K$ in Eq.(3) is replaced by a $2\times 2$ matrices
${\mathbf{K}}_{ij}$ with elements that might depend also on the angles
$\phi_{i}$ and $\phi_{j}$. The equations now read

$\frac{d\vec{\sigma_{i}}}{dt}=\mathbf{W}_{i}\vec{\sigma_{i}}+\frac{1}{N}\sum_{j}[{\mathbf{K}}_{ij}\vec{\sigma_{j}}-(\vec{\sigma_{i}}\cdot{\mathbf{K}}_{ij}\vec{\sigma_{j}})\vec{\sigma_{i}}]$
(6)

and can be interpreted as a generalized frustrated model, as $\mathbf{K}_{ij}$
rotates $\vec{\sigma}_{j}$ hindering its alignment with $\vec{\sigma}_{i}$ and
inhibiting synchronization.

Defining the auxiliary vectors
$\vec{q}_{i}=\sum_{j}{\mathbf{K}}_{ij}\,\vec{\sigma}_{j}$ we can rewrite
theses equations as

$\frac{d\vec{\sigma_{i}}}{dt}=\mathbf{W}_{i}\vec{\sigma_{i}}+[\vec{q}_{i}-(\vec{\sigma_{i}}\cdot\vec{q}_{i})\vec{\sigma_{i}}].$
(7)

Norm conservation, $|\vec{\sigma_{i}}|=1$, is guaranteed for any set of
regular matrices ${\mathbf{K}}_{ij}$, as can be seen by taking the scalar
product of Eqs.(6) or (7) with $\vec{\sigma_{i}}$. In this paper we shall only
consider the case where ${\mathbf{K}}_{ij}={\mathbf{K}}=constant$ and
$\vec{q}_{i}=\vec{q}=\mathbf{K}\vec{p}$. Writing

$\mathbf{K}=\left(\begin{array}[]{cc}a&b\\\ c&d\end{array}\right)$ (8)

Eq.(6) can be written back in terms of the phases $\theta_{i}$ as

$\dot{\theta}_{i}=\omega_{i}+\frac{1}{2N}\sum_{j=1}^{N}\left[(a+d)\sin{\theta_{ji}^{-}}+(d-a)\sin{\theta_{ji}^{+}}+(c-b)\cos{\theta_{ji}^{-}}+(b+c)\cos{\theta_{ji}^{+}}\right]$
(9)

where ${\theta_{ji}^{-}}=\theta_{j}-\theta_{i}$ and
${\theta_{ji}^{+}}=\theta_{j}+\theta_{i}$. From this representation we
immediately recognize two special cases: (i) $a=d$, $b=c=0$, corresponding to
the usual Kuramoto model and, (ii) $a=d=K\cos\alpha$, $b=-c=K\sin\alpha$ (a
rotation matrix). In this case the terms inside the brackets in Eq.(9)
simplify to $\sin(\theta_{j}-\theta_{i}-\alpha)$, corresponding to the
Kuramoto-Sakaguchi model Sakaguchi and Kuramoto [1986], Yue et al. [2020].
This motivates us to re-parametrize $\mathbf{K}$ as

$\mathbf{K}=K\left(\begin{array}[]{cc}\cos\alpha&\sin\alpha\\\
-\sin\alpha&\cos\alpha\end{array}\right)+J\left(\begin{array}[]{cc}-\cos\beta&\sin\beta\\\
\sin\beta&\cos\beta\end{array}\right)\equiv\mathbf{K}_{R}+\mathbf{K}_{S}$ (10)

where $\mathbf{K}_{R}$ is a rotation and $\mathbf{K}_{S}$ a symmetric matrix.
In these variables Eq.(9) reads

$\dot{\theta}_{i}=\omega_{i}+\frac{1}{N}\sum_{j=1}^{N}\left[K\sin(\theta_{j}-\theta_{i}-\alpha)+J\sin(\theta_{j}+\theta_{i}+\beta)\right].$
(11)

and the eigenvalues of $\mathbf{K}$ are given by

$\lambda_{\pm}=K\cos\alpha\pm\sqrt{J^{2}-K^{2}\sin^{2}\alpha}$ (12)

and are independent of $\beta$.

###  II.1 Continuity equation

In the limit of infinitely many oscillators we define $f(\omega,\theta,t)$ as
the density of oscillators with natural frequency $\omega$ at position
$\theta$ in time $t$. It satisfies the continuity equation

$\frac{\partial f}{\partial t}+\frac{\partial(fv_{\theta})}{\partial\theta}=0$
(13)

with velocity field

$\vec{v}=\vec{\omega}\times\hat{r}+\vec{q}-(\hat{r}\cdot\vec{q})\hat{r}=(\omega+q_{\theta})\hat{\theta}\equiv
v_{\theta}\hat{\theta}$ (14)

where $\vec{\omega}$ points perpendicular to the plane of rotation. Eq.(5) for
the order parameter becomes

$\vec{p}(t)=\int\hat{r}(\theta)f(\omega,\theta,t)d\theta\,d\omega.$ (15)

###  II.2 OA ansatz for density function

The density of oscillators is a periodic function of $\theta$ and, therefore,
can be expanded in Fourier series. Ott & Antonsen Ott and Antonsen [2008]
showed that if Fourier coefficients are chosen as $\rho^{|m|}e^{-im\phi}$ the
solution is self-consistent, in the sense that it preserves this form at all
times, remaining in this restricted subset of density functions. The ansatz,
parametrized by $\rho(\omega,t)$ and $\phi(\omega,t)$ is:

$\displaystyle
f(\omega,\theta,t)=\frac{g(\omega)}{2\pi}\sum_{m=-\infty}^{\infty}\rho^{|m|}e^{im(\theta-\phi)}=\frac{g(\omega)}{2\pi}\frac{(1-\rho^{2})}{1+\rho^{2}-2\rho\cos(\theta-\phi)}.$
(16)

It is convenient to define the vector $\vec{\rho}=\rho(\cos\phi,\sin\phi)$ so
that $\cos(\theta-\phi)=\hat{r}\cdot\hat{\rho}$. Substituting Eq.(16) into
(15) we find:

$\vec{p}(t)=\int\vec{\rho}(\omega,t)g(\omega)\,d\omega.$ (17)

Also, substituting Eq.(16) into Eq.(13) we obtain Barioni and de Aguiar
[2021a]

$\displaystyle\dot{\vec{\rho}}=\vec{\omega}\times\vec{\rho}+\frac{1}{2}(1+\rho^{2}){\mathbf{K}}\vec{p}-({\mathbf{K}}\vec{p}\cdot\vec{\rho})\vec{\rho}.$
(18)

###  II.3 Complex variables

The Ott-Antonsen ansatz is particularly useful for Lorentzian distributions
$g(\omega)$. In order to derive the equation satisfied by the order parameter
in this case we define the complex variables

$z=pe^{i\psi},\qquad w=\rho e^{i\phi},\qquad u=qe^{i\xi}.$ (19)

Eq.(18) then becomes

$\dot{w}=i\omega w+\frac{u}{2}-\frac{u^{*}}{2}w^{2}.$ (20)

where

$u=Kze^{-i\alpha}-Jz^{*}e^{-i\beta}$ (21)

is obtained from the definition $\vec{q}=\mathbf{K}\vec{p}$. Substituting (21)
into (20) we obtain

$\dot{w}=i\omega
w+\frac{1}{2}\left(Kze^{-i\alpha}-Jz^{*}e^{-i\beta}\right)-\frac{w^{2}}{2}\left(Kz^{*}e^{i\alpha}-Jze^{i\beta}\right).$
(22)

Equation (17), on the other hand becomes

$z=\int w(\omega)g(\omega)d\omega.$ (23)

We note that if we define $F=-Jz^{*}e^{-i\beta}$, Eq.(22) can be written as

$\dot{w}=i\omega
w+\frac{1}{2}\left(Kze^{-i\alpha}+F\right)-\frac{w^{2}}{2}\left(Kze^{-i\alpha}+F\right)^{*}$
(24)

which is identical to the forced Kuramoto model discussed in Childs and
Strogatz [2008], with the difference that now the ‘force’ $F$ is generated by
the system itself, instead of being applied externally.

##  III Dynamics and phase diagram

For the Lorentzian distribution

$g(\omega)=\frac{1}{\pi}\frac{\Delta}{(\omega-\omega_{0})^{2}+\Delta^{2}}$
(25)

Eq.(23) can be integrated in the complex plane using a close contour from
$\omega=-R$ to $\omega=R$ and back to $-R$ through the upper half circle
$\omega=Re^{i\phi}$, enclosing the pole at $\omega=\omega_{0}+i\Delta$. Taking
$R\rightarrow\infty$ we obtain $z=w(\omega_{0}+i\Delta)$ Ott and Antonsen
[2008]. Calculating Eq.(22) at $\omega_{0}+i\Delta$ we can replace $w$ by $z$
to get

$\dot{z}=i(\omega_{0}+i\Delta)z+\frac{1}{2}\left(Kze^{-i\alpha}-Jz^{*}e^{-i\beta}\right)-\frac{z^{2}}{2}\left(Kz^{*}e^{i\alpha}-Jze^{i\beta}\right).$
(26)

Separating real and imaginary parts we obtain

$\dot{p}=-\Delta
p+\frac{p}{2}\left(1-p^{2}\right)\left[K\cos\alpha-J\cos(2\psi+\beta)\right]$
(27)

and

$\dot{\psi}=\omega_{0}-\frac{1}{2}\left(1+p^{2}\right)\left[K\sin\alpha-J\sin(2\psi+\beta)\right].$
(28)

###  III.1 Equilibrium solutions

Strict equilibrium solutions of Eqs.(27) and (28) are only possible if
$J=\alpha=\omega_{0}=0$. In this case $\psi$ is constant and Eq.(27) is
autonomous. However, a less stringent definition of equilibrium can be
obtained by setting only $J=0$. In this case Eq.(27) has two stationary
solutions, $p=0$ and $p=\sqrt{1-K_{c}/K}$, where $K_{c}\equiv
2\Delta/\cos\alpha$. The latter bifurcates at $K=K_{c}$ and corresponds to the
synchronized state of the Kuramoto-Sakaguchi model. At this solution $\vec{p}$
rotates with constant angular velocity
$\dot{\psi}=\omega_{0}-K\sin\alpha+\Delta\tan\alpha$. It is then possible to
change to a frame of reference that rotates with $\psi$, making the system
stationary. This is possible because $\mathbf{K}$ is itself a rotation matrix
and commutes with the operation that changes reference frames.

###  III.2 Phase tuned states - real eigenvalues

If $|J|>|K|\sin\alpha$ the eigenvalues of $\mathbf{K}$ are real (see Eq.(12)).
Let $\vec{v}=(\cos\gamma,\sin\gamma)$ be the eigenvector corresponding to the
largest eigenvalue of $\mathbf{K}$. From
$\mathbf{K}\vec{v}=\lambda_{+}\vec{v}$ it follows that
$K\sin\alpha-J\sin(2\gamma+\beta)=0$. This implies that when $\vec{p}$ is an
eigenvector of $\mathbf{K}$ ($\psi=\gamma$) we find $\dot{\psi}=\omega_{0}$.
We call this region phase tuned, as the phase $\psi$ can be tuned with the
choice of $\mathbf{K}$. Eq. (7) then seems to recover the Kuramoto model with
$\lambda_{+}$ playing the role of the scalar coupling constant $K$, as
$\vec{q}=\mathbf{K}\vec{p}=\lambda_{+}\vec{p}$. This, however, is not so. To
see why, let us consider the trivial solution $p=0$ of Eq.(27). Linear
stability analysis leads to the equation

$\delta\dot{p}=-\delta
p\left(\Delta-\frac{K}{2}\cos\alpha+\frac{J}{2}\cos(2\psi+\beta)\right)$ (29)

which becomes unstable if

$K\cos\alpha-J\cos(2\psi+\beta)\geq 2\Delta$ (30)

Using $K\sin\alpha=J\sin(2\psi+\beta)$ (as $\vec{p}$ is eigenvector of
$\mathbf{K}$) we can rewrite this condition as

$K\cos\alpha+\sqrt{J^{2}-K^{2}\sin^{2}\alpha}\geq 2\Delta$ (31)

which corresponds to set $\lambda_{+}\geq\Delta$.

The solution that branches off from $p=0$ at $\lambda_{+}=\Delta$ is

$p=\sqrt{1-\frac{2\Delta}{K\cos\alpha-J\cos(2\psi-\beta)}}$ (32)

but this is only stationary if $\psi$ is constant ($\omega_{0}=0$) or $J=0$.
For $\omega_{0}\neq 0$, there is no nontrivial equilibrium solution. The trick
of changing to a new frame of reference rotating with $\omega_{0}$ does not
work, as $\mathbf{K}$ does not commute with rotations if $J\neq 0$. Such
operation would change $\mathbf{K}$ to
$\mathbf{\tilde{K}}=\mathbf{R}^{T}\mathbf{K}\mathbf{R}$ which would be itself
time dependent. For $\omega_{0}\neq 0$ the order parameter branches off to an
active state where it rotates while its module oscillates, as if an external
force were acting on the system (see next subsection). We recall the
interpretation of Eq.(24) where $J$ indeed acts as an external drive to the
system. For $\omega_{0}=0$ the order parameter remains stationary and pointing
in the direction of the eigenvector of $\mathbf{K}$, breaking the rotational
symmetry present in the original Kuramoto model. Some of these results can be
demonstrated for any symmetric distribution of frequencies $g(\omega)$, not
just the Lorentzian (see Appendix A).

###  III.3 Active states - complex eigenvalues

Figure 1: Phase diagram for $\alpha=0.5$ and $\Delta=1$. The red line
separates synchronized motion (to the right) from disordered motion (shaded
bluled area to the left). Diagonal lines with slope equal to $\sin\alpha$
separate regions with complex eigenvalues (enclosing the line $J=0$) from real
eigenvalues (enclosing $K=0$). In the phase tuned region, the phase of the
order parameter points in the direction of the eigenvector of $\mathbf{K}$
with the largest eigenvalue if $\omega_{0}=0$. In the oscillatory region
$\vec{p}$ rotates and its module oscillates even if $\omega_{0}=0$ (active
states). Exactly at the diagonals the eigenvalues of $\mathbf{K}$ are
degenerated and the system behaves like the original Kuramoto model. On the
K-axis ($J=0$) we recover the Kuramoto-Sakaguchi model. Figure 2: Phase
diagram in the plane $J-\omega_{0}$ for fixed $K$ and $\alpha$. The x-axis is
plotted in units of $K\sin\alpha$ and the y-axis in units of
$K\sin\alpha/\zeta$. For $\omega_{0}=0$ and $|J|<K\sin\alpha$ (red line) the
system is in the oscillatory region. As $\omega_{0}$ is increased and go
inside the gray cone generated by the black lines (Eq.(35)) $p$ stops
oscillating (Fig. 4(a)). For $|J|>K\sin\alpha$ and $\omega_{0}=0$ (blue lines)
the system is phase tuned region with $\vec{p}$ pointing in the direction of
the eigenvector of $\mathbf{K}$ with largest eigenvalue. As $\omega_{0}$
increases or decreases and go outside the gray cone, $p$ starts oscillating
(Fig. 4(b)).

If $|J|<|K|\sin\alpha$ the eigenvalues of $\mathbf{K}$ are complex and the
relation $K\sin\alpha-J\sin(2\gamma+\beta)=0$ does not hold, as the
eigenvectors are also complex. Still, the trivial solution $p=0$ satisfies
Eq.(27). For $\psi$ we find

$\dot{\psi}=\omega_{0}-\frac{1}{2}\left[K\sin\alpha-J\sin(2\psi+\beta)\right]$
(33)

whose solution is

$(2\omega_{0}-K\sin\alpha)\tan(\psi+\beta/2)=-J+2\Omega\tan[\Omega(t-t_{0})]$
(34)

where $\Omega=\sqrt{(2\omega_{0}-K\sin\alpha)^{2}-J^{2}}/2$. For
$\omega_{0}=0$ this corresponds to the imaginary part of $\lambda$,
generalizing the Kuramoto-Sakaguchi dynamics. For $\omega_{0}\neq 0$,
oscillations in $\psi$ (and therefore in the modulus of $\vec{p}$) require
$\Omega$ to be real, otherwise $\psi$ goes to a constant that can be computed
from Eq.(34) replacing $\Omega\tan[\Omega(t-t_{0})]$ by
$|\Omega|\tanh[|\Omega|(t-t_{0})]$. Therefore, oscillatory solutions can be
‘stabilized’ by choosing $\omega_{0}$ in the window

$\frac{K\sin\alpha-|J|}{\zeta}<\omega_{0}<\frac{K\sin\alpha+|J|}{\zeta}.$ (35)

where $\zeta=2$. The boundaries of this region correspond to phase transitions
from oscillatory to constant behavior of $\vec{p}$. Notice that these
equations for $\psi$ were derived at the trivial solution $p=0$. However, we
can assume that they are also valid for $p\neq 0$, as this would only change
factor $\zeta=2$ in Eq.(33) to $2/(p+1)$, which can be approximated by its
time-averaged value $\zeta=2/(\langle p\rangle+1)$. A similar effect occurs in
the phase-tuned region, where $\omega_{0}$ can induce oscillations in $p$.

The solution of Eq.(29) for $J\neq 0$ is

$\delta p(t)=\delta
p(0)\exp{\left\\{\left(\frac{K}{2}\cos\alpha-\Delta\right)t+\frac{J}{2}\int_{0}^{t}\cos(2\psi(t^{\prime})+\beta)dt^{\prime}\right\\}}.$
(36)

Because $\cos(2\psi(t)+\beta)$ oscillates around zero, the integral in the
exponent remains finite and the condition for $p=0$ to be unstable is

$K\cos\alpha\geq 2\Delta$ (37)

which is independent of $J$. Figure 1 shows the complete phase diagram in the
$J-K$ plane for $\omega_{0}=0$ displaying the regions where the three
different types of solutions occur. A diagram illustrating the transitions as
a function of $\omega_{0}$ is shown in Fig.2.

Figure 3: (a) Phase diagram for $\alpha=0.5$ and Gaussian distribution of
frequencies with $\omega_{0}=0$. Stars show the values of $K$ and $J$
displayed on the plots in panel (b): (b1) $K=2.5;J=0$ (Kuramoto-Sakaguichi);
(b2) $K=2.5;J=0.5$ (active oscillatory state); (b3) $K=2.5;J=1.6$ (phase tuned
state); (b4) $K=-1;J=3$ (phase tuned state with negative $K$. The red lines in
panels (b3) and (b4) show $\cos\psi$ and correspond to the direction of the
eigenvector of $\mathbf{K}$ with largest eigenvalue.

##  IV Simulations

The theoretical diagrams described in Figs. 1 and 2 can be confirmed by
numerical simulations. In this section we show results for $N=10000$
oscillators, $\alpha=0.5$ and Gaussian distribution of natural frequencies, as
the system converges faster in this case than with Lorentzian distributions.
Fig. 3(a) shows the phase diagram in the $J-K$ plane adapted to Gaussian
distributions. Red stars show parameter values selected for the simulations
shown in Fig. 3(b). In the active region $p$ oscillates (panel b2) unless
$J=0$ (panel b1). In the phase tuned region (panels b3 and b4) $p$ converges
to a constant value, and so does $\psi$ (red lines), as $\vec{p}$ becomes
parallel to the eigenvector of $\mathbf{K}$ with largest eigenvalue.

Fig. 4 shows the transitions induced by $\omega_{0}$. Panel (a) displays
results for the active region with $K=2.5$ and $J=0.5$. For $\omega_{0}=0.3$
the oscillations in $p$ persist. For $0.6<\omega_{0}<1.4$ the oscillations are
dumped and $p$ converges to a constant. However, for $\omega_{0}>1.4$ the
system oscillates again (see Fig. 2). Panel (b) shows similar results for the
phase tuned region with $K=2.5$ and $J=1.6$. In this case $p$ remains constant
for $-0.35<\omega_{0}<2.4$. Outside this range $p$ acquires oscillations
induced by $\omega_{0}$.

Figure 4: (a) Transition from oscillatory to static and back to oscillatory
behavior as $\omega_{0}$ is varied in the active oscillatory region for
$K=2.5$ and $J=0.5$. Values of $\omega_{0}$ are indicated in the panels. (b)
Similar transitions for the phase tuned region for $K=2.5$ and $J=1.6$ .

##  V Conclusions

We have extended the Kuramoto model by writing the dynamic equations in vector
form Chandra et al. [2019a] and promoting the scalar coupling constant $K$ to
a $2\times 2$ matrix $\mathbf{K}$. Splitting the matrix into a rotation
$\mathbf{K}_{R}$ plus a symmetric matrix $\mathbf{K}_{S}$ we were able to
write down the equations governing the evolution of the order parameter in the
limit of infinite oscillators and Lorentzian distributions of natural
frequencies. We showed that $\mathbf{K}_{R}$ gives rise to the Kuramoto-
Sakaguchi model whereas the symmetric part $\mathbf{K}_{S}$ acts as internal
driving force that breaks the rotational symmetry.

We constructed the complete phase diagram for the case where the average
frequency of the oscillators, $\omega_{0}$, is zero. Solutions are divided
into disordered (no synchrony), phase tuned (where module and phase of the
order parameter converges to specific values) – and active states (where
module and direction of the order parameter oscillate in time). The original
Kuramoto model and the frustrated model of Kuramoto and Sakaguchi are special
cases of our model. Finally, we showed that non-zero values of $\omega_{0}$
induce further transitions in the behavior of the order parameter, creating a
cone in the $J$-$\omega_{0}$ plane (Fig. 2) inside of which the solutions are
static (constant $p$) and outside are oscillatory.

This novel generalization of the Kuramoto model opens a broad range of
modifications to be explored and studied. Examples include oscillators of
higher dimensions, networks replacing the all-to-all connections of the model
and applications to swarmalators O’Keeffe et al. [2017]. We hope the
contributions discussed in this paper can lead to new interpretations and to
deeper comprehension of the Kuramoto model and its applications.

###### Acknowledgements.

This work was partly supported by FAPESP, grants 2019/24068-0 (AEDB),
2021/10709-3 (GLB), 2016/01343‐7 (MAMA, ICTP‐SAIFR) and CNPq, grant
301082/2019‐7 (MAMA).

##  Appendix A General distribution of natural frequencies

In this appendix we show that the critical line separating asynchronous and
synchronous states can be obtained for general distributions of natural
frequencies $g(\omega)$ using more qualitative arguments. For this we go back
to Eq.(18) and take the scalar product with $\hat{\rho}$ to get

$\displaystyle\dot{\rho}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1-\rho^{2})({\mathbf{K}}\vec{p}\cdot\hat{\rho}).$
(38)

Equilibrium requires that either $\rho=1$, and we only need to find the
direction $\hat{\rho}$, or $\vec{\rho}$ is perpendicular to
$\mathbf{K}\vec{p}$. For $\rho=1$ we obtain

$0=\vec{\omega}\times\hat{\rho}+{\mathbf{K}}\vec{p}-({\mathbf{K}}\vec{p}\cdot\hat{\rho})\hat{\rho}.$
(39)

Taking the cross product with $\hat{\rho}$ we get
$\vec{\omega}=-\hat{\rho}\times{\mathbf{K}}\vec{p}$ and
$\omega=q\sin(\phi-\xi)$, which requires $|\omega|<q$.

The tree important vectors are: the order parameter
$\vec{p}=(p\cos\psi,p\sin\psi)$, the ansatz vector
$\vec{\rho}=(\rho\cos\phi,\rho\sin\phi)$ and the auxiliary vector
$\vec{q}={\mathbf{K}}\vec{p}=(q\cos\xi,q\sin\xi)$. Using
$d\omega=q\cos(\phi-\xi)d\phi$ we obtain, for the x-component of Eq.(17)

$p\cos\psi=\int_{-q}^{+q}\cos\phi\,g(\omega)d\omega=\int_{\xi-\pi/2}^{\xi+\pi/2}q\cos\phi\cos(\phi-\xi)g(\omega(\phi))d\phi$
(40)

as the solutions for $|\omega|>q$ (corresponding to
${\mathbf{K}}\vec{p}\cdot\hat{\rho}=0$) do not contribute to $\vec{p}$. The
trivial solution corresponds to $p=0$, as this implies $q=0$. The critical
condition for synchronization corresponds to the non-trivial solution
calculated at $p=0$. Evaluating $g(\omega)$ at $q=0$ and doing the integral
over $\phi$ we get

$p\cos\psi=(g(0)\pi/2)q\cos\xi$ (41)

and, similarly,

$p\sin\psi=(g(0)\pi/2)q\sin\xi.$ (42)

Dividing one equation by the other we find that $\tan\xi=\tan\psi$, implying
that $\vec{q}$ is parallel to $\vec{p}$, or
${\mathbf{K}}\vec{p}=\lambda\vec{p}$. The direction of the order parameter is
therefore fixed by coupling matrix if its eigenvalues are real. Moreover,
replacing $q$ by $\lambda p$ in Eqs. (41) and (42) we find the critical
condition for synchronization as

$\lambda_{c}=2/(\pi g(0))$ (43)

which is identical to the condition found for the Lorenztian distribution,
where $g(0)=(\Delta\pi)^{-1}$.

## References

  * Kuramoto [1975] Yoshiki Kuramoto.  Self-entrainment of a population of coupled non-linear oscillators.  In _International Symposium on Mathematical Problems in Theoretical Physics_ , pages 420–422. Springer-Verlag, Berlin/Heidelberg, 1975\.  doi: 10.1007/BFb0013365.  URL http://www.springerlink.com/index/10.1007/BFb0013365.
  * Kuramoto [1984] Yoshiki Kuramoto.  Chemical Waves.  In _Chemical Oscillations, Waves, and Turbulence_ , pages 89–110. Springer Berlin Heidelberg, 1984.  doi: 10.1007/978-3-642-69689-3˙6.  URL http://www.springerlink.com/index/10.1007/978-3-642-69689-3{_}6.
  * Hong and Strogatz [2011] Hyunsuk Hong and Steven H Strogatz.  Kuramoto model of coupled oscillators with positive and negative coupling parameters: an example of conformist and contrarian oscillators.  _Physical Review Letters_ , 106(5):054102, 2011\.
  * Yeung and Strogatz [1999] MK Stephen Yeung and Steven H Strogatz.  Time delay in the kuramoto model of coupled oscillators.  _Physical Review Letters_ , 82(3):648, 1999.
  * Breakspear et al. [2010] Michael Breakspear, Stewart Heitmann, and Andreas Daffertshofer.  Generative models of cortical oscillations: neurobiological implications of the kuramoto model.  _Frontiers in human neuroscience_ , 4:190, 2010.
  * Rodrigues et al. [2016] Francisco A. Rodrigues, Thomas K D M Peron, Peng Ji, and J??rgen Kurths.  The Kuramoto model in complex networks.  _Physics Reports_ , 610:1–98, 2016.  ISSN 03701573.  doi: 10.1016/j.physrep.2015.10.008.  URL http://dx.doi.org/10.1016/j.physrep.2015.10.008.
  * Climaco and Saa [2019] Joyce S. Climaco and Alberto Saa.  Optimal global synchronization of partially forced kuramoto oscillators.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 29(7):073115, 2019.  doi: 10.1063/1.5097847.  URL https://doi.org/10.1063/1.5097847.
  * Gomez-Gardenes et al. [2011] Jesus Gomez-Gardenes, Sergio Gomez, Alex Arenas, and Yamir Moreno.  Explosive synchronization transitions in scale-free networks.  _Physical Review Letters_ , 106(12):1–4, 2011\.  ISSN 00319007.  doi: 10.1103/PhysRevLett.106.128701.
  * Ji et al. [2013] Peng Ji, Thomas K Dm Peron, Peter J. Menck, Francisco A. Rodrigues, and J??rgen Kurths.  Cluster explosive synchronization in complex networks.  _Physical Review Letters_ , 110(21):1–5, 2013\.  ISSN 00319007.  doi: 10.1103/PhysRevLett.110.218701.
  * Acebrón et al. [2005] Juan A. Acebrón, L. L. Bonilla, Conrad J Pérez Vicente, Félix Ritort, and Renato Spigler.  The Kuramoto model: A simple paradigm for synchronization phenomena.  _Reviews of Modern Physics_ , 77(1):137–185, 2005\.  ISSN 00346861.  doi: 10.1103/RevModPhys.77.137.
  * Dörfler and Bullo [2011] Florian Dörfler and Francesco Bullo.  On the critical coupling for kuramoto oscillators.  _SIAM Journal on Applied Dynamical Systems_ , 10(3):1070–1099, 2011.
  * Olmi et al. [2014] Simona Olmi, Adrian Navas, Stefano Boccaletti, and Alessandro Torcini.  Hysteretic transitions in the kuramoto model with inertia.  _Physical Review E_ , 90(4):042905, 2014.
  * Childs and Strogatz [2008] Lauren M. Childs and Steven H. Strogatz.  Stability diagram for the forced Kuramoto model.  _Chaos_ , 18(4):1–9, 2008.  ISSN 10541500.  doi: 10.1063/1.3049136.
  * Moreira and de Aguiar [2019a] Carolina A Moreira and Marcus AM de Aguiar.  Global synchronization of partially forced kuramoto oscillators on networks.  _Physica A: Statistical Mechanics and its Applications_ , 514:487–496, 2019a.
  * Moreira and de Aguiar [2019b] Carolina A Moreira and Marcus AM de Aguiar.  Modular structure in c. elegans neural network and its response to external localized stimuli.  _Physica A: Statistical Mechanics and its Applications_ , 533:122051, 2019b.
  * O’Keeffe et al. [2017] Kevin P O’Keeffe, Hyunsuk Hong, and Steven H Strogatz.  Oscillators that sync and swarm.  _Nature communications_ , 8(1):1–13, 2017.
  * O’Keeffe et al. [2022] Kevin O’Keeffe, Steven Ceron, and Kirstin Petersen.  Collective behavior of swarmalators on a ring.  _Physical Review E_ , 105(1):014211, 2022.
  * Chandra et al. [2019a] Sarthak Chandra, Michelle Girvan, and Edward Ott.  Continuous versus discontinuous transitions in the d-dimensional generalized kuramoto model: Odd d is different.  _Physical Review X_ , 9(1):011002, 2019a.
  * Ott and Antonsen [2008] Edward Ott and Thomas M. Antonsen.  Low dimensional behavior of large systems of globally coupled oscillators.  _Chaos_ , 18(3):1–6, 2008.  ISSN 10541500.  doi: 10.1063/1.2930766.
  * Chandra et al. [2019b] Sarthak Chandra, Michelle Girvan, and Edward Ott.  Complexity reduction ansatz for systems of interacting orientable agents: Beyond the kuramoto model.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 29(5):053107, 2019b.
  * Barioni and de Aguiar [2021a] Ana Elisa D Barioni and Marcus AM de Aguiar.  Complexity reduction in the 3d kuramoto model.  _Chaos, Solitons & Fractals_, 149:111090, 2021a.
  * Barioni and de Aguiar [2021b] Ana Elisa D Barioni and Marcus AM de Aguiar.  Ott–antonsen ansatz for the d-dimensional kuramoto model: A constructive approach.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 31(11):113141, 2021b.
  * Sakaguchi and Kuramoto [1986] Hidetsugu Sakaguchi and Yoshiki Kuramoto.  A soluble active rotater model showing phase transitions via mutual entertainment.  _Progress of Theoretical Physics_ , 76(3):576–581, 1986.
  * Yue et al. [2020] Wenqi Yue, Lachlan D Smith, and Georg A Gottwald.  Model reduction for the kuramoto-sakaguchi model: The importance of nonentrained rogue oscillators.  _Physical Review E_ , 101(6):062213, 2020.