Multi-Population Phase Oscillator Networks with Higher-Order Interactions

The classical Kuramoto model consists of finitely many pairwise coupled oscillators on the circle. In many applications a simple pairwise coupling is not sufficient to describe real-world phenomena as higher-order (or group) interactions take place. Hence, we replace the classical coupling law with a very general coupling function involving higher-order terms. Furthermore, we allow for multiple populations of oscillators interacting with each other through a very general law. In our analysis, we focus on the characteristic system and the mean-field limit of this generalized class of Kuramoto models. While there are several works studying particular aspects of our program, we propose a general framework to work with all three aspects (higher-order, multi-population, and mean-field) simultaneously. Assuming identical oscillators in each population, we derive equations for the evolution of oscillator populations in the mean-field limit. First, we clarify existence and uniqueness of our set of characteristic equations, which are formulated in the space of probability measures together with the bounded-Lipschitz metric. Then, we investigate dynamical properties within the framework of the characteristic system. We identify invariant subspaces and stability of the state, in which all oscillators are synchronized within each population. Even though it turns out that this so called all-synchronized state is never asymptotically stable, under some conditions and with a suitable definition of stability, the all-synchronized state can be proven to be at least locally stable. In summary, our work provides a rigorous mathematical framework upon which the further study of multi-population higher-order coupled particle systems can be based.

patterns. Synchrony in its many forms often is also relevant for the function of a network dynamical system: In neural networks, synchrony is linked to cognitive function as well as disease [51].
A classical question in network dynamical systems is how the network structure and functional interactions shape the collective network dynamics. The Kuramoto model has been instrumental in understanding synchronization in coupled oscillators [32,45,1]. In the case of N identical oscillators, the state θ k ∈ S := R/2πZ of oscillator k ∈ {1, . . . , N } evolves according tȯ where ω ∈ R is the intrinsic frequency of the oscillators. To understand (global) synchronization, the Kuramoto model makes several simplifying assumptions: The network is homogeneous and interactions are pairwise and all-to-all, mediated by the sine of the phase differences. The dynamics of (1.1) are particularly simple: The dynamics are effectively two-dimensional and for generic initial conditions all oscillators synchronize [54]. However, many real-world network dynamical systems are more complex and cannot be captured using the Kuramoto model. First, one needs to consider interaction functions that contain more than a single harmonic [17], for example if the interactions are state-dependent [2]. Considering a general coupling function g : S → R breaks the degeneracy and chaotic dynamics are possible even for small oscillators [13]. Second, many networks that arise in real-world systems are not homogeneous or allto-all coupled but have some form of modularity or community structure [22,40]: There are different communities or populations that are characterized by the property that the coupling within a population is different from the coupling between populations. Even if the populations are identical, multi-population oscillator networks give rise to a variety of synchrony patterns; see [10] for a recent review. Third, classical network dynamical systems, such as the Kuramoto model (1.1), assume that interactions take place in terms of pairs: The influence of two nodes on a third is simply the sum of two individual contributions. For example, Skardal and Arenas [43] consider the phase oscillator dynamics given bẏ with nonadditive interactions that involve triplets and quadruplets of oscillators. Indeed, nonpairwise higher-order interactions arise naturally when considering phase reductions (cf. [3,34] and the appendix) and recent work has highlighted the dynamical importance of higher-order interactions [6,11]. In all these cases, going beyond the Kuramoto model leads to new dynamics. Further examples of new dynamics are possible if these features are combined, e.g., modular networks of phase oscillator networks with higher-order interactions allow for heteroclinic structures between different synchrony patterns [7,8,12]. So far, rigorous insights that relate network structure and dynamics have relied on specific assumptions on the size of the network and the type of interactions. Classical dynamical systems techniques can be applied if the networks have relatively few nodes: The heteroclinic structures between synchrony patterns have been shown to exist in networks of up to nine phase oscillators [8]. Since many real-world networks have many nodes, it is often instructive to consider these networks in the limit of infinitely many nodes. In this limit, the Ott-Antonsen reduction [41,10] has been instrumental to elucidate the dynamics of coupled oscillator populations. The main drawback is that this reduction requires network interactions that are mediated by a single harmonic, whether it is pairwise or contains higher-order terms. Moreover, the limit itself has only been rigorously justified in a restricted setup.
Here, we give general results on the mean-field limit of multi-population phase oscillator networks with higher-order interactions. We rigorously describe the evolution of the mean-field limit where each population is represented by a probability measure that describes the distribution of oscillators on the circle. Moreover, we identify how phase space is organized for a class of coupled oscillator networks with higher-order interactions and calculate the (local) stability properties of synchrony patterns. Our results contribute to the understanding of coupled oscillator networks in several ways. First, we generalize the results of [33] to transport equations involving finitely many oscillator populations with higher-order interactions. Second, our stability analysis complements the work in [15] and [19] by considering a general setting with multiple oscillator populations and general coupling with nonsinusoidal pairwise and nonpairwise higher-order interactions. Third, our explicit stability results on a measure-valued evolution shows in this setting one cannot expect asymptotic stability of full phase synchrony but just a weaker form of stability. Finally, our results provide a first step towards understanding global dynamical phenomena in the general mean-field limit: Our stability results outline necessary conditions to prove that the heteroclinic structures for small networks exist in the mean-field limit of large networks.
This work is organized as follows. In the remainder of this section we fix some notation that will be used throughout and introduce a general system of equations that describe the network dynamics. In Section 2 we establish the existence and uniqueness of the equations in the space of probability measures. This section also covers the relation of the mean-field limit to the continuity equation and networks of finitely many oscillators. The main results regarding synchronization are given in Section 3. We apply our results to three explicit examples of coupled phase oscillator networks in Section 4. Finally, in Section 5 we give some concluding remarks and an outlook on future research.

Notation
We first fix some notation that will be used throughout this paper. Let P(X) denote the set of all Borel probability measures on the set X. If S = R/2πZ is the unit circle then P(S) represents the set of all Borel probability measures on the circle. The symbol P ac (S) is used to denote the set of absolutely continuous probability measures, i.e., those which have a density, on the unit circle. Whenever we write α 1 − α 2 for two points α 1 , α 2 ∈ S on the circle, we refer to the value of α 1 − α 2 ∈ [0, 2π). Further, let us define open intervals on the circle as (α 1 , α 2 ) := {α 1 + t(α 2 − α 1 ) : t ∈ (0, 1)}.
To compare two measures µ, ν ∈ P(S), we use the Wasserstein-1 distance [52], which is also referred to as the bounded-Lipschitz distance

