Discontinuous phase transitions in the q-voter model with generalized anticonformity on random graphs

We study the binary q-voter model with generalized anticonformity on random Erdős–Rényi graphs. In such a model, two types of social responses, conformity and anticonformity, occur with complementary probabilities and the size of the source of influence \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_c$$\end{document}qc in case of conformity is independent from the size of the source of influence \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_a$$\end{document}qa in case of anticonformity. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_c=q_a=q$$\end{document}qc=qa=q the model reduces to the original q-voter model with anticonformity. Previously, such a generalized model was studied only on the complete graph, which corresponds to the mean-field approach. It was shown that it can display discontinuous phase transitions for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_c \ge q_a + \Delta q$$\end{document}qc≥qa+Δq, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta q=4$$\end{document}Δq=4 for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_a \le 3$$\end{document}qa≤3 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta q=3$$\end{document}Δq=3 for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_a>3$$\end{document}qa>3. In this paper, we pose the question if discontinuous phase transitions survive on random graphs with an average node degree \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle k\rangle \le 150$$\end{document}⟨k⟩≤150 observed empirically in social networks. Using the pair approximation, as well as Monte Carlo simulations, we show that discontinuous phase transitions indeed can survive, even for relatively small values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle k\rangle$$\end{document}⟨k⟩. Moreover, we show that for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_a < q_c - 1$$\end{document}qa<qc-1 pair approximation results overlap the Monte Carlo ones. On the other hand, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_a \ge q_c - 1$$\end{document}qa≥qc-1 pair approximation gives qualitatively wrong results indicating discontinuous phase transitions neither observed in the simulations nor within the mean-field approach. Finally, we report an intriguing result showing that the difference between the spinodals obtained within the pair approximation and the mean-field approach follows a power law with respect to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle k\rangle$$\end{document}⟨k⟩, as long as the pair approximation indicates correctly the type of the phase transition.

For researchers outside of statistical physics, it may seem strange why do physicists try so hard to distinguish between different types of phase transitions and why do they search so eagerly for discontinuous ones in the models of opinion dynamics [1][2][3][4][5][6][7][8][9][10] . In the context of social systems, such efforts can be justified by the empirically confirmed existence of phenomena strictly related to discontinuous phase transitions, such as tipping points, critical mass, and hysteresis [11][12][13][14][15] .
Several factors may facilitate discontinuous phase transitions, such as: (a) increasing the number of states that the system's entities (agents, particles, etc.) can take 1,16-18 , (b) introducing an aversion to change, which can be related to the size of the group needed to influence the agent 19,20 or another kind of inertia 3,5 , (c) higher dimensionality of the system 21 , higher system degree 5 , or the larger number of layers in the multiplex 22 . Finally, the presence of the annealed instead of the quenched disorder usually promotes discontinuous phase transitions [23][24][25][26] .
Until recently, it seemed that in models of opinion dynamics, a type of nonconformity could also be a factor influencing the type of the phase transition. For example, in the case of the q-voter model 27 , the presence of independence (noise) was needed to induce discontinuous phase transitions and the anticonformity itself was not sufficient 26,28,29 . Analogous observation was done for the symmetrical threshold model 30 . However, two years ago it was shown that discontinuous phase transitions can occur also in the q-voter model with anticonformity if we assume that the size of the group of influence needed for conformity q c is independent of the size of the group of influence needed for anticonformity q a 7 . Such a generalized model, which reduces to the basic q-voter model with anticonformity for q c = q a = q , was studied so far only on the complete graph. It was shown that in the case of the complete graph, which corresponds to the mean-field approach (MFA), it displays discontinuous phase transitions for q c ≥ q a + q , where q = 4 www.nature.com/scientificreports/ for q a ≤ 3 and q = 3 for q a > 3 . However, real social networks, especially the friendship ones, do not have a complete graph's structure. It was predicted within the social brain hypothesis and then confirmed empirically that humans have an average of about 150 relationships and this is true for both offline and online (Facebook, Twitter) networks [31][32][33] . It means that the average degree k of a node in the social network is substantially smaller than the size N of the entire social network, whereas for the complete graph �k� = N − 1 . On the other hand, the higher network degree facilitates discontinuous phase transitions 5 . Therefore, the natural question that arises here is if the model with generalized anticonformity can display discontinuous phase transitions for �k� ≪ N.
To answer the question, we study the model on random Erdős-Rényi (ER) graphs with an average node degree k that corresponds to those that were found empirically in social networks 32,33 . We decided on such a simple structure because in such a case the pair approximation (PA) should be more valid than for networks with a higher clustering coefficient 26,[34][35][36][37] . However, to validate the analytical results, we conduct also Monte Carlo simulations, because it was reported recently that PA can give invalid results even on ER graphs if the mean degree of nodes is small and comparable with the size of the influence group q 10,29,38 . In such a case, the results can be wrong not only quantitatively, but even qualitatively, indicating discontinuous phase transitions, which are not observed in the computer simulations.

