Unstable diffusion in social networks

How and to what extent will new activities spread through social ties? Here, we develop a more sophisticated framework than the standard mean-field approach to describe the diffusion dynamics of multiple activities on complex networks. We show that the diffusion of multiple activities follows a saddle path and can be highly unstable. In particular, when the two activities are sufficiently substitutable, either of them would dominate the other by chance even if they are equally attractive ex ante. When such symmetry-breaking occurs, any average-based approach cannot correctly calculate the Nash equilibrium - the steady state of an actual diffusion process.


Introduction
Social activities, new ideas, and innovative technologies spread through networks of social ties formed by friends, colleagues, and followers in social media, such as Facebook, Twitter, and Instagram. Social ties (in physical and online spaces) are not only a channel through which information flows, but also a channel through which influence is propagated. For example, an individual's decision regarding which activities to join often depends on the fraction and/or 40% adopts activity B, because there is no factor that differentiates between the two, at least on average. In fact, this is not necessarily the case in the simulated diffusion processes. The analytical solutions would indeed be correct if the topological properties of the initially active players (i.e., seed nodes) are symmetric, but such a situation rarely occurs in complex networks. Our numerical experiments reveal that diffusion processes on complex networks generated from the same degree distribution can still reach totally different Nash equilibria when the two activities are substitutes. This suggests that equally-attractive yet substitutable activities do not necessarily gain equal popularity, and either activity may even dominate the other through the cascade of peer effects.
The possibility of "symmetry breaking" also raises an important issue for theoretical studies of network games. With such an unstable equilibrium, any deterministic equilibrium would fail to predict the "true" Nash equilibrium since one of the possible equilibria will inevitably be achieved by chance, depending on the details of network structure that cannot be captured by the degree distribution. Symmetry breaking has long been recognized as a source of diversity in economic development (Matsuyama, 2002;Acemoglu et al., 2017), financial globalization (Matsuyama, 2004), and international trade (Matsuyama, 2013;Chatterjee, 2017). To the best of our knowledge, this is the first study that shows why almost equally attractive activities (or technologies, products, etc) can gain totally different levels of popularity.
In the theoretical analysis, we use stylized networks such as Erdős-Rényi random graphs (Erdős and Rényi, 1959) and random regular graphs, with which we can obtain exact degree distributions. While this greatly facilitates the analysis, these networks are not necessarily realistic (Barabási, 2016). To check how well our model would explain diffusion processes on real social networks, we also examine the goodness of fit using empirical data. Here, we construct a social network of economists based on the acknowledgments of articles published in American Economic Review between 2019 and 2020, where nodes represent authors and edges correspond to social ties created by giving and receiving comments.

