Opinion Dynamics Systems on Barabási–Albert Networks: Biswas–Chatterjee–Sen Model

A discrete version of opinion dynamics systems, based on the Biswas–Chatterjee–Sen (BChS) model, has been studied on Barabási–Albert networks (BANs). In this model, depending on a pre-defined noise parameter, the mutual affinities can assign either positive or negative values. By employing extensive computer simulations with Monte Carlo algorithms, allied with finite-size scaling hypothesis, second-order phase transitions have been observed. The corresponding critical noise and the usual ratios of the critical exponents have been computed, in the thermodynamic limit, as a function of the average connectivity. The effective dimension of the system, defined through a hyper-scaling relation, is close to one, and it turns out to be connectivity-independent. The results also indicate that the discrete BChS model has a similar behavior on directed Barabási–Albert networks (DBANs), as well as on Erdös–Rènyi random graphs (ERRGs) and directed ERRGs random graphs (DERRGs). However, unlike the model on ERRGs and DERRGs, which has the same critical behavior for the average connectivity going to infinity, the model on BANs is in a different universality class to its DBANs counterpart in the whole range of the studied connectivities.


Introduction
The interest in sociophysics has greatly increased in the past two decades, mainly when considering the dynamics that are present in social systems or networks [1][2][3][4][5][6][7][8]. Stauffer [5] and Galam [9], who proposed models that use local majority rule arguments, are considered the predecessors of the sociophysics emerging field, also termed as "the dynamics of opinions". In fact, the dynamics of opinions are treated in the same way that researchers treat the usual real world [10]. As a result, opinion-dynamics systems have turned out to be one of the most studied subjects in sociophysics [7,[11][12][13][14].
Biswas, Chatterjee, and Sen [8] proposed a continuous opinion dynamics model that is nowadays called the BChS model. This non-equilibrium system has pair interactions that can be either positive or negative, modeled by a single noise parameter q, that represents the fraction of negative interactions. The continuous version, on fully connected networks, has been studied through numerical simulations. It has been shown that a continuous phase transition takes place at a critical value q = q c with mean-field critical exponents [8]. On the other hand, this same model has been treated on regular lattices in two and three dimensions by using Monte Carlo (MC) simulations [15]. Although a continuous phase transition is also obtained, the criticality is now the same as the Ising model in the same dimensions.
The BChS model can also be defined on a discrete version and it has been studied by using Monte Carlo simulations on several scale-free networks and random graphs. The discrete model undergoes a second-order phase transition and, from the results that have been reported in the literature, we can say that: (i) its critical behavior is the same as the continuous version on fully connected networks and on regular two-and three-dimensional lattices [8,15]; (ii) on Apollonian networks (ANs), the critical behavior is different from the Ising model on regular lattices [16]; (iii) on directed Barabási-Albert networks (DBANs), the critical exponent ratios change with the connectivity of the networks [17], but the model belongs to the same universality class as the majority vote model (MVM) [18] on the same DBANs; (iv) on Erdös-Rènyi random graphs (ERRGs) and directed ERRGs random graphs (DERRGs) the critical exponents are different from the Ising model and also the MVM model on the same random graphs [19]. For additional details on the related MVM see, for instance, Refs. [20][21][22][23] It is then well understood that it is not possible to know, a priori, whether a given model, on a particular scale-free network, will present a phase transition and, if a secondorder phase transition occurs, what will be its universality class. For this reason, it will then be quite interesting to understand what should be the behavior of the discrete BChS model, when defined on BANs and, at the same time, compare the new behavior with the previous results of the model on other networks (and also with the MVM). Moreover, except for the Voronoi-Delaunay random lattices, the BANs seem to be the lacking most common scale-free topology to be studied within the context of the discrete BChS model [24]. In fact, as we will see below, there are some quite different and unexpected critical behaviors of the discrete BChS model on BANs, when compared to the DBANs and also the same model on other random network topologies.
Thus, in this work, the social consensus formation in the non-equilibrium discrete version of the Biswas-Chatterjee-Sen dynamic system on BANs, for several values of the noise parameter q, has been studied through Monte Carlo simulations. The scope of the paper is as follows. The next section presents the model, the MC simulations that have been employed, together with the evolution of the physical quantities from which the transition has been characterized. In Section 3, the obtained results are discussed and, in the last section, conclusions and some final remarks are addressed.

Model and Simulation
In order to study the BChS model [8,16,17,19] defined on Barabási-Albert networks, we will closely follow the procedure outlined in Ref. [19]. However, for questions of completeness, the model and simulations will be shortly summarized below.
Agents (or individuals) are set on each i node of a BAN with N sites. At time step t, the agents have opinion variables o i (t), which can assume three different values say −1, 0, or +1 . The rules for updating o i (t) according to the BChS model are as follows: (i) The initial configuration is constructed by randomly assigning one of the three opinion states for each site i of the BAN; (ii) A site i is then randomly select to be updated; (iii) One bound of the site i is also randomly selected and an affinity µ ij is given for this bond (j is the corresponding site sharing the bond with site i). This affinity parameter is another discrete variable that assumes a value +1, but can be turned negative with a probability q. The parameter q acts, in fact, as an external noise, modeling local discordances; (iv) The opinion variable of both sites sharing the selected bond are now updated following the rules where o i (t) and o j (t) are the opinion states at time t, while o i (t + 1) and o j (t + 1) are the updated opinion states of the two sites i and j, respectively; (v) When the opinion state is out of the interval [−1, +1], for example being larger than +1, it is automatically made equal to +1. The same happens when the opinion state is smaller than −1, when it is made equal to −1.
An order parameter O for this dynamical system can be defined by averaging the opinion variables o i (t) over all individuals, i.e., where t is large enough for the system reaching the stationary state. In the thermodynamic limit (infinite size networks), there is a critical value q = q c , where for q < q c an ordered phase with O = 0 is present, while for q > q c one has a disordered phase with O = 0 instead. Exactly at q = q c , both ordered and disordered phases become equal at a secondorder phase transition. It is interesting that, in this case, rather than a thermal driven critical behavior as in usual magnetic systems, one has indeed a kind of a random configuration driven phase transition. This dynamic order-disorder transition can be studied, via MC simulations, using similar magnetic-like variables, so common at second-order phase transitions. For instance, for a given value of q, Equation (3) gives us a natural order parameter O(q). The fluctuation of the order parameter, O f (q), is in this way the analogous to the magnetic susceptibility. We can even define the equivalent of the reduced Binder cumulant of the order parameter, namely O 4 (q). These quantities are easily computed from MC simulations as where · · · t stands for time averages, that are computed after the system has reached the stationary state, and · · · av means the averages over different initial configurations. Some details of the simulations are described below. In this work, the above quantities, as a function of q, have been computed through extensive MC simulations on BANs of finite sizes ranging from N = 250, 500, 1000, 2000, 4000, 8000, up to 16, 000. In order to let the system reach its stationary state, the initial 10 5 MC steps (MCS) have been discarded. The corresponding time averages have then be computed by taking the next 2 × 10 5 MCS. Here, one MCS consists of randomly choosing N sites of the network, as described above. For each set of network size N and parameter q, 10 3 to 10 4 different configurations have been considered to obtain the configurational averages. Now, close to the critical noise parameter q = q c , and for large system sizes, quantities (4), (5) and (6) should have a power-law behavior of the form (for further details see, for instance, Ref. [25]) where c is a non-universal constant and q c the corresponding critical noise for an infinite size graph. The critical exponents β, ν, and γ come from the behavior of the order parameter, correlation length, and fluctuation of the order parameter, respectively, when taken as a function of the noise parameter q close to q c . In the above equations, however, we have the size behavior instead, where f k (x), with x = N 1/ν (q − q c ) and k = {O, O f , O 4 }, are scaling functions. As we shall see below, for each finite graph of N sites, we can estimate its pseudo-critical noise parameter q c (N) by looking, e.q., at the peak presented by the fluctuation of the order parameter as a function of q.
We can now follow the standard procedure to obtain, numerically, the phase transition properties of the model in the thermodynamic limit. For a second-order phase transition, which will be the present case, the reduced Binder cumulant O 4 , in Equation (9), should be independent of the system sizes for |x| << 1 (in fact, for large systems and close to q c ). This means that the crossing points of O 4 , as a function of q and for different values of N, will give a good estimate of the critical noise parameter q c [25]. With q c in hands, and using Equations (7) and (8), one can compute the exponents ratio β/ν and γ/ν, respectively. Now, by locating the value of the noise parameter q c (N), where the maximum of O f occurs as a function of q, we can explore Equation (10) and further obtain an estimate of the critical exponent 1/ν.
Within this scheme, the main critical behavior of the system can be quantitatively specified. In addition, using the hyper-scaling hypothesis [26] 2β/ν + γ/ν = D eff (11) it is also possible to compute the effective dimension D eff .

Results and Discussion
The dependence of the reduced Binder cumulant O 4 on the noise parameter q, for several finite sizes N, is shown in Figure 1a, for the smallest connectivity z = 2, and in Figure 1b, for a quite higher value z = 100. From these figures, one can see that: (i) the cumulant crossings attest indeed that we have a second-order phase transition in the system, where q c = 0.173(1) for z = 2 and q c = 0.248(2) for z = 100; (ii) from these values of q c , we can see that the critical noise parameter is dependent on the connectivity z. Similar results are obtained for different connectivities and one can thus compute q c for several values z. Table 1 gives the results so obtained for other selected connectivities.   (9) 0.917 (12) With q c in hand for different connectivities z, we can now compute the critical exponent ratios from the power-law relations of Equations (7), (8), and (10). The log-log plot of the quantities defined in Equations (4)-(6) and (10), as a function of N, on BANs with different connectivities z, are shown in Figure 2. All quantities have been computed at the extrapolated infinite network critical value q c for the corresponding connectivity z. The alignment of the data in those figures strongly corroborates the transition being indeed of second order, with the slope from a linear fit giving the desired critical exponent ratio. The critical exponent ratios so obtained are given in the legend of Figure 2 and also in Table 1, together with additional values of the connectivity z. at q max c for each BANs size N , can also be used in Eq. (8) to get an additional estimate o the critical exponent ratio γ/ν. Proceeding in this way, we obtain exponent ratios that ar At this point, the universal behavior of the present finite-size-scaling functions, which are given in Equations (7)-(9), can be further tested in a still wider range of q. Note that we have, up to now, been working only with noises in the region close to the transition q c . Figure 3 depicts the collapse of the present data of the scaled order parameter ON β/ν (a), the scaled fluctuation of the order parameter O f N −γ/ν (b), and the scaled reduced Binder cumulant O 4 (c), as a function of the scaled displacement (q − q c )N 1/ν for the particular connectivity z = 20. We can see from the collapse of the data that we are not only dealing with a second-order phase transition, but also the computed critical exponents from the simulations seem to be indeed accurate. The fluctuation of the order parameter O f (q) has a clear peak close to the connectivity transition, as is depicted in Figure 3b for z = 20. The value of this peak O f (q max c ), located at q max c for each BANs size N, can also be used in Equation (8) to obtain an additional estimate of the critical exponent ratio γ/ν. Proceeding in this way, we obtain exponent ratios that are similar to those obtained from the same scaling relation using the previously computed q c , for all values of z. This is clearly seen in Table 1 by comparing the numbers in the fifth and sixth columns, and also by inspecting Figure 2b,d.
Regarding the γ/ν quantity, it is worthwhile to stress that on DBANs [17], for all values of z, a different critical exponent ratio is obtained when considering the corresponding maximum value of the fluctuation of the order parameter O f (q max c ). This can be seen in the fifth and sixth columns on Table 2, that reproduce the previous data of the same model on DBANs from Ref [17].
From Tables 1 and 2, it is also possible to have a numerical comparison between the present results on BANs and those earlier obtained on DBANs, for several connectivities z. However, a global view of the data conveyed in Tables 1 and 2 is better seen in the left panels of Figure 4. They plot the effective dimension D e f f , and the ratios β/ν and γ/ν, for the selected values of the connectivity z, on BANs (full circles) and on DBANs (open circles). From those plots we see that: (i) the effective dimension value D e f f ≈ 1.0 for all z on both networks; (ii) after a slight variation for small values of z, the critical exponents become almost constant for z > 20 on both networks; (iii) the critical exponent ratios on BANs are different from those on DBANs, meaning that they belong to different universal classes.
Although the effective dimension is close to one, the model has indeed a clear phase transition, in contrast to the same model on regular one-dimensional lattice, where no transition is observed. Nevertheless, we have to keep in mind that on these networks there are longer-range interactions, and not only short-range finite interactions that allows us to prove that one-dimensional classical models are free from any phase transition. Table 2. Results of the model on DBANs according to Ref. [17]. The legend of Table 1 also applies in this case, with the exception of the column regarding the exponent 1/ν, which has not been computed in [17]. It is also interesting to compare the critical behavior of the discrete BChS model on BANs and DBANs to the critical behavior of the model previously studied on ERRGs and DERRGs [19] by using similar MC simulations. To make such comparison, we have the same quantities on ERRGs and DERRGs from Ref. [19] in the right panels of Figure 4. One can see that, on ERRGs and DERRGs, after an oscillatory behavior of the critical exponents for small values of z, the exponents are almost the same in both graphs, even in the limit of z → ∞. This is due to the fact that in this limit, ERRGs and DERRGs have the same small world networks behavior, which does not happen to be the case on BANs and DBANs. Finally, from the data of Table 1, one can notice that the value of the critical noise parameter increases with the connectivity, and approaches a limiting value as z → ∞. This behavior is better seen in the phase diagram depicted in Figure 5 for the model on BANs. We also note, from Tables 1 and 2, that q c , for large values of z, is different for both networks. This is due to the fact that, unlike the ERRGs and DERRGs, in the limit z → ∞ the BANs and DBANs do not have the same behavior as the small-world lattice. However, as for the model on ERRGs and DERRGs, the critical transition line does obey a similar inverse power law of the form where A and B are non-universal constants. The inverse power-law of the transition line can be clearly seen in Figure 5b. From a linear fit to the data we obtain the constant values given in Figure 5b, from which we have the critical noise in the z → ∞ limit as q c = 0.2488 (1). It is also useful to have Equation (12), since one can now estimate the critical noise q c for any value of z.

Concluding Remarks
The non-equilibrium BChS model, in its discrete version and defined on BANs, has been studied through Monte Carlo simulations for several values of the local consensus controlling parameter q. The results show that the model undergoes a second-order phase transition with critical exponents in a different universality class as the same model on regular twoand three-dimensional lattices [15], ERRGs and DERRGs [19], as well as on DBANs [17]. Despite of belonging to another universality class, the hyper-scaling [26] relation is always valid, independent of the connectivity z, and provides an effective dimension close to 1.0. It should be stressed that the critical exponents of the discrete BChS model on BANs are not the mean-field ones, as in the continuous version of the model on the same BANs [27].
Unlike the model on DBANs, the critical exponent on BANs, coming from the peak of the fluctuation of the order parameter, is comparable to the exponent directly obtained from the critical noise parameter. In addition, in the z → ∞ limit, the results on BANs are different from those on DBANs, a result of the non-equivalence of both networks for large values of the connectivity. We believe that the present results on the BANs, in some sense complete the study of the discrete BChS model on scale-free networks and random graphs. Table 3 has a summary of the universality class of the BChS model on several networks and regular lattices. Table 3. Summary of the universality class of the discrete BChS model on several networks. All transitions are second order and the class of its own means that we are not aware yet of another member sharing the same critical exponents.