Conformity, anticonformity and polarization of opinions: insights from a mathematical model of opinion dynamics

Understanding and quantifying polarization in social systems is important because of many reasons. It could for instance help to avoid segregation and conflicts in the society or to control polarized debates and predict their outcomes. In this paper we present a version of the $q$-voter model of opinion dynamics with two types of response to social influence: conformity (like in original $q$-voter model) and anticonformity. We put the model on a social network with the double-clique topology in order to check how the interplay between those responses impacts the opinion dynamics in a population divided into two antagonistic segments. The model is analyzed analytically, numerically and by means of Monte Carlo simulations. Our results show that the systems undergoes two bifurcations as the number of cross-links between cliques changes. Below the first critical point consensus in the entire system is possible. Thus two antagonistic cliques may share the same opinion only if they are loosely connected. Above that point the system ends up in a polarized state.


Introduction
What do affirmative action and gun control [1], same-sex marriage and sexual minority rights [2], abortion [3] stem cell research [4], global warming [5], attitudes toward political candidates [6] or the recent refugee crisis in Europe [7] have in common? All these keywords are examples of topics known to ignite polarized debates in the society. Studying them could thus shed more light on the phenomenon of polarization, which is one of the central issues in the recent opinion dynamics' research. Polarization is understood here as a situation in which a group of people is divided into two opposing parties having contrasting positions [8]. It is sometimes referred to as bi-polarization [9] to distinguish it from the so called group polarization, i.e. the tendency for a group to make decisions that are more extreme than the initial inclination of its members [10,11].
The paper is organized as follows. In the next section we introduce the model and shortly describe how it differs from the one presented in [55]. Then, we will investigate the model both analytically and numerically. Finally, some conclusions will be presented.

Basic assumptions
We begin with a brief description of the assumptions of the model analysed in [55]. We choose the q-voter model [34] as our modelling framework. Within the original model, q randomly picked neighbours influence a voter to change his opinion. The voter conforms to their opinion if all q neighbours agree. If they are not unanimous, the voter can still flip with a probability . The unanimity rule is in line with a number of social experiments. For instance, it has been observed that a larger group with a nonunanimous majority is usually less efficient in terms of social influence than a small unanimous group [57,58]. Moreover, Asch found that conformity (i.e. matching attitudes, beliefs and behaviors to group norms) is reduced dramatically by the presence of a social supporter: targets of influence having a partner sharing the same opinion were far more independent when opposed to a seven person majority than targets without a partner opposed to a three person majority [59].
From social networks analysis it follows that the existence of segments within a network may be correlated to polarization [60][61][62]. We will thus assume that our social network is already modular. For the sake of simplicity we will put the model on the so called double-clique topology consisting of two complete graphs (cliques) connected by some cross links with each other [63]. An example of such network is shown in Figure 1. It should be emphasized that this is a rather strong assumption, because one actually cannot rule out the opposite possibility that segmentation is induced or intensified by polarization. However, analysis of the casual connection between the network segmentation and the polarization is beyond the scope of this work and will be addressed in a forthcoming paper. The original q-voter model uses conformity, i.e. the act of matching attitudes and opinions to group norms, as the only response to social influence. Other possible types of responses has been described for instance within the diamond model [64][65][66][67] and are shown in Figure 2. The anticonformity (i.e. challenging the position or actions of a group) representing negative ties is of particular interest, because from the social balance theories it follows that both positive and negative ties are needed for the polarization to emerge and to prevail [40].Thus we will add this type of response to our model to check how the interplay between conformity and anticonformity impacts its dynamics.
It is worth to note here that the anticonformity as a type of social response was introduced into binary models of opinion dynamics for the first time probably by Serge Galam, who used the notion of 'contrarians' in order to describe agents always adopting opinions opposite to the prevailing choice of others [68]. In his seminal paper he showed that there exist a critical density of contrarians above which a system always ends up in a bi-polarized state.
We are however aware that the assumption on negative influence is still a subject of an intense debate. There are several models able to explain polarization without any kind of negative influence, for example the argument-communication theory of bi-polarization [9] or the bounded-confidence model [16]. Some empirical studies on negative influence do not provide convincing support for this assumption as well [69]. Here presented within a q-voter model framework with q = 4 [34]. The source of influence is a group consisting of unanimous agents (schematically pictured as a cloud). The "up" and "down" spins (arrows) represent agents with two different opinions on a single issue.
There is an important point while speaking about anticonformity: it is relative. It turns out that in many settings multiple sources of norms are possible. As a consequence, conformity to one source can at the same time be anticonformity to another. For instance, a teenager's conformity to peers is very often anticonformity to his parents [66]. Therefore we will assume within our model, that an agent strives for agreement within his own clique and simultaneously anticonforms to individuals from the other clique.
Within the q-voter model all individuals are characterized by a single dichotomous variable, i.e. the model belongs to the class of binary models. At first glance this approach may seem unrealistic, because the opinions of individuals on specific subjects are expected to vary gradually. Therefore they should be rather represented by continuous variables [9,13,14,16]. But from empirical findings it follows that the distribution of opinions on important issues measured on some multivalued scale is often bimodal and peaked at extreme values [70,71]. Moreover, many data on social networks are characterized by a semantic unicity, meaning that opinions and interactions of networks' members are restricted to a single domain or topic [72]. Very often those opinions may be interpreted as simple "yes"/"no", "in favour of"/"against" or "adopted"/"not adopted" answers [73]. In other words, in some situations the most important characteristics of the system under investigation may be already captured by the relative simple models of binary opinions.

