Glassy states and super-relaxation in populations of coupled phase oscillators

Large networks of coupled oscillators appear in many branches of science, so that the kinds of phenomena they exhibit are not only of intrinsic interest but also of very wide importance. In 1975, Kuramoto proposed an analytically tractable model to describe these systems, which has since been successfully applied in many contexts and remains a subject of intensive research. Some related problems, however, remain unclarified for decades, such as the existence and properties of the oscillator glass state. Here we present a detailed analysis of a very general form of the Kuramoto model. In particular, we find the conditions when it can exhibit glassy behaviour, which represents a kind of synchronous disorder in the present case. Furthermore, we discover a new and intriguing phenomenon that we refer to as super-relaxation where the oscillators feel no interaction at all while relaxing to incoherence. Our findings offer the possibility of creating glassy states and observing super-relaxation in real systems.

T he Kuramoto model (KM) 1 was introduced and developed to provide an analytically tractable description of the populations of coupled phase oscillators that so often appear in real life. It has been applied successfully in many fields 2,3 , for example, to describe the collective behaviour of lasers 4,5 , neurons in the brain 6,7 , Josephson junction arrays 8 and even humans 9 .
The widespread applications of the KM to real problems has ensured that the basic model, together with its variants and modifications [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] , has been studied very extensively over the past few decades. The practical usefulness of these studies stems from widespread abundance of user-defined systems that can be described by the KM, for example, the above-mentioned laser arrays 4,5 and Josephson junction circuits 8 , where one can adjust the coupling between the oscillators quite freely. Therefore, if some new behaviour has been uncovered theoretically, it can immediately be implemented and observed in practice. Furthermore, for those systems where the exact configuration is not known precisely and cannot be changed, such as interacting neurons in the brain 6,7 , it is obviously useful to know what kinds of dynamics the different forms of KM can demonstrate, so that the observed behaviour can be matched with an appropriate existing model and thus explained.
The KM has been found to exhibit a remarkable diversity of interesting phenomena and states and, by now, most of them have been thoroughly investigated. There is, however, one important exception-the so-called 'oscillator glass' state 25 . Thus, many systems, including networks of interacting spins 26,27 , dipoles 28 , and electrons 29 , have been found to display behaviour reminiscent of a glass structure 30 , and it has been suggested 25 that populations of coupled oscillators can also demonstrate glassy states of some kind. Studies of such states are expected to be very fruitful, in the same way as were the studies of spin glasses, which led to many new techniques and applications in other fields (biology, computer science, economics and so on 27 ). However, the provenance and properties of the oscillator glasses are still subject to debate 25,[31][32][33][34][35][36] .
Here we report the analysis of a very general, yet analytically tractable, form of KM. We derive a set of equations describing its steady-state behaviour and show that, for a particular class of distributions, the consideration may be reduced to a much simpler model. Thus, in the limit t-N, the system can demonstrate the same macroscopic behaviour for different coupling configurations, enabling one to select and study the simplest one. This is, in itself, an interesting theoretical result. It also implies that there may be cases where it is, in principle, impossible to infer the underlying coupling structure from the observed data.
We also discuss when and how the full-time evolution of the system parameters can be obtained and find some interesting features related to this question. Among the phenomena that the model can exhibit, we find states with a glassy structure that can be studied analytically within the framework presented, and we discuss these states and their properties in detail. Finally, we describe a completely new phenomenon which we refer to as super-relaxation where, under certain conditions, the coupling between the oscillators effectively disappears during their relaxation to incoherence.

