# Synchronization of Spin-Torque Oscillators via Continuation Method

Denis Nikitin, Carlos Canudas de Wit, Paolo Frasca, Ursula Ebels

► To cite this version:

Denis Nikitin, Carlos Canudas de Wit, Paolo Frasca, Ursula Ebels. Synchronization of Spin-Torque Oscillators via Continuation Method. IEEE Transactions on Automatic Control, 2023, 68 (11), pp.6621-6635. 10.1109/TAC.2023.3298288. [hal-03315718](https://hal.science/hal-03315718)

HAL Id: hal-03315718

https://hal.science/hal-03315718

Submitted on 5 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Synchronization of Spin-Torque Oscillators via Continuation Method

Denis Nikitin, Carlos Canudas-de-Wit, Paolo Frasca and Ursula Ebels

***Abstract***—In this paper we study synchronization phenomenon for spin-torque oscillators coupled on a ring. Spin-torque oscillators are nanoelectronic devices which promise efficient microwave generation provided they are synchronized in large arrays. Due to their nonlinear and non-isochronous nature, their synchronization is difficult to analyse explicitly. We employ a recently developed continuation method that transforms the network of coupled oscillators (each described by an ordinary differential equation) into a single nonlinear partial differential equation (PDE). We then analyse the synchronization of this PDE in two cases: when all the oscillators are identical and when there are two different types of oscillators. In the case of identical oscillators we reconstruct all possible synchronous solutions in the system and provide explicit conditions for stability. For non-identical oscillators we derive and solve a differential synchronization condition which allows us to reconstruct the shape of the equilibrium profiles. All the presented results, which are derived for the PDE, are validated by numerical simulations of the original network of ODEs.

***Index Terms***—Synchronization of Large-Scale Networks, Spin-Torque Oscillators, Non-Isochronous Oscillators, Kuramoto

## I. INTRODUCTION

SYNCHRONIZATION is an astonishing phenomenon that can occur in large systems with many interacting nonlinear agents. Manifestations of this effect are usually associated with a significant increase in the energy efficiency of the system because many agents behave as a whole, oscillating at a single common frequency. In the control community, the Kuramoto model was mainly investigated as a prototypical example of synchronization. It is used to express a behaviour of a network of coupled oscillators whose dynamics are described using phase angles. The Kuramoto model and its second-order generalization can be applied to study various applications such as chemical oscillators [1], synchronization of smart grids [2], power networks [3] and even crowd synchrony on London's Millennium Bridge [4]. The issue of synchronization in the Kuramoto model depending on different network topologies and different parameters is extensively studied, see [5] for the review and [2], [6] for recent comprehensive synchronization conditions.

The Kuramoto model describes only the phase dynamics of oscillators, assuming that the oscillation amplitude is constant for each oscillator. In many practical applications, however, the relationship between amplitude and phase dynamics cannot be neglected. The class of oscillators in which the oscillation frequency can vary with amplitude is called non-isochronous.

D. Nikitin, C. Canudas-de-Wit and P. Frasca are with Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, GIPSA-lab, 38000 Grenoble, France. U. Ebels is with Univ. Grenoble Alpes, CEA, CNRS, Grenoble INP, Spintec, 38000 Grenoble, France.

Many well-known models, such as van der Pol oscillator [7] or FitzHugh-Nagumo neuron model [8] have this property.

This class of non-isochronous oscillators includes spin-torque oscillators (STO), nanoelectronic devices that are based on the spin-transfer torque effect discovered by [9] and [10]. It appears that an electric direct current which passes through a magnetized layer can become spin-polarized, and moreover this spin-polarized current can further transfer angular momentum to another magnetized layer. This transfer induces torque on the magnetization of the second layer, which can lead to switching of the magnetization direction. A steady state magnetization precession can be induced under appropriate conditions of external field and for current densities larger than a critical value. Due to the fast precession of magnetization in a ferromagnetic layer, STOs produce microwave output voltage signals. Thus large arrays of STOs can theoretically serve as very efficient microwave generators. This is why the question of synchronization of STOs is very important: synchronous oscillations of many oscillators amplify each other due to constructive interference, while asynchronous oscillations exhibit destructive interference and thus produce less power. Different magnetic configurations facilitating synchronization were studied in [11], [12], while the impact of different device properties were investigated in [13], [14]. However, analytic studies were mostly limited to investigating properties of amplitude-frequency coupling, synchronization to an external oscillating force and phase-locking effects, as in [15], [16].

In this paper we study the synchronization of a large array of STOs with a ring interaction topology. In order to completely describe synchronization on a ring, we are going to reconstruct different synchronization profiles possible in the system and to derive particular conditions on the system parameters which can be checked to ensure that these synchronized solutions are stable. One could then use these conditions as a guidance for realization and deployment of large arrays of synchronously operating STOs. We approach this problem by utilizing a recently developed continuation method, which transforms a network of coupled ODEs into a single PDE. The general continuation method was introduced by the authors in [17] and was applied in [18] to derive an urban traffic PDE model and in [19] to stabilize a Stuart-Landau laser chain. To demonstrate applicability of the method to the problem of network synchronization, in Section II we start with a toy example of analysis of the Kuramoto model synchronization on a ring.

We then proceed to introduce the STO model and its continuation in Sections III and IV, respectively. Next, the analysis of STO synchronization is split into two cases. In the first case we assume that all oscillators are the same, i.e. the
system parameters do not change along the ring. In the second case, we extend our analysis to the case when there are several different types of oscillators in the system. This setup allows us to investigate how the synchronization properties change due to possible parameter variations that may arise because of fabrication inaccuracies among other reasons.

In the case of identical oscillators in Section V, our main results provide a characterization of all synchronized solutions and the explicit conditions on the system parameters which can be checked to ensure that a chosen synchronized solution is stable. Due to the complexity of STOs, these conditions depend on all of the system's parameters, but we found that the coupling phase between oscillators is one of the most important ones. We validated our results using numerical simulations and found that the derived conditions are able to reconstruct the behavior of the original system.

The analysis of the system with non-identical oscillators is much more challenging. In Section VI we assume that the system consists of two different types of oscillators (e.g. they were produced in two batches). We then derive a differential equation serving as a condition for synchronization and then solve it to obtain an implicit function describing synchronous solution's profile. Thus we are successfully able to analytically reconstruct solutions arising in the inhomogeneous system, which we demonstrate using numerical simulations.

## II. TOY EXAMPLE: RING OF KURAMOTO OSCILLATORS

To illustrate how a PDE approximation of an oscillator network derived via the continuation method can be used, in this section we focus our attention on a network of Kuramoto oscillators with local interactions, namely coupled on a 1D ring. Deriving PDE representation of Kuramoto system, we show that this model can be more appropriate for analysis (in the same way as continuous dynamical systems can be more tractable than the discrete ones). As a toy example, here we present one possible application of this representation, namely we analytically find a synchronization threshold for a 1D ring topology. Problem of computation of a general synchronization threshold for different topologies and frequency distributions was recently solved by [2] and [6] with the help of graph theory. However, the method presented in these papers is sophisticated and not straightforward to extend to other types of oscillators. Contrary, the idea based on continuation which is presented in this section helps to find a synchronization condition in a very natural way. Moreover, this method will be extended in the next sections to a more general class of non-isochronous oscillators in complex domain. Apart from deriving synchronization threshold we demonstrate that the continuation method produces an accurate representation of the Kuramoto network by performing numerical simulations comparing the original ODE system with the obtained PDE.

### A. Continuum description of the Kuramoto model

We start by analysing a system of Kuramoto oscillators

$$ \dot{\phi}_i = \omega_i + F(\sin(\phi_{i+1} - \phi_i) - \sin(\phi_i - \phi_{i-1})), \quad (1) $$

where $\phi_i$ is a phase angle of i-th oscillator, $\omega_i$ is its natural frequency and $F > 0$ is a coupling strength. Each oscillator is coupled with its two neighbours, forming a closed ring. We assume that there are $n$ oscillators and that each oscillator has a position on a ring defined by $x_i \in [0, 2\pi)$, with $x_{i+1} - x_i = \Delta x$ and $x_1 - x_n + 2\pi = \Delta x$, meaning that the oscillators are spaced equally on the ring. Using these positions, we can further define a natural frequency function $\omega(x_i) = \omega_i$ and then a state function $\phi(x_i) = \phi_i$.

The main idea now is to use a PDE approximation of the original ODE system (1). A general description of our method to derive PDE approximations for nonlinear systems of ODEs is given in [17]. Here we propose a self-contained example application of this method to obtain a Kuramoto PDE. Since we assume continuity of an underlying space on a ring, let us define a function $s(x)$ by the rule

$$ s(x_{i-1/2}) = \sin(\phi(x_i) - \phi(x_{i-1})), \quad (2) $$

which leads to the system

$$ \dot{\phi}(x_i) = \omega(x_i) + F(s(x_{i+1/2}) - s(x_{i-1/2})). $$

Now we can calculate the Taylor expansion of the function $\phi(x_i)$ at the point $x_{i-1/2}$:

$$ \phi(x_i) = \phi(x_{i-1/2}) + \frac{\Delta x}{2} \frac{\partial \phi}{\partial x}(x_{i-1/2}) + \frac{\Delta x^2}{8} \frac{\partial^2 \phi}{\partial x^2}(x_{i-1/2}) + \dots $$

The same expansion can be used for $\phi(x_{i-1})$. Thus the argument of sine in (2) is

$$ \phi(x_i) - \phi(x_{i-1}) = \Delta x \frac{\partial \phi}{\partial x}(x_{i-1/2}) + O(\Delta x^3). $$

In general these series are infinite, but we want to obtain a low-order approximation, therefore we assume that

$$ \phi(x_i) - \phi(x_{i-1}) = \Delta x \frac{\partial \phi}{\partial x}(x_{i-1/2}). $$

Note that this coincides with the usual second-order finite difference method for the discretization of first-order spatial derivatives in PDEs. As a result, we finally obtain that

$$ s(x) = \sin \left( \Delta x \frac{\partial \phi}{\partial x}(x) \right). \quad (3) $$

Moreover, using the same procedure applied to the function $s(x)$, we can write

$$ \frac{\partial \phi}{\partial t} = \omega(x) + F \Delta x \frac{\partial s}{\partial x}, $$

and using the approximation of $s(x)$ in (3), the obtained PDE approximation is

$$ \frac{\partial \phi}{\partial t} = \omega(x) + F \Delta x \frac{\partial}{\partial x} \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right). \quad (4) $$

