Explosive cooperation in social dilemmas on higher-order networks

Understanding how cooperative behaviours can emerge from competitive interactions is an open problem in biology and social sciences. While interactions are usually modelled as pairwise networks, the units of many real-world systems can also interact in groups of three or more. Here, we introduce a general framework to extend pairwise games to higher-order networks. By studying social dilemmas on hypergraphs with a tunable structure, we find an explosive transition to cooperation triggered by a critical number of higher-order games. The associated bistable regime implies that an initial critical mass of cooperators is also required for the emergence of prosocial behavior. Our results show that higher-order interactions provide a novel explanation for the survival of cooperation.

Introduction.The pervasiveness of cooperation in our world has long puzzled researchers [1,2].After all, the natural world (and human society is not an exception) obeys Darwinian selection, which is driven by the self-interest of individuals.In such a competitive world, costly altruistic cooperative behaviours seem inappropriate, since they do not bring any immediate advantage to the cooperators in the brutal fight for the survival of the fittest [3][4][5][6].It is instead more profitable for a self-interested individual to defect (i.e.not participating to the costly altruistic behaviours), taking advantage of the benefits from the actions of cooperators who, in turn, see their sustainability jeopardized by the higher profits of free-riders [7,8].
A well-known theoretical framework for studying the problem of the survival of cooperative traits in human societies is that of social dilemmas or collective action problems in which, given a set of actors, each of them can choose between two strategies, either to cooperate or to defect [9,10].In this context, a defector receives a higher payoff than a cooperator when the two interact, but if everyone adopts the more profitable selfish strategy of defection, the payoff of the agents vanishes [11][12][13][14].Social dilemma scenarios are typically studied in evolutionary game theory [3,7,[15][16][17][18] by implementing games, such as the Prisoner's Dilemma (PD), on structured populations [19][20][21].The underlying structure of a population is usually modeled as a network, where links represent the interactions between pairs of agents [22][23][24].In some cases, the structure of the network has been shown to promote prosocial behaviours through, e.g., mechanisms of network reciprocity [4,25,26], the heterogeneity of the nodes [27][28][29] and the presence of clustering [30].Despite their contributions to our comprehension of social dilemmas, these attempts to consider realistic interaction structures are limited in their representation of real-world systems.The links of a network can indeed only describe pairwise interactions, while the units of a complex system can also interact in groups of more than two.Thus, networks do not allow to accommodate more realistic and general forms of higher-order social interactions.
In the last years, different higher-order mathematical structures, such as hypergraphs and simplicial complexes, have started to be used to represent interactions among three or more units [31][32][33].From contagion processes [34] to synchronization dynamics [35][36][37] and ecological competition [38], many studies have illustrated that higher-order interactions can give rise to the emergence of novel collective behaviors and dynamical patterns not observable in pairwise networks.The first steps have been moved to consider higher-order interactions also in the context of evolutionary games.However, all the works dealing with n-person social dilemmas [39] either still rely on pairwise networks to define group interactions, without the flexibility and generality of real higher-order networks [40], or make too strong assumptions (e.g.regarding the payoff structure) that can be justified only in the particular scenario under investigation [41][42][43][44][45].For example, when hyperedges have been used to describe group interactions at the microscale level [44,45], the payoff associated to each group is a linear function of the strategies of the group members, hence seriously limiting the general representation of the dynamics of social dilemmas.Conversely, when more general payoff structures have been adopted in an extended framework of n-person games, only unstructured (well-mixed) populations composed of groups of the same size have been considered, and with the main focus on games with more than 2 strategies, thus not addressing the case of social dilemmas [46].
In this Letter, we introduce a general framework to extend social dilemmas to structured populations with the presence of interactions in groups of variable size.We model the interactions structure of a population of players as a hypergraph where players are involved in both pairwise and higher-order games represented as hyperedges of different sizes.In this way, the payoff of each player is determined by the strategies of all the players involved in the interaction at once.By comparing extensive numerical simulations of the evolutionary dynamics of the game on random hypergraphs to the analytical results in well-mixed approximation, we find that the higher-order interactions can dramatically change the Nash Equilibria (NE) of the game, allowing cooperators to survive in the Prisoner's Dilemma (PD).In fact, above a critical number of higher-order interactions which depends on the parameters of the game, the dynamics shows an explosive transition to a bistable state, where besides full defection (the only NE in case of just pairwise interactions) a cooperative stable state emerges.Moreover, we found that an initial critical mass of cooperators is also needed to sustain cooperation in the long term: below this critical mass every player becomes a defector, even if the number of higher-order interactions is above the critical threshold.
The model.We consider a population of N players taking part in M strategic interactions, which can either be pairwise or in groups of three or more players.Such interactions are described by a hypergraph H(V, E), where V is the set of N vertices or nodes representing players, and E is the set of M hyperedges [31,32].Each hyperedge e g , with g ∈ 1, • • • , M , is a group (a subset of V) of two or more players interacting in game g.The hypergraph can be represented by a N × M incidence matrix B, whose entry b ig is equal to 1 if the player i is playing the game g, or is zero otherwise.The number of games in which a player i takes part is given by the hyperdegree k i = M g=1 b ig , while the number of players in a game g is the size of the hyperedge q g = |e g | = N i=1 b ig .We focus here on the case of hypergraphs with hyperedges of size two (2-hyperedges, or simply edges) and three (3-hyperedges), respectively corresponding to classical pairwise games (2-games) and games played in groups of three players (3-games).Concerning the payoffs, given that we have q g players involved in a game g, if we indicate as n s the number of different strategies available, the identical players (or symmetry) requirement reduces the total number of different payoffs n qg s q g to just n s ns+(qg−1)−1 (qg−1) payoffs (see SM).As in the classical pairwise symmetric games, here we consider only n s = 2 possible strategies, cooperation (C) and defection (D).This means that for symmetric 2-games there are 4 possible different payoffs while for 3 players there are 6.As usual, the payoffs for 2-games can be displayed as a 2 × 2 matrix Π, whose element π sisj = π si (s j ), π sj (s i ) is the pair of payoffs for player i and j respectively, when the first player plays strategy s i and the second s j .Generalizing this approach to the case of interactions in groups of three players, the payoffs for a 3-game can then be represented as a 2 × 2 × 2 tensor T , whose element τ sisj s k = τ si (s j , s k ), τ sj (s i , s k ), τ s k (s i , s j ) is now a 3-tuple with the value of the payoff for each of the three players i, j and k, playing strategies s i , s j , s k .Fig. 1 illustrates how to implement social dilemmas on higher-FIG. 1. Higher-order games on a hypergraph.The orange triangular areas are hyperedges of size qg = 3, corresponding to games played by three players (3-games), while the purple segments are pairwise interactions, hyperedges of size qg = 2, representing standard games played by two players (2-games).The payoff structures of symmetric 2-games and 3-games are reported in the two boxes.
order systems.The complete payoff structure for both 2-games (q g = 2) and 3-games (q g = 3) is shown, using different symbols to denote different values of payoffs.As commonly done in the study of social dilemmas, without loss of generality we choose the payoff for mutual cooperation equal 1, while the payoff for mutual defection is equal 0, for both 2-games and 3-games [42].In a similar manner, i.e. independently from the number of players (2 or 3) in the game, with and we indicate the payoffs received for unilaterally deviating from mutual cooperation and defection respectively.In this way it is immediate to identify in and in the payoffs usually denoted, in pairwise social dilemmas, as the temptation T and the sucker 's payoff S. Identifying T and S is a fundamental step for the characterization of the game.According to the values of T and S, classical pairwise games are classified into four different types, each characterized by a different set of Nash Equilibria (NE): the Prisoner's Dilemma (T > 1, S < 0), the Chicken game (T > 1, S > 0), the Stag Hunt game (S < 0, T < 1) and the Harmony game (S > 0, T < 1) (see SM). Hence, we can now extend the same classification to 3-games.In this case there are two additional payoffs, namely for defection against a cooperator and a defector ( namely W), and for cooperation against a cooperator and a defector ( namely G).According to the relative value of these two additional payoffs (if G > W or G < W ) each type of 3-games is divided in two disjoint subsets with different Nash Equilibria (see SM).
Results.To investigate the effects of higher-order interactions on the equilibria of a system with N players, we considered the following evolutionary game dynamics.We start from a population with an initial fraction ρ 0 = ρ(t = 0) of cooperators.At each time step, one player (namely the focal) is selected at random, and a second player (the model) is chosen at random among the neighbouring nodes of the focal player on the hypergraph, i.e. those nodes which are connected to the focal player by hyperedges of any size.Each of the two selected players plays a 2-game with all its neighbors connected through a 2-hyperedge, and a 3-game for each 3-hyperedge it takes part in.A 2-game is completely defined by the values of the payoff matrix entries T and S, while the 3-game has the same T and S of the 2-game, but is also defined by the payoffs G and W .For each game, the focal (respectively model) player earns a payoff depending on its strategy and on the strategies of the other players involved in that particular game (i.e. of the other players in the hyperedge).The sum of all the game payoffs defines the total payoff π f of the focal player and the total payoff π m of the model player.The focal player has then the possibility to adopt the strategy of the model player s m , with a probability which is a non-decreasing function of the total payoff difference π m − π f , modelled as a Fermi function [9,17,47]: where w represents the strength of selection.Since we are interested in the Nash Equilibria of the game we iterate the dynamics to compute the quasistationary (QS) probability distribution [48,49] of the fraction of players adopting strategy C (cooperators), whose maxima correspond to the Evolutionary Stable States (ESS) [50], ρ * , of the evolutionary dynamics [41,51,52].As for the underlying structure of interactions, we have considered random hypergraphs with different numbers of higher-order interactions.We have constructed hypergraphs of order N with tunable average hyperdegree k = N i=1 k i /N and fraction of 3-hyperedges δ = n ∆ /M , where M = n ∆ + n / is the sum of the total number of n / 2-player and n ∆ 3-player interactions in the hypergraph (see SM). Fig. 2 shows the results for the case of the Prisoner Dilemma (PD), the most relevant game in the study of social dilemmas.We recall that the pairwise PD is defined by payoff values T > 1 and S < 0. In particular, for our simulations we chose T = 1.5, S = −0.5 and for the strength of selection w = 1/ k (see SM).For the 3-game PD we consider G and W such that (G − W ) > 0, since in this case the one-shot 3-game has 4 different NE: full defection (D,D,D) and all the permutations of 2 cooperators and 1 defector (see SM).In the two top panels we show ρ * , the ESS of the Replicator dynamics (RD), as a function of the fraction δ of 3-hyperedges, for different values of a := 2(G − W ) and of k , the average hyperdegree of the hypergraph.The colored symbols represent the numerical results for ρ * , obtained from the peaks of the QS distribution p QS (ρ) in panels (c-f).We observe a bifurcation in the stable points of the dynamics when the fraction δ of 3-hyperedges exceeds a critical value δ c (a).In particular, while for δ < δ c the only stable NE is full defection ρ * D = 0, as in the standard pairwise PD, for δ > δ c we observe the emergence of a bistable behaviour where cooperation survives: besides the full defection ρ * D , a new stable state 0.5 ≤ ρ * + ≤ 1 appears due to the effect of the higher-order interactions.Fig. 3(a) illustrates the typical time evolution of the system.It reports the fraction of cooperators ρ(t) as a function of time for 20 different initial conditions characterized by different initial values ρ 0 = ρ(0).We notice that when ρ 0 is smaller than a given threshold ρ * − , the dynamics will converge to the full defection state.Conversely, when ρ 0 > ρ * − , it will converge to the stable state ρ * + where a finite fraction of the population are cooperators.In other words, ρ * − represents the initial critical mass of cooperators needed for cooperation to survive in the long term.Fig. 3(b) shows that ρ * − is a decreasing function of δ for any value of the parameter a.This implies that smaller initial densities of cooperators are sufficient to sustain stable cooperation in systems with a larger fraction δ of 3-games interactions (see SM for further details on the numerical results).
To better characterize the role of higher-order interactions on the outcome of the game, we have solved analytically the game in the well-mixed population case, where each player can interact with all the others.Each interaction is considered to be a 3-game, with probability δ, and a 2-game with a probability 1−δ.The evolutionary dynamics of the fraction ρ of cooperators for a well-mixed population in the thermodynamic limit is given by the mean-field Replicator Equation (RE) [7,53,54]: where π C and π D are respectively the expected payoff of a cooperator and of a defector and are both functions of the density of cooperators ρ and of the fraction δ of 3-games (see SM). Hence, the expected payoff difference is also a function of ρ and δ: It follows immediately that when ∆ = (cδ − b) 2 + 4S(b + S) ≥ 0, then ρ * ± are real-valued for every a, b, δ, S. In particular, given that (cδ − b) 2 is always positive, a sufficient condition for the existence of real-valued solutions is 4S(b + S) = 4S(T − 1) > 0, which is always satisfied for the Stag Hunt game and Chicken game.For the Prisoner's Dilemma and the Harmony game instead ∆ > 0 holds only for certain values of the parameters.In particular, for the game we are focusing on in this Letter, namely the PD, we have T > 1 and S < 0, hence b = T − S − 1 > 0.Moreover, c = a + b > 0, given that we are considering 3-games with a = 2(G − W ) > 0. In this case, we find that ρ * ± are real-valued for: It is easy to prove (see SM) that if δ < δ th 2 the real-valued solutions ρ * ± are negative, while if δ > δ th 1 , ρ * ± are always positive and such that 0 ≤ ρ * − ≤ 0.5 ≤ ρ * + ≤ 1.The stability analysis of the solutions yields that, while ρ * D = 0 and ρ * + are stable, ρ * − and ρ * C = 1 are unstable stationary states.Therefore, Eq. ( 4) gives us the mean field critical threshold δ th 1 of 3-player interactions for cooperation to survive in the higher-order PD.In fact, if the number of 3-player interactions is below this critical threshold the only stable stationary state of the higher-order PD is full defection ρ * D = 0, as in the pairwise case.If instead the fraction of 3-games δ exceeds δ th 1 , an explosive transition to a bistable state emerges, where both ρ * D = 0 and 0.5 ≤ ρ * + ≤ 1 are stationary stable states.In Fig. 2 the analytical mean-field results are reported as dashed lines.In particular, the analytical predictions for the stable states ρ * + and ρ * D are in perfect agreement with the peaks of the quasistationary distributions in Fig. 2(c-f) and with the symbols in panels (a-b) reporting the ESS obtained numerically on random hypergraphs.At the same time, also the critical fraction of 3-games δ th 1 (vertical dashed lines in Fig. 2(a-b) ) correctly marks the discontinuous transition to bistability observed numerically.Fig. 3 displays the unstable solution ρ * − , which defines the basins of attraction of the two stationary stable states ρ * D and ρ * + , showing again a good agreement between the meanfield predictions (dashed lines) and the numerical results (trajectories in Fig. 3(a) and symbols in Fig. 3

(b)).
Conclusions.In this Letter, we have introduced a general game theory framework to study social dilemmas in systems where not only pairwise but also higher-order interactions are possible.The main finding of our work is that cooperation can survive even in cases, such as the PD, where pairwise interactions would lead to full defection.Moreover, the observed transition to a state with a stable fraction of cooperators is explosive when the number of higher-order interactions of the system is above a critical threshold that depends on the parameters of the game.However, the observed bistability implies that, even when possible, the survival of cooperators is not guaranteed: a critical mass of initial cooperators is in this case needed to achieve stable pro-social behaviour.This is in agreement with empirical observations regarding the critical mass of initiators required to trigger social and cultural changes [55,56].Our findings demonstrate that higher-order interactions can promote cooperation in competitive environments, showing a new way out of social dilemmas.While in this Letter we have been focusing on the PD, our higher-order framework can be easily applied to any other game.So we hope, our work will inspire new research on the investigation of higherorder interactions and their effects in different strategic scenarios.