Network structure
For theoretical analysis, we construct synthetic networks formed by N players, where N is assumed to be sufficiently large. Player i is connected to k i other players by undirected and unweighted edges. k i is called the degree of player i (or node i). Players at the end of the edges emanating from i are called the neighbors of player i. We consider a whole ensemble of many possible networks in a given class, where a particular network structure is realized with a certain probability. That is, we do not focus on a particular single network, rather we specify the distribution of all possible network structures that would appear in a given network model. Any realized network is, therefore, an instance drawn from the ensemble uniformly at random. Examples of Erdős-Rényi networks (Erdős and Rényi, 1959) and random regular graphs are shown in Fig. 1a.
In Erdős-Rényi networks, any two nodes are connected with probability q, and the model is expressed as G(N, q). Let M denote the total number of (undirected) edges. In the G(N, q) model, the probability that a realized network has x edges is given by where the probability of observing each particular network with x edges is q x (1 − q) ( N 2 )−x . Let q k denote the degree distribution, namely, the probability that a randomly chosen node has degree k. In the Erdős-Rényi model, q k is given by a binomial distribution, so it is approximated by a Poisson distribution for large N : where z = k q k k denotes the average degree. Practically, we need to determine the maximum value for degree k, denoted by k max . Throughout the analysis, we set a sufficiently large k max such that it covers at least 99.9% of the entire distribution, i.e., k max = min{k : k k=0 q k ≥ 0.999}.
In random z-regular networks, nodes are connected uniformly at random subject to the degree constraint that k i = z for every i = 1, . . . , N . The number of edges is thus prespecified as M = zN/2 since every node has degree z. While we mainly use the Erdős-Rényi model in the following analysis, we will also use, when necessary, random z-regular networks to maintain analytical tractability. 4 While we use Erdős-Rényi networks and random z-regular graphs in our analysis, the only network property required for our approach is a locally tree-like structure with which the presence of local cycles can be ignored. For instance, Erdős-Rényi networks are locally tree-like if the connecting probability q = z/(N − 1) is sufficiently small (i.e., z N ) (Newman, 2018). If a network is locally tree-like, the average number of second neighbors, i.e., the nodes at distance two from a starting node, is given by k, q kq k = k q k k q +1 ( + 1)/z = z q ( − 1) /z = z 2 , whereq ≡ ( + 1)q +1 /z denotes the excess degree distribution 5 , and we used the fact that the variance of degrees is given by z when the degree distribution is Poissonian. This indicates that each of the z first neighbors of a starting node connects to z second neighbors, on average, and the first and second neighbors of the same starting node will not overlap (Fig. 1b). To give an intuition, consider a representative situation in which a starting node has z (first) neighbors and one of the neighbors has an excess degree z. Since nodes are connected uniformly at random, the probability of a first neighbor being connected to at least one of the other first neighbors (e.g., dotted line in Fig. 1b denotes the probability that one of the other z − 1 first neighbors is chosen among N − 2 nodes, excluding the starting node (black node) and the focal first neighbor (gray node). This suggests that the chance that a local triangle is formed will be vanishingly small if z N , where it will be unlikely that the first neighbors of the starting node are connected with each other. Since this argument holds for any starting node, the neighbors of a first neighbor are also unlikely to be connected with each other, and the neighbors of a second neighbor are unlikely to be connected with each other, and so on. 6 A more general class of network models, called configuration models, in which the degree distribution is prespecified while nodes are connected at random subject to the degree constraint, also exhibit a locally tree-like property (Molloy and Reed, 1995;Newman, 2018). Clearly, Erdős-Rényi networks and z-regular networks are special cases of the configuration models such that the degree distributions are specified by a Poisson distribution and q z = 1, respectively.

Coordination games with a bilingual option
Throughout the paper, we consider two types of activities, denoted by A and B. Players will benefit from an activity if they enjoy the same activity as their neighbors. The two activities may be complements (e.g., drinking and smoking), substitutes (e.g., committing a crime and taking higher education), or neutral. For each of the two activities, players face a binary problem for which they adopt either action 0 ("don't do it") or action 1 ("do it"). The strategy set is thus given by {00, 01, 10, 11} ≡ S. Strategy 11 is called the bilingual option where the player engages in both activities (Oyama and Takahashi, 2015;Arigapudi, 2020). Each player selects a pure strategy s ∈ S, taking all the neighbors' actions as given.

Pure strategy equilibria in a bilateral game
The payoff matrix of the 4 × 4 bilateral coordination game is given in Table 1. A player receives payoff a > 0 (resp. b > 0) of coordinating with a neighbor if they both join activity A (resp. B). The cost of participating in an activity is given by c, where 0 < c < a and c < b, so the net payoff of coordinating on activity A (resp. B) leads to a − c > 0 (resp. b − c > 0). If two players successfully coordinate on both activities, both players receive a + b plus extra payoffs δ, for a 5 The excess degree distributionq is the probability that the total degree of a randomly selected neighbor is + 1, i.e., the probability that a randomly selected neighbor of a starting node has edges except the edge emanating from the starting node (Newman, 2018). is called the excess degree.
7 In Oyama and Takahashi (2015), they consider an extra cost of taking a bilingual option, which needs to be incurred independently of the other player's response. In the current payoff structure, we do not consider such an extra cost, where the cost of taking a bilingual option is given by the sum of the costs for each action.
Instead, we introduce an extra payoff, δ, of coordinating on the bilingual option.
The conditions (4)-(6), combined with the four assumptions, define pairwise risk-dominant equilibria depending on the relative sizes of the parameters, as presented in (3).
When δ is negative and large in absolute value such that δ < c − a or δ < c − b, the action pair (11, 11) is no longer a Nash equilibrium since joining both activities is not beneficial. In the following analysis, we employ the assumption that δ > c − a and δ > c − b to focus on the situation in which the bilingual option can be an equilibrium strategy.

Multilateral games with k neighbors
In games with k neighbors, the total payoffs of a player are given by the sum of the payoffs received in the k bilateral games. Let v(s, m) denote the total payoffs of a player who adopts strategy s ∈ S and faces the neighbors' strategy profile m = (m 00 , m 01 , m 10 , m 11 ) , where m s ∈ Z ≥0 denotes the number of neighbors that adopt strategy s ∈ S. Note that we have s∈S m s = k for nodes with degree k. For a given m, the payoff of each strategy leads to v(00, m) = 0, v(11, m) = −2ck + am 01 + bm 10 + (a + b + δ)m 11 .
The optimal strategy s * for a given m, is then given by Since v(00, m) = 0, players adopt strategies other than 00 if the total payoffs of those strategies are positive. 8 The conditions for v(01, m) > 0 and v(10, m) > 0 are given by the following simple threshold rules: It should be noted that if there were only one activity, say activity A, the players' behavior would be ruled by a fractional threshold rule, i.e., m 01 /k > c/a, as in the well-studied models of contagion (Morris, 2000;Watts, 2002;Jackson, 2008), in which each player adopts either strategy 00 or 01. Since we have a bilingual option here, the optimal strategy is not determined simply by the two conditions (12) and (13).
The conditions v(11, m) > v(01, m) and v(11, m) > v(01, m) are respectively rewritten as v(11, m) > v(01, m) iff m 10 k Suppose for the moment that δ > 0. Conditions (13) and (14) indicate that the fraction of neighbors joining activity B needed for strategy 11 to be preferable to strategy 01 will be less than that required for strategy 10 to be preferable to strategy 00. This is because when δ > 0, the payoff of engaging in both activities is greater than the sum of the payoffs of each activity. However, if the two activities are substitutes and thereby δ < 0, the benefit of joining an additional activity is diminished, so the condition for joining activity B in addition to A will be more stringent than condition (13), which is the condition for deciding whether to participate in B or do nothing. The same argument also holds for the relationship between v(11, m) and v(10, m) (Eq. 15). A crucial difference from the bilateral coordination games is that in this multilateral environment, coordinating with neighbors may not necessarily be the best response. For example, suppose that m 11 = 0, while m 01 and m 10 are sufficiently large such that Eqs. (12) and (13) are satisfied (e.g., most friends use Windows or Mac OS, but none of them use both). Since m 11 = 0, conditions (14) and (15) reduce to m 10 /k > c/b and m 01 /k > c/a, respectively, and these conditions are also satisfied since Eqs. (12) and (13) hold for m 11 = 0. Therefore, the best strategy turns out to be s * = 11 (e.g., using both Windows and Mac), regardless of the fact that no neighbor adopts s = 11.
It is also straightforward to show that v(01, m) > v(10, m) iff am 01 − bm 10 + (a − b)m 11 > 0, which suggests that a player's strategy may switch from 01 to 10 or vice versa in the process of diffusion, depending on the relative size of m 01 and m 10 . The diffusion process is thus generally non-monotonic, unlike the standard (irreversible) binary-state cascade models (Morris, 2000;Watts, 2002;Jackson, 2008;Unicomb et al., 2021).

p-dominance and its relation to the threshold conditions
In the previous subsection, we obtained the fractional threshold conditions for a strategy to be the best response. Here, we argue that a commonly used equilibrium selection criterion, p-dominance (Morris et al., 1995;Kajii and Morris, 1997), can be interpreted as a sufficient condition for the corresponding threshold condition.
Definition (Morris et al. 1995, Kajii andMorris 1997). Letũ(s i , s j ) denote the (i, j)th element of the payoff matrix, and let π on S be a probability distribution where s∈S π(s) = 1. The strategy pair (s i , s i ) is p-dominant if for every probability distribution π on S such that π(s i ) ≥ p, and for every s ∈ S, A Nash equilibrium (s i , s i ) is said to be p-dominant if s i is a best strategy as long as the neighbors adopt strategy s i with a probability greater than or equal to p. While p-dominance is originally an equilibrium characterization for games with incomplete information, this concept is closely related to the threshold conditions shown in the previous section. From Eqs. (12), (14) and (16), one can show that s = 01 will be the best response if the fraction of 01-neighbors, m 01 /k, satisfies all of the following three conditions: where we used the inequality m 01 + m 10 + m 11 ≤ k. Now, consider more stringent conditions for which the inequalities (18)-(20) will hold for any configuration of (m 00 /k, m 10 /k, m 11 /k). Such conditions are given by The combined condition of (21) is thus a sufficient condition for s = 01 to be the best response. Indeed, the strategy pair (s, s) is p-dominant for any p ∈ [p s , 1] and s ∈ S such that Note that p 01 is equivalent to the maximum threshold in (21). This implies that if the fraction m 01 /k is interpreted as the probability of a randomly selected neighbor adopting s = 01 (i.e., π(01)), then p 01 would coincide with the relevant threshold value in (21). This argument obviously holds true for the other strategies as well. Fig. 2 shows the values of p s for different parameter combinations.
In section 4, we calculate the dynamical paths for the aggregate share of each strategy. There, we will examine whether a strategy s would increase its popularity once its share exceeds p s indicated in Fig. 2. If the probability distribution of neighbors' strategies is well approximated by the aggregate shares of each strategy in the population, then p-dominance would indeed become a sufficient condition for a strategy to be adopted by the majority of players in the steady state.  3 Analytical framework for describing the dynamics Now we present analytical approaches to describing the dynamics of diffusion. Throughout the analysis, we assume that strategy 00 is the status quo for all players, except for a certain fraction of players who are initially "active." Then, the neighbors of the initially active players (or the "seed players") may respond by adopting optimal strategies other than 00, which may initiate cascades through edges in the network.
Since the strategic choice of each player can be expressed as a function of the neighbors' profile m, we introduce a probability F m (s → s ) in which the strategy changes from s to s : 9 Note that the choice of new strategy does not depend on the player's current strategy s. This indicates that, unlike the standard binary-state models (Morris, 2000;Watts, 2002), the strategic choice is fully reversible and is determined solely in response to the neighbors' profile m. We also call F m the response function.

Approximate master equations
Let V s |m|=k (t) denote the set of k-degree players that belong to the (s, m) class at time t, i.e., all the players in V s |m|=k adopt strategy s ∈ S and have k neighbors. The fraction of k-degree players belonging to the (s, m) class is given by where q k N is the total number of players with degree k. Then, the expected fraction of players adopting strategy s ∈ S (i.e., s-players) at time t is given as where |m|=k denotes the sum over m = (m 00 , m 01 , m 10 , m 11 ) such that s∈S m s = k. Our interest in this section is to calculate the dynamical path of ρ s for each strategy s ∈ S. To this end, we need to describe the dynamics of ρ s k,m , which capture changes in the population in a given (s, m) class.
There are four factors that change ρ s k,m over time. Players will leave the (s, m) class if i) their strategy changes from s to s ( = s), or ii) their neighbor profile changes from m to m ( = m). On the other hand, players will newly enter the (s, m) class if iii) the players' strategies shift from s ( = s) to s, or iv ) the neighbor profiles shift from m ( = m) to m. To take into account all of these four factors that will affect the behavior of ρ s k,m , we employ an approximate master equations (AMEs) approach (Gleeson, 2011(Gleeson, , 2013Fennell and Gleeson, 2019). The dynamics of ρ s k,m is given by the following differential equation: for all s ∈ S and m such that s∈S m s = k ∈ {1, . . . , k max }. We assume that d dt ρ s k,m = 0 for k = 0, meaning that isolated players will not be influenced by other players. φ s (r → r ) denotes the probability that a neighbor of an s-player changes the strategy from r to r for r, r ∈ S: where the denominator k q k |m|=k m s ρ r k,m represents the expected number of (s)-(r) edges that connect s-players and r-players at a given point in time. The expected number of (s)-(r) edges that change to (s)-(r ) in a small time interval dt is then given as k q k |m|=k m s ρ r k,m F m (r → r )dt. The probability of a (s)-(r) edge being changed to a (s)-(r ) edge in the interval dt, denoted by φ s (r → r )dt, is calculated as the ratio of the expected number of edges that changes from (s)-(r) to (s)-(r ) and the expected number of (s)-(r) edges, which leads to Eq. (30).
The first and second terms in Eq. (29) respectively capture the aforementioned factors i) and ii). The first term captures the rate at which the strategy of a player in the (s, m) class shifts from s to s ( = s) in an infinitesimal time interval dt. Similarly, the second term denotes the rate at which the neighbors' profile differs from m. The third and fourth terms correspond to the factors iii) and iv), respectively. The third term captures the rate at which the strategy of the (s , m)-class players changes from s ( = s) to s. The fourth terms shows the rate that the neighbors' profile newly becomes m. e r denotes the 4 × 1 vector with 1 in the r-th element and 0 in the other elements. The expression m − e r + e r thus represents the neighbor profile that has m r + 1 in the r -th element and m r − 1 in the r-th element. It should be noted that Eq. (29) describes the dynamics under asynchronous updates in which only a fraction dt of players can change their strategies in response to their neighbor profiles in a small time interval. Therefore, while the optimal strategy is given as a function of m (Eq. 11), the current strategy of a player does not necessarily have a one-to-one correspondence with m in the process of diffusion. s and m recover the one-to-one correspondence defined by Eq. (11) in the steady state at which players have no incentive to change their strategies, i.e., a Nash equilibrium.
The system of differential equations can be solved by providing initial values ρ s k,m (0) for all s ∈ S and m ∈ {m : |m| = k, k = 0, . . . , k max }. The number of differential equations in the system is calculated as the total number of ways of picking k balls in an urn with replacement and without ordering. In the urn there are balls of four different colors, so the total number of color patterns when one picks k balls is given by 4+k−1 k = k+3 k . Since a player selects one of the four strategies and the degree k ranges from 0 to k max , the total number of differential equations leads to 4 kmax k=0 k+3 k . 10

