Describing synchronization and topological excitations in arrays of magnetic spin torque oscillators through the Kuramoto model

The collective dynamics in populations of magnetic spin torque oscillators (STO) is an intensely studied topic in modern magnetism. Here, we show that arrays of STO coupled via dipolar fields can be modeled using a variant of the Kuramoto model, a well-known mathematical model in non-linear dynamics. By investigating the collective dynamics in arrays of STO we find that the synchronization in such systems is a finite size effect and show that the critical coupling—for a complete synchronized state—scales with the number of oscillators. Using realistic values of the dipolar coupling strength between STO we show that this imposes an upper limit for the maximum number of oscillators that can be synchronized. Further, we show that the lack of long range order is associated with the formation of topological defects in the phase field similar to the two-dimensional XY model of ferromagnetism. Our results shed new light on the synchronization of STO, where controlling the mutual synchronization of several oscillators is considered crucial for applications.

Scientific RepoRts | 6:32528 | DOI: 10.1038/srep32528 only a few oscillators has been demonstrated 30,32 . Theoretically, the magnetization dynamics of STO is modeled with the Landau-Lifshitz-Gilbert-Slonzewski (LLGS) equation 41,42 , but large number of STO lead to challenging computations caused by the non-local dipolar fields. It is important to consider that in these non-linear systems "more is different", and that the collective behavior can not be derived simply from the behavior of its individual elements. Thus, a theoretical framework capable to capture the essential dynamics would be ideal to explore those systems.
Here, we show that two-dimensional arrays of STO coupled via dipolar fields can be modeled by a variant of the Kuramoto model. We begin with describing two coupled STO with the Thiele equation 43 and show that for small-amplitude oscillations the system can be described as a simple phase oscillator model. Next, we model the interactions for the case of a two-dimensional array of oscillators based on the dipolar coupling and obtain a modified Kuramoto model. Finally we compare the results from our model to the micromagnetic solution of the LLGS equation.
We find that the synchronization in two-dimensional arrays of dipolar coupled STO is purely a finite size effect and the critical coupling strength for obtaining a globally synchronized state scales with the number of oscillators N as λ crit ∝ log(N). Using realistic values of the dipolar coupling strength between STO we show that this imposes an upper limit for the maximum number of STO that can be synchronized. Further, we study the synchronization transition between the initial formation of locally synchronized clusters and the globally synchronized and phase coherent state and correlate it with a transition in the local order of the system. We also observe the emergence of topological defects and the formation of patterns in the phase field similar to the two-dimensional XY-model of magnetism-suggesting a connection between arrays of STO, systems described by a 2d Kuramoto model and the 2d XY model of statistical mechanics.

