Synchronization of Spin-Torque Oscillators via Continuation Method

In this article, we study synchronization phenomena of spin-torque oscillators coupled on a ring. Spin-torque oscillators are nanoelectronic devices that promise efficient microwave generation provided they are synchronized in large arrays. Due to their nonlinear and nonisochronous nature, their synchronization properties are difficult to analyze explicitly. To address this challenge, we employ a recently developed continuation method and transform the network of coupled oscillators (each described by an ordinary differential equation) into a single partial differential equation (PDE). We then analyze 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 characterize all possible synchronous solutions and provide necessary and sufficient conditions for stability. For nonidentical oscillators, we derive and solve a differential synchronization condition, which allows us to accurately reconstruct the shape of the equilibrium profiles. All the presented results are derived for the PDE and then validated by numerical simulations of the original network of ordinary differential equations (ODEs).

example of synchronization.The Kuramoto model describes a network of coupled oscillators whose dynamical variables are phase angles.The Kuramoto model and its second-order generalization can be applied to study various applications, such as chemical oscillators [1], 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 topologies of the interconnection network and different parameters, has, thus, been extensively studied, see [5] for a broad literature review, and [2] and [6] for recent comprehensive synchronization conditions.
The Kuramoto model describes only the phase dynamics of oscillators, while assuming that the oscillation amplitude is constant for all oscillators.In many practical applications, however, the relationship between amplitude and phase dynamics cannot be neglected.The oscillators for which the oscillation frequency can vary with the amplitude are called nonisochronous.Many well-known models, such as the van der Pol oscillator [7] or the FitzHugh-Nagumo neuron model [8], belong to this class.
The class of nonisochronous oscillators also includes spintorque oscillators (STOs), nanoelectronic devices that are based on the spin-transfer torque effect discovered by [9] and [10].By this physical effect, a direct electric current that passes through a magnetized layer can become spin polarized, and this spin-polarized current can further transfer angular momentum to another magnetized layer.This transfer induces a torque on the magnetization of the second layer, which can lead to the switching of the magnetization direction.A steady-state magnetization precession can be induced under appropriate conditions on the external field and for current densities larger than a critical value.Due to the fast precession of magnetization, STOs produce microwave output voltage signals.Therefore, large arrays of STOs can theoretically serve as very efficient microwave generators.This potential application explains why the question of synchronization of STOs is very important: synchronous oscillations of many oscillators amplify each other due to constructive interference, whereas asynchronous oscillations exhibit destructive interference and, thus, produce less power.Different magnetic configurations that facilitate synchronization have been studied in [11] and [12], whereas the impact of different device properties have been investigated in [13] and [14].However, analytic studies were mostly limited to the properties of amplitude-frequency coupling, synchronization to an external oscillating force and phase-locking effects, as in [15] and [16].
In this article, we study the synchronization of an array of STOs with a ring interaction topology.In order to completely describe synchronization on a ring, we are going to reconstruct the different synchronization profiles that are possible in the system, and to derive necessary and sufficient conditions on the system parameters that can be checked to ensure that these synchronized solutions are stable.These conditions could then offer guidance for the realization and deployment of large arrays of synchronously operating STOs.
We approach this challenging problem by utilizing a recently developed continuation method, which transforms a network of coupled ODEs into a single partial differential equation (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.The continuation method provides continuous approximations that are guaranteed to become arbitrarily accurate as the order of the PDE grows [17], in contrast to other "mean-field" methods that require the number of ODEs to be large for the approximation to be meaningful.
In order to demonstrate the applicability of the method to the problem of network synchronization, in Section II, we start with an illustrative but significant example, which is the analysis of Kuramoto model synchronization on a ring.We then proceed to introduce the STO model and its continuation in Sections III and IV, respectively, again in the case of ring interconnection.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 (Section V), our main results provide a characterization of all synchronized solutions and the explicit necessary and sufficient conditions on the system parameters that 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 nonidentical 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 derive a differential equation serving as a condition for synchronization and solve it to obtain an implicit function, which describes the synchronous solution's profile.In this way, we are successfully able to analytically reconstruct solutions arising in the inhomogeneous system, as we demonstrate using numerical simulations.

II. 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 1-D ring.By deriving PDE representation of Kuramoto system, we show that the PDE model can be more convenient for analysis, in the same way as continuous dynamical systems can be more tractable than the discrete ones.Networks of Kuramoto oscillators synchronize when their coupling is stronger than a certain threshold.This has been known for a long time for complete interaction topologies, whereas the general problem of computing this synchronization threshold for different network topologies and frequency distributions was recently solved in [2] and [6] with the help of graph theory.However, the method presented in these papers is not straightforward to extend to other types of oscillators.On the contrary, the idea based on continuation, which is presented in this section, allows for finding a synchronization condition in a very natural way.Moreover, this method will be extended in the next sections to a more general class of nonisochronous 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 analyzing a system of Kuramoto oscillators where φ i is a phase angle of ith oscillator, ω i is its natural frequency, and F > 0 is a coupling strength.Each oscillator is coupled with its two neighbors, forming a closed ring.We assume that there are n oscillators and that each oscillator has a position on the ring defined by x i ∈ [0, 2π), with x i+1 − x i = Δx and x 1 − x n + 2π = Δx, meaning that the oscillators are spaced equally on the ring.Using these positions, we can naturally define a frequency function ω(x i ) = ω i and a state function 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 which leads to the system Now we can calculate the Taylor expansion of the function φ(x i ) at the point The same expansion can be used for φ(x i−1 ), so that the argument of sine in (2) is In general, these Taylor series are infinite, but we want to obtain a low-order approximation; therefore, we assume that Note that this formula 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 Moreover, by applying the same procedure to function s(x), we can write ∂φ ∂t = ω(x) + F Δx ∂s ∂x and by using the approximation of s(x) in (3), we obtain the PDE approximation 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 ω(x) = 1 + x sin(2x) for x ∈ [0, 2π) as a working example (in general any function can be used, but for the future analysis we choose one that can be integrated in closed form).The coupling strength is set to F = 4.We numerically simulate the approximated PDE (4) on a grid with 500 points.The results of the simulation are shown in Fig. 1: while the ODE solution splits into several clusters, the PDE solution continuously connects these clusters (see Fig. 1, left) and remains along time a rather accurate approximation of the former (see 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 φ * of (4) in the case of complete synchronization.In such a case, there exists ω = 1 2π 2π 0 ω(x)dx such that all oscillators share the same frequency Therefore, the equilibrium solution should satisfy Let us integrate this equation from x 0 to x 1 , where both are chosen arbitrary where Ω(x) is some primitive function of ω − ω(x).Rearranging terms, we obtain where C is a 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 Ω(x) written in the form If such Ω(x) exists, φ * 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 Ω(x) such that (7) is possible, meaning that the sine value lies in the interval with C being some constant.Recalling that Ω(x) is a primitive function of ω − ω(x) and that in general is defined up to a constant, this is equivalent to the condition for any Ω(x).Recovering the synchronization threshold F * only requires to replace the inequality with the equality sign The synchronization threshold (8) provides a condition on the existence of equilibrium solutions.One can interpret F * as a measure of the cumulative deviation of the natural frequencies of the oscillators from the mean frequency: the larger this deviation is, the stronger coupling is required in order to keep oscillators together.Consistently, this formula can be seen as a continuous counterpart of the critical coupling max i ω i − min i ω i that has been established for classical Kuramoto networks; see, e.g., [5,Eq. (46)].Furthermore, it is possible to prove that for all F > F * , there will exist a stable equilibrium solution.
Proof: Without loss of generality we assume that C = 0 (because the primitive function Ω(x) is defined up to a constant).Note that Ω(x) ∈ [−F * Δx, F * Δx].Then, the equilibrium solution can be recovered from (7) (up to a constant) as Now let us assume that the equilibrium solution is slightly perturbed: φ = φ * + φ.Then, by ( 4) Rewriting sine, we get where the fact that φ is a small perturbation was used.Now (5) cancels the natural frequencies, therefore we arrive at This is a standard linear diffusion equation with the diffusion coefficient cos(Δx ∂φ * ∂x ).For stability, it remains to prove that this coefficient is always positive.Indeed, ) and the linearized system is locally asymptotically stable.
Remark 1 (Necessary synchronization condition): Theorem 1 shows that the condition F > F * is a sufficient condition for the existence of a stable equilibrium solution.In fact, it is possible to show that the condition F F * is a necessary one: indeed, if an equilibrium solution exists, then it should satisfy (9) in order to render (5) true.But this is possible only if the arcsine function in ( 9) is well defined for all x, which means that F F * .The case of F = F * , thus, leads to the existence of an equilibrium solution, even though it is not possible to prove its asymptotic stability in general.To validate this analysis, we use the same parameters of the simulation as for Fig. 1: the length of the ring is L = 2π, the number of ODE nodes n = 50, and the natural frequency ω(x) = 1 + x sin(2x) (which is an integrable function, thus can be analytically treated).By definition of the positions of nodes, Δx = L/n = 2π/50.Further, thus ω = 1/2.The primitive function Ω(x) can be taken as Substituting these values in (8) gives This value is the smallest F for which the equilibrium solution exists.To verify the result (13) for the original system (1), we simulated the latter for F ∈ [0, 20] and calculated ω = max φi − min φi , which we call desynchronization frequency.In the case of complete synchronization, ω should be zero.Indeed, Fig. 2 shows that ω is zero for F > 18.9, and it increases when F becomes smaller, thus behaving in accordance with the derived value of F * in (13).We would like to stress that, having chosen ω(x) to be analytically integrable, the PDE approximation gives benefits in terms of computation complexity: we have been able to simply integrate it to find Ω(x) and obtain an analytic answer to the synchronization threshold problem, thus requiring only O(1) computations.Instead, solving this problem directly for the ODE system requires at least O(N ) computations where N is the number of oscillators.

III. MODEL FOR THE RING OF STOS
It was shown in the previous section that synchronization analysis can be performed for Kuramoto oscillators using the PDE representation.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 nonisochronous oscillators.As a particular example within this class, we focus on STOs.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 3. Schematic representation of a possible geometry of STO (one single element of the ring interconnection).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 P. The magnetization M of the "free" magnetic layer then oscillates under the effect of polarized current and the effective magnetic field H. Fig. 4. Close view on the dynamics of the magnetization M of the "free" layer, governed by (14).Damping and current-induced spintransfer torque compensate each other, stabilizing steady oscillations caused by precession around the magnetic field H.
A typical STO consists of two ferromagnetic layers, a thick one called "fixed" and a thin one called "free," see 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, which manifests itself as a torque acting on the free layer magnetization.This effect creates a precession, depicted in Fig. 4. Denote the magnetization of the "free" magnetic layer by vector M, the magnetization of the "fixed" magnetic layer by vector P, and the effective magnetic field by vector H; then the dynamics of the "free" layer magnetization is governed by the Landau-Lifshitz-Gilbert equation spin transfer (14) where the parameters γ, α, and σ depend on the system's geometry and materials, and I is a current density that is applied to the system.For a more complete review of the spin-transfer torque effect and of the physics of STOs, we refer the readers to [15], [16], and [20].Equation ( 14) can be simplified for analysis.The magnetization vector M oscillates around the effective magnetic field vector H.Let us project M on a plane orthogonal to 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 ( 14) of an STO can be modeled through where p = |c| 2 represents a squared amplitude of oscillations, ω is a linear frequency, N is a nonlinear frequency coefficient, Γ G is a linear damping, Q is a nonlinear damping coefficient, and I and σ are the same as in (14).Model ( 15) is nonlinear since the oscillations' frequency depends on the amplitude through the frequency gain N .In the case of STOs, this amplitude-related frequency shift happens to be very strong, thus these oscillations cannot be described by simpler linear models.
If σI Γ G , then the origin c = 0 is a stable equilibrium point, whereas oscillations will occur if σI > Γ G .Assuming the latter to be true, we define the linear part of the dissipative terms Γ = σI − Γ G > 0 and the nonlinear gain of the dissipative terms S = Γ G Q + σI, thus rewriting system (15) as System (16) will oscillate with amplitude √ p = Γ/S and with frequency φ = ω + N Γ/S, where φ is a phase of an oscillator.
For the amplitude of oscillations to be well defined, we also require S > 0.
In physical systems, the phase of an oscillator usually changes much faster than its amplitude.Therefore, model ( 16) is often studied in amplitude-phase representation c = √ pe iφ , where √ p is the amplitude of the oscillations and φ is their phase.Instead of writing two separate equations for these variables, we will write model ( 16) in logarithmic representation.Define z = ln c.Then, the real part of z will represent the amplitude, namely exp{2 Re z} = p.Let us denote r := Re z = 1 2 ln p.The imaginary part of z is the phase of the oscillator, φ := Im z, thus this transformation allows to track phase information directly.Since dc = c • dz, the model ( 16) now becomes 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 neighbors.As in the previous section, let n denote the number of oscillators and let x i ∈ [0, 2π) be a position on a ring of the ith oscillator.The distance between oscillators is x i+1 − x i = Δx and x 1 − x n + 2π = Δ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 Here, F i is a (possibly complex) coupling constant, with an amplitude representing the coupling strength and a phase representing the coupling phase.Using the logarithmic representation, the model ( 18) reads as

IV. CONTINUATION AND SYNCHRONIZATION CONDITION
It is now possible to perform the continuation for the coupled system in the same way it was done for Kuramoto oscillators in Section II.By following the method illustrated in [17], the continuation of ( 19) is performed in several steps: 1) By combining these continuation formulas, we finally get Thus, system (19) can be transformed into the PDE model  (20) and conditions for their existence and stability.A synchronized solution is a solution to (20) such that ∂z/∂t = iω, where ω is a synchronization frequency.Thus, we are interested in the question when such a solution z = z(x) exists for some ω.The condition for synchronization is or in terms of r(x) and φ(x) Note that we divided the equation by F before splitting real and imaginary parts so that the hyperbolic functions take the simplest form.
The exponential term e 2r in (24) can be removed by combining two equations together.Using amplitude-phase notation, we can introduce where Therefore, the synchronization condition is equivalent to (25) combined with one of the equations in (24) to determine the connection between r and φ.Equilibrium solutions to (25) represent synchronized states of the original PDE (20).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 must exist a solution where all oscillators share the same amplitude r and the same phase φ.However, it will appear that, depending on the number of oscillators and their parameters, there can be more solutions, and that their stability properties are not trivial.
In order to search for equilibrium solutions, let us assume that in the synchronized state the amplitudes of oscillators r are slowly varying in space, namely, ∂r/∂x ≈ 0. We will further show that in the identical oscillators case, the amplitude is indeed constant in space.Due to the assumption ∂r/∂x ≈ 0, we can use sinh(Δx ∂r/∂x) ≈ 0 and cosh(Δx ∂r/∂x) ≈ 1.With these simplifications, (25) depends only on φ(x) and, thus, can be solved independently If the parameters of oscillators were nonidentical, (27) 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 constants A and B, (27) is separable.We can notice that it depends only on the derivative of φ(x), not on the phase itself.Define θ = Δx ∂φ/∂x.The physical meaning of θ is the difference in phases between two consecutive oscillators.With this definition, (27) becomes The general solution to (28) is described in Section VI.For now, let us note that by the structure of (28) any nonconstant continuous solution θ(x) would be monotonic.Further, apart from being a solution to (27), synchronization means that the solution φ(x) is a continuous angle, thus φ(x + 2π) − φ(x) = 2πk for some k ∈ Z.This implies two conditions that θ should satisfy In the case of identical oscillators with constants A and B, it appears that the only possible solution to ( 27) is a constant one.Indeed, nonconstant solutions should be monotone with respect to coordinate; however, the first condition in (29) requires θ to be periodic, which is not possible if θ is not constant.Thus, looking at (27), we see that for constants A and B, there is a simple solution cos θ = −B/(2 A), which means that all possible synchronized solutions for (27) are given by