We validate this PDE approximation by the simulation of an ODE system with $n = 50$ oscillators, placed on a ring. Having freedom to choose any nontrivial natural frequency function, we take $\omega(x) = 1 + x \sin(2x)$ for $x \in [0, 2\pi)$ (in general any function can be used, but for the future analysis we choose an integrable one). The coupling strength is set to $F = 4$. We numerically simulate the approximated PDE (4) on a grid with
<figure><img src="image_2.png" /><figcaption>Figure 1. Comparison between simulations of a Kuramoto ODE network (1) with $n = 50$ oscillators and its PDE approximation (4). **Left:** Snapshot of profiles of both systems at time $T = 500$. **Right:** Evolution of a mean-square absolute divergence between solutions.</figcaption></figure>

500 points. The results of simulation are shown in Fig. 1. It is clear that while the solution to the ODE system splits into several clusters, solution to the PDE model continuously connects these clusters (Fig. 1, left), remaining rather accurate at the positions of the oscillators of the original ODE system (Fig. 1, right).

## B. Synchronization threshold for Kuramoto oscillators

The main advantage of describing the system in terms of partial derivatives is that now the space becomes a continuum, thus integrals can be taken (and in general integrals are much more tractable than series). We will show how the obtained PDE (4) can be used to find the parameter $F^*$ for which phase transition from the complete synchronization to the emergence of clusters occurs. Namely, let us try to find an equilibrium solution $\phi^*$ of (4) in case of complete synchronization. It is clear that then there exists $\bar{\omega} = \frac{1}{2\pi} \int_0^{2\pi} \omega(x)dx$ such that all oscillators share the same frequency:

$$ \frac{\partial \phi^*}{\partial t} = \bar{\omega}. $$

Therefore the equilibrium solution should satisfy

$$ F \Delta x \frac{\partial}{\partial x} \sin \left( \Delta x \frac{\partial \phi^*}{\partial x} \right) = \bar{\omega} - \omega(x). \quad (5) $$

Let us integrate this equation from $x_0$ to $x_1$, where both are chosen arbitrary:

$$ \begin{aligned} \sin \left( \Delta x \frac{\partial \phi^*}{\partial x} (x_1) \right) - \sin \left( \Delta x \frac{\partial \phi^*}{\partial x} (x_0) \right) &= (6) \\ &= \frac{1}{F \Delta x} (\Omega(x_1) - \Omega(x_0)), \end{aligned} $$

where $\Omega(x)$ is some primitive function of $\bar{\omega} - \omega(x)$. Rearranging terms, we obtain

$$ \begin{aligned} \sin \left( \Delta x \frac{\partial \phi^*}{\partial x} (x_1) \right) - \frac{1}{F \Delta x} \Omega(x_1) &= (7) \\ &= \sin \left( \Delta x \frac{\partial \phi^*}{\partial x} (x_0) \right) - \frac{1}{F \Delta x} \Omega(x_0) =: C, \end{aligned} $$

where $C$ is therefore some constant independent of the choice of $x_0$ and $x_1$. We obtained that the existence of an equilibrium solution is equivalent to the existence of the primitive function $\Omega(x)$ written in the form

$$ \frac{1}{F \Delta x} \Omega(x) = \sin \left( \Delta x \frac{\partial \phi^*}{\partial x} (x) \right) + C. \quad (8) $$

If such $\Omega(x)$ exists, $\phi^*$ can be recovered by taking arcsine and then integrating. Therefore, a complete synchronization for a given $F$ is possible if and only if there exists $\Omega(x)$ such that (8) is possible, in a sense that the sine value lies in the interval $[-1, 1]$. Essentially this means that

$$ \Omega(x) \in [-F\Delta x + C, F\Delta x + C] \quad \forall x \in [0, 2\pi], $$

with $C$ being some constant. Recalling that $\Omega(x)$ is a primitive function of $\bar{\omega} - \omega(x)$ and that in general it is defined up to a constant, this is equivalent to the condition

$$ \max_{x \in [0, 2\pi]} \Omega(x) - \min_{x \in [0, 2\pi]} \Omega(x) \le 2F\Delta x $$

for any $\Omega(x)$. To recover synchronization threshold $F^*$ it requires only to replace inequality with equality sign:

$$ F^* = \frac{1}{2\Delta x} \left( \max_{x \in [0, 2\pi]} \Omega(x) - \min_{x \in [0, 2\pi]} \Omega(x) \right). \quad (9) $$

Synchronization threshold (9) provides a condition on the existence of equilibrium solutions. It is possible to prove that for all $F > F^*$ there will exist a stable equilibrium solution.

**Theorem 1.** For all $F > F^*$, there exists an equilibrium solution $\phi^*$, satisfying (8), which is locally asymptotically stable.

*Proof.* Without loss of generality we assume that $C = 0$ (because the primitive function $\Omega(x)$ is defined up to a constant). Note that $\Omega(x) \in [-F^*\Delta x, F^*\Delta x]$. Then the equilibrium solution can be recovered from (8) (up to a constant) as

$$ \phi^*(x) = \frac{1}{\Delta x} \int_{0}^{x} \arcsin \left( \frac{1}{F \Delta x} \Omega(x) \right) dx. \quad (10) $$

Now let us assume that the equilibrium solution is slightly perturbed: $\phi = \phi^* + \tilde{\phi}$. Then, by (4),

$$ \frac{\partial \tilde{\phi}}{\partial t} = \omega(x) - \bar{\omega} + F \Delta x \frac{\partial}{\partial x} \sin \left( \Delta x \frac{\partial \phi^*}{\partial x} + \Delta x \frac{\partial \tilde{\phi}}{\partial x} \right). \quad (11) $$

Rewriting sine, we get

$$ \sin \left( \Delta x \frac{\partial \phi^*}{\partial x} + \Delta x \frac{\partial \tilde{\phi}}{\partial x} \right) \approx \sin \left( \Delta x \frac{\partial \phi^*}{\partial x} \right) + \cos \left( \Delta x \frac{\partial \phi^*}{\partial x} \right) \Delta x \frac{\partial \tilde{\phi}}{\partial x}, $$

where the fact that $\tilde{\phi}$ is a small perturbation was used. Now (5) cancels the natural frequencies, therefore we arrive at

$$ \frac{\partial \tilde{\phi}}{\partial t} = F \Delta x^2 \frac{\partial}{\partial x} \left[ \cos \left( \Delta x \frac{\partial \phi^*}{\partial x} \right) \frac{\partial \tilde{\phi}}{\partial x} \right]. \quad (12) $$

This is a standard linear diffusion equation with the diffusion coefficient $\cos(\Delta x \frac{\partial \phi^*}{\partial x})$. For stability it remains to prove that this coefficient is always positive. Indeed,

$$ \begin{aligned} \cos \left( \Delta x \frac{\partial \phi^*}{\partial x} \right) &= \cos \left( \arcsin \left( \frac{1}{F \Delta x} \Omega(x) \right) \right) = (13) \\ &= \sqrt{1 - \left( \frac{\Omega(x)}{F \Delta x} \right)^2} > \sqrt{1 - \left( \frac{\Omega(x)}{F^* \Delta x} \right)^2} \ge 0, \end{aligned} $$

and the linearised system is locally asymptotically stable. $\square$

To validate this analysis we use the same parameters of the simulation as for Fig. 1: the length of the ring is $L = 2\pi$,
<figure><img src="image_3.png" /><figcaption>Figure 2. Desynchronization frequency $\tilde{\omega}$ depending on the coupling strength $F$. Other parameters as in Fig. 1. Left: $F \in [0, 20]$. Right: zoom in, $F \in [18, 20]$.</figcaption></figure>

the number of ODE nodes $n = 50$, and the natural frequency $\omega(x) = 1 + x \sin(2x)$ (which is an integrable function, thus can be analytically treated). By definition of the positions of nodes, $\Delta x = L/n = 2\pi/50$. Further,

$$\int_{0}^{x} \omega(s) ds = x - \frac{1}{2} x \cos(2x) + \frac{1}{4} \sin(2x),$$

thus $\tilde{\omega} = 1/2$. Primitive function $\Omega(x)$ can be taken as

$$\Omega(x) = \int_{0}^{x} (\tilde{\omega} - \omega(s)) ds = \frac{1}{2} x \cos(2x) - \frac{1}{4} \sin(2x) - \frac{1}{2} x,$$

with $\max \Omega(x) = \Omega(3.06) = 0.0203$ and $\min \Omega(x) = \Omega(4.765) = -4.726$. Substituting these values in (9) gives

$$F^* \approx 18.88. \qquad (14)$$

The value (14) is the smallest $F$ for which the equilibrium solution exists. To verify the result (14) for the original system (1), we simulated it for $F \in [0, 20]$ and calculated $\tilde{\omega} = \max \dot{\phi}_i - \min \dot{\phi}_i$, which we call desynchronization frequency. In case of complete synchronization $\tilde{\omega}$ should be zero. Indeed, Fig. 2 shows that $\tilde{\omega}$ is zero for $F > 18.9$, and it increases when $F$ becomes smaller, thus behaving in accordance with the derived value of $F^*$ in (14).

### III. MODEL FOR THE RING OF SPIN-TORQUE OSCILLATORS

It was shown in the previous section that synchronization analysis can be easily performed for Kuramoto oscillators using PDE representation. The same synchronization conditions were already derived for ODE representation by [2] and [6], however these results were based on an extensive use of nontrivial graph theory and linear algebra. We can further show that PDE-based models allow for a more natural analysis of systems by applying the continuation method to a more complex class of oscillators, namely non-isochronous oscillators. As a particular example of this class of oscillators we focus on spin-torque oscillators (STO).

A typical spin-torque oscillator consists of two ferromagnetic layers, a thick one called "fixed" and a thin one called "free", see the left panel of Fig. 3. The "fixed" layer will spin polarize an electrical current that passes through it. This polarized current, when traversing the second "free" layer, will transfer a spin angular momentum, that manifests itself as

<figure><img src="image_4.png" /><figcaption>Figure 3. Schematic representation of a possible geometry of spin-torque oscillator. Red blocks represent ferromagnetic layers with their magnetization directions denoted by black arrows. Electrons flow from bottom to top, first passing through the "fixed" magnetic layer which induces spin polarization coinciding with its magnetization direction $\mathcal{P}$. The magnetization $\mathcal{M}$ of the "free" magnetic layer then oscillates under the effect of polarized current and the effective magnetic field $\mathcal{H}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Close view on the dynamics of the magnetization $\mathcal{M}$ of the "free" layer, governed by equation (15). Damping and current-induced spin-transfer torque compensate each other, stabilizing steady oscillations caused by precession around the magnetic field $\mathcal{H}$.</figcaption></figure>