NUMBER OF DIFFERENT PAYOFFS
We consider a general q g -person game, where each player can choose among n s different strategies.If the players are distinguishable (i.e.not identical) then there are n qg s possible different elements of the payoff tensor and for each element, there are q g possible different individual payoffs (since each player is different).That is, in the case of distinguishable players the maximum number of different payoffs N max π is: Instead, if the players are identical, the payoff of a player depends on its strategy and on the unordered sample of the strategies of the other q g − 1 players.Unordered because, since the players are identical, it does not matter which player plays which strategy.Given that there are n s possible different strategies, by applying the formula for unordered sampling with replacement of q g − 1 items picked at random from n s choices, we find that for identical players (i.e. for symmetric games) the maximum number of different payoffs is: Substituting, q s = 3 and n s = 2, we get N max π = 6.In our model we set the payoff for mutual cooperation R = 1 and that for mutual defection P equals 0. The remaining four payoffs are then denoted as , , , and .

CLASSIFICATION OF SOCIAL DILEMMAS
In a social dilemma each player can choose between two strategies, either to cooperate (strategy C) or to defect (strategy D) [9,10,19].Defecting brings a higher individual payoff than cooperating when facing one or more cooperators (i.e. it is convenient to free-ride the cooperative efforts of the other players).However, if all players defect, everyone (defectors included) suffers, since the collective payoff vanishes.For pairwise (i.e.2-player) social dilemmas, the payoffs for mutual cooperation (namely Reward, R) and mutual defection (namely Penalty, P ) can be set to 1 and 0 respectively, without loss of generality.Moreover, the payoff associated with unilaterally deviating from mutual cooperation is T (Temptation), while a player receives the payoff S (Sucker) for deviating from mutual defection.It follows that if S > 0 it is convenient for a rational player to deviate from mutual defection, while if T > 1 it is preferable to deviate from mutual cooperation.Therefore, depending on the combination of values of T and S (i.e, above or below the threshold values 1 and 0), we get four scenarios that depict four possible different games.These games are characterized by different Nash equilibria and can be conveniently represented as a square of games as shown in Fig. S1.In particular, we have a Prisoner's Dilemma for T > 1, S < 0, a Chicken game for T > 1, S > 0, a Stag Hunt game for T < 1, S < 0, and T < 1, S > 0 define a Harmony game.A pairwise game can be represented using the so-called payoff matrix representation, where the element of the matrix π sisj = π si (s j ), π sj (s i ) is the pair of payoffs for player i and j respectively, when the first player plays strategy s i and the second s j [5].For 3-player games, the payoff matrix is then substituted by a 2 × 2 × 2 payoff tensor T , whose element τ sisj s k = τ si (s j , s k ), τ sj (s i , s k ), τ s k (s i , s j ) is now a 3-tuple with the value of the payoff for each of the three players i, j and k, playing strategies s i , s j , s k .As seen in the first section of the SM, the number of possible different payoffs for 3-player symmetric games is equal to 6. Consistently with the pairwise social dilemmas, we choose the payoff for full cooperation (i.e.strategy profile (C, C, C)) equal to (1, 1, 1) and the payoff for mutual defection (strategy profile (D, D, D)) equal to (0, 0, 0).As shown in the manuscript, in a 3-player game the payoff for unilaterally deviating from mutual cooperation (respectively the payoff for deviating from mutual defection) is analogous to the temptation payoff T (respectively sucker's payoff S) in the pairwise social dilemma.Figure S2 shows the full 2 × 2 × 2 payoff tensor.In fact, as in pairwise games, in 3-person games, it is advantageous for a rational player to deviate from FIG. S2.Payoff tensor for 3-player social dilemmas.Consistently with the notation for pairwise games, we assume the payoffs for mutual cooperation and defection respectively equal to (1, 1, 1) and (0, 0, 0).The two matrices represent the two levels of the 2 × 2 × 2 payoff tensor, whose elements are the triplets of payoffs (τ1, τ2, τ3).The top matrix shows the payoffs when player 3 adopts cooperation, i.e. for s3 = C.The matrix on the bottom reports the payoffs for the case s3 = D.It is worth noticing that despite the game being symmetric it would be not obvious to reconstruct the whole payoff tensor just from the payoffs of player 1, as usually done in the case of pairwise symmetric games.
mutual cooperation if > R = 1, while it is beneficial to deviate from mutual defection if > P = 0.However, unlike the pairwise games, in 3-person games there are two additional payoffs ( and ) that define a new threshold.In particular, if > , it is favorable to be a defector when playing against a cooperator and a defector, while if < , it is convenient to side with the cooperator.

