{'fraction_non_alphanumeric': 0.03684731920651412, 'fraction_numerical': 0.006862058134196728, 'mean_word_length': 4.565147236789028, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-body Interactions and Non-Linear Consensus Dynamics on Networked Systems

Multi-body interactions can reveal higher-order dynamical effects that are not captured by traditional two-body network models. In this work, we derive and analyse models for consensus dynamics on hypergraphs, where nodes interact in groups rather than in pairs. Our work reveals that multi-body dynamical effects that go beyond rescaled pairwise interactions can only appear if the interaction function is non-linear, regardless of the underlying multi-body structure. As a practical application, we introduce a specific non-linear function to model three-body consensus, which incorporates reinforcing group effects such as peer pressure. Unlike consensus processes on networks, we find that the resulting dynamics can cause shifts away from the average system state. The nature of these shifts depends on a complex interplay between the distribution of the initial states, the underlying structure and the form of the interaction function. By considering modular hypergraphs, we discover state-dependent, asymmetric dynamics between polarised clusters where multi-body interactions make one cluster dominate the other.


I. INTRODUCTION
In recent years, researchers have developed a variety of higher-order models for complex systems in order to enrich the standard network formalism when it is not sufficient to capture their structure and function [1]. Among those, multi-body models, sometimes called combinatorial models, focus on the importance of group interactions, that is situations when the basic unit of interaction involves more than two nodes. Multi-body interactions are observed in different areas in nature [2], society [3] and technology [4], with examples ranging from collaborations of authors [5] to neuronal activity [6,7]. Such systems may be represented as hypergraphs or simplicial complexes, and a substantial body of work has characterised their structural properties. However, a proper understanding of how multi-body interactions affect spreading dynamics in networked systems is still nascent [3,[8][9][10][11].
When considering a dynamical process driven by multibody interactions, it is important to distinguish between an accumulation of two-body dynamics and genuine multi-body dynamics. In the former case, the influence on an agent can be fully explained by its pairwise relationships to other group members. In contrast, multi-body dynamical systems additionally account for the effect of the group as a whole. An area where this distinction is often blurred is the modelling of complex diffusion in social networks. The adoption of norms and opinion between agents is known to be a non-trivial process, which may drastically differ from a simple model of epidemic spreading. Experiments in social psychology such as the conformity experiment [12] indicate that multiple exposures might be necessary for an agent to adopt a certain state. This phenomenon may be described with threshold models, in which the state of agents switches if a certain fraction of their neighbours agrees. Observe that these models are based on independent, two-body interactions that are linearly accumulated and therefore do not account for multi-body effects resulting from those edges being part of a closed group [13]. Yet, it is well known in sociology that the dynamics in a social clique are determined not just by the pairwise relationships of its members, but often by complex mechanisms of peer influence and reinforcement [14][15][16].
There are different ways of encoding the multi-body structure of a networked system. They range from set systems and general hypergraphs [10,17] to approaches adapted from algebraic topology based on simplicial complexes [9,14,[18][19][20][21][22]. Hypergraphs generalise standard networks by allowing nodes to interact through hyperedges of different sizes. Size two corresponds to two-body edges, size three to three-body edges and, in general, size k to k-body edges, associated to k-cliques. Importantly, different types of triangles may be found in hypergraphs, those obtained by combining three two-body edges, and those corresponding to one three-body edge, and they have a different nature. In simplicial complex approaches, a fully connected group of d + 1 nodes is formalised as a d-simplex. The unity of all simplices that are formed by the underlying edge structure of a graph is called a simplicial complex and can be used as the underlying topology [9,14]. These complexes can be broken down into different simplex skeletons that depend on the dimension of simplices they include. The 1-skeleton of a simplicial complex just includes (two-body) edges and is equivalent to the conventional network structure. The 2-skeleton includes all simplices up to dimension two (nodes, edges and triangles) and a general d-skeleton all simplices up to dimension d [9,19].
Hypergraphs and simplicial complexes are natural languages to describe the structure of multi-body networked systems. However, understanding how these substrates affect the dynamics requires an additional modelling step. Multi-body dynamical models that use simplicial complexes may, for instance, exploit an algebraic structure called the Hodge-Laplacian, analogous to the graph Laplacian [9]. However, the resulting diffusive process becomes more complicated to study, as diffusion is now defined between simplices of any dimension, and not only between nodes as is usually the case on networks. This generalisation leads to challenges in both analysis and interpretability, as state variables can usually be measured only on nodes, and not on other network entities, in empirical data. In this paper, we explore a different modelling approach based on hypergraphs where, more realistically, the states are defined on the nodes, and their evolution depends on the combined values of all the nodes inside each multi-body interacting unit. For the sake of simplicity, we focus mostly on three-body dynamical systems as a first extension to two-body dynamical systems [23], that is we consider the dynamics taking on an hypergraph where all hyperedges have size three. We show that linear dynamics on hypergraphs can always be rewritten as a dynamics on a standard, two-body network. This observation emphasises that non-linear interactions are necessary for the multi-body dynamical system not to be reducible to a two-body dynamical system. As a next step, we propose and analyse a minimal non-linear model for consensus dynamics including reinforcing group effects, which we study analytically and by means of numerical simulations.
The remainder of this article is outlined as follows. In Section II we describe the model for three-body dynamical systems. We show further in Section III the role of the interaction function and differentiate between the linear and non-linear case. Section IV is devoted to the analysis of a specific non-linear interaction function which captures reinforcement behaviour in triangles. Finally in Section V we discuss the consequences of incorporating multi-body interactions and propose future research directions.

II. FROM TWO-BODY TO THREE-BODY INTERACTIONS
Two-body dynamical systems describe general, possibly non-linear, dynamics on conventional edge-based networks [23]. Let G be a network consisting of a set V (G) = {1, . . . , N } of N nodes connected by a set of edges E(G) = {(i, j) : i, j ∈ V (G)}, described by ordered tuples of nodes. The structure of the network can be represented by the adjacency matrix A ∈ R N ×N with entries For undirected networks, which we consider in this paper, A is a symmetric matrix. Each node i ∈ V (G) is endowed with a dynamical variable, x i ∈ R for i ∈ {1, . . . , N }, whose evolution is determined by the underlying graph G, defining the edges over which interactions take place, and the set of interaction functions quantifying how the states of neighbouring nodes affect each other. The time-evolution of the system is then defined by where the total interaction on a node is the sum, over its neighbours, of possibly non-linear interactions. Particular examples of (3) include the Kuramoto model [24], continuous-time random walks and linear consensus on networks [25]. For three-body dynamical systems, the topology is encoded by a hypergraph H consisting of a set V (G) = {1, . . . , N } of N nodes connected by a set of three-body interactions, or triangles. We describe undirected triangles T (G) = {{i, j, k} : i, j, k ∈ V (G)} as unordered triples of nodes. The structure of the network can be represented by the adjacency tensor A ∈ R N ×N ×N with entries The adjacency tensor is symmetric in all dimensions. Moreover, we assume that the graph has no self-loops, to avoid classifying edge-combinations involving self-loops as triangles. Thus, A ijk = 0 if i = j, j = k, or i = k. By this definition, we consider a system that consists only of three-body interactions, since A cannot describe twobody edges. A three-body dynamical system generalising Eq. (3) is defined by its structure, encoded by H, and by the set of interaction functions quantifying, for each triangle {i, j, k}, the combined influence of the states of nodes j and k on the state i. The dynamics of the variables x i ∈ R on the hypergraph is then given by In an undirected setting, it is natural for the symmetries of f (jk) i to be aligned with the symmetries of A. Therefore, we assume that f for each nodes. The differences between two-body or three-body interactions are visualised in Figure 1.

A. Symmetries and quasi-linearity
Our main goal is to examine if, and under which conditions, a three-body dynamical system (6) can not be rewritten as a two-body dynamical system (2), that is when a multi-body formalism is truly necessary to capture the complexity of a system. We tackle this problem for classes of interaction functions satisfying desirable symmetries. As is often the case for models of non-linear consensus or synchronisation on standard networks, we require the process to be invariant to translation and rotation. This is a reasonable assumption for physical and sociological interaction processes ensuring independence on the global reference frame. A function is rotational and translational invariant if it is invariant under application of elements from the special Euclidean group SE(N ), which is defined as the symmetry group of all translations and rotations around the origin. As we restrict the scope to scalar values x i on nodes, and do not consider vectors here, the rotational invariance simply means invariance under a change of signs of the values. In the case of twobody dynamical systems, it is known that a necessary and sufficient condition for these symmetries to be satisfied is the quasi-linearity of the interaction [23], that is where k ij is an arbitrary function from R to R. This form implies that the interaction function is, for each edge, an odd function of (x j − x i ), which is a popular choice in the study of non-linear consensus [26]. Within the language of non-linear consensus, this model belongs to the family of relative non-linear flow. Also note that, as long as this function is positive-definite, the dynamics result in a time rescaled process, which will finally lead, for any initial condition, to a consensus in which all node states have the same value. While we cannot generally transfer these results to three-body dynamical systems, they will provide us a guide on how to define a 'minimal non-linear' model in that case.

B. Linear Dynamics and Motif Matrices
We first investigate the relations between two-body and three-body dynamical systems in the case of linear interaction functions. Linear dynamics is obviously of critical importance as it often serves as a first approximation and determines the linear stability of critical points for non-linear systems. With two-body dynamical systems, the interaction function is given by is a scaling constant, and the resulting dynamics readṡ where L ij = D ij − A ij is the network Laplacian. Here the degree matrix D ij = δ ij k i is a diagonal matrix of the degrees k i = j A ij . Eq. (8) naturally arises when modelling continuous-time random walks on networks [25], but also for opinion-formation and decentralized consensus, as in the continuous-time DeGroot model [27]. For undirected, connected networks, the dynamics asymptotically converges to consensus, with a rate determined by the second dominant eigenvalue of the Laplacian. With three-body interaction systems, the linear interaction function is a linear combination of the diffusive couplings on the two edges (i, j) and (i, k). After accounting for the symmetry in j and k, one finds the linear interaction function where c is again a constant. The three-body dynamical system simplifies enormously in this case, aṡ is the triangle motif Laplacian, defined as the standard Laplacian of a network whose adjacency matrix is This rescaled network is thus obtained by weighting each edge by the number of triangles to which it belongs. The diagonal degree matrix counts the number of triangles around a node. In other words, a three-body dynamical system can be rewritten as a two-body dynamical system in the case of linear dynamics, after a proper rescaling of the adjacency matrix. This observation reveals that a genuine multi-body dynamics on hypergraphs requires a non-linear interaction function, and that multi-body interactions are thus not sufficient to produce dynamics that can not be reduced to two-body dynamical systems. It is therefore essential to consider the interplay between the interaction function and the multi-body topology of the system. Before exploring further the dynamics, let us note the connection between the operation (11) and the 'motif matrix' used to uncover communities in higher-order networks [28]. Motifs are an important object of study in network science [29] and are defined as follows. A motif on k nodes is defined by a tuple (B, P ), where B is a k×k binary matrix and P ⊂ {1, 2, . . . . , k} is a set of anchor nodes. The matrix B encodes the edge pattern between the k nodes. In [28], the authors searched to define a generalisation of conductance and of the cut, where the basic unit of interaction would be a specific motif and not an edge. Essentially, their analysis led to the analysis of a motif adjacency matrix (W M ) ij = number of instances of motifs in M (13) containing i and j, from which a motif Laplacian could be defined. The results from (11) provide a dynamical interpretation of this quantity for triangles, which can readily be generalised to arbitrary k-cliques, and could also open the way to the use of random-walk based community detection techniques, such as the Map equation [30] and Markov stability [31], for higher-order networks.

A. Modelling Consensus Dynamics with Group reinforcement
In this section, we introduce and study a specific form of non-linear interaction function, aiming at modelling consensus dynamics with group reinforcement on hypergraphs. Note that other choices of non-linear interaction functions, akin to Watts threshold models [32], have been considered recently for information spreading [10]. As previously mentioned, multi-body group effects that cannot be reduced to pairwise interactions can appear in various contexts. In the area of sociology, reinforcing group effects such as peer pressure are a long-standing area of study, for instance in social psychology [12]. In the case of binary opinions, for instance, these phenomena emerge when the influence of one opinion on the other in a group depends non-linearly on its popularity. Another important mechanism for opinion dynamics is based on homophily [33], and the fact that only sufficiently similar nodes tend to interact with each other [34]. It is important to develop models that capture these multibody mechanisms to better understand phenomena such as hate communities [35], echo chambers and polarisation [36] in society.
Motivated by these observations, and inspired by Eq. (3), we introduce a three-body consensus model (3CM), with a non-linear interaction function of the form where we assume the function on each triangle is the same, for the sake of simplicity, and with a resulting dynamics for each node i given bẏ This expression models, for each triangle {i, j, k}, the multi-body influence of nodes j and k on node i by the by an influence function s(|x j − x k |) of their state differences. If s(x) is monotonically decreasing, j and k reinforce their influence on i if they have similar states and hinder each other if they are very different. This property is reminiscent of non-linear voter models in the case of discrete dynamics [37][38][39], where the voters change opinion with a probability p f that depends non-linearly on the fraction f of disagreeing neighbours. It is also akin to the Sznajd model modelling the proverb that 'United we stand, divided we fall' in sociophysics [40]. In Eq. (14), the interaction function is non-linear for non-constant s(x) and captures multi-body effects, as the interactions on a triangle can no longer be split into pairwise interaction functions. If s(x) is constant, we recover the linear case discussed in Section III B. Note that our choice of interaction function is clearly symmetric in i and j and therefore congruent with the underlying undirected triangle topology. As rotations are isometries preserving the norm, s(|x j − x k |)(x j − x i ) is rotational and translational invariant. Therefore, this property is also true for f which is a sum of two functions of this form. The value xi ∈ [0, 1] of the nodes is visualised by a colour gradient between white and black. We consider a monotonically decreasing influence function (e.g. s(x) = exp(λx) for λ < 0). As a result, nodes p and k reinforce each other as they have similar values, and hence a large s(|xp − x k |). Nodes k and j have instead distinct values, which leads to a smaller scaling function s(|xj − x k |). As the edge {i, k} is part of both triangles, the diffusion on this edge is scaled by both processes.
A natural choice for the function s(x) is where the sign of the parameter λ determines if the function monotonically decreases or increases and, therefore, if similar or disparate values of j and k accelerate the dynamics on node i. When λ = 0, we recover linear interactions with a constant s(x) = 1. Figure 2 shows the influences on node i for λ < 0, i.e. where similar node states reinforce each other. Another option is the Heaviside function, given by which switches between a zero interaction and linear diffusion when the difference of the neighbouring triangle nodes becomes smaller than a threshold φ ∈ (0, 1). Note that this function is not positive-definite, so that the dynamics do not necessarily converge to consensus, with all nodes have the same value, asymptotically. We will examine the dynamics governed by these two different forms of interaction functions in more detail in our numerical experiments. Until then, we will assume that s(x) is an arbitrary scalar function.

B. Derivation of a Weighted, Time-dependent Laplacian
In Section III B, we showed that, in the case of linear interactions, a three-body dynamical system can be rewritten as a standard two-body dynamical system, defined on a network where the weight of an edge is the number of triangles to which it belongs. Let us explore how this result extends for 3CM. Recall that we assumed that the adjacency tensor A ijk is symmetric. We define I ij as the index-set of nodes that form a triangle with nodes i, j. Note that I ij = ∅ if no triangle exists or, equivalently, if the edge {i, j} does not exist. We now define the weighted adjacency matrix W as The entries (W) ij measure the three-body interactions on node i over edge {i, j}. The matrix verifies diag(W) = 0, because I ii = ∅ as we excluded triangles formed by selfloops. The corresponding degree matrix measuring the total three-body influence on node i is defined as and the corresponding Laplacian is then given by Using (20), we can rewrite the dynamics aṡ The 3CM can thus also be rewritten in terms of the Laplacian of a network, with the essential difference that this network is now time-dependent, through its dependence of the adjacency matrix W = W(t) on the dynamical node states x i = x i (t) in (18). We drop this notation for simplicity, simply using W from now on, but note that the weighted Laplacian is indeed constant in time in the linear case. The weighted Laplacian is the matrix representation of the non-linear dynamics and therefore the analogue of the motif Laplacian, introduced in Section III B for linear dynamics. Even if this reformulation suggests that the dynamical system may be rewritten as a pairwise system, the entries of the weighted Laplacian are dependent on the node states of the neighbouring triangle nodes, reflecting the genuine three-body interactions of 3CM.

C. Conservation of the Average Node State
Our objective is to determine if, and how, 3CM asymptotically reaches consensus. The identification of conserved quantities is usually an essential step to understand the properties of dynamical systems. In the case of consensus on networks, it is well-known that the average state at time t,x is conserved under general conditions. Consider a twobody dynamical system described bẏ The initial averagex(0) is conserved if the derivativė is zero for all times. This is true if the adjacency matrix A ij of the underlying graph is symmetric and the interaction function g(x) is odd, which is verified for quasi-linear dynamics.
Consider now a three-body dynamical system with adjacency tensor A ijk and the interaction function f (x) where the form of the interaction function g ensures that the dynamics are rooted in i and symmetric in j and k, in congruence with our model. Let Π(i, j, k) be the set of all permutations of the three indices. Using this notation, we can conclude that the derivativė ) is zero for all times if A π = A τ for all permutations π, τ ∈ Π(i, j, k) and g(x j (t) + x k (t) − 2x i (t)) + g(x i (t) + x k (t) − 2x j (t)) + g(x j (t) + x i (t) − 2x k (t)) = 0. This is the case for an undirected three-body interactions (such that A is symmetric in all dimensions) if we choose g(x) to be a linear function. We thus conclude that in threebody dynamical systems, linear dynamics conserves the average state of the system.
For non-linear dynamics, the conservation of the average state is in general not verified. For 3CM, where the interaction function takes the specific form , the change in the average state can be written aṡ Trivially, the average state is conserved if W is symmetric. For all i, j, which is only true, for all times, if the influence function is constant, s(x) = c, that is when the dynamics is linear. The weighted matrix W ij = c(W T ) ij is then the motif adjacency matrix scaled by the constant c, in agreement with the results in Section III B. Other symmetries are visible in the weighted matrix: if every entry of the i-th row of W ij is the same, node i has the same incoming multi-body effect from all its neighbours. If the same holds for the i-th column, node i has equal outgoing multi-body effect over all triangle edges. When the dynamics is non-linear and s(x) is not a constant, W(t) and the weighted Laplacian are timedependent. If, at a certain time t, W(t) is symmetric, the derivative of the average state vanishes locally in time. This happens whenever the multi-body effects on the triangles are in balance with each other, as in Eq. (27). This balance is determined by how many triangles a node is a part of, what states the nodes in these triangles are in, and the form of the influence function s(x). Let us assume that the system is initially such that W(0) is symmetric, and examine how the symmetry is affected for later times t > 0. As we are working with undirected and thus symmetric three-body interactions, the dynamics are symmetric on each triangle. The triangle topology and the interaction function itself play a role in the evolving process, as they are encoded in the weighted matrix. Therefore, we need a three-body topology that is symmetric concerning the form of s(x) in addition to a symmetric initial adjacency matrix for this property to be conserved. Otherwise, the average state shifts according to the interplay of influence function s(x), initialisation, and network topology.

D. Influence of Initialisation in the Mean-Field
In order to investigate how the non-linearity of 3CM affects the conservation of the average state, we first consider a mean-field approach. We consider a fully-connected hypergraph H with N nodes V (G) = {1, . . . , N }, such that each triplet of distinct nodes is connected by a triangle. The generalised adjacency tensor is thus We have that I ij = V (G)/{i, j} with |I ij | = (N − 2) for all i, j with i = j and I ij = ∅ for i = j. The dynamics at time t are then given bẏ The average state of those dynamics is represented bȳ x i (t) and is conserved if W(t) is sym-metric. In the fully connected case this condition holds if (29) for all nodes i and j in the network. This means that the multi-body effects of the interactions must be balanced for all nodes in a fully connected system. This condition is the same as in the general case, see Eq. (27), although the mean-field ignores the additional effect of topological distributions of the triangles on the graph. The equality only holds in general for linear dynamics if s(x) is constant. Otherwise, a shift of the initial state may be observed, depending on the initialisation. As an illustration, consider a situation where the number of nodes is even and when the initial values on the nodes is binary, that is either zero or one. Condition (29) will be satisfied only if the initial configuration is balanced, that is when x(0) = 0.5, in which case that quantity is conserved in time. In contrast, if the initial configuration is unbalanced, there will necessarily exist edges over which (29) does not hold, and the average state is expected to evolve in time. If s(|x j − x k |) is given by a decreasing function, that is when similar nodes reinforce each other, the deviation from 0.5 is expected to grow in time, with a drift towards the majority. In contrast, if the influence function is such that dissimilar nodes reinforce each other, one expects to observe a drift to the balanced state 0.5.
To validate these findings, we have performed numerical simulations of 3CM on a fully connected hypergraph of 100 nodes, and used the exponential influence function s(x) = exp(λx) with λ = −1 for a decreasing function, λ = 1 for a growing function and λ = 0 for a constant. Considering different initial distributions of the node states, we have compared: In the first two cases we do not observe any shift in the average state, as expected. However, we do see a shift for the asymmetric initialisation, as shown in Figure 3. The simulations confirm that the average state is conserved for linear dynamics (λ = 0) and multi-body effects only occur for non-linear interaction functions with λ = 0. For λ < 0 we observe a shift towards the majority, resulting in an asymptotic average smaller than the initial value of 0.2 and growing imbalance. For λ < 0 we see the opposite phenomenon, with a shift of the average opinion towards balance. Finally, let us propose a simple method to approximate the asymptotic value of the average state based on the initial configuration. To do so, we estimate the dynamical importance w i of a node i based on the initial configuration as w i (t) = influence of node i total weight in the system The asymptotic value ofx is then obtained by one explicit Euler step of the dynamics from the initial configuration The simulations in Figure 3 also display the predicted value (31), which correctly identifies the direction of the shift.

E. Influence of the Triangle Network Topology
As a next step, let us go beyond the mean-field and consider non-trivial topologies. A critical aspect is the presence of asymmetries in the topology, in particular of two connected clusters. Consider two fully connected clusters that consist of the same number of nodes, and the clusters are connected by a set of triangles. Since triangles are made of three nodes and three edges, each triangle defines an asymmetric connection between the clusters. Two edges connect nodes of different clusters (outer-edges) and one edge connects two nodes within one cluster (inner-edge). We call a connecting triangle between cluster A and cluster B 'directed towards cluster B' if A contains the inner-edge.
The effect of this topological construction becomes clear if we additionally take the initial node states into account. Consider the case of a binary initialisation, shown in Figure 4, with the nodes in cluster A in the initial state x A (0) = 0 and the nodes in cluster B in the initial state x B (0) = 1. We consider a positive-definite, decreasing influence function s(x), so that similar states reinforce each other. Moreover, the configuration is such that there is one single triangle between A and B, and it is directed towards cluster B. Due to the consensus in cluster A and the fact that it contains an inner-edge, the diffusion on the two outer-edges of the connecting triangle is accelerated. On the contrary, the influence is inhibited in the opposite direction, as the two outer-edges damp the diffusion because of the large state difference between the clusters. For this reason, one expects the initial value in A to dominate that in B and thus to determine the asymptotic consensus value. Despite the lack of direction of edges, the process leads to an asymmetric flow of influence from one cluster to the other [41].
In order to analyse this mechanism more quantitatively, we perform numerical simulations on two fully In contrast, the node-state difference between the clusters, and thus on the outer-edges, is maximal, which slows down the effect of cluster B on A.
connected clusters, each consisting of 10 nodes, with the binary initialisation specified above. We then connect the clusters with 80 randomly placed triangles, such that a fraction p ∈ [0, 1] of triangles are directed towards cluster A and the rest towards cluster B. We first examine the influence of the directedness parameter p. For that purpose, we take the interaction function s(x) = exp(λx) with λ = −100, so that pairs of similar nodes exert an overwhelmingly strong influence on other nodes. We show the results of averages over 20 simulations in Figure 5. In (a), we observe a shift in the final consensus value towards the initial value in cluster A or in cluster B, depending on the percentage of triangles directed away from that cluster. This is expected for the above argument, as the connection is maximally directed towards cluster B for p = 0 and towards A for p = 1. The rate of convergence is significantly faster when the initial configuration is very asymmetric, that is extreme values of p.
The asymmetry also influences the rate of convergence towards consensus, as shown in (b). More asymmetric configurations lead to a faster rate of convergence. The simulations also reveal higher fluctuations in the asymptotic state for values close to p = 0.5. This result indicates that the process is sensitive to even small deviations from balance in the initial topology, which can lead to large differences in the consensus value.
In a second set of experiments, we examine the impact of the choice of interaction function. We first consider a maximally directed connection towards cluster B (p = 0), in order to amplify possible differences, and the influence function s(x) = exp(λx). In Figure 6, we plot the consensus value as a function of λ, showing a transition from the initial value in cluster A to that in cluster B, as expected. For λ < 0 (λ > 0), the consensus is biased towards the initial value in A (B). For λ = 0, the dynamics is linear and the initial average statex(0) = 0.5 is conserved. The final consensus value of the two cluster system (with p = 0), dependent on the parameter λ. As the connection is fully directed towards cluster B, it depends on λ if the nodes in cluster A are reinforcing each others' influence which leads to directed dynamics towards cluster B (λ < 0) or inhibit each other (λ > 0) which leads to the contrary effect. Therefore, the consensus value shifts towards the mean of cluster B with growing λ. For λ = 0 we have linear dynamics and the initial average 0.5 is conserved.
We also consider a more extreme choice of interaction function, the Heaviside function, s(|x j − x k |) = H(|x j − x k | − φ) with the threshold φ ∈ (0, 1). This function switches on and off the impact of f {jk} i (x i , x j , x k ) on node i depending on the node state difference |x j − x k |. If the interaction is switched on, it is linear. This property is reminiscent of the bounded confidence model [42], which makes the current framework an extension of it to threebody interactions. In Figure 7, we show the simulation results for φ = 0.2. As the difference between the two clusters is initially larger than φ, the diffusion is only switched on in the direction of the connecting triangle (so towards B for p = 0). Therefore, only nodes of cluster B change their value initially as shown in Figure 7 (a). As soon as the distance between the two clusters is smaller than φ, the dynamics becomes linear and the asymmetry of the dynamics disappears. For p = 0.5 the dynamics are symmetric as the direction of the triangles is balanced. Note that there is no shift in the average in this case, as we previously observed for s(x) = exp(λx). Since the dynamics are simultaneously switched on and are then linear, the different realisations of the directing triangles between the clusters do not impact the dynamics.
In a third experiment, we consider the effect of initialisation on the dynamics. As seen in Section IV, conservation of the average node state can be achieved if the multi-body effects are balanced on each node, which depends on node states and the three-body structure simultaneously. In the mean-field, we had a fully connected, symmetric topology and therefore an initial state average ofx(0) = 0.5 was sufficient. To achieve a balanced initialisation for the clustered topology, the initialisation has to be chosen to reflect the topological structure, and thus the parameter p. This is a difficult task as the system is sensitive to small deviations from a symmetric topol- . We connect the two clusters with 400 triangles, with p = 0 (left), p = 0.5 (right). As before, the value of p quantifies the asymmetry of the topology. We observe that the initial imbalance of the node-states in the clusters is reinforced for p = 0.5, even for very small deviations from an initial average ofx(0) = 0.5.
ogy for certain influence functions, as shown in Figure 5. We would expect to get close to a balance if we initialise the system randomly. We examine the roles of fluctuations in the initialisation in Figure 8, where each node is initially assigned a random number in [0, 1], so that x(0) = 0.5. We use the influence function s(x) = exp(λx) with λ = −100. The value of p determines if we have an asymmetric topology (p = 0) and therefore deviations from the initial mean inside clusters are reinforced (Figure 8, left), or if we have almost symmetric topology (p = 0.5), which leads to a conservation of the mean for small deviations because the system is sufficiently balanced (Figure 8, right). Note that, for p = 0, the system first reaches a local consensus inside each cluster, before attaining the global consensus asymptotically. As a last set of experiments, we investigate the emergence of different time scales in the dynamics, that is a fast convergence of states inside clusters, followed by a slower convergence towards global consensus. To do so, we initialise nodes in different clusters uniformly randomly over separate intervals, such that nodes of cluster initial averages in the two clusters are far apart. The clusters are internally fully-connected and inter-clusterdynamics are generally more inhibited due to the difference of the initial cluster means. As a result, we first observe a fast symmetric dynamics within the clusters towards their means (Figure 9, bottom) and then a slower dynamics between the two clusters ( Figure 9, top). However, the outcome of this process critically depends on the interaction function. For s(x) = exp(λx), with λ = −100, we observe an asymmetric shift towards cluster B for p = 1, as shown in Figure 9 (left). If we consider the Heaviside function instead, the dynamics between the clusters show a similar asymmetry as in the exponential case until the two cluster means are less separated than the Heaviside threshold φ = 0.2. As shown in Figure 9 (right), the dynamics between the clusters then become linear and symmetric.
Together, these results show the importance of the initial distribution of the node states, their alignement with the modular structure of the hypergraph, and also the specific shape of the non-linear interaction function.

V. DISCUSSION
Network science provides a powerful framework for the modelling and description of interacting systems. Its strength comes from its minimalism and generality, dissecting the notion of connectivity into core elements, nodes and edges, that may then be combined to form indirect connections. Yet, interactions may not always be decomposed into a sum of pairwise edges, and the atomic unit of interaction may thus involve more than two nodes. In this paper, we have explored a simple model, 3CM, for consensus in order to identify the impact of three-body interactions on dynamics.
As a first step, we have clarified the difference between the model of the multi-body structure of a system, and the model of its multi-body dynamics. The distinction is most apparent in the case on linear consensus models, whose dynamics can be reduced to a two-body dynamical system even when they are defined on a higher-body hypergraph. In other words, the interaction needs to be non-linear for genuine, non-reducible multi-body dynamical phenomena to emerge. As a second step, we have introduced a non-linear interaction function inspired by models in opinion dynamics, in which two nodes in a triangle either reinforce or inhibit each other's effect on the third node depending on the own similarity. Note that 3CM is quadratic in the node states and that it may be seen as a first-order correction to linear, two-body dynamics. We have shown that the resulting dynamics may lead to a shift of the average state in the system. In the mean-field, we found that this shift only depends on the initial states of the nodes but, in general hypergraph topologies, this dependency is associated to the dominance of certain sub-graphs over others, due to a complex interplay between the form of the interaction function, the topology of the cluster connection and the initialisation of the node states.
This work opens different research directions. First, we found an interesting relation between linear consensus dynamics on hypergraphs and the so-called motif Laplacian [28] proposed for community detection in higherorder networks. This relations suggests a way to generalise dynamics-based community detection methods for hypergraphs [30,31]. Another important extension to this work would be to go beyond three-body interactions and thus consider hyperedges of any cardinality. A natural way to achieve this goal would be, for each hyper-edge, to weight the influence of a node j on node i by its distance to the average opinion in the hyperedge (excluding node i). For three-body interactions, as considered here, this construction reduces to a scaled 3CM aṡ Other generalisations include the study of directed three-body interactions, which may include asymmetric roles to their constituents, but also of stochastic models of opinion dynamics like the majority rule [43], instead of the deterministic models considered here. Another limitation of our work comes from the very constrained hypergraph structures that we used, e.g. connected clusters, and the fact that the initialisation of the node states was aligned with the topology. To better understand the dynamical properties of 3CM in a more practical setting, it would be essential to understand how, given a fixed hypergraph structure for instance obtained from empirical data, the initial configuration of the states would affect the asymptotic consensus. These findings could then be exploited in order to modify balances or imbalances in the dynamical system by seeding (i.e. changing the initial states) or eliminating components (i.e. changing the topology). In a linear model of a two-body dynam-ical system, it is known that asymptotic properties are dominated by the structure, e.g. the mixing time is determined by the spectral gap. In the case of non-linear dynamics, the structure alone is not sufficient, and it becomes essential to understand how the states of the nodes affect notions such as convergence time and lead to an asymmetry between hypergraph components.