Collective dynamics of swarmalators with higher-order interactions

Higher-order interactions shape collective dynamics, but how they affect transitions between different states in swarmalator systems is yet to be determined. To that effect, we here study an analytically tractable swarmalator model that incorporates both pairwise and higher-order interactions, resulting in four distinct collective states: async, phase wave, mixed, and sync states. We show that even a minute fraction of higher-order interactions induces abrupt transitions from the async state to the phase wave and the sync state. We also show that higher-order interactions facilitate an abrupt transition from the phase wave to the sync state by bypassing the intermediate mixed state. Moreover, elevated levels of higher-order interactions can sustain the presence of phase wave and sync state, even when pairwise interactions lean towards repulsion. The insights gained from these findings unveil self-organizing processes that hold the potential to explain sudden transitions between various collective states in numerous real-world systems.


I. INTRODUCTION
The dual interplay between swarming and synchronization is in the heart of swarmalation phenomena which is commonly used to delineate the collective behaviors of entities called swarmalators.With ever-growing advancements and discoveries in multi-agent system studies, it has been observed that there are numerous systems where entities aggregate in space and synchronize over time.Such systems, both natural and man-made, are appropriate examples of swarmalator systems that exhibit simultaneous swarming and synchronization effects.Japanese tree frogs [1], magnetic domain walls [2], swarming robots [3], magnetotactic bacteria [4], vinegar eels [5], Quincke rollers [6], Janus particles [7] are few examples where swarmalation effects are encountered.
Although the two fields synchronization [8,9] and swarming [10][11][12] have been extensively explored over the last few decades, studies on their combined effect began not very long ago.In the seminal work of Vicsek et.al [13], particles moved inside a bounded region, and their directions were influenced by the neighboring particles lying inside a unit radius.Despite being novel to illustrate the phase transition from a desynchronized to a synchronized state, the Vicsek model places little focus on the spatial position and structures of the particles.Later, depending on the spatial movement of the particles, synchronization phenomena were studied by the introduction of mobile agents or moving oscillators [14][15][16].Here also, the effect of spatial position and internal dynamics is unidirectional: the position of the particles influences their internal dynamics, but not the other way round.
The pivotal works that laid the platform for the swarmalator systems were carried out by Tanaka et al. [17] and Isawa et al. [18,19] while studying the movement and dynamics of chemotactic oscillators.The movements of these oscillators are mediated by the surrounding chemical.In 2017, O'Keeffe et al. [20] proposed a simple mathematical model of swarmalators where they move in the two-dimensional plane with Kuramoto-like oscillator dynamics.Five long-term collective states for position aggregation and phase synchronization, viz., static sync, static phase wave, splintered phase wave, active phase wave, and static async were reported in this study.The spatial attraction between two swarmalators was affected by their relative phase, and the spatial distance between them influenced the phase coupling.Adopting this central idea of mutual influence of the spatial position and phase, this model has further been studied with different interaction functions [21][22][23][24], coupling schemes [25][26][27], external forcing [28], large particle limit [29,30], etc. [31,32] A plethora of new collective states has been found which in turn made the researchers interested in delving deeper into the study of such systems.However, mathematical analysis in terms of the solvability of the models or analytical properties of the emerging states was lacking in most of the cases.
To address this problem, O'Keeffe et al. [33,34] came up with a 1D swarmalator model which looked like a pair of coupled Kurmaoto equations, where x i and θ i represent the position on a 1D ring and the phase of the i th swarmalator, respectively for arXiv:2309.03343v1[nlin.AO] 6 Sep 2023 i = 1, 2, . . ., N .v i , ω i are the velocity and internal frequency, and J, K are the inter-element coupling strengths.Being similar to the Kuramoto model, this model was analytically tractable.Analyses of the emerging states had been possible when the model is studied with nonidentical velocities and frequencies [34], distributed couplings [35], random pinning [36,37], thermal noise [38], phase lags [39], etc.All these prior research on swarmalators has predominantly focused on a thorough exploration of their behavior within the framework of pairwise interactions among the constituent entities of the system.More precisely, the spatial and phase dynamics that dictate the interactions among swarmalators have been exclusively governed by the presence of pairwise connections that link them together.Nevertheless, the reliance on a hypothesis rooted solely in pairwise interactions proves inadequate in capturing a wide array of pertinent scenarios [40][41][42].An example includes the group interactions among microorganisms in microbial communities [40].This limitation becomes especially apparent in cases where the interplay of phase and spatial dynamics among individual swarmalators is not solely influenced by pairwise connections among swarmalators but rather by the simultaneous influence of multiple interconnected swarmalators.This intricate interdependence cannot be adequately decomposed solely into the framework of pairwise connections and thus necessitates the introduction of higher-order (group) interactions [40,[43][44][45][46].
Recent advances in physics and other communities have drawn specific attention to the significance of interactions among dynamic units that extend beyond the pairwise realm.Notably, three-and four-way interactions have come to the forefront, revealing their pivotal role in shaping collective behaviors [47][48][49].As a result, the field of network science has shifted its focus toward comprehending higher-order structures to more accurately capture the diverse interactions that exist beyond conventional pairwise connections [43][44][45][46].These intricate interactions are frequently encoded within simplicial complexes [50][51][52], delineating various levels of simplex structures within the network.An assemblage of 1-simplices (edges/links), 2-simplices (filled triangles), and so on constitute the intricate framework of the simplicial complex, reflecting the essence of these higher-order interactions.
In this context, the impact of higher-order interactions on the domain of synchronization has been the subject of thorough investigation in recent years [53][54][55][56][57][58][59][60][61][62].These studies have unveiled that the incorporation of higherorder interactions among dynamic units has the potential to give rise to a plethora of new collective phenomena.However, the influence of higher-order interactions on the realm of swarmalators, which is indissolubly linked to the field of synchronization, remains an unexplored territory to date, and therefore, it is imperative to investigate this uncharted territory.
Motivated by this, in this paper, we propose a model of swarmalators that encompasses both pairwise and higher-order interactions, notably three-body interactions among the swarmalators.These interactions are intricately woven into a simplicial complex framework at the microscopic level.Our proposed model extends the 1D swarmalator model [Eq.(1)] on a ring to the framework that incorporates higher-order interactions among the phase and space dynamics of the swarmalators and thus analytically tractable using the generalized Ott-Antonsen (OA) ansatz [63] in the thermodynamic (N → ∞) limit.Similar to the pairwise model, the present model also displays a diverse range of dynamics, featuring four distinct collective states: async, phase wave, mixed, and sync states.We aim to understand how higher-order interactions influence the formation and characteristics of these distinct collective states.Due to the introduction of higher-order interactions, several significant phenomena arise that are absent when swarmalator interactions are limited to only pairwise connections.Our observations highlight that the inclusion of higher-order interactions leads to abrupt transitions from the async state to the phase wave state and the sync state, contingent on the specific configurations of coupling strengths.We also observe that stronger higherorder interactions can give rise to the persistence of phase wave and sync states, even in cases where pairwise couplings are negative (i.e., repulsive).Furthermore, our findings also reveal that substantial higher-order couplings can facilitate a direct emergence of the synchronized state from the phase wave state, bypassing the intermediate mixed state.

