{'fraction_non_alphanumeric': 0.0437796652510617, 'fraction_numerical': 0.009638604380048297, 'mean_word_length': 4.498740842490842, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ott-Antonsen attractiveness for parameter-dependent oscillatory networks

The Ott-Antonsen (OA) ansatz [Chaos 18, 037113 (2008), Chaos 19, 023117 (2009)] has been widely used to describe large systems of coupled phase oscillators. If the coupling is sinusoidal and if the phase dynamics does not depend on the specific oscillator, then the macroscopic behavior of the systems can be fully described by a low-dimensional dynamics. Does the corresponding manifold remain attractive when introducing an intrinsic dependence between an oscillator's phase and its dynamics by additional, oscillator specific parameters? To answer this we extended the OA ansatz and proved that parameter-dependent oscillatory systems converge to the OA manifold given certain conditions. Our proof confirms recent numerical findings that already hinted at this convergence. Furthermore we offer a thorough mathematical underpinning for networks of so-called theta neurons, where the OA ansatz has just been applied. In a final step we extend our proof by allowing for time-dependent and multi-dimensional parameters as well as for network topologies other than global coupling. This renders the OA ansatz an excellent starting point for the analysis of a broad class of realistic settings.

The Ott-Antonsen (OA) ansatz [ Chaos 18, 037113 (2008), Chaos 19, 023117 (2009)] has been widely used to describe large systems of coupled phase oscillators. If the coupling is sinusoidal and if the phase dynamics does not depend on the specific oscillator, then the macroscopic behavior of the systems can be fully described by a low-dimensional dynamics. Does the corresponding manifold remain attractive when introducing an intrinsic dependence between an oscillator's phase and its dynamics by additional, oscillator specific parameters? To answer this we extended the OA ansatz and proved that parameter-dependent oscillatory systems converge to the OA manifold given certain conditions. Our proof confirms recent numerical findings that already hinted at this convergence. Furthermore we offer a thorough mathematical underpinning for networks of so-called theta neurons, where the OA ansatz has just been applied. In a final step we extend our proof by allowing for time-dependent and multi-dimensional parameters as well as for network topologies other than global coupling. This renders the OA ansatz an excellent starting point for the analysis of a broad class of realistic settings. Coupled phase oscillators are being widely used to describe synchronization phenomena. The study of their collective dynamics has experienced a major breakthrough by the results by Ott and Antonsen [1][2][3] . The asymptotic behavior of the mean field of infinitely many coupled oscillators can be cast into a reduced, low-dimensional system of ordinary differential equations. The evolution is hence captured by the so-called Ott-Antonsen (OA) manifold.
Very recently, the OA ansatz has been applied to networks of theta neurons, see, e.g., Refs. 4-10. A particular property of coupled, inhomogeneous theta neurons is that both the phase of a single neuron as well as its dynamics depend on a parameter, which establishes an intrinsic relation between them. While numerical results suggest the attractiveness of the OA manifold in the presence of such a parameter dependence, it has as to yet not been proven whether the dynamics really converges to it. For a certain class of parameter dependencies we here extend the existing theory of the OA ansatz and show that the OA manifold continues to asymptotically attract the mean field dynamics.
Parameter-dependent systems and their description through the OA ansatz have been considered by, e.g., Strogatz and co-workers 11 , Wagemaker and coworkers 12 , and So and Barreto 13 . There, parameters seemingly did not yield a correlation between an oscillator's phase and its dynamics but a rigorous proof for this is still missing. We explicitly address this last point. In particular, we prove a conjecture later formulated by Montbrió and co-workers 7 on the attractiveness of the OA manifold for parameter-dependent systems. The case of parameters serving as mere auxiliary variables readily follows from our result -we will refer to this as "weak" parameter-dependence 14 . By showing that a network of theta neurons can be treated as a parameter-dependent oscillatory system, our result establishes an immediate link to networks of quadratic integrate-and-fire (QIF) neurons: That is, the so-called Lorentzian ansatz as an equivalent approach to the OA ansatz is analytically substantiated. By this we may exert an important impact in mathematical neuroscience.
Finally, we extend the parameter-dependence for more general classes of networks. First, we address non-autonomous systems and show that our proof can be applied to time-varying parameters. An important example here is a biologi-cally realistic approach to oscillatory systems proposed by Winfree 15 . Second, we include multiple distributed parameters illustrated by coupled limit-cycle oscillators with shear. Third, we apply our proof to networks with different coupling topologies including non-local coupling by using an heterogeneous mean field approach.

I. INTRODUCTION
The Kuramoto model can be considered the most seminal description of globally coupled networks of phase oscillators. It has been investigated in great detail but its various extensions still make it the model-to-work-with when it comes to the study of network dynamics 16,17 . We adopt the notion of Montbrió, Pazó, and Roxin 7 and write the Kuramotolike model asθ where the phase dynamics of the j-th oscillator (j = 1, . . . , N) depends on its natural frequency ω j and a driving complex-valued field H. The latter can depend on time t, on the mean field z(t) = N j=1 e iθ j (t) , and on other auxiliary variables, but not on the (index of) oscillator, i.e. it remains identical for all oscillators j = 1, . . . , N. Given the right-hand side of (1), the oscillators are sinusoidally coupled.
In the thermodynamic limit (N → ∞) the OA ansatz yields solutions for the dynamical evolution of the corresponding distribution function (of all the oscillators), which are attracted towards a reduced manifold of states 1,2 . Central to this is the description of the system via its distribution density ρ(θ, ω, t). The quantity ρ(θ, ω, t) dθ dω is the fraction of oscillators whose phases are in the range [θ, θ+dθ] and have natural frequencies in [ω, ω+dω] at time t. The distribution function ρ obeys the continuity equation with velocity field The latter can equivalently be written as 11,12 In agreement with the assumptions on H we require that the functions f and h may explicitly depend on time t, on the (now continuum form of the) mean field z(t) = ∞ −∞ 2π 0 ρe iθ dθdω, and on other auxiliary variables, but not on the the phase θ itself.
Asymptotic attractiveness of the OA manifold, given by distribution functions of the form that satisfy the normalization condition has been proven for continuous frequency distribution functions g(ω) of non-zero width and for H being independent of θ; c.c. stands for complex conjugate. Other requirements include |α(ω, t)| ≤ 1, and some analytic continuity conditions. 1,2 In what follows we extend this approach by rigorously proving the asymptotic attractiveness of the OA manifold in the case of H and ω depending on an additional parameter η that may also influence θ. Equivalently, we include a time-and η-dependence of f and h in (3b).
By this, we allow for an intrinsic relation between θ, H, and ω, or θ, f , and h, respectively.
As of today, the attractiveness of the OA manifold in the (time-and) parameter-dependent case has only been hypothesized 11,18 but not proven.

II. PARAMETER-DEPENDENT SYSTEMS
When including additional parameters at the oscillator level, the dynamics (1) becomeṡ The natural frequency Ω of oscillator j may therefore deviate from ω j , which promotes further heterogeneity among oscillators. Moreover the driving field H may depend on η j .
The right-hand side of (6) expresses a certain dependence on the (index of the) j-th oscillator.
Hence, such a dependence is no longer exclusive to the sinusoidal coupling, but also affects the natural frequency Ω(ω j , η j ) and the driving field H(η j , t) .
When considering η a random variable, we may regard η j to be drawn from a distribution function g(η). Likewise ω j may be drawn from a (different) distribution function. The oscillator-specific parameter η j may change this distribution function in the oscillator's favor.
Therefore, we here incorporate a joint distribution g(ω, η) in the normalization condition (5). In general, ω and η are not independent and the joint distribution consists of two nested distributions. We hence replace Ω(ω j , η j ) by ω(η j ). Then, in the continuum limit (6) reads: The relation through η becomes now even more evident as the temporal derivative of θ has become partial.
Again, one can introduce a distribution function ρ(θ, ω, η, t), which now additionally depends on η. And again, this distribution function satisfies the continuity equation (2) with velocity field (7). In line with the parameter-independent case, in which the distribution function g(ω) of the natural frequencies ω had non-zero width 1,2 , we assume that the distribution function g(η) of the parameter η also has non-zero width. The frequency ω, thus, cannot be constant but depends on η. Likewise, the driving field H depends on η. Importantly, these two terms exhibit so an implicit dependence on θ, such that the proof for the attractiveness of the OA manifold as has been derived in Ref. 2 may no longer hold.
However, there is strong numerical incentive that the OA manifold fully covers the long-term behavior of the dynamics of the population of parameter-dependent phase oscillators; see, e.g., Refs. 4-10, 18-23.
In the following we demonstrate the proof of this conjecture for a particular class of parameter-dependent systems. We consider η to follow a Lorentzian distribution and assume that ω depends linearly on η, i.e. ω(η, t) = a · η + c, where, without loss of generality, we set a = 1 and consider c = c(t) ∈ L 1 (R) an integrable, and in particular piecewise smooth, function. Our line of argument follows closely that of Ott and Antonsen 2 but we extend their results whenever necessary. We would like to note that our findings remain valid for a larger class of distribution functions as has been depicted in detail in Ref. 3. We will comment on this and consider more general η-dependencies of ω in Sections IV and V.
Ifρ + is analytic in the lower half η-plane and decays to zero as Im(η) → −∞, one can multiply (9) by g(η)dη and integrate the result by employing the residue theorem. Hence, the integrals can be evaluated at the residue of the enclosed pole of g(η) at η = −i∆. We find ∂ tρ+ (θ, −i∆, t) + ∂ θ −i∆ ·ρ + (θ, −i∆, t) The two remaining integrals can be determined provided that H and H * have no singularities in the lower half η-plane and do not increase "too" fast for Im(η) → −∞. Since g is a Schwartz function, we only need H to diverge at most sub-exponentially. For common choices of H, as listed in Ref. 2, these requirements are met indeed, which yields Here we substituted f + (θ, t) =ρ + (θ, −i∆, t) and H(t) = H(−i∆, t). These equations agree exactly with Eqs. (17 & 18) in Ref. 2. Hence, following the same reasoning around Eqs. (19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30)(31) in Ref. 2 one can conclude that (11) is fulfilled. To underscore the line of argument, we would like to give a short sketch of the proof. First, by introducing a conformal transformation of the upper half complex θ-plane into the unit disc via w = e iθ , one can rewrite (12 & 13) as wheref + andṽ are the transformed functions from (12 & 13), and d/dt = ∂/∂t +ṽ∂/∂w. (14) can be integrated using the method of characteristics for linear and homogeneous partial differential equations 25 . Here we requiref + ∈ C 2 (R) butṽ does not need to be continuous.

This yieldsf
as solution with and the characteristics are given by with final condition W (w, t) = w. Finally, in order to show thatf + (w, t) → 0 for t → ∞, which, by (15), we prove that The details for the rather lengthy computation can be found in Ref. 2. We here we would only like to mention that the integral in (16) is split into three distinct parts, each of which is evaluated and while two of them remain bounded, the third diverges at the rate ∆t, presuming ∆ > 0. This eventually completes the proof and underlines the importance that the distribution function g(η) must have non-zero width ∆. We would also like to note that in the final step of the proof the continuity of v is required, i.e. H in (13) must be continuous.
If one includes, e.g., square functions in the time-dependent parts of the frequency term and/or driving field, one is confronted with jump discontinuities, which become present in the right-hand side of (13) either directly or indirectly via the order parameter z(t). A closer look at Ref. 2, however, confirms that for small jumps the reasoning can be guaranteed and for proper choices of a time constant T their Eq.(31) holds. Thus, we can argue that OA attractiveness will be maintained even in the case of discontinuities, which also confirms our rather long assumption for c(t) to be in L 1 (R) in the linear dependence of ω(η) = aη + c.
So far we only considered a Lorentzian distribution and some linear dependence of ω on η. However, our result can be extended to a much broader class of distribution functions g(η), non-linear dependencies ω(η), or even joint distributions g(ω, η) in the case of Ω(ω, η); see Section IV below. Hence, it is proper to say that the asymptotic attractiveness of the OA manifold for parameter-dependent systems of coupled phase oscillators is generic. Note that the proof remains identical if θ = θ(t) does not depend on the parameter η, that is, when there is no correlation between specific oscillators and their dynamics. We call this case "weak" parameter-dependence, which has been coined in several earlier studies, e.g., Refs. 11-13, and 18, where parameters were introduced as auxiliary variables. Our result therefore confirms the attractiveness of the OA manifold also in this case, as has simplifyingly been taken for granted in the afore-cited studies.

III. NETWORKS OF QIF AND THETA NEURONS
As mentioned above, there is a variety of recent papers that showed numerically how the dynamics of networks of theta neurons is time asymptotically attracted by the OA manifold 4-6 . Recently, Montbrió and co-workers studied how the macroscopic dynamics of a network of quadratic integrate-and-fire (QIF) neurons is described by a low-dimensional system by using a so-called Lorentzian ansatz 7 . By transforming the QIF neurons into a network of theta neurons, their Lorentzian ansatz does resemble the OA ansatz with parameter-dependent frequency and driving field, as considered in Section II.
To be more precise, the dynamics of the membrane potential V j of a QIF neuron may be described byV for j = 1, . . . , N. Here, I j denotes an input current, V p a peak value, and V r a reset value.
Once the membrane potential V j reaches V p , the neuron emits a spike, and V j will be reset to V r . Commonly, the limit V p = −V r → ∞ is considered. The input current I j consists of a neuron-specific quenched component η j , a common time-dependent input I(t) and a coupling term Js(t), combining the synaptic weight J and a smooth mean synaptic activation s(t), resulting in The latter two time-dependent components are identical for all neurons in the network. In order to describe the macroscopic behavior of the network, Montbrió and co-workers used the Lorentzian ansatz with center y(η, t) and time-dependent half-width x(η, t), which turns out to exhibit the long-term solution for the distribution of the membrane potentials. The properties x(η, t) and y(η, t) that define the distribution function (21) are also closely linked to the firing rate of the neuronal population and to the mean membrane potential, respectively. While the Lorentzian ansatz applies to the (membrane voltage) dynamics of QIF neurons, we are here primarily interested in the phase dynamics. Using V j = tan(θ j /2) one can transform (19 & 20) into theta neurons, 26 In (22) the time-independent injected current η j is drawn from a distribution function g(η).
For the sake of legibility we abbreviate the non-autonomous part of (22) as Rearranging terms and considering the thermodynamic limit, one can rewrite (22) as with To apply our result from above, one has to show that H does not diverge exponentially when Im(η) → −∞, and that c(t) possesses an antiderivative. On the one hand, for the components of c(t) with s(t) being smooth and I(t) piecewise smooth and integrable, there will always exist an antiderivative of c(t). On the other hand, we have H(η) = iη + const, such that H grows only linearly for Im(η) → −∞. Hence, we find that the OA manifold does asymptotically attract the macroscopic behavior of a network of coupled theta neurons.
Due to the existence of a conformal mapping between the quantity w(η, t) = x(η, t) + iy(η, t) and the function α(η, t) defining the OA manifold (4)

IV. GENERAL PARAMETER DISTRIBUTIONS
As already mentioned in Section II, the assumptions of a linear relation between ω and η and of η being drawn from a Lorentzian can be loosened in many respects. We first consider g(η) to still be a Lorentzian centered around η = η 0 with width ∆, i.e. g(η) ∼ L(η 0 , ∆). The linear dependency ω(η, t) = a · η + c may be generalized by considering both a = a(t) and c = c(t) time-dependent. Then, by the common transformation properties for Lorentzian (Cauchy) distributions, ω follows a Lorentzian of the form g(ω) ∼ L(aη 0 + c, ∆|a|). Let a = 0 be constant. Then a similar change of variables,θ = θ − (aη 0 t + C(t)), with C(t) being the antiderivative of c(t), keeps the distribution function centered around 0. Without loss of generality we set a = 1; even if a = a(t) and a(t) > 0 or a(t) < 0 for all t > 0, the rescaling of θ retrieves that we can stick to our assumption a = 1. If, however, a changes sign at, e.g., t = t 0 , then the scale parameter ∆|a| tends to zero for t → t 0 . Due to (4) also ρ(θ, ω, t) will exhibit a δ-peak at t = t 0 . In this case our results are not readily applicable 3 . However, iḟ a(t 0 ) = 0, then we can shift the initial time to zero, t 0 → 0. Whenever ρ + (θ, ω, t 0 ) satisfies the necessary initial conditions, the OA manifold will remain attracting for all t > t 0 , given We proceed with more general cases of frequency and parameter distributions. In Ref. 3, the authors elegantly extend the original proof, which considers only Lorentzian frequency distributions: Instead of demanding analytic continuity of both the frequency distribution g(ω) and the initial condition into the whole lower ω-half plane, it suffices that g and the initial condition have analytic continuations into a strip S defined by 0 ≥ Im(ω) > −σ and −∞ ≤ Im(ω) ≤ +∞ with σ > 0, where neither of them has singularities and both approach zero as |ω| → ∞. Thereby the class of applicable distribution functions includes Gaussians, sech-distributions, and many more, and even multimodal distributions can be incorporated as long as these functions have finite non-zero widths; see references in Ref. 3.
This approach can be adopted and used in our η-parameter-dependent case. For this let us assume again individual oscillators given by (6). As mentioned in Section II, we might be confronted with a nesting of the distributionsg(ω) and g(η) for ω and η. In particular, the latter may determine the first in an oscillator-specific way. That is the reason why the resulting distribution functionĝ(η) can become arbitrarily complicated. However, as long as the analytic continuations ofg and g into the strip S (for some σ > 0 as defined above) do not have singularities, and neitherg nor g features a δ-peak in their time-evolutions, alsoĝ will behave as required. An additional requirement is that the product H(η, t)ĝ(η) satisfies these conditions, too. This means that we have to find a strip S ′ ⊂ S, defined by 0 < σ ′ ≤ σ, in which Hĝ has an analytic continuation, does not have singularities, its time evolution does not feature δ-peaks (if necessary we have to reset the initial time point after such a peak), and that we require |H(η r + iη i , t)ĝ(η r + iη i )| → 0 for |η r | → ∞ and 0 > η i > −σ ′ . In particular, H must not grow faster thanĝ decays, such that the OA manifold continues to capture the long-term dynamics of the system.
Revisiting the example from Section III, where H(η) = i(η + c − 2) andĝ(η) ∼ L(η 0 + c, ∆), we find thatĝ decays exponentially for |η r | → ∞ such that H must not increase at an exponential rate. In fact, H dot not have any singularities in the whole complex Consequently, for large |η r |, the product Hĝ will be dominated byĝ such that all assumptions are fulfilled. Hence, we can confirm again the attractiveness of the OA manifold. 28 We would like to remark that initial conditions on the oscillator distribution function, ρ(θ, η, 0), play an important role. If they fail to be satisfied, this may hinder the OA manifold to attract the dynamics. For an example we would like to refer to Appendix C of Ref. 29, in which the specific time point has to be determined appropriately in order to set up promising initial conditions.

V. APPLICATIONS -REALISTIC SETTINGS
So far, we only considered non-independent frequency and parameter distributions,g(ω) and g(η), respectively. In general, however, one cannot take this "simple" dependence for granted. The additional parameter might be multi-dimensional, i.e. η ∈ R n with n > 1. When considering the thermodynamic limit of infinitely many coupled oscillators, the dynamics (6) may obey Employing the OA ansatz for this system one has to encounter distribution functions given the joint distribution g(ω, η) is a major modification to the setting considered before. Does the OA manifold remain attracting? (24) suggests the phase θ = θ(η, t) to depend on the parameter η in line with our notion of parameter-dependent systems. But it is unclear whether the OA manifold is attracting even without this particular correlation between phase, natural frequency, and driving field. If, however, the OA attractiveness can be proven for systems with generalized natural frequency Ω and driving field H as in (24), this will allow for a further and even broader extension of the existing theory. In the following we first list a few examples for which numerical simulations have been reported and that give strong incentive that the OA ansatz may indeed be valid. We will show how our proof can be adopted, thereby confirm the OA attractiveness, and set the numerical results on solid ground. Last, we provide some general properties of Ω and H for which the OA ansatz holds.
We start with the Winfree model 15 which is an early mathematical description of syn- where h(t) is a smooth function depending only on the mean field z(t) but not on the phase itself 30 . In particular, this model contains time-dependent parameters, see also Ref. 31.
Next, we consider reaction-diffusion systems with heterogeneous, self-oscillating elements.
In particular, we study the mean-field version of the complex Ginzburg-Landau equation, whose equation describes a population of globally coupled limit-cycle oscillators. Hence, we can rewrite the dynamics in form of (24). By introducing a shear (or nonisochronicity) parameter η as an additional random variable and transforming the system through a phase reduction, the governing equations in the continuum limit read [19][20][21] : where K denotes the coupling strength and z = z(t) is the order parameter. The frequency ω and the shear η are drawn from a joint distribution g(ω, η). In contrast to Section II, we explicitly allow the additional parameter η to be drawn from another frequency distribution.
For the joint distribution one has to address two scenarios. Either, the random variables are independent, such that the joint distribution can be split into g(ω, η) = g 1 (ω)g 2 (η), Last but not least, the upcoming branch of heterogeneous mean fields 35 falls in a category whose mean field dynamics can be described along the OA ansatz. The heterogeneous mean field approach deals with networks that are not all-to-all coupled but they exhibit some particular (and sparse) network topology, and therefore can barely be studied analytically.
Given a network with a particular degree distribution, however, it is possible to introduce so-called degree-block variables, whose dynamics govern the evolution of all nodes which have the same degree k. This approach reveals the same equations as the annealed networks approximation 17,36 , which can hence be considered equivalent. Recent studies considered the heterogeneous mean fields of the Kuramoto model, e.g., on scale-free 37-39 and random Erdös-Rényi networks 37 . The starting point is a specifically coupled Kuramoto network with coupling strength K and adjacency matrix A = (a ij ) with i, j = 1, . . . , N, We can cluster various node dynamics by replacing the adjacency term with an expectation value for their node degree η j . Ideally, the underlying topology exhibits some well-defined degree distribution P (η). In the continuum limit N → ∞, these node degrees are substituted in the phase dynamics as weighted, distributed coupling strengths, so that the governing dynamics read where ω and η are drawn from a joint distribution g(ω, η) = P (η)g 1 (ω). This setup is amenable to, e.g., random fields, as has been presented in Ref. 39 where oscillators are enforced through local fields, which find their way into the specific forms for Ω and H.
In all these different classes of parameter-dependent networks, we will show how the OA attractiveness can be regained.

A. Winfree model
As said, the Winfree model describes macroscopic synchronization phenomena of large oscillator systems whose individual nodes are naturally pulse-coupled with one another.
The introduction of phase response curves (PRC) allows for quantifying how the phase of an oscillator responds to the pulse-like perturbations from the other oscillators. The general form of model reads at the single node level where ε denotes the coupling strength, Q is the PRC and P is a pulse-like signal. Following the notation of Pazó and Montrbrió in Ref. 30, we consider PRCs with sinusoidal shape, with an offset parameter σ, and a phase-lag β. Moreover, we assume the pulse-like signal to be smooth, P (θ) = P n (θ) = a n (1 + cos θ) n , with n ∈ N ≥1 controlling the width of the pulses, and a n is a normalizing constant. In the thermodynamic limit, we regain (26) as where the coupling function incorporates the smooth mean field with z the common (Kuramoto) order parameter (10). The frequency Ω(ω, t) = ω + c(t) with c(t) = εσh(t) has a form identical to Section II, where ω follows a Lorentzian frequency distribution g(ω). Since the order parameter z(t) is bounded with |Re(z)| ≤ 1, we have h(t) ≥ 0 for all t ≥ 0. Furthermore, the driving field does not depend on additional parameters, so that our proof can be directly applied, confirming that the OA ansatz holds and the OA manifold indeed captures the long-term dynamics of the Winfree model.
An alternative proof for the case of time-dependent frequency and driving field can be found in Ref. 31. However, as we have depicted in Section IV, our proof generalizes their findings and extends them to a broader class of frequency distribution functions g(ω). Of particular interest in the non-autonomous extension is also the matter of discontinuities.
Recall that in Section III we introduced a time-dependent input current I(t), see (20), which can, e.g., take the form of a square function with jump-discontinuities. Our proof applies to this specific feature and confirms existing numerical results 7 .

B. Limit-cycle oscillations with shear
Investigating collective synchronization usually addresses networks of coupled elementary oscillatory units. The dynamics of these units may be described as normal forṁ where ̺ denotes the radius and ω determines the frequency of rotation on the stable limit cycle with ̺(t) ≡ 1. The parameter η quantifies the shear, or non-isochronicity, of the flow, i.e. how strongly perturbations away from the limit cycle modify the phase dynamics. When we consider an all-to-all coupled population of N ≫ 1 of these oscillatory units, we arrive at the mean-field version of the complex Ginzburg-Landau equation with dissipative couplinġ z j = ̺ j e iφ j . Heterogeneity among the population is promoted by having the frequency ω j and shear parameters η j drawn from a distribution function g(ω, η). In the weakly coupled case, i.e. the coupling strength |K| is small, a phase reduction allows us to describe the dynamics of the system by their phases only. In the continuum limit N → ∞, we can introduce the phase distribution function ρ(θ, ω, η, t). Note that ω and η are independent, so that neither of them is redundant. Accordingly, the order parameter z takes now the form Thus, the phase dynamics reads and the phase distribution function satisfies the continuity equation with v the right-hand side of (38), see also Refs. 19 and 21. Using the notion of (24), the frequency and the driving field are both time-varying and depend on the additional shear parameter η: To assure that the OA manifold indeed exhibits the mean field dynamics of this system with shear, we have to adapt our proof from Section IV for the joint distribution g(ω, η).
The general idea is again to decompose the distribution function ρ in Fourier space into and use symmetry assumptions to focus on ρ + , which again will be decomposed into ρ + = ρ + +ρ ′ + . Whileρ ′ + lies on the OA manifold and has Fourier coefficientsρ ′ +,n = [α(ω, η, t)] n , ρ + solves The assumptions on the analytic continuation properties of Section IV hold -in particular we need analytic continuations with respect to both ω and η into strips S ω and S η . Hence we have to show that Discussing general solutions of (43) given an arbitrary joint distribution function are beyond the scope of this paper. However, for particular g(ω, η) we can affirm the attractiveness of the OA manifold for these parameter-dependent systems. To begin with, we use the assumption of Montbrió and Pazó that the joint distribution can be written as the product of two Multiplying (42) with g(ω, η) and integrating over (ω, η), we can use Fubini's theorem (on the assumption of integrability of Ωgρ + and Hgρ + ) and compute the double integral by changing the order of integration. First, we can evaluate the integral over ω by applying the residue theorem as in Section II and then move on to the second integral, which reads While the term Ωg 2ρ+ can be evaluated at the pole η = η 0 ±iγ (± depending on the contour of integration, which again depends on the coupling K, see also Ref. 19) , we have to assure that the product H(η, t)g 2 (η) vanishes for Im(η) → ±∞. Indeed, the linear growth of H in η, see (40), will be dominated by the exponential decay of g 2 , such that the residue theorem can be applied here, too, which results finally in (12)&(13), from which the claim follows as presented in Section II. As has been shown in Section IV, the restrictions to unimodal Lorentzians can be dropped and the OA attractiveness is sustained. Here we can even handle δ-functions as long as one of the partial distribution functions has finite width: due to the special form of Ω(ω, η, t), the OA ansatz holds for homogeneous frequencies ω j = ω while the shear is heterogeneous and the coupling K > 0 does not vanish.
The case in which the joint distribution g(ω, η) is no longer uncorrelated, i.e. if the first equality in (44) with (ω, K) following a joint distribution g(ω, K). Given the strong resemblance between their numerical simulations and the predictions via the OA ansatz, the authors realized that the latter "formulas were derived on the assumption of at least asymptotic validity of the OA ansatz." 22 They also investigated necessary initial conditions with respect to their analytic continuation and applicability to the OA ansatz. Unfortunately, they did not prove this their system dynamics (45) does not belong the classes of systems considered in the proofs by Ott and Antonsen 1-3 . Recall, a general characterization of correlated joint distribution g(ω, K) = g 1 (ω)g 2 (K) that are applicable for the extended OA ansatz is hardly feasible.
However, for three examples used in literature we can prove that the OA manifold defines the asymptotic evolution of the whole system.
Here, L(ω; ω n , γ n ) denotes a Lorentzian of width γ n > 0 and centered around ω = ω n , and Γ(K) is a multimodal-δ-function. For properly chosen q 1,2 (47) can be regarded a generalization of (46) so that it is sufficient to deal with the former. For simplicity, let us consider N q = 2, i.e. g(ω, K) to be a bimodal joint distribution. Employing g(ω, K) in the definition of the order parameter (37), we see that we can decompose it into z(t) = q 1 z 1 (t) + q 2 z 2 (t) with q 1 + q 2 = 1. That is, we can view our system as two all-to-all coupled populations with population-specific coupling strengths K 1,2 . Given that the frequency distributions are Lorentzians of finite width γ 1,2 we can apply the results for two-population/bimodal Kuramoto models as in Refs. 40 and 41, which confirms the attractiveness of the OA manifold for this kind of joint distributions. The case of multiple Kuramoto populations with specific coupling strengths can be approached by transforming the system into one global system whose oscillators' frequencies follow a multimodal distribution consisting of weighted inhomogeneous unimodal distributions, which can mirror the underlying coupling topology across populations 42 .
It is true that the examples mentioned above are not exhaustive but rather represent a small set of a broad variety of joint distribution functions. Nevertheless, we believe that our results may be a major breakthrough for the applicability of the OA ansatz for systems with more intricate distribution functions.

C. Heterogeneous mean field models
While the general case of uncorrelated joint distributions has already been covered in the preceding Section V B, we would like to concentrate on the specific derivation of the heterogeneous mean field model. Recall the standard Kuramoto model on a given network, where K is the coupling strength and the adjacency matrix is given by A = (a ij ) i,j=1,...,N .
We substitute the adjacency values a jk ∈ {0, 1} by their expectation values a jk ∈ [0, 1], which are given by Introducing the complex order parameter as the dynamics for all nodes with the same degree η k reaḋ In this special form, in which the single nodes are replaced by block-degree variables, we returned to the all-to-all coupling. For a given degree distribution P (η) property (49) also holds in the continuum limit N → ∞ where the governing dynamics read with ω and η being drawn from a joint distribution g(ω, η) = P (η)g 1 (ω). As before we can introduce a phase distribution function ρ(θ, ω, η, t), which fulfills the continuity equation ∂ t ρ + ∂ θ (vρ) = 0 with v the right-hand side of (50). Note, however, that depending on the underlying network topology and its degree distribution P (η), one has to choose the domain of η properly. In the case of a scale-free network, the degree distribution follows P (η) ∝ η −γ with γ > 1. Hence the normalization conditions for the distribution function ρ obey ∞ 1 2π 0 ρ(θ, ω, η, t) dθdη = g 1 (ω) and ∞ −∞ 2π 0 ρ(θ, ω, η, t) dθdω = P (η) .
We can apply the OA ansatz as before. By the same reasoning as in Section V B, we can so prove the OA attractiveness for heterogeneous mean field models, rendering also nonglobally coupled oscillator networks applicable to have their mean field dynamics evolved on a low-dimensional manifold.

D. Non-local coupling
Two months before Ott and Antonsen published their ansatz, Ko and Ermentrout investigated the creation of partially locked states in a network of identical all-to-all coupled oscillators due to inhomogeneous coupling 43 . Instead of heterogeneity of the oscillators' frequencies, it was the coupling heterogeneity that led to partial synchronization. Carlo Laing analytically investigated this network of globally coupled oscillators with coupling strengths drawn from a power-law distribution 44 along the line of the OA ansatz -recall the resemblance to the heterogeneous mean field approach for scale-free networks. Assuming "nearly" identical oscillators, i.e. the frequencies ω were drawn from a Lorentzian with width 0 < ∆ ≪ 1, he could verify the earlier results that were derived via a self-consistency argument 43 , and extend them by including a thorough bifurcation analysis. Our findings in Section V B put these results on a solid mathematical ground.
Of particular interest is Laing's work on a ring of oscillators 44,45 . For a given ring topology, the typical coupling scheme is neither local neighbor-to-neighbor, nor global coupling.
Instead, the oscillators are non-locally coupled via a coupling kernel G. We assume that each oscillator k = 1, . . . , N has some fixed spatial position x k ∈ [−π, π], a natural frequency ω k drawn from a continuous distribution function g(ω) with non-zero width, and interacts with the others depending on the distance between their sites modulo periodic boundary conditions. The governing dynamics reaḋ where α is a phase-lag parameter and G : R → R a continuous even and 2π-periodic coupling function 46 . We retrieve global coupling, if G = 0 is constant. Commonly used coupling functions G are of exponential form G(x) ∼ e −κ|x| with κ > 0, or of trigonometric form The reflection symmetry of G is lost for B = 0. In the continuum limit, the velocity field (7) becomes While the inner two integrals have the form of a local complex order parameter z(y, t), measuring the synchronization degree of oscillators around y, we can interpret the last integral as a convolution of the local order parameter with the (spatial coupling) kernel G. In particular, we can regard the dynamics ∂ t θ(x, t) of an oscillator at position x as being controlled by the local mean field H(x, t). Unlike the case of global coupling, the order parameter has become space-dependent and thus the driving field. However, a similar "physical picture" as for global coupling is valid: practically we deal with an assembly of independent oscillators under the control of a common forcing field 46,47 . We now go a step further and interpret the space variable x as a subpopulation index 48 . Equivalent to the block-degree variables in the heterogeneous mean field approach, we consider the subpopulation index as a parameter that follows a particular, in this case a uniform, distribution function. Hence, (52) represents the governing dynamics of a parameter-dependent system, for which we proved the OA attractiveness in the preceding sections.

VI. RELAXATION DYNAMICS
As discussed, we allow time-varying parameters to affect the oscillator dynamics. The change of parameters comes with its time scale(s). The change can be periodic. This periodicity may also influence the evolution of the mean field and thereby the OA manifold.
Therefore, the relation between this periodicity and the characteristic time of the system to approach the manifold needs to be investigated. If the relaxation dynamics onto the manifold is way slower than the characteristic time scale of the time-varying manifold itself, then our findings will remain true for the limit t → ∞. They are, however, of minor interest for describing the transient behavior of the mean field. Several numerical results 7,22,23,31 suggest that the relaxation to the OA manifold is reasonably fast, in some cases even instantaneous.
To address this analytically, we briefly recall the proof for the attractiveness from Section II. After having Fourier expanded the phase distribution function ρ(θ, η, t), and then decomposed the positive Fourier modes into a part that already lies on the manifold,ρ ′ + , and a residual partρ + , we showed how the latter converged to zero in a weak sense, cf. (11). We can extract the relaxation time to the OA manifold from out of the proof: From (12)& (13) we obtain a solution f + (θ, t) =ρ ′ + (θ, −iσ, t), with σ ′ > σ > 0 whereρ ′ + (θ, η, t) admits an analytic continuation into the strip the solution (15) obeysf hence the relaxation time τ is by definition . showed that for frequency distributions g(ω) supported on the whole real axis, the decay towards the incoherent state is exponentially fast for coupling strengths below the critical threshold, K < K c . If g(ω) has compact support, i.e. g is non-zero only on a compact interval [−γ, γ] ⊂ R, 0 < γ < ∞, the rate may be considerably slower, even polynomial. In the example they used to illustrate their result, the authors assumed the frequencies ω to be distributed uniformly on I = [−γ, γ], i.e. g(ω) = 1/2γ if ω ∈ I, and 0 otherwise. The jump discontinuities of g on ∂I, however, prohibited an analytic continuation of g into a strip S in the lower complex ω-plane, contradicting the required conditions for applying the OA ansatz 3 . That is why the proofs above cannot be applied here, and our argumentation about the relaxation times remains unaffected.
Last but not least, we would like to add that decay times typically depend on initial conditions. Pikovsky and Rosenblum pointed out that for identical macroscopic, i.e. mean field, initial conditions the microscopic initial states can lead to very different transient dynamics towards the OA manifold, see Section 3.2 in Ref. 18. A more thorough investigation about this specific topic has not been undergone, yet, but might shed light on the underlying dynamics of the microscopic variables of large oscillatory systems in contrast to its mean field behavior.

VII. DISCUSSION AND CONCLUSION
The OA ansatz has proven rather fruitful for investigating the macroscopic behavior of systems of coupled phase oscillators in terms of a low-dimensional system. Although parameter dependence has already been mentioned in Ott and Antonsen's original work, parameters were merely considered auxiliary variables and the velocity field was required to incorporate the phase only through a sinusoidal coupling term.
Our main result was to prove that the η-dependence sustains the time-asymptotic attractiveness of the OA manifold for systems of coupled oscillators. For this we required that the driving field H does not have singularities in the complex η-plane and that it diverges at most sub-exponentially for Im(η) → −∞, next to the conditions in the original Ott & Antonsen formulation 1,2 . Furthermore, we assumed the frequency ω(η, t) to be linear in η.
We were able to depict the proof step by step. Subsequently we loosened the restrictive assumptions and showed that our results remain valid for a much broader class of distribution functions g(η) as well as more complex dependencies of the driving field H(η) and the natural frequencies ω(η) on the parameter η.
Although the main idea of introducing a common parameter η was to correlate the driving field and the natural frequency with their specific oscillator, our proof is identical for the case when η does only influence the mean field dynamics. By this, we have proved the claim in Ref. 3 that the OA manifold remains attractive in the "weak" parameter-dependent case when H depends on "other non-phase-oscillator variables obeying auxiliary dynamical systems." Common choices of H and ω usually fulfill the aforementioned assumptions as stated in Section II. That is, our result can be immediately applied in a variety of circumstances. Here, we highlighted an application in mathematical neuroscience. By this, our findings strengthen the theory of coupled theta neurons: The many recent numerical findings of Ref. 7 and the references therein are finally set in a solid mathematical framework. Moreover, the link between QIF neurons and theta neurons has been underscored by proving the attractiveness of the Lorentzian ansatz.
We generalized and extended existing proofs for non-autonomous systems. In particular we addressed the Winfree model, which is biologically more realistic than the Kuramoto model and therefore closer to applications. We also addressed coupled oscillatory systems with an additional shear parameter, another important tool to render the Kuramoto model more realistic. The major novelty was our rigorous proof of the OA attractiveness for systems with uncorrelated joint distribution functions when more parameters than only the natural frequencies are treated as a random variable. This finding opened the way for networks with specific underlying coupling topologies other than the restrictive global coupling. Using the heterogeneous mean field approach, we showed how these networks can be treated along the OA ansatz. First steps were also taken in the direction of correlated joint distributions.
All in all, we consider the explicit dependence on an additional parameter η of both the oscillator's phase and the (non-sinusoidal) components an important extension introducing an intrinsic relation between phase, frequency, and driving field of an oscillator. The latter two are correlated with the phase so that the η-dependence does not allow for applying the original theory.
Still, there are several open problems concerning the mean field dynamics of an oscillatory system and its description by a low-dimensional system. A first urgent one is the case of δ-peaked frequency distributions. Numerical simulations 54 and heuristic arguments hint at convergence of the OA manifold, where a proper mathematical derivation is omitted under the pretence of "nearly identical oscillators" 45,55,56 . A thorough proof would render the OA ansatz rigorously applicable to "chimera states", a topic that is particularly en vogue; see, e.g., the recent review paper by Panaggio and Abramscite 57 . Importantly, such a proof has to circumvent the main argument of Ott and Antonsen's original proof, where the width ∆ > 0 of the distribution g(ω) allowed for a consequent evaluation of the mean field dynamics. On the other hand, Pikovsky and Rosenblum 58 already showed that more complicated dynamics can emerge from the OA manifold when describing the system along Another intriguing open problem is whether the mean field dynamics is attracted by a low-dimensional manifold when the parameter-dependence of the frequency and driving field is extended by an explicit dependence on the phase. A recent example is given by Laing 60 , who considered the driving field H to follow a dynamics that explicitly depends on the phase θ. This system exhibits partial synchronization patterns, which are also covered by the OA ansatz, but any attempt to apply the OA ansatz has been avoided "due to the dynamics of the extra variables." 60 When the coupling term incoporates higher harmonics, see, e.g., Refs. 61 and 62, no lowdimensional analytic solution for the mean field evolution has been found. This is another open question whether further generalizations of the work of Ott and Antonsen [1] can be rigorously manifested. We believe that our current proof for parameter-dependent networks is a good starting point for tackling these important issues.

ACKNOWLEDGMENTS
This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sk lodowska-Curie grant agreement #642563 (COS-MOS). We would like to thank Ernest Montbrió, Oleh Omel'chenko, Spase Petkoski and