A. Equilibrium Points
Recall that θ = Δx ∂φ ∂x .Since φ(x) is a phase, it is defined up to a constant.Assuming φ(x) = 0 at x = 0, using (30) and the definitions of A and B, the solution for φ(x) is a linear function Note that ω is a synchronization frequency and is still unknown in this equation.Position x is itself defined on a ring, thus x ∈ [0, 2π).Moreover, since the equilibrium solution is a periodic function, φ(2π) should also be a multiple of 2π.We can define k ∈ Z + such that φ(2π) = 2πk.Since k determines the number of phase turns along the ring, we shall call it the winding number.Therefore, the solution can exist for any ω such that The case k = 0 corresponds to an in-phase synchronized system, meaning that the phases of all oscillators coincide, whereas 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 neighbors is Note also that the system is symmetric for simultaneous substitution k → −k and x → −x, thus phases can turn both clockwise and counter-clockwise along the ring.The principal branch of arccos has a range of values [0, π], therefore k should satisfy |k| π/Δx (other solutions will just copy the ones included in this range due to periodicity).Since Δx is defined as the distance between two oscillators and is assumed to be constant, Δx = 2π/n, where n is the number of oscillators in the system.Thus, |k| n/2, with |k| = n/2 corresponding to the case when two neighbor oscillators are in antiphase.
The synchronization frequency is, thus, given by for k ∈ {0, . .., n/2}.In particular, depending on the sign of 2f • sin(β + γ)/ cos γ, the in-phase synchronized state is either the fastest or the slowest one.