Results
Outline of the results. To make what follows clearer and easier to understand, we start by summarizing the main results of the work. First, we derive equations (18)(19)(20)(21)(22) describing the stationary states (SSs) of the system (1). Then we reduce the macroscopic steady-state behaviour of the system (1) (five distributed parameters) to a much simpler model (29) (two distributed parameters) for a large class of distributions (equation (23)). We discover new states in the coupled oscillator populations, which are in some sense similar to physical glasses, and we establish the conditions for their appearance (equations (36) and (37)). And finally, we discover the emergence of a super-relaxation phenomenon, where the oscillators evolve interaction-free, irrespectively of the couplings between them, if the conditions (41), (42), (43) are fulfilled. After reviewing the background and terminology, each of these results is obtained and considered in detail below.
The model. We consider the KM in the form where N is the number of oscillators, y i (t) and o i are the ith oscillator's phase and natural frequency, respectively, and k i q j (b i þ g j ) represent the coupling strengths (phase lags) between the ith and jth oscillators; all parameters The case considered earlier 10 corresponds to (1) with The KM has not previously been treated in such a general form (1), though it actually includes as special cases the many KM modifications and extensions studied earlier. Thus, not much is known about the possible behaviour of the system except in those particular cases. To study it, we first generalize the recently presented framework 10 to encompass (1), and then we proceed to a consideration of the phenomena that it can exhibit.
Main equations. The oscillators' collective behaviour in (1) can be described by two complex parameters where Z is the mean field whose amplitude R quantifies the extent of the agreement between the oscillators' phases y i , while Y represents the weighted mean field, with amplitude W reflecting the agreement between y i þ g i ( þ p for q i o0), weighted by |q i |. With the use of the definition (2), the model (1) can be rewritten as As can be seen, the dynamics of the system is governed mathematically by the weighted mean field Y, which determines the effective interaction between the oscillators. The ordinary mean field Z, on the other hand, represents a more physical macroscopic variable, being the one that is most often observed in practice (for example, it quantifies the overall output current of a Josephson junction array 8 ). Note that, unlike the mean fields' phases F and C individually, their difference F À C is invariant under the phase shifts y i -y i À j(t), and thus represents another meaningful parameter. In the continuum limit N-N, the system (1) is treated using the probability density function f(y,C,t), which reflects the probability that the oscillator has parameters C and phase y at time t. It can be further factorized as f ðy; C; tÞ rðy; t j CÞGðCÞ ð 4Þ where r(y,t|C) is the conditional probability density function (CPDF), reflecting the probability that at time t the oscillator has a phase y, given its parameters C. By definition, the CPDF should satisfy R p À p rðy; t j CÞdy ¼ 1, which leads to the continuity equation @ t rðy; t j CÞ þ @ y f½o À kW sinðy À b À FÞrðy; t j CÞg ¼ 0; ð5Þ where we have used expression (3) for _ y. The CPDF can usually [37][38][39][40] be represented using the ansatz introduced by Ott and Antonsen (OA-ansatz) 41 rðy; t j CÞ ¼ 1 2p 1 þ 2Re aðC; tÞe À iy 1 À aðC; tÞe À iy ! ; j aðC; tÞ j 1 ) Z p À p e iy rðy; t j CÞdy ¼ aðC; tÞ; ð6Þ where we use a(C,t) ¼ a*(C,t) instead of the original OAvariable 41 a(C,t), and the * denotes complex conjugate.
Finally, by substituting equation (6) into equations (5) and (2), one obtains the full system of equations describing the dynamics of the system (1): where here and below, unless otherwise specified, the integration over dCdodkdqdbdg is taken over the whole domain (( À N, N) for o, k, q and ( À p, p) for b, g).
Parameter redundancy. It should be clarified that parametrizing the KM (1) with both b i and g j is, in a strict sense, mathematically redundant, since one of these phase shifts can be removed by a change of variables. Thus, in terms ofỹ i ¼ y i À b i , the system equation (1) becomes with the new distribution of parameters e Gð e CÞ being given as e Gðo; k; q; e gÞ ¼ R Gðo; k; q; b; gÞdðe g À b À gÞdbdg. Taking into account the relationship between the variables in equations (10) and (1), it is clear that the conditional distribution of the new phases e rðỹ; t j o; k; q; e gÞ is related to the original one as rðy; t j o; k; q; b; gÞ ¼ e rðy À b; t j o; k; q; b þ gÞ. Substituting this into equation (2), it can be shown that the weighted mean fields in terms of y i andỹ i are equal, e YðtÞ ¼ YðtÞ, but that the ordinary mean fields can, in general, change in a non-trivial fashion, e ZðtÞ 6 ¼ ZðtÞ.
Nevertheless, both b i and g j can sometimes be practically relevant, that is, have physical meanings. For example, such a coupling structure can be designed artificially in the case of userdefined systems (for example, laser arrays 4,5 ), but we expect it to appear in the observation-only systems as well (for example, interacting neurons 6,7 ). Hence, in some cases, the full model (1) will actually provide a more straightforward and meaningful description of the system, while equation (10) will represent only a mathematical formulation. Thus, for example, Z might be related to the real physical quantity, such as the total power output of the laser array 4,5 , while e Z will then be a purely mathematical characteristic. To preserve generality, therefore, and to avoid complicating the discussion by introducing different variable transformations, we study the model (1) including both b and g. Furthermore, as will be seen below, the distributed b also has significant consequences in terms of the non-equilibrium dynamics, which cannot be straightforwardly obtained from the transformed system equation (10).
It should also be noted, that the model (1) is invariant under (k i ,b i )-( À k i ,b i þ p) or (q i ,g i )-( À q i ,g i þ p) or rescalings (k i ,q i )-(k i /r,rq i ); to preserve generality, we retain this ambiguity, as the normalization can be fixed at any time, and there is no universal choice (see Methods).
Terminology and notation. We adopt terminology similar to that introduced earlier 10,42 . The KM form (1) is invariant under y 0 ¼ y À Ot, which just changes the natural frequencies to . Thus, one can consider the KM in different frames rotating at frequency O with respect to some reference frame. For the latter, we select the frame with zero mean frequency /oS R oG(C)dC ¼ 0 and call it the natural frame; the distribution G(C) is defined in this frame. We define a SS as a state with time-independent CPDF q t r(y,t|C) ¼ 0, which also implies q t Z ¼ q t Y ¼ 0. The state can be stationary only in a particular rotating frame, so it is characterized by its frame frequency O and mean fields Z,Y. We refer to SSs with O ¼ 0 as natural states, while SSs with Oa0 are travelling wave states. For later convenience we define where 0rrr1; the distribution P(f,r) and the way to simulate it are discussed in Methods. Note, that P(f,0) ¼ 1/2p and As will be seen below, I and J represent an effective complex distribution of o,k to be used for the determination of Z and Y, respectively.
Validity of the OA-ansatz. In what follows, we will make frequent use of the OA-equations (7)(8)(9). For the case k i ¼ const, and for a very large class of frequency distributions, their validity has been proven 39,40 in the asymptotic limit t-N, corresponding to system's steady-state behaviour; though not justified rigorously, this result seems to hold for distributed k, q, b, g as well (based on numerical simulations). As the next step, Pikovsky and Rosenblum 37,38 derived more general equations and showed that the full dynamics of the model-for all t and almost any parameter distribution-obeys equations (7)-(9) only if the initial phase configuration also belongs to the OA-manifold rðy; 0 j CÞ ¼ Pðy À f 0 ðCÞ; r 0 ðCÞÞ ð13Þ with any f 0 (C) and r 0 (C). Although not given explicitly by Pikovsky and Rosenblum 37,38 , equation (13) can be deduced from their equation (3) in ref. 38, which generates (13) in the case of a uniform distribution of the constants of motion c k (for which the OA-description was proven to hold 37,38 ); see also the work of Marvel et al. 43 Note that the OA-ansatz often provides a good approximation to the system dynamics even when equation (13) is not satisfied (for example, for r(y, In addition, it might sometimes be useful to analytically continue a(C,t) into the complex plane over some of C, for example, to consider o to be complex. This is required to apply the conventional OA-reduction procedure 41 (see also refs 13-15) which, where possible, allows one to obtain a finite-dimensional system of equations for Y(t),Z(t). However, one should always have |a(C,t)|r1 as, otherwise, the OA-ansatz (6) becomes illdefined. Hence, a(C,t) can be considered only in the region of C for which j aðC; 0Þ j ¼ Z e iy rðy; 0 j CÞdy 1; ð14Þ Re½kðe ib e À ij À e À ib e ij Þ 0 at any W and jarg[a(C,t)] À F. The condition (14) establishes that |a|r1 is satisfied at t ¼ 0, while (15) then guarantees that it holds at all other times too. Some manifestations of the issues related to (14) and (15) (14) is always satisfied when there is no correlation of initial phases with the system parameters C: Unless otherwise specified, we will assume that the system starts from such a configuration.
Stationary states. Having reviewed the background and related issues, we now proceed to the derivations. We start by finding possible configurations into which the system (1) can settle as t-N. Although in specific cases it can converge to some inherently time-dependent solution, such as a standing wave 43 or an oscillating p-state 44 , we restrict the consideration to SSs only, treating them in frames where they are stationary. The effective parameter distribution therefore becomes G(C)-G O þ (C), with O denoting SS frame frequency. By definition and equation (6), the SSs satisfy q t a ¼ 0. Using this in equation (7) and taking account of the OA validity condition |a|r1, one finds that all possible SSs in their own rotating frames correspond to a s ðCÞ ¼ e iðF þ bÞ Â ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi The stationary distribution of the oscillators' phases is then recovered by substituting equation (17) into equation (6), which gives where H( À k) denotes a Heaviside function. As can be seen, the oscillators with |o|o|k|W (in a rotating frame) are frozen around the positions determined by b, o, k and the weighted mean field Y ¼ We iF , while the others are incoherent.
It should be noted that, in addition to a s (C) given by equation (17), there exists one more stationary solution of equation (7): the same as equation (17), but with a minus before the ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k 2 W 2 À o 2 p for |o|r|k|W; however, this solution is never realized in reality because it corresponds to the unstable position on the phase circle, as can be seen by recovering the respective CPDF (which will be the same as equation (18), but with y-y þ p for |o|r|k|W).
Self-consistency and stability conditions. While the solution (18) gives a qualitative picture of how the phases are distributed in the stationary regime, one cannot infer directly from r s (y|C) the stability of such a configuration or the macroscopic parameters W,O for which it can be realized. To find the latter, we first note that the equation (17) was derived from equation (7) only, without taking account of equation (8), which it should also satisfy. Therefore, substituting equation (17) into equation (8), one obtains the self-consistency conditions (SCCs) where J AE O are as defined in equation (12). Taking the real and imaginary parts of equation (19) yields two equations from which the SS parameters W,O can be determined. These can then be used to determine Z as well, for which, using equation (17) in equation (9), one gets The SS stability, on the other hand, is generally hard to study analytically. For the incoherent state (W ¼ 0), however, it can be analysed using the approach of ref. 45. Thus, performing a linear stability analysis of (5) above incoherence and using a selfconsistency argument (see Methods), one can show that incoherence changes stability when there exists a solution x to the equations Interestingly, from equations (21) and (12) it follows that in the case G(C) ¼ G(o, q, b, g)p(k) the stability of incoherence does not depend on a particular form of p(k), but only on the mean coupling strength /kS R kp(k)dk, as noticed previously for a simpler model 10,15 . Note also that, if phase shifts g or b are present then, with increasing coupling strength, the incoherence can not only lose, but also gain 46,47 stability at the transition points determined from equation (21).
To estimate (at least approximately) the stability of the SSs with W40, one can utilize the approach of ref. 10 and devise the empirical stability conditions (ESCs) which, for the considered model (1), take the form (see Methods) where e F e F½J AE O ; W (see equation (19)). Despite being empirical and thus approximate, the ESCs (22) work well in the majority of cases, though not in all, for example, they can fail in the presence of standing waves 10 ; their performance also becomes less good when arg[J(o,k)]a0. Nevertheless, ESCs seem to be exact if arg[J(o,k)] ¼ 0 and the distribution G(C) is unimodal over o.
Uncoupled distributions. As discussed above, for the system (1), the parameters of its possible SSs can be found from the SCCs (19) and (20), their stability can be deduced from equation (21) for the incoherent state, and (approximately) from equation (22) for other SSs, while the associated phase distributions are given by equation (18). The corresponding expressions, however, are generally quite complicated, but it turns out that they simplify greatly if we consider the distribution of q, b, g to be uncorrelated The SCCs (19) then simplify to are fully analogous to F R,O in ref. 10; the expression (20) for R and F À C reduces to the incoherence stability conditions (21) become System reduction. From equations (25)(26)(27)(28) it is evident that, for uncoupled distributions (23), all the macroscopic properties of the SSs are completely characterized by g(o, k), |J| and f J (see equation (24)), irrespectively of the particular form of h(q, b, g), while |I|, f I serve merely to specify Z (27). Instead of (1), therefore, one can consider the system with the same distribution of o i , k i defined by g(o, k). Then each SS of (29) corresponds to an SS of (1), and their parameters are related as where the subscript SK denotes the states of the model (29); the stability of the corresponding states will also be the same, as follows from equation (28) for incoherence, and from the ESCs (22) and numerical evidence for other SSs. Hence, for any distribution of q, b, g obeying (23), one can reduce the consideration of the steady-state behaviour of (1) to the much simpler Sakaguchi-Kuramoto model (29). This elegant and unexpected result is illustrated in Fig. 1. Interestingly, something similar was noted earlier 21 , but for a significantly less general model than (1). The reduction (30), however, relates only to the macroscopic properties of SSs (t-N), whereas the full evolutions Z(t),Y(t) and the microscopic properties of the resultant states cannot be obtained in this way. Note that equation (30) implies Wr|J|,Rr|I|, that is, the distribution of q, b, g imposes an upper bound on the SSs' mean field strengths, which they cannot exceed however strong the coupling is.
Definition of glassy states. Ordinary glass 30 is a state of matter which, in contrast to the liquid state, is solid; but, in contrast to the crystalline state, has a structure without any long-range translational order. States that are in some sense similar were reported in several other systems [26][27][28][29] , but whether or not there can exist glassy analogues in networks of coupled oscillatorsoscillator glasses-has remained unclear 25,31-36 .
However, a very basic question is how to define an oscillator glass, that is, what properties the system should possess in order to qualify for this title. Obviously, one should draw an analogy with known glassy states in order to illuminate this issue. Thus, the defining feature of the structural glasses 30 -frozen disordercan be specified as: (i) the particles, for example atoms, lack any long-range translational order, similarly to the liquid state; (ii) the particles are frozen, that is, do not move with respect to each other, similarly to the solid state. Later, states with frozen spins lacking long-range orientational order were discovered in spin systems, and by analogy were called spin glasses 26,27 . The canonical spin glasses, in addition to properties (i) and (ii), are characterized by: (iii) a significantly redundant ground state, so that the same macroscopic behaviour can be realized by many microscopic spin configurations, not related to each other by any simple symmetry transformation; (iv) frustration, that is, no spin configuration can simultaneously satisfy all energy bounds; (v) slow non-exponential relaxation of the order parameter (magnetization for spin systems), as well as other exotic dynamical features.
Relating the phases y i of the oscillators to the positions of the particles in space, or the spin orientations, the analogies of the above properties for populations of oscillators will be (all for t-N) Uniform distribution of phases at any time : rðy; tÞ R rðy; t j CÞGðCÞdC ¼ 1=2p; Existence of a ðfrozenÞ group of oscillators not moving with respect to each other : Q  (25)(26)(27)(28). The distribution for the original system was: gðo; kÞ $ dðk À KÞ½o 2 þ e À o 2 À 1 ; hðq; b; gÞ $ e À ðq À 1Þ 2 =0:02 ½Pðb À p=8; 0:8Þ þ Pðb þ p=2; 0:9ÞPðg À p=3; 0:7Þ, but the same picture appears for any h(q, b, g) with identical f J , |J|, |I| (24). The simulations were performed for 500 s, and the values presented are averages over the last 100 s. NATURE COMMUNICATIONS | DOI: 10.1038/ncomms5118 ARTICLE Nonuniqueness of lim t!1 ½y i ðtÞ À y j ðtÞ for the frozen oscillators; which might therefore end up in different places; depending on the initial conditions; Frustration; that is; ðarguablyÞ nonexistence of stable phase configurations with all _ y i ðtÞ ¼ 0; even if all o i ¼ 0; ð34Þ Slow nonexponential relaxation of RðtÞ: Generally, it is quite a subtle issue to decide what is a 'true' oscillator glass, that is, should it satisfy all the criteria (31)(32)(33)(34)(35), or only some of them? Thus, the state reported by Daido 25 has properties (31), (34) and (35), but not (32) and (33), since the oscillators there, although adjusting their mean phase velocities to some extent, perform a diffusive motion and thus, are never frozen (Q ¼ 0). This state, therefore, bears more analogy to spin glasses than to structural glasses.
Here, we utilize a different approach, considering state to be glassy if it satisfies the first two criteria (31), (32). To avoid possible terminological issues, we refer to such states as quasiglassy. They are thus defined as the states with a uniform distribution of phases y i (31), indistinguishable from incoherence but where, in contrast to the latter, some macroscopic portion of oscillators are frequency-locked (32), that is, have equal phase velocities _ y i . It should be noted that the condition for the absence of phaselocking is sometimes taken as R ¼ 0. However, although being implied by equation (31), it is not as rigorous as the latter. Thus, there might be many configurations for which R ¼ 0 but where the oscillators actually adjust their phases. Examples include two phase-locked populations of the same size, which are in the antiphase with each other. States with R ¼ 0 satisfying condition (32), but not condition (31), will be called spurious glassy, according to ref. 31 where they were first discovered. All SSs with R40, including the usual synchronized states, p-states 19,20 and travelling waves, will be classified as coherent.
As a simplified picture, the distinctions between the states can be understood in terms of a group of people doing cyclical exercises, each with their own tempo and other parameters. The incoherent state is when everyone proceeds independently; the coherent state is when they all move synchronously, at any given time having similar poses; spurious glassy is when they exercise at the same tempo but remain pairwise in the opposite poses, and quasi-glassy is when they adjust their tempos, but always remain in the independent random poses, thus representing a kind of synchronous disorder.
Realization of glassy states. Considering model (1) with the uncoupled distribution (23), it can be shown (see Methods) that if the marginal distribution of b is uniform and |J|40, then all SSs except incoherence satisfy equations (31) and (32). Therefore, the quasi-glassy states appear when Z hðq; b; gÞdqdg ¼ 1=2p; Representing h(q, b, g)h 1 (b)h 2 (q, g|b), equation (36) becomes h 1 (b) ¼ 1/2p. Thus, to satisfy also condition (37), the distribution of q, g should be specifically correlated with b.
The simplest examples are h(q, g|b) ¼ d(q À q 0 )P(g þ b,r) and h(q, g|b) ¼ L(q À q 0 cosb, D)P(g À g 0 , r). For spurious glassy states, the condition (36) is not satisfied, but R ¼ 0, which in the present case, is equivalent to |I|| R e ib h 1 (b)db| ¼ 0, as can be seen from equation (27). There are many possible h 1 (b) satisfying this condition, for example, Examples of the four types of states occurring in the model (1) are shown in Fig. 2, presented in the natural frame where the population does not move as a whole (/oS ¼ 0). In the (y,o)plane, the only difference between the quasi-glassy state and incoherence is that, in the former, phases within the glassy cluster (|o|o|k|W) are frozen around random angles b (static disorder), as seen in (y, b)-plane; this is in contrast to their asynchronous  P(b, 0.8). The pictures will be the same for any h 2 (q, g|b) satisfying |J| ¼ 1, f J ¼ 0 in each case. All states have O ¼ 0 and are presented in the natural frame (/oS ¼ 0); depending on the distribution G(C) of parameters in (1), similar states may also appear in a rotating frames (Oa0), corresponding to quasi-glassy, spurious glassy or coherent travelling waves. movement (dynamic disorder), observed for incoherence. Drawing an analogy between the oscillators' phases and the particles' positions in space, the incoherent state of the system (1) can be related to a liquid or gaseous state of matter; coherent to the crystalline state; spurious glassy to a crystal with a specific structure; and quasi-glassy to a glass.
As discussed previously, the model (1) can be reduced to a form without b (10) by a change of variables. Thus, any quasiglassy or spurious glassy state in terms of y i (1) can be transformed to a coherent state in terms ofỹ i ¼ y i À b i . Contrariwise, glassy states of the system (1) can always be viewed as coherent states but in a specifically disordered coordinate system, being therefore redundant in a strict mathematical sense. In practice, however, the coordinates are chosen based on their physical meaning, so that these states are realizable when the actual interaction between the oscillators contains phase shifts b, which can be achieved, for example, in a real systems with programmable coupling 48,49 . A similar situation occurs for the Mattis spin glass 50 : it can be transformed to a ferromagnetic state by change of variables 26,27 . Nonetheless, the Mattis model found a number of applications [51][52][53]  Relaxation dynamics. Up to now, we have concentrated mainly on the behaviour of the system (1) in the asymptotic limit t-N, restricting our consideration to a set of its possible SSs. The relaxation to these states is determined by the parameter distribution and can take various forms. It is hard to analyse in general, but for particular G(C) the full-time evolutions Y(t),Z(t) can be obtained analytically by using the OA-reduction procedure 41 , which can be extended to incorporate the distributed phase shifts. Thus, consider the distribution In this case, a(C,t) can be analytically continued inside the unit circle of z g ¼ e ig , as the condition (15) is satisfied there. Hence, in equations (8) and (9), one can integrate over g by changing the integration to the unit circle of z g and taking the residue at the pole z g ¼ r 0 e if 0 of P(g À f 0 ,r 0 ) (see equation (11)). The integration over o is performed in the usual way 41 , that is, by taking the residue at o ¼ iD. Then, from equations (8) and (9) one obtains ZðtÞ ¼ aðiD; e ig ¼ r 0 e if 0 ; tÞ, YðtÞ ¼ r 0 e if 0 ZðtÞ. Substituting this into equation (6) and solving the resultant equations, one gets finally where we have represented all in terms of I, J (see equation (24)), which in the present case are j I j Obviously, a similar procedure to that outlined here can be applied if the distribution of g has a few poles inside the unit circle. Now consider the KM (1) with distributed b: in which case, one has Z(t) ¼ Y(t) by definition (2). By the change of variablesỹ i ¼ y i À b i , the present case can be reduced to the already-studied KM with distributed g (see equation (38)). Then, because ZðtÞ ¼ YðtÞ ¼ e YðtÞ (see discussion of equation (10)), the macroscopic parameters should still evolve according to equation (39), although one should now use j I j e if I ¼j J j e if J ¼ r 0 e if 0 in the latter. The same result may be obtained by applying the OA-reduction, with the integrals (8) and (9) over b being evaluated by taking residues inside the unit circle of e ib , similarly to what was done in the previous case.
However, the behaviour predicted by equation (39), and the actual behaviour of Z(t),Y(t), agree only for distributed g, but not for distributed b, as demonstrated in Fig. 3. The reasons for this are, first, that one cannot continue a(C,t) inside the unit circle of b, as (15) is not satisfied there, which makes OA-reduction impossible when b is distributed. Second, the evolution of Y(t),Z(t) for (40) cannot be obtained from the system of y i ¼ y i À b i , at least if starting from the standard initial conditions (16), as assumed here. This is because such transformation introduces a correlation of the initial phases with system parameters. For example, if all y i (0) ¼ 0, then one will haveỹ i ð0Þ ¼ À b i ¼ Àe g i , so that e að e C; 0Þ ¼ e À ig and the OAreduction cannot be applied because now (14) is not satisfied. Therefore, although the cases (38) and (40) can be transformed into each other by changes of variables, the time evolution for distributed b appears to be much more complex than for distributed g, and cannot be obtained in a simple form.
Super-relaxation. Astonishingly, for a class of distributions G(C), and a very large family of initial configurations, the oscillators do not feel any interaction at all while relaxing to incoherence. This behaviour occurs when the following three conditions are satisfied: Z qe ig Gðo; k; q; b; gÞdkdqdbdg ¼ 0; 8o; ð41Þ Initial conditions rðy; t ¼ 0 j CÞ are not correlated with q or g; either directly or indirectly ðthrough the other parametersÞ; Incoherence is the only stable state of the system: Thus, conditions (41) and (42) imply W(0) ¼ R e iy r(y,0|C)qe ig G(C) dydC ¼ 0. Based on numerical evidence, the weighted mean field then stays at zero during the whole evolution W(t) ¼ 0, leading to an effective disappearance of interaction between the oscillators, as follows from equation (3). As a result, the oscillators evolve freely ( _ y i ¼ o i ) and so the relaxation (both its rate and form) depends This phenomenon, which we refer to as super-relaxation, is illustrated in Fig. 4. Note that for uncoupled distributions, equation (23), the condition (41) reduces to R qe ig h(q, b, g) dqdbdg ¼ 0.
The claim that W(t) ¼ 0 if all the conditions (41), (42) and (43) are satisfied can be proven rigorously for particular distributions, for example, G(C) ¼ G(o, k, q, b þ g)/2p (see Methods), though we were unable to prove it in general; based on simulations, however, it seems to be satisfied in most (if not all) cases. Thus, it is clear that the system always has a solution with W(t) ¼ 0 if conditions (41) and (42) are fulfilled: in this case, to have W(t)40, the distribution r(y, t|C) should acquire a particular dependence on the system parameters; but evolving according to equation (5) with W(t) ¼ 0, it can become dependent only on o, which in turn cannot increase W(t) above zero if condition (41) holds. The stability of this solution is rather hard to prove, but it appears to be stable if the incoherence is the only stable state (condition (43)). Nevertheless, even if W(t) is not exactly zero, below some relatively high threshold its value seems to have a negligible effect on the dynamics of R(t), as will be seen in Fig. 5 below.
Super-relaxation can appear only if the final state is incoherence (W ¼ 0), while relaxation to SSs with W40 will generally be coupling dependent, even if conditions (41) and (42) are satisfied. In the latter case, relaxation occurs in two stages as shown in Fig. 5 for the example of quasi-glassy states. The phases first begin to disorder in the same way as for incoherence, but are slowly entrained while passing their equilibrium positions. When the field of the entrained oscillators, characterized by W, becomes strong enough, they begin to force the unentrained ones to take their positions, so the relaxation switches to a faster, coupling-dependent regime. This switch occurs sooner for stronger coupling.