a torque acting on the free layer magnetization. This effect creates precession, depicted in the right panel of Fig. 3. Denote the magnetization of the "free" magnetic layer by vector $\mathcal{M}$, the magnetization of the "fixed" magnetic layer by vector $\mathcal{P}$ and the effective magnetic field by vector $\mathcal{H}$. Then the dynamics of the "free" layer magnetization is governed by the Landau-Lifshitz-Gilbert equation:

$$\frac{\partial \mathcal{M}}{\partial t} = \underbrace{-\gamma(\mathcal{M} \times \mathcal{H})}_{\text{precession}} + \underbrace{\frac{\alpha}{|\mathcal{M}|} \left( \mathcal{M} \times \frac{\partial \mathcal{M}}{\partial t} \right)}_{\text{damping}} + \underbrace{\frac{\sigma I}{|\mathcal{M}|} (\mathcal{M} \times (\mathcal{M} \times \mathcal{P}))}_{\text{spin transfer}} \qquad (15)$$

where parameters $\gamma$, $\alpha$ and $\sigma$ depend on the system's geometry and materials, and $I$ is a current density which is applied to the system. For a review of the spin-transfer torque effect and STOs see [16], [20], [15].

Equation (15) can be simplified for analysis. Magnetization vector $\mathcal{M}$ oscillates around the effective magnetic field vector
Let us project $\mathcal{M}$ on a plane orthogonal to $\mathcal{H}$ and denote the resulting projection via a complex variable $c$. Then, with some additional transformations (see [16] for details) it is possible to show that the magnetization dynamics (15) of an STO can be modelled through:

$$ \dot{c} = i(\omega + Np)c - \Gamma_G(1 + Qp)c + \sigma I(1 - p)c, \quad (16) $$

where $p = |c|^2$ represents a squared amplitude of oscillations, $\omega$ is a linear frequency, $N$ is a nonlinear frequency coefficient, $\Gamma_G$ is a linear damping, $Q$ is a nonlinear damping coefficient, and $I$ and $\sigma$ are the same as in (15). Model (16) is nonlinear since the oscillations' frequency depends on the amplitude through the frequency gain $N$. In case of spin-torque oscillators this amplitude-related frequency shift happens to be very strong, thus these oscillations cannot be described by simpler linear models.

If $\sigma I \le \Gamma_G$, the origin $c=0$ is a stable equilibrium point. Oscillations will occur if $\sigma I > \Gamma_G$. Assuming it is true, define a linear part of sum of dissipative terms $\Gamma = \sigma I - \Gamma_G > 0$ and further a nonlinear gain of sum of dissipative terms $S = \Gamma_G Q + \sigma I$, thus the system (16) can be written as

$$ \dot{c} = i(\omega + Np)c + (\Gamma - Sp)c. \quad (17) $$

System (17) will oscillate with amplitude $\sqrt{p} = \sqrt{\Gamma/S}$ and with frequency $\dot{\phi} = \omega + N\Gamma/S$, where $\phi$ is a phase of an oscillator. For the amplitude of oscillations to be well defined, we also require $S > 0$.

In physical systems a phase of an oscillator usually changes much faster than its amplitude. Therefore the model (17) is often studied in amplitude-phase representation $c = \sqrt{p}e^{i\phi}$, where $\sqrt{p}$ is an amplitude of oscillations and $\phi$ is a phase of an oscillator. Instead of writing two separate equations for them, we will write model (17) in logarithmic representation. Define $z = \ln c$. Then the real part of $z$ will represent the amplitude, namely $\exp\{2\operatorname{Re}z\} = p$. Let us denote $r := \operatorname{Re}z = \frac{1}{2}\ln p$. The imaginary part of $z$ is a phase of an oscillator, $\phi := \operatorname{Im}z$, thus such transformation allows to track phase information directly. Since $dc = c \cdot dz$, the model (17) now becomes

$$ \dot{z} = \Gamma + i\omega - (S - iN)e^{2\operatorname{Re}z}. \quad (18) $$

Now let us move to a system of coupled oscillators. We assume the oscillators are placed on a ring, and each oscillator is coupled with its two neighbours. As in the previous section, let $n$ denote the number of oscillators and let $x_i \in [0, 2\pi)$ be a position on a ring of the $i$-th oscillator. The distance between oscillators is $x_{i+1} - x_i = \Delta x$ and $x_1 - x_n + 2\pi = \Delta x$, meaning that the oscillators are spaced equally on the ring. Coupling between oscillators means that each oscillator has its neighbors' states as an external force:

$$ \dot{c}_i = i(\omega_i + N_i p_i)c_i + (\Gamma_i - S_i p_i)c_i + F_i(c_{i-1} + c_{i+1}). \quad (19) $$

Here $F_i$ is a (possibly complex) coupling constant, with an amplitude representing coupling strength and a phase representing coupling phase.

Using logarithmic representation, the model (19) reads as

$$ \dot{z}_i = \Gamma_i + i\omega_i - (S_i - iN_i)e^{2\operatorname{Re}z_i} + F_i(e^{z_{i-1}-z_i} + e^{z_{i+1}-z_i}). \quad (20) $$

## IV. CONTINUATION AND SYNCHRONIZATION CONDITION

It is now possible to perform continuation for the coupled system in the same way it was done for Kuramoto oscillators in the previous section. According to [17], the continuation of (20) is performed in several steps:

$$ 1) \quad z_{i-1} - z_i \rightarrow -\Delta x \frac{\partial z}{\partial x_{i-1/2}} $$

$$ 2) \quad z_{i+1} - z_i \rightarrow \Delta x \frac{\partial z}{\partial x_{i+1/2}} $$

$$ 3) \quad e^{-\Delta x \frac{\partial z}{\partial x_{i-1/2}}} \rightarrow e^{-\Delta x \frac{\partial z}{\partial x_i}} - \Delta x \frac{1}{2} \frac{\partial}{\partial x} e^{-\Delta x \frac{\partial z}{\partial x_i}} $$

$$ 4) \quad e^{\Delta x \frac{\partial z}{\partial x_{i+1/2}}} \rightarrow e^{\Delta x \frac{\partial z}{\partial x_i}} + \Delta x \frac{1}{2} \frac{\partial}{\partial x} e^{\Delta x \frac{\partial z}{\partial x_i}} $$

Using these continuations, we finally get

$$ e^{z_{i-1}-z_i} + e^{z_{i+1}-z_i} \rightarrow \left( e^{\Delta x \frac{\partial z}{\partial x_i}} + e^{-\Delta x \frac{\partial z}{\partial x_i}} \right) + \Delta x \frac{\partial}{\partial x} \left( \frac{e^{\Delta x \frac{\partial z}{\partial x_i}} - e^{-\Delta x \frac{\partial z}{\partial x_i}}}{2} \right), $$

or simply

$$ e^{z_{i-1}-z_i} + e^{z_{i+1}-z_i} \rightarrow 2 \cosh \left( \Delta x \frac{\partial z}{\partial x} \right) + \Delta x \frac{\partial}{\partial x} \sinh \left( \Delta x \frac{\partial z}{\partial x} \right). $$

Thus, system (20) can be written using PDE model as

$$ \frac{\partial z}{\partial t} = \Gamma + i\omega - (S - iN)e^{2\operatorname{Re}z} + F \left[ 2\cosh \left( \Delta x \frac{\partial z}{\partial x} \right) + \Delta x \frac{\partial}{\partial x} \sinh \left( \Delta x \frac{\partial z}{\partial x} \right) \right], \quad (21) $$

where parameters $\Gamma$, $\omega$, $S$, $N$ and $F$ are (possibly) varying functions of space, determined by approximating sampled values $\Gamma_i$, $\omega_i$, $S_i$, $N_i$ and $F_i$ at points $x_i$.

Separating (21) into a system of two equations for $r = \operatorname{Re}z$ and $\phi = \operatorname{Im}z$, one gets

$$ \left\{ \begin{aligned} \frac{\partial r}{\partial t} &= \Gamma - Se^{2r} + \\ & \quad + \operatorname{Re}F \left[ 2\cosh \left( \Delta x \frac{\partial r}{\partial x} \right) \cos \left( \Delta x \frac{\partial \phi}{\partial x} \right) + \right. \\ & \quad \left. + \Delta x \frac{\partial}{\partial x} \left( \sinh \left( \Delta x \frac{\partial r}{\partial x} \right) \cos \left( \Delta x \frac{\partial \phi}{\partial x} \right) \right) \right] - \\ & \quad - \operatorname{Im}F \left[ 2\sinh \left( \Delta x \frac{\partial r}{\partial x} \right) \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right) + \right. \\ & \quad \left. + \Delta x \frac{\partial}{\partial x} \left( \cosh \left( \Delta x \frac{\partial r}{\partial x} \right) \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right) \right) \right], \\ \\ \frac{\partial \phi}{\partial t} &= \omega + Ne^{2r} + \\ & \quad + \operatorname{Re}F \left[ 2\sinh \left( \Delta x \frac{\partial r}{\partial x} \right) \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right) + \right. \\ & \quad \left. + \Delta x \frac{\partial}{\partial x} \left( \cosh \left( \Delta x \frac{\partial r}{\partial x} \right) \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right) \right) \right] + \\ & \quad + \operatorname{Im}F \left[ 2\cosh \left( \Delta x \frac{\partial r}{\partial x} \right) \cos \left( \Delta x \frac{\partial \phi}{\partial x} \right) + \right. \\ & \quad \left. + \Delta x \frac{\partial}{\partial x} \left( \sinh \left( \Delta x \frac{\partial r}{\partial x} \right) \cos \left( \Delta x \frac{\partial \phi}{\partial x} \right) \right) \right]. \end{aligned} \right. \quad (22) $$

It is interesting to note that (22) includes a standard Kuramoto PDE derived in Section II-A as a particular case. Indeed,
assuming $r = r_0 = \text{const}$ both in space and time and assuming $F \in \mathbb{R}$, one gets an equation for $\phi$ as

$$ \frac{\partial \phi}{\partial t} = \omega + N e^{2r_0} + F \Delta x \frac{\partial}{\partial x} \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right), \quad (23) $$

which exactly coincides with (4) by replacing $\omega$ with $\omega + Ne^{2r_0}$.

Similar to Section II-B, we are interested in possible synchronized solutions of (21) and conditions for their existence and stability. A synchronized solution is a solution to (21) such that $\partial z/\partial t = i\bar{\omega}$, where $\bar{\omega}$ is a synchronization frequency. Thus we are interested in the question when such a solution $z = z(x)$ exists for some $\bar{\omega}$. Then the condition for synchronization is