II. MODEL
We consider an ensemble of N swarmalators subjected to two-and three-body interactions, embedded in a simplicial complex at the microscopic level, where the instantaneous position and phase of the i th swarmalator are represented by (x i , θ i ) ∈ (S 1 , S 1 ).When decoupled, each swarmalator is characterized by a set of natural velocity and frequency (v i , ω i ), drawn from a specific distribution g v,ω .For the sake of pedagogy, we here choose the intrinsic frequencies to be drawn from a Lorentzian distribution, g v,ω (x) = ∆v,ω π(x 2 +∆ 2 v,ω ) , with zero mean and half-width ∆ v,ω .The evolution of the swarmalators under the impression of pairwise and triadic interactions is then given by, where (J 1 , K 1 ) and (J 2 , K 2 ) are the pairwise and triadic coupling strengths associated with the spatial and phase interactions, respectively.Note that when (J 2 , K 2 ) = (0, 0), it coincides with the conventional evolution equation of swarmalators over a ring [34].Therefore, analogous to the pairwise swarmalators model over a ring, Eq. (2b) incorporates position-dependent synchronization in swarmalators.To achieve synchronization, we employ the well-known Kuramoto sine terms [48,64] which minimize the pairwise and triadic phase differences between the swarmalators, and the associated distance-dependent cosine terms amplify the coupling intensity among the interacting swarmalators.On the other hand, Eq. (2a) serves as a counterpart of Eq. (2b) and captures phasedependent swarming behavior.In this case, the sine terms minimize the distances between the swarmalators, leading them to aggregate or swarm together, while the cosine terms, similar to Eq. (2b), bolster the coupling between the swarmalators, but this time based on their phase similarity.One can also interpret both Eqs.(2a) and (2b) as models that represent synchronization on the unit torus with both pairwise and three-body interactions.Hence, the model provides a more general framework for swarmalators by considering beyond pairwise interactions.
In order to simplify the model, we convert the trigonometric function to complex exponentials and introduce new variables where i = √ −1 and with and Here Z ± 1 (S ± 1 ) indicates the order parameters associated with the conventional swarmalator model that quantifies the space-phase order of the system [20,33,34].When the correlation between phase and space is perfect (i.e., x i = ±θ i +C 0 , for some constant C 0 ), the value of the order parameter S ± 1 is equal to 1. While, when θ i and x i are uncorrelated, the value of S ± 1 is 0. Therefore the order parameters S ± 1 measure the degree of correlation between the space (x i ) and phase (θ i ) variables, with S ± 1 ranging from 0 (no correlation) to 1 (perfect correlation).On the other hand, Z ± 2 (S ± 2 ) can be interpreted as new order parameters that come up as a result of higher-order interactions, analogous to the higher-order Kuramoto phase models [48,65].In our present study, we focus only on the evolution of conventional order parameters Z ± 1 (S ± 1 ).

III. RESULTS
Now, the coupling dependency of order parameters Z ± 1 (S ± 1 ) in Eqs. ( 3) and ( 5) indicate that depending on the values of coupling strengths, several combinations for order parameters (S + 1 , S − 1 ) can be achieved, which eventually leads to the emergence of different collective states.We, therefore, start by integrating the Eqs.(2a) and (2b) for the calculation of the order parameters S + 1 , S − 1 with N = 10 5 swarmalators and the half-widths of Lorentzian distribution ∆ v = ∆ ω = 1.The results show that depending on the values of coupling strengths, four distinct stable states emerge in the system, which can be classified by the dyad (S + 1 , S − 1 ) as follows: (i) The first state is referred to as the "Async" state, denoted by (S + 1 , S − 1 ) = (0, 0) state.In this state, the swarmalators are uniformly distributed in phase and space, as demonstrated in Fig. 1

