Hopf bifurcations of twisted states in phase oscillators rings with nonpairwise higher-order interactions

Synchronization is an essential collective phenomenon in networks of interacting oscillators. Twisted states are rotating wave solutions in ring networks where the oscillator phases wrap around the circle in a linear fashion. Here, we analyze Hopf bifurcations of twisted states in ring networks of phase oscillators with nonpairwise higher-order interactions. Hopf bifurcations give rise to quasiperiodic solutions that move along the oscillator ring at nontrivial speed. Because of the higher-order interactions, these emerging solutions may be stable. Using the Ott–Antonsen approach, we continue the emergent solution branches which approach anti-phase type solutions (where oscillators form two clusters whose phase is π apart) as well as twisted states with a different winding number.


Introduction
Coupled phase oscillators on a network provide essential models to understand synchronization phenomena [1].Apart from global phase synchrony, where the phase of all oscillators coincides, twisted states in nonlocally coupled networks have attracted attention [2].For example, consider N Kuramoto oscillators on a ring network, whose phase θ k ∈ T := R/2πZ evolves according to θk := d dt where ω ∈ R is the intrinsic oscillator frequency and G : R → R is a smooth one-periodic function that determines the strength of the oscillator interaction depending on their relative position on the ring.For q ∈ N, the q-twisted states with β ∈ T, Ω ∈ R, k = 1, . . ., N are periodic solutions of (1); see Figure 1(a) for an example.These solutions are also known as rotating wave solutions or, for q = 1, as splay phase configurations.While the stability of such solutions has been analyzed explicitly [2,3], the specific form of Kuramoto phase coupling (a single harmonic without phase shift) imposes a gradient structure, which prevents the emergence of bifurcations to periodic solutions.
Richer dynamics and bifurcation behavior are possible for more generic phase interactions that one would expect from phase reductions [4].Indeed, phase reductions also give rise to nonpairwise phase interaction terms [5,6,7], which can give rise to Hopf bifurcations of splay phase solutions and more complicated dynamical behavior [8].A generalization of the Kuramoto model (1) for ring-like networks to include nonpairwise interaction terms is where K 2 ∈ R determines the strength of the pairwise interactions with a phase shift α 2 ∈ T and K 3 ∈ R determines the strength of nonlinear interactions between three phase variables.Such phase oscillator networks with "higher-order" nonpairwise interactions have their intrinsic interest [9,10] as dynamical systems on (weighted) hypergraphs where the pairwise phase interaction function sin(θ j −θ k +α 2 ) corresponds to interactions along edges and the triplet phase interaction function sin(2θ j − θ l − θ k + α 3 ) corresponds to interactions along hyperedges.The typical setup to analyze q-twisted states is the continuum limit of N → ∞ oscillators; see already [2].In the continuum limit the phase of oscillator x ∈ S := [0, 1]/(1 ∼ 0)-the unit interval with end points identified-at time t is given by a phase θ(x, t) = θ x (t) that, for the ring network (3), evolves according to θx = ω + K 2 S G(x − y) sin(θ y − θ x + α 2 ) dy + K 3 S S G(x − y)G(x − z) sin(2θ y − θ z − θ x + α 3 ) dzdy.(4) Now q-twisted states Θ q (x, t) := 2πqx + Ωt + β (5) are periodic solutions with a smooth phase profile (in x).Their stability has been analyzed for networks with Kuramoto coupling (α 2 = 0) without higher-order interactions [2,11] and more recently also with higher-order interaction [12].
Here we analyze Hopf bifurcations of q-twisted states in phase oscillator networks (4) with higher-order interactions.To identify Hopf bifurcation points, we linearize the system at q-twisted states and find conditions for a complex conjugate pair of eigenvalues to cross the imaginary axis; we focus primarily on 1-twisted states.Computing higher-order derivatives allows to determine whether the Hopf bifurcation is subcritical or supercritical and estimate amplitude and period of the bifurcating solution.The emergent periodic solutions are quasiperiodic solutions with nontrivial rotation along the spatial domain S-see Figure 1(b) for an example of the corresponding solution in a finite network.Because of the higher-order interactions, these emergent solutions can be stable.To continue the solutions further from the bifurcation point, we consider the dynamics on the Ott-Antonsen manifold [13,14].Adapting recent approaches to continue periodic solutions on the Ott-Antonsen manifold [15,16], we compute bifurcation branches for varying parameters; cf. Figure 1(c) for an example of branch of periodic solutions bifurcating from a 1-twisted state.Further bifurcations involve (traveling and stationary) antiphase solutions and some branches appear to come close to −1-twisted solution.
This paper is organized as follows.In Section 2 we specify the dynamical equations as well as relevant parameters and discuss them in the context of recent related work.In Section 3, we linearize the equations and give the eigenvalues that determine the bifurcations as well as the higher-order derivatives that determine the bifurcation type.We consider the system on the Ott-Antonsen manifold in Section 4 and outline the continuation technique.In Section 5 we describe the bifurcations that arise and conclude in Section 6 with some remarks.