Classification of 3-player games
We saw that, given their definitions, the payoffs and can be effectively regarded as the 3-player extension of T and S, and hereafter to avoid confusion we will refer to them with the same letters of the pairwise case (i.e. as T and as S).Hence, we can extend to 3-games the same classification based on the values of T and S of pairwise social dilemmas.That is, for the 3-player games, according to the values of T and S we have the four social dilemmas we saw for the pairwise case.However for 3-person games, depending on whether > or < , each of these four games is now further divided into two disjoint subsets with different NE as shown in the following list.

3-player Stag hunt game
The Stag hunt game is defined by T < 1 and S < 0. In this case, the values of the payoffs and do not change the two NE, (C, C, C) and (D, D, D).However, depending on which payoff between and is higher, the ways in which is possible to reach these two NE changes, and one NE is favored over the other.In Game theory notation, this threshold influences the basin of attraction of the two NE, without changing the NE themselves, i.e. it makes one or the other NE risk dominant: • > : there are more strategic moves leading to (D, D, D) (it has a larger basin of attraction, i.e. it is risk dominant) than to (C, C, C); defection is promoted over cooperation.
• < : cooperation is promoted since there are more strategic paths bringing to (C, C, C).

GENERATING RANDOM HYPERGRAPHS
As for the underlying structure of interactions, we considered random hypergraphs with different numbers of higher-order interactions.We have constructed hypergraphs with tunable numbers n / and n ∆ of 2-and 3-hyperedges, respectively.Let δ = n ∆ /M , where M = n ∆ + n / is the total number of interactions in the hypergraph, be the fraction of 3-hyperedges.For fixed values of N , δ and M , we start with N nodes and first connect each of the possible N (N − 1)/2 pairs of distinct nodes with a probability where k = M/N is the desired average hyperdegree of each node.We then connect each of the N (N − 1)(N − 2)/6 triplets of distinct nodes with a probability Given that every time we add a pairwise interaction the total hyperdegree of the network increases by 2, while when we add a 3-hyperedge it increases by 3, we obtain a random hypergraph with the desired k and δ.If the final hypergraph is not connected, we take the largest connected component.Fig. S3 shows that the networks obtained from this algorithm correctly reproduce the desired numbers of 2-hyperedges, 3-hyperedges, and average hyperdegree.In particular, in Fig. S3.a we notice that the numbers of 2-hyperedges and 3-hyperedges in which each player takes part in are distributed as binomial distributions centered around k = 10, as expected for random hypergraphs given the chosen parameters k = 20 and δ = 0.5 [23].