The "old" (quenched disorder) model
In this section we recall the model introduced in [55]. We consider a set of 2N agents, each of whom has an opinion on some issue that at any given time can take one of two values: S i = −1 or S i = 1 for i = 1, 2, . . . , 2N. We will sometimes call these agents spinsons to reflect their dichotomous nature originating in spin models of statistical physics and humanly features and interpretation (spinson = spin+person) [35,74].
We put our agents on a double-clique network. It consists of two complete graphs of N nodes connected with L × N 2 cross links ( Figure 1). The parameter L is the fraction of the existing cross links and N 2 -their maximum number. The cross links between the cliques are chosen randomly and the resulting network does not evolve in time during simulation. Thus, from the statistical physics' point of view the model belongs to the class of models with quenched disorder [75].
We will assume that a spinson strives for agreement within his own clique (conformity) and simultaneously challenges the opinions of individuals from the other clique (anticonformity). In other words we link the type of social response of agents with their group identity. To account for the possibility of acting as both conformist and anticonformist at the same time within the q-voter model, we introduce the notion of signals and slightly alter the concept of unanimity of the influence group. A signal is just a state of the neighbour when coming from spinson's clique or its inverted state otherwise. The target of influence changes its opinion only if all members of the influence group emit the same signal. No other modification of the q-voter framework are needed to account for anticonformity.
We use Monte Carlo simulation techniques with a random sequential updating scheme to investigate the model. Each Monte Carlo step in our simulations consists of 2 × N elementary events, each of which may be divided into the following substeps: 1. Pick a target spinson at random (uniformly from 2N nodes). 2. Build its influence group by randomly choosing q neighboring agents. 3. Convert the states of the neighbors into signals that may be received by the target. Assume that the signals of the neighbors from target's clique are equal to their states. Invert the states when from the other clique. 4. Calculate the total signal of the influence group by summing up individual signals of its members. 5. If the total signal is equal to ±q (i.e. all group members emit the same signal), the target changes its opinion accordingly. Otherwise nothing happens.
Thus, our model is nothing but a modification of the q-voter with = 0 and an additional social response of spinsons. You may refer to [55] for further details of the model.