Twisted States in Ring Networks with Higher-Order Interactions
In the following we analyze q-twisted states in the continuum limit equation (4).Without loss of generality we will henceforth set ω = 0 by exploiting the phase-shift symmetry θ x → θ x +β, β ∈ T, that shifts all oscillator phases by a constant phase angle β.Thus, the dynamical equations are and the q-twisted states (5) are actually equilibria relative to the phase-shift symmetry action.The 0-twisted state corresponds to full phase synchrony in the system and the ±1-twisted states to the classical (anti)splay phase configuration.
The network has a ring structure (with orientation) since the system has a spatial rotational symmetry, where s ∈ S acts on S by s : x → x + s.Specifically, we consider a network with a coupling kernel function (7) G(x) = 1 + A cos(2πx) + B sin(2πx).
with only the first nontrivial harmonic being present; see also [17,18,19].A nonzero parameter B breaks the reflectional symmetry of the ring network: For B = 0, the reflection x → −x is also a symmetry of the system.The higher-order triplet interactions are determined by the triplet coupling kernel that is a product of the same coupling kernel function G that determines the pairwise interactions.Such a product structure is natural if the phase equations are obtained through a phase reduction [7].However, it is distinct from other triplet interaction kernels considered in the literature in the context of ring networks.Specifically, in [12] the authors considered a generalized top-hat coupling kernel W (x, y, z) = W r (z + y − 2x) where W r (v) = 1 if min(|v|, 1 − |v|) ≤ r and W r (v) = 0 otherwise.This interaction function lacks a product structure but generalizes coupling with a finite coupling range considered in the context of twisted states on rings with pairwise coupling [3].By contrast, with the triplet coupling kernel (8) and the kernel function (7), the interactions in (6) only consists of finitely many Fourier modes; this facilitates analytic computations, as we will see in Section 3.
As for the phase interaction function in (6), it only depends on the first harmonic of the state of the oscillator at x. Hence, the system can be reduced to the Ott-Antonsen manifold [13,14] as for globally coupled networks [20,21].