(a).
There is no space-phase order among the swarmalators and so (S   2a) and (2b) using Julia Tsit5 adaptive differential equation solver [66] with N = 10 5 swarmalators whose intrinsic velocity and frequency have been drawn from the Lorentzian distribution with zero mean and half-width ∆ω = ∆v = 1.
x i and phase angles θ i are related as x i ≈ ∓θ i , depending on whether it is (S, 0) or (0, S) state, respectively.In the (ξ, η) coordinate system, the swarmalators are partially locked in either ξ i or η i and drift in the other variable.
(iv) The fourth state is known as "Sync" state, denoted by (S, S), where S ̸ = 0 [Fig.1(h)].In this state, the swarmalators are partially locked in both ξ i and η i .For most initial conditions, two clusters of locked swarmalators emerge spontaneously, as depicted in Fig. 1(d).However, one can also observe a single cluster of locked swarmalators for some initial conditions.
Next, in order to analyze all these four states, we employ the Ott-Antonsen (OA) ansatz [63] in the thermodynamic (N → ∞) limit and derive the expressions for order parameters in each state.In the N → ∞ limit, the collective states of the swarmalators can be defined by a continuous function ρ(v, ω, ξ, η, t) as, where ρ(v, ω, ξ, η, t) is the probability to have a swarmalator at time t with intrinsic frequency ω, intrinsic velocity v, and coordinates η and ξ.Differentiating (6) with respect to t, one can obtain the continuity equation Given that our model is a higher-order Kuramoto model occurring on a torus, we are in pursuit of a "torodoidal" OA ansatz [34,63,67], which can be described as a multiplication of Poisson kernels, where α(v, ω, t) and β(v, ω, t) are undetermined and need to be solved in a self-consistent manner, and "c.c" refers to the complex conjugate of its preceding terms.Now plugging the expression for ρ [given by Eq. ( 8)] into the continuity Eq. ( 7), we obtain that the Fourier modes α(v, ω, t) and β(v, ω, t) are subsequently constrained to adhere to the identical conditions for all harmonics p and q, which leads to the fulfillment of a coupled complexvalued differential equation as follows, in the submanifold ∥α∥ = 1 = ∥β∥.Subsequently, the order parameters Equations ( 9a)-(10b) contain a set of self-consisting equations for the order parameters Z ± m in the N → ∞ limit.

A. Stability of the async state
In the async state, the order parameters Z ± m are zero.When Z ± m = 0, Eqs.(9a) and (9b) give the solu- Clearly, these solutions are selfconsistent, as substituting them into the Eqs.(10a) and (10b) one can obtain , which converges to zero in the t → ∞ limit.Now, to investigate the stability of the async state, we introduce a small perturbation around the solutions α 0 and β 0 given by This eventually gives the perturbed order parameters as where we assume that 1 ≫ ∥α 1 (v, ω, t)∥, ∥β 1 (v, ω, t)∥ and Z ± 11 .Substituting the expressions for α 1 , β 1 and Z ± 11 into the Eqs.(9a) and (9b), and considering the terms up to first order, we obtain the following set of evolution equations One can notice that the terms with α 0 and β 0 as a function of v and ω oscillate rapidly for t ≫ 1.These rapidly oscillating terms barely contribute when integrating over v and ω, and thus to approximate the integrals, we can neglect these small terms.Now, introducing a new variable τ = v + ω, α and β can be expressed as α(v, ω, t) = α(y, t), and β(v, ω, t) = β(y, t), respectively.Then, integrating the Eqs.(13a) and (13b) over v and ω, we eventually obtain where To evaluate the integration, we use the fact that it has a residue τ = i(∆ v + ∆ ω ) in the upper half plane where α * 0 (τ ) and β * 0 (τ ) are analytic.Now, the async state becomes stable when the perturbed order parameters Z ± 11 die out in time.Hence, from Eqs. (14a) and (14b), one can conclude that the async state becomes unstable when J + 1 − 2(∆ v + ∆ ω ) > 0, or in other words the async state sustains its stability for J + 1 −2(∆ v +∆ ω ) < 0. Therefore, the curve satisfying signifies the critical curve above which the async state loses its stability, and the swarmalators form a phase wave state.Eq. ( 16) reveals that the transition from async state to phase wave state depends solely on the pairwise interactions by means of the pairwise coupling strengths (J 1 , K 1 ).

