Chimeras and complex cluster states in arrays of spin-torque oscillators

We consider synchronization properties of arrays of spin-torque nano-oscillators coupled via an RC load. We show that while the fully synchronized state of identical oscillators may be locally stable in some parameter range, this synchrony is not globally attracting. Instead, regimes of different levels of compositional complexity are observed. These include chimera states (a part of the array forms a cluster while other units are desynchronized), clustered chimeras (several clusters plus desynchronized oscillators), cluster state (all oscillators form several clusters), and partial synchronization (no clusters but a nonvanishing mean field). Dynamically, these states are also complex, demonstrating irregular and close to quasiperiodic modulation. Remarkably, when heterogeneity of spin-torque oscillators is taken into account, dynamical complexity even increases: close to the onset of a macroscopic mean field, the dynamics of this field is rather irregular.

prototype model for such a coupling has been suggested, where N STOs are connected in series and are subject to a common dc current, with a parallel resistive load. The coupling is due to the giant magnetic resistance (GMR) effect, as the resistance of an STO depends on the orientation of its magnetization, so that the redistribution of the ac current between the STO array and the load depends on the the average (over the array) value of this resistance. This setup thus corresponds to the general scheme of mean-field coupling, as discussed above. However, application of the standard Kuramoto approach here 43,47 is rather questionable, because the STOs are highly nonlinear (see refs 44-46, 48-50) and are especially sensitive close to the homoclinic gluing bifurcation 45,51,52 .
In this paper we study an array of electrically coupled STOs with an RC-load. This setup is close to that considered in refs 19, 20. We will show that the STOs demonstrates a plethora of nontrivial dynamical regimes, including partial synchronization, pure and clustered "chimera" states, and clustering. Furthermore, the mean field appears to be highly irregular in a large range of parameters, making the transition to synchrony in this ensemble quite different from the usual transitions observed, e. g., in the Kuramoto model.
Before proceeding to the results, we make some remarks about the terminology used in this paper and its relation to the terms used in the literature. Most controversial is the definition of chimera states. Commonly, this term refers to situations where elements or their groups differ with respect to a certain characteristics or property. In a traditional approach, this characteristics is the individual frequency of oscillations, averaged over time: if frequencies of oscillators in a symmetric setup differ, one calls this state a chimera (see, e.g. ref. 53). Another widespread approach is based on spatial organization: one defines chimera as a profile with a coexistence of smooth and discontinuous in space patches (see, e.g. refs 54, 55). Each approach, along with merits, has its obvious limitations: the former, being restricted to the cases of continuous-time oscillators with well-defined mean frequency, is applicable neither to the ensembles of coupled maps 54,55 , nor to the spin dynamics 56 ; the latter approach assumes spatial ordering and cannot be applied to interactions in globally coupled arrays that are insensitive to distances in physical space. Below we treat a situation where none of the above approaches is applicable: we discuss globally coupled units in the regime of strong coupling where existence of a well-defined mean frequency cannot be guaranteed. Instead, below we adopt a definition of "chimera" state, based on the clustering property (the respective characteristics is the instantaneous position of the unit in its state space). This approach is applicable to globally coupled populations both of discrete maps and of continuous in time systems. We call a state in a population of identical units a "chimera" if it consists of a macroscopic cluster of units, the states of which coincide, and of a cloud of units, the states of which are different. Notice, that we put in this paper the term "chimera" into quotation marks, to emphasize the usage of a particular definition which may differ from other definitions existing in the literature.

Results
Basic equations. In this paper we consider an array of spin-torque oscillators, connected to the external driving current via an RC load. The external, generally time-dependent current I(t), is therefore divided between the current through the STO-array J and the current through the RC-load (capacitance C, resistance r) Here V = JR is the voltage on the array, depending on the (generally time-dependent) resistance R of the STO elements. This relation gives the first differential equation of the system The remaining equations are those for N STO oscillators. Each STO is described by its free-layer magnetization M i . This vector has constant unit length, and its orientation varies according to the Landau-Lifshitz-Gilbert-Slonczewski equation where γ is the gyromagnetic ratio; α is the Gilbert damping constant; β contains material parameters; J is the current through the STO defined above; the effective magnetic field H eff contains an external magnetic field, an easy-axis field, and an easy-plane anisotropy field; M 0 is magnetization of the fixed layer. Following 44 we assume . The Landau-Lifshitz-Gilbert-Slonczewski equations can be rewritten in terms of the spherical coordinates (φ i , θ i ) determining orientation of M i , these equations are coupled via the current J = V/R (see Supplementary Materials for details): Finally, the interaction between the dynamics of the STOs and the load voltage V is due to the dependence of the resistance R on the magnetization vectors M i , according to the mechanism described in ref. 19: sin cos 1 sin cos (5) The main parameters that determine the coupling are the level of magnetoresistance variations ε, and the dimensionless time constant of the load τ. The main parameter determining the dynamics of the STOs is the external current I. Below we will study the dynamics of the array in the cases of identical and non-identical STOs. The mean magnetization X will be used for the visualization of the mean field behavior.