$$ 0 = \Gamma + i(\omega - \bar{\omega}) - (S - iN)e^{2\mathrm{Re}z} + F \left[ 2\cosh\left(\Delta x \frac{\partial z}{\partial x}\right) + \Delta x \frac{\partial}{\partial x} \sinh\left(\Delta x \frac{\partial z}{\partial x}\right) \right], \quad (24) $$

or in terms of $r(x)$ and $\phi(x)$

$$ \left\{ \begin{array}{l} 0 = \mathrm{Re} (F^{-1} [\Gamma + i(\omega - \bar{\omega}) - (S - iN)e^{2r}] ) + \\ \quad + \left[ 2 \cosh \left( \Delta x \frac{\partial r}{\partial x} \right) \cos \left( \Delta x \frac{\partial \phi}{\partial x} \right) + \right. \\ \quad \left. + \Delta x \frac{\partial}{\partial x} \left( \sinh \left( \Delta x \frac{\partial r}{\partial x} \right) \cos \left( \Delta x \frac{\partial \phi}{\partial x} \right) \right) \right], \\ \\ 0 = \mathrm{Im} (F^{-1} [\Gamma + i(\omega - \bar{\omega}) - (S - iN)e^{2r}] ) + \\ \quad + \left[ 2 \sinh \left( \Delta x \frac{\partial r}{\partial x} \right) \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right) + \right. \\ \quad \left. + \Delta x \frac{\partial}{\partial x} \left( \cosh \left( \Delta x \frac{\partial r}{\partial x} \right) \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right) \right) \right]. \end{array} \right. \quad (25) $$

Note that we divided the equation by $F$ before splitting real and imaginary parts such that the hyperbolic functions take the simplest form.

The exponential term $e^{2r}$ in (25) can be removed by combining two equations together. Using amplitude-phase notation we can introduce $f = |F|$, $\beta = \arg F$, $G = |S + iN|$ and $\gamma = \arg(S + iN)$. With this notation

$$ \mathrm{Re} (F^{-1} (S - iN)) = \frac{G}{f} \cos (\gamma + \beta), $$

$$ \mathrm{Im} (F^{-1} (S - iN)) = -\frac{G}{f} \sin (\gamma + \beta), $$

therefore defining $A = \tan(\gamma + \beta)$, multiplying the first equation in (25) by $A$ and summing it with the second one we obtain

$$ \begin{aligned} & A \Delta x \frac{\partial}{\partial x} \left[ \cos \left( \Delta x \frac{\partial \phi}{\partial x} \right) \sinh \left( \Delta x \frac{\partial r}{\partial x} \right) \right] + \\ & + \Delta x \frac{\partial}{\partial x} \left[ \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right) \cosh \left( \Delta x \frac{\partial r}{\partial x} \right) \right] + \\ & + 2 \left[ A \cos \left( \Delta x \frac{\partial \phi}{\partial x} \right) \cosh \left( \Delta x \frac{\partial r}{\partial x} \right) + \right. \\ & \qquad \left. + \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right) \sinh \left( \Delta x \frac{\partial r}{\partial x} \right) \right] + B = 0, \end{aligned} \quad (26) $$

where

$$ B = \frac{1}{f \cos(\gamma + \beta)} [\cos \gamma (\omega - \bar{\omega}) + \sin \gamma \Gamma]. \quad (27) $$

Therefore, the synchronization condition is equivalent to (26) combined with one of the equations in (25) to determine connection between $r$ and $\phi$. Equilibrium solutions to (26) represent synchronized states of the original PDE (21). In the following sections we will characterize these solutions for different configurations and provide stability conditions.

## V. IDENTICAL OSCILLATORS CASE

In this section we will focus on the case when the ring consists of oscillators having identical parameters. Intuitively it is clear that in this case there exists a solution where all oscillators share the same amplitude $r$ and the same phase $\phi$. However it appears that depending on the number of oscillators and their parameters there can be more solutions, and that their stability properties are not trivial.

For the search of equilibrium solutions let us assume that in the synchronized state the amplitudes of oscillators $r$ are slowly varying in space, namely $\partial r/\partial x \approx 0$. We will further show that in the identical oscillators case the amplitude is indeed constant in space. Due to the assumption $\partial r/\partial x \approx 0$, we can use $\sinh(\Delta x \partial r/\partial x) \approx 0$ and $\cosh(\Delta x \partial r/\partial x) \approx 1$. With these simplifications the equation (26) depends only on $\phi(x)$ and thus can be solved independently:

$$ \Delta x \frac{\partial}{\partial x} \sin \left( \Delta x \frac{\partial \phi}{\partial x} \right) + 2A \cos \left( \Delta x \frac{\partial \phi}{\partial x} \right) + B = 0. \quad (28) $$

If the parameters of oscillators would be non-identical, the equation (28) would be very difficult to solve analytically since $A$ and $B$ are varying functions of space (it can be equivalently converted to the Abel equation of the second kind which has no closed-form solution). Therefore in this section we assume $A$ and $B$ being constant. A more general scenario of a piecewise constant functions $A$ and $B$ will be covered in the next section.

For constant $A$ and $B$ equation (28) is separable. We can notice that it depends only on the derivative of $\phi(x)$, not on the phase itself. Define $\theta = \Delta x \partial \phi / \partial x$. A physical meaning of $\theta$ is a difference in phases between two consecutive oscillators. With this definition, (28) becomes

$$ \frac{\cos \theta}{-B - 2A \cos \theta} d\theta = \frac{1}{\Delta x} dx. \quad (29) $$

The general solution to the equation (29) is described in Section VI. For now let us note that by the structure of (29) any non-constant continuous solution $\theta(x)$ would be monotonic. Further, apart from being a solution to (28), synchronization means that the solution $\phi(x)$ is a continuous angle, thus $\phi(x + 2\pi) - \phi(x) = 2\pi k$ for some $k \in \mathbb{Z}$. This implies two conditions which $\theta$ should satisfy:

$$ \begin{array}{l} 1). \quad \theta(x) = \theta(x + 2\pi) \quad \forall x \in \mathbb{R}, \\ 2). \quad \displaystyle\int_{0}^{2\pi} \theta(x) \, dx = 2\pi \Delta x k \quad \text{for some } k \in \mathbb{Z}. \end{array} \quad (30) $$

In the case of identical oscillators with constant $A$ and $B$ it appears that the only possible solution to (28) is a constant
one. Indeed, non-constant solutions should be monotone with respect to coordinate, however the first condition in (30) requires $\theta$ to be periodic, which is not possible if $\theta$ is not constant. Thus looking at the equation (28) we see that for constant A and B there is a simple solution $\cos \theta = -B/(2A)$, which means that all possible synchronized solutions for (28) are given by

$$ \theta = \arccos\left(-\frac{B}{2A}\right). \qquad (31) $$

### A. Equilibrium points

Recall that $\theta = \Delta x \frac{\partial \phi}{\partial x}$. Since $\phi(x)$ is a phase, it is defined up to a constant. Assuming $\phi(x) = 0$ at $x = 0$, using (31) and the definitions of A and B, the solution for $\phi(x)$ is a linear function

$$ \phi(x) = \frac{x}{\Delta x} \arccos \left( - \frac{\cos \gamma (\omega - \bar{\omega}) + \sin \gamma \Gamma}{2f \sin(\beta + \gamma)} \right). \quad (32) $$

Note that $\bar{\omega}$ is a synchronization frequency and is still unknown in this equation.

Position $x$ is itself defined on a ring, thus $x \in [0, 2\pi)$. Moreover, since the equilibrium solution is a periodic function, $\phi(2\pi)$ should also be a multiple of $2\pi$. We can define $k \in \mathbb{Z}_+$ such that $\phi(2\pi) = 2\pi k$. Since $k$ determines the number of phase turns along the ring, we call it the winding number. Therefore, the solution can exist for any $\bar{\omega}$ such that

$$ k = \frac{1}{\Delta x} \arccos \left( - \frac{\cos \gamma (\omega - \bar{\omega}) + \sin \gamma \Gamma}{2f \sin(\beta + \gamma)} \right), \quad k \in \mathbb{Z}_+. $$

The case $k = 0$ corresponds to an in-phase synchronized system, meaning phases of all oscillators coincide, while the case $k = 1$ corresponds to the state where the phases of the oscillators do a round turn along the ring. It is clear that in general the phase difference between neighbours is

$$ \theta^* = k\Delta x. $$

Note also that the system is symmetric for simultaneous substitution $k \to -k$ and $x \to -x$, thus phases can turn both clockwise and counter-clockwise along the ring.

The principal branch of arccos has a range of values $[0, \pi]$, therefore $k$ should satisfy $|k| \le \pi/\Delta x$ (other solutions will just copy the ones included in this range due to periodicity). Since $\Delta x$ is defined as the distance between two oscillators and is assumed to be constant, $\Delta x = 2\pi/n$, where $n$ is the number of oscillators in the system. Thus $|k| \le n/2$, with $|k| = n/2$ corresponding to the case when two neighbor oscillators are in anti-phase.

The synchronization frequency is thus given by

$$ \bar{\omega} = \omega + \tan \gamma \Gamma + 2f \frac{\sin(\beta + \gamma)}{\cos \gamma} \cos(k\Delta x) \quad (33) $$

for $k \in \{0, ..., n/2\}$. In particular, depending on the sign of $2f \cdot \sin(\beta + \gamma)/\cos \gamma$, the in-phase synchronized state is either the fastest or the slowest one.

### B. Stability analysis

Assume the equilibrium solution is given by (32) with the frequency (33) for $k \in \{0, ..., n/2\}$. We want to study for which of these $k$ the solution is stable.

Define $z^*(x,t) = r^* + i\phi^*(x) + i\bar{\omega}t$ to be an equilibrium solution for (21). Thus $\phi^*(x)$ is defined by (32) for a chosen $k$, $\bar{\omega}$ is a frequency of synchronized solution (33), and a constant $r^*$ can be found from (24) by taking its real part:

$$ e^{2r^*} = \frac{\Gamma + 2 \operatorname{Re} F \cos(k\Delta x)}{S} \quad (34) $$

Note that exponential should be positive to be well defined, therefore we require

$$ \Gamma + 2f \cos \beta \cos(k\Delta x) > 0. $$

Now let us define a deviation from the equilibrium solution $\bar{z}(x,t) = z(x,t) - z^*(x,t)$. It is governed by a difference of (21) for $z(x,t)$ and for $z^*(x,t)$, taking into account (24). Assuming $\bar{z}(x,t)$ is small, linearization of (21) around $z^*(x,t)$ is given by