Discussion
We have generalized our earlier approach 10 to make it applicable to the more general KM (1) so that, using equations (19)(20)(21)(22), one can immediately obtain a macroscopic characteristics of possible SSs. The generalized KM (1) encompasses a variety of KM variants studied earlier 10,[13][14][15]19,21,24 , allowing one readily to reproduce and extend many of the previous results. Remarkably, the steady-state behaviour of (1) with any distribution h(q,b,g) (see equation (23)) can be obtained from the simple Sakaguchi-Kuramoto model (29), (30). It should be noted, however, that (17)(18)(19)(20)(21)(22) describe only SSs, being inapplicable to inherently non-stationary solutions such as standing waves 55 or oscillating p-states 44 . Note also that most expressions can straightforwardly be extended to the case G(C)-G(C,R,W) (see Methods); for examples when this might be relevant see refs 56-59. Most interestingly, we have found, that the model (1) can exhibit exotic behaviour, such as glassy states and superrelaxation, thereby opening new horizons for KM-related investigations, both theoretical and practical. These discoveries have a far-reaching implications. For example, it should now be possible to create, observe and study glassy behaviour in real systems of coupled oscillators, where a variety of novel phenomena may be anticipated. As one possible application, if some physical quantity can be associated with the weighted mean field W (equation (2)) in laser arrays 5  (R ¼ 0) but for which the other effects are nonvanishing (W40); similar considerations apply to a wide range of different KM applications. Furthermore, the phenomenon of super-relaxation might be used to design systems whose dynamics remains highly stable in the face of different perturbations and parameter changes.