Mean-field approximation
In the AME approach, it is assumed that the transition rate of a neighbor's strategies, φ s (r → r ), is independent of the states of the other neighbors, and their profile m is used to characterize each player's state. In the mean-field (MF) approach, we impose a stronger assumption that the strategies of neighbors are independent and randomly distributed following a multinomial distribution. Thus, the neighbor profile m is ignored and each player's state is characterized by a combination of strategy and degree, (s, k), rather than (s, m).
To calculate the average fraction of players belonging to the (s, k) class, we define ρ s k as the sum of ρ s k,m over m: The average probability that a neighbor of a player adopts strategy r ∈ S, denoted by ω r , is given by where kq k /z is the probability that a randomly selected neighbor has degree k. The dynamics of ρ s k are then expressed as d dt 10 We solve the system of differential equations using an ODE solver, ode45, for Matlab.
The first term in Eq. (33) captures the rate at which a player leaves the (s, k) class by changing the strategy from s to s ( = s). The second term denotes the rate at which (s , k)-class players newly employ strategy s. Note that since we have four strategies and the degree ranges from 0 to k max , there are 4(k max + 1) differential equations in the MF scheme, where d dt ρ s 0 = 0.

Dominant strategies in Nash equilibria
Let us examine how the steady state of diffusion processes is affected by the payoff parameters (a, b, c), substitutability parameter δ, and network connectivity z. Fig. 3a shows the cascade region calculated using the AME. We find that within the cascade region in which 1 − ρ 00 > 0, there are some distinct subregions characterized by different dominant strategies, where we define a dominant strategy as the strategy employed by more than 50% of players in the Nash equilibrium (Fig. 3a). For a given δ, the dominant strategy varies with the relative attractiveness, a/b, and the average degree, z. When δ = 0, for instance, strategy 01 (resp. strategy 10) prevails only when the payoff value a (resp. b) is relatively large. Which strategy dominates the others also depends strongly on the degree of substitutability δ. Intuitively, strategy 11 will be dominant in most of the cascade region when both activities are complements (Fig. 3a, middle), whereas either strategy 01 or 10 will be adopted when the two activities exhibit strong substitutability (Fig. 3a, right).
It should be noted that shifts in dominant strategies, or phase transitions, occur in a discontinuous manner at the boundary of the dominant regions. A small change in a parameter may cause the current dominant strategy to be discarded in favor of another. (Fig. 3b). In the class of binary-state threshold models, the possibility of global cascades can arise only when the network connectivity (i.e., the average degree z) is neither too weak nor too strong, and the boundaries of the average degree at which phase transitions occur are called critical points (Watts, 2002;Gai and Kapadia, 2010). In our "multistate" cascade model, there are possibly six types of boundaries at which dominant strategies switch: (00, 01), (00, 10), (00, 11), (01, 10), (01, 11), and (10, 11). The critical points for these transition patterns are characterized not only by the network connectivity z, but also by the relative attractiveness a/b and the degree of substitutability δ.
Note that while the theoretical path of (ρ 01 , ρ 10 ) suggested by the MF equations approaches the saddle point, the actual (or simulated) path could deviate from the theoretical path because seed players are not necessarily located in symmetric positions. The extent to which a seed player affects the other players' behavior would depend not only on the number of their direct neighbors, but also on the number of neighbors at two or more steps away. Thus, there is no guarantee that the realized fraction of s-players is equal to the theoretical average ρ s , as there are some fluctuations in the network structure that could not be captured by the current "averagebased" approximation methods. If the realized path in a given network deviates from the MF path at least to some extent, then the subsequent path will move further away from the saddle point. However, because ρ 11 is increasing, the feasible region continues to shrink, resulting in lim t→∞ ρ 01 (t) = lim t→∞ ρ 01 (t) = 0 and lim t→∞ ρ 11 (t) = 1 (Fig. 4c). This indicates that the observed symmetry breaking, if any, is transient, and strategy 11 will always be dominant in the steady state. Fig. 5 illustrates the phase diagrams for substitutes (δ = −2). We see that there is still a saddle point as seen in Fig. 4, but there arises a region in which dρ 11 /dt < 0 (shaded in pale red). In the diffusion of neutral activities, the feasible region always shrinks and thus any deviation from the theoretical path will be diminished, leading to the unique equilibrium (ρ 01 , ρ 10 ) = (0, 0). In the diffusion of substitutes, the same mechanics will still hold if the deviation from the MF path is not sufficiently large (Fig. 5c). However, the feasible region will expand if the deviation from the theoretical path is sufficiently large such that dρ 11 /dt < 0 (Fig. 5b). If this is the case, the observed symmetry breaking will no longer be a transient phenomenon, where the equilibrium for (ρ 01 , ρ 10 ) will be given by (1, 0) or (0, 1) (Fig. 5d). In the following section, we will show that such a persistent symmetry breaking can indeed occur not only in z-regular random networks, but also in more complex Erdős-Rényi networks.
1. For a given z and N , generate an Erdős-Rényi network with connecting probability q = z/(N − 1).
2. Select seed players at random so that there are ρ 01 0 N players adopting strategy 01 and ρ 10 0 N players adopting strategy 10. The other players employ strategy 00 as the status quo.
3. Choose a fraction dt ∈ (0, 1) of players uniformly at random and update their strategies to maximize their payoff v.
4. Repeat step 3 until convergence, where further updates would not change the strategy of any player.