$$ \begin{aligned} \frac{\partial \bar{z}}{\partial t} = & -2(S - iN)e^{2r^*} \operatorname{Re} \bar{z} + 2F \sinh \left( \Delta x \frac{\partial z^*}{\partial x} \right) \Delta x \frac{\partial \bar{z}}{\partial x} + \\ & + F \Delta x \frac{\partial}{\partial x} \left[ \cosh \left( \Delta x \frac{\partial z^*}{\partial x} \right) \Delta x \frac{\partial \bar{z}}{\partial x} \right]. \end{aligned} \quad (35) $$

Using

$$ \Delta x \frac{\partial z^*}{\partial x} = i \Delta x \frac{\partial \phi^*}{\partial x} = i \theta^* = i k \Delta x, $$

we get

$$ \cosh \left( \Delta x \frac{\partial z^*}{\partial x} \right) = \cos(k\Delta x), \quad \sinh \left( \Delta x \frac{\partial z^*}{\partial x} \right) = i \sin(k\Delta x), $$

which can be substituted in (35), resulting in

$$ \begin{aligned} \frac{\partial \bar{z}}{\partial t} = & -2(S - iN)e^{2r^*} \operatorname{Re} \bar{z} + 2iF \Delta x \sin(k\Delta x) \frac{\partial \bar{z}}{\partial x} + \\ & + F \Delta x^2 \cos(k\Delta x) \frac{\partial^2 \bar{z}}{\partial x^2}. \end{aligned} \quad (36) $$

Separating (36) into real and imaginary parts $\bar{z} = \tilde{r} + i\tilde{\phi}$ and using $F = f e^{i\beta}$:

$$ \left\{ \begin{aligned} \frac{\partial \tilde{r}}{\partial t} = & -2Se^{2r^*} \tilde{r} - 2f \Delta x \sin \beta \sin(k\Delta x) \frac{\partial \tilde{r}}{\partial x} - \\ & - 2f \Delta x \cos \beta \sin(k\Delta x) \frac{\partial \tilde{\phi}}{\partial x} + \\ & + f \Delta x^2 \cos \beta \cos(k\Delta x) \frac{\partial^2 \tilde{r}}{\partial x^2} - \\ & - f \Delta x^2 \sin \beta \cos(k\Delta x) \frac{\partial^2 \tilde{\phi}}{\partial x^2}, \\ \frac{\partial \tilde{\phi}}{\partial t} = & 2Ne^{2r^*} \tilde{r} + 2f \Delta x \cos \beta \sin(k\Delta x) \frac{\partial \tilde{r}}{\partial x} - \\ & - 2f \Delta x \sin \beta \sin(k\Delta x) \frac{\partial \tilde{\phi}}{\partial x} + \\ & + f \Delta x^2 \sin \beta \cos(k\Delta x) \frac{\partial^2 \tilde{r}}{\partial x^2} + \\ & + f \Delta x^2 \cos \beta \cos(k\Delta x) \frac{\partial^2 \tilde{\phi}}{\partial x^2}. \end{aligned} \right. \quad (37) $$

System (37) is a system of linear equations, thus the method of separation of variables can be applied to solve it. Moreover, it
is homogeneous, thus the basis functions should be exponen-
tial. Therefore stability of (37) can be checked by substituting
exponential basis functions

$$
\tilde{r} = r_0 e^{\lambda t} e^{imx}, \quad \tilde{\phi} = \phi_0 e^{\lambda t} e^{imx} \quad (38)
$$

for some $\lambda \in \mathbb{C}$ and $m \in \mathbb{Z}$, since basis should be periodic in $x$ along the ring. For asymptotic stability there should exist no solution of (37) with $\text{Re}\,\lambda > 0$. Substituting (38) in (37) one gets

$$
\left\{
\begin{array}{l}
\lambda r_0 = -2Se^{2r^*}r_0 - \\
\quad -2f\Delta x \sin\beta \sin(k\Delta x)imr_0 - 2f\Delta x \cos\beta \sin(k\Delta x)im\phi_0 - \\
\quad -f\Delta x^2 \cos\beta \cos(k\Delta x)m^2r_0 + f\Delta x^2 \sin\beta \cos(k\Delta x)m^2\phi_0, \\
\\
\lambda \phi_0 = 2Ne^{2r^*}r_0 + \\
\quad + 2f\Delta x \cos\beta \sin(k\Delta x)imr_0 - 2f\Delta x \sin\beta \sin(k\Delta x)im\phi_0 - \\
\quad -f\Delta x^2 \sin\beta \cos(k\Delta x)m^2r_0 - f\Delta x^2 \cos\beta \cos(k\Delta x)m^2\phi_0.
\end{array}
\right.
\tag{39}
$$

Define

$$
\begin{align*}
P &= f \Delta x^2 \cos \beta \cos(k \Delta x) m^2 + 2i f \Delta x \sin \beta \sin(k \Delta x) m, \\
Q &= -f \Delta x^2 \sin \beta \cos(k \Delta x) m^2 + 2i f \Delta x \cos \beta \sin(k \Delta x) m.
\end{align*}
$$

Note that as $m \to \infty$, first terms become dominating. With
the help of these functions and with $\bar{S} = 2Se^{2r^*} > 0$ and $\bar{N} =
2Ne^{2r^*}$, (39) becomes

$$
\lambda \begin{pmatrix} r_0 \\ \phi_0 \end{pmatrix} = \begin{pmatrix} -P - \bar{S} & -Q \\ \bar{N} + Q & -P \end{pmatrix} \begin{pmatrix} r_0 \\ \phi_0 \end{pmatrix}, \quad (40)
$$

thus we are interested in the eigenvalues of the matrix in (40).
It is trivial to show that they are given by

$$
\lambda = \frac{1}{2} \left( -2P - \bar{S} \pm \sqrt{(2P + \bar{S})^2 - 4P(P + \bar{S}) - 4Q(Q + \bar{N})} \right). \quad (41)
$$

Taking $m = 0$, one of the eigenvalues becomes zero, corre-
sponding to the fact that the phase is defined up to a constant,
and the other eigenvalue is $-\bar{S}$.

Further assume $m \neq 0$ and thus $P, Q \neq 0$. Condition for stability Re $\lambda < 0$ translates as

$$
\mathrm{Re}(2P + \bar{S}) > \mathrm{Re} \sqrt{(2P + \bar{S})^2 - 4P(P + \bar{S}) - 4Q(Q + \bar{N})}. \quad (42)
$$

In particular, as $m \to \infty$, ReP should be positive. Now for simplicity define

$$
H = (2P + \bar{S})^2, \quad D = 4P(P + \bar{S}) + 4Q(Q + \bar{N}).
$$

Using complex relation for any $c \in \mathbb{C}$

$$
2(\mathrm{Re}c)^2 = \mathrm{Re}(c^2) + |c|^2 \quad (43)
$$

and taking square of inequality (42), it becomes

$$
|H| + \mathrm{Re}D > |H - D|, \quad (44)
$$

in particular Re $D > -|H|$. Taking square once more we get

$$
(\mathrm{Im}D)^2 - 2 \mathrm{Im} H \mathrm{Im} D - 2|H|(\cos \nu + 1) \mathrm{Re} D < 0, \quad (45)
$$

where $\nu = \arg H$. By (43) $|H|(\cos \nu + 1) = 2(\mathrm{Re}\sqrt{H})^2$, and thus (45) means that

$$
(\mathrm{Im} D)^2 - 2 \mathrm{Im} H \mathrm{Im} D < 4 \mathrm{Re} D (\mathrm{Re} \sqrt{H})^2. \quad (46)
$$

Defining sequences $h_m$, $d_m$ for $m \in \mathbb{Z} \setminus \{0\}$ as

$$
h_m = f \Delta x^2 \cos(k \Delta x) m^2 \quad \text{and} \quad d_m = 2f \Delta x \sin(k \Delta x) m, \quad (47)
$$

we can express P, Q, D and H as

$$
P = h_m \cos \beta + i d_m \sin \beta, \quad Q = -h_m \sin \beta + i d_m \cos \beta,
$$

$$
D = 4(h_m^2 - d_m^2 + \bar{S}h_m \cos \beta - \bar{N}h_m \sin \beta) + \\
\qquad + 4i(\bar{S}d_m \sin \beta + \bar{N}d_m \cos \beta),
$$

$$
H = (4h_m^2 \cos^2 \beta - 4d_m^2 \sin^2 \beta + \bar{S}^2 + 4h_m \bar{S} \cos \beta) + \\
\qquad + 4i(\bar{S}d_m \sin \beta + 2h_m d_m \sin \beta \cos \beta),
$$

which then being inserted in (46) results in

$$
\begin{align}
& (\bar{S}d_m \sin \beta + \bar{N}d_m \cos \beta) \cdot \nonumber \\
& \cdot (\bar{N}d_m \cos \beta - \bar{S}d_m \sin \beta - 4h_m d_m \sin \beta \cos \beta) < \tag{48} \\
& < (\bar{S} + 2h_m \cos \beta)^2 (h_m^2 - d_m^2 + \bar{S}h_m \cos \beta - \bar{N}h_m \sin \beta). \nonumber
\end{align}
$$

Finally, defining $\bar{G} = 2e^{2r^*}G$ and using $\bar{S} = \bar{G}\cos\gamma$ and $\bar{N} = \bar{G}\sin\gamma$, we get

$$
\begin{equation}
\begin{aligned}
d_m^2 \bar{G} \sin(\gamma + \beta) (\bar{G} \sin(\gamma - \beta) - 4h_m \sin \beta \cos \beta) < \\
< (\bar{G} \cos \gamma + 2h_m \cos \beta)^2 (h_m^2 - d_m^2 + \bar{G}h_m \cos(\gamma + \beta)).
\end{aligned}
\tag{49}
\end{equation}
$$

Recall that by (42) Re $P > 0$, and substituting $P$ and $h_m$ from (47) we imply also that $\cos(k\Delta x)\cos\beta > 0$. Thus we just proved a theorem, which can be seen as the main result of this paper:

**Theorem 2.** *Necessary and sufficient condition for stability of the synchronized solution (32) with the frequency (33) for $k \in \{0, ..., n/2\}$ for the system (21) with constant parameters is given by the inequality (49) for $m \in \mathbb{Z} \setminus \{0\}$ together with the requirement $\cos(k\Delta x)\cos\beta > 0$.*

Due to the dependence of (49) on $m$ it is difficult to
check this condition explicitly. Therefore we will state several
corollaries for particular values of the winding number $k$,
providing explicit inequalities to check. These corollaries can
be then directly used for particular systems establishing a
guaranteed stability of synchronized solutions.

**Corollary 1.** *Necessary and sufficient conditions for in-phase synchronization are given by*

$$
\cos \beta > 0, \quad \cos(\gamma + \beta) > -\frac{f \Delta x^2}{2e^{2r^*} G}.
$$

*Proof.* Indeed, in-phase equilibrium solution satisfies $k=0$, thus by (47) $d_m = 0$ and $h_m = f\Delta x^2 m^2 > 0$. From the second condition of Theorem 2 we recover $\cos\beta > 0$. Finally, (49) with $d_m = 0$ requires right-hand terms to be greater than zero, which is just $h_m(h_m + \bar{G}\cos(\gamma + \beta)) > 0$. Since this is always true as $h_m \to \infty$ with $m \to \pm\infty$, it is enough to satisfy this inequality for $m = \pm 1$, leading to $\cos(\gamma + \beta) > -f\Delta x^2/\bar{G}$. $\square$

Notice that conditions required in Corollary 1 as well as in
all other corollaries below immediately ensure the existence
of exponential representation of the amplitude of oscillations
defined in (34).
**Corollary 2.** Necessary and sufficient conditions for anti-phase synchronization are

$$ \cos \beta < 0, \quad \cos(\gamma + \beta) < \frac{f \Delta x^2}{2e^{2r^*} G}. $$

*Proof.* The proof follows the same steps as the previous one, switching the sign of $h_m$. □

**Corollary 3.** Sufficient conditions for synchronization with $\sin(k\Delta x) \neq 0$ are given by

$$ \cos(k\Delta x) \cos \beta > 0 \qquad (50) $$

together with

$$ \begin{gathered} \Upsilon < \frac{\Delta x^2}{4} \cot(k\Delta x)^2 + \frac{Ge^{2r^*} \cos(\gamma + \beta)}{2f} \frac{\cos(k\Delta x)}{\sin(k\Delta x)^2} - 1, \\ \text{where } \Upsilon = \begin{cases} 0, & \cos^2 \beta \le \cos^2 \gamma, \\ \frac{\cos^2 \beta}{\cos^2 \gamma} - 1, & \cos^2 \beta > \cos^2 \gamma. \end{cases} \qquad (51) \end{gathered} $$

Condition (50) is also a necessary condition for stability.

*Proof.* First, condition (50) repeats the second condition of Theorem 2. Further, since $\sin(k\Delta x) \neq 0$, $d_m$ is non-zero. Divide (49) by $d_m^2$ and by $(\bar{G}\cos\gamma + 2h_m\cos\beta)^2$, obtaining

$$ \frac{\bar{G}\sin(\gamma + \beta)(\bar{G}\sin(\gamma - \beta) - 4h_m \sin \beta \cos \beta)}{(\bar{G}\cos \gamma + 2h_m \cos \beta)^2} < \frac{h_m^2 - d_m^2 + \bar{G}h_m \cos(\gamma + \beta)}{d_m^2}. \qquad (52) $$

Inserting the definitions of $h_m$ and $d_m$ from (47) we see that right-hand side of (52) is strictly increasing with $m^2$, therefore it can be simplified by setting $m^2 = 1$ as in the worst-case, thus obtaining right-hand side of (51).

Now, to find sufficient conditions for satisfaction of (52), let us bound the left-hand side from above. For this we will use the following Lemma:

**Lemma 1.** Function $f(x)$, defined as

$$ f(x) = \frac{V + \mu x}{(U + x)^2} \qquad (53) $$

with $U > 0$ and $x > 0$ is bounded from above by

$$ f(x) \le \begin{cases} 0, & V \le 0 \text{ and } \mu \le 0, \\ \frac{V/U^2}{4\mu U - 4V}, & V > 0 \text{ and } U\mu \le 2V, \\ \frac{\mu^2}{4\mu U - 4V}, & \mu > 0 \text{ and } U\mu > 2V. \end{cases} \qquad (54) $$

The proof of this lemma can be found in Appendix A. We apply this lemma to the left-hand side of (52), with $U = \bar{G}\cos\gamma > 0$, $x = 2h_m\cos\beta > 0$, $V = \bar{G}^2\sin(\gamma + \beta)\sin(\gamma - \beta)$ and $\mu = -2\bar{G}\sin(\gamma + \beta)\sin\beta$, obtaining (51). Note that due to the trigonometric properties the conditions $U > 0$, $\mu > 0$ and $U\mu > 2V$ are contradicting by definitions of variables, thus only the first two cases of (54) are present in (51). Further, $V \le 0$ and $U > 0$ implies $\mu \le 0$, while $V > 0$ and $U > 0$ implies $U\mu < 2V$, thus it is sufficient to check only $V$ in (54). □

Assume $\cos\beta > 0$ such that the in-phase solution is stable. Then the second condition of Theorem 2 requires that

<figure><img src="image_6.png" /><figcaption>$k=0$, in-phase stable</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>$k=1$, stable</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>$k=2$, stable</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>$k=3$, unstable</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>$k=4$, unstable</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>$k=5$, anti-phase unstable</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d)</figcaption></figure>

<caption>Figure 6. Synchronized solutions for system of $n = 50$ coupled identical spin-torque oscillators. **Top row:** analytic results for PDE (21). **(a):** diagram of possible regimes by Corollaries 1 and 3. Color code denotes the highest guaranteed existing regime, *chaotic* means that no stable solution exists. **(b):** Synchronization frequency $\omega$ by (33) for different $k$ depending on $\beta$ for $f = 0.75$. **Bottom row:** numerical simulation of (19). **(c):** diagram of numerically established regimes. **(d):** numerically measured synchronization frequency.</caption>

$\cos(k\Delta x) > 0$. Thus for the stability $k\Delta x$ should be smaller than $\pi/2$. This means $k < n/4$, where $n$ is the number of oscillators. In particular, the phase difference between two neighbouring oscillators should be smaller than $\pi/2$. Also this means that to observe a state with $k = 1$ one needs at least 5 coupled oscillators, and to observe higher-order states one needs at least 9 oscillators. As an example, all possible states in the system with 10 oscillators are shown in Fig. 5.
### C. Numerical simulation

To compare predictions from the previous section we performed a numerical simulation of system of $n = 50$ coupled spin-torque oscillators. Simulation parameters were chosen according to [15], namely, we set $\omega = 6.55 \cdot 2\pi$, $N = -3.82 \cdot 2\pi$, $\Gamma_G = 0.375 \cdot 2\pi$ (all those measured in radians per nanosecond), $Q = -0.24$ and $\sigma = 5.48 \cdot 10^{-13} \cdot 2\pi$ rad/nanosec/A*m$^2$ for (16). In this case the critical current density which is required to start oscillations is $I_c = \Gamma_G/\sigma = 684.3 \cdot 10^9$ A/m$^2$. In our experiments we use a larger current density $I = 1.5I_c$ to observe steady state oscillations. With this setup the parameters of oscillator (17) are $\Gamma = 1.1781$ and $S = 2.9688$. Further, using definitions $G = |S + iN|$ and $\gamma = \arg(S + iN)$, we get $G = 24.1847$ and $\gamma = -82.95^\circ$.

Due to large negative $\gamma$ conditions in Corollaries 1 and 3 are not easy to satisfy. We can check which stable synchronized solutions are admitted by the coupled system depending on different coupling parameter $F$. Comparison between analytic predictions and numerical simulation results is shown in Fig. 6. We take different couplings $F = fe^{i\beta}$ with $f$ changing from 0 to 10 and $\beta$ changing from -0.2 to 0.2 radians. For each set of parameters we check the highest $k$ for which conditions in Corollaries 1 and 3 are satisfied. These results are depicted in the diagram Fig. 6a. Further, we compare them with experimental results by simulating the original ODE system (19). We initialize all oscillators in this system using an amplitude $\sqrt{p_i} = \sqrt{\Gamma}/S$ and a phase $\phi_i = ik\Delta x$ for the $i$-th oscillator, such that the phase makes $k$ turns along the ring. Finally a small gaussian noise with a standard deviation of 0.05 is added to phases. The system is simulated for 5000 nanoseconds (corresponding roughly to 15000 periods of oscillation for $f = 0.75$). When simulation ends, we check if the system remained stable or it diverged from the corresponding equilibrium solution. The obtained highest possible stable regimes are depicted in the diagram Fig. 6c. Comparing it with the diagram Fig. 6a, we see that the analytic prediction almost perfectly reconstructs the experimental diagram, with deviations probably being attributed to the inaccuracies in numerical stability check.

Finally we compare synchronization frequency $\omega$ predicted by (33) with the one measured in simulation. To measure synchronization frequency in simulation we first notice that for every agent oscillating with constant amplitude its immediate frequency can be found as $\omega \approx \text{Im}(\dot{c}/c)$. Then we average this frequency over all agents and over last 1000 nanoseconds. The measured synchronization frequency for $f = 0.75$ and for $\beta \in [-0.1, 0.2]$ is depicted in Fig. 6d. It is clear that for the higher regimes for $k = 1$ and $k = 2$ stable solutions exist only for sufficiently high values of $\beta$. Comparing measured frequency with analytically predicted by (33) in Fig. 6b one can see that the trends and relative frequency differences between different regimes are reproduced correctly and that the measured frequency is about 0.1 rad/nanosec higher than the predicted one. This effect vanishes for higher values of $f$. This mismatch can originate from the fact that the analytic prediction was found for the PDE model (21), while the simulation was performed for the original ODE system (19).

<figure><img src="image_16.png" /><figcaption>Figure 7. Function $g(\theta, J)$ defined in (56) with respect to $J$ for different values of $\theta$. Colors denote different branches in (56).</figcaption></figure>

## VI. NON-IDENTICAL OSCILLATORS

In the previous section we assumed that all oscillators are identical and that the solutions' magnitude is constant in space. In this section we relax a requirement on homogeneity but keep the assumption that $\partial r/\partial x \approx 0$. Instead of analysing systems with constant parameters $A$ and $B$ along the ring we assume that these parameters are piecewise constant.

It was shown in Section V that in the small magnitude variation case the synchronization condition is equivalent to the equation (28), and for constant parameters $A$ and $B$ it is represented by a separable equation (29). Let us define $J = B/A$. Performing integration of (29) it is trivial to show that all solutions to (29) are given in a form

$$A \frac{x}{\Delta x} + C = g(\theta, J), \quad (55)$$

where $C$ is an integration constant, and $g(\theta, J)$ is a special function parametrized by $J$ which is given by

$$g(\theta, J) = \begin{cases} \frac{J}{2\sqrt{4-J^2}} \ln \left| \frac{1 + \left( \frac{2-J}{\sqrt{4-J^2}} \tan \frac{\theta}{2} \right)}{1 - \left( \frac{2-J}{\sqrt{4-J^2}} \tan \frac{\theta}{2} \right)} \right| - \frac{1}{2}\theta, & |J| < 2, \\ \frac{J}{\sqrt{J^2 - 4}} \arctan \left( \frac{J-2}{\sqrt{J^2-4}} \tan \frac{\theta}{2} \right) - \frac{1}{2}\theta, & |J| > 2, \\ \frac{1}{2} \tan \frac{\theta}{2} - \frac{1}{2}\theta, & J=2, \\ -\frac{1}{2} \cot \frac{\theta}{2} - \frac{1}{2}\theta, & J=-2. \end{cases} \quad (56)$$

It is interesting to note that there is a complex relation between arctangent and logarithm functions

$$\arctan s = -\frac{i}{2} \ln \left( \frac{1+is}{1-is} \right), \quad (57)$$

which means that the first two cases in (56) are essentially the same. In fact, the definition (56) defines a piecewise continuous function with at most two singularities with respect to $J$, see Fig. 7. In particular (55) means that the constant solution (31) is captured by the singularity approaching $J = -2$ from the left. Note further that $g(\theta, J)$ is an odd function with respect to $\theta$.

Using the solution (55) it becomes possible to analyse systems with several different types of oscillators. Here for simplicity we will focus on the case of two types of oscillators. The first type of oscillators has a set of parameters $\omega_1, N_1, \Gamma_1,$
<figure><img src="image_17.png" /><figcaption>Figure 8. Examples of schematic representations of a ring with $n = 20$ oscillators with two different oscillator types placed periodically. K is a number of periods.</figcaption></figure>

$S_1$ and $F_1$, and similarly the second type has a corresponding set of its own parameters. We further assume that oscillators' types are repeated $K$ times along the ring, and that every continuous chunk of a particular oscillators' type consists of a fixed number of oscillators depending on its type (evidently this implies $K$ is a divisor of the number of oscillators $n$). This means that the type of oscillators is a periodic function on the ring with period $2\pi/K$. For example if $K=1$ this setup corresponds to one large set of oscillators of the first type followed by only one large set of oscillators of the second type, while if $K=n/2$ the types of oscillators alternate. We can define a set of switching points as $y_j$ for $j \in \{0, ..., 2K-1\}$, with $y_0 = 0$ and $y_j = j/2 \cdot 2\pi/K$ for even $j$. Finally, for odd $j$ we require $y_j - y_{j-1} = \text{const}$, thus the proportion of types is preserved. Oscillators placed in $[0, y_1) \cup [y_2, y_3) \cup ...$ are of the first type, and oscillators placed in $[y_1, y_2) \cup [y_3, y_4) \cup ...$ are of the second type. In particular this means that oscillators of the first type occupy proportion $y_1/y_2$ of the whole ring. Some possible examples of such distributions are schematically presented in Fig. 8.

Since oscillators are of different types, aggregated parameters $A$ and $B$ will have different values $A_1, A_2, B_1$ and $B_2$, leading to two different decision parameters $J_1$ and $J_2$. However an unknown synchronization frequency $\bar{\omega}$ should be common for both types, therefore by definition of $B$ in (27) we can write $J_1 = \bar{J}_1 + \tau_1 \bar{\omega}$ and $J_2 = \bar{J}_2 + \tau_2 \bar{\omega}$, where

$$ \bar{J}_1 = \frac{\cos \gamma_1 \omega_1 + \sin \gamma_1 \Gamma_1}{f_1 \sin(\gamma_1 + \beta_1)}, \quad \tau_1 = -\frac{\cos \gamma_1}{f_1 \sin(\gamma_1 + \beta_1)}, \quad (58) $$

with $\bar{J}_2$ and $\tau_2$ being defined in a similar way.

We are now interested in particular solutions $\theta(x)$ to (28). By (30) $\theta$ should be periodic. Since intervals of types of oscillators are equal, symmetry leads to the fact that $\theta$ should be periodic with period being equal to two intervals of different types of oscillators, namely $\theta(y_0) = \theta(y_2) = \theta(y_4) = ... = \theta(y_{2K-2})$. Further, one could expect to obtain continuous solutions, however performing numerical simulations of such systems we made an observation regarding possible synchronized solutions:

**Observation 1.** Solution $\theta(x)$ behaves continuously and monotonically in the first type domain and is constant with discontinuity in the interior in the second type domain. More-

over, solution endpoints are symmetric about zero, namely $\theta(y_0) = -\theta(y_1)$.

The set of all possible solutions is not covered only by those proposed by Observation 1, however each particular class of solutions heavily depends on the properties of the function (56) and thus requires special treatment. Further in this section we will stick to the class of solutions in agreement with Observation 1.

Defining $\theta^* = \theta(0)$ and assuming $\theta(y_1) = -\theta^*$ by Observation 1, we can compute (55) in points $x = y_0 = 0$ and $x = y_1$ for the first type of oscillators and subtract one from another, obtaining

$$ A_1 \frac{y_1}{\Delta x} = 2g(\theta^*, J_1), $$

where we used the fact that the function $g(\theta, J)$ is odd with respect to $\theta$. Substituting $J_1$ as in (58), we get a condition which should be satisfied for the first type of oscillators

$$ 2g(\theta^*, \bar{J}_1 + \tau_1 \bar{\omega}) - A_1 \frac{y_1}{\Delta x} = 0, \quad (59) $$

which have two unknowns: $\theta^*$ and $\bar{\omega}$. The second condition comes from the assumption that for the second type domain the solution is constant and thus it is determined by (31). Using it for the second type domain we get

$$ \theta^* = \arccos \left( -\frac{J_2}{2} \right). \quad (60) $$

Note that both $\theta^*$ and $-\theta^*$ are solutions to (31), which is consistent with Observation 1. Now, substituting $J_2$ by (58) in (60) and then substituting the result in (59) we obtain an equation with a single unknown $\bar{\omega}$:

$$ 2g \left( \arccos \left( -\frac{\bar{J}_2 + \tau_2 \bar{\omega}}{2} \right), \bar{J}_1 + \tau_1 \bar{\omega} \right) - A_1 \frac{y_1}{\Delta x} = 0. \quad (61) $$

This equation can be solved for $\bar{\omega}$ using numerical methods such as Newton method for example. Once $\bar{\omega}$ is known, we can find $J_1$ and $J_2$ by (58) and then compute $\theta^*$ by (60). The full solution on the first domain is then reconstructed by (55).

To determine the shape of solution $\theta(x)$ it remains only to find an exact position denoted by $y^* \in (y_1, y_2)$ where a discontinuous jump from $\theta^*$ to $-\theta^*$ happens in the second type domain. This position can be obtained if one recalls that $\theta = \Delta x \partial\phi/\partial x$ and thus integral of $\theta$ should have fixed value by (30) for some $k \in \mathbb{Z}$. In particular due to the periodic nature of the problem with $K$ periods we have

$$ \int_{0}^{y_2} \theta(x) dx = \frac{2\pi\Delta x}{K} k. \quad (62) $$

Since on the first type domain $\theta(x)$ is symmetric, its contribution to the integral is zero. Further, $\theta(x) = \theta^*$ on $x \in [y_1, y^*)$ and $\theta(x) = -\theta^*$ on $x \in (y^*, y_2]$, therefore (62) is just

$$ (2y^* - y_1 - y_2) \theta^* = \frac{2\pi\Delta x}{K} k, $$

which leads to

$$ y^* = \frac{\pi \Delta x}{K \theta^*} k + \frac{y_1 + y_2}{2}. \quad (63) $$

Thus the solution's shape $\theta(x)$ is fully reconstructed.
<figure><img src="image_18.png" /><figcaption>Figure 9. Comparison of numerical and analytical synchronized solutions of systems with $n = 500$ oscillators separated into two classes. Horizontal axis: index of oscillator. Vertical axis: phase difference between two consecutive oscillators in degrees. Yellow line denotes solution obtained by numerical simulation of (19), black line denotes analytic solution by (60)-(63). Parameters: (a): $K = 1$, $y_1/y_2 = 0.2$, $\Gamma_2 = 1.05 \cdot \Gamma_1$, $k = 3$. (b): $K = 4$, $y_1/y_2 = 0.6$, $N_2 = 1.03 \cdot N_1$, $k = -2$.</figcaption></figure>

**Remark 1.** Observation 1 assumes the first part of the solution behaves continuously and the second part is piecewise constant. In real system these parts can be interchanged, which depends on the obtained values of $J_1$ and $J_2$: for the continuous part $|J| > 2$, while for the piecewise constant part $|J| < 2$ (while they are both usually negative and close to -2).

**Remark 2.** Other types of solutions except those presented in Observation 1 are also possible if parameter variations are very high. In this case there is no piecewise constant domain and all solution's parts behave according to (55). It is then possible to formulate a system of nonlinear equations with several unknown variables which should be solved numerically. However we found that solutions to this system lie very close to singularities of $g(\theta, J)$, thus they cannot be found reliably by numerical methods without additional problem reformulation.

### A. Numerical simulation

To demonstrate how solutions to the synchronization condition (28) found by (60)-(63) approximate synchronized solutions of the original system (19) we performed numerical simulations of (19) with $n = 500$ oscillators being split into two types as it was described earlier in this section. Parameters of the first type of oscillators were taken the same as in Section V-C, and for the second type slight deviations in parameters were added. Oscillators were placed periodically on the ring with $K$ periods, thus there were $2K$ groups of oscillators as it was shown in Fig. 8. Each group of oscillators of the first type occupies $y_1/y_2$ proportion of the period of the length $y_2$, and each group of oscillators of the second type occupies $(y_2 - y_1)/y_2$ proportion. Numerical simulation was initialized in the same way as in Section V-C with $k$ denoting initial shift in phases of consecutive oscillators such that the phase makes $k$ turns along the ring.

We performed two simulations:

1. In the first simulation we altered damping parameter $\Gamma$ for the second type of oscillators such that $\Gamma_2 = \Gamma_1 \cdot 1.05$. We used only two types of oscillators, thus $K = 1$. The first type occupies only 20% of the whole ring, thus $y_1/y_2 = 0.2$. Finally, oscillators were initialized such that the phase makes $k = 3$ turns along the ring.

2. In the second simulation we changed frequency gain parameter $N$ for the second type of oscillators such that $N_2 = N_1 \cdot 1.03$. We used eight groups of oscillators, four of each type, thus $K = 4$. The first type occupies 60% of every period, thus $y_1/y_2 = 0.6$. In this simulation oscillators were initialized such that the phase makes $k = -2$ turns along the ring, rotating in opposite direction.

Results of the simulation are presented in Fig. 9. Simulation was performed for 2000 nanoseconds and then phase differences between consecutive oscillators were computed. The result was then compared with analytic predictions by (60)-(63). It is clear that the shape of solutions is reconstructed almost perfectly even though our analysis was based on the continualized PDE model of the network and a small magnitude variation assumption. Also, numerical simulations have shown that synchronized solutions are very fragile in a sense that small deviations in parameters result in very large differences in phases between consecutive oscillators, although the system still remains stable and the predictions are still correct.

## VII. CONCLUSION

Analysis of synchronization of spin-torque oscillators has a big practical importance since synchronous oscillations produce much more energy, therefore it is very important to realize when synchronized solutions do exist and what deviations in fabrication (which result in deviations in parameters) they do tolerate. In this paper we have shown how the continuation method can help in the analysis of this problem and then we derived some results which could be useful in practical applications. In particular we completely treated the case of identical oscillators, providing explicit formulas for equilibrium solutions and their stability conditions. For non-identical oscillators we analysed one particular class of possible equilibrium solutions, showing that its shape can be analytically reconstructed and thus opening new possibilities for more efficient modeling and future analysis of the system. Still, there are many questions that could be investigated in details regarding the system (19), its PDE approximation (21) and the synchronization condition (24). First, Corollary 3 for the general winding number $k$ in the case of identical oscillators gives only sufficient conditions on stability and probably more rigorous statements could be made based on Theorem 2. Second, a practically more important case of non-identical oscillators, discussed in Section VI, covers only a search for specific equilibrium solutions. Therefore more general equilibrium analysis and investigation of stability conditions could be performed. Finally, this work in general and the model (19) in particular were devoted to the analysis of the ring topology of oscillators. Study of more general topologies, especially 2-dimensional, would be of a great importance for practical applications.

# APPENDIX A

## PROOF OF THE LEMMA 1

The function $f(x)$ is defined for $x \in [0, +\infty)$, thus its supremum is achieved either at $x = 0$, $x = +\infty$ or at $f'(x) = 0$.
If $V \le 0$ and $\mu \le 0$, then the function is nonpositive with asymptotic value $f(+\infty) = 0$, thus we use 0 as a bound in this case. Let us now find its extremum:

$$f'(x) = \left( \frac{V + \mu x}{(U + x)^2} \right)' = \frac{\mu(U - x) - 2V}{(U + x)^3} = 0, \quad (64)$$

thus it is achieved at $x_{extr} = U - 2V/\mu$. Substituting it back in (64) we obtain

$$f(x_{extr}) = \frac{\mu^2}{4\mu U - 4V}. \quad (65)$$

Finally we notice that the extremum (64) is indeed maximum only if $\mu > 0$ and if $x_{extr} > 0$, otherwise the maximum is achieved at zero, $f(0) = V/U^2$. Therefore, combining the bounds together we get

$$f(x) \le \begin{cases} 0, & V \le 0 \text{ and } \mu \le 0, \\ \frac{V}{U^2}, & V > 0 \text{ and } U\mu \le 2V, \\ \frac{\mu^2}{4\mu U - 4V}, & \mu > 0 \text{ and } U\mu > 2V. \end{cases} \quad (66)$$

## ACKNOWLEDGMENT

The Scale-FreeBack project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement N 694209).