Bifurcations of Twisted States
In this section, we study the bifurcation that occurs when a q-twisted state gains or loses its stability under variation of parameters.As noted above, the system (6) has a T-symmetry that maps a solution θ x to a solution θ x + β for a given constant β ∈ T. Thus, the linearization of the right-hand side of (6) always has a zero eigenvalue, which makes a rigorous bifurcation analysis tedious.In order to avoid this zero eigenvalue, we change to the system of phase differences and define Ψ x (t) := θ x (t) − θ 0 (t).Then, the function Ψ x (t) satisfies and Ψ 0 (t) = 0 (10) for all times t.The process of transitioning from (6) to the system of phase differences (9) reduces the T-symmetry.In particular, the T-symmetry is present in the system (6) but not in (9), since every function in the system of phase differences has to satisfy Ψ 0 (t) = 0 and thus cannot be shifted by a constant.After this symmetry reduction, the q-twisted states (5) are represented by the function Ψ q x = Ψ q (x) with Ψ q (x) := 2πqx, (11) which does not depend on the time t, satisfies Ψ q (0) = 0, and cannot be perturbed by a constant function.Consequently, when linearizing the righthand side of ( 9), which we denote by G, around a q-twisted state (11) there is no trivial zero eigenvalue.To make this precise, we define a function G q : X × P → X as where q ∈ Z is the winding number, p = (K 2 , K 3 , A, B, α 2 , α 3 ) summarizes all parameters and X = H 1 0 (S, R) is the space of once weakly differentiable functions on S with zero boundary conditions.Since H 1 (S, R) ⊂ C(S), these boundary conditions can be imposed in the classical sense.The function G q consequently gives the local behavior of the right-hand side of (9) around the q-twisted state, v ∈ X can be seen as a perturbation of the twisted state, and the condition v(0) = 0 ensures (10).Conducting a bifurcation analysis of (9) instead of (6) simplifies the setting.
3.1.Linearization.The bifurcation analysis of twisted states is based on a stability analysis, which can be conducted using the eigenvalues of the linearization of the right-hand side around the twisted state.More precisely, we linearize G q (v, p) around v = 0, consider this linearization as an operator D v G q (0, p) from X to itself and determine the eigenvalues of this operator.
To obtain D v G q (0, p), we take a function η • ∈ X and h ∈ R, and calculate Using the definition of G q we obtain • cos(2Ψ q y + Ψ q z + α 3 ) dzdy.Even though this computation was very formal, it can be rigorously shown that D v G q (0, p), as calculated here, is indeed the Fréchet derivative of G q (v, p) at v = 0, see [12].
Henceforth we focus on the twisted state with the smallest winding number, namely q = 1; the case q = −1 is analogous.We note that the functions u k (x) = sin(2πkx) and w k (x) = 1 − cos(2πkx) for k ≥ 1 form a Schauder basis of X, see [12,Lemma B.1].Consequently, we evaluate D v G 1 (0, p 0 ) on these basis functions.For k = 1, this yields which implies a complex conjugated pair of eigenvalues where Whenever Im (λ + ℓ ) ̸ = 0 but Re (λ + ℓ ) = 0, for some ℓ, we expect a Hopf bifurcation to occur under variation of parameters.The transverse stability of this bifurcation is then determined by the other eigenvalues.If there exists k ∈ N with k ̸ = ℓ such that Re (λ + k ) > 0, all equilibria and periodic orbits that emanate from the Hopf bifurcation are not transversely stable.If on the other hand Re (λ + k ) < 0 for all k ̸ = ℓ, the equilibria and periodic orbits are at least transversely stable.As we are particularly interested in transversely stable bifurcations, we assume from now on that this is the case.As one can see in Figure 2, there are parameter regions where the critical eigenvalue is attained for ℓ = 1 and other regions where ℓ = 2 is the dominant eigenvalue.Here, a Hopf bifurcation can occur.However, there also exist parameter regions, where ℓ = 4 is the critical eigenvalue.Since Im (λ ± 4 ) = 0, there is no generic Hopf bifurcation for these parameter values.
3.2.Bifurcations.Having investigated the linear stability of the 1-twisted state, we now analyze the periodic solutions that emanate from the 1-twisted state in a Hopf bifurcation.We use the notation C T (R, X) for all continuous T -periodic functions with values in X.Moreover, we assume for simplicity that B is the main bifurcation parameter, i.e., we fix all other parameters and vary only B to initiate the bifurcation, which then occurs at some value B 0 .We can then write the eigenvalues of linearization around the 1-twisted state as functions of B. At the critical value B 0 we have Re (λ ± ℓ (B 0 )) = 0 and we additionally assume that Re (λ ± k (B 0 )) < 0 for all k ̸ = ℓ and Im (λ ± ℓ (B 0 )) ̸ = 0.In particular, we also denote λ ℓ (B) for the critical eigenvalue that has positive imaginary part, i.e., λ ℓ (B) : Similarly, we denote v ℓ for the eigenfunction that corresponds to this critical eigenvalue.
A generic theorem, that relies on some technical assumptions [22, Theorem I.8.2], then guarantees that periodic solutions bifurcate from the 1twisted state when B passes through B 0 .In particular, there is a continuously differentiable function κ : U → R, where U is a small neighborhood of 0 ∈ R, such that Im (λ ℓ (B 0 )) = κ(0) > 0.Moreover, there is a continuous Figure 3. Illustration of the Hopf bifurcation.Note that the vertical axis represents the space C 2π/κ(r) , which is different from classical bifurcation diagrams.Thus, there is also no flow around the paraboloid, but each point on the paraboloid represents a solution.Each point on the black parabola can be reached via the curve r → (B(r), v(r)).Every other periodic solution can be obtained by a phase shift (B(r), S τ v(r)) for a suitable τ ∈ R.
Additionally to the first derivatives of G 1 , second and third derivatives are required to approximate the curve (v(r), B(r)), see [22, Theorem I.9.1].In particular, the curve of 2π/κ(r)-periodic solutions can be approximated as d dr v(r) where ζ ∈ C can be computed from the first, second and third derivative of the right-hand side G at the 1-twisted state; see Appendix A. Further, Re ( d dB λ ℓ (B)| B=B 0 ) is the speed with which the real part of the critical eigenvalue passes through zero.Basically, (13a) helps to approximate the amplitude and profile of bifurcating periodic solutions, (13b) determines their period, and (13c) connects the parameter B with r and thereby determines for which parameter value of B these periodic solutions occur.In particular, if (13c) is positive, these periodic solutions are only existent when B > B 0 .Conversely, if (13c) is negative, they exist for B < B 0 .Using the principle of exchange of stability of the equilibrium with the periodic orbits in a Hopf bifurcation, one can even determine the stability of the periodic orbits.Specifically, if Re (ζ) > 0 one can consider two cases depending on the sign of Re ( d dB λ ℓ (B)| B=B 0 ).If it is positive, then the 1-twisted state is unstable for B > B 0 .Moreover, (13c) is positive and thus the periodic solutions only exist when B > B 0 .Since the twisted states are unstable in that parameter regime, the periodic solutions have to be stable.If on the other hand Re ( d dB λ ℓ (B)| B=B 0 ) is negative, the twisted state is unstable for B < B 0 .Since (13c) is then also negative, the periodic solutions also exist in the same parameter region and are stable.In both cases the bifurcation is supercritical.Repeating the argument for Re (ζ) < 0 shows that the bifurcation is then subcritical, meaning that the periodic solutions that emanate from the bifurcation are then unstable.In conclusion, the sign of Re (ζ) can be used to distinguish between sub-and supercritical Hopf bifurcations.