New ("annealed") version of the model
In the model described in the previous section the cross links between the cliques are generated randomly at the beginning of a simulation and remain fixed while the system evolves in time. If the number of cross links is smaller than their maximum number N 2 (i.e. L < 1), some agents may have no connections to the other clique, some others -multiple ones. In other words, the agents may differ from each other because of the distribution of links between the cliques. While it can be handled with ease within a computer simulation, this feature constitutes usually a challenge for mathematical modelling due to the necessity to perform a quenched average over the disorder [56] . Thus, we decided to modify the model slightly for simpler mathematical treatment.
Most of the assumptions presented earlier in this section hold, i.e. we consider 2N agents living on a double-clique network. Each agent may be in one of the states {+1, −1} representing its opinion on some issue. It seeks for agreement within his own clique (conformity) and simultaneously anticonforms to individuals from the other clique (anticonformity). And it changes its opinion if the members of an influence group emit the same signal. Just to recall, a signal is just a state of the neighbour when coming from agents's clique or its inverted state otherwise.
If the size N of each clique is large enough, then the number of links within a single clique is given by With the number of cross connections equal to L × N 2 the quantity gives us the probability of choosing one cross link out of all edges in the double-clique network.
Assuming that every agent from one clique is connected with probability p with an agent from the other clique, and with probability 1 − p with an agent from its own clique we arrive at the new version of our model. Technically this approximation is nothing but an average of the original (quenched disorder) model over different configurations of cross links in the network. Thus it corresponds to annealed models from statistical physics [76]. The step 2 from the update rules of the model defined in paragraph 2.2 requires some adjustments: 1. Pick a target spinson at random (uniformly from 2N nodes). 2. Build its influence group by randomly choosing q agents. In the quenched disorder model we simply followed 4 randomly chosen links of the target to achieve that. Due to the setup of that model some targets usually had no cross connections, some others -multiple ones. Now the situation is different -each target has the same probability of being cross-connected and the actual links to other agents have to be built first. Thus, for each member of the influence group we decide first which clique it will belong to (with probability 1 − p to target's clique, with p to the other one). Then we choose the member randomly from the appropriate clique (see Figure 3). 3. Convert the states of the group members into signals. 4. Calculate the total signal of the influence group. 5. If the total signal is equal to ±q (i.e. all group members emit the same signal), the target changes its opinion accordingly. Otherwise nothing happens. It should be noted that a similar model, but with a broken symmetry between the cliques, was used in both the quenched and the annealed version to explain recurring fashion cycles [77]. Let us denote the state of the i-th agent at the discrete time τ by S i (τ). There are two natural quantities that fully describe the state of the system: the concentration of agents in state +1 and the average opinion (magnetization in physical systems) [35]. The concentration at time τ is defined as where N ↑ (τ) is the number of agents in state +1 at time τ. The average opinion is given by Please note that there is a simple relation between these two quantities: Thus it actually does not matter which one will be chosen for representation of the state of the system. For the sake of convenience we will usually work with the concentration below. However, some of the results will be transformed to average opinions to allow for comparisons with the findings from [55].
In order to easily detect polarized states in the system (i.e. all agents in state +1 in one clique and in state −1 in the other) we will often calculate the concentration separately for each clique: where A and B are the labels of the cliques. The interpretation of c X is as follows: • c X = 1 -positive consensus in clique X, i.e. all agents in that clique are in state +1, • 1 2 < c X < 1 -partial positive ordering in clique X, i.e. the majority of agents is in state +1, • c X = 1 2 -no ordering in clique X, i.e. the numbers of agents in state +1 and −1 are equal, • 0 < c X < 1 2 -partial negative ordering in clique X, i.e. the majority of agents is in state −1, • c X = 0 -negative consensus in clique X, i.e. all agents in that clique are in state −1.

Transition probabilities
In each elementary time step the number of agents in state +1 in one clique -say A -can increase by 1 only if: 1. a target from clique A is chosen (probability 1/2), 2. the target is in state −1 (probability 1 − c A ), 3. it flips, i.e. an influence group emitting signal +q is chosen.
We can immediately write down the transition probability for such an event: We have introduced here a scaled time t = τ 2N and a scaled time step ∆ N = 1 2N . We will use them below to derive a limiting dynamical system for the model. Moreover, one can easily check that the term of the form (u + v) q in the above equation is just the generating function for the probabilities of those compositions of q members of an influence group which can cause an opinion switch event (see Figure 4).
Similarly, the number of spinsons in state +1 in clique A decreases by 1 if: 1. a target from clique A is chosen (probability 1/2), 2. the target is in state +1 (probability c A ), 3. it flips, i.e. an influence group emitting signal −q is chosen.
These conditions lead to the following transition probability: It is also possible that the number of agents in state +1 remains unchanged in elementary time step. The probability of this event is simply given by: (9) After repeating analogous considerations for clique B we get: Thus, given the states of the cliques at time t, the expectations for the numbers of agents in state +1 at time t + ∆ N are given by the following expressions: The abbreviationp = 1 − p was used in the above formulas.

Asymptotic dynamical system
We would like to derive from Eqs. (11) a limiting dynamical system for N → ∞ in scaled time t = τ 2N . Let us first divide the above equations by N: It is very likely that in the limit N → ∞ the random variables c i = N ↑ i N localize and hence become almost surely equal to their expectations. We get Taking the limit N → ∞ and denoting the limiting variables c A and c B by x and y we arrive at

Annealed model as a birth-death process
According to Eqs. (7)-(10) we have only two types of state transitions in each clique: 'births', which increase the state variable by one, i.e. N ↑ X → N ↑ X + 1 (X = A, B), and 'deaths', which decrease it by one, N ↑ X → N ↑ X − 1. Thus our model may be seen as two coupled birth-death processes [78]. Since such a process is relatively easy to simulate, we will use it as an additional benchmark while comparing the results for the quenched disorder and annealed models.

Results
All results presented in this section were obtained via symbolic and numerical calculations by making use of the Python's scientific stack [79]. Python codes needed to reproduce some of them may be found in the supplementary materials.
Although we will often use the case q = 3 for presenting the results, there is no particular reason for choosing this value. We considered in our analysis influence groups of sizes ranging from 1 to 6. The upper bound of the group size was motivated by the conformity experiments by Asch [59]. Qualitatively, the results turned out to be independent on the actual value of q. However, with increasing q the critical points were shifted towards higher values of the interaction parameter p (see Fig. 10).

Direction fields and stationary points
Our goal is to investigate the dynamical system given by Eq. (14), It is customary to start such an analysis by plotting direction fields in the state space of the system [80]. Examples of the fields for q = 3 and two values of the parameter p are shown in Figure 5. As a remainder: the solution trajectory through a given initial state is a curve in the state space which at every point is tangent to the field at that point. Several things are immediately clear from the picture shown in the figure. At p = 0.2 the flows in the state plane suggest that there are nine stationary points (already marked with red dots). Some of these points are easy to classify. For instance, there are two attractors (i.e. points toward which the system tends to evolve for a wide variety of initial conditions) at (0, 1) and (1, 0) corresponding to a polarized state of the system, i.e. all agents in one clique are in state +1 and in the other -in state −1. Moreover, there are two other symmetric attractors close to the coordinates (0, 0) and (1, 1). It seems that the (almost) complete consensus is possible in our system for some initial configurations, at least for that particular value of p. The point (0.5, 0.5) is a repeller (the system tends to evolve away from it) and the remaining four points seem to be hyperbolic (near such points the orbits of a two-dimensional, non-dissipative system resemble hyperbolas).
At p = 0.4 (right plot in Figure 5) all hyperbolic points and the symmetric attractors disappear. The point (0.5, 0.5) becomes hiperbolic. The only remaining attractors are (0, 1) and (1, 0). Hence for higher values of p the polarization of the system is the only possible asymptotic state. It should be noted that this findings recapture the results from [55].
To find the exact coordinates of the system we just set x and y equal to zero in Eq. (15) and solve the resulting set of equations with respect to x and y, For p = 0.2 and q = 3 we get: The linear stability analysis of these points reveal that indeed P 1 , P 2 , C 1 and C 2 are stable equilibria, R is a repeller and the remaining points are hyperbolic ones, in agreement with our analysis of the direction field in Figure 5. Similarly, for p = 0.4 and q = 3 we have: As before, P 1 and P 2 are stable, but R 1 is now hyperbolic. The remaining points disappeared (they became complex). If we repeat the above calculations for other values of the parameter p and plot the results, we get a bifurcation diagram showing how the dynamics of the system changes with p, i.e. with increasing degree of anticonformity in the system. The plot of the x coordinates of the fixed points as functions of p is shown in Figure 6. The picture for the y coordinates would look the same (but with rearranged labels of the unstable points U i ). Figure 6. Bifurcation diagram of the system given by Eq. (15) for q = 3. Solid lines indicate stable fixed points, the remaining ones -repellers and hyperbolic equilibria. At small values of p the system has four attractors: P 1 and P 2 correspond to asymptotic polarization, C 1 and C 2 -to consensus. There is also a repeller R 1 and four hyperbolic points U 1 -U 4 . As p increases a bifurcation happens at p * 0.26. The four hyperbolic non-symmetric points U 1 -U 4 disappear and the consensus equilibria C 1 and C 2 become hyperbolic. With further increase of p the symmetric points eventually disappear and the repeller R 1 becomes hyperbolic.
Stable equilibria are indicated with solid lines. We see that the system has two attractors P 1 = (0, 1) and P 2 = (1, 0) as well as the unstable point R 1 = (0.5, 0.5) as stationary solutions independently of p. However, the system undergoes two phase transitions. At p 1 0.26 the hyperbolic points U 1 -U 4 disappear and the nontrivial symmetric fixed points C 1 and C 2 become hyperbolic. At p 2 0.32 these symmetric points disappear as well.
The above results were obtained numerically. However, in the special case of q = 2 one can find both critical values of p analytically. For a general value of q an analytical computation of the second critical point is possible as well. We present corresponding calculations in Appendices B and C.

Time evolution of the asymptotic system
The time evolution of our dynamical system for two different values of p is shown in Figure 7. The set (15) of ordinary differential equations was solved numerically in Python by making use of the odeint function from the SciPy package [79].
As already known from Figs. 5 and 6, at p = 0.2 our system may end up either in the polarized state (i.e. P 1 or P 2 ) or in the consensus one (C 1 or C 2 ) depending on the initial conditions for the concentrations of +1 spinsons. We see that the results shown in the left column of Figure 7 are in line with these findings. If for instance the starting point is the total positive consensus in clique A (x 0 = 1.00) and almost total consensus in clique B (y 0 = 0.99), then the system ends up in state C 2 = (0.981, 0.981) representing consensus in the entire system (top left plot in Figure 7). In this case the anticonformistic links between cliques, the number of which is represented by p, lead only to a tiny decrease of the initial concentrations of +1 agents in both cliques. If the initial conditions are close to the repeller R 1 and symmetric (i.e. we start from a point on the diagonal in the state plane, which means that there is (almost) no ordering in each clique), then again the system reaches the consensus state (middle left plot). However, a small deviation from the diagonal pushes the system towards the polarized state (bottom left).
As we see, the behaviour for p = 0.4 is different. The system usually ends up in the polarized state (top and middle right plots in Figure 7). The only exception are initial conditions along the diagonal in the state plane (x, y). Since R 1 is now not a repeller, but a hyperbolic fixed point, the system is pushed towards it in this case (bottom right plot). Again, this is in agreement with the direction field shown in Figure 5.

Basins of attraction
Attractors of every dynamical system are surrounded by a basin of attraction representing the set of initial conditions in the state space whose orbits approach the attractor as time goes to infinity. In the previous paragraph we have seen already examples of initial conditions belonging to the basins of the polarization equilibria (P 1 or P 2 ) and the consensus ones (C 1 or C 2 ) (see Figure 7). Now, we would like to quantify the shapes of the basins of different attractors of our system. Of particular interest in a model with segmentation and negative ties are the basins of the consensus states, since in such a setting one intuitively expects polarization as the natural asymptotic state.
Since it is not possible to calculate the shapes analytically, we will resort to numerical methods again. For that purpose we create first a grid in the state plane (x, y) with both coordinates varying from 0.0 to 1.0 with step 0.01. The points on the grid represent different initial conditions uniformly distributed in the whole state plane. Then we solve the dynamical system (15) for each grid point and check what attractor the solution is converging to at long evolution times. Results for p = 0.2 and p = 0.4 are presented in Figure 8. As often in this paper, the parameter q was set to 3 in the whole procedure.  17) and (18). Boundaries between the basins are formed by stable manifolds of the hyperbolic points U 1 -U 4 (left plot) or of the point R 1 (right plot).
For p = 0.2 (left plot in Figure 8) the whole state plane is divided into the basins of four attractors: P 1 , P 2 , C 1 and C 2 . Their coordinates are defined in Eq. (17). Surprisingly the basins of both positive (C 2 ) and negative (C 1 ) consensus are relatively large. Thus, if the number of negative connections between the cliques is small, the consensus may still by reached in a double-clique network for a range of initial conditions. It is worth to mention that the boundaries between the basins observed in the plot are formed by the stable manifolds of non-symmetric hyperbolic points U 1 -U 4 .
The picture at p = 0.4 (right plot in Figure 8) is simpler. Consensus is no longer possible and the whole state plane is divided into the basins of two polarization points P 1 and P 2 (see Eq. (18) for their coordinates). The boundary between them corresponds to the stable manifold of R 1 , which is now hyperbolic.

