Modeling public opinion control by a charismatic leader

We study the average long-time behavior of the binary opinions of a social group with peer-to-peer interactions under the influence of an external bias and a persuadable leader, a strongly-biased agent with a dynamic opinion with the intention of spreading it across the system. We use a generalized, fully-connected Ising model, with each spin representing the binary opinion of an agent at a given time and a single, super spin representing the opinion of the leader. External fields and interaction constants model the opinion bias and peer-to-peer interactions, respectively, while the temperature $T$ models an idealized social climate, representing an authoritarian regime if $T$ is low or a liberal one if $T$ is high. We derive a mean-field solution for the average magnetization $m$, the"social mood", and investigate how $m$ and the super spin magnetization vary as a function of $T$. We find that, depending on the initial conditions, due to the presence of metastable states, the sign of the average magnetization depends on the temperature. Finally, we verify that this effect is also present even if we consider only nearest-neighbor interactions within the social group.


Introduction
Social interactions between humans are an example of a complex multi-component system.While the behavior of a single individual is quite complex, the collective behavior emerging from the interactions among peers might produce global patterns, which are independent of the characteristics of each individual [1,2,3].The Ising model has found application in social science, since it was first purposed by Weidlich [4] 50 years ago to analyze the statistics of opinion polarization.Later, the model was also applied to study critical properties in the collective behavior in a strike [1], and, more recently, a generalized version of this model considering clustered networks was used to study the evolution of binary opinions in society [5].A number of other statistical models have been developed to analyze collective social phenomena, which include the models proposed by Schelling [6] and Axelrod [7], and more physical ones, such as the voter model [8,9,10], the Sznajd model [11,12] and the majority-voter model [13], which has been shown to have the same critical exponents as the Ising model [14].Within the conceptual framework of social dynamics, a spin i represents an individual or agent, with its state or spontaneous magnetization σ i being its current (binary) opinion.An external field h i quantifies its bias to a certain opinion, while the interaction constants J ij quantify the peer-to-peer conditioning on personal opinions [5].The temperature T can be interpreted as a "climate parameter" [4], quantifying the random sources of external conditioning that may influence an individual opinion.In a simplified view, large T could be interpreted as a "liberal" regime, where opinions tend to fluctuate, while low T corresponds to a "totalitarian" regime, where opinion changes are more difficult to achieve [4].The Ising model can easily be modified to account for extra features.For instance, heterogeneity can be introduced in the form of a so-called "zealot" [15] or, sometimes, "stubborn" [16] agent: an agent with a fixed opinion [17,18] or with an opinion that evolves independently from its peers [15], with the intention of spreading it across the system.It is our purpose to investigate the effect of such a charismatic leader on the evolution of the global opinion as function of the degree of liberalism and initial conditions.
In this work, we study a generalized ferromagnetic Ising model which represents an idealized social group biased towards some specific opinion with peer-to-peer interactions among all possible spin pairs, as well as a super spin, a zealot-like agent with a large but finite bias to a contrarian opinion interacting with all other spins.A possible consequence of the introduction of zealots is the appearance of steady states where consensus (all spins in the same state) is not achieved [19,20], enriching the phenomenology seen in this kind of models.A crucial novelty in our model is that the super spin, unlike traditional zealots spins, also has the possibility to change its state or opinion over time due to the interactions with the social group.We will focus on the equilibrium behavior of the "social mood" [11], i.e. the average magnetization of the Ising model, as well as the particular average magnetization of the super spin, as a function of the climate parameter T and the initial distribution of opinions.The presence of heterogeneity in the model might give rise to metastable states, in which the system might reside for long times, giving rise to opinion distributions that do not minimize the energy of the system but are nonetheless important to consider due to their persistence [5].We also study the behavior of the model in systems with finite size, as it is usually the case in opinion models [21,3,18].This is because, in the context of social science, where the basic elements are humans, the system sizes considered are much smaller (at most of the order ∼ 10 9 if we consider the current global scale [22]) when compared to the ones considered in material science, where the elements are instead atoms or molecules, being of the order of the Avogadro constant ∼ 10 23 .Finally, we also consider a case where interactions in the social group are short-ranged.
The rest of the manuscript is organized as follows: section 2 describes the model in more detail and presents a mean field solution for the average magnetization.In section 3 we study finitesized systems using Monte Carlo simulations and compare numerical solutions with the mean field solution.Section 4 then summarizes the main results and presents possible future directions for the study of this model.