Continuation on the Ott-Antonsen Manifold
It is easy to verify that the dynamics described by Eq. ( 6) is equivalent to the dynamics of the Ott-Antonsen equation with a convolution-type integral operator after insertion of the ansatz z(x, t) = e iΘ(x,t) .
Eq. ( 14) is useful for two reasons.First, it can be used to perform a linear stability analysis of twisted states in an alternative way; see Appendix B. Second, this equation can be used to compute the solution branches emanating from the Hopf bifurcations found above.In this section, we outline how to compute these solution branches; the approach is based on recent results presented in [15] and adapted to our setting.
4.1.Self-consistency equation for traveling nonuniformly twisted waves.All periodic solutions emerging at the Hopf bifurcations of 1-twisted states have form of quasiperiodically evolving synchronization patterns with some Θ 0 (x) and s, Ω ∈ R. For Eq. ( 14) they correspond to traveling wave solutions of the form where |a(x)| = 1 for all x ∈ S. Such solutions can be efficiently computed, using the self-consistency equation derived below.By inserting ansatz (15) into Eq.( 14) we find that the profile function a(x) in ( 15) is a periodic solution of the integro-differential equation If s ̸ = 0, Eq. ( 16) is equivalent to the complex Riccati equation contains the network coupling terms.Rather than finding the unknown function a(x), our strategy is to determine the corresponding mean field w(x) using a self-consistency equation.First, we recall some facts about Eq. ( 17), which were previously proved in [15, Section 2] by one of the authors of this paper.The most important fact is that for an arbitrary periodic function w(x) and an arbitrary real coefficient Ω/s, the complex Riccati equation ( 17) usually has two periodic solutions.The initial conditions of these solutions are determined by the fixed points of the Poincaré map of Eq. (17).Moreover, due to the special structure of Eq. ( 17), its Poincaré map has the form of a Möbius transformation M(z) = e iθ (z + b) bz + 1 with some b ∈ {z ∈ C : |z| < 1} and θ ∈ 2πS determined by the choice of w(x) and Ω/s.If |b| > | sin(θ/2)|, then the both fixed points of M(z) lie on the unit circle and therefore Eq. ( 17) has two periodic solutions satisfying |a(x)| = 1.One of these solutions is asymptotically stable, while the other is asymptotically unstable.In contrast, if |b| < | sin(θ/2)|, then the fixed points of M(z) lie inside or outside of the unit circle.More specifically, in this case, Eq. ( 17) has one solution that satisfies the inequality |a(x)| < 1 and another solution that satisfies the inequality |a(x)| > 1.
Suppose that the Ott-Antonsen equation ( 14) has a stable traveling wave solution of the form (15) with |a(x)| = 1.For this solution, we can calculate its mean field w(x) through (18).Then, considering w(x) and Ω/s as given, we can try to solve the periodic boundary value problem for Eq. ( 17) with the constraint |a(x)| = 1.Above we have shown that such a problem can either have two solutions (one stable and one unstable) or none.Moreover, in the first case, depending on the sign of the speed s, only one of the two solutions can be relevant for a stable traveling wave (15).Indeed, recalling that in (15) we use the moving coordinate frame ξ = x − st, we easily conclude that if s > 0 (s < 0) then a stable traveling wave (15) corresponds to an unstable (stable) solution of Eq. ( 17).Altogether these facts allow us to say that using the periodic boundary value problem for Eq. ( 17) we can always reconstruct the profile a(x) of a stable traveling wave (15) if the corresponding mean field w(x) and the ratio Ω/s are known.Denoting the resulting solution operator by U we can write (19) a(x) = U(w(x), Ω/s).
Obviously, the last expression agrees with the definition of w(x) if and only if which is an integral self-consistency equation for w(x).
To complete the definition of Eq. ( 20), we show a simple and fast way to calculate the operator U.The justification of this method, consisting of five steps, is given in [15, Section 2]: (i) Given w(x) and Ω/s, solve Eq. ( 17) starting from the initial condition a(0) = 1, and denote ζ 1 = a(1).
(iii) Calculate the coefficients b and e iθ of the Möbius transformation M(z) representing the Poincaré map of Eq. ( 17), (Note that by construction |ζ 1 | = 1.) (iv) Check that |b| > | sin(θ/2)| (otherwise the operator U is not welldefined) and then calculate the initial value a * of the periodic solution of interest, namely (v) Solve Eq. ( 17), starting from the initial condition a(0) = a * .This yields the periodic solution of Eq. ( 17).
Remark 4.1.If in step (iv) of the above algorithm we choose the formula of a * with '−' for s < 0 and the formula with '+' for s > 0, we obtain another solution operator for the periodic boundary value problem associated with Eq. (17).But this operator, by construction, gives profile functions a(x) of traveling waves (15), which are unstable with respect to Eq. ( 14).4.2.Algebraic self-consistency equation.The self-consistency equation ( 20) is a nonlinear integral equation that can be difficult to solve.However, in the case of coupling kernel (7), it can be reduced to a finite-dimensional nonlinear system.To see this, let us denote Then using trigonometric identities we can write where is the usual L 2 inner product.Moreover, for complex numbers a 0 , This relation together with Eq. ( 20) implies ( 21) with some complex coefficients ŵk .Inserting ansatz (21) into Eq.( 20) and equating the terms proportional to different ψ k separately, we reformulate Eq. ( 20) as a system of five complex equations with and System (22) needs to be solved with respect to s, Ω and ŵk , k = 1, . . ., 5.But it is obviously underdetermined for this.The problem can be resolved by recalling that the original Eq. ( 14) has two continuous symmetries.Therefore, we may add two pinning conditions to (22).For the sake of convenience, we choose these pinning conditions in the form (23) Im ŵ2 = Re ŵ3 = 0. Now, if we find a solution of the extended system ( 22), (23), we can use formula (21) to calculate w(x) and then (19) to calculate the corresponding a(x).Altogether, two scalars s and Ω and the profile function a(x) allow us to determine the traveling wave solution (15) of Eq. ( 14).Finally, if we want to show a typical snapshot of the corresponding solution of Eq. ( 6), we can use Θ(x) = arg a(x).