Results
From the Thiele equation to the Kuramoto model. We are considering STO whose free layer ground state configuration is a magnetic vortex. The vortex state is characterized by in-plane curling magnetization, and a small (~10 nm) region of the vortex core with out-of-plane magnetization 44 . The gyrotropic motion of the vortex core is driven by the injection of a DC spin polarized current through the STO stack, and can be described by a gyration radius r and phase θ, as illustrated in Fig. 1a. We first model the interaction of two vortices with the Thiele equation with an extra term that accounts for the vortex interaction. These equations describe the vortices motion given by their coordinates X 1,2 in their self induced gyrotropic mode, and include the spin-transfer-torque (STT) as well as a coupling term 43,45 .  Here, G is the gyroconstant, k(X 1,2 ) the confining force, D 1,2 the damping coefficient and F STT the STT. The interaction between the neighboring STO illustrated in Fig. 1b is summarized by a dipolar coupling term given by F int = − μ(d)X 2,1 , where μ(d) describes the interaction strength as a function of the separation d between the two STO. Assuming a small difference in the nominal frequencies of two coupled STO described by Eq. (1), one can linearize the set of equations following the approach by Belanovsky et al. 27 , showing that the dynamics of the phase difference between the STO can be described by Adler's equation 46 . Following these approximations, the set of equations reduce to that of two coupled phase oscillators θ 1 and θ 2 : (see 'Supplementary information' for details).
To check the validity of the approximations, one can compare the results obtained using a simplified phase oscillator model to a numerical solution of Eq. (1) as well as a micromagnetic solution of the full system using the LLGS equation. This was done by Belanovsky et al. 27 , where they found that for a small difference in nominal frequencies, in their case given by a difference in STO disc diameter Δ D/D 0 ≤ 5%, the synchronization can be qualitatively described using the simplified model. Assuming the error in state-of-the-art fabrication processes is below this limit, the simplified equations are a valid description of the system.
The functional form of Eqs (2-3) is the same as that of the well known Kuramoto model 8,9 , which is a generalization for the case of an ensemble of weakly coupled phase oscillators. Considering the interaction between several STO, we obtain a Kuramoto model where the single oscillator state is described through the dynamic equation of its phase θ i due to the interaction with its surrounding oscillators θ j : The coupling term is here generalized to include the interaction between several oscillators, determined by the interaction strength λ ij between oscillators θ i and θ j . This determines the nature of the interaction, ranging from a global all-to-all coupling where λ ij = λ for all oscillators, to a local interaction where λ ij = 0 for all but the nearest neighbors. Here, we are considering the intermediate case of a non-local coupling to mimic the dipolar interaction between neighboring STO. Starting from a macrodipole approximation for the dipolar energy between two magnetic dipoles μ 1 and μ 2 , the average interaction strength is found to decay as µ 47 , where d ij is the distance between oscillators θ i and θ j . We thus set the coupling strength to where we include interactions within a coupling radius R, indicated by the blue circle in Fig. 1c. The network model for the STO array is implemented with bi-periodic boundary conditions and the time evolution of the oscillator phases given by Eq. (4) is solved numerically. A small random disorder in the oscillator eigenfrequencies is included by setting ω i = ω 0 ± δω i . Here, ω 0 = 1 GHz and δω i represents a uniformly distributed random disorder where δω i /ω 0 ≤ 2.5%. The interaction strength is determined by the STO spacing and size, as well as the magnetic material properties 31 . Here, the interaction strength has been varied in the range λ = 1-20 MHz, and is in the same range as the interaction strength extracted from micromagnetic simulations for similar STO 31 .
In order to evaluate the Kuramoto model as a valid description for arrays of STO, we compare it to a micromagnetic solution of the complete system, accounting for all dynamic dipolar terms (see 'Methods' section). To compare the Kuramoto model and the micromagnetic solution, we define a suitable order parameter to distinguish disordered and synchronized states. The phase of the individual oscillators θ i is used to define the order parameter ρ, describing the phase coherence in a system of N oscillators: The case ρ = 0 corresponds to the maximally disordered state, whereas ρ = 1 represents the state where all oscillators are perfectly synchronized and phase coherent. In addition to the global order parameter ρ, we define a local correlation function β: i j i j where the brackets indicate a summation over neighboring oscillators, and n is the number of neighbors. β i is a measure of the phase correlation of oscillator θ i and its neighbors, indicated by the blue and red oscillators in Fig. 1c respectively. If oscillator θ i is located within a synchronized cluster, β i → 1. Calculating β thus allows for investigating the formation of locally synchronized clusters and the emergence of patterns of synchronized states, which can not be obtained simply from the global order parameter ρ.

