Kinetic Models of Discrete Opinion Dynamics on Directed Barabási–Albert Networks

Kinetic models of discrete opinion dynamics are studied on directed Barabási–Albert networks by using extensive Monte Carlo simulations. A continuous phase transition has been found in this system. The critical values of the noise parameter are obtained for several values of the connectivity of these directed networks. In addition, the ratio of the critical exponents of the order parameter and the corresponding susceptibility to the correlation length have also been computed. It is noticed that the kinetic model and the majority-vote model on these directed Barabási–Albert networks are in the same universality class.


Introduction
There has been, recently, a great interest in the study of networks [1][2][3], which are different from the regular crystalline Bravais lattices, and are also frequently called scale-free networks. The treatment of such complex lattices was originally motivated by social organizations and computer connections, ranging from networks in nature to networks of people as well. However, spin systems defined on such complex networks have also been considered in the purpose to determine the character of its phase transition, if present, and the corresponding universality class in the case of critical behavior (for a recent review see reference [4]).
It should be mentioned that not only spin models have been considered on such networks. In the particular case of the directed Barabási-Albert networks (DBAN), it has been shown that the nearest-neighbor spin-1/2 Ising model has no phase transition [4][5][6]. Nevertheless, on these same networks the non-equilibrium majority-vote model (MVM) [7] presents a well-defined order-disorder dynamical phase transition [8].
Regarding models on networks, in 2012 Biswas, Chatterjee, and Sen [9] introduced a very interesting kinetic model of continuous (or discrete) opinion dynamics. The Biswas-Chatterjee-Sen (BCS) model has mutual interactions that can be both positive and negative, and a single parameter p that represents the fraction of negative interactions. Numerical simulations of the continuous version of the BCS model have indicated the existence of a universal continuous phase transition at a critical value p = p c with mean field exponents ν = 2.00(1) for the correlation length, β = 0.50(1) for the order parameter and γ = 1.00(1) for the order parameter fluctuation (susceptibility). The original BCS model has been defined on a fully connected graph of infinite range and, as said above, it can be treated in its continuum or discrete opinion dynamics versions. However, Mukherjee and Chatterjee [10] extended the model to square and cubic lattices, and the numerical results have indicated that the critical behavior of the BCS model on these lattices is the same as that of the Ising model in the same dimension. Thus, it should be interesting to study the BCS model on the complex DBAN and to compare the results to those obtained by the majority-vote model on the same network [8]. For a better comparison we will treat here the model within the discrete opinion dynamics version.
The plan of the paper is the following. In the next section, the BCS model in the discrete opinion dynamics, with the corresponding Monte Carlo simulations, and the thermodynamic quantities used to obtain the critical behavior are presented. The results are discussed in Section 3 and some conclusions are summarized in the last section.

Biswas-Chatterjee-Sen Model
The discrete version of the BCS model [9] on a DBAN can be described as follows. Consider a set of agents (individuals) which have, at time t, opinion variables o i (t), where i denotes the particular node of the DBAN with the network having a total of N agents. In the discrete version of the model, the opinion of an agent or individual i assumes only values ±1. In addition, the opinions change out of pairwise interactions via mutual influences/couplings µ ij in such a way that we have where t is the discrete time and the pairwise interactions between sites i and j, µ ij , and between sites j and i, µ ji , are quenched and integer variables as well, since they are associated with the links of the network, so they do not change in time. Due to the directed character of network links, the interactions µ ij and µ ji are not symmetrical and, in the most general case, can also assume zero value, stemming from the network construction. It can be seen, in the above dynamics, given by Equations (1) and (2), that an agent i updates his opinion by interacting with agent j, and it is influenced by the mutual coupling term µ ij . As in the Ising model and MVM this discrete version of the BCS model presents the up-down symmetry with o i (t) = ±1. In the present study, µ ij is a discrete variable which takes the value −1 or +1 with probability p or 1 − p, respectively. In other words, the disorder parameter 1 − p denotes the fraction of positive pairwise interactions.
The ordering in this system can be determined by the average opinion over all individuals in the network, defined by By changing the fraction p of negative interactions in the fully connected graph of infinite range, one can observe, for p < p c , a symmetry breaking ordered phase with a non-zero value of the order parameter O. On the other hand, for p > p c , a disordered phase exists with O = 0. At p = p c a continuous transition between the ordered and disordered phases takes place.
To study the critical behavior of the present model in the DBAN we can consider the same Ising-like quantities for the second-order phase transition. However, instead of a thermal driven criticality, here we have a configurational disorder driven transition. Thus, for the present model we are interested in the order parameter as a function of the disorder p, O(p), the order parameter fluctuations χ(p), which is the analogous of the magnetic susceptibility in spin models, and the corresponding reduced fourth-order Binder cumulant U 4 (p), where the last two quantities can be defined as where · · · stands for configurational averages (and sample averages in our simulations, as is discussed below) computed at the steady states. The above-mentioned quantities are functions of the disorder parameter p and, close to the transition point p = p c , they obey the additional finite-size scaling relations where ν, β, and γ are the usual critical exponents of the correlation length, order parameter and susceptibility, respectively, and f O (x) and f χ (x) are the finite-size scaling functions with being the scaling variable. The fourth-order Binder cumulant U 4 should be independent of (large) system sizes for |x| << 1, which allows one to estimate the critical value p c . Therefore, from the size dependence of O and χ at p c we can obtain the exponents ratio β/ν (O) and γ/ν (χ).
It is also possible to evaluate the effective dimensionality, D eff , from the hyperscaling hypothesis