Model
We consider a system made of N bulk spins, under the influence of a weak external, negative field h < 0, with ferromagnetic connections J > 0 between every pair of bulk spins.An additional super spin s, under the influence of a local, strong positive field h s ≫ |h| is connected to all bulk spins through a strong ferromagnetic connection J s ≫ J.The Hamiltonian H of this generalized Ising model can be decomposed into two main components, with H s and H bulk referring to the super spin s and the bulk spins, respectively, where 1] is the state of the super spin s and the 1/2 in the last term of (3) takes into account the double counting of the pairs ij.The probability P (σ) to observe a particular state σ = {σ 1 , σ 2 , ..., σ N , σ s } at equilibrium is given by the Boltzmann distribution, where β = 1/(k B T ), with k B being the Boltzmann constant and T the temperature, and Z = σ e −βH is the partition function, where σ denotes a sum over all possible states σ.
Within the conceptual framework of social dynamics, the state σ i represents the current (binary) opinion of agent i.The bulk is the social group, with the external field h < 0 quantifying the intrinsic tendency or the exposure of this group towards a specific opinion [5], with the negative sign being just a convention, while J models the peer-to-peer interactions between bulk spins [5], with a tendency for consensus since J > 0. The super spin, on the other hand, represents a leader with a strong (h s ≫ |h|) and contrarian ideology (h s with opposite sign to h) to the rest of the social group, and has the means to be persuasive or indoctrinate efficiently (J s ≫ J), seeking to align the opinion of its peers with its own (J s > 0).
To solve this generalized Ising model, we will use the mean field approximation and assume that the sum of all the fluctuations around the average magnetization are negligible, where m ≡ ⟨σ i ⟩, ∀i, is the average magnetization per spin.To proceed, we can start by adding and subtracting the average magnetization m in the last term of Eq. ( 3), where in the last step we expanded the product and wrote explicitly the results of the sums Using approximation (5), we can disregard the last term in (6), and replace the result into Eq.( 3), approximating the Hamiltonian (1) as where /2 is a term that sets the energy scale.Notice that, in the thermodynamic limit N → ∞, the last term of Eq. ( 6) is effectively zero when compared to the first two, which are both proportional to N .Therefore the mean field Hamiltonian (7) becomes exact when N → ∞.We can then define a "mean-field" partition function To proceed it is useful to expand the summation over the two super spin σ s = ±1 states.We can then obtain the following expression for Z MF (see Appendix A), with where is the effective field felt by the bulk spins when the super spin is up (down), and we introduced the definitions h ′ s ≡ h s /N and J ′ ≡ (N − 1)J.From (8) we can define the probability p(+) or p(−) for the super spin to be up or down, respectively related by p(+) + p(−) = 1.As usual, the important physical observables of the system can be calculated from the derivatives of the logarithm of the partition function (8), such as the bulk magnetization m (per spin), Performing the derivative in Eq. ( 16) we find the mean field solution for the average bulk magnetization, m = p(+) e βA − e βB e βA + e βB + p(−) e βC − e βD e βC + e βD .(17) This is a self-consistent, transcendental equation for m, since the exponential terms A, B, C and D depend on m.It does not seem possible to find a closed form expression for m.Clearly, the first term corresponds to the average bulk magnetization when the super spin is up (p(+) = 1), while the second one is the average bulk magnetization when the super spin is down (p(−) = 1).Eq. ( 17) can be approximated, however, by simplifying the expression for the partition function (8) in the thermodynamic limit.First, notice that, the ratio (e βC +e βD )/(e βA +e βB ) < 1, and in the limit N → ∞, the term (e βC + e βD )/(e βA + e βB ) N goes to zero.Consequently, Applying this approximation into Eq.( 16) yields m ≈ e βA − e βB e βA + e βB = tanh [βh + ] , (21) where in the last step we cancelled out the common term e βh ′ s .Interestingly, comparing result (21) with Eq. ( 17) implies p(+) = 1 and p(−) = 0, i.e. in the thermodynamic limit, the super spin is expected to be always aligned with its local field h s if condition (19) is satisfied.On the other hand, if instead e βA + e βB < e βC + e βD , ( we have and This result, in turn, implies that p(−) = 1 and p(+) = 0, namely the super spin is expected to be always anti-aligned with its local field in the thermodynamic limit if condition ( 22) is satisfied instead.
In conclusion, to find the solutions of m for fixed values of the external parameters, one has to solve Eqs. ( 21) and ( 24) depending if temperature T and average bulk magnetization m satisfy, respectively, conditions (19) or (22).The behavior of the bulk, according to the mean field description in the thermodynamic limit, is consistent with a ferromagnetic system with N fully-connected spins, each pair interacting with strength J = J ′ /(N − 1), under the influence of an effective "external" field (±J s + h).The sign of J s depends on the state of the super spin which, as a function of T and m, may flip, acting like a switch.This is a consequence of the cooperative behavior emerging from the interaction between the super spin and the bulk.
From the partition function (8) one can also calculate the mean field expression for the free energy, Using the respective N → ∞ approximations ( 20) and ( 23), the (intensive) free energy is given by where + or − hold if condition (19) or ( 22) is satisfied, respectively.For the bulk parameters, we define h = −1 and J = 1/(N − 1).Notice that these definitions, i.e. h independent of N and J ∝ N −1 , are important to have the Hamiltonian H ∝ N and the free energy F ∝ N , and therefore well-defined extensive energies [23].For the super spin parameters, we define h s = N to achieve field neutrality, i.e. h ′ s = h s /N = |h| = 1.For the interaction constant, first notice that each bulk spin interacts with N − 1 other bulk spins with strength J, therefore, if all other spins are aligned, the energy contribution of the bulk interactions is of the order J • (N − 1) = J ′ = 1.We set J s = 1.01, in order to study a regime where the super spin interacts with the bulk with a slightly stronger interaction than the bulk interactions, even when global opinion accordance is achieved within the bulk.
In Fig. (1) we plot the mean-field solutions for m, Eqs. ( 21) and ( 24), as well as the free energy (26), for three values of the temperature T = {0.01,0.80, 1.50}, in the thermodynamic limit N → ∞.For low temperature, there are two solutions for m (Fig. In this regime, the super spin is expected to always be aligned with the bulk, as seen from the dark green curve (super spin always down, anti-aligned with its field), which intersects the bisector x = y at m ≈ −1, while the light green one (super spin always up, aligned with its field) intersects the bisector at m ≈ 1.At intermediate temperatures, the situation is quite different, as, independently of the bulk average magnetization, the super spin is expected to always be up, aligned with its field h s .The solutions of m at positive and negative magnetization are now |m| < 1 (Fig.   21) and ( 24), and the free energy FMF (26) as a function of m (d-f), for three values of T = {0.01,0.80, 1.50}, in the thermodynamic limit N → ∞.The light and dark green curves correspond to the regimes where conditions ( 19) and ( 22) are satisfied, respectively, and therefore the regimes where the super spin is predicted to be always up (p(+) = 1) or down (p(−1) = 1) in the thermodynamic limit.The vertical red dashed lines in (a-c) indicate the mean-field solutions of m (both stable and unstable ones), for which Eqs. ( 21) or ( 24) intersect the bisector x = y.The parameters are: h = −1, J ′ = 1, h ′ s = 1, Js = 1.01.

Numerical results
Next, we will study finite systems through the numerical estimation of the average bulk m and super spin m s magnetization from Monte Carlo simulations as a function of the temperature T , and compare with the mean field predictions ( 21) and (24), depending on which condition (19) or (22) is satisfied.
For the Monte Carlo simulations, we use the Metropolis algorithm [24].One Monte Carlo step corresponds to a full lattice sweep.For each temperature T and a given initial spin configuration, we thermalize by disregarding the first N therm = 2•10 4 steps, and then average over N MC = 10 4 spin configurations, storing spin configurations for averaging only every two steps, to reduce correlations between configurations.We repeat this procedure and average over N IC = 10 4 independent initial spin configurations.For low T , energy barriers between metastable states might be large compared to the thermal energy [23].In finite systems, it is understood that there exists always a finite time for the system to overcome these barriers and reach thermal equilibrium, but this time might be much longer than the Monte Carlo simulation times considered.A possible approach to check if we are averaging the bulk magnetization m over states which have different time averaged magnetization, when considering several independent initial spin configuration, is to split the measure of the average m into two components [23], where m + is the bulk magnetization averaged over several independent initial spin configurations which end up with a positive time-averaged bulk magnetization, while m − is the average over initial spin configurations which end up in turn with a negative time-averaged bulk magnetization.g + and g − correspond to the fraction of independent initial configurations that end up with either positive or negative average bulk magnetization, respectively.Naturally, they are constrained by the normalization condition g + + g − = 1.We split the average super spin magnetization in an analogous way, where the split components m s+ and m s− refer now to the average super spin magnetization over several independent initial configurations when its time average ends up positive or negative, respectively, and g s+ (g s− ) is the fraction of initial configurations that end up with positive (negative) average super spin magnetization.In Fig. (2) we plot the values of both m and m s , as well as their respective split components, as defined in Eqs. ( 27) and ( 28), as a function of temperature T , for a system with N ∼ 10 4 spins, for three different initial conditions: starting with a uniform distribution of σ i (no consensus), with all spins σ i = 1, aligned with the super spin field h s (consensus with the bias of the leader), and all spins σ i = −1, aligned with the bulk field h (consensus with the bulk bias), with the same respective initial conditions applied to the super spin σ s .
Looking first at the case starting with all σ i randomly assigned (Fig. (2)a-f), a configuration corresponding to no consensus, the average bulk magnetization (Fig. (2)a) exhibits different behavior depending on the temperature T .For low T ≲ 1, the system can evolve towards one of two stable states (Fig. while m − is a local one.Despite this, numerical results show that the state m − is a persistent metastable state for low T and, curiously, has a higher probability for the system to evolve to, as seen from Fig. (2)c, where g − > g + for low T .The mixed average of the metastable m < 0 and stable m > 0 states with corresponding uneven weights g − and g + therefore leads to a negative value in the average bulk magnetization seen in Fig. (2)a.This interestingly suggests that in an authoritarian regime, starting from a state with no consensus, it is more likely the bulk reaches consensus with the opinion of the peers and the leader adapts to the majority choice (Fig. (2)f).For T ≳ 1, the system undergoes a phase transition, and only the positive m + state is observed, as g + = 1 in Fig.  Results from Monte Carlo simulations of the average bulk m and super spin ms magnetization (first column), their partition into positive or negative average magnetization (middle column) and the fractions of Monte Carlo initial configurations which evolve towards a state with positive or negative average magnetization (last column), as a function of T ∈ [0.01, 1.50], for N ∼ 10 4 and three different initial conditions: no consensus (uniform distribution of {σi} and random σs, (a-f)), consensus with the super spin bias (all {σi}, σs = 1, (g-l)), and consensus with the bulk bias (all {σi}, σs = −1, (m-r)).For each point in T , results are averages over NIC = 10 4 independent Monte Carlo simulations.In plots (b), (h) and (n) we also show the mean-field solutions for m in the thermodynamic limit (colored dashed lines), as predicted by Eqs.(21) and (24) .The horizontal grey dashed lines indicate m = 0 or ms = 0.The parameters are: in the bulk magnetization, suggesting that more fluctuations are necessary to observe the change in opinion for the bulk.Notice also that, except for T ≈ 0.4, m s+ = 1 and m s− = −1, meaning that the super spin is either always up or down depending on which state the system evolves to, consistently with the mean field prediction in the thermodynamic limit.
For the case of initial configurations with all σ i = 1 (Fig.
(2)g-l), a configuration corresponding with consensus with the super spin bias, the long-time behavior of the system is quite different.Now, only the state with positive average bulk magnetization m s is observed (Fig. Lastly, for the case with spins starting all in the σ i = −1 state, corresponding to consensus with the bulk bias, the behavior is again quite different.The average bulk magnetization (Fig. (2)m-o) is negative for low T ≲ 1, with m = −1 indicating that all bulk spins are aligned with the bulk field h.For T ≳ 1, the state with negative m becomes unstable, and the average m becomes positive even when starting with all σ i = −1.The average super spin magnetization m s (Fig. (2)p-r) also changes sign with temperature.For low T ≲ 1, it remains aligned with the bulk field h (m s = −1), while for high T ≳ 1 it flips, remaining aligned with its field h s (m s = 1).Near the phase transition T ≈ 1, the state of the super spin fluctuates over time, as evidenced by the value of |m s | < 1 (Fig.
(2)p and q).Interestingly, unlike the case starting from a configuration with no consensus (Fig. (2)a-f), the transition temperature for the super spin is similar to that of the bulk, at T ≈ 1.It is also noteworthy that the numerical value for the partial average magnetization m − (Fig. ( 2)n) remains close to m − ≈ −1 until T ≈ 1, despite not being a minimum of the mean field free energy (blue dashed line).Similar results are obtained if we consider only nearest neighbour interactions in the bulk (see supplementary Fig. S1), where the effect of the changing sign of m with T seems to become less relevant as the system size increases (see supplementary Fig. S2).

Conclusions
We have studied the opinion dynamics in a social system under the influence of a charismatic leader by means of a generalized Ising model where all spins interact with a super spin via strong couplings.The interesting feature of the model is the appearance of a change in sign of both the bulk and leader average opinions as function of the social climate, behavior not usually found in Ising models.In summary, the mean-field approach (Fig. (1)) suggests that the state of total (m = 1) or partial (0 < m < 1) consensus with the leader is predicted to be the stable state, for any social climate T .For finite systems, numerical results (Fig. (2)) depend on the initial conditions for a totalitarian regime (low T ) due to metastable states corresponding to total (m = −1) or partial (−1 < m < 0) consensus with the bulk bias.When there is disagreement with the leader, i.e. no consensus (m ≈ 0) or consensus with the bulk (m = −1), the sign of the time-averaged social mood m depends on the social climate T .If the system starts from a state of no consensus (m ≈ 0, Fig. (2)a-f), the effect that the social mood m changes sign as function of temperature is due to an ensemble average over both metastable and stable states.Interestingly, the probabilities to evolve to either state are unequal (Fig. (2)c), with the metastable state being more probable than the stable one for low T .For a more liberal regime T ≳ 1, only the stable state of partial consensus with the leader (m > 0) is observed (Fig. (2)b-c) in the long-time behavior of the system.These results suggest that, if a social group is initially in a state of no consensus, agreement with the ideology of a leader is difficult, because it is more likely for the bulk to align its opinion with its own bias, when the social climate is authoritarian (low T ).For the initial condition corresponding to consensus with the bulk bias (m = −1, Fig. (2)m-q), m changes sign with T because the system always remains stuck in a persistent metastable state, aligned with the bulk ideology (m = −1), and, as the social climate changes to a more liberal regime, the system transitions into the equilibrium state of partial agreement with the leader (m > 0).This behaviour is also observed if only nearestneighbor interactions are considered in the bulk, corresponding to a scenario of limited social contact within the social group.Possible future directions for the study of this model include analyzing how the system behaves as a function of the super spin parameters h s and J s and investigating non-equilibrium properties, such as the transition rates, as is usually done in the context of opinion dynamics [3,14,10,5,25,26,27].

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
(1)a) with similar free energy (Fig. (1)d), corresponding to consensus aligned with the bulk field h, m ≈ −1, and opposite to it, m ≈ 1.
(1)b), indicating that consensus is no longer stable, as expected from a more liberal social climate.A new solution also appears near m ≈ 0, but the respective free energy plot (Fig.(1)e) indicates that this solution is unstable.Finally, at high temperatures, the system has only a single stable solution at m ≈ 0 (Fig.(1)c and f), which indicates indicates a social climate where opinions are highly varied and any kind of consensus is unstable.The super spin, as before, is expected to be always aligned with its field h s .

Figure 1 :
Figure 1: Mean field solution.Graphical solutions of the mean field self-consistent equation for for the average bulk magnetization m (a-c), Eqs.(21) and(24), and the free energy FMF (26) as a function of m (d-f), for three values of T = {0.01,0.80, 1.50}, in the thermodynamic limit N → ∞.The light and dark green curves correspond to the regimes where conditions (19) and (22) are satisfied, respectively, and therefore the regimes where the super spin is predicted to be always up (p(+) = 1) or down (p(−1) = 1) in the thermodynamic limit.The vertical red dashed lines in (a-c) indicate the mean-field solutions of m (both stable and unstable ones), for which Eqs. (21) or (24) intersect the bisector x = y.The parameters are: h = −1, J ′ = 1, h ′ s = 1, Js = 1.01.
(2)b), with either positive m + or negative m − average bulk magnetization.The numerical values of m + and m − are consistent with the mean field solutions for m (dashed lines in Fig. (2)b), where m + is predicted to be the global minimum of the free energy (see Fig.(1)b) (2)c, again consistent with the mean field solution (see Fig.(1)c and f), leading to a positive m in Fig. (2)a.Considering the average super spin magnetization m s (Fig. (2)d), for low T ≲ 0.4, two states can be observed, with either positive m s+ or negative m s− super spin magnetization (Fig. (2)e), with the negative one m s− being again more likely to be observed for low T , as g s− > g s+ in Fig. (2)f.For higher T ≳ 0.4, only the positive m s+ is observed, with g s+ ≈ 1 in Fig. (2)f.Notice that this occurs for a much lower temperature than the phase transition observed

Figure 2 :
Figure 2: Numerical results.Results from Monte Carlo simulations of the average bulk m and super spin ms magnetization (first column), their partition into positive or negative average magnetization (middle column) and the fractions of Monte Carlo initial configurations which evolve towards a state with positive or negative average magnetization (last column), as a function of T ∈ [0.01, 1.50], for N ∼ 10 4 and three different initial conditions: no consensus (uniform distribution of {σi} and random σs, (a-f)), consensus with the super spin bias (all {σi}, σs = 1, (g-l)), and consensus with the bulk bias (all {σi}, σs = −1, (m-r)).For each point in T , results are averages over NIC = 10 4 independent Monte Carlo simulations.In plots (b), (h) and (n) we also show the mean-field solutions for m in the thermodynamic limit (colored dashed lines), as predicted by Eqs.(21) and(24) .The horizontal grey dashed lines indicate m = 0 or ms = 0.The parameters are: h = −1, J ′ = 1, h ′ s = 1, Js = 1.01.
(2)h and i), corresponding to the global free energy minimum (dashed line in Fig.(2)h).Consequently, the average bulk magnetization (Fig.(2)g) no longer changes sign with temperature.For this case, the super spin is always up (see Fig.(2)j-l), aligned with its field h s .