Model
We consider a system of N interconnected voters (also called agents or spins). Each of them is characterized by a single dynamical binary variable s(x, t) = j , where j = +1 (↑) or j = −1 (↓) , x = 1 , . . . , N, and t denotes time. From the social point of view, s(x, t) represents an opinion of an agent placed at node x at time t on a given topic measured in the two-point psychometric scale (yes/no, agree/disagree). We use j to describe the state of a single spin for consistency with earlier papers in which PA was used for the q-voter model 26,39 .
Agents are placed in the vertices of an arbitrary graph and can interact only with those agents that are directly linked to them. Following 40 we use here the term source for the group of influence and the term target for the voter, who is influenced (the active one). We are taking into account two types of social response, conformity and anticonformity, and we use the annealed (situation) approach 41 . It means that each agent can behave in two distinct ways with complementary probabilities: it can copy the opinion of the source or take the opposite one.
The social influence in this model is imposed only by the unanimous group of q agents that are randomly chosen (without repetitions) out of k x nearest neighbors of a target at node x. Theoretically, it may happen that a given voter at node x cannot be influenced at all, because its degree k x < q . However, in social reality and our study, such a situation appears rarely. The typical size of the influence group q varies from a few to a dozen people 42 . On the other hand, the real egocentric personal social networks consist on average of 150 relationships of successively inclusive layers at 5, 15, 50 and 150 alters 32,33 . Therefore, here we will also consider �k� ≤ 150 , with the condition k > q.
In the original model, the source of influence consists of q agents for both types of social response 19,29,41 . However, recently it was proposed that the size of the source needed for conformity q c may not be the same as the size q a needed for anticonformity 7 . Such a generalization led to discontinuous phase transitions, which were not observed before within the q-voter model with anticonformity, neither on the complete graph 28,41 , nor on the random graph 29 .
The generalized model was studied so far only via the mean-field approach, which corresponds to the complete graph 7 . Here, to answer the question posed in the Introduction, we study the model on random graphs with �k� ≤ 150 . As in all previous studies on the q-voter model, we use the random sequential updating (RSU) scheme, which is used to mimic the continuous time, and the algorithm of an elementary update is the following (see also Fig. 1): 1. a single vertex is chosen randomly x ∼ U{1, N}, 2. a random number from the uniform distribution r ∼ U(0, 1) is drawn to decide about the type of response: (a) if r < p , where p denotes probability of anticonformity, then a source of q a agents is randomly drawn (without repetitions) out of k x target's direct neighbors; if the source is unanimous, the target takes the opposite state to the one of the source, (b) otherwise, a source of q c agents is randomly drawn (without repetitions) out of k x target's direct neighbors; if the source is unanimous, the target takes the same state as the source.
As usually, a unit time, i.e., Monte Carlo step (MCS), consists of N elementary updates.