Simulations
We give below some details of the computer simulations. First, we initialize a directed Barabási-Albert network with only z neighbors which are connected with each other. Then, additional nodes are connected to this core following the Barabási-Albert prescription. In the end, we will reach a complete graph with a maximum value of N + z nodes. The links are directed by building a "Kertesz list", where the node j can be in the neighbor list of node i, whereas the node i does not necessarily belong to the neighbor list of node j. After, the Monte Carlo simulations have been performed on DBAN with various values of the connectivity z, namely ranging from z = 2 to z = 100. In the present case, the Monte Carlo method is used in the evolution of the network together with the acceptance of the agents update. However, differently from the MVM that uses a master equation for its time evolution, here one uses the acceptance rules of the BCS model evolving on a Markov chain, which is a usual tool of the Monte Carlo approach.
For each connectivity z we have used system sizes with N = 250, 500, 1000, 2000, 4000, 8000, and 16000 nodes. Although the relaxation time, and the number of Monte Carlo steps (MCS), as well as the number of independent runs to generate the statistics depend on the system size, we have taken, in general, for the larger networks, 1.5 × 10 5 (MCS) for the system to reach its steady state, and then the time averages have been calculated over the next 3 × 10 5 MCS. One MCS here is accomplished after N attempts to update the opinions of agents i and j, considering the evolution Equations (1) and (2). For all sets of parameters p and z, we have generated 100 distinct networks and, for each distinct network, we have simulated 100 independent runs. These quantities have shown to be enough to get data with reasonable error bars (as will become clearer below) and within our available computational time. Figure 1 displays the fourth-order Binder cumulant U 4 as a function of the disorder parameter p for several DBAN with different number of nodes. In (a) the connectivity number is z = 2 and in (b) we have z = 100. It is clear from these figures that the system undergoes a second-order phase transition, since the cumulants tend to cross at the same value at the critical disorder parameter p c [11]. We can also note that the smaller the connectivity is, the more apparent the crossings are. Despite still having finite-size effect in the cumulants, the estimate of p c has been made by taking only the cumulant crossings of the larger nodes lattices (pairs of networks ranging from 2000 to 16,000 nodes).