B. Stability Analysis
Assume that the equilibrium solution is given by (31) with the frequency (32) for k ∈ {0, . .., n/2}.We want to study for which of these k the solution is stable.
Let z * (x, t) = r * + iφ * (x) + iωt be an equilibrium solution for (20).Thus, φ * (x) is defined by (31) for a chosen k, ω is the frequency of the synchronized solution (32), and the constant r * can be found from (23) by taking its real part Note that the exponential should be positive to be well defined, therefore we require Γ + 2f cos β cos(kΔx) > 0.
Now let us define a deviation from the equilibrium solution z(x, t) = z(x, t) − z * (x, t).It is governed by a difference of (20) for z(x, t) and for z * (x, t), taking into account (23).Assuming that z(x, t) is small, the linearization of (20) around z * (x, t) is given by Using which can be substituted in (34), resulting in Then, we can separate (35) into real and imaginary parts z = r + i φ and use F = fe iβ to obtain System (36) 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 exponential.Therefore, the stability of (36) can be checked by substituting exponential basis functions r = r 0 e λt e imx , φ = φ 0 e λt e imx (37) for some λ ∈ C and m ∈ Z, since basis should be periodic in x along the ring.For asymptotic stability, there should exist no solution of (36) with Re λ > 0. Substituting (37) in (36), one gets For notational convenience, define and note that as m → ∞, the first terms become dominating.
With the help of these functions and with S = 2Se 2r * > 0 and N = 2Ne 2r * , (38) becomes Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
thus we are interested in the eigenvalues of the matrix in (39).It is trivial to show that they are given by Taking m = 0, one of the eigenvalues becomes zero, corresponding to the fact that the phase is defined up to a constant, and the other eigenvalue is − S. Further, assume that m = 0 and, thus, P, Q = 0. Condition for stability Re λ < 0 translates as Re(2P + S) > Re (2P + S) 2 − 4P (P + S) − 4Q(Q + N ).
(40) In particular, as m → ∞, Re P should be positive.Now for simplicity, define Using the complex identity which then being inserted in (44) results in Finally, defining Ḡ = 2e 2r * G and using S = Ḡ cos γ and N = Ḡ sin γ, we get Recall that Re P > 0 by (40) and, substituting P and h m from (45), we also imply that cos(kΔx) cos β > 0. Thus, we just proved a theorem, which can be seen as the main result of this article.
Due to the dependence of ( 46) 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 directly used for particular systems to establish a guaranteed stability of synchronized solutions.
Corollary 1 (k = 0): Necessary and sufficient conditions for in-phase synchronization are given by Proof: The in-phase equilibrium solution satisfies k = 0; thus, by (45), d m = 0 and h m = f Δx 2 m 2 > 0. From the second condition of Theorem 2, we recover cos β > 0. Finally, (46) with d m = 0 requires right-hand terms to be greater than zero, which is just h m (h m + Ḡ cos(γ + β)) > 0. Since this is always true as h m → ∞ with m → ±∞, it is enough to satisfy this inequality for m = ±1, leading to cos(γ + β) > −f Δx 2 / Ḡ.
Notice that the conditions in Corollary 1, as well as in all other corollaries in the following, immediately ensure the existence of the exponential representation of the amplitude of oscillations defined in (33).
Corollary 2 (k = n/2): Necessary and sufficient conditions for antiphase synchronization are Proof: The proof follows the same steps as the previous one, switching the sign of h m .
Corollary 3 (k ∈ {1, . . ., n/2 − 1}): Sufficient conditions for synchronization with sin(kΔx) = 0 are given by cos(kΔx) cos β > 0 (47) together with Condition ( 47) is also a necessary condition for stability.Proof: First, condition (47) repeats the second condition of Theorem 2. Further, since sin(kΔx) = 0, d m is nonzero.Divide Inserting the definitions of h m and d m from (45), we see that the right-hand side of ( 49) is strictly increasing with m 2 , therefore it can be simplified by setting m 2 = 1 as in the worst case, thus obtaining the right-hand side of (48).Now, to find sufficient conditions for the satisfaction of (49), let us bound the left-hand side from above.For this, we will use the following Lemma, whose proof can be found in Appendix A.
Remark 3 (Winding number and number of oscillators): Assume that cos β > 0 such that the in-phase solution is stable.Then, the second condition of Theorem 2 requires that cos(kΔx) > 0. Thus, for the stability, kΔx should be smaller than π/2.This means k < n/4, where n is the number of oscillators.In particular, the phase difference between two neighboring oscillators should be smaller than π/2.This also means that in order to observe a state with k = 1, one needs at least five coupled oscillators, and to observe higher order states, one needs at least nine oscillators.As an example, all possible states in the system with ten oscillators are shown in Fig. 5.