## REFERENCES

[1] I. Z. Kiss, Y. Zhai, and J. L. Hudson, "Emerging coherence in a population of chemical oscillators," *Science*, vol. 296, no. 5573, pp. 1676-1678, 2002.

[2] F. Dörfler, M. Chertkov, and F. Bullo, "Synchronization in complex oscillator networks and smart grids," *Proceedings of the National Academy of Sciences*, vol. 110, no. 6, pp. 2005-2010, 2013.

[3] G. Filatrella, A. H. Nielsen, and N. F. Pedersen, "Analysis of a power grid using a Kuramoto-like model," *The European Physical Journal B*, vol. 61, no. 4, pp. 485-491, 2008.

[4] S. H. Strogatz, D. M. Abrams, A. McRobie, B. Eckhardt, and E. Ott, "Crowd synchrony on the Millennium Bridge," *Nature*, vol. 438, no. 7064, pp. 43-44, 2005.

[5] F. Dörfler and F. Bullo, "Synchronization in complex networks of phase oscillators: A survey," *Automatica*, vol. 50, no. 6, pp. 1539-1564, 2014.

[6] S. Jafarpour and F. Bullo, "Synchronization of Kuramoto oscillators via cutset projections," *IEEE Transactions on Automatic Control*, vol. 64, no. 7, pp. 2830-2844, 2018.