Kuramoto model versus micromagnetic simulations.
We now compare the Kuramoto model given by Eq. (4) and the full micromagnetic solution of the LLGS equation (see 'Methods' section for details). Starting from a disordered initial state, we investigate the synchronization dynamics by calculating the time evolution of the phase distribution θ i and local correlation β i . As an example, we show in Fig. 2 snapshots of θ i and β i for the Kuramoto model and the micromagnetic solution for a system of 45 × 45 oscillators at times t 1 and t 2 > t 1 , with random initial phases. At time t 1 , one notices the initial formation of small locally synchronized clusters, as seen through both the phase distribution and the bright areas in the correlation maps, where β → 1. As time progress to t 2 , these clusters grow in size and merge with neighboring clusters. Comparing the two models, we find that they both show the same behavior. For weak interaction strengths, the system tends to be in a disordered state with no correlation between neighboring oscillators. By increasing the interaction strength above a certain threshold, synchronized clusters begin to form. The oscillators within each cluster are synchronized, but might not be phase coherent with other clusters. This can be seen in the phase maps in Fig. 2, where the individual clusters have different phases. As time progresses, the transition from a disordered to a synchronized state is governed by the growth and merging of neighboring clusters, reaching a globally synchronized and phase coherent state for sufficiently strong interactions.
Depending on the interaction strength, the system ends up in either a disordered, partially synchronized or globally synchronized state. Controlling the interaction strength is thus the key parameter to determine the system behavior. The interaction strength needed to obtain synchronization will depend on the differences in the nominal frequencies of the oscillators 31 . Another important consideration, is whether the critical interaction strength also depends on the number of oscillators. Lee et al. 48 have studied the synchronization in a 2d Kuramoto model with a nearest neighbor interaction. They showed that the transition to a synchronized state depends strongly on system size, and that the critical coupling strength needed to synchronize scales with the number of oscillators N as λ crit ∝ log(N). This raises the question if such a scaling law can also be observed in our models: observing the same scaling laws in both the Kuramoto model and the micromagnetic solution would strengthen the suggestion of the Kuramoto model as a valid description of arrays of STO.
We first consider our Kuramoto model, which has a non-local interaction to mimic the dipolar interaction in arrays of STO. Due to the increased complexity compared to the nearest neighbor model studied by Lee et al. 48 , an analytical derivation of the scaling behavior with system size is to our knowledge still an open question. To investigate the scaling behavior we thus resort to a numerical solution. We performed simulations with the number of oscillators ranging from N = 9 to N = 2500, gradually increasing the interaction strength between each simulation until the system reaches a synchronized state at a critical coupling strength, λ crit . 100 simulations were performed for each system size, with different initial oscillator phases and eigenfrequencies. In Fig. 3a we show a plot of λ crit vs. number of oscillators, N. The results indicate that the critical coupling strength scales as λ crit ∝ log(N), same as the nearest neighbor Kuramoto model investigated by Lee et al. 48 . The main difference compared with the nearest neighbor model is that we include interactions within a coupling radius R, as indicated by the blue circle in Fig. 1c. The imposed cutoff radius has a physical justification when considering a realistic system, which would inevitably include thermal noise. As the STO we consider are weakly coupled, and the dipolar interaction decay with distance, there will be a limiting spacing where the thermal noise level is comparable to the coupling strength. In our model we observe that the results do not depend qualitatively on the range of the cutoff radius, and that the large scale behavior is dominated by the diffusive coupling. This suggests that insight from studying the analytically tractable nearest neighbor model of Lee et al. 48 might provide valuable insight into the behavior of STO arrays.
The coupling in the Kuramoto model is defined simply as an interaction strength given by λ in Eq. (4). In the micromagnetic model the coupling comes from dipolar interactions, determined by the magnetic material properties and the STO spacing. In order to investigate the scaling behavior in the micromagnetic solution, one thus needs to relate the critical interaction strength to a critical STO spacing. Following the aforementioned macrodipole approximation that the effective dipolar interaction decay as 1/d 3