Stable states
In order to estimate ρ * + , we simulated 1000 runs of the evolutionary dynamics.For each run, we start with a randomly chosen fraction of defectors 0 < ρ 0 = ρ(t = 0) < 1 and we use a different instance of the random hypergraph (generated using the algorithm described in the SM).We use the quasistationary (QS) method [48,49] to evolve the system allowing sufficient time for thermalization.In particular, for our simulations on hypergraphs of size N = 1000, we chose a thermalization time of 10 6 time steps and a total simulation length of 10 7 time steps.We recall from the manuscript that we are focusing on the case of the Prisoner's Dilemma and to define the game we arbitrarily chose the payoff values T = 1.5 and S = −0.5.As for the strength of selection we chose w = 1/ k , however we have verified that the results are consistent for at least one order of magnitude above and below this choice for w.We chose w proportional to 1/ k in order to have a comparable strength of selection among different hypergraphs with different k , since the average payoff in the hypergraph increases with the average hyperdegree.It has been proved that the peaks of the QS probability distribution of players with a given strategy (in our case cooperators) correspond to the evolutionary stable state (ESS) of the system [51].We, therefore, find the peak(s) of the QS probability distribution obtained by averaging the distribution of cooperators over all the 1000 runs.Through this method, if δ < δ c we obtain one peak (corresponding to ρ * D ) of the QS distribution, instead if δ > δ c we find two peaks (ρ * D and ρ * + ) .To estimate the error on our measurements of the ESS, we first found the peaks of the QS probability distribution of each of the 1000 runs, obtaining in this way one peak, (ρ * D ) i , or two peaks, (ρ * D ) i and (ρ * + ) i , for each run i.We then compute the absolute deviations of these peaks from the measured ESS (i.e., from the corresponding peak of the QS distribution averaged over all 1000 runs).The median of these absolute deviations is taken as the error ∆ on the estimate of the ESS, that is: Unstable state (critical mass of cooperators) We measured numerically also ρ * − , the critical mass of cooperators needed to observe stable cooperation.Since ρ * − is the unstable solution of the replicator equation, we need to employ a different approach than the one used for finding the stable solution.First, we divide the 1000 simulation runs (see the SM section on the stable states' numerical results for details on the simulations) into 25 batches consisting of 40 runs each.For each batch i, we found the point (ρ min ) i corresponding to the minimum value of the quasistationary (QS) probability distribution between the two peaks (ρ * D ) i and (ρ * + ) i .To estimate (ρ * − ) i we then integrate the QS probability distribution up to (ρ min ) i .The idea behind this approach is that, if a system starts in an initial condition with fewer cooperators than ρ * − , on average it will end up in (or close to) ρ * D , while if it starts with more cooperators than ρ * − , it will end up in (or close to) ρ * + .Thus, the area under the QS distribution until the minimum (ρ min ) i is proportional to the fraction of initial conditions while if c < 0: It can be shown that δ th 2 < δ th + < δ th 1 and therefore we have positive real-valued stationary solutions 0 < ρ * ± < 1 only for δ > δ th 1 , since if δ < δ th 2 < δ th + the real-valued solutions are negative.