Solution Theory
In this section we first introduce the system of characteristic equations stated below as (2.1)-(2.2). Then, we establish its well-posedness, its continuous dependence on initial conditions, and finally relate the solution of the characteristic system to mean-field limit equations of Vlasov-Fokker-Planck type and finite-dimensional generalized Kuramoto models with multiple populations and higher-order coupling.

The System of Characteristic Equations
In this paper, we consider the dynamics of M ∈ N coupled phase oscillator populations. We now introduce a general set of equations that describes the network evolution, where the state of population σ ∈ [M ] is given by a probability measure µ σ . The network interactions are determined by a multi-index s σ ∈ [M ] Rσ for each population together with Lipschitz continuous coupling functions G σ : S |s σ | × S → R. Specifically, these coupling functions are supposed to be L-Lipschitz when S |s σ | × S is considered with the metric d(α, β) = . . , µ in M ) ∈ P(S) M denotes the initial state of the network, # denotes the push-forward operator and µ = (µ 1 , . . . , µ M ), then the evolution of µ(t) = (µ 1 (t), . . . , µ M (t)) is determined by the characteristic equations for σ ∈ [M ] and the evolution operator where ω σ ∈ R is the instantaneous frequency of all oscillators in population σ. We remark that the general idea of using a mean-field formulation involving probability measures instead of densities is quite classical [24]. Our equations (2.1)-(2.2) provide a very general variant of this principle allowing for multipopulation higher-order coupled systems. Yet, the formulation also naturally reduces to classical cases, as we illustrate in Section 4.
1a) together with (2.1b),(2.1c) and (2.2), they are referred to as the mean-field characteristic flow. In this case, µ ∈ C M P(S) , given by (2.1b), is a solution of the system (2.1)-(2.2). Remark 2.2. For a given mean-field characteristic flow Φ σ , the solution µ of the characteristic system is uniquely given by (2.1b). Conversely, for a given solution µ ∈ C M P(S) , the mean-field characteristic flow Φ σ is unique, as one can see by integrating (2.1a) and (2.1c).

Existence, Uniqueness and Continuous Dependence on Initial Conditions for the Characteristic System
We start with existence and uniqueness building upon ideas by Neunzert [39] developed in the context of more classical single-population kinetic models, which lead to similar mean-field limits in comparison to our system (2.1)-(2.2). Since the proof given in [39] makes abstract assumptions about the coupling function and only deals with one population, there are quite some differences to the existence and uniqueness proof of our system, which is why we decided to include the full details here. First, we have to define a suitable space for solutions.
is continuous. Let C P(S) be the set of all weakly continuous functions µ : [0, T ] → P(S).
Remark 2.4. As T > 0 is arbitrary, we do not explicitly include T in the notation C P(S) for functions mapping from [0, T ] to P(S). Since the following existence and uniqueness result as well as results regarding continuous dependence on initial conditions are valid for all T > 0, they can be extended to hold on the half-open interval [0, ∞).
Lemma 2.5. The space C P(S) together with the metric is a complete metric space.
Proof. Let (γ n ) n∈N be a Cauchy sequence in P(S). Then, we obtain from [25, Satz 3] the convergence of (γ n ) n∈N to an arbitrary positive measure. Testing with the constant 1-function, i.e., choosing f ≡ 1 in the representation (1.3b), yields that the Cauchy sequence is even converging to a probability measure. This shows the completeness of (P(S), W 1 ). Now, let (µ n ) n∈N be a Cauchy sequence in (C P(S) , d). The completeness of (P(S), W 1 ) causes the existence of measure-valued functions µ ∞ : [0, T ] → P(S). To show weak continuity of µ ∞ , let f ∈ C(S) first be 1-Lipschitz continuous and calculate is a complete metric space.
For a given µ ∈ C M P(S) , let T σ t,s [µ] denote the flow induced by the velocity field (K σ µ(t))(φ), i.e., We proceed to verify several assumptions and constructions in the arguments in [39] within our new setting.
Lemma 2.7 (cf. [39], Assumption on page 236). For all σ ∈ [M ], φ, ψ ∈ S and µ ∈ P(S) M we have Proof. Using the assumption about Lipschitz continuity of G σ , we can estimate This completes the proof.
In particular, Proof. This is a direct consequence of Lemma 2.7 and Gronwall's Lemma.
Lemma 2.11. For all µ, ν ∈ C M P(S) and Proof. Integrating (2.4) and then using Lemmas 2.7 and 2.10, we can estimate Gronwall's Lemma gives us v σ (t) ≤ Ls σ e Lt t 0 w(τ )e −τ L dτ, which was the claim. . Let h 1 , h 2 : S → S be bijective measurable mappings and let µ ∈ P(S). Then, For given µ in ∈ P(S) M we now consider the mapping Lemma 2.13. The mapping A, defined in (2.7), is indeed a self mapping.
Proof. It is clear, that (Aµ) σ (t) is again a probability measure, so it is left to show that t → (Aµ) σ (t) is weakly continuous. Note that for each φ ∈ S, the map t → T σ t,0 [µ]φ is continuous (even differentiable). Thus, for any f ∈ C(S), the composition t → f (T σ t,0 [µ]φ) is continuous as well and uniformly bounded. Consequently, by the change of variables rule for push forward measures and the dominated convergence theorem, we have which proves weak continuity.
Finally, we can establish existence and uniqueness: Theorem 2.14 (cf. [39], Theorem 2). For given µ in = (µ in 1 , . . . , µ in M ) ∈ P(S) M there exists a unique solution µ ∈ C M P(S) for the system (2.1)-(2.2). Proof. Let us prove, that for α large enough, A is a contraction with respect to the metric d α , defined in (2.3). To estimate d α (Aµ, Aν), we use Lemmas 2.11 and 2.12: This shows that A is a contraction for suitable values of α. Therefore, by the contraction mapping principle, A has a unique fixed point µ. To conclude, if µ is the fixed point of A, µ satisfies (2.1)- 2), µ is also a fixed point of A, which completes the proof.
As a next step, we want to establish continuous dependence on initial conditions for the characteristic system. Since this system links to the finite-dimensional generalized Kuramoto models via the empirical measure as well as to the mean-field limit, we may later conclude from a suitable continuous dependence that the mean-field limit is a good approximation on finite time scales for the generalized Kuramoto system of oscillators. Lemma 3]). Let T : X → X be a surjective and Lipschitz continuous mapping with Lipschitz constant L and let µ, ν ∈ P(X). Then, The next result is sometimes also referred to as a Dobrushin-type estimate: Theorem 2.16. Let µ in , ν in ∈ P(S) M be two different initial measures and let µ, ν : [0, T ] → P(S) M be the solutions to system (2.1)- , π σ (t)) + W 1 (π σ (t), ν σ (t)) .
By using Lemma 2.11, Lemma 2.12 and the notation introduced in (2.6) we see, that By Lemma 2.15, we obtain where L T is the Lipschitz constant of T σ t,0 [ν]. By Corollary 2.8, we have L T = e Lt . Putting everything together, we get Dividing by e Lt and then applying Gronwall's Lemma yields that so the result follows.