Scaling of output power in arrays of STO.
To investigate the implications of the scaling with system size, we consider a model calculation of the output power as we increase the number of STO. For applications of STO as e.g. nanoscale microwave generators, the power output of a single STO is not competitive. Decisive improvement is expected from the synchronization and phase locking of several STO, as this would result in a quadratic increase of the output power, P ∝ N 2 for N synchronized oscillators.
The output of a single STO can be described by its amplitude and phase, θ a e j i j . In our model we have assumed a constant amplitude for all STO and the total amplitude A for an array of STO is given by: , from the definition of the order parameter ρ in Eq. (6). A quadratic scaling in the power output, P ∝ N 2 , implies a perfectly synchronized and phase coherent state, given by ρ = 1. However, as the coupling strength to obtain a synchronized state scales with the number of STO, this will affect the power output when scaling up to large arrays.
As an example, we consider a system of STO composed of 200 nm diameter spin valve nanopillars with 15 nm thick Permalloy as the ferromagnetic layer. The average interaction energy can be extracted from micromagnetic simulations, and for an edge-edge spacing of 150 nm the interaction strength is found to be λ≈  7-8 MHz 31 . To increase the output power, we are now interested in scaling up to a large number of STO. We keep the same STO spacing when scaling up, e.g. keeping λ fixed. From Fig. 3a, we see that for an interaction strength λ  , there is a corresponding number of STO, ∼ N , where λ λ >  crit . This can be illustrated by calculating the total power output as the number of STO is increased, as shown in Fig. 3c. For a small number of STO, we see that the power output follows close to the ideal N 2 scaling. However, as the number of STO is increased, the scaling with system size becomes increasingly important. For a certain number of STO, indicated by ∼ N in Fig. 3a, the interaction is no longer strong enough to obtain a synchronized state. This is also illustrated in the inset of Fig. 3c, where we plot the order parameter ρ for arrays of 3 × 3 and 13 × 13 STO respectively, showing the transition from a synchronized state to chaotic behavior as the number of STO is increased above ∼ N . This illustrates the importance of our findings that synchronization in such 2d arrays is purely a finite size effect. The interaction strength is limited by the material properties and STO spacing, and using realistic values of the coupling strength we start to see significant deviations from the ideal P ∝ N 2 scaling for array sizes larger than 10 × 10 STO (see Fig. 3c). This means that in a physical realizable system, the scaling with system size imposes an upper limit for the maximum number of STO that can be synchronized.
Topological defects. That the synchronized and phase coherent state is purely a finite size effect, is similar to that of the classical 2d XY model of magnetism. The Kuramoto model is indeed similar to the 2d XY model 49 , where the direction of spin in the XY model corresponds to the oscillator phase in the Kuramoto model. In the 2d XY model a long range ordered phase is absent due to the presence of spin wave fluctuations and topological defects. The lack of long range order is a specific case of the Mermin-Wagner theorem in spin systems 50 , stating that continuous symmetries cannot be spontaneously broken in systems with sufficiently short range interactions in dimensions d ≤ 2. The fluctuations preventing long range order in the 2d XY model diverge logarithmically with system size 49 , in agreement with the logarithmic scaling observed in our system of STO.
Similar topological defects at the boundaries between locally synchronized clusters have previously been observed in the nearest neighbor Kuramoto model 48 as well as in other two-dimensional oscillator network models 17,[51][52][53] . In oscillator networks this is associated with the appearance of topological defects in the oscillator phase field, θ i . In the continuum limit this is expressed as: where dl is an integration path enclosing the defect, and n is the topological charge. Such topological features are observed also in our Kuramoto model for arrays of STO. The presence of vortices in the phase field is more pronounced as the system size increases. As an example we here consider an array of 50 × 50 oscillators. The disorder in the system is kept constant (given by the difference in the nominal frequencies of the oscillators) and the interaction strength λ is varied, acting as the inverse temperature: as coupling increases, the system becomes more ordered. Starting from a disordered initial state and varying the interaction strength between each simulation, we observe 4 different regimes: For weak interaction strengths we observe the formation of locally synchronized clusters, where cluster sizes increase with interaction strength. Apart from the localized clusters there is no long range order in the system (indicated as regime 1 in Fig. 4). Increasing λ above a certain threshold, we enter regime 2. Here we observe the formation of vortices in the phase field, and as an example we show in Fig. 4c a state with 4 vortices. The topological charge is conserved in the system, and two vortices of charge ± 1 respectively is present.
In both regime 1 and 2 long range order in the system is absent and ρ ≈ 0, as indicated in Fig. 4a. (that ρ > 0 here is a result of fluctuations due to the finite array size). Increasing λ further we enter regime 3, where the transition from regime 2 → 3 is governed by vortex annihilation processes (see 'Supplementary information'). Here there are no topological defects in the phase field, and the lack of global phase coherence is due to spin waves in the phase field where the oscillator phases change smoothly across the array (regime 3 in Fig. 4b). For sufficiently strong interactions we enter regime 4. Increasing the interaction strength is analogous to increasing the exchange coupling in a Ferromagnetic system, resulting in a more ordered state. The result here is a gradual suppression of the spin waves in the phase field observed in regime 3 as the interaction strength is increased. Regime 4 is thus characterized as the globally phase coherent state where all oscillators are synchronized and phase coherent (regime 4 in Fig. 4b).
The growth of the order parameter ρ with increasing coupling strength λ in Fig. 4a resembles that of a phase transition. Previous work have shown that the synchronization transition in the globally coupled Kuramoto model can be described as a phase transition, where the nature of the transition can be of first or second-order depending on the frequency distribution and coupling topology 8,54 . The Kuramoto model with finite range couplings is less studied, as these systems are difficult to analyze and solve analytically. A study of the locally coupled Kuramoto model on a d-dimensional lattice have shown that the synchronization transition depends strongly on the lattice dimensionality, and indicates d = 4 as the lower critical dimension for phase synchronization 55 . This is in agreement with the observed scaling with system size in our model, which indicates that the synchronization transition is purely a finite size effect.
In order to investigate the synchronization transition in our model further, we calculate the spatial correlation function for the array of oscillators. The correlations decay with distance, and asymptotically the correlation function is given by: 〈 θ(r) ⋅ θ(R)〉 ∝ e −|r−R|/ξ /|r − R| η . This describes the correlation between oscillators at positions r and R respectively, and the correlation length ξ is obtained by averaging over all positions r and R in the array (an example is shown in 'Supplementary information' Fig. S2). From the decay of the correlation function, we Scientific RepoRts | 6:32528 | DOI: 10.1038/srep32528 obtain the correlation length ξ as a function of the interaction strength λ. Conventional phase transitions are accompanied by a diverging correlation length close to the transition. Here, we do not observe a diverging ξ going from the disordered to the phase coherent state (1 → 4 in Fig. 4), and the correlation length remains finite. As inset in Fig. 4a we show a log-log plot of the correlation length normalized to the system size, ξ/L, (N = L × L).
The results indicate a power law relating the correlation length and the interaction strength as ξ ∝ λ ν , where the exponent for this case was found to be ν = 2.1 ± 0.1. This means that the correlation length simply scales with the coupling strength, and the transition between regimes 1-4 in Fig. 4 correspond to structures of ever increasing length scales. The transition to the phase coherent state (ρ → 1) occurs when the correlation length approaches the system size L, underlining the finite size effects on the synchronization transition and that the system is not undergoing a conventional phase transition. Further investigations of finite size effects, the lack of long ranger order in the Kuramoto model and the connection to the 2d XY model will be the subject of future work.