C. Numerical Simulation
We now want to test the predictions about synchronization of coupled STOs that can be drawn from the PDE results in the previous section.To this aim, we perform a numerical simulation of n = 50 coupled STOs.The comparison between analytic predictions and numerical simulations is shown in Fig. 6.The simulation parameters are chosen to be physically sound according to [15], namely, we set ω = 6.55 • 2π, N = −3.82• 2π, Γ G = 0.375 • 2π (all these measured in radians per nanosecond), Q = −0.24,and σ = 5.48 • 10 −13 • 2π rad/nanosec/A * m 2 for (15).In this case, the critical current density, which is required to start oscillations, is I c = Γ G /σ = 684.3• 10 9 A/m 2 .In our experiments, we use a larger current density I = 1.5I c that enables  First, we notice that the condition cos β < 0 in Corollary 2 is unlikely to be satisfied in real systems.We, thus, focus the simulation on testing the conditions in Corollaries 1 and 3, which are also quite restrictive due to the negative value of γ in our Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.example.We can check which stable synchronized solutions are admitted by the coupled system depending on different coupling parameter F .We take different couplings F = fe iβ with f changing from 0 to 10 and β 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 Fig. 6(a).Further, we compare them with experimental results by simulating the original ODE system (18).We initialize all oscillators in this system using an amplitude √ p i = Γ/S and a phase φ i = ikΔx for the ith 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 ns (corresponding roughly to 15 000 periods of oscillation for f = 0.75).When the simulation ends, we check if the system remained stable or diverged from the corresponding equilibrium solution.The system is considered to be stable if the amplitude of the deviation from the equilibrium solution in the end of the simulation is not larger than the deviation at the start.Notice that due to the nature of numerical simulations, this stability check can be imprecise near the regime boundaries.Fig. 6(c) depicts the highest stable regimes that are observed.Comparing it with Fig. 6(a), we see that the analytic prediction almost perfectly reconstructs the experimental diagram, with small deviations that should probably be attributed to inaccuracies in the aforementioned numerical stability check.Finally, we compare the synchronization frequency ω predicted by (32) 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 ω ≈ Im( ċ/c).Then, we average this frequency over all agents and over the last 1000 ns.The measured synchronization frequency for f = 0.75 and for β ∈ [−0.1, 0.2] is depicted in Fig. 6(d).It is clear that for the higher regimes for k = 1 and k = 2, stable solutions exist only for sufficiently high values of β.Comparing measured frequency with the analytical prediction by (32) in Fig. 6(b), 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/ns 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 (20), while the simulation was performed for the original ODE system (18).