Results and Discussion
In the examples of Figure 1 we have p c = 0.439(3) for z = 2 and p c = 0.388(4) for z = 100. In the above data the errors are statistical ones. These, and other values of the critical disorder parameter for different connectivities, are given in the first row of Table 1.  Table 1. For clarity, only the general trend of the cumulant behavior is shown without the error bars.
In Figure 2 we have the ln-ln plot of the average opinion at the critical disorder, O(p c ), as a function of the number of nodes N, for several values of z. In this case, it is easy to see that, from Equation (6a), the slope of the linear fit gives the critical exponent ratio β/ν. The corresponding results are reproduced in Table 1 together with the results for additional values of the connectivity z. From Figure 2 we can also notice that the critical exponent ratio β/ν changes as the connectivity z changes, implying in a sort of non-universal behavior. As we will discuss below, this is the general trend for the susceptibility exponent ratio as well. Figure 3 displays the ln-ln plot of the order parameter susceptibility at the critical disorder value, χ(p c ), as a function of the number of nodes N, for several values of the connectivity z. For the susceptibility, from Equation (6b), the slope of the linear fit gives the critical exponent ratio γ/ν. One can clearly see that the exponent ratio γ/ν, contrary to what happened to the β/ν exponent, now increases as the number of nodes increases, implying, nevertheless, in a non-universal behavior for different values of z, as before. In this case, a stronger finite-size dependence of the critical exponent ratio can also be noticed for higher values of the number of nodes. However, instead of taking into account corrections-to-finite-size scaling, which further needs the knowledge of the (unknown) correction-to-scaling exponent, we opted to analyze the behavior of the fitted exponent by systematically dropping data from the smaller lattices until the large node sizes regime has been achieved. That is why, for z greater than 50, only the larger number of nodes (N ≥ 1000) have been considered in the linear fit. The corresponding exponents are also shown in the first row of Table 1, including additional values of different parameters z.  The susceptibility, given by Equation (6b), can still be used, at least in principle, to obtain the value of the ratio γ/ν by using its maximum value χ max as a function of p, in this case at p max . It turns out, however, that p max is of the order p max ≤ 0.2 for the connectivity in the range 2 ≤ z ≤ 100. This is a quite smaller value when compared to the p c ≤ 0.44 obtained above from the cumulant crossings U 4 . In addition, even if we consider the maximum of the susceptibility in Equation (6b) regardless the value of p max , a different ratio exponent is obtained, as is depicted in Figure 4 for several values of z. One can see that γ/ν = 0.128(5) obtained from the data with p c is much smaller than γ/ν = 0.890 (10) computed from the maximum value of the susceptibility.
One can clearly notice that not only p c , coming from the maximum of the susceptibility, but also the corresponding critical exponent ratio γ/ν, are completely different from those values obtained from the fourth-order cumulant. It should be stressed that this very behavior is qualitatively the same as that obtained for the majority-vote model on the same DBAN [8]. However, the MVM on undirected Barabási-Albert networks (UBAN) exhibits a critical value of the connectivity and a corresponding critical exponent ratio which are of the same order whether estimated from χ max or χ(p c ) [8,12]. Although we do not have results of the BCS model on UBAN, the above results may be an indication that, independent of the dynamic model, the DBAN itself promotes this rather strange behavior of the order parameter fluctuation. This is further seen for other values of the connectivity z, as shown in Figure 5, and given in the first row of Table 1.  Open symbols correspond to the susceptibility evaluated at p c , while full symbols are evaluated at the maximum value of the susceptibility p max . The lines are linear fits with the corresponding exponents ratio given in Table 1. It is apparent from the results of Table 1 that both models have similar exponents ratio. In fact, in Figure 6 it is shown the critical exponent ratios β/ν and γ/ν, together with the effective dimensionality D e f f /2, as function of the connectivity z, with the data taken from Table 1 for both BCS model and MVM on the DBAN. It is clear from this figure that while the exponents ratio changes as z changes, the effective dimensionality remains the same and close to D e f f ≈ 1. Although the critical exponents for both models are close to each other, as expected by the Grinstein criterion [13], the agreement is better for the effective dimensionality D e f f . Figure 7 shows the phase diagram in the connectivity z and disorder p parameters plane for the BCS model compared to the MVM on the DBAN. In that scale we can notice that the non-universal critical disorder is indeed different for both models. Table 1. The critical connectivity parameter p c , the critical exponents ratio β/ν and γ/ν, and the effective dimension D eff on directed BA networks for different connectivities z. For each entry of the connectivity z, the upper row corresponds to the results for the discrete BCS model, while the lower row corresponds to the results for the MVM. The ratio γ/ν is obtained from two estimates: at p c and at p max . The former one comes from the values of χ at the critical concentration and the latter one from the maximum value of χ as a function of p. Error bars are statistical only.

Conclusions
We have studied a discrete version of the non-equilibrium kinetic Biswas-Chatterjee-Sen (BCS) model through extensive Monte Carlo simulations on directed Barábasi-Albert networks (DBAN). In contrast to the Ising model, which has no phase transition on DBAN, the BCS model shows indeed a continuous phase transition. Our Monte Carlo simulations indicate that the critical exponent ratios β/ν and γ/ν change as the connectivity of the network changes, and are in the same universality class as the majority-vote model (MVM). However, the effective dimensionality D eff is equal to unit, for all values of z, implying that the hyperscaling relation 2β/ν + γ/ν = 1 may be valid. A similar behavior has been found in references [8,12,14] for the MVM, even regarding the non-expected behavior of the maximum value of the order parameter susceptibility and its corresponding critical exponent ratio. This fact also indicates that the DBAN itself may be responsible for this strange behavior regarding the susceptibility quantity on this dynamical model. Of course, the study of other dynamic models defined on the DBAN would be very welcome to check this important point.