Continuation of Periodic Traveling Solutions
The approach outlined in the previous section now allows to continue periodic solutions that emanate from a Hopf bifurcation of a splay solution.The linear stability analysis of the 1-twisted state (cf.Section 3) indicates the location of the Hopf bifurcations; cf.Fig. 2 for the stability diagram of the 1twisted state in the (B, K 3 )-plane for fixed A = 0.9, Where the Hopf bifurcation is supercritical new types of stable time-dependent synchronization patterns can emerge.We tested this theoretical prediction in numerical simulations for finite networks (3) of N = 512 phase oscillators and found that the new solutions take the form of spatially modulated traveling waves.After identifying such solutions, we used the self-consistency equation (22) with the pinning condition (23) to perform their arc-length continuation.terms of the drift speed s and the asymmetry parameter B or the strength of the higher-order interactions K 3 .We do not compute stability along the branches explicitly but summarize stability properties of the branches indicated by direct numerical simulations of finite networks.Note that for small values of K 3 and B, the supercritical Hopf bifurcation of 1-twisted state is mediated by the eigenvalue λ + 1 with the eigenfunction v + 1 (x) = 1 − e −2πix , see Fig. 2(b).For example, Figs. 4 and 5 show the solution branches for K 3 = 0.1 and B = 0.05, respectively.(These values correspond to the bottom horizontal and left vertical lines in Fig. 2.) The branch for K 3 = 0.1, as shown in Fig. 4, starts at the point (B, s) ≈ (0.05, 0.079) in which the drift speed s coincides with the value Im(λ + 1 )/(2π) predicted by Remark 3.1.The branch consists of two parts with different slopes that meet at the fold point (b).Numerical simulations of the finite system (3) suggest that the upper part is stable whereas the lower part is unstable.Moving along the solution branch, we see that as the drift speed s decreases, the original straight profile Θ(x) bends down and up in its left and right sections.Moreover, it seems likely that the solution branch extends asymptotically to B → ∞.In this limit, the speed s vanishes and the phase profile Θ(x) converges to a horizontal line with a phase-slip discontinuity at one point.The shape of the branch for fixed B = 0.05, shown in Fig. 5, is more complex.It is characterized by the non-monotonic dependence of s on K 3 such that at least four different slope parts separated by three fold points can be found in the corresponding diagram.Numerical simulations of the finite system (3) suggest that the negative slope parts of the branch are stable whereas those with positive slope are unstable.Moreover, for large negative values of K 3 we observe similar asymptotic behavior as in the B → ∞ limit in Fig. 4. Now we describe two examples of Hopf bifurcation for larger values of K 3 and B, when this bifurcation is mediated by the eigenvalue λ + 2 with the eigenfunction v + 2 (x) = 1 − e −4πix .Note that in this case the drift speed is equal to Im (λ + 2 )/(4π) as expected for the linear approximation; cf.Remark 3.1.Fig. 6 shows the solution branch for K 3 = 0.5; see also the top horizontal line in Fig. 2. The branch starts at the point (B, s) ≈ (0.27, 0.024), folds at the point (c) and numerical continuation terminates at the point (e).Using numerical simulations of the finite system (3), we find that the lower part of the branch is unstable, whereas the upper part is only stable from the right-most point to some point between (b) and (c).This indicates that before approaching the fold point (c), the solution is destabilized by some dynamical bifurcation, most likely a secondary Hopf bifurcation, although we have not checked this hypothesis rigorously.At (e) the numerical To complete the description of the solution branch shown in Fig. 6, we note that in this case all solutions of the self consistency system (22), (23) satisfy the identities ŵ1 = ŵ4 = ŵ5 = 0. Therefore, we have On the other hand, the profiles a(x) determined by (19) satisfy a symmetry relation a(x + 1/2) = −a(x).With Θ(x) = arg a(x) this means that Θ(x + 1/2) = Θ(x) + π mod 2π.The latter relation is clearly seen in the panels (a)-(e) of Fig. 6.
The last example of this section is the solution branch of ( 22), ( 23) for B = 0.2 as shown in Fig. 7.It shows the nonmonotonic dependence of the drift speed s on the parameter K 3 .Starting from the Hopf bifurcation point, the speed s decreases up to point (c) and then increases for larger values of K 3 .Although there are no fold points on the branch, numerical simulations of the finite system (3) show that the corresponding solutions become unstable for large values of K 3 .This is another sign of the existence of secondary bifurcations the system.Remarkably, the solution branch in Fig. 7 seems to have an asymptotic limit for K 3 → ∞.In this limit, the speed s tends to some nonzero value and the phase profile Θ(x) approaches a two-cluster anti-phase state (cf.Panel (e)) with a phase-slip discontinuity at one point.