Correlation between cliques
All results presented up to this point indicate that some sort of a competition between conformity and anticonformity may be responsible for substantial changes in the dynamics of the model. To elaborate on that issue we will look at the product of the final states of the cliques, as a function of p at different values of q. We will call this quantity a correlation between the cliques. To allow for comparisons with the results presented in [55] we will focus our attention on the total positive consensus x 0 , y 0 = 1.0 as the initial condition. It is referred to as the Scenario I in [55] and corresponds to the following situation: two cliques with a natural tendency to disagree with each other evolve at first independently. They get in touch by chance and establish some cross-links to the other group once they both reached consensus on a given issue. Results for the correlation between the cliques are presented in Figure 9. For values of p smaller than a critical value ( 0.267 for q = 3) both cliques always end up in positive consensus. In other words, in this regime the intra-clique conformity wins with the inter-clique anticonformity and both communities are able to maintain their initial positive consensus. If the value of p is larger than the critical one, the anticonformity induced effects take over and the whole system ends up in a polarized state. Moreover, the critical point shifts with increasing q towards higher values of p (see Table 1 for more details) . Thus, the bigger the influence group, the more cross-links are needed to polarize the society. Let us compare the above results with the model presented in [55]. For that purpose we have to convert the concentrations of agents in state +1 in each clique into average opinions according to Eq. (5) and to transform the probability p of being connected to other clique into the number of cross-links L. We obtain the transformation formula immediately from Eq. (2): The comparisons for two different values of q are shown in Figure 10. The agreement of the results is quite good, despite the differences between the models and the fact that the numerical solution for the annealed model was obtained for an infinite system, whereas the simulations were performed for only N = 100 agents in each clique. Most notably, the transition from consensus to polartization sets in at similar values of L in both models.
Some of the discrepancies between the numerical results for the annealed model and the simulation results for the quenched disorder one are due to the finite system size of the quenched disorder model and the stochastic nature of its simulation. Indeed, the agreement between models is even better, if we leave the numerical solution out of consideration for a moment and concentrate only on the simulation results in both cases. Moreover, as it follows from Figure 11, the differences between the models decrease with increasing system size. Thus, for N → ∞ both variants of the model will probably converge to each other.  Figure 10 for explanation of the labels. The differences between simulation results for both models decrease with increasing system size.