Identical oscillators.
For identical spin-torque oscillators, different cluster regimes, where the states (i.e. variables (θ I , φ)) of different groups of units coincide identically, are possible. The simplest of these regimes is that of full synchrony (one-cluster state), where The full system of 2N + 1 equations then reduces to a three-dimensional one; for large enough values of the external current > = + + α β 2 dz the reduced system possesses a limit cycle solution, born in the supercritical Hopf bifurcation at I = I H . Mean field in this case coincides with the field of a unit. This regime has a chance to be observed if it is stable towards "evaporation" of units from the cluster. Quantitatively, this stability is measured by the "evaporation" or "split" Lyapunov exponent 57 , see section Methods for details. A negative evaporation Lyapunov exponent means that the cluster is stable, a positive one indicates instability. We show the stability diagram in Fig. 1. One can see the large range of the parameters of the external load τ and of the current I, where the full synchrony is stable. This result is confirmed by direct numerical calculations: when the initial states are prepared sufficiently close to each other, evolution of the array ends up in the fully synchronous periodic regime.
Direct calculations however show that the full synchrony is not globally stable. Starting the simulations of Eqs (5) from random initial conditions, we never observed the fully synchronized state after transients. Instead, different complex configurations occur. These configurations can be classified as follows (for a practical implementation of this classification see section Methods below): 1. Partially synchronous states. In this regime the states of all oscillators are different (i.e. not a single persistent cluster is observed), however the distribution is not uniform so that the mean field X(t) performs macroscopic oscillations. These states are of the same type, as described in refs 10-12, 58 for other models of globally coupled oscillators. Noteworthy, in the partially synchronous states described in these references, the mean frequencies of all oscillators are the same. However, the notion of partial synchronization can be also applied to systems without a well-defined frequency.  The dynamics here is high-dimensional, although the dimension is significantly reduced compared to the dimension 2N + 1 of the original system. 4. Clustered "chimera" states. Here several clusters are formed, but also a significant number of non-clustered units in a cloud is present. This regime can be considered as a combination of regimes 2 and 3.
We stress here that this classification is purely compositional one, it can be inferred from the instantaneous snapshot and does not rely on the particular dynamics of the STOs or of the mean fields; rather we discuss below how/whether dynamics is related to the compositional properties.
In Fig. 2 we show the statistics of these states in a population of N = 200 spin-torque oscillators (this number is not very high because we needed here a large set of trials; below, when illustrating different states, we use ensembles with sizes up to N = 1000). Here the parameter of the external load is fixed at τ = 6. The computational procedure was as follows. For every checked value of I we localized the limit cycle corresponding to the fully synchronous state. Given a point θ φ v ( , , ) on this limit cycle, we randomly chose ca. 1800 initial conditions v v v 0 1, and determined the proportions of those conditions that, after a transient time of 2000 characteristic periods, ended up in one of the four states described above.
One can see in Fig. 2 that there are regions where the partial synchronous states dominate, and a large region   .
. I 0 0063 0 0115 where partial synchrony is never observed, but instead clustered and "chimera" states are dominant. Remarkably, in this region the fully synchronized state is stable, so that possibly this stability is inherited by sufficiently big clusters.
Above we characterized compositional properties of the complex states, by classifying them according to cluster structure; next we describe dynamical properties of the evolution of the mean field. For the most of the values of the current I presented in Fig. 2, nontrivial aperiodic dynamics of the mean field X(t) has been observed. We show three typical examples in Fig. 3. To make the time dynamics clear, we present it on two panels, at a long and at a relatively short time scales. Characterization of irregularity of the mean field in a system of many interacting units is not an easy task. Indeed, here the usual methods suitable for low-dimensional systems, such as calculation of Lyapunov exponents and of dimension are hardly applicable. Therefore we characterize regularity of the mean field by calculating the autocorrelation function (see section Methods for details). For small values of the current I, a very slow irregular modulation of the oscillations is observed (Fig. 3, panels (a,b)). The autocorrelation function (panel (c)) slowly decays in this case. For larger values of I the modulation is quite ordered, close to a quasiperiodic regime with two frequencies (Fig. 3, panels (d,e)). Here the autocorrelation function (panel (f)) returns nearly to 1, as one expects for quasiperiodic processes. For even larger values of the current, the dynamics is strongly irregular (Fig. 3, panels (g,h)). Here the autocorrelation function (panel (i)) rapidly decays nearly to zero.
A natural question is whether the dynamical properties of the array of STOs are related to the compositional ones. We have found that this relation is very weak . For example, regimes (a,b) and (e,f) in Fig. 3, while being both compositionally partial synchronous states, demonstrate in one case quite regular dynamics and in another case rather irregular one. The case, where for the same value of the parameters different structures are possible (like in panels (c,d) in Fig. 3) requires special consideration. Here, for I = 0.008, according to diagram Fig. 2, clusters, clustered "chimera", and "chimera" states can occur. In Fig. 4 we compare the dynamical regimes for these three compositional states. Together with the time course of the mean field, we present in the right panels the Poincaré maps. Recall that a quasiperiodic attractor with two incommensurate frequencies would be seen as a smooth closed curve on the map. One can see that in all the cases the dynamics is close to a quasiperiodicity, with deviations much more noisy in "chimera" (a,b) and clustered "chimera" (c,d) regimes compared to the clustered state (e,f). This has quite natural explanation: clustered state has a rather small dimension compared to situations where there is a large cloud of unclustered oscillators, thus the effective noise is smaller.
Nonidentical oscillators. Above we focused on the properties of an array of identical STOs. Here we study a more realistic situation, where parameters of the STOs are different. Following ref. 19, we take parameter H k in the effective magnetic field → H eff to be uniformly dispersed in a range H k0 − ΔH k < H k < H k0 + ΔH k . According to general theory of the synchronization transition, e.g. from the exact solution of the Kuramoto model 1, 2 , one expects, that for a fixed coupling strength and strong enough diversity, no synchronization is observed, and this state turns into a synchronous regime with typically periodic behavior of the mean field, through a bifurcation at a certain critical value of diversity. The bifurcation is perfect in the thermodynamic limit (infinite number of oscillators) and is spoiled by finite-size fluctuations for finite ensembles.
Our numerical simulations only partially support this scenario, see Fig. 5. For large diversity in the array, indeed the oscillators do not synchronize and the mean field vanishes (see panels e,f). However, we do not observe, for a fixed level of coupling strength, a transition to a regular, periodic mean field like in the Kuramoto model: the mean field appears to be highly irregular close to the transition (panels c,d) and even for a small diversity (panels a,b). This holds for the case where identical oscillators demonstrate erratic behavior (case I = 0.012, panels b,d,f), as well as for the case where for identical oscillators a regime close to quasiperiodicity is observed (case I = 0.008, panels a,c,e). This effect is opposite to "taming spatiotemporal chaos with disorder" reported in ref. 62. Here disorder enhances collective chaos, as is especially pronounced by comparing Fig. 3(c) with Fig. 5(a,c). A more detailed exploration of the range of parameters is presented in the Supplementary Material, Figs 3 and 4. There we show the dynamics of the mean field for different levels of disorder and for different coupling strengths. One can see there that only for weak disorder and weak coupling above the synchronization threshold, the mean field is regular; otherwise irregular variations are observed. We have to stress, however, that finite-size effects may hide or destroy regularity of oscillations, therefore more detailed studies of very large arrays are needed to clarify this observation.

