{'fraction_non_alphanumeric': 0.06692610988785908, 'fraction_numerical': 0.014849710686671104, 'mean_word_length': 4.14760147601476, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Long Time Evolution of Phase Oscillator Systems

It is shown, under weak conditions, that the dynamical evolution of an important class of large systems of globally coupled, heterogeneous frequency, phase oscillators is, in an appropriate physical sense, time-asymptotically attracted toward a reduced manifold of system states. This manifold, which is invariant under the system evolution, was previously known and used to facilitate the discovery of attractors and bifurcations of such systems. The result of this paper establishes that attractors for the order parameter dynamics obtained by restriction to this reduced manifold are, in fact, the only such attractors of the full system. Thus all long time dynamical behavior of the order parameters of these systems can be obtained by restriction to the reduced manifold.

Systems consisting of many coupled phase oscillators have been used to model diverse situations ranging from Josephson junction circuits, to circadian rhythms, to synchronization of cardiac pacemaker cells. In previous work by us [15], it was shown that a large class of such models possess solutions on an invariant manifold M. It has since proved possible to simply obtain various attractors of the dynamics on M. A remaining open question is that of whether such attractors for the dynamics on M are also attractors for the dynamics of the full system, and, if so, whether all of the attractors of the full system lie on M. In this paper we prove, under very general conditions, that, in an appropriate sense, the answer to these questions is yes. This result establishes that restriction of consideration to the manifold M can be used as an effective computational and analysis method for obtaining all the typical, long-time dynamical behavior of these systems.

I. INTRODUCTION
Large systems of coupled phase oscillators with heterogeneous frequency distributions are of general interest and are the essential modeling tool in past analyses of a variety of interesting situations in physics, chemistry, biology, etc. Perhaps the simplest and best known such system is the Kuramoto model [1], which treats the synchronization of globally (all-to-all) coupled phase oscillators for which the coupling between pairs of oscillators appears as the sine of the phase difference between the oscillators. Examples where this basic framework has been extended to more complex situations include Josephson junction circuits [2], pedestrian induced oscillation of footbridges [3,4], systems with time-dependent coupling [5], driven systems describing circadian rhythm in mammals [6,7], the effect of time-delay in oscillator interactions [8,9], the effect of non-unimodal distribution of the natural frequencies of the phase oscillators [10,11], "communities" of phase oscillators interacting with multiple other phase oscillator communities [12,13], the "chimera" model of certain mammals that experience sleep with only one of their two brain hemispheres at a time [13,14], etc.
The large number of interesting applications of phase oscillator models motivates the attempt to find general analysis tools applicable to these models. In this vein, it has recently been shown that, in the continuum limit (i.e., the number of oscillators approaches infinity), such models possess solutions on a reduced manifold of system states [15]. Furthermore, for the case of a Lorentzian distribution of oscillator frequencies, the dynamics on the reduced manifold is typically describable by a finite number of ordinary differential equations. This finding has been utilized to determine attractors and their bifurcations on the reduced manifold for all the applications previously mentioned (see Refs. [2,3,5,7,9,11,13,14]). Two basic questions remain: (i) are attractors for the dynamics restricted to the reduced manifold also attractors of the full system; and (ii) are there attractors of the full system that do not lie on the reduced manifold? Indications of results so far are mixed. On the one hand, numerical results from Refs. [6,9], and especially [11], are consistent with the supposition that all attractors of the full system lie on the reduced manifold. On the other hand, Refs. [13,16] find long-time asymptotic behavior that is not on the reduced manifold. The result of our paper is that, in an appropriate sense (that we specify later in this paper), the reduced manifold is globally attracting provided that the spread ∆ in the distribution of oscillator frequencies is nonzero. In particular, for ∆ > 0 all attractors of the full system lie on the reduced manifold, and all attractors of the dynamics on the reduced manifold are attractors of the full system. This greatly facilitates the task of finding the attractors of the full system, since they now can be sought using the reduced system. The result also resolves the puzzle posed by the previous results, since the finding by Refs. [13,16] of long time motion not on the reduced manifold was for the case of ∆ = 0, while the opposite indication from the numerical results of Refs. [6,9,11] treated situations in which ∆ > 0.

II. FORMULATION
We begin by noting that the models in the class of problems in which we are interested all involve the determination of a distribution function F (θ, ω, t), where θ is the phase of an oscillator, and ω is the natural frequency an oscillator would have in isolation from the outside world (e.g., from other oscillators); F dθdω is the fraction of oscillators at time t whose phases and natural frequencies lie in the range [θ, θ + dθ] and [ω, ω + dω]. Since the natural frequency ω of an oscillator is assumed not to change with time, the marginal frequency distribution, is time independent. The key quantity characterizing the macroscopic behavior of the distribution function F is the "order parameter" r(t) originally introduced by Kuramoto [1] and defined by Since the number of oscillators is conserved, F obeys an oscillator continuity equation, and, for all of the problems previously mentioned (Refs. [2] to [14]), v θ (θ, t) is expressible in the form [2], Equations (3) and (4) constitute an ω-dependent partial differential equation in the two real variables (θ, t) to be solved subject to an initial condition F (θ, ω, 0). The problem is ostensibly complicated by the fact that the time dependence of the quantity H may, in general (Refs. [1] to [15]), depend on F through the complex order parameter r(t) defined by (2), as well as through other, non-phase-oscillator variables obeying auxiliary dynamical equations, which themselves may [2,3] depend on r(t), or (as in Ref. [5]) through explicit external time dependence of system parameters. Here are some examples: for the classical Kuramoto [1] problem, H = kr(t), where k is the strength of the coupling between oscillators; for the circadian rhythm problem [6,7], H = kr(t)+Λ, where Λ is a constant reflecting the strength of the diurnal drive of the day-night sunlight cycle, and Λ might be given an explicit time dependence, Λ = Λ(t), to model variation between sunny and cloudy days; for the case of time delay in the response of oscillators to other oscillators in the system [8,9], is the distribution function [9] of delays along the links between oscillators; in the cases treated in Refs. [12,13] (communities of oscillators), [11] (nonunimodal frequency distribution g(ω)), and [13,14] (the chimera model), there are several distribution functions, i.e., F and H in (3) and (4) are replaced by F σ and H σ (σ = 1, 2, · · · , s, where s is the number of distributions) and each H σ is a function of all the order parameters r 1 , r 2 , · · · , r s ; in the case of pedestrian induced oscillation of footbridges [3,4], H = kÿ(t), whereÿ(t) is the side-to-side acceleration of the bridge, which obeys a damped oscillator equation driven by r(t), where r(t) represents the effects of the pedestrians.
For a general problem of the type described above, as the system evolves, H(t) will change self consistently with t. We will show in what follows that, whatever is the evolution of H(t), in the long time limit, solutions for the order parameter evolution r(t) (or r σ (t)) obey the differential equation that applies for evolution on the reduced manifold of Ref. [15]. Because the precise time dependence of H(t) will not matter in the derivation of our result, it suffices to consider H(t) as some general function of time without regard to how this time dependence is determined.
Expanding the distribution F as a Fourier series in θ, we write F in the form, We note that 2π 0 F ± dθ = 0 and that, assuming absolute convergence of the Fourier series, the analytic continuation of F + (F − ) into Im(θ) > 0 (Im(θ) < 0) has no singularities and decays exponentially to zero as Im(θ) → +∞ (Im(θ) → −∞). As will soon become evident, the decomposition of F given by (5) is a key step. We note that since F − = F * + for real (θ, ω), it suffices to consider only F + . Substituting (5) into Eqs. (3,4) and projecting the result onto the function space spanned by the basis functions, e iθ , e 2iθ , e 3iθ , · · · , we obtain As previously noted, our result will not depend on the precise time dependence of H(t). Thus, whatever is the dependence of H(t), we can formally regard it as given. Adopting this viewpoint, Eq.(7) is linear in F + with an inhomogeneous driving term on the right hand side (namely, 1 2 H * exp(iθ)). As such, we can write F + as whereF + is a homogeneous solution to (7) andF + ′ is an inhomogeneous solution. An inhomogeneous solution is given by taking the Fourier coefficients ofF + ′ to be given byF n ′ (ω, t) = [α(ω, t)] n , as proposed in Ref. [15]. When this ansatz is used in (7), it is found that (7) is indeed satisfied if α(ω, t) satisfies which, for each value of ω, is an ordinary differential equation in time t. We note further that, as shown in [15], |α(ω, t)| < 1 so that the summation of the Fourier series forF + ′ converges and yieldŝ This form ofF + ′ specifies the reduced manifold found in Ref. [15]. Thus, at long time, F + would tend toF + ′ (i.e., F would tend to the reduced manifold) if lim t→∞F+ = 0. However, a simple counter examples shows that this cannot always be true. In particular, if H = 0, Eq. (7) has homogeneous solutions, for any set {A n } for which this series converges. Since the magnitude of each term of the summation in (11) is time-independent (ω is real),F + does not go to zero at t → ∞. On the other hand, we note that, as t increases, the individual terms, e in(θ−ωt) , oscillate more and more rapidly in ω. Thus for any such term A n e in(θ−ωt) dω decays exponentially in time for sufficiently smooth g(ω). For example, our subsequent considerations will be for the case of a Lorentzian frequency distribution, for which which decays exponentially to zero as t → ∞ provided that ∆ > 0 (the case ∆ = 0 corresponds to g(ω) being a delta function in ω). [We remark that the mean value of ω has been taken to be zero in (12), but that no generality is lost by this, as a mean value can be restored by the change of variables, θ ′ = θ − Ωt, ω ′ = ω + Ω.] Thus, while we cannot expect to show thatF + → 0 as t → ∞, Eq. (13) suggests that this may not be necessary to obtain order parameter dynamics that tend to the order parameter dynamics that applies on the reduced manifold. In particular, noting from (2), (5) and (6) that what we require is that where In what follows we will demonstrate that Eq.(15) indeed holds under very general conditions.

III. DEMONSTRATION OF THE MAIN RESULT
To show that (15) applies, we now assume that the analytic continuation ofF + (θ, ω, t) into Im(ω) < 0 has no singularities in Im(ω) < 0 and approaches zero as Im(ω) → −∞. To show that this last assumption is a consistent one, let |ω| be very large, |ω| ≫ H. Then the homogeneous version of Eq.(7) forF + is approximately which has solutions for its Fourier θ-componentsF n ∼ exp[in(θ − ωt)] which go to zero as Im(ω) → −∞. (It was to achieve this that we have introduced the decomposition of F given by Eqs. (5) and (6).) We will further discuss this analyticity assumption at the end of this paper.
We now specialize to the case of Lorentzian g(ω), Eq. (12). We multiply the homogeneous version of Eq.(7) by g(ω)dω, integrate the result from ω = −R to ω = +R, analytically continue into the complex ω-plane, close the integration path with a semicircle of radius R in the lower half ω-plane, and let R → ∞. Using our assumption that F + (θ, ω, t) is analytic in the lower half ω-plane and decays to zero as Im(ω) → −∞, the integral along the large semicircle approaches zero as R → ∞, and the integrals from ω = −∞ to ω = +∞ along the real ω-axis may thus be evaluated as the residue of the enclosed pole of g(ω) at ω = −i∆ (see Eq. (12)). This yields where f + (θ, t) =F + (θ, −i∆, t).
We now introduce a conformal transformation of the upper half complex θ-plane into the unit disc, z = e iθ . Equations (17) and (18) then become whereṽ (z, t) = ∆z + 1 2 Noting that (19) can be written as where d/dt = ∂/∂t +ṽ∂/∂z, we can integrate (22) along the characteristics of this equation to obtaiñ where dt ′ , and the characteristics are given by the orbit equation, with the final condition Z(z, t) = z. Thus Z(z, t ′ ) for t ′ < t represents the location Z of the orbitṽ(z ′ , t ′ ) that winds up at point z at time t. It is useful to rewrite (24) by introducing with h, β, ρ and φ real. The real and imaginary parts of (24) then give We note from (25) that when ρ = 1, we have dρ/dt ′ = ∆ > 0. Thus for final conditions on ρ = 1, the orbits backward in time move into ρ < 1. Thus |Z(z, t ′ )| < 1 for |z| ≤ 1 and t ′ < t (i.e., Z(z, t ′ ) is in the unit disc). We wish to show thatf + (z, t) → 0 as t → +∞. From (23) we see that this will be the case if In order to show this, we first note that the real part of ∂ṽ(z ′ , t ′ )/∂z ′ is simply one half the divergence of the two dimensional flowṽ =ṽ ρ (ρ, φ)ρ 0 +ṽ φ (ρ, φ)φ 0 (whereṽ ρ andṽ φ are given by (25) and (26), and ρ 0 and φ 0 are unit vectors in the ρ and φ directions), i.e., Equation (27) is most easily demonstrated in rectangular co-ordinates: z ′ = x + iy,ṽ(z ′ , t ′ ) =ṽ x (x, y, t ′ ) + iṽ y (x, y, t ′ ), whereṽ x ,ṽ y , x and y are real. Then (27) immediately follows by settingṽ =ṽ x x 0 +ṽ y y 0 , and using the Cauchy-Riemann condition, ∂ṽ x /∂x = ∂ṽ y /∂y, in the expression for the divergence in rectangular coordinates, ∇ ·ṽ = ∂ṽ x /∂x + ∂ṽ y /∂y. Now evaluating ∇ ·ṽ in polar coordinates (ρ, φ), we have Solving (28) for h cos(φ−β) in terms of ∇·ṽ and inserting the result in Eq.(25) for dρ/dt ′ , we obtain after some rearrangement Inserting (29) into the integral for η(z, t) and choosing a fixed reference time T satisfying 0 < T < t, we have We are interested in final (t ′ = t) conditions on the unit circle, Z(z, t) = z = e iθ for θ real, and their continuation into the unit disc |z| ≤ 1, corresponding to ρ ≤ 1 at the final time t ′ = t. For ρ sufficiently near one, Eq.(25) shows that dρ/dt ′ ∼ = ∆. Thus by the continuity of the right hand side of (25), there is an annulus in the Z-plane, 1 ≥ ρ ≥ ρ 0 , in which dρ/dt ′ > 0, implying that as t ′ is reduced from t (i.e., t−t ′ is increased), ρ moves uniformly from ρ = 1 at time t to smaller values. Thus any final point in the annulus eventually enters the disc ρ < ρ 0 < 1 and never leaves it. We can therefore choose the time T such that for all final conditions |Z(z, t)| = |z| ≤ 1, we have where the first inequality applies for 0 ≤ t ′ < t − T . We consider T to be held fixed, and we ask how η(z, t) behaves as t → +∞. By (28) the integrand in the first of the three terms of (30) is bounded, and, since the integration range in t ′ (namely, T ) for this term is fixed, we conclude that the first term is bounded. By (31) the second term in (30) is also bounded. Again by (31) the integrand of the third term of (30) is positive and greater than one. Thus the third term exceeds (t − T )∆. Hence Re[η(z, t)] → +∞ for t → +∞, if ∆ > 0, thus demonstrating that the exponential factor in (23) goes to zero for large time [17]. We therefore conclude from (30) and (23) that Eq. (15) is satisfied if ∆ > 0, which is the desired result.
In particular, at large t the order parameter will approach the quantity ∞ −∞ α(ω, t)g(ω)dω, where α(ω, t) evolves by Eq. (9). This implies that r(t) will satisfy the differential equation, which follows from multiplying (9) by g(ω)dω, integrating from ω = −∞ to ω = +∞ and, as done previously, using the residue method to evaluate the integrals. Hence, for ∆ > 0, the long time dynamics of the order parameter r(t) is governed by the ordinary differential equation (Eq.(32) and Ref. [15]) that describes its dynamics for distribution functions F on the reduced manifold. This is our main result.

IV. DISCUSSION
The principal conditions for the applicability of our result are that the initial condition is such that, when F + (θ, ω, t) is continued into the complex ω-plane, it is analytic in Im(ω) < 0 and decays to zero as Im(ω) → −∞. As discussed in Ref. [15], if these conditions are satisfied initially, then they are also satisfied for all t > 0. What happens if the condition at Im(ω) → −∞ is not satisfied initially? Here, a simple example [18] may be instructive. We again consider the case H = 0. Say the initial condition on F + has a component exp(inθ + iγω) with γ real and positive. This initial condition violates our assumption of decay to zero as Im(ω) → −∞. However, use of this initial condition in Eqs. (3,4) with H = 0 yields the solution exp(inθ + i(γ − t)ω), and, for large enough time, t > γ, the result satisfies the required condition that it approaches zero as Im(ω) → −∞ [19]. Thus, even if our desired condition at Im(ω) → −∞ is not satisfied initially, in many cases, the result that the long time dynamics of r(t) is described by Eq.(32) may still apply.
We now connect our result with the concept of an inertial manifold. An inertial manifold M with respect to a distance metric µ satisfies the condition that, for any initial condition in the state space, the subsequent system evolution is such that the distance between the evolved orbit and the manifold M as measured by the metric µ approaches zero as t → +∞. What we have shown in this paper is that, in the space of distribution functions F (θ, ω, t), our reduced manifold (Eq.(10)) is inertial with respect to the proper distance metric µ. In particular, this is so if we take the distance between two distribution functions F 1 (θ, ω, t) and F 2 (θ, ω, t) to be defined by (33) (For F 1 not on the reduced manifold M , the distance from F 1 to M is µ(F 1 , F 2 ) minimized over all F 2 on M .) Note that, by this choice of distance metric, the problem associated with the example of Eq.(11) is avoided.
Finally, we remark that, while our result is for the special case of a Lorentzian distribution of oscillator frequencies (Eq.(12)), we believe that this restriction does not greatly limit the usefulness of the resulting formulations for discovering typical system behavior. Indeed, past numerical experiments [6,11] comparing results obtained using Lorentzian g(ω) and using Gaussian g(ω) were found to yield qualitatively identical bifurcation structures.
We are grateful to B.R. Hunt for an extremely useful comment. We also thank W.S. Lee and S.H. Strogatz for comments on a preliminary draft of this paper. This work was supported by the ONR(N00014-07-0734) and by NSF (PHY 0456249).