Discussion
Due to the computational complexity of Monte Carlo simulation approach to agent-based models we were able to investigate only few distinct initial conditions in [55]. In this paper, thank to a slight modification of the model, we could analyse it mathematically and therefore explore the whole state plane.
Again, our results indicate that the interplay between conformity and anticonformity may lead to a polarized state of the system. We now have however much better understanding of the conditions necessary to arrive at consensus and we have determined regimes in which polarization takes over. Thus the present results complete the analysis started in [55]. The most important message of the study is that consensus between two antagonistic communities is possible only if they are loosely connected with each other.
It should be noted that similar results have been achieved by Shin and Lorenz [42] within the information accumulation system model of continuous opinion dynamics. The authors found that the convergence of two internally highly connected communities with a comparably low number of cross-links to the same opinion is less possible the more overall interaction between agents takes place.
Explosive growth in Internet mediated communication facilitates the exchange of opinions between people, both passively and actively [81]. Within the language of our model it means that modern communication tools like Facebook and Twitter increase the number of links between people in general and between members of different minded cliques in particular. Our results imply that if the fraction of cross links (i.e. the value of p) between such cliques exceeds some critical value, then the polarized state is the only attractor of the system. In real social systems the situation is of course more complicated, because the cliques interact not only with each other, but with other actors as well. However, our results may indicate one of the possible mechanisms for the omnipresence of polarization in social media. Paradoxically, the often criticized "filter bubbles" on Facebook or Google [82], which separate users from information (people) that disagrees with their viewpoints may help to weaken the problem with bi-polarization and to maintain overall consensus, because they reduce the fraction of cross links between cliques (the discussion of negative effects of "filter bubbles" on the society is beyond the scope of this work).
As far as potential extensions to our model are concerned, the suggestions given in our previous paper still hold. For instance, it could be very informative to check how robust the model is to the introduction of noise, because it is already known that including noise to models of opinion dynamics may significantly change their predictions [83]. Another interesting aspect worth to address in future studies is the casual connection between the network segmentation and the polarization. From the analysis presented in the main text (see Sec. 3.1) it follows that the system undergoes two phase transitions. At the first critical point four hyperbolic nonsymmetric fixed points disappear and the nontrivial symmetric equilibria become hyperbolic. At the second critical point these symmetric points disappear as well.
Our goal is to find both critical values of p analytically in the special case q = 2. We will proceed by first computing the coordinates of the nontrivial symmetric fixed points as functions of p. Then we find coordinates of a symmetric point, at which the largest eigenvalue of the Jacobian is equal to zero. Combining these two results will give us the critical value of the first transition. The second phase transition occurs at a value of p, for which both symmetric points disappear (i.e. they become complex).
For computing the nontrivial symmetric fixed points we put x = y and q = 2 into Eq. (15). Setting x to zero yields Factoring out and simplifying the above equation gives which may be written as and after some calculations we arrive at Solving the last equation gives Note that for p > 1 4 both solutions are complex (i.e. the symmetric points disappear). Hence the second transition occurs at We proceed by computing the Jacobian at a symmetric point (x, x). Due to the symmetry of the dynamical system (15) the Jacobian at such a point is of the form Hence the eigenvalues are (we omit the explicit dependency on x for the sake of simplicity) After some straightforward but lengthy computation we get Since b is negative for p < 1 2 , the largest eigenvalue is always a − b. Setting it to zero yields Combining the last expression with Eq. (A5) gives which is equivalent to We have to solve the above equation in order to get the critical value for the first transition. The relevant solution is p 1 = The critical value of p for the second phase transition may be calculated analytically in case of a general q. At that critical point the equilibrium R 1 = (1/2, 1/2) changes its character from a repeller to a hyperbolic one. Hence, we have to evaluate the Jacobian at R 1 and look for a value of p, for which the smallest eigenvalue becomes zero.
The Jacobian will again have the form given by Eq. (A9). Taking the derivatives of the right-hand side of Eq. (15) with respect to x and y and setting x = y = 1/2 in the resulting expressions give The smallest eigenvalue of the Jacobi matrix is given by Equating it to zero yields the critical value of p, For q = 2 we obtain which agrees with the result obtained in the previous section.

Appendix C. Supplementary materials
A Jupyter notebook with Python code used for the numerical analysis of the dynamical system (15) may be found in [84].