Discussion
In this paper we studied complex collective regimes in an array of spin-torque oscillators coupled via a common RC load. Our main focus was on the compositional properties of the ensemble, and on the dynamical properties of the mean field. In almost all aspects, these properties are very much different from those in the standard models of globally coupled periodic oscillators, like the Kuramoto model. We did not use for classification of the states the properties of the mean frequencies, because the STOs in some regimes demonstrate mixed mode oscillations (see Fig. 2 of Supplementary Material), for which the unique frequency cannot be defined straightforwardly.
Compositional properties can be most clearly described for the identical oscillators. Here, although stability theory predicts a possibility of a fully synchronized state, practically only complex regimes between full synchrony and full asynchrony are encountered. In a large range of parameters we observed a "chimera" state, where  part of the population builds a cluster and the rest remains disperse. While this regime has been recently reported for several models [13][14][15] , a more complex clustered "chimera" state, with several clusters and dispersed oscillators, appears to be a novel one. Together with these two types of "chimera", a clustered state where all oscillators belong to several clusters, is also possible. Alternatively to appearance of clusters, the array can demonstrate a partially synchronous state where all oscillators are dispersed but remain correlated forming a relatively strong mean field. We reiterate that our definition of "chimera" is not based on the comparison of mean frequencies of oscillations, as has been suggested in some recent studies. Indeed, the frequency is a convenient observable for phase oscillators (or, e.g., for those obeying differential equations on a cylinder, where rotations around the cylinder can be straightforwardly counted). In our case, every individual oscillator lives on a sphere, and it cannot be guaranteed that trajectories do form in this projection a "band attractor" allowing to count oscillations and determine the mean frequency. Proper definition of a mean frequency and its analysis in different regimes may be a subject of future work. Dynamical properties are related to the regularity of the mean field, which is well defined for nonidentical oscillators as well. One expects that for strong diversity of the oscillator's properties, the correlation between them will be minimal and they should sum up incoherently into a constant mean field. We confirm this property for the spin-torque oscillators as well. Close to the desynchronization threshold, usually the appearing mean field is dynamically simple (periodic in time). For example, in the case where for identical oscillators one observes complex behavior related to a heteroclinic cycle 63 , by following the transition from large diversity to the small one, in ref. 64 first a transition to small periodic oscillations was detected; these oscillations undergo secondary bifurcations to more complex states only when diversity becomes sufficiently small. For the spin-torque oscillators we observe a different scenario: here close to the transition to partial coherence where the mean field is small, we observe rather irregular dynamics. For smaller disorder, the dynamics becomes more regular, close to a quasiperiodic one, for clustered states, and remains rather irregular for the partial synchronous state. This interesting observation deserves more detailed study, where especially the role of the finite size effects should be clarified.
Finally, we discuss how general are the results described above. In this paper we focused on a particular setup, where the common load for the array of STOs is of RC-type. In refs 20, 51 an array with an RCL-load has been i i at the moments where X reaches local maxima. We stress here that the states are classified ("chimera", clustered "chimera", clustered) not according to the behavior of the mean field, but through the analysis of the instantaneous states of the units in the population, as defined above and described in section Methods.
considered. Preliminary calculations for this case show, that all four types of the composition are observed in this setup as well; however the frequency of their occurrence and the dynamical properties are different.