Repeat steps 1-4.
To conduct simulations in a manner consistent with the continuous-time framework, we implement an asynchronous update in step 3, where a randomly chosen fraction dt of the players updates their strategies in an infinitesimally small interval dt (Gleeson, 2011(Gleeson, , 2013. This is consistent with the AME and MF formulations because changes in the fraction dρ s k,m in Eq. (29) (resp. dρ s k in Eq. 33) are explained by a fraction dt of the possible shifts in the players' states captured by the RHS of Eq. (29) (resp. Eq. 33). We set dt = 0.01 in all simulations.

Dynamics of diffusion: symmetric payoffs
Let us first consider the case of symmetric activities, where a = b and ρ 01 (0) = ρ 10 (0). Since a = b, the two activities A and B are equally attractive, so it is expected that ρ 01 (t) = ρ 10 (t) for all t, other things being equal. We find that when the two activities are neutral (δ = 0) or complements (δ > 0), the path of ρ s (t) obtained by the AME method well matches the simulated path for any strategy s, while the MF approximation is generally less accurate (Fig. 6, left and  middle). The standard deviations of the simulated ρ s (t) over 100 runs become vanishingly small as t → ∞. In contrast, when the two activities are substitutes (δ < 0), the standard deviations of the simulated ρ s (t) increase over time, while the average values are still well approximated by the corresponding theoretical values obtained by the AME (Fig. 6, right). The discrepancy between theory and simulation observed in the propagation of substitutes reflects the fact that one of the two activities occasionally dominates the other even if the payoff parameters and the seed fractions are perfectly symmetric. This symmetry breaking occurs through the mechanism described in section 3.4; when the activities are substitutable, players would not have an incentive to choose strategy 11, so once cascade occurs, either strategy 01 or 10 would prevail. The AME and MF solutions suggest the existence of a symmetric equilibrium, but such an equilibrium would not be achieved in numerical simulations where players are connected in a heterogeneous way. Fig. 6 also shows the values of p such that strategy pair (s i , s i ) is p-dominant (denoted by p s i in Fig. 6, middle). As argued in section 2.2.3, if the probability distribution of neighbors' strategies is approximated by the aggregate shares of each strategy, i.e., π(s) ≈ ρ s for all s ∈ S, then p s i would be roughly interpreted as a (sufficient) threshold of ρ s i above which the strategy s i will prevail. We see from Fig. 6 that in the cases of δ = 0 and 2, in which s = 11 prevails, ρ 01 does not exceed p 01 while strategy 01 gained some extent of popularity when t is between 5 and 8. On the other hand, ρ 11 exceeds p 11 around those times and continued to gain popularity until it reaches the steady state. When δ = −2, in contrast, no strategy exceeds the corresponding p to be adopted by the majority of players. Figs. 7a and 7b show the distributions of the difference ρ 01 (T ) − ρ 10 (T ) for different values of δ. When the two activities are complements, we always have ρ 01 (T ) = ρ 10 (T ) (Fig. 7a), in which case the AME solution is accurate for every instance of simulated contagion processes. In contrast, when they are substitutes, we have either (ρ 01 (T ), ρ 10 (T )) = (1, 0) or (0, 1) in most of the simulation runs. This suggests that either activity A or B dominates with probability ≈ 0.5, depending on the details of the structure of the network that the AME or MF methods do not capture. Note that the AME solution will still describe the simulated equilibrium on average, but it does not necessarily mean that the AME solution is accurate for every instance of the propagation processes. In general, there is a negative relationship between δ and the standard deviation of ρ 01 (T )−ρ 10 (T ), and the essential result also holds true for random regular networks (Fig. 7c).

Dynamics of diffusion: asymmetric payoffs
As we saw in the previous section, symmetry-broken diffusion is occasionally observed when the degree of substitutability δ is low (Fig. 7), in which case the approximation methods do not accurately predict the path of simulated ρ s (t). However, when there is a certain extent of intrinsic asymmetry between the two activities, the AME method works quite well even for substitutable activities. If a > b, for example, activity A is always preferred to activity B, other things being equal, and accordingly the number of 10-players will be smaller than the number of 01-players. This makes it difficult to have a situation where the total payoff v(10) is comparable to v(01), which would be achieved if a certain fraction of neighbors adopt strategy 10 while a smaller fraction of neighbors adopts strategy 01. Thus, the unstable nature of the diffusion process that we see for equally attractive activities will not materialize when the intrinsic attractiveness is different. Fig. 8 shows the dynamical path of ρ s (t) when a > b (The other parameters are the same as  6 and b = 4). See caption of Fig. 6 for details. those in Fig. 6). When the activities are neutral (δ = 0, Fig. 8, left), strategy 01 initially spreads to more than half of the players, and then most of the players begin to shift their strategies to s = 11, which leads to a decay in the fraction of 01-players. In the context of the phase diagram in Fig. 4, this corresponds to a situation where the feasible region ρ 01 + ρ 10 + ρ 11 ≤ 1 continues to shrink, and thereby ρ 01 is bounded from above by 1 − ρ 10 − ρ 11 . Note also that ρ 01 did not exceed p 01 , meaning that (01, 01) was not p-dominant.
When the activities are substitutes (δ < 0), strategy 01 dominates the other strategies (Fig. 8, right). Since the total payoff of choosing activity B cannot be comparable to that of activity A in the process of diffusion, the simulated paths are quite stable over different runs, and thus the standard deviations are fairly small. Note that ρ 01 surpasses p 01 during the diffusion process, which implies that the strategy pair (01, 01) is p-dominant.

Alternative model of diffusion
We have seen that players' interactions through 4 × 4 coordination games provide a mechanism by which an activity spreads over a network. To check the robustness of the results, here we consider an alternative model of diffusion in which players select their strategies to maximize utility functions.
Let x i ∈ {0, 1} denote player-i's binary action on activity ∈ {A, B}. The strategy of player i is expressed by vector , and the strategies of the other players are summarized as Note that there is a one-to-one correspondence between x i ∈ {(0, 0), (0, 1), (1, 0), (1, 1)} and s i ∈ {00, 01, 10, 11}. Following Chen et al. (2018), we consider the following quadratic utility function: where g ij ∈ {0, 1} denotes the (i, j)th element of the adjacency matrix G = (g ij ); g ij = 1 if there is an edge between i and j and g ij = 0 otherwise. α is a parameter that captures the benefit of enjoying an activity itself, which entails costs (= 1/2). γ and β respectively denote the benefit of cooperating with neighbors and the degree of substitutability. Thus, the larger β is, the greater the complementarity between the two activities.

Absolute-threshold rules
The utility for each strategy is given by Since the strategy of a player can be considered as a function of the neighbors' strategy profile m = (m 00 , m 01 , m 10 , m 11 ) , we can redefine the utility function asũ i (s, m) ≡ u(s, s −i ). The optimal strategy is then given by where we drop the subscript i for simplicity.
These assumptions guarantee thatũ(s, (k, 0, 0, 0) ) ≤ 0 ∀s ∈ S, which prohibits players from adopting strategies other than 00 when no neighbors enjoy activities. This means that neighbors' influence is necessary for an activity to spread over the network. Sinceũ(00, m) = 0, a player selects a strategy other than 00 as long as the utility is positive. 12 The conditions forũ(01, m) > 0 andũ(10, m) > 0 are respectively given by the following threshold conditions:ũ u(10, m) > 0 iff m 10 + m 11 > 1 γ The most important difference from the threshold conditions in coordination games, Eqs. (12) and (13), is that conditions (41) and (42) are independent of the player's degree k. This indicates that a player is more likely to join activity A as more neighbors join activity A, regardless of the fraction of active neighbors. This also implies that the larger the degree, the more likely it is for an activity to spread over a network. This type of diffusion has been studied within a class of absolute-threshold models (Granovetter, 1978;Karimi and Holme, 2013;Unicomb et al., 2021).
u(11, m) >ũ(10, m) iff m 01 + m 11 > 1 γ Suppose for the moment that β > 0. Then Eq. (43) indicates that the threshold of the number of neighbors enjoying activity B (= m 10 +m 11 ) above which strategy 11 is preferable to strategy 01 is smaller than that required for strategy 10 to be more preferable than 00 (Eq. 42). This is because when β > 0, the utility of enjoying both activities is greater than or equal to the sum of the utilities for each activity. However, if the two activities are substitutable (i.e., β < 0), the incentive to enjoy both activities is discouraged, so the condition for enjoying activity B in addition to A becomes more stringent than condition (42), which is the condition to enjoy activity B or do nothing. The same argument also holds for the relationship betweenũ(11, m) andũ(10, m) (Eq. 44).
Regarding the choice between the two activities, we havẽ Whether this inequality holds depends on the relative number of active neighbors, m 01 − m 10 , which will be time-varying. Thus, a player's strategy may switch from 01 to 10 or vice versa in the process of diffusion, depending on the neighbors' status. Therefore, diffusion processes in this utility-based model are generally non-monotonic, similar to the contagion processes through coordination games.
On the other hand, as the substitutability parameter β becomes negative, there arises a wider region within which dρ 11 /dt ≤ 0 (Figs. 9b and c). In particular, when β = −1, we see that dρ 11 /dt is non-positive for the entire feasible region of (ρ 01 , ρ 10 ). This suggests that the steady state obtained by the AME method is practically not attainable since any path deviating from the stable path (i.e., stable arm) would cause symmetry breaking, leading to (ρ 01 (T ), ρ 10 (T ), ρ 11 (T )) = (1, 0, 0) or (0, 1, 0). Again, this is a situation in which one of the activities would dominate the other purely by chance, depending on the realization of an instance drawn from an ensemble of random networks. When the degree of substitutability is moderate (Fig. 9b), there is a wide region of (ρ 01 , ρ 10 ) in which dρ 11 /dt > 0, so the feasible region may shrink over time. This will make the upper bounds for ρ 01 and ρ 10 well below 1, resulting in a smaller degree of asymmetry compared to the case of β = −1. Fig. 10a illustrates how the accuracy of the AME solutions is related to the degree of substitutability β and the average degree z. We find that when β = −1, each simulation run will be likely to deviate from the AME solution, while the average over multiple simulations appears to be close to the value predicted by the AME. These deviations from the AME values can be interpreted as an outcome of incidental diffusion due to symmetry breaking, whose dynamics are captured by the phase diagrams in Fig. 9c.
Naturally, we find that there is a negative relationship between the degree of substitutability β and the standard deviations of ρ 01 (T )−ρ 10 (T ) (Fig. 10b). When the two activities are strongly substitutable, it is highly likely that either of the two dominates the other (i.e., ρ 01 (T ) = 1 or ρ 10 (T ) = 1), resulting in the standard deviations being close to 1. When the substitutability is moderate, symmetry breaking still occurs, but the difference between ρ 01 (T ) and ρ 10 (T ) would be well below 1 because strategy 11 also prevails to some extent. In the phase diagram in Fig. 9b, this corresponds to an expansion of the area in which dρ 11 /dt > 0.

Asymmetric preferences
Now we analyze the case of asymmetric preferences where there is no possibility of symmetry breaking by definition. In particular, we investigate the relationship between network connectivity and the steady-state share of each strategy. Throughout this section, we assume that activity B is a little less attractive than activity A, where α A = 0.4 and α B = 0.3.
When the two activities are neutral (i.e., β = 0), we see that lower connectivity would lead to a higher chance of strategy 01 being propagated over a network (Fig. 11, left). Under our parameter configurations, strategy 01 is the dominant strategy up to z ≈ 2.5, but for more densely connected networks, strategy 11 dominates and both activities A and B are likely to prevail. Thus, the dominant strategy switches around z ≈ 2.5, at which the simulated cascade sizes can be highly volatile and their average would not match the AME solution unless the number of runs is large enough. 13 When the activities are complements (e.g., β = 0.1), strategy 11 would be dominant in well connected networks such that z > 1 (Fig. 11, middle). Indeed, both activities will spread to most of the players located in the giant component of a network, as is seen in the standard binary cascade models (Watts, 2002). When the activities are strongly substitutable (e.g., β = −1), only activity A would gain popularity, while there would be no chance for B to propagate (Fig. 11, right).
It should be noted that there are some differences between the current utility-based model and the model of coordination games concerning the role of network connectivity. First, in the diffusion mechanism through coordination games, there is an upper bound of z above which propagation decays (Fig. 3). This is because players' thresholds are given by the fractions of neighbors already participating in certain activities, where an increase in the number of neighbors would dilute the influence of each neighbor. However, such a dilution effect does not exist in the diffusion process through utility-based games since the threshold conditions (41)-(45) depend only on the absolute numbers of neighbors adopting certain strategies. Therefore, there is no upper bound of z for the activities to propagate in the utility-based model. Second, while there is a common property that a larger degree of substitutability will lead players to select a more attractive activity, we see an essential difference when the two activities are neutral. In the diffusion through utility-based games, lower connectivity would lead to a propagation of a single activity, as shown in Fig. 11. In contrast, in the diffusion through coordination games, a lower degree of connectivity will enhance the propagation of both activities, while a higher degree of connectivity would allow either activity to spread (Fig. 3). In the fractional threshold model, the influence of a single neighbor will be diluted as the number of neighbors increases, and accordingly, it becomes harder for an activity to gain popularity. In contrast, in the absolute threshold model, a rise in the number of neighbors just facilitates the transmission of peer effects, as it increases the number of routes through which information arrives.

Empirical social networks
Now we examine how accurately the AME approach can describe the dynamics of diffusion in a real-world social network. For this purpose, we construct a social network of economists based on the acknowledgments of the papers published in American Economic Review (AER).

Data: a social network of AER authors
To construct a network, we first create a list of authors who published at least one paper in AER between January 2019 and December 2020. During this period AER published 244 papers. Two authors are connected by an undirected and unweighted edge if one of them gives credit to the other in the acknowledgments of a published paper, assuming that there is a social tie between them. 14 The constructed network consists of some small isolated components, so we focus on the largest connected component on which a global diffusion could occur. A visualization of the network structure is presented in Fig. 12. We have N = 447 with mean degree z = 4.99 and maximum degree k max = 29.

Diffusion on a social network of economists
We examine how well the AME will predict the dynamical paths for each strategy. We consider asymmetric activities in the same way as in sections 4.2 and 5.3 (i.e., a > b and α > β) since the AME works well on synthetic networks when there is intrinsic asymmetry. To calculate the AME solution, we obtain the empirical degree distribution p e k as p e k = # of nodes with degree k N .
We find that in both types of games, the simulated diffusion dynamics are well replicated by the AME approach (Fig. 13). A difference from the simulations based on synthetic networks is that we see larger standard deviations for the fraction of each strategy, indicated by the error bars, while the two activities are intrinsically asymmetric. This would be because the size of 14 To extract author names from the PDF file of each paper that was downloaded manually from the AER web site, we used a Python module Stanza (Qi et al., 2020)  the author network is relatively small, thereby deteriorating the accuracy of the approximation method which is expected to work well for sufficiently large networks. To mitigate the small-size problem, we imposed a relatively high value for the seed fractions, ρ 01 0 = ρ 10 0 = 0.1, yet we still see some fluctuations in the simulated results. Arguably, the standard deviations would diminish if we use larger empirical networks, such as networks aggregated over longer periods. However, the maximum degree k max of such a large network tends to be so high that the computational cost of the AME would be prohibitively expensive.

Conclusion and discussion
We studied the dynamics of diffusion based on two different models of network games. In both models, we showed that the diffusion dynamics exhibit instability due to phase transitions and symmetry breaking. When the activities are complements or neutral, the AME approach is highly accurate in predicting the popularity of each activity at a given point in time. On the other hand, when the activities are substitutable, "average-based" analytical equations such as MF and AME may not correctly describe the Nash equilibria attained in simulated propagation processes.
There are some issues to be addressed in future research. First, while we considered two activities and four pure strategies {00, 01, 10, 11}, one could extend the models to study more than two activities since, in principle, the AME method can be applied for any number of states (or strategies) n. However, the difficulty is that the number of differential equations will increase to n kmax k=0 n+k−1 k , so the computational cost would easily become prohibitive. For instance, if we have three types of activities and a degree distribution over k = 0 to 10, then the number of possible strategies is given by n = 2 3 , and the number of differential equations will be 8 10 k=0 8+k−1 k = 350, 064. Second, it would be useful to introduce nonrandom structural properties, such as clustering and assortativity, which are frequently observed in real social networks (Barabási, 2016;Newman, 2018). Since the AME system is already  complex and there are a number of equations to be solved, introducing additional network properties would be a difficult task. Nevertheless, incorporating structural properties would be important not only to improve the accuracy of the AME running on real data, but also enhances the predictability of diffusion. For example, firms would aim to enhance the popularity of new products through online social networks (Watts and Dodds, 2007), and governments may need to monitor the spread of misinformation and manipulations that could spread through social media (Badawy et al., 2019). Third, it would also be useful to take into account the time-varying aspects of real-world networks. In many social networks, nodes and edges frequently appear and disappear, leading to changes in the network structure. The dynamic nature of networks has been extensively investigated within the framework of "temporal networks" in the field of network science (Holme and Saramäki, 2013). We hope that our study will stimulate further research in these directions.