VI. NONIDENTICAL 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 the requirement on homogeneity but we keep the assumption that ∂r/∂x ≈ 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 (27), and for constant parameters A and B, it is represented by a separable (28).Let us define J = B/A.By performing the integration of (28), it is trivial to show that all solutions to (28) are given as follows: where C is an integration constant, and g(θ, J) is a special function parametrized by J, which is given by (53) It is interesting to note that there is a complex relation between arctangent and logarithm functions which means that the first two cases in (53) are essentially the same.In fact, (53) defines a piecewise continuous function with at most two singularities with respect to J, see Fig. 7.In particular, (52) means that the constant solution (30) is captured by the singularity approaching J = −2 from the left.Note further that g(θ, J) is an odd function with respect to θ.

A. Case of Two Types of Oscillators
Using the solution (52), it becomes possible to analyze 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 ω 1 , N 1 , Γ 1 , 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π/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, whereas if K = n/2, the types of oscillators alternate.We can define a set of switching points as y j for j ∈ {0, . .., 2K − 1}, with y 0 = 0 and y j = j/2 • 2π/K Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.for even j.Finally, for odd j, we require y j − y j−1 = const, thus the proportion of types is preserved.Oscillators placed in [0, y 1 ) ∪ [y 2 , y 3 ) ∪ . . .are of the first type, and oscillators placed in [y 1 , y 2 ) ∪ [y 3 , y 4 ) ∪ . . .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 ω should be common for both types; therefore, by definition of B in (26), we can write J 1 = J1 + τ 1 ω and J 2 = J2 + τ 2 ω, where (55) with J2 and τ 2 being defined in a similar way.
We are now interested in particular solutions θ(x) to (27).By (29), θ should be periodic.Since intervals of types of oscillators are equal, symmetry leads to the fact that θ should be periodic with period being equal to two intervals of different types of oscillators, namely, θ(y 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 θ(x) behaves continuously and monotonically in the first type domain and is constant with discontinuity in the interior in the second type domain.Moreover, solution endpoints are symmetric about zero, namely, θ(y 0 ) = −θ(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 (53) and, thus, requires special treatment.Further, in this section, we will stick to the class of solutions in agreement with Observation 1.
Defining θ * = θ(0) and assuming θ(y 1 ) = −θ * by Observation 1, we can compute (52) in points x = y 0 = 0 and x = y 1 for the first type of oscillators and subtract one from another, obtaining where we used the fact that the function g(θ, J) is odd with respect to θ. Substituting J 1 as in (55), we get a condition that should be satisfied for the first type of oscillators which have two unknowns: θ * and ω.The second condition comes from the assumption that for the second type domain, the solution is constant, and thus, it is determined by (30).Using it for the second type domain, we get Note that both θ * and −θ * are solutions to (30), which is consistent with Observation 1. Now, substituting J 2 by (55) in (57) and then substituting the result in (56), we obtain an equation with a single unknown ω This equation can be solved for ω using numerical methods, such as Newton method, for example.Once ω is known, we can find J 1 and J 2 by (55) and then compute θ * by (57).The full solution on the first domain is then reconstructed by (52).
To determine the shape of solution θ(x), it remains only to find an exact position denoted by y * ∈ (y 1 , y 2 ), where a discontinuous jump from θ * to −θ * happens in the second type domain.This position can be obtained if one recalls that θ = Δx ∂φ/∂x and, thus, integral of θ should have fixed value by (29) for some k ∈ Z.In particular, due to the periodic nature of the problem with K periods, we have Since, on the first type domain, θ(x) is symmetric, its contribution to the integral is zero.Further, θ(x) = θ * on x ∈ [y 1 , y * ) and θ(x) = −θ * on x ∈ (y * , y 2 ], therefore (59) is just Thus, the solution's shape θ(x) is fully reconstructed.Remark 4: 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 5: 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 (52).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(θ, J), thus they cannot be found reliably by numerical methods without a suitable problem reformulation.

B. Numerical Simulation
To demonstrate how solutions to the synchronization condition (27) found by ( 57)-( 60) approximate synchronized solutions of the original system (18), we performed numerical simulations of (18) with n = 500 oscillators being split into two types as it was described earlier in this section.The parameters of the first type of oscillators were taken the same as in Section V-C, and for the second type, slight deviations in the parameters were added.The 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 a y 1 /y 2 proportion of the period of the length y 2 , and each group of oscillators of the second type occupies a (y 2 − y 1 )/y 2 proportion.The numerical simulation was initialized in the same way as in Section V-C, with k denoting the initial phase shifts between consecutive oscillators such that the phase makes k turns along the ring.We performed two simulations, with two rather different sets of parameters in order to demonstrate that our analytic results can be applied to different kinds of STO systems.
1) In the first simulation, we altered the damping parameter Γ for the second type of oscillators, so that Γ 2 = Γ 1 • 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 in such a way that the phase makes k = 3 turns along the ring.2) In the second simulation, we changed the frequency gain parameter N for the second type of oscillators such that N 2 = N 1 • 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 in such a way that the phase makes k = −2 turns along the ring, rotating in the opposite direction.The results of the simulation are presented in Fig. 9.The simulation was performed for 2000 ns, and then, the phase differences between consecutive oscillators were computed.The result was then compared with analytic predictions by (57)-(60).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 an assumption of small magnitude variation.Also, numerical simulations have shown that synchronized solutions are very fragile, meaning that small deviations in the parameters result in very large phase differences between consecutive oscillators, although the system still remains stable and the PDE model predictions are still correct.