Discussion
To summarize, we have shown that the Kuramoto model provides a good description of arrays of STO. It provides a simple theoretical model to study large populations of coupled STO, which were previously unaccessible due to the long computation time for a full micromagnetic solution. By investigating the collective dynamics in large arrays of STO, we observed a scaling with system size indicating that the synchronization in arrays of dipolar-coupled STO is purely a finite size effect. The critical coupling strength to obtain a globally synchronized state scales with the number of oscillators, as λ crit ∝ log(N), preventing global synchronization for large system sizes. As a consequence of the scaling with system size, we showed that for realistic values of the dipolar coupling strength between STO this imposes an upper limit for the maximum number of oscillators that can be synchronized. Further, we showed that the lack of long range order and scaling with system size is associated with the emergence of topological defects and the formation of patterns in the phase field, similar to that of the 2d XY-model of magnetism.
In the present study we considered dipolar-coupled STO, where the short time delay in the coupling between neighboring oscillators compared to the oscillator frequency means that phase delay in the couplings can be neglected. However, for other coupling mechanisms, time delay can become significant. Interaction mediated by spin waves provide a different mechanism to obtain synchronization of STO 32 , where the finite propagation speed of the spin waves results in a phase offset in the couplings. Another recent proposal includes the use of non-local electrical couplings, where the coupling phase can be externally tuned through an electrical delay line 56 .
From a dynamical systems point of view, the study of time delay induced modifications to the couplings is of fundamental interest, as well as of practical relevance for modeling of physical, biological and chemical systems. In such systems, time delay is associated with finite propagation velocity of the couplings via e.g latency times of neuronal excitations, reaction times in chemical systems etc. (see e.g. ref. 9 and references therein). The possibility of designing STO arrays with a defined phase offset in the couplings suggests a real world analog to the more general Sakaguchi-Kuramoto model on a 2d lattice, which allows for a phase lag in the couplings 57 .
Our study suggests, on the one hand, that the use of models from non-linear dynamics can be useful for describing synchronization of magnetic oscillators and, on the other hand, arrays of STO as a a physical realizable model system for the Kuramoto model on a 2d lattice.