Discussion
In this paper, we performed a detailed analysis of Hopf bifurcations of twisted states in a ring network of phase oscillators with nonlocal higherorder interaction.We were able to conduct not only the linear stability analysis of twisted states, but also to describe the global properties of new periodic solutions emanating from the Hopf bifurcation.For the numerical continuation, we focused on the strength of the nonpairwise interactions K 3 -in line with previous work for α 2 = 0 [12]-as well as the asymmetry parameter B of the coupling kernel that facilitates bifurcations to traveling solutions (see [23,24]).While we restricted the bifurcation analysis to compute one-parameter bifurcation diagrams in Section 5, computing these sheds light on potential codimension two bifurcations.Fig. 8 shows one parameter bifurcation diagrams in the asymmetry parameter B for different strength K 3 of the nonpairwise interactions (cf.also Fig. 6).For K 3 ≈ 1.4 it appears that there is a cusp point where the fold point "turns over".Inspecting the phase profile at the fold point for varying parameter K 3 reveals that they are very similar to Fig. 6(c).Note that the phase profile at the end of the numerically computed solution branch in Fig. 6 for small but finite s resembles a −1-twisted (antisplay) phase configuration up to a single twist around the torus as shown in Panel (e); similar phase profiles are obtained at the end of the branches in Fig. 8 (not shown).Indeed, the point in parameter space where numerical continuation terminates is close to a bifurcation of the (stationary) −1-twisted phase configuration at B ≈ 0.09: First, note that the linear stability of the −1-twisted phase configuration is given by ( 12) with B replaced by −B.Now the real eigenvalues λ ± ℓ for ℓ ≥ 4 pass through zero at B = A tan(α 2 ) −1 ≈ 0.09 for the parameter chosen independent of K 3 .Together with the fact that s is close to zero where the numerical continuation terminates, one may speculate that these branches relate to this degenerate bifurcation of the −1-twisted (antisplay) phase configuration.While this may be possible for networks of finitely many oscillators, in the limit of N → ∞ the solutions are topologically different due to distinct winding numbers and the bifurcation likely involves the essential spectrum.Clarifying the nature of the singularity requires further investigation and is beyond the scope of the current work.
In our analysis we focused on specific examples of higher-order interactions and a sinusoidal coupling kernel that facilitated the analysis.Many variations of the model are possible, such as considering other interactions in oscillator rings beyond Kuramoto-type pairwise interactions [25,26] or rings of phase oscillators with nonidentical intrinsic frequencies [27].Furthermore, turbulence is one of the complex dynamical behavior that is observed in the phase oscillator networks [28] and whether this can be understood in terms of the bifurcation scenarios outlined here is an open question.Of course, new phenomena related to twisted states can also be expected for more complicated higher-order interactions, such as adaptive higher-order interactions [29].
The dynamics of phase oscillator networks naturally relate to the dynamics of more general, nonlinear oscillator networks through phase reduction (see, for example [5,6,7]).More specifically, higher-order phase interactions can arise through higher-order corrections to the phase reductions even when the coupling of the nonlinear oscillators is pairwise.While rings of nonlinear oscillators beyond weak coupling can also be analyzed directly [30], we anticipate that the bifurcation mechanisms uncovered here can shed light on the dynamics of nonlinear oscillator networks such as those discussed in [31,32].Making this explicit in a specific model leaves interesting directions for future research.for all other eigenfunctions of A. According to [22,Formula I.9.11]Inserting the ansatz u(x, t) = 1 + v(x, t) into Eq.( 27) and linearizing the result with respect to the small perturbation v, we obtain where η = iΩ + K 2 e −iα 2 ĝ−1 + K 3 e −iα 3 ĝ−2 ĝ1 .
To investigate the decay of different spatial Fourier modes, we insert the ansatz v(x, t) = v + e 2πikx e λt + v − e −2πikx e λt with v + , v − ∈ C and k ∈ Z into Eq.( 28) and equate separately the terms at e λt and e λt .Thus, we obtain a spectral problem