VII. CONCLUSION
Understanding the synchronization of STOs has a big practical importance since synchronous oscillations produce energy much more efficiently.Therefore, it is very important to realize when synchronized solutions do exist and to which extent they tolerate deviations in the parameters (which originate from imperfections in the fabrication of the oscillators).In this article, we have shown how the continuation method can help in the analysis of this problem, and we derived results that can 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 nonidentical oscillators, we analyzed 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.Despite these initial insights, there remain many questions that could be investigated in details regarding system (18), its PDE approximation (20), and the synchronization condition (23).First, Corollary 3 for the general winding number k in the case of identical oscillators gives only sufficient conditions on stability: likely, tighter sufficient conditions could be derived from Theorem 2. Second, the important case of nonidentical oscillators, discussed in Section VI, covers only a search for specific equilibrium solutions.Beyond this initial insight, more general equilibrium analysis and investigation of stability conditions could be performed.
Furthermore, this work has been devoted to the analysis of oscillators that are arranged on a ring, per model (18), as a prototypical case of study.The study of more general topologies, especially 2-D, would be of great importance for practical applications: we are confident that continuation methods can be useful for these extensions.Indeed, continuation methods do extend to geometries that are nonperiodic, irregular, and in higher dimensions: 2-D and 3-D geometries.
Finally, we would like to comment on the scope and validity of application of our continuation method.In that respect, it should be noted that, even though continuation is most natural for large arrays, the precision of the continuous approximation is guaranteed when the order of approximation goes to infinity, irrespective of the number of oscillators [17,Th. 1 and 2].This fact suggests that the method can be effective to approximate small arrays of oscillators as well.