Methods
Integration of equations. Simulation of the system (5) of ODEs governing the array of spin-torque oscillators is performed by virtue of the 4th-order Runge-Kutta method. In the case of identical oscillators, there is a numerical trap due to a finite precision of computer representation of real numbers: if the states of two units are closer than the accuracy of representation of numbers in double precision (16 decimal digits), these states "merge" and remain henceforth indistinguishable. This pitfall may lead to appearance of clusters which otherwise (e.g. with quadruple precision calculations) would not exist. To avoid this effect, we add a very small diversity (of the order 10 −12 ) to parameters of the oscillators, thereby preventing exact coincidence of their dynamics.

Evaporation Lyapunov exponent.
To determine the stability of the cluster solution (6), one imposes a lin- to the permutation symmetry such a perturbation can be imposed to any pair of oscillators). This perturbation does not influence the mean field X and thus also the global variable v, it just corresponds to a small deviation of two oscillators from the cluster. The equations for δθ and δφ thus read Integration of these linear equations together with the full system (5)  which determines whether the cluster is transversely unstable (λ ev > 0) or stable (λ ev < 0). Distinguishing different compositional states. We used the following algorithm to distinguish different compositional states. In a population of STOs, we first identified clusters. Two units with indices k, m are assigned to the same cluster if their states are sufficiently close: d(k, m) < d cl , where d(k, m) = |x k − x m | + |y k − y m |, d cl is a small number (we used in calculations d cl = 10 −10 ), and the local observables x k , y k are defined in the same way as the mean fields x k = sinθ k cosφ k , y k = cosθ k sinφ k . If no clusters were detected in the population, (but the macroscopic mean field was present, that has always been the case) the state was classified as a partially synchronous one. The "chimera" state was the state with just one cluster which however did not include all units. We defined the clustered "chimera" state as one with more than one cluster, and a sufficient (larger that 10% of all units) number of units not belonging to the clusters. Finally, a state with many clusters and possibly a small amount (less than 10% of all units) of units outside the clusters was identified as a multiclustered state. This algorithm was used to produce Fig. 2.
Autocorrelation function and Poincaré map. The autocorrelation function of the mean field X(t) is . Practically, the autocorrelation function was calculated via time averaging; the length of a time series was 16000 characteristic periods. Another way to characterize the complexity of the dynamics is a Poincaré map. In our case we have chosen a two-dimensional map, by plotting all the points at which the mean field X(t) reaches a maximum, on the plane (X, Y) where the mean field Y is defined as For the periodic motion attractor of this map reduces to a point, while for a quasiperiodic motion with two incommensurate frequencies it yields a closed curve. Poincaré maps presented in Fig. 4 are close to this structure. More complex states (quasiperiodic ones with many incommensurate frequencies and chaotic ones) result in a cloud of points lying on a fractal set (for a low-dimensional chaos) or without any structure.