Methods
Numerical simulations. Except where otherwise specified, all simulations were performed by integrating the full system equation (1)  Parameter normalization. Since the model (1) , there exist some ambiguity in the parameter definitions. To remove this ambiguity, one can fix the normalization of q i and specify a rule for choosing the sign of k i , q i ; but different choices might be preferred in different situations. For example, the normalization P q i ¼ 1 is inapplicable when q i ¼ cos(g i ) with g i being uniform; P |q i | ¼ 1, on the other hand, fails for a Lorenzian distribution of q i . The most universal choice seems to be R | R qe i(b þ g) G(C)dqdbdg|dodk ¼ 1, which additionally assures Wr1. However, when q i ¼ 1 and b i are distributed (so Y ¼ Z), this normalization will lead to a rescaling of q i (so YaZ), which might be inconvenient. Therefore, because the normalization can be fixed at any time, we retain the associated ambiguity in order to preserve generality. Note that W (equation (2)) but not R and not |k|W, changes under the (k,q)-rescalings and thus can be higher than unity.
Phase distribution. The distribution P(f, r), equation (11), represents a Poisson kernel for the unit disc and is very common for the KM. For example, it can be shown that the OA-ansatz (equation (6)) can be rewritten 43 as r(y, t|C) ¼ P(y À arg[a(C, t)], |a(C, t)|). Because P(f, r) has one pole inside and one outside the unit circle of z f e if , it is very convenient in relation to analytic derivations, for example, one has R p À p e if Pðf À f 0 ; rÞdf ¼ re if 0 . To simulate this distribution, one sets f i ¼ 2 arctan½ 1 À r 1 þ r tanðpðp i À 1=2ÞÞ, where p i A[0,1] are uniformly distributed random numbers.
Derivation of the incoherence stability conditions. To perform a linear stability analysis of incoherence, we examine the maximum growth rate for small perturbation eZ(y, t|C) added to the incoherent solution r(y, t|C) ¼ 1/2p. Following ref. 45, this perturbation is expressed as Zðy; t j CÞ ¼ cðCÞe À iy e lt þ c Ã ðCÞe iy e l Ã t þ Z ? ðy; t j CÞ; ð45Þ where the Fourier expansion of Z > (y, t|C) contains only terms e iny with |n|41 (so that R (e ±iy )Z > (y, t|C)dy ¼ 0). Substituting r(y,t|C) ¼ 1/2p þ eZ(y,t|C) and (45) into (5), and collecting the terms proportional to e À iy , to the first order in e one obtains Considering the part Be iy of equation (5), one will get the complex conjugate of equation (46), while for terms Be in y with |n|41 the growth rate of the associated perturbations q t [eZ > (y,t|C)] can be shown to be O(e 2 ), so that they are of no interest for linear analysis. From (46) it follows that c(C) ¼ (ke ib Y(t)e À lt )/(2e(l À io)). Substituting this into Y(t) ¼ R e iy r(y,t|C)qe ig G(C)dydC ¼ ee lt R qe ig c(C)G(C)dC and dividing both sides by Y(t)/2, one self-consistently obtains where we have used the definition of J(o,k) (see equation (12)). The incoherence changes stability when Re[l] crosses zero. Hence, denoting l ¼ l r þ il i , taking the limit l r -0 in equation (47) and using lim lr!0 ½l r À iðo À l i Þ À 1 ¼ pdðo À l i Þ þ i½o À l i À 1 , one obtains (21) with x ¼ l i .
Derivation of the ESCs. The ESCs (22) cannot be derived rigorously but, rather, are based on empirical assumptions. Thus, we first assume that the perturbations to the SS obey where dY ¼ (dW þ iWdF)e iF and dO are the deviations of the SS parameters from their stationary values Y and O, respectively; A is a constant, G þ O is defined in (11) and a s (C,Y) is the OA stationary solution (17). The latter can be represented in unified form for all o, k as Although equation (48) represents a purely empirical assumption, some of the motivation behind it is that for incoherence it gives the correct transition points (equation (21)), as can be checked by performing a linear stability analysis of equation (48) with Y ¼ 0.
For Y40, based on trial and error we set A40, dO ¼ W 2 dF in equation (48). Then using equation (49) and retaining in (48) only the terms of the first order in dW, dF, one obtains  (50) is stable only if the trace and determinant of the corresponding matrix are lower and higher than zero, respectively, which gives (22).
Glassy conditions for the model considered. For the system equation (1) with the uncoupled distribution, equation (23), it is easy to see that the conditions (31,32) for observing quasi-glassy states can be rewritten as (36,37). Thus, the phases of the oscillators with |o i |r|k i |W in a stationary regime are frozen at in the rotating frame), as follows from equation (18). Therefore, for states with W40, an absence of phase locking between the oscillators (condition (31)) is equivalent to a uniform distribution of b i (equation (36)). Next, W40 (37) is the necessary and sufficient condition for (32): necessary because, otherwise, all oscillators have _ y i ðtÞ ¼ o i ¼ 0, as implied by equation (3); and sufficient because it establishes frequency locking of the oscillators with |o i |r|k i |W (since their phases are constant). Finally, as is clear from equation (25), for the uncoupled distributions (23), SSs with W40 can appear only if |J|40, as additionally indicated in equation (37).
Glassy states exhibiting nearly all glassy properties. One can modify the model (1) to observe glassy states with most of the desired properties (31)(32)(33)(34)(35), and not only the first two of them. Thus, frustration (34) can be introduced, for example, by considering the van Hemmen type of interactions 32,33 , while the high multiplicity of states (33) can be achieved by using the higher-order coupling 22,54 . As an example, it can easily be checked numerically that a state satisfying conditions (31-34) appears in the system _ y i ¼ o i À 10 N X j ðk i q j þ q i k j Þ sin½3ðy i À y j À b i À g j Þ; ð51Þ where k i and q i are independent random variables taking values ±1 with equal probability, b i ¼ À g i are uniformly distributed and the frequencies o i are drawn from g (o)B[1 þ o 2 ] À 1 . Regarding the last glassy property (35), it might be possible to satisfy it with an appropriate choice of the parameter distribution, which, to a large extent, determines the relaxation of R(t). However, most interestingly, it turns out that for systems of the type (51), there exist many configurations with R(t-N)40, and the latter is proportional to the initial R(0). Thus, in many cases the phases do not disorder completely, so the basic criterion (31) becomes violated if starting from initial conditions other than incoherence.
Super-relaxation claim W(t) ¼ 0. In some cases, one can rigorously prove that W(t) remains zero at all times if the conditions (41), (42) and (43) are fulfilled. Thus, consider the distribution Changing the variables toỹ i ¼ y i À b i , one obtains the system (10) with e Gð e CÞ ¼ Gðo; k; q; e gÞ, characterized by the same weighted mean field e WðtÞ ¼ WðtÞ (see discussion of (10)). Under such a change, however, one subtracts from each y i a uniform random variable b i , so that any initial conditions for y i satisfying (42) will be mapped to a uniform initial condition forỹ i ; moreover, since G(C) depends only on the sum g þ b, these initialỹ i will not be correlated with e g i ¼ g i þ b i , as would be otherwise. Therefore, one has e rðỹ i ; 0 j e CÞ ¼ 1=2p ) e að e C; 0Þ ¼ 0. Then, taking into account that e Wð0Þ ¼ Wð0Þ ¼ 0, it follows from equations (7) and (8) that e að e C; tÞ ¼ 0 ) e WðtÞ ¼ WðtÞ ¼ 0. Finally, because the incoherent state for y i is stable (condition (43)), it will obviously be stable in terms ofỹ i as well, so that the perturbations to the weighted mean field will decay, ensuring that it stays at zero. from which the SS parameters are determined, will also be the same except Note, however, that if the parameter distribution depends on R (q R G(C,R,W)a0), and q or g are distributed (so that WaR), then equations (19) and (20) become coupled and should be solved at the same time to get all the R,W,O,F À C. Such situation occurs because, in this case, the effective interaction between the oscillators, although explicitly determined by W (3), has additional implicit dependence on R; as a result, the dynamics of the model becomes more sophisticated, and its analysis appears to be an interesting topic for future research.
Assuming that lim R!0;W!0 R n W m @ n R @ m W GðG; R; WÞ ¼ 0; 8n; m ! 0; n þ m40, the conditions for incoherence stability can be derived in the usual way and take the same form (21), except that J(o, k)-J(o, k, 0, 0) (53). The ESCs (22), on the other hand, cannot so easily be generalized to field-dependent distributions, at least if q R G(C, R, W)a0 and WaR; otherwise, that is, when G(C)-G(C, W), they preserve the original form of equation (22), but with a modified J given by equation (53) (though it is not clear how well they work in this case). Having obtained the general equations, their simplification for the uncoupled distributions G(C, R, W) ¼ g(o, k, R, W)h(q, b, g, R, W) to analogues of equations (24)(25)(26)(27)(28) is straightforward.