Figure 2 .
Figure 2. Eigenvalues of the linearization of a 1-twisted state.Panel (a) shows max k Re (λ + k ).In the green regions, the 1-twisted state is stable, whereas it is unstable in the red regions.At the boundary a bifurcation occurs.Panel (b) depicts which eigenvalue is the critical one, i.e., for which ℓ we have Re (λ + ℓ ) = max k Re (λ + k ).Finally, Panel (c) shows |Im (λ + ℓ )|, i.e., the modulus of the imaginary part of the critical eigenvalue.Parameter values: K 2 = 1, A = 0.9, α 2 = π 2 − 0.1, α 3 = 0. Dashed lines in all panels indicate the parameter values for the numerical continuation shown in Figs.4-7 below.
Figs. 4-7 show typical solution branches in

Figure 4 .
Figure 4.The drift speed s of the solution of Eq. (20) versus the parameter B for K 3 = 0.1.HB indicates the position of Hopf bifurcation.Panels (a)-(e) show the arguments Θ(x) = arg a(x) of solutions corresponding to the black dots in the main diagram.

Figure 5 .
Figure 5.The drift speed s of the solution of Eq. (20) versus the parameter K 3 for B = 0.05.HB indicates the position of Hopf bifurcation.Panels (a)-(e) show the arguments Θ(x) = arg a(x) of solutions corresponding to the black dots in the main diagram.

Figure 6 .
Figure 6.The drift speed s of the solution of Eq. (20) versus the parameter B for K 3 = 0.5.HB indicates the position of Hopf bifurcation.Panels (a)-(e) show the arguments Θ(x) = arg a(x) of solutions corresponding to the black dots in the main diagram.

Figure 7 .
Figure 7.The drift speed s of the solution of Eq. (20) versus the parameter K 3 for B = 0.2.HB indicates the position of Hopf bifurcation.Panels (a)-(e) show the arguments Θ(x) = arg a(x) of solutions corresponding to the black dots in the main diagram.