[7] A. Kuznetsov, N. Stankevich, and L. Turukina, "Coupled van der Pol-Duffing oscillators: Phase dynamics and structure of synchronization tongues," *Physica D: Nonlinear Phenomena*, vol. 238, no. 14, pp. 1203-1215, 2009.

[8] A. Schmidt, T. Kasimatis, J. Hizanidis, A. Provata, and P. Hövel, "Chimera patterns in two-dimensional networks of coupled neurons," *Physical Review E*, vol. 95, no. 3, p. 032224, 2017.

[9] J. C. Slonczewski, "Current-driven excitation of magnetic multilayers," *Journal of Magnetism and Magnetic Materials*, vol. 159, no. 1-2, pp. L1-L7, 1996.

[10] L. Berger, "Emission of spin waves by a magnetic multilayer traversed by a current," *Physical Review B*, vol. 54, no. 13, p. 9353, 1996.

[11] S. Sani, J. Persson, S. M. Mohseni, Y. Pogoryelov, P. Muduli, A. Eklund, G. Malm, M. Käll, A. Dmitriev, and J. Åkerman, "Mutually synchronized bottom-up multi-nanocontact spin-torque oscillators," *Nature communications*, vol. 4, no. 1, pp. 1-7, 2013.

[12] M. Carpentieri and G. Finocchio, "Spintronic oscillators based on spin-transfer torque and spin-orbit torque," in *Handbook of Surface Science*. Elsevier, 2015, vol. 5, pp. 297-334.

[13] J. Persson, Y. Zhou, and J. Åkerman, "Phase-locked spin torque oscillators: Impact of device variability and time delay," *Journal of applied physics*, vol. 101, no. 9, p. 09A503, 2007.

[14] P. Muduli, Y. Pogoryelov, F. Mancoff, and J. Åkerman, "Modulation of individual and mutually synchronized nanocontact-based spin torque oscillators," *IEEE transactions on magnetics*, vol. 47, no. 6, pp. 1575-1579, 2011.

[15] C. Dieudonné, "Synchronization of a spin transfer oscillator to a RF current: mechanisms and room-temperature characterization," Ph.D. dissertation, Université Grenoble Alpes, 2015.

[16] A. Slavin and V. Tiberkevich, "Nonlinear auto-oscillator theory of microwave generation by spin-polarized current," *IEEE Transactions on Magnetics*, vol. 45, no. 4, pp. 1875-1918, 2009.

[17] D. Nikitin, C. Canudas-de Wit, and P. Frasca, "A continuation method for large-scale modeling and control: from ODEs to PDE, a round trip," *Conditionally accepted to IEEE Transactions on Automatic Control*, arXiv preprint arXiv:2101.10060, 2021.

[18] L. Tumash, C. Canudas-De-Wit, and M. L. Delle Monache, "Multi-directional continuous traffic model for large-scale urban networks," *HAL preprint hal-03236552*, 2021.

[19] D. Nikitin, C. Canudas-De-Wit, and P. Frasca, "Boundary control for stabilization of large-scale networks through the continuation method," *HAL preprint hal-03211021f*, 2021.

[20] M. D. Stiles and J. Miltat, "Spin-transfer torque and dynamics," *Spin dynamics in confined magnetic structures III*, pp. 225-308, 2006.

<figure><img src="image_19.png" /></figure>

Denis Nikitin received both the B.Sc. degree in 2016 and the M.Sc. degree in 2018 in Mathematics and Mechanics Faculty of Saint Petersburg State University, Saint Petersburg, Russia, specializing on the control theory and cybernetics. He won several international robotics competitions while being student and was a teacher of robotics in the Math&Phys Lyceum 239 in Saint Petersburg. He is currently a doctoral researcher at CNRS, GIPSA-Lab, Grenoble, France. His current research mainly focuses on control of large-scale systems.

<figure><img src="image_20.png" /></figure>

Carlos Canudas-de-Wit (F'16) was born in Villahermosa, Mexico, in 1958. He received the B.S. degree in electronics and communications from the Monterrey Institute of Technology and Higher Education, Monterrey, Mexico, in 1980, and the M.S. and Ph.D. degrees in automatic control from the Department of Automatic Control, Grenoble Institute of Technology, Grenoble, France, in 1984 and 1987, respectively.

Dr. Canudas-de-Wit is an IFAC Fellow and an IEEE Fellow. He was an Associate Editor of the IEEE Transactions on Automatic Control, Automatica, the IEEE Transactions on Control Systems Technology, the IEEE Transactions on Control of Network Systems. He is an Senior Editor of the Asian Journal of Control and the IEEE Transactions on Control of Network Systems. He served as the President of the European Control Association from 2013 to 2015 and as member of the Board of Governors of the IEEE Control System Society from 2011 to 2014. He acted as General Chair of the 58th IEEE Conference on Decision and Control, Nice 2019. He holds the ERC Advanced Grant Scale-FreeBack from 2016 to 2022.
<figure><img src="image_21.png" /></figure>

**Paolo Frasca** (M'13-SM'18) received the Ph.D. degree from Politecnico di Torino, Turin, Italy, in 2009. After Postdoctoral appointments in Rome and in Turin, he has been an Assistant Professor with the University of Twente, Enschede, The Netherlands, from 2013 to 2016. Since October 2016, he has been CNRS Researcher with GIPSA-lab, Grenoble, France.

His research interests cover the theory of network systems with applications to infrastructural and social networks. He has (co)authored more than 40 journal publications and the book *Introduction to Averaging Dynamics Over Networks* (Springer). He has been an Associate Editor for the Editorial Boards of several conferences and journals, including the IEEE Control Systems Letters, the *Asian Journal of Control*, and *Automatica*.

<figure><img src="image_22.png" /></figure>

**Ursula Ebels** received her PhD in Physics in 1995 from the University of Cambridge UK, and held subsequently two postdoctoral positions (Ohio State University, USA and IPCMS in Strasbourg, France). In 2002 she joined SPINTEC laboratory as a CEA research engineer to develop the RF magnetization dynamics group. In 2009 she obtained her habilitation. She has co-authored more than 100 articles and 4 book chapters, and presented her work at numerous international conferences and workshops. She was coordinator of several French national ANR projects

as well as one FP7 project. Her research work covers different aspects of spin torque driven magnetization dynamics of nanostructured magnetic systems in view of microwave applications as well as for defining oscillator based computing schemes.