APPENDIX A PROOF OF THE LEMMA 1
The function f (x) is defined for x ∈ [0, +∞), thus its supremum is achieved either at x = 0 and x = +∞ or at f (x) = 0.If V 0 and μ 0, then the function is nonpositive with asymptotic value f (+∞) = 0, thus we use 0 as a bound in this case.Let us now find its extremum thus it is achieved at x extr = U − 2 V /μ.Substituting it back in (61), we obtain Finally, we notice that the extremum (61) is indeed maximum only if μ > 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 (63)

Fig. 1 .
Fig. 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 meansquare absolute divergence between solutions.

Fig. 5 .
Fig. 5. Six possible equilibrium solutions (31) for the ring of ten STOs.Assuming cos β > 0, the first three are stable and the second three are unstable.

Fig. 6 .
Fig. 6.Synchronized solutions for a system of n = 50 coupled identical STOs.Top row: analytic results for PDE (20).(a) Diagram of possible regimes by Corollaries 1 and 3. Color code denotes the highest guaranteed existing regime, and chaotic means that no stable solution exists.(b) Synchronization frequency ω by (32) for f = 0.75 and different k depending on β.Bottom row: numerical simulation of (18).(c) Diagram of numerically established regimes.(d) Numerically measured synchronization frequency.

Fig. 8 .
Fig. 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.
, and F are functions of space, determined by approximating sampled values Γ i , ω i , S i , N i , and F i at points x i .The analysis of this PDE system equation will be our focus for the rest of this article.
Separating (20) into a system of two equations for r = Re z and φ = Im z, one gets