Pair approximation
In the previous paper, we considered the model on the complete graph, for which the MFA gives exact results 7 .
Here we consider the model on random ER graphs with the average node degree �k� < N − 1 and thus we use another analytical approach, so-called pair approximation 26,34,43 .
Within MFA, we assume that the global concentration of agents with positive opinion, defined as approximates correctly the local one. Therefore, the probability to choose randomly an up-spin at time t in the neighborhood of node x is just equivalent to c(t). However, this assumption is strictly valid only for the complete graph. A step further is to distinguish between the probability of choosing randomly an up-spin in the www.nature.com/scientificreports/ neighborhood of an up-spin and in the neighborhood of a down-spin. To make this distinction, the concentration of active bonds b(t) is introduced, with an active bond defined as a bond between opposite spins. As shown in 26,44 the set of equations that describes the evolution of the system within PA is universal for all single-spin flip binary dynamics: where k is a degree of a given node, j is the state of this node, i ∈ {0 , 1, . . . , k} is the number of active edges attached to this node, P(k) is the degree mass probability function, θ j is the conditional probability of selecting a node that is in the opposite state to its neighbor in a state j, and c = −j and � b = 2(k − 2i)/�k� are rescaled elementary changes in corresponding quantities per one Monte Carlo step. Moreover, for all single-flip binary dynamics θ ↑ = b/2c and θ ↓ = b/2(1 − c) 26,44 . Only the function f (. . .) depends on a given model and its parameters, since it describes the flipping probability of a node in a given state. Therefore, the only task in a broad class of binary opinion dynamics is to calculate f (. . .) 44 .
In our case, an opinion change happens only if q c neighbors with opposite opinions are drawn in case of conformity or q a neighbors with matching opinions are selected in case of anticonformity. Therefore, as the total number of neighbors with opposite opinions is equal to the number of active links i of a given voter and the group of influence is chosen without repetition, the flipping probability takes the following form www.nature.com/scientificreports/ Thus, in the limit of an infinite graph, the time evolution of this model is defined by a system of differential equations: As always for the q-voter model, evolution equations (4) does not depend on the degree distribution P(k) but only on the average degree k 26,29 . Hence, from now on we denote k by k for brevity. Unfortunately, for q a = q c we were not able to obtain the analytical solution of (4), even for the stationary state. Therefore, we solved it numerically.