Vlasov-Fokker-Planck Mean-Field Equation for Absolutely Continuous Measures
As a first special case, we can state the mean-field limit PDE, which is given by the initial value problem Here, µ ρ(t) ∈ P(S) M is the collection of measures whose densities are given by (ρ σ (t, ·)) σ∈ [M] . The initial conditions are given such that We already remark that we are eventually going to show below the natural interpretation that ρ σ (t, φ) describes the probability of finding an oscillator of population σ at time t at a position φ. • For every f ∈ C(S) and for all σ ∈ [M ], the maps t → S f (α)ρ σ (t, α)dα are continuous.
Proof. We consider the measures µ in σ given by µ in σ (A) = A ρ in σ (φ) dφ. Theorem 2.14 yields the existence of a unique solution µ(t) = (µ 1 (t), . . . , µ M (t)) for the system (2.1)- So µ σ (t) is absolutely continuous with density To show that (ρ σ ) σ∈[M] is a weak solution of the initial value problem, we take By using the chain rule, it turns out, that h σ (t, Putting these two calculations together, we obtain This shows that the collection (ρ σ ) σ∈[M] is a weak solution of the initial value problem (2.8)-(2.9). Remark 2.20. In case the initial conditions ρ in σ (φ) and the coupling functions G σ are sufficiently smooth, it follows from standard ODE theorems, that the flow T σ t,0 [µ](φ) and thus also the densities ρ σ (t, φ) at later times are smooth with respect to φ.

Discrete Initial Measures
As a further special case of the characteristic system (2.1)-(2.2), one is often interested in choosing µ in as a discrete measure. In particular, initial measures of the form describe the discrete states of N oscillators in each of the M populations. It can easily be seen, that if functions φ σ,k : R ≥0 → S solve the generalized Kuramoto ODE systeṁ with initial condition φ σ,k (0) = φ in σ,k , then the measures solve the characteristic system (2.1)-(2.2).

Relation between Absolutely Continuous and Discrete Initial Measures
In order to find a relation between absolutely continuous initial measures and discrete initial measures in the characteristic system (2.1)-(2.2) we apply Theorem 2.16 for the case of µ in ∈ P(S) M consisting of discrete measures and the components of ν in ∈ P(S) M being absolutely continuous measures. In fact we even consider a whole family of discrete initial measures given by (2.13) where N parameterizes the family (µ in,N σ ) N ∈N and Ω N,σ denotes the amount of discrete oscillators in population σ for µ in,N σ . Suppose that this family of measures converges for each σ ∈ [M ] as N → ∞ in the Wasserstein-1 distance to measures ν in σ , which can be represented by densities ρ in σ . The solution of the original system (2.1) for initial conditions given by (2.13) can be obtained by solving the ODE (2.10)-(2.11), whereas the solution of this system for initial condition given by the densities ρ in σ can be obtained by solving the mean-field limit (2.8)-(2.9). The mean-field limit is analytically much easier than systems of the form (2.11) with large N . Therefore, one typically tries to find results for the mean-field limit (2.8) first. Then, one employs Theorem 2.16 with µ in,N σ for µ in and ν in being the limit of µ in,N σ as N → ∞. This theorem gives an upper bound of how well the mean-field limit approximates solutions to the finite N system (2.11). In fact, Theorem 2.16 tells us that whenever a sequence of discrete initial measures (2.13) converges to absolutely continuous initial measures, represented by densities ρ in σ , then the solutions of system (2.11)-(2.12) will also converge at a fixed later time t to densities ρ σ (t, ·) which solve the initial value problem (2.8)-(2.9).

Specific coupling functions
Now that we have set up a general framework for multi-population oscillators with higher-order coupling, it is important to stress that the additional terms beyond the classical Kuramoto model, can also induce new dynamical effects. This is one main motivation to keep the theory as general as possible. Yet, to see that new effects occur, it is interesting to consider specific coupling functions, so that the dynamics of the characteristic equation (2.1) are more restricted.
To illustrate the emergence of new dynamics, we first consider finite networks of N phase oscillators with sinusoidal coupling. As shown in [54], a trajectory of a coupled oscillator system of the forṁ where φ k ∈ S is the phase of each oscillator and f, g, h : S N → R are functions of φ = (φ 1 , . . . , φ N ), does not explore the whole N -dimensional phase space, but is rather constrained to a three-dimensional manifold. In particular, when employing the transformation where γ(t), Θ(t) and Ψ(t) are functions satisfying the 3-dimensional ODE systeṁ the coefficients ψ 1 , . . . , ψ N in (2.15) are actually constant in time [54]. This degeneracy is a consequence of the sinusoidal coupling and can be understood in terms of Möbius group actions [35,44]. Even though (2.14) seems quite restrictive, it is actually contains a large class of systems going beyond the Kuramoto [16]; see also [10]. For example, as one can see by applying trigonometric identities, the system (1.2) from the introduction belongs to this class. When converting an initial condition φ(0) = (φ 1 (0), . . . , φ N (0)) of (2.14) into initial conditions γ(0), Ψ(0), Θ(0) of (2.16) and constants ψ 1 , . . . , ψ N , there are N equations (2.15), but N + 3 variables to choose. Therefore, three degrees of freedom are left, and one typically chooses the constants ψ 1 , . . . , ψ N such that It was shown in [54], that one can find suitable parameters ψ 1 , . . . , ψ N satisfying (2.17) and γ(0), Ψ(0) and Θ(0) for almost all initial conditions φ(0) = (φ 1 , . . . , φ N ). In particular, the only exception is when there are so called "majority clusters", when the phases of at least N/2 oscillators coincide. By imposing the restrictions (2.17), one can think of the parameters ψ 1 , . . . , ψ N as the ones foliating the N -dimensional phase space into a N − 3 dimensional continuum of 3-dimensional invariant manifolds. The dynamics inside these manifolds is then described by γ(t), Ψ(t) and Θ(t).
Let us now apply this theory to the system (1.2), that we considered in the introduction, with K 2 = 0: Here, k = 1, . . . , N , ω ∈ R is the intrinsic frequency of all oscillators and K 1 , K 3 ∈ R are coupling constants. By introducing the order parameter this system can also be written aṡ The equations (2.16), determining the evolution of γ(t), Ψ(t) and Θ(t), then turn intȯ Since the order parameter fulfills as one can see by using trigonometric identities (see [54,Eq. (3.3)]), it is a function of γ and Ψ only. Therefore, Θ does not appear on the right-hand side of (2.20) and consequently, the dynamics of (2.20) is essentially two dimensional. Analogously to [54], we can now define a potential function As calculated in [54], the partial derivatives of H are given by .
This allows us to rewrite the (γ, Ψ)-dynamics of (2.20) in terms of partial derivatives of H: Using this, we can now calculatė However, since as calculated in [54], we obtainḢ = r 2 (K 1 + K 3 r 2 ). At this point it is important to understand the meaning of H. As explained in [54], on the one hand and on the other hand Consequently, whenever K 1 + K 3 r 2 > 0 and r > 0, H increases until either r = 1 or Example 2.21. Let K 1 = 1, K 3 = −4 and the initial condition φ be given such that there are no majority clusters and 0 < r < 1/2. H is stationary only if r = 0 or r = −K 1 /K 3 = 1/2. Otherwise H is increasing. Thus, there are only two possibilities: Firstly, H can diverge to infinity or secondly lim t→∞ H exists in R. In the first case, r → 1, by (2.21) and consequently r has to achieve the value 1/2 at an intermediary time. This, however, cannot happen, since at r = 1/2, all oscillators of (2.19) are only rotating around the circle with a common frequency. So r = 1/2 is invariant and thus r → 1 is impossible. Consequently lim t→∞ H exists in R, which in turn can only happen if r → 1/2. To conclude this example, almost all initial conditions display r = 1/2 in the limit as t → ∞, and the oscillators consequently partially synchronize. This occurs even though the intrinsic frequency ω in (2.18) does not depend on k as it is typically the case for the classical Kuramoto model when one speaks of partial synchronization. Moreover, this observation is in line with the emergence of new stationary phase configurations when considering "nonlinear coupling" through the order parameter [16] beyond the Kuramoto model.
While the above analysis applies only to discrete initial measures, a similar analysis can be conducted when the initial measure has a density ρ in (φ). In this case, the evolution of ρ is governed by the continuity equation. Then, there are no constants ψ 1 , . . . , ψ N , but there is a density χ(ψ) on the circle that does not depend on time. The transformation to get from φ to ψ is exactly the same as described in (2.15) and the requirement for the density χ(ψ) to be independent of time is given by (2.16), as well. Since the initial density ρ cannot have any majority clusters, it is always possible to find initial conditions for γ, Ψ, Θ such that the density χ satisfies In general, sums that depend on (ψ k ) k=1,...,N and appear when dealing with discrete measures have to be replaced by integrals with respect to the density χ(φ) in the absolutely continuous case. Therefore, the potential function H reads as still holds true. Consequently, the explanation in Example 2.21 can be adopted to absolutely continuous measures.
As we have seen in the previous examples, the order parameter can give insights about the degree of synchrony in an oscillator system. The dynamics of this order parameter has interesting bifurcations at r = 0 and r = 1 but there can also be further equilibria of r ∈ (0, 1). For instance, in Example 2.21, r = 0 and r = 1 were repelling equilibria whereas r = 1/2 was attracting. It is therefore interesting under which conditions certain invariant states of an oscillator system are stable or unstable. Since r = 1 in each population is the only invariant state of the general characteristic oscillator system (2.1) and under minor assumptions, r = 0 includes a general invariant state, as well, the next section is devoted to a stability analysis of these two states.

Synchrony and Synchrony Patterns and their Stability
Of course, instead of working with the simpler mean-field limit, or with a particular finite-dimensional Kuramoto model, an alternative route is to directly study the characteristic system. This system links finite-dimensional oscillator systems to the mean-field. In particular, it contains the full information about both systems, so we start our dynamical analysis here using the characteristic system.
From now on we assume the multi-indices s σ to have the form for some L σ ∈ N and a σ j ∈ [M ], i.e., |s σ | = 2L σ + 1. Further, the coupling functions G σ : S |s σ | × S → R are supposed to be of the form G σ (α, φ) = g σ (α 1 − α 2 , . . . , α 2Lσ−1 − α 2Lσ , α 2Lσ+1 − φ), for functions g σ : S Lσ × S → R. These conditions can be summarized by requiring the velocity field K σ µ to be given by for new multi-indices r σ = (a σ 1 , . . . , a σ L σ ) with |r σ | = L σ and coupling functions g σ : S |r σ | × S → R. In the remainder of this paper we study the system (2.1) with (3.2). We adopt the notation in [8] and write D if a population is in splay configuration and S if it is phase synchronized. That is, we write µ 1 · · · µ σ−1 Sµ σ+1 · · · µ M = µ ∈ P(S) N µ σ ∈ S , (3.3a) to indicate that population σ is fully phase synchronized or in splay phase. We extend the notation to intersections of the sets (3.3). Consequently, S · · · S (M times) is the set of cluster states where all populations are fully phase synchronized and D · · · D the set where all populations are in splay phase. Proof. Without loss of generality we can assume that m = 1, since the case for general m then follows by repeatedly applying this proposition. After reindexing populations, we can also assume the M th population to be fixed in synchronized or splay state. For a convenience of notation we suppose the multi-indices r σ to be sorted in ascending order, which can easily be achieved by changing the order of integration. Let us now denote χ σ = |{i : r σ i = M }| and write v σ ∈ [M ] p σ with p σ = |r σ | − χ σ for the multi-index having the same entries as r σ except that the last χ σ entries, i.e., all entries valued M , are missing. In case we fix the M th population to be synchronized, the other M − 1 populations yield a system of the form (2.1),(3.2) with coupling functionŝ

Stability of the All-Synchronized State
As we have seen, the state in which one population is synchronized, is invariant. However, as the synchronized state is not really a single point, but rather a whole set of phase configurations it makes more sense to study the stability of this set. Similarly, for the multi-population model, it makes sense to analyze the stability of the set S M := S · · · S = {µ ∈ P(S) M : µ = (δ ξ1 , . . . , δ ξM ), ξ 1 , . . . , ξ M ∈ S}.
Notation 3.5. For an ease of notation, we will write where B(δ, ǫ) denotes the ball centered at a measure δ of radius ǫ in the Wasserstein-1 metric.
The following definition is the natural variant of (Lyapunov) stability in our setting: Definition 3.6. The set S M is stable if for all σ ∈ [M ] and all neighborhoods U σ ⊂ P(S) of S there exist neighborhoods V σ of S such that for any µ in = (µ in 1 , . . . , µ in M ) ∈ V 1 × · · · × V M , the solution µ(t) of (2.1),(3.2) satisfies µ(t) ∈ U 1 × · · · × U M for all t ≥ 0.
Of course, one often does not only want to show stability of solutions staying near an invariant set but also that the solutions tend towards the invariant set: Definition 3.7. The set S M is asymptotically stable if it is stable and, additionally, there exists a neighborhood V = V 1 × · · · × V M ⊂ P(S) M such that for all µ in ∈ V the solution of (2.1),(3.2) satisfies µ σ (t) → S as t → ∞ for all σ ∈ [M ].
Remark 3.8. The two Definitions 3.6 and 3.7 are formulated in terms of the topology created by the distance max σ∈[M] W 1 (S, µ σ ). However, instead of taking the maximum, one can also sum over W 1 (S, µ σ ), to get a definition that resembles the metric used to prove existence and uniqueness of (2.1) more closely. However, as these topologies are equivalent, we use Definitions 3.6-3.7 in the topology generated by max σ∈[M] W 1 (S, µ σ ), as this setting seems easier to work with. Remark 3.9. The two Definitions 3.6 and 3.7 also make sense when considered with the more general coupling (2.2). However, we only work with them in the context of the velocity fields (3.2).

No Generic Asymptotic Stability
This section aims to illustrate that the all-synchronized state can not be asymptotically stable under a generic condition, which is in this case ∂ ∂γ g σ (0, γ) = 0 for at least one σ ∈ [M ]. Assume that we do not have generic asymptotic stability for M = 1. Generalizing this to M populations can then be easily done by applying the results for one population to an invariant subset of the form S · · · Sµ σ S · · · S, with σ ∈ [M ] chosen such that ∂ ∂γ g σ (0, γ) = 0. No asymptotic stability in this subset with one free population then yields no asymptotic stability in the whole system with M free populations.
Thus, we only consider the case M = 1, so we assume ∂ ∂γ g 1 (0, γ) = 0. The strategy of the proof is to construct a sequence of steady states converging to the synchronized state. Along this family, no asymptotic convergence can take place. To accomplish this construction, consider a perturbation of the synchronized state of the form for φ in 1 , φ in 2 , ∈ S and n ∈ N. Now, µ(t) := µ 1 (t) obeys the equations ∂ t Φ(t, ξ, µ in ) = (Kµ(t))(Φ(t, ξ, µ in )) and χ = r 1 .

Stability
Next, we are going to show that large classes of generic systems do admit at least stability of the synchronized state for large parameter regions. To begin, we use the following abbreviations: Let us now assume that a σ > 0 for all σ ∈ [M ] and for all σ ∈ [M ] and all α ∈ (−η, η) |r σ | , φ ∈ (−η, η). We will later impose conditions on κ > 0 and then choose η > 0 accordingly to (3.5). We want to remark that results from this section rely on the notation defined in Section 1.
Lemma 3.11. For all σ ∈ [M ] and any two particles φ σ , which is by Lemma 3.10 independent of t. Then, there exists a constant C > 0 such that whenever 0 < Ψ σ (t) < η for all σ ∈ [M ], they satisfẏ Proof. Let us consider a decomposition of the probability measures µ σ (t) into two measures µ inside σ (t) and µ outside and supp(µ inside σ (t)) ⊂ ((φ σ 1 (t), φ σ 2 (t)), supp(µ outside σ (t)) ⊂ S \ ((φ σ 1 (t), φ σ 2 (t)). Then, a calculation showṡ Here, the equality ( * ) was achieved by decomposing each measure µ i into its components according to (3.6) and rearranging terms such that every integrand with an integral running over at least one measure of the type µ outside is contained in the part "integrals over µ outside ". We can easily estimate these integrals from above by first combining the integrals into a single one running over µ outside σ (with the integrand still consisting of integrals) and then taking the supremum norm C of the integrand. As the total mass of µ σ (t) equals 1 and the µ σ -mass inside the interval (φ σ 1 (t), φ σ 2 (t)) is m inside σ , the mass outside this interval evaluates to 1 − m inside σ . Consequently, the terms summarized in "integrals over µ outside " can be bounded from above by C(1 − m inside σ ). The inequality ( * * ) is based on linear approximation and the fact that Ψ σ (t) < η. Lemma 3.12. Let µ ∈ P(S) be a probability measure on the circle, ξ 1 , ξ 2 ∈ S and m inside = (ξ1,ξ2) dµ(α).
Proof. Using (1.3a), we calculate We can not put the previous lemmas together and formulate our main theorem: Theorem 3.13. If the coupling functions g σ are continuously differentiable, i.e. g σ ∈ C 1 (S |r σ | × S), and they satisfy a σ > 0 for all σ ∈ [M ] then, the set of all-synchronized states S M is stable.

Almost Asymptotic Stability
One might now hope that although we do not have asymptotic stability, we can expect asymptotic stability of large classes of initial conditions as the family of steady states constructed above is a rather small part of phase space. Before stating with theorems regarding asymptotic stability, we need to introduce the concept of phase differences, as this concept becomes important in the subsequent proofs. Similarly to the original system (2.1), the system of phase differences describes the temporal evolution of oscillators, which can be grouped into populations, on the circle. Unlike the original system (2.1), in the system of phase differences, the position of the oscillators is not given in absolute coordinates but instead with respect to reference oscillators. The system of phase differences is given by Notation 3.14. Let ζ ∈ S. When using the notation m ζ , we refer to the function m ζ : Lemma 3.15. Let ζ 1 , . . . , ζ M ∈ S, µ in ∈ P(S) M , suppose that µ(t) solves the system (2.1), (3.2) and let Φ σ (t, ξ in σ , µ in ) be its mean-field characteristic flow. Now define and . Then, ν(t) and Ψ σ (t, ξ in σ , ν in ) solve the system (3.12), (3.13). Before proving this lemma, we remark that Φ σ (t, ζ σ , µ in ) can be seen as the position of reference oscillators we have talked about at the beginning of this section.
Proof of Lemma 3.15. It is easy to verify (3.12c): To check (3.12b), take a measurable set A ⊂ S and calculate To finally show (3.12a), we use the notation ). Then, a rather lengthy calculation confirms This completes the proof.
Remark 3.16. This Lemma is especially useful because the only operation used to create the measures ν σ (t) from the measures µ σ (t) is a rotation by Φ σ (t, ζ σ , µ in ) around the circle. Therefore, W 1 (S, ν σ (t)) = W 1 (S, µ σ (t)) and µ σ (t) → S if and only if ν σ (t) → S. is Lipschitz continuous with constant L 1 then, the coupling operator F σ satisfies Proof. Follows from (3.13) and Lemma 2.9. . Then, initial configurations in the space of densities µ in ∈ P ac (S) M , which are close enough to the all-synchronized state, converge to the all-synchronized state as t → ∞.
Remark 3.19. Note that the assumption on absolute continuity of the measures eliminates the counterexamples from Section 3.2.1 as only perturbations into measures with densities are allowed. This clearly shows, why studying the mean-field Vlasov-Fokker-Planck equation for densities is often easier in comparison to our goal of deriving directly the maximum information from the characteristic system.
Proof of Theorem 3.18. As the assumptions of this theorem include the assumptions of Theorem 3.13, we can choose ǫ V > 0 such that for any µ in ∈ B(S, ǫ V ) M , µ(t) ∈ B(S, ǫ U ) M for all t ≥ 0. To prove asymptotic stability in the space of absolutely continuous measures, let µ in ∈ B(S, ǫ V ) M ∩ P ac (S) M , with ǫ V specified later, and proceed analogously to the proof of Theorem 3.13 until we get to the point when φ σ 2 (t) − φ σ 1 (t) ≤ 2ζ for all σ ∈ [M ] and t ≥ 0. Next, we use Lemma 3.15 in order to switch to the system of phase differences with ζ σ = φ σ,in 1 = φ σ 1 (0). The reference oscillators in this system of phase differences are consequently given by φ σ for all t ≥ 0, ν σ (t, dγ) Therefore, a computation similar to the one in the proof of Theorem 3.13 shows so the solution ν(t) always stays close to the all-synchronized state located at the origin. In the system of phase differences, individual particles then follow the flow However, by (3.15) and Lemma 3.17, (3.16) can be rewritten in the forṁ for a small perturbation p σ ∈ C 1 (S × R). Next, we fix ǫ U and with that also ǫ V such that for all ν 1 , . . . , ν M ∈ B(δ 0 , ǫ U ), p σ (·, t) C 1 is small enough such that the flow induced by the dynamical system (3.16) is equivalent to the one induced byΨ σ = (F σ δ M 0 )(Ψ σ ) for all σ ∈ [M ]. So we consider (3.17) as a perturbed one-dimensional autonomous ODE. The existence of such an ǫ U is guaranteed by Lemma 3.17. Specifically, this lemma states that p σ (·, t) C 1 ≤ (4L |r σ | + 2L + 4L 1 |r σ | + 2L 1 )ǫ U =: ǫ σ f , uniformly in t. A particular choice of ǫ U can be constructed as follows: Asĝ ′ σ (0) = a σ > 0 and there are only two roots with non-vanishing derivative on the circle, b σ < 0. Further, let us write η σ 1 , η σ 2 > 0 for radii of intervals such that inf α∈(−η σ Now, ǫ U can be chosen such that for all f σ ∈ C 1 (S) with f σ C 1 < ǫ σ f the following criteria are satisfied: . While (C1) ensures that F σ ν(t) has no zeros away from the roots −ψ σ 0 and 0 for all t ≥ 0, (C2) guarantees that F σ ν(t) is strictly monotonic in the two neighborhoods around the roots. This monotonicity also causes the existence of at most one zero of F σ ν(t) near the two roots. Even though the two zeros of F σ ν(t) may be varying over time, (C2) ensures that the flow of (3.16) is still exponentially contracting in (−η σ 1 , η σ 1 ) and exponentially expanding in (−ψ 0 σ − η σ 2 , −ψ 0 σ + η σ 2 ). Thus, at least one of two distinct test particles starting in (−ψ 0 σ − η σ 2 , −ψ 0 σ + η σ 2 ) leaves this region and eventually ends up in (−η σ 1 , η σ 1 ). Since Ψ σ (t, 0, ν in ) = 0 for all t ≥ 0 and the contracting property of Ψ σ (t, 0, ν in ) around 0, the particle even converges to the origin at 0. So there exists only one trajectory, which starts at an arbitrary point Ψ in σ ∈ S, that does not converge to 0. By assumption, ν in σ ({Ψ in σ }) = 0 and hence all the mass concentrates around 0. Therefore, W 1 (δ 0 , ν σ (t)) → 0 as t → ∞. Because this holds true for all σ ∈ [M ], the all-synchronized state is asymptotically stable for absolutely continuous perturbations.
Remark 3.20. Theorems 3.13 and 3.18 per se only apply to perturbations in all populations. However, if we exemplarily want to analyze the stability of DSSD in a network of M = 4 populations, Proposition 3.2 allows us to reduce the system of four populations to equations describing only the evolution of population #2 and #3 while we keep population #1 and #4 fixed in splay state. Applying Theorems 3.13 and 3.18 consequently yields criteria for the (asymptotic) stability of DSSD with respect to perturbations in the second and third population.

A Review of Linear Stability with Pairwise Coupling in One Population
In the easiest case we consider only one population and no higher-order coupling. Then, the velocity field (3.2) is given by A well known way [45,48] for analyzing stability of the splay state is to look at the mean-field (or continuity) equation which describes the evolution of the density ρ(t, φ). Next, one typically inserts the ansatz ρ(t, φ) = 1 2π + ǫη(t, φ) into the continuity equation and collects terms of order ǫ. Assuming Fourier representations where c.c. denotes the complex conjugate of the previous term, one can derive differential equations for the evolution of the coefficients c k (t): Fortunately, these equations are uncoupled and linear stability of the splay state can thus simply be infered if Im(a k ) > 0 for all k ≥ 1. In other words, when writing the coupling function g(γ) as linear combinations of sin(γ), cos(γ) and trigonometric monomes of higher order, the prefactors of sin(γ), sin(2γ), . . . have to be negative. A similar analysis yields the linear instability of the splay state if Im(a k ) < 0 for at least one k.

Non-Pairwise Coupling in Multi-Population Systems
Let us now consider the more general case of multi-population systems and higher-order interactions.
Assuming the initial measures to be represented by densities ρ in σ (φ), the velocity field is given by and the densities ρ σ (t, φ) solve the continuity equation (2.8). Here, ρ (r σ ) is the shorthand notation for In this section, we extend the formal calculations from Section 3.3.1 to the case of multi-population systems and higher-order interactions. Such a formal derivation of criteria for linear stability of the all-splay state can be done with the same techniques as those used in Section 3.3.1. As in this section, we therefore consider a small perturbation around the all-splay state, i.e., with Fourier decompositions c σ k (t)e ikφ + c.c. (3.19) Further, the coupling functions g σ : S |r σ | × S → R are supposed to be given in terms of its Fourier expansion as well: The requirement a σ 0,0 = 0 is not really a limitation as possible non-zero values of a σ 0,0 = 0 can be absorbed into ω σ .
Given these representations, we formally insert (3.18) into the continuity equation (2.8) to obtain ∂ ∂t with the usual abbreviation Collecting terms of order O(ǫ) yields Let us now interchange the sums with the integrals and evaluate each of the summands individually to obtain The same computations holds true if we replace η (r σ j ) (t, α j ) with η (r σ j ) (t, β j ). Therefore, the sum in (3.21) vanishes and as we see next, the only term that does not vanish inside the brackets of this equation is η σ (t, γ). Multiplying it with the coupling function g σ , integrating and subsequently simplifying yields Combining these two results with (3.21) we get Thus, after having taken the derivative and having collected e ikφ -terms, it is easy to see that c σ k (t) obeys the differential equation Therefore, small perturbations of population σ in direction of e ikφ + c.c. with k ≥ 1 decay on a linear level if Re(−(ā σ 0,k + ω σ )ik) = −k Im(a σ 0,k ) < 0. Similarly, they grow if −k Im(a σ 0,k ) > 0. However, it is important to note that the equations (3.22) are only based on a formal derivation. Assuming nonetheless that our formal calculations can be made rigorous in this way, we can summarize our results by claiming linear stability of the all-splay state if Im(a σ 0,k ) > 0 for all σ ∈ [M ], k = 1, 2, . . . and linear instability if Im(a σ 0,k ) < 0 for one σ ∈ [M ] and one k = 1, 2, . . . . Remark 3.21. There are several challenges when trying to obtain rigorous (linear) stability results. First, we have to rigorously linearize by constructing a suitable function space in which the operator F defined by is Fréchet differentiable. Then, we have to check that the formal calculation above holds within this function space, and that we have described the spectrum completely. Finally, one has to invoke a suitable result that linear stability entails local nonlinear stability. Carrying out this full stability analysis program is beyond the scope of the current work.

Mean-Field Dynamics of Phase Oscillator Networks
In this section, we give several examples to illustrate the theory and results we have developed so far.

The Kuramoto Model for Identical Oscillators
In the easiest case, M = 1 and s 1 = {}, so the function G 1 in (2.2) is mapping only from S to R. The measure µ 1 (t) then simply gets transported along the time-independent velocity field However, this case is not really interesting, which is why we now consider the case M = 1, s = s 1 = (1), G 1 (α, φ) = sin(α − φ). For an initial measure µ in = µ in 1 ∈ P(S), the characteristic system (2.1)-(2.2) simplifies to with the coupling function where ω ∈ R is the common oscillator frequency of the single population. So we have recovered the classical characteristic system of the Kuramoto model for identical oscillators. Indeed, if one assumes that the initial measure has the form of an empirical measure given by for some N ∈ N, then it is easy to see that the solution µ(t) is of the form µ(t) = 1 The equations (4.2) are the classical finite-dimensional Kuramoto model [31]. In fact, one may prove that as N → ∞, then the evolution of the empirical measures due to (4.2) is well-approximated by a mean-field limit [39] as described in Section 2.2. This model can be extended by replacing the sinusoidal coupling function G 1 by a more general coupling function G 1 (α, φ) = f (α − φ). The velocity field in this generalized Kuramoto model for identical oscillators is then given by Theorem 3.13 yields the stability of the synchronized state in this system if f ∈ C 1 (S) and f ′ (0) > 0. By Theorem 3.18, the synchronized state is asymptotically stable in the space of densities if furthermore the function ψ → f (ψ) − f (0) has only one root with non-vanishing derivative around the circle except 0 and f ′ is Lipschitz continuous.

One Population with Higher-Order Interactions
Let us now consider a system which still consists of only one population but involves higher-order interactions. Specifically, we reconsider system (1.2): Here, N denotes the amount of discrete oscillators, i ∈ [N ] and K 1 , K 2 , K 3 ∈ R. Putting this system into our framework (2.1)-(2.2), the coupling function G := G 1 : S 3 × S → R is given by and the multi-index s 1 is trivially given by s 1 = (1, 1, 1). Consequently, by Theorem 2.14, there even exists a measure valued solution of the mean-field limit of the system (4.3). However, to put this system into the more restrictive form of (3.2), we need to assume K 2 = 0, as we have also assumed in Section 2.4. In this case, the function g := g 1 : S × S → R from (3.2) reads as g(α, γ) = K 1 sin(γ) + K 3 sin(α + γ), (4.4) with multi-index r 1 = (1). Note that g(0, γ) = (K 1 + K 3 ) sin(γ). Theorem 3.13 thus tells us that the synchronized state in the system (2.1),(3.2),(4.4) is stable if K 1 + K 3 > 0. Furthermore, in this case, by Theorem 3.18, it is even asymptotically stable in the space of densities.
Example 4.1 (Example 2.21 revisited). Again, let K 1 = 1 and K 3 = −4. Then, K 1 + K 3 = −3 and thus, the synchronized state is not stable. This is consistent with Example 2.21, since there we found r = 1/2 to be attractive. In particular, almost all initial conditions globally converge to r = 1/2 and so it makes sense that the synchronized state is not stable.
Note that the network interactions in (4.3) are quite specific: The coupling functions are purely sinusoidal and thus the authors use the Ott-Antonsen reduction to understand the system dynamics [43]. Note that our analysis here is not limited to such networks but applies also to phase oscillator networks with general (higher-order) interactions where the reduction methods cease to apply. These networks arise naturally, for example, through phase reductions in oscillator networks that generically contain multiple harmonics; cf. [4,9,34].

Multiple Coupled Populations with Higher-Order Interactions
In [8] networks of M = 3 finite phase oscillator populations coupled by higher-order interactions have been considered. In particular, the main equations from [8] are given bẏ where the index σ ±1 for the population has to be understood modulo M if σ ±1 ∈ {1, 2, 3}. Furthermore, φ σ,k refers to the phase of the kth oscillator in population σ for k = 1, . . . , N , h 2 : S → R is a Lipschitzcontinuous intra-population coupling function and with a Lipschitz-continuous inter-population coupling function h 4 : S → R. It can be seen that if one population, say the first, is initially synchronized, i.e., φ 1,1 (0) = · · · = φ 1,N (0) then it is also synchronized at later times. Similarly, if the the oscillators of the first population are initially in splay state, they exhibit this property also at later times. If two populations are fixed to be either synchronized or in splay state, the oscillators in the remaining free populationσ evolve according tȯ where the coupling function g is made up of h 2 and h 4 and depends on the state in which the fixed populations are locked. More interestingly, by choosing appropriate coupling functions h 2 , h 4 and coupling constants K + , K − , one can achieve that the (de)synchronization of one population can cause another population to start synchronizing or desynchronizing. In particular, existence of a heteroclinic cycle for small populations was shown in [8]. Within these heteroclinic cycle populations alternatingly synchronize and desynchronize. Our numerical simulations have confirmed that such a cycle continues to exist if the size of the populations grows to infinity. With little change to the system (4.5),(4.6), this system can be written in the form of (2.1),(2.2). In fact, when summing over all j from 1 to N in (4.5) instead of over all j except j = k, the coupling in this system is of the form (2.2) when s 1 = (3, 3, 2, 2, 1), s 2 = (1, 1, 3, 3, 2), s 3 = (2, 2, 1, 1, 3) and the coupling functions are defined by G 1 (α, φ) = G 2 (α, φ) = G 3 (α, φ) := G(α, φ) with The velocity field (2.2) then evaluates to Our results allow to analyze the stability of invariant sets of these networks in the mean-field limit.
Note that all of the multi-indices are of the special form (3.1). Thus, with r 1 = (3, 1), r 2 = (1, 3), r 3 = (2, 1) this system can also be put into the form (3.2). Then, the coupling functions are given by Having put the system into the form (3.2) allows us to apply results from Section 3. For example Theorem 3.13 yields the stability of the all-synchronized state if the function Let us now try to investigate the stability of SDD with respect to perturbations in the first population. Unfortunately, we cannot apply Theorem 3.13 immediately but we have to do some preparatory steps first. So we first choose µ in 2 (A) = µ in 3 (A) = 1 2π λ S (A). Then, Proposition 3.3 causes the second and third population to stay in splay state for all t ≥ 0. The velocity field according to which µ 1 (t) is transported is given by Therefore, the dynamics in µ 1 DD can be described by a single coupling function g(γ) := h 2 (γ) + K + − K − 2π S h 4 (α, γ)dα.
Only now, we can apply Theorem 3.13 to see that SDD is stable with respect to perturbations in the first population if the functionĝ satisfiesĝ ′ (0) > 0. In order to investigate the linear stability of DDD we need to assume Fourier expansions ξ k e iγk + c.c., ζ l,k e iαl e iγk + c.c.
By inserting them into the representation of the coupling function g, we see that a l,k e iα1l1 e iα2l2 e iφk + c.c.
for Fourier coefficients a l,k that satisfy a 0,k = ξ k − K − ζ 0,k + K + ζ 0,k . By the results obtained in Section 3.3.2, the all-splay state is linearly stable if Im(ξ k + (K + − K − )ζ 0,k ) > 0 for all k = 1, 2, . . . and linearly unstable if one of these coefficients is less than 0. These results give necessary conditions for the emergence of heteroclinic cycles involving the invariant sets SDD, SSD, . . . to exist not only in networks of finitely many oscillators but also in the mean-field limit of these systems. Note that a similar analysis is possible for the mean-field limit of networks that consist of M = 4 coupled oscillator populations that support heteroclinic networks with multiple cycles in [12].

Discussion and Outlook
In this paper, we have proposed a new general framework for the evolution of coupled phase oscillator populations with higher-order coupling. First, we provided a general solution theory. We clarified existence and uniqueness of weakly continuous solutions to the characteristic system and justified the mean-field limit. Then, we studied some dynamical properties of the characteristic system. We showed that the subspaces, which are characterized by one or more populations being either in synchronized or in splay state, are dynamically invariant. Next, we proved the stability of the all-synchronized state under concrete conditions on the coupling functions. Finally, we also analyzed developed linear stability analysis for the splay state via the mean-field limit equation.
Although we have provided the general mathematical foundations for studying large-scale multipopulation oscillator networks with higher-order coupling, there are still many open questions for future work. Here, we considered the case of populations with identical oscillators. This means that in the mean-field limit there are atomic measures that are invariant under the flow. There are two ways to break this degeneracy. First, one can assume that the intrinsic frequencies of the oscillators follow a distribution with a density as in the classical Kuramoto model; cf. [10]. Second, adding noise to the evolution leads to a diffusive terms in the Fokker-Planck equation [50]. In either case, the synchronized phase configuration is not invariant anymore and deforms to a near-synchronous stationary solution. Insights into how the stability properties derived here change through these perturbations would be desirable. Unlike the synchronized state, the splay state would stay invariant for non-identical oscillators but investigating its linear stability would be more complicated due to the frequency dependence. One expects bifurcation structures to be affected generically by higher-order coupling, e.g., being able to change super-to sub-critical transitions [29,43].
If the network interactions contain just a single harmonic, the Watanabe-Strogatz reduction applies. Its application in the mean-field limit has so far been heuristic and one typically assumes the existence of densities [42] and it would be interesting to understand Watanabe-Strogatz theory for the characteristic equation (2.1) that describes the evolution of general measures. Together with nonidentical frequencies within a populations, this would be the first step towards a rigorous description of Ott-Antonsen theory in a measure theoretic sense; see also [20].
While we discussed network dynamical systems with higher-order interactions, we did not assign an explicit algebraic structure to the dynamical equations. Recently, dynamical systems on higherorder networks-whether hypergraphs or simplicial complexes-have attracted attention. While the assignment of higher-order network structure may not be unique, this perspective has its advantages: It naturally leads to limiting systems involving hypergraph variants of graphons [37,36,30], or more generally hypergraph variants of graphops [5,28,23]. In this context, one could also aim to link the stability analysis via the Vlasov-Fokker-Planck equation for hypergraphs better with direct methods on the level of the finite-dimensional ODEs for hypergraphs such as master stability functions [38,6]. Since the uncoupled system (A.2) has a supercritical hopf bifurcation at λ = 0 and the equilibrium point at the origin is stable for λ < 0, we expect the coupled system to have an invariant attracting torus whenever λ > 0 and ǫ are close to 0. In fact, a theorem from [4] shows exactly that and furthermore states that the flow on the invariant torus can be approximated by a higher-order coupled oscillator system.
Specifically, as stated in [4], over a longer timescale a system of the form better approximates the dynamics on the invariant torus than a system without higher-order coupling.