B. Analysis of phase wave state
In the phase wave state, swarmalators develop a phase wave or band with x i = ∓θ i for (S, 0) and (0, S) states, respectively.Therefore, we seek a solution to Eqs. (9a)-(10b) that satisfies α = 0, β ̸ = 0, Z + m ̸ = 0, and Z − m = 0, (m = 1, 2).Substituting these relations into the Eqs.(9a) and (9b), we have Notice that α(v, ω, t) depends on v + ω, which is distributed according to a Lorentzian distribution with spread ∆ v + ∆ ω .This allows us to evaluate the integral for the order parameter Z + m explicitly as Z + 1 = α * (i∆ v + i∆ ω , t) using the Cauchy's residue theorem by closing the contour to an infinite-radius semicircle in the upper half-plane.Similarly, we can obtain Substituting the relation between Z + 1 and Z + 2 into the Eqs.(17a) and (17b), and assuming ψ + 1 = 0, we obtain the expressions for α and β as follows, and where we introduce a function F as and the term r is given by, Equation ( 19) results in Z − m = 0, as anticipated.On the other hand, Eq. ( 18) implies that Solving Eq. ( 22) for S + 1 gives the expression of S + 1 as where ∆ = ∆ v + ∆ ω .The plus and minus sign inside the square root corresponds to the stable and unstable solutions if they exist.Equation ( 23) suggests the presence of a bistable region upon the selection of coupling strengths.In this scenario, a particular combination of coupling values leads to the coexistence of asynchronous and phase wave states depending on the initial conditions chosen.Saying differently, with varying coupling strengths, there exist two different transition scenarios: one corresponds to the transition from the async state to the phase wave state, referred to as forward transition, and the other is backward transition, which is associated with the transition from phase wave state to async state.
For the forward transition case, S + 1 bifurcates from 0 (i.e., async state) at which is consistent with Eq. ( 16), where forward transition from async state to phase wave state emerges and only the stable branch of S + 1 exists.This once again guarantees that the transition from async state to phase wave state is solely dependent on the pairwise coupling strengths.On the other hand, in the case of backward transition, both the stable and unstable branches coexist for a region of parameter values called the hysteresis region.But as soon as the backward critical coupling condition J + 1 = J + 1,b is reached, both the stable and unstable branches clash and destroy each other.As a result, the stability of the phase wave is completely shattered, and S + 1 = 0 is the sole viable solution.From Eq. ( 23), we obtain that both the stable and unstable branches of S + 1 exist if the following coupling condition satisfies, Eq. ( 25) thus corresponds to the critical curve for the backward transition.
To illustrate this, we numerically integrate the Eqs.(2a) and (2b) with N = 10 5 swarmalators for the calculation of S ± 1 .In Fig. 2(a) we plot the variation of order parameter S ± 1 as a function of pairwise coupling strength K 1 .K 1 is first increased adiabatically from K 1 = 0 to an adequately large value and then decreased back with the other couplings fixed at J 1 = 1, J 2 = 5, and K 2 = 9.The results reveal that S − 1 remains zero all the time, while an abrupt transition from S + 1 ≈ 0 to S + 1 ≈ 0.7 occurs at K 1 = 7 [obtained from Eq. ( 24)] as the coupling strength K 1 is increased.Another abrupt transition from S + 1 ≈ 1 to S + 1 ≈ 0 occurs at K 1 = 6 [obtained from Eq. ( 25)] as K 1 is decreased starting from the phase wave state.Therefore, the system supports bistability behavior (where both phase wave and async states are stable) for K 1 ∈ [6,7].This is further substantiated by our theoretical predictions regarding the order parameters (illustrated using continuous and dashed magenta lines for the stable and unstable branches, respectively), demonstrating a commendable concurrence with the numerical findings (represented by solid circles).To emphasize the bistability nature, in Fig. 3 we demonstrate the async (upper row) and phase wave (lower row) states, respectively, for K 1 = 6.5 while keeping the other couplings fixed at specific values as earlier.
The outcomes given above highlight a significant finding that the introduction of higher-order interactions can cause a sudden transition from the async state to phase wave state and vice-versa, which was not the case with only pairwise interactions among the swarmalators [34].Thereafter, to better understand the effect of higher-order interactions in promoting bistability behavior, we plot the complete stability profile for the system in Fig. 2(b).It shows that for adequately small values of higher-order coupling (J + 2 < 4), the transition from async state (I) to phase wave state (III) is continuous and takes place through a supercritical pitchfork bifurcation at J + 1 = 4.However, for larger values of higherorder coupling (J + 2 > 4), the pitchfork bifurcation at J + 1 = 4 becomes subcritical and a saddle-node bifurcation emerges at a lower value of J + 1 , given by Eq. ( 25) (depicted in blue).These two bifurcations correspond to the sudden transition observed in Fig. 2(a), and the region bounded by them refers to the region of bistability (II) between async and phase wave state.Figure 2(b) also reveals an interesting observation that for relatively larger values of higher-order coupling (J + 2 ≥ 16), the region of bistability stretches into the negative region J + 1 < 0, demonstrating the fact that higher-order interactions can stabilize the phase wave state even when the pairwise interactions are repulsive.

C. Analysis of sync state
In the sync state we have (S + 1 , S − 1 ) = (S, S), with S ̸ = 0. Therefore, we seek solutions to the Eqs.(9a)-(10b) 1 as a function of pairwise coupling strength K1 for J1 = 1 and fixed three-body coupling strengths J2 = 5, K2 = 9.Solid and dashed magenta curves represent the stable and unstable solutions given by Eq.( 23), respectively.The dashed vertical lines correspond to the critical couplings for forward (in red) and backward (in black) transitions obtained from Eqs. ( 24) and ( 25 and higher-order coupling J + 2 = J 2 +K 2

2
. Two distinct types of bifurcations, saddle-node and pitchfork, are represented by blue and red curves, respectively.These two curves intersect each other at (J + 1 , J + 2 ) = (4, 4).When J + 2 < 4, the pitchfork bifurcation is deemed supercritical, whereas when J + 2 > 4, the pitchfork bifurcation becomes subcritical.3: Bistablilty between async and phase wave state.Scatter plot for async (0, 0) and phase wave (S, 0) states at K1 = 6.5, J1 = 1, J2 = 5, and K2 = 9 [drawn from the region (II) in Fig. 2(a)] is depicted in the upper and lower row, respectively.The order parameter S + 1 (S − 1 ) is indicated by the black (magenta) circle in the left column, where the values of S + 1 and S − 1 are represented by the length of the line joining the center with the respective circles.Clearly, in the upper row, both the values of order parameters are almost zero, indicating the async (0, 0) state, while in the lower row, the value of S + 1 is non-zero and S − 1 is zero, characterizing the phase wave (S, 0) state.In the right panel, the corresponding scatter plots in the (x, θ) plane are displayed.Swarmalators are uniformly distributed in both phase and space, characterizing the async state (in the upper row), whereas in the lower row, swarmalators display a correlation between phase and space and thus correspond to the phase wave state.
that beyond K 2 = 4, a bistable dynamics emerges in the system due to the interplay between saddle-node and subcritical pitchfork bifurcations, and the associated region of bistability (II) is bounded by the curves of these bifurcations.On the other hand, for K 2 < 4, a continuous transition from async state (I) to sync state (III) occurs at K 1 = 4 through a supercritical pitchfork bifurcation.The stretch of the bistability region (II) in the regime K 1 < 0 reveals the fact that for sufficiently larger higher-order coupling strengths, the system can achieve a stable sync state even when the pairwise interactions are repulsive (i.e., K 1 = J 1 < 0).

General case: J1
In contrast to the earlier scenario, in this case, the coupled complex-valued differential equations (9a) and (9b) cannot be separated from each other.As a result, it's not possible to directly express the order parameters Z ± 1 in terms of α and β, respectively.Therefore, we look for a solution of Eqs.(9a) and (9b) that satisfy α ̸ = 0, β ̸ = 0, Z ± m ̸ = 0. Assuming ψ ± m = 0, we find where (31) Now, we solve the integrals for the order parameters given by Eqs.(10a) and (10b) using residue theorem.Further, considering the fact that in the sync state S + 1 = S − 1 = S, we obtain where ). Eq. ( 32) eventually provides the expression for the order parameter, which can be acquired by solving the following equation for S where y = S 2 .From Eq. ( 33), it is clear that S bifurcates from zero at ∆ 1 = D, i.e., Notice that for J 1 = K 1 , the condition (34) coincides with the critical coupling condition (29) and thus once In (a), the behavior of the order parameter S ± 1 is presented in relation to the pairwise coupling strength K1.With a fixed setting of J1 = K1 and constant three-body coupling strengths J2 = K2 = 9, both stable and unstable solutions (obtained from Eq. ( 28)) are showcased through solid and dashed magenta curves, respectively.The dashed vertical lines correspond to critical coupling values for forward (red) and backward (black) transitions.These values are derived from Eqs. ( 29) and (25), with J1 = K1, J2 = K2.The outcome obtained from direct simulations of Eqs.(2a) and (2b) is depicted in solid circles for N = 10 5 swarmalators with Lorentzian distribution half-widths ∆ω = ∆v = 1.The findings distinctly reveal an abrupt transition from the async (I) state to the sync (III) state, establishing a bistable domain (II) where both states are stable.Moving to (b), a comprehensive stability diagram is presented, illustrating the states of async (I), phase wave (III), and bistability (II) as functions of pairwise coupling K1 and higher-order coupling K2.The diagram encompasses two distinct types of bifurcations: the blue curve representing a saddle-node bifurcation and the red curve depicting a pitchfork bifurcation.Notably, these curves intersect at (K1, K2) = (4, 4).When K2 < 4, the pitchfork bifurcation takes a supercritical nature, whereas, for K2 > 4, the pitchfork bifurcation becomes subcritical.Scatter plot for async (0, 0) and sync (S, S) states at K1 = J1 = 3.5, and K2 = J2 = 9 [drawn from the region (II) in Fig. 4(a)] is depicted in the upper and lower row, respectively.The order parameter S + 1 (S − 1 ) is indicated by the black (magenta) circle in the left column, where the values of S + 1 and S − 1 are represented by the length of the line joining the center with the respective circles.Clearly, in the upper row, the values of both order parameters are almost zero, indicating the async (0, 0) state, while in the lower row, both the order parameters take non-zero equal value, characterizing the sync (S, S) state.In the right panel, the corresponding scatter plots in the (x, θ) plane are displayed.Swarmalators are uniformly distributed in both phase and space, characterizing the async state (in the upper row), whereas in the lower row, swarmalators are locked in both phase and space and thus correspond to the sync state.
again guarantees that the sync state can bifurcate directly from the async state without passing through any intermediate state [Fig.4].However, in general circumstances (i.e., J 1 ̸ = K 1 , J 2 ̸ = K 2 ), the transition from async state to sync state is not direct but occurs by passing through intermediate states, namely phase wave state and mixed state.To describe this, in Fig. 6 we plot the order parameters S ± 1 for varying K 1 while keeping the other couplings fixed at nominal values J 1 = 9, J 2 = 6.5, and K 2 = 5.5.The solid magenta and blue lines correspond to the analytical predictions of stable solutions of phase and sync states obtained from Eqs. ( 23) and (33), respectively, which are in favorable alignment with the results obtained through direct simulation (shown in solid and void markers) for the finite-size system.As observed, with increasing K 1 , the system passes from the async state (I) to the sync state (V ) through the intermediate phase wave state (III) and mixed state (IV ).Here, the evolution of the order parameters (S + 1 , S − 1 ) shows that with increasing K 1 , the system passes from (0, 0) (i.e., async) state to (S, 0) (i.e., phase wave) state through an abrupt transition, resulting into an interval of bistability (II) where both the async and phase wave states are stable.On the other hand, the system goes through a continuous transition from (S, 0) (i.e., phase wave) state to (S, S) (i.e., sync) state, resulting in the intermediate mixed state (S ′ , S ′′ ) where S ′ > S ′′ .

D. Mixed State
In the mixed state we have (S + 1 , S − 1 ) = (S ′ , S ′′ ), with S ′ ̸ = S ′′ .This state resides as an intermediary between the phase wave and sync states and bifurcates from the (S, 0), or (0, S) state when S ′ > S ′′ or S ′′ > S ′ , through a continuous transition to the (S, S) state.In Fig. 6, the shaded black region (IV ) depicts the interval of K 1 (keeping the other coupling fixed at a nominal value) for which the system passes through the mixed state.Unfortunately, we are unable to find the analytical boundaries in terms of the coupling strengths for the mixed state, and thus obtain the domain of mixed state through the numerical investigations by evaluating the order parameters S ± 1 with S + 1 > S − 1 .The distinctive characteristic of the mixed state lies in the fact that while both S ′ and S ′′ remain time-independent, the functions α and β exhibit time dependency [34].This contrasts with the timeindependent equations ( 18) and ( 30) associated with the phase wave and sync states, respectively.

E. Abrupt transition from phase wave state to sync state
As discussed above, in general, with increasing coupling strengths, a continuous transition from the phase wave state to the sync state emerges in the system through an intermediate mixed state.This smooth transition phenomenon is consistent with previously acquired results for only pairwise interactions between the swarmalators [34].When the higher-order coupling strengths are relatively small, a comparable transition phenomenon is still evident in our present system with higher-order interactions [Fig.6].Conversely, when sufficiently significant higher-order couplings come into play, a critical observation comes to light.In this scenario, the influence of higher-order interactions leads to an abrupt and noteworthy transition from the (S, 0)/(0, S) state to the (S, S) state, bypassing the intermediate mixed state.To illustrate this, in Fig. 7 we plot the order parameters (S + 1 , S − 1 ) as a function of K 1 for adequately large higherorder couplings J 2 = 8 and K 2 = 9, keeping J 1 fixed at J 1 = 7.The solid magenta and blue lines represent the analytically predicted stable solutions for the phase and synchronized states, respectively, derived from Eqs. ( 23) and (33).Impressively, these analytical predictions closely match the outcomes obtained through direct simulations on the finite size system, depicted by the solid 1 as a function of pairwise coupling strength K1 for J1 = 9 and fixed three-body coupling strengths J2 = 6.5, K2 = 5.5.Solid magenta and blue curves represent the stable solutions for phase wave and sync states, given by Eqs.( 23) and ( 33), respectively.Solid circles depict the result obtained from the direct simulation of Eqs.(2a) and (2b) for N = 10 5 oscillators with half widths of the Lorentzian distribution ∆v = ∆ω = 1.An abrupt transition from the async (I) to phase wave (III) state is observed, which creates a domain of bistability (II), where both the states are stable.On the other hand, from the phase wave state (III) to sync state (V ), a smooth (continuous) transition occurs, resulting in the shaded region (IV ), which corresponds to (S ′ , S ′′ ) state, with S ′ > S ′′ , and called as mixed state.and void markers.As K 1 is being first adiabatically increased to a large value and then decreased, we observe two different abrupt transitions, resulting in two distinct bistable domains.The first one corresponds to a sudden transition from the (0, 0) state (I) to the (S, 0) state (III), inducing a region of bistability (II), where both the async and phase wave states are stable (shaded grey region).This bistable phenomenon has already been discussed in Sec.III B. The second one is associated with an abrupt transition from (S, 0) state (III) to (S, S) state (V ), giving rise to a bistable domain (IV ), where both stable phase wave and sync state can emerge (shaded red region).To elucidate this bistability nature, in Fig. 8, we demonstrate the two distinct states: phase wave state (upper row) and sync state (lower row) for fixed K 1 = 1.5, drawn from the region of bistability (IV ).

IV. DISCUSSIONS
Summing up, here we have introduced an analytically tractable model of swarmalators that incorporates both pairwise and higher-order interactions (specifically threebody interactions), embedded in a simplicial complex at the microscopic level.This proposed model exhibits a high degree of complexity with four distinct collective states, namely async, phase wave, mixed, and sync states.The higher-order interactions introduce supplementary layers of nonlinearity into the behavior of the microscopic system dynamics.As a result, a few pivotal phenomena emerge that remain absent when interactions between the swarmalators are confined to only pairwise connections and lack the influence of higher-order interactions.We observe that the inclusion of higher-order interactions leads to abrupt transitions from the async state to either the phase wave state or the sync state, depending on the specific coupling strength configurations.Our observations also reveal that when the higher-order interactions are sufficiently strong, phase wave and sync state can emerge and persist, even in scenarios where pairwise couplings are repulsive.This implies that despite the potential decay of specific coupling types, the existence of alternative forms of coupling can play a pivotal role in upholding the regimes of bistability between async and phase wave (sync) states.Furthermore, our findings extend to the discovery that substantial higherorder couplings can facilitate a direct emergence of the sync state from the phase wave state without passing through the mixed state.This distinct behavior stands in contrast to situations involving solely pairwise interactions, where the synchronized state bifurcates from the phase wave state through the intermediate mixed state.
Thus, our study makes a substantial contribution to understanding the influence of higher-order interactions on shaping the collective dynamics of swarmalators, although numerous avenues for further exploration remain open.In this context, we specifically focus on higherorder interactions up to the third order, employing an all-to-all coupling arrangement.Consequently, investigating interactions beyond the three-body and exploring localized coupling configurations becomes a particularly 1 as a function of pairwise coupling strength K1 for J1 = 7 and fixed three-body coupling strengths J2 = 8, K2 = 9.Solid magenta and blue curves represent the stable solutions for phase wave and sync states, given by Eqs.( 23) and (33), respectively.Solid circles depict the result obtained from the direct simulation of Eqs.(2a) and (2b) for N = 10 5 oscillators with half widths of the Lorentzian distribution ∆v = ∆ω = 1.Two distinct abrupt transition is revealed.First, an abrupt transition from the async (I) to phase wave (III) state is observed, which creates a domain of bistability (II), where both the async and phase wave states are stable.Second, from the phase wave state (III) to sync state (V ), generating a bistable domain (IV ), where both the phase wave and sync states are stable.
intriguing prospect for future research.It is also worth noting that our current investigation is limited to analyzing the effect of higher-order interactions solely within a 1D swarmalator model on a ring.As such, an equally captivating avenue for future exploration lies in examining how higher-order interactions impact the collective behavior within 2D and other higher-dimensional swarmalator systems. 1 is non-zero, and S − 1 is zero, representing the phase wave (S, 0) state, while in the lower row, both the order parameters take non-zero equal value, characterizing the sync (S, S) state.In the right panel, the corresponding scatter plots in the (x, θ) plane are displayed.Swarmalators display a correlation between phase and space in the upper row, characterizing the async state, whereas, in the lower row, swarmalators are locked in both phase and space and thus correspond to the sync state.

Figure 2 :
Figure 2: Abrupt transition from the async state to the phase wave state.(a) Order parameter S ±1 as a function of pairwise coupling strength K1 for J1 = 1 and fixed three-body coupling strengths J2 = 5, K2 = 9.Solid and dashed magenta curves represent the stable and unstable solutions given by Eq.(23), respectively.The dashed vertical lines correspond to the critical couplings for forward (in red) and backward (in black) transitions obtained from Eqs. (24) and (25), respectively.Solid circles depict the result obtained from the direct simulation of Eqs.(2a) and (2b) for N = 10 5 oscillators with half widths of the Lorentzian distribution ∆v = ∆ω = 1.The results reveal an abrupt transition from async (I) state to phase wave (III) state, which results in a bistable domain (II) where both the referred states are stable.(b) The comprehensive stability diagram illustrating the states of async (I), phase wave (III) and bistability (II) as a function of pairwise coupling J + 1 = J 1 +K 1 Figure 2: Abrupt transition from the async state to the phase wave state.(a) Order parameter S ±1 as a function of pairwise coupling strength K1 for J1 = 1 and fixed three-body coupling strengths J2 = 5, K2 = 9.Solid and dashed magenta curves represent the stable and unstable solutions given by Eq.(23), respectively.The dashed vertical lines correspond to the critical couplings for forward (in red) and backward (in black) transitions obtained from Eqs. (24) and (25), respectively.Solid circles depict the result obtained from the direct simulation of Eqs.(2a) and (2b) for N = 10 5 oscillators with half widths of the Lorentzian distribution ∆v = ∆ω = 1.The results reveal an abrupt transition from async (I) state to phase wave (III) state, which results in a bistable domain (II) where both the referred states are stable.(b) The comprehensive stability diagram illustrating the states of async (I), phase wave (III) and bistability (II) as a function of pairwise coupling J + 1 = J 1 +K 1 2

Figure
Figure3: Bistablilty between async and phase wave state.Scatter plot for async (0, 0) and phase wave (S, 0) states at K1 = 6.5, J1 = 1, J2 = 5, and K2 = 9 [drawn from the region (II) in Fig.2(a)] is depicted in the upper and lower row, respectively.The order parameter S + 1 (S − 1 ) is indicated by the black (magenta) circle in the left column, where the values of S + 1 and S − 1 are represented by the length of the line joining the center with the respective circles.Clearly, in the upper row, both the values of order parameters are almost zero, indicating the async (0, 0) state, while in the lower row, the value of S + 1 is non-zero and S − 1 is zero, characterizing the phase wave (S, 0) state.In the right panel, the corresponding scatter plots in the (x, θ) plane are displayed.Swarmalators are uniformly distributed in both phase and space, characterizing the async state (in the upper row), whereas in the lower row, swarmalators display a correlation between phase and space and thus correspond to the phase wave state.

Figure 4 :
Figure 4: Abrupt transition from async to sync state.In (a), the behavior of the order parameter S ± 1 is presented in relation to the pairwise coupling strength K1.With a fixed setting of J1 = K1 and constant three-body coupling strengths J2 = K2 = 9, both stable and unstable solutions (obtained from Eq. (28)) are showcased through solid and dashed magenta curves, respectively.The dashed vertical lines correspond to critical coupling values for forward (red) and backward (black) transitions.These values are derived from Eqs. (29) and (25), with J1 = K1, J2 = K2.The outcome obtained from direct simulations of Eqs.(2a) and (2b) is depicted in solid circles for N = 10 5 swarmalators with Lorentzian distribution half-widths ∆ω = ∆v = 1.The findings distinctly reveal an abrupt transition from the async (I) state to the sync (III) state, establishing a bistable domain (II) where both states are stable.Moving to (b), a comprehensive stability diagram is presented, illustrating the states of async (I), phase wave (III), and bistability (II) as functions of pairwise coupling K1 and higher-order coupling K2.The diagram encompasses two distinct types of bifurcations: the blue curve representing a saddle-node bifurcation and the red curve depicting a pitchfork bifurcation.Notably, these curves intersect at (K1, K2) = (4, 4).When K2 < 4, the pitchfork bifurcation takes a supercritical nature, whereas, for K2 > 4, the pitchfork bifurcation becomes subcritical.

Figure 5 :
Figure 5: Bistablilty between async and sync state.Scatter plot for async (0, 0) and sync (S, S) states at K1 = J1 = 3.5, and K2 = J2 = 9 [drawn from the region (II) in Fig.4(a)] is depicted in the upper and lower row, respectively.The order parameter S + 1 (S − 1 ) is indicated by the black (magenta) circle in the left column, where the values of S + 1 and S − 1 are represented by the length of the line joining the center with the respective circles.Clearly, in the upper row, the values of both order parameters are almost zero, indicating the async (0, 0) state, while in the lower row, both the order parameters take non-zero equal value, characterizing the sync (S, S) state.In the right panel, the corresponding scatter plots in the (x, θ) plane are displayed.Swarmalators are uniformly distributed in both phase and space, characterizing the async state (in the upper row), whereas in the lower row, swarmalators are locked in both phase and space and thus correspond to the sync state.

Figure 6 :
Figure 6: Transition from async to sync state through intermediate mixed state.Order parameter S ±1 as a function of pairwise coupling strength K1 for J1 = 9 and fixed three-body coupling strengths J2 = 6.5, K2 = 5.5.Solid magenta and blue curves represent the stable solutions for phase wave and sync states, given by Eqs.(23) and(33), respectively.Solid circles depict the result obtained from the direct simulation of Eqs.(2a) and (2b) for N = 10 5 oscillators with half widths of the Lorentzian distribution ∆v = ∆ω = 1.An abrupt transition from the async (I) to phase wave (III) state is observed, which creates a domain of bistability (II), where both the states are stable.On the other hand, from the phase wave state (III) to sync state (V ), a smooth (continuous) transition occurs, resulting in the shaded region (IV ), which corresponds to (S ′ , S ′′ ) state, with S ′ > S ′′ , and called as mixed state.

Figure 7 :
Figure 7: Abrupt transitions from async to phase wave state, and phase wave to sync state.. Order parameter S ±1 as a function of pairwise coupling strength K1 for J1 = 7 and fixed three-body coupling strengths J2 = 8, K2 = 9.Solid magenta and blue curves represent the stable solutions for phase wave and sync states, given by Eqs.(23) and (33), respectively.Solid circles depict the result obtained from the direct simulation of Eqs.(2a) and (2b) for N = 10 5 oscillators with half widths of the Lorentzian distribution ∆v = ∆ω = 1.Two distinct abrupt transition is revealed.First, an abrupt transition from the async (I) to phase wave (III) state is observed, which creates a domain of bistability (II), where both the async and phase wave states are stable.Second, from the phase wave state (III) to sync state (V ), generating a bistable domain (IV ), where both the phase wave and sync states are stable.

Figure 8 :
Figure 8: Bistablilty between phase wave and sync state.Scatter plot for phase wave (S, 0) and sync (S, S) states at K1 = 1.5, J1 = 3.5, K2 = 9, and J2 = 8 [drawn from the region (IV ) in Fig. 7] is depicted in the upper and lower row, respectively.The order parameter S + 1 (S − 1 ) is indicated by the black (magenta) circle in the left column, where the values of S + 1 and S − 1 are represented by the length of the line joining the center with the respective circles.Clearly, in the upper row, the value of S +1 is non-zero, and S − 1 is zero, representing the phase wave (S, 0) state, while in the lower row, both the order parameters take non-zero equal value, characterizing the sync (S, S) state.In the right panel, the corresponding scatter plots in the (x, θ) plane are displayed.Swarmalators display a correlation between phase and space in the upper row, characterizing the async state, whereas, in the lower row, swarmalators are locked in both phase and space and thus correspond to the sync state.