Monte Carlo results
To compare PA with MC results, we conducted simulations on sufficiently large graphs N = 5 · 10 5 . We used ensemble averaging, which means that for each set of parameters ( q a , q c , p, k, N) we constructed M independent graphs. It occurred that for a graph of size N = 5 · 10 5 the number of samples M = 10 was sufficient. As usual, the most problematic simulation parameter to choose was the "thermalization" time τ needed to reach the stationary state. It is well known that it increases dramatically near the critical point and therefore initially we used, similarly as for N and M, different values of τ to check which one would be sufficient. Finally, we decided to use τ = 5 · 10 5 . The main goal of this paper is to check if the discontinuous phase transitions appear for realistic values of k. Hence, having in mind the condition q c , q a < k and realistic values of q and k found in the empirical studies 32,33,42 , we use k = 50 , 150 to show the dependence between the stationary value of the up-spins' concentration and the probability of anticonformity p. Because we want to check the type of the phase transition, we start from two different types of initial conditions, ordered c 0 = c(0) = 1 and disordered c 0 = c(0) = 0.5 ones. Such an approach allows us to see the hysteresis, which appearance is one of the indicators of a discontinuous phase transition.
As shown in Fig. 2, starting from the ordered initial state c 0 = 1 the system remains in the ordered state for p < p * upper , i.e. below the upper spinodal, and it reaches a disordered state for p > p * upper , i.e. above the upper spinodal. On the other hand, for c 0 = 0.5 the system remains in the disordered state for p > p * lower , whereas for p < p ( www.nature.com/scientificreports/ The fact that p * lower � = p * upper , denotes the existence of the hysteresis: for p ∈ (p * lower , p * upper ) the stationary state depends on the initial one. Additionally, the jump of c st , which can be easily converted to the order parameter φ = 2c st − 1 , is seen. Both signatures of discontinuous phase transition, the hysteresis and the jump of the order parameter, are visible in PA and MC results. Moreover, as shown in Fig. 2, MC overlaps PA results almost perfectly and approaches MFA with the increasing k, as expected.
However, it occurs that not for all values of parameters q a , q c the agreement between PA and MC results is equally good, which is clearly seen in Fig. 3. Within PA the hysteresis is seen, i.e., p * lower < p * upper , whereas within MC simulation no hysteresis is observed, i.e. p * lower = p * upper . As previously, PA results approach MFA ones with the increasing k, and the hysteresis vanishes. To evaluate the values of the parameters for which PA gives correct predictions, we measure the width of the hysteresis, defined here as the distance between the upper and the lower spinodal p * upper − p * lower , as a function of q a for the fixed value of q c within all three methods: PA, MFA and MC. In Fig. 4 we see results for q c = 10 but we have checked also other values of q c , which allows to draw general conclusions.
Let us recall first the results obtained analytically within MFA 7 : a hysteresis indicating a discontinuous phase transition has been found for q a ≤ q c − 3 . Indeed, we see in Fig. 4 that for q c = 10 and q a ≤ 7 , the hysteresis appears within MFA, PA, and MC. Moreover, for q a < q c − 1 MC overlap PA results if only k is sufficiently large (as in many other studies PA fails if the size of the influence group approaches k 9,10,26,29 ).
For q a ≥ q c − 1 the hysteresis starts to increase within PA, then it reaches the maximum at a certain value q * a = q * a (q c ) , for example q * a (10) = 13 , and then it decreases. According to MFA, no hysteresis should be seen for q a > q c − 3 and indeed no hysteresis is seen for these values of parameters within MC. Even at q * a = q * a (q c ) , in which PA predicts the maximum width of the hysteresis, it does not occur in Monte Carlo simulations, which is clearly seen in Fig. 5.
In all presented results on Figs. 2, 3, 4, 5, it is seen that the difference between PA and MFA results decreases with increasing k, which is an expected result. However, it is not that obvious how it decreases. Therefore, we decided to check how the spinodals p * lower , p * upper obtained within PA deviates from the mean-field's ones. Results for q c = 10 are shown in Fig. 6. First of all, we see that for both spinodals the power-law decay with k is observed in case of q a < 9 (panels in two upper rows in Fig. 6). We checked that the same phenomenon appears for other values of q c . As long as q a < q c − 1 PA predicts correctly the type of phase transition and simultaneously the power-law decay is observed. The small deviation from the power-law is seen only for small values of k, for which PA also fails. Obviously, we cannot prove that the power-law holds for arbitrarily large k, but it exists for all examined values up till k = 109350 , which is a huge number from the social point of view.
For q a ≥ q c − 1 PA predicts discontinuous phase transitions, which are neither predicted by the MFA nor seen in the simulations, and simultaneously the power-law breaks down. It is probably worth to recall that the condition q a = q c − 1 corresponds to the situation in which the same number of unanimous agents have to be selected to change the state of the system for both types of social response: www.nature.com/scientificreports/ • in case of anticonformity, the state of the system changes if a voter and the source of size q a are in the same state, so in total q a + 1 agents need to be unanimous, • in case of conformity, the state of the system changes if a voter has an opposite state to the source of size q c , so in total q c agents need to be unanimous.  www.nature.com/scientificreports/ Unfortunately, the above remark neither explains why PA fails nor why the power-law breaks down. We decided to report here these intriguing results although unfortunately we are not able to explain them. We never saw this kind of analysis before, so we are not able to compare our results with others to gain some additional intuition. However, we hope that maybe some readers will be able to explain this intriguing phenomenon.

Summary
In this paper, we investigated the q-voter model with generalized anticonformity on random graphs via Monte Carlo simulations, as well as pair approximation. It occurs that, similarly as within MFA, discontinuous phase transitions appear only if the size of the influence group q c needed for conformity is sufficiently larger than the size of the influence group q a needed for anticonformity, precisely for q a ≤ q c − 3 . For these values of parameters PA results overlap MC ones. Moreover, PA gives results in agreement with Monte Carlo simulations in the case of q a < q c − 1 for any value of q a . It predicts properly both the type of phase transition, as well as the values of the spinodals, which are overestimated by MFA for k ≪ N . However, for q a ≥ q c − 1 PA wrongly predicts the type of phase transition and shows an unexpected non-monotonic behavior of the size of hysteresis. As expected, PA results approach MFA ones with the increasing average node degree k. However, the difference between the results obtained within PA and MFA decreases with k as power-law. Results are presented in the log-log scale (the natural logarithm is used here). For the lower spinodal, the power-law appears in the whole examined range of k. On the other hand, the upper spinodal deviates from the power-law for q a ≥ q c − 1 for small values of k. It should be noticed that simultaneously for these values of parameters PA wrongly predicts discontinuous phase transition, i.e. p * upper � = p * lower , whereas MFA indicates the continuous one, i.e. p * upper = p * lower .