ALTERNATIVE BIFURCATION
In the previous section of the SM, we found that the non-trivial stationary states ρ * ± described by Eq.S13 exist (i.e., are real-valued and positive) iff δ > δ th 1 , where the critical threshold of 3-player interactions δ th 1 is a function of a, b and S.However, the condition given by inequality Eq. (S14) can also be expressed as a critical threshold on one of the other variables a, b, and S. For example, we can easily get a critical threshold on a as a function of δ, b and S: Fig. S4 shows the bifurcation curve as a function of a for various values of δ.We observe a bifurcation in the stable points of the dynamics when a = 2(G − W ) exceeds a critical value a c .In particular, while for a < a c the only stable NE is full defection ρ * D , as in the standard pairwise PD, for a > a c we observe the emergence of a bistable behaviour where cooperation survives: besides the full defection ρ * D , a new stable state 0 < ρ * + < 1 appears due to the effect of the payoffs associated with higher-order interactions.

FIG. 2 .
FIG. 2. Stable stationary states for the PD on random hypergraphs with N = 1000 and tunable ratio δ of three-body interactions.(a) Fraction of cooperators at equilibrium as a function of δ for average hyperdegree k = 20 and different values of a, and (b) for a = 5 and different values of k .(c-f ) Quasistationary distributions for a = 5 and k = 20 and four different values of δ.Symbols represent the numerical results averaged over 1000 independent runs (the error bars are smaller than the symbols), while dashed lines are the analytical mean-field predictions.

FIG. 3 .
FIG. 3. Basins of attraction and critical mass of cooperators for the PD on random hypergraphs.(a) Temporal evolution of the fraction of cooperators for various initial conditions and δ = 0.5, a = 5, k = 20.(b) Unstable stationary state ρ * − as a function of δ for average hyperdegree k = 20 and different values of a. Symbols show the numerical results, while the dashed lines are the analytical mean-field predictions.The shaded areas represent the errors.

where a = 2 (
G − W ), b = T − S − 1 and c = (a + b).Therefore, besides the two trivial stationary absorbing states of the RE, namely full-defection ρ * D = 0 and full cooperation ρ * C = 1, Eq. 1 has other two stationary states ρ * ± for which π C − π D = 0: FIG. S1.Games square showing the four different games defined by the combination of values of T and S above or below the thresholds R = 1 and P = 0.

3 -
player Harmony game Pairwise Harmony games are defined by T < 1 and S > 0. Depending on the values of and we now have the following Nash Equilibria: • > : the game has 4 different pure NE, (C, C, C), (C, D, D), (D, C, D), (D, D, C). • < : (C, C, C) is the only NE of the game.

3 -
player Chicken game The pairwise Chicken game (CG) is defined by T > 1 and S > 0. The values of and characterize two different subsets of CG with different Nash Equilibria as: • > : the NE are (D, D, C), (D, C, D) and (C, D, D). • < : the NE are (C, C, D), (C, D, C) and (D, C, C).
FIG. S3.(a) Probability mass function (PMF) of the hyperdegree k, distinguishing between the contribution to k of 2-hyperedges and 3-hyperedges.(b) PMF of the fraction of 3-hyperedges δ.(c) PMF of the average hyperdegree k of an hypergraph.The PMF are computed over 100 instances of a random hypergraph of size N = 1000.The dotted lines denote the mean of the distributions.The desired fraction of 3-hyperedges and average hyperdegree are δ = 0.5 and k = 20, respectively.

4 FIG
FIG. S4.Fraction of cooperators at equilibrium as a function of a for average hyperdegree k = 20 and different values of δ.Symbols represent the numerical results averaged over 1000 independent runs (the error bars are smaller than the symbols), while dashed lines are the analytical mean-field predictions.For these results we choose T = 1.5 and S = −0.5.