# Equilibrium and dynamics of a three-state opinion model

Irene Ferri  Albert Díaz-Guilera  Matteo Palassini<EMAIL_ADDRESS>Departament de Física de la Matèria Condensada and Institute of Complex
Systems (UBICS)
Universitat de Barcelona, 08028 Barcelona, Spain

###### Abstract

Abstract: We introduce a three-state model to study the effects of a neutral
party on opinion spreading, in which the tendency of agents to agree with
their neighbors can be tuned to favor either the neutral party or two
oppositely polarized parties, and can be disrupted by social agitation
mimicked as temperature. We study the equilibrium phase diagram and the non-
equilibrium stochastic dynamics of the model with various analytical
approaches and with Monte Carlo simulations on different substrates: the
fully-connected (FC) graph, the one-dimensional (1D) chain, and Erdös-Rényi
(ER) random graphs. We show that, in the mean-field approximation, the phase
boundary between the disordered and polarized phases is characterized by a
tricritical point. On the FC graph, in the absence of social agitation,
kinetic barriers prevent the system from reaching optimal consensus. On the 1D
chain, the main result is that the dynamics is governed by the growth of
opinion clusters. Finally, for the ER ensemble a phase transition analogous to
that of the FC graph takes place, but now the system is able to reach optimal
consensus at low temperatures, except when the average connectivity is low, in
which case dynamical traps arise from local frozen configurations.

##  I Introduction

Within the field of complex systems, social questions are perhaps the most
elusive, as the agents involved (humans) exhibit a sophisticated individual
behavior, not easily reducible to a few analyzable parameters. Nevertheless,
many models have been proposed to capture different aspects of societal
interaction, such as bipartidism Abelson (1964); Deffuant  _et al._ (2002);
Banisch and Olbrich (2018); Yang  _et al._ (2020), gerrymandering Jacobs and
Walch ; Stewart  _et al._ (2019); Bergstrom and Bak-Coleman (2019), or echo
chambers formation Gaisbauer  _et al._ (2020); Kureh and Porter (2020);
Baumann  _et al._ (2020).

The consensus problem on a given social question, such as which kind of energy
is the most suitable for subsistence or which political party should govern,
has been addressed using a variety of agent-based models, both discrete, such
as the voter Holley and Liggett (1975); Suchecki  _et al._ (2005); Kureh and
Porter (2020), Ising Li  _et al._ (2019) or Potts models Schulze (2005), and
continuous, such as the Deffuant model Deffuant  _et al._ (2000). Continuous
models often predict the formation of opinion clusters Weisbuch (2004);
Sobkowicz (2015), thereby reinforcing a discrete description of the opinion
space.

A common goal in many of these works is to understand the transition between
an initial disordered state, in which opinions are random and uncorrelated, to
a state in which agents exhibit some kind of local or global consensus. The
simplest case occurs when there are only two opinions, as in polarized
situations with a clear bipartidist scenario. In other situations, however, it
is more realistic to consider at least one intermediate or neutral state
representing, for example, centrists or undecided voters. Several three-state
models have been proposed, using various approaches for introducing the
neutral state. Some of these models prevent agents that hold extreme opinions
from interacting directly, forcing them to pass through the neutral opinion.
For instance, Vazquez and Redner Vazquez and Redner (2004) study a stochastic
kinetic model in the mean-field approximation, and find that the final
configuration depends strongly on the initial proportion of agents in each
state. Along similar lines, the authors of Ref.Svenkeson and Swami (2015)
incorporate temperature and find a phase transition analogous to that of the
Ising model.

In this paper we propose a three-state Hamiltonian agent-based model for
opinion spreading in which agents interact in a pairwise manner that tends to
promote consensus, with a tunable neutrality parameter that controls the
relative preference for the neutral state over the polarized states. Agents
can change their opinion according to a stochastic dynamics in which the
effects of social agitation are taken into account by mimicking them as a
temperature.

The model can be mapped to a special case of the Blume-Emery-Griffiths (BEG)
model Blume  _et al._ (1971) from condensed matter physics. Other variants of
the BEG model have been applied to sociophysics before Yang (2010); Fernandez
_et al._ (2016), but our model allows to study directly the role of the
neutral state in the dynamics of the opinion formation.

The geometry of social structures is crucial in opinion formation and other
contemporary questions such as pandemic spreading, economics Souma  _et al._
(2003); Li  _et al._ (2003); Silva and Matsushita (2022) and smart cities
design Lu  _et al._ (2015); AlSonosy  _et al._ (2018). In order to understand
the role of the network geometry in the dynamics of opinion formation, we
embed our model on different types of networks: the fully connected (FC)
graph, the one-dimensional (1D) chain, and Erdös-Rényi (ER) random graphs. We
study with different analytical approaches the equilibrium phase diagram of
the model, and use Monte Carlo (MC) simulations at zero and non-zero
temperature to investigate the stochastic evolution of the population starting
from random configurations.

The paper is organized as follows. In Section II we introduce the model,
placing it in a social context and discussing its general features. In Section
III we determine the equilibrium phase diagram in the mean-field
approximation, which is characterized by a phase boundary with a tricritical
point. In Section IV we analyze the zero-temperature dynamics on the FC graph
and identify the basins of attraction of the absorbing states. We predict, and
confirm via MC simulations, that random configurations always evolve towards
an all-neutral state, and a finite temperature is necessary to overcome the
barriers towards equilibrium. In Section V we derive an exact solution for the
1D chain, and obtain the scaling of the magnetization from a domain-growth
argument, which is validated by our MC results. In Section VI, devoted to ER
graphs, we determine the transition temperature using the annealed mean-field
approximation, and show that for low connectivities the system gets stuck in
dynamical traps, while for high connectivities it is always able to reach the
optimal configuration. Finally, in Section VII we present our conclusions. The
Appendices contain details of the analytical calculations and numerical
methods.

##  II The Model

We propose an agent-based model with a discrete opinion space. Agents live on
the $N$ nodes of an undirected graph, and their opinions are represented by
two-dimensional vectors ${\bf S}_{i}$ ($i=1,\dots,N$) that can take three
orientations:

  * •

${\bf S}_{i}=(1,0)$; positive opinion / rightist

  * •

${\bf S}_{i}=(0,\alpha)$; neutral opinion / centrist

  * •

${\bf S}_{i}=(-1,0)$; negative opinion / leftist

where $\alpha$ is the neutrality parameter. We assume that the agents prefer
to agree with their neighbors so as to minimize, in the absence of social
agitation, the following cost function, or Hamiltonian:

$H=-J\displaystyle\sum_{\langle i,j\rangle}{\bf S}_{i}\cdot{\bf S}_{j},$ (1)

where $J>0$ and the sum runs over all the undirected edges of the graph. In
the optimal configuration, i.e. the ground state of the Hamiltonian, agents
reach consensus on the neutral opinion if $\alpha>1$, or on one of the two
polarized opinions if $\alpha<1$. The possible values of the interaction
energy between two agents, $-J{\bf S}_{i}\cdot{\bf S}_{j}$ are shown in Fig.
1. The quantity $J(\alpha^{2}-1)$ thus measures the reward of the neutral
opinion over the polarized ones. We assume that the system is in contact with
a thermal bath at temperature $T$, which should be understood as a coarse-
graining of all the sources of noise (e.g. social agitation) that affect
individual opinions.

An equivalent way to express the above Hamiltonian is to replace the state
vectors with scalar variables taking the values $\sigma_{i}\in\\{1,0,-1\\}$,
which gives

$H=-J\displaystyle\sum_{\langle
i,j\rangle}\left[\sigma_{i}\sigma_{j}+\alpha^{2}(\sigma_{i}^{2}-1)(\sigma_{j}^{2}-1)\right]\,.$
(2)

Figure 1: Contribution of a pair of interacting agents to the energy of the
system.

In this form, the model can be seen as a special case of the BEG model Blume
_et al._ (1971), which in its general form is described by the Hamiltonian

$H^{BEG}=-J\displaystyle\sum_{\langle
i,j\rangle}\sigma_{i}\sigma_{j}-K\displaystyle\sum_{\langle
i,j\rangle}\sigma_{i}^{2}\sigma_{j}^{2}+\displaystyle\sum_{i}\Delta_{i}\sigma_{i}^{2}+C\,.$
(3)

In fact, we recover Eq.(2) by setting

$K=J\alpha^{2},\quad\Delta_{i}=k_{i}K,\quad C=-NzK\mathbin{/}2,$ (4)

where $k_{i}$ is degree of node $i$, i.e. the number of nodes connected to it,
and $z=\sum_{i}k_{i}/N$ is the average degree of the graph. Note that the
standard BEG model has a unique value $\Delta_{i}=\Delta$ for all $i$.
Therefore, for graphs with constant degree $k_{i}=z$, we can obtain the phase
diagram of our model by projecting that of the BEG model, which is defined in
the three-dimensional parameter space ($K/J$, $\Delta/J$, $T/J$), onto the
two-dimensional semiplane ($\alpha,T/J$) defined by Eq.(4).

We assume that the system evolves stochastically via either one of two common
discrete-time Markov processes, the Metropolis and the Glauber dynamics, both
described in Appendix A. For $T>0$, both dynamics converge, given sufficient
time, to a stationary state in which the probability of a configuration
$\sigma=\\{\sigma_{i}\\}_{i=1}^{N}$ is given by the Boltzmann distribution
$p(\sigma)\propto\exp[-H(\sigma)/T]$. At $T=0$, depending on the graph and
initial conditions, they can either converge to the ground state of the
Hamiltonian, or get trapped forever in metastable configurations, as we will
discuss later. We perform MC simulations with both dynamics starting, unless
otherwise specified, from a random configuration in which the agents
independently take one of the three opinion states with uniform probability
(i.e. sampling from the infinite-temperature Boltzmann distribution), and then
quenching the system instantaneously to the desired temperature $T$ and
letting it evolve.

The order parameters of the model are the Ising-like magnetization
$m=\sum_{i}\langle\sigma_{i}\rangle/N$, namely the difference between the
fractions of rightists and leftists, and the fraction of neutral agents,
$n_{0}=1-\sum_{i}\langle\sigma_{i}^{2}\rangle/N$. Here $\langle\dots\rangle$
denotes the expectation with respect to the Boltzmann distribution.

In the rest of the paper we will use dimensionless units for the energy and
temperature such that $J=1$. We can anticipate some general features of the
equilibrium phase diagram in the plane $(\alpha,T)$. For $\alpha<1$, the
ground state is ferromagnetic, namely the agents achieve a global consensus
either in the positive or in the negative state ($m=\pm 1,n_{0}=0$). For
$\alpha>1$, all agents assume the neutral opinion ($m=0,n_{0}=1$) in the
ground state. Hence, moving along the zero-temperature axis we encounter a
discontinuous phase transition at $\alpha=1$.

For $\alpha=0$, the model can be thought of as an Ising model with vacancies.
Hence, in the thermodynamic limit it will generically display a continuous
phase transition in the Ising universality class at a critical temperature
$T_{c}^{0}>0$, between a low-temperature polarized (ferromagnetic) phase,
characterized by $m\neq 0$, and a high-temperature disordered (paramagnetic)
phase, in which $m=0$.

We thus generally expect a phase boundary in the plane $(\alpha,T)$ connecting
the two points $(\alpha=0,T=T_{c}^{0})$ and $(\alpha=1,T=0)$. Since the
transition is discontinuous at one end of the boundary and continuous at the
other, we also expect a tricritical point $(\alpha_{tc},T_{tc})$ at some point
along the boundary, separating a line of continuous transitions at
$T=T_{c}(\alpha)$ for $0\leq\alpha\leq\alpha_{tc}$ from a line of
discontinuous transitions at $T=T_{d}(\alpha)$ for $\alpha_{tc}<\alpha\leq 1$.
Such a phase boundary was indeed observed for a different projection of the
BEG model, the Blume-Capel model, on various types of random graphs, including
some cases in which the phase boundary is reentrant Leone  _et al._ (2002).
Exceptions to this scenario are represented by 1D systems with short-range
interaction, in which $T_{c}^{0}=0$, since no long-range order can survive at
finite temperature, and by graphs with a degree distribution falling more
slowly than $k^{-3}$ for large $k$, in which the system is known to remain
ferromagnetic at all temperatures for $\alpha=0$ Leone  _et al._ (2002) (we
expect this to remain true for all $\alpha<1$).

##  III Mean-field phase diagram

The mean-field approximation provides a useful first understanding of the
equilibrium behavior of the model. In this approximation, the free-energy per
spin is given by the minimum of a free-energy function ${\cal L}(m,n,\beta)$
with respect to $m$ and $n=1-n_{0}$, where $\beta=1/T$ is the inverse
temperature. As shown in Appendix A, the stationarity conditions
${\partial{\cal L}}/{\partial m}={\partial{\cal L}}/{\partial n}=0$ give two
coupled self-consistent equations (SCEs),

$\displaystyle m$ $\displaystyle=$ $\displaystyle\frac{2e^{\beta
z\alpha^{2}(n-1)}\sinh(\beta zm)}{1+2e^{\beta z\alpha^{2}(n-1)}\cosh(\beta
zm)},$ (5) $\displaystyle n$ $\displaystyle=$ $\displaystyle\frac{2e^{\beta
z\alpha^{2}(n-1)}\cosh(\beta zm)}{1+2e^{\beta z\alpha^{2}(n-1)}\cosh(\beta
zm)}.$ (6)

By expanding them for small $m$, we find a line of continous transitions
between the disordered phase and the polarized phase, at an inverse critical
temperature $\beta_{c}(\alpha)$ given by

$\alpha^{2}=\frac{1}{\beta_{c}(\alpha)z-1}\ln\left[2(\beta_{c}(\alpha)z-1)\right]\,.$
(7)

In particular, we have $\beta_{c}(0)^{-1}=2z/3$, which is below the critical
temperature $T_{c}=z$ of the Ising model in the mean-field approximation. This
is because, even if at $\alpha=0$ the neutral state does not contribute to the
energy, it brings an additional entropy that destabilizes the polarized phase.
The line of critical points, shown in Fig. 2, ends at a tricritical point
$(\alpha_{tc},\beta_{tc})$ determined by the condition

$2\ln[2(\beta_{tc}z-1)]=3-\beta_{tc}z$ (8)

By solving this numerically we obtain $\beta_{tc}^{-1}=0.532573\,z$, and
substituting this value into Eq.(7) gives $\alpha_{tc}=0.800354$.

Figure 2: Mean-field phase diagram. The line of continuous transitions
$T=\beta_{c}(\alpha)^{-1}$ for $\alpha<\alpha_{tc}$ is given by Eq.(7). For
$\alpha>\alpha_{tc}$, the points are obtained numerically as explained in the
text, the lines being only a guide to the eye.

For $\alpha>\alpha_{tc}$, the phase diagram displays a line of discontinuous
transitions at inverse temperature $\beta_{d}(\alpha)$, which we locate by
finding numerically the values of $m$, $n$, $\beta$ that satisfy
simultaneously Eqs.(5) and (6), together with the condition of equality
between the free energies of the ferromagnetic and paramagnetic phases. The
resulting phase boundary is shown in Fig. 2. Also shown in the figure is the
limit of metastability of the ferromagnetic phase, namely the value of $T$
above which the SCEs no longer admit a solution with $m\neq 0$.

The equilibrium values $\langle m\rangle$ and $\langle n_{0}\rangle=1-\langle
n\rangle$, obtained by solving numerically Eqs.(5) and (6), are shown in Fig.
3 as a function of temperature. For $\alpha<\alpha_{tc}$, the magnetization
vanishes continuously at the critical temperature, and is described by the
usual mean-field critical and tricritical exponents as $T\to
T_{c}(\alpha)^{-}$, namely $\langle m\rangle\sim(T_{c}(\alpha)-T)^{1/2}$ for
$\alpha<\alpha_{tc}$ and $\langle m\rangle\sim(T_{tc}-T)^{1/4}$ for
$\alpha=\alpha_{tc}$. The fraction of neutral agents $\langle n_{0}\rangle$
increases with $T$ up to the critical point $T_{c}(\alpha)$, then it decreases
monotonically towards $\langle n_{0}\rangle=1/3$ in the $T\to\infty$ limit, in
which the three states become equiprobable. We note that at $T=T_{c}(\alpha)$
we have $\langle n_{0}\rangle=1-T_{c}(\alpha)/z$.

For $\alpha_{tc}<\alpha<1$, $\langle m\rangle$ decreases monotonically with
$T$ and jumps to zero at $T=T_{d}(\alpha)$, while $\langle n_{0}\rangle$ has a
strongly non-monotonic temperature dependence: starting from $\langle
n_{0}\rangle=0$ at $T=0$, it increases slowly with $T$, then it jumps to a
large value at the discontinuous transition, before decreasing monotonically
towards $1/3$.

Finally, for $\alpha\geq 1$, we have $\langle m\rangle=0$ at all temperatures,
and $\langle n_{0}\rangle$ decreases monotonically with $T$, starting from
$\langle n_{0}\rangle=1$ at $T=0$, since in the ground state all agents are in
the neutral state.

From a social point of view, we see that agents agree on one of the polarized
opinions when coupling dominates over temperature, while at intermediate
levels of upheaval above the discontinuous transition they have a neutral
preference.

Figure 3: Mean-field results for the equilibrium values of the magnetization
(a) and the fraction of neutral sites (b), as a function of the temperature
and for different values of $\alpha$. The thick red line corresponds to
$\alpha=\alpha_{tc}$, separating the continuous transition from the
discontinuous transition. In b) the thick green curve corresponds to
$\alpha=1.2$, and the other curves are as in a).

The mean-field approximation is exact for the FC graph in the $N\to\infty$
limit, formally setting $z=1$ (see Section IV). In addition, the qualitative
features of the phase diagram and of the behavior of $\langle m\rangle$ and
$\langle n_{0}\rangle$ described above will hold in a large class of graphs.
In particular, for $d$-dimensional regular lattices we expect the same
qualitative phase diagram when $d\geq 2$ (no finite-temperature phase
transition can exist for $d=1$), with mean-field critical exponents for $d>4$
and non-mean-field exponents of the Ising universality class for $d=2,3$.

For random graphs, the mean-field approximation is not exact and various other
theoretical approaches have been developed. For the Ising model, it was found,
using the replica method, that if the degree distribution falls off as
$k^{-5}$ for large $k$ or faster (which includes the case of random-regular
and ER graphs), a continuous transition with mean-field exponents takes place.
On the other hand, for $\gamma\leq 5$ different scenarios depending on
$\gamma$ are observed Leone  _et al._ (2002). In the Blume-Capel model, a
tricritical point was found using an annealed mean-field approximation (see
Section C) De Martino  _et al._ (2012).

##  IV Fully-connected graph

Next, we discuss the stochastic dynamics of the model on the FC graph, in
which each agent interacts with all other agents. In this case, in order for
the Hamiltonian to be extensive (i.e. proportional to $N$), we must replace
the coupling constant $J$ in Eq.(1) by $J/N$. Since the degree is $z=N-1$, for
large $N$ we have $zJ/N=J$, which is equivalent to setting $z=1$ in the mean-
field solution. The Hamiltonian can then be written, up to terms of order
$1/N$ and recalling that $J=1$ in our units, as

$H=-\frac{N}{2}(m^{2}+\alpha^{2}n_{0}^{2})\,.$ (9)

The dynamics can then be represented as the evolution of a point $(m,n_{0})$
inside the triangle of vertices $(1,0),(-1,0),(0,1)$ shown in Fig. 4. We will
refer to any such point as the macrostate of the system, to distinguish it
from the microscopic configuration of the $N$ agents.

###  IV.1 Zero-temperature dynamics

If we sample the initial configuration of the $N$ agents independently and
uniformly at random among the three states, the probability distribution of
the initial macrostate is

$p(m,n_{0})=3^{-N}\frac{N!}{(Nn_{+})!\,(Nn_{-})!(Nn_{0})!},$ (10)

where $n_{\pm}=(1-n_{0}\pm m)/2$ is the fraction of agents in the positive and
negative state, respectively. For large $N$, one has
$p(m,n_{0})\sim\exp[Nf(m,n_{0})]$, where $f(m,n_{0})$ is a large-deviation
function that has a maximum at $m=0,n_{0}=1/3$, which corresponds to
equiprobable opinions, and is shown by the red point in Fig. 4. The
probability of any macrostate different than the equiprobable one decreases
exponentially with $N$.

It is nevertheless interesting to analyze the fate of the system prepared in
an arbitrary macrostate $(m,n_{0})$, even if exponentially rare, as it might
be relevant from a social viewpoint. At zero temperature, with both the
Glauber and the Metropolis dynamics the system can only evolve from any given
configuration by moves that do not increase $H$. In the following we will
focus on the Metropolis dynamics, in which an elementary move consists in
choosing an agent at random and proposing to change its state with probability
1/2 to either of the two states different from the current one, and accepting
the proposal if the energy change $\Delta H$ is negative or zero, in which
case a new macrostate $(m^{\prime},n_{0}^{\prime})=(m+\Delta m,n_{0}+\Delta
n_{0})$ is reached.

Since there are two possible proposals starting from each of the three opinion
states, there are six possible moves, except at the edges of the triangle
where some transitions are forbidden, and at the vertices where they are all
forbidden. The six moves are listed in Table 1 together with their proposal
probability $w_{r}$ ($r=1,\dots,6$), the displacement vector $(\Delta
m^{r},\Delta n_{0}^{r})$ multiplied by $N$, and the energy change $\Delta
H_{r}$ (up to terms of order $1/N$) upon accepting the move $r$. For example,
move $r=1$ consists in picking at random an agent $i$ that is currently in the
state $\sigma_{i}=-1$ and proposing to change its state to $\sigma_{i}=+1$.
The probability to propose this move is $n_{-}/2$ (since the fraction of
agents with $\sigma_{i}=-1$ is $n_{-}$ and we can choose to go to
$\sigma_{i}=1$ or 0 with probability 1/2). The total magnetization changes by
$+2$, the number of neutral sites remains unchanged, and the energy changes by
$-[(M+2)^{2}-M^{2}]/(2N)=-(2m+2/N)$.

All macrostates are unstable with respect to downhill moves (i.e. moves with
$\Delta H_{r}<0$) except the vertices of the triangle, which are absorbing
states. The lines $m=\pm\alpha^{2}n_{0}$ are the separatrices of the basins of
attractions of the absorbing states, as they discriminate between different
sets of allowed moves.

A full analysis of the master equation associated to the stochastic dynamics
would allow to determine the probability to end up in each of the absorbing
states for a given initial macrostate, but is beyond the scope of the paper.
In Fig. 4 we display instead the flux lines corresponding to the average
displacement, $(\langle\Delta m\rangle,\langle\Delta n_{0}\rangle)$, given by

$\langle\Delta m\rangle=\sum_{r=1}^{6}w_{r}\Delta m_{r}\theta(-\Delta
H_{r})\,,$ (11)

and similarly for $\langle\Delta n_{0}\rangle$, where $\theta$ is a Heaviside
function ensuring that only downwill moves contribute. We can see that below
the separatrices the average flux flows to the polarized $m=\pm 1$ states, and
above the separatices it flows to $n_{0}=1$.

In more detail, we observe that if $\alpha<1$ and $m>0$ (the case $m<0$
follows by symmetry) all the allowed moves from a point below the separatrix
(moves 1, 3 and 5, as shown in Fig. 4) displace the point further away from
the separatrix. Thus all the points below the separatix belong to the basin of
attraction of the $m=1$ absorbing state and have probability zero to flow to
$n_{0}=1$. In contrast, a point above the separatrix can move away from the
separatrix with probability $n_{+}/2$ (move 4) or move towards it with total
probability $n_{-}$ (the sum of the probabilities of moves 1 and 3). If
$n_{0}<1-3m$, the former probability is smaller than the latter, thus the
preferred direction is towards the separatrix. Therefore, a starting point in
the triangle defined by $m/\alpha^{2}<n_{0}<1-3m,m>0$ will have a non-zero
probability to end up in the “wrong” absorbing state $m=1$. However, since
there are many more paths towards $n_{0}=1$ than towards $m=1$, this
probability will in fact decrease exponentially with $N$.

Figure 4: Evolution of the Metropolis dynamics on the FC graph. The red lines are the separatrices $m=\pm n_{0}\alpha^{2}$ of the three basins of attraction (shown here for $\alpha=0.8$). The red point shows the equiprobable macrostate $m=0,n_{0}=1/3$. The solid black arrows represent the allowed moves ($\Delta H\leq 0$) from a macrostate below (green point) and above (blue point) the separatrix for $m>0$. The dashed orange arrows represent the forbidden moves ($\Delta H>0$). The light grey lines show the average direction of the evolution at each point. $r$ | Transition | $w_{r}$ | $N(\Delta m^{r},\Delta n_{0}^{r})$ | $\Delta H_{r}$
---|---|---|---|---
1 | $-\longrightarrow+$ | $n_{-}/2$ | $(2,0)$ | $-2m$
2 | $+\longrightarrow-$ | $n_{+}/2$ | $(-2,0)$ | $2m$
3 | $-\longrightarrow 0$ | $n_{-}/2$ | $(1,1)$ | $-m-\alpha^{2}n_{0}$
4 | $+\longrightarrow 0$ | $n_{+}/2$ | $(-1,1)$ | $m-\alpha^{2}n_{0}$
5 | $0\longrightarrow+$ | $n_{0}/2$ | $(1,-1)$ | $-m+\alpha^{2}n_{0}$
6 | $0\longrightarrow-$ | $n_{0}/2$ | $(-1,-1)$ | $m+\alpha^{2}n_{0}$
Table 1: Allowed transitions for the Metropolis algorithm. The second column
shows the transition between two spin values. $w_{r}$ is the probability of
proposing the transition, $N\Delta m^{r}$, $N\Delta n_{0}^{r}$, $\Delta H_{r}$
are the changes in total magnetization, number of neutral spins, and
Hamiltonian upon executing the transition $r$.

We verified the above predictions by performing repeated MC simulations with
the Metropolis dynamics at $T=0$ and estimating, for every possible initial
macrostate, the probability to reach each of the three aborbing states. In
Fig. 5 (left column) we display our results for $\alpha=0.8$. Indeed, points
below the separatrices end up in the $m=\pm 1$ absorbing states, while above
the separatrices we observe regions of “mixed fate” points, which we define as
those with a probability larger than 1.5% to end up in a state different than
their “natural” absorbing state (for example, points displayed in green can
end up in $m=1$ instead of $n_{0}=1$). The width of the mixed fate regions
decreases with $N$, approximately as $1/N$ as implied by the previous
argument. We also verified that the probability of a point above the
separatrix to end up in $m=1$ decreases exponentially with its distance from
the separatrix (not shown). An analogous analysis for the Glauber dynamics
shows that, since there are fewer allowed downhill moves, the probability to
cross the separatrix is even smaller and thus the mixed-fate region is
thinner, as confirmed numerically in Fig. 5 (right column).

Figure 5: Basins of attraction at $T=0$ and $\alpha=0.8$ on the FC graph. (a)
$N=50$, Metropolis dynamics; (b) $N=50$, Glauber dynamics; (c) $N=200$,
Metropolis dynamics; (d) $N=100$, Glauber dynamics. The lines are
separatrices. Symbols at each point indicate the fate of the system starting
from that point. Results are obtained with $1000$ MC runs for each point above
the separatrix, and $40$ for each point below it.

It is also interesting to observe that a deterministic downhill dynamics that
follows the direction opposite to the gradient of $H$, $-\nabla
H=N(m,\alpha^{2}n_{0})$, also has the same separatrices.

From a social perspective, the above results show that for $\alpha<1$, when we
start from equiprobable opinions ($m=0,n_{0}=1/3$), the population evolves
towards neutrality ($n_{0}=1$), despite the optimal configurations are the
polarized states. As we show below, social agitation is necessary in order to
overcome the energy barrier and reach optimality.

###  IV.2 Finite-temperature dynamics

We performed MC simulations at $T>0$ with the Glauber dynamics with $N=500$
agents, again starting from a random configuration. When $\alpha<\alpha_{tc}$
or $\alpha>1$, the system is able to equilibrate at all temperatures $T>0.1$
in less than $10^{3}$ Monte Carlo steps (MCS), where 1 MCS = $N$ elementary
moves. Our MC estimates for $\langle m\rangle$ and $\langle n_{0}\rangle$
agree with the exact mean-field results displayed in Fig. 3, with negligible
finite-size effects.

For $\alpha_{tc}<\alpha<1$, at temperatures below the phase boundary in Fig.
2, we expect metastability effects due to the discontinuous transition:
starting from a random configuration, the system gets trapped in a region of
the configuration space with $m\simeq 0$ for a time that diverges
exponentially with $N$ and with $1/T$.

This is confirmed by inspecting the probability distribution of the
magnetization obtained from many MC runs, shown in Fig. 6. The right column
shows results for $\alpha=0.9$, which is well into the discontinuous region
and for which the transition temperature is $T_{d}=0.47$. We see that for
$T=0.4$, after $10^{3}$ MCS the distribution is bimodal, with one peak at the
mean-field equilibrium magnetization $\langle m\rangle=0.83$, and the other
near $m=0$, corresponding to runs stuck in the metastable state. After
$10^{6}$ MCS, the system is able to overcome the free-energy barrier and the
peak at $m=0$ disappears. In contrast, for $T\leq 0.3$ even after $10^{6}$ MCS
the distribution remains peaked around $m=0$, while the equilibrium value is
close to $m=1$. This shows that almost all agents remain in the neutral state
(similarly to the $T=0$ case analyzed in the subsection IV.1), and social
agitation is not enough to overcome the free energy barrier towards polarized
consensus.

This contrasts with the results for $\alpha=0.75$, shown in the left column of
Fig. 6, at which the continuous transition takes place at $T_{c}=0.56$. We now
see that for $T\geq 0.3$, already after $10^{3}$ MCS the distribution is
peaked around the mean-field equilibrium value ($\langle
m\rangle=0.95,0.85,0.63$ for $T=0.3,0.4,0.5$ respectively), which shows that
in this case thermal fluctuations are able to bring the system towards
consensus. Only at $T=0.1$ we observe that the system is more likely to reach
$m\simeq 0$ than the equilibrium value $\langle m\rangle\simeq 1$, even after
$10^{6}$ MCS. At such low temperature, the probability to flip a polarized
agent surrounded by agents of opposite sign is exponentially small.

Figure 6: Probability distribution of the magnetization for two values of
$\alpha$ below (left columb) and above (right column) the tricritical point
for different temperatures obtained with Glauber MC simulations with $N=500$
agents and $10^{4}$ repetitions. The first row show results for a
computational time of $10^{3}$ MCS and the second row for $10^{6}$ MCS.

##  V One-dimensional lattice

We now study the model on a 1D chain in which each agent interacts only with
its two nearest neighbors. The main interest of this case is that it is
exactly solvable and it is the furthest from the FC graph. In Appendix D we
solve the model using the transfer matrix method, adopting periodic boundary
conditions (it is straightforward to consider other boundary conditions, such
as open or free, and the results for large $N$ are not affected by the
choice). In the large-$N$ limit, we obtain the free energy per site

$\displaystyle f(\beta,\alpha)$ $\displaystyle=$
$\displaystyle-\frac{1}{\beta}\ln\left[\cosh\beta+e^{\beta\alpha^{2}}/2\right.$
(12) $\displaystyle+$
$\displaystyle\left.\sqrt{(\cosh\beta-e^{\beta\alpha^{2}}/2)^{2}+2}\,\,\right].$

This is an analytic function of its arguments, thus it does not display any
finite-temperature phase transition, as expected in any 1D models with short-
range interactions. Only in the zero-temperature limit there is a phase
transition at $\alpha=1$. In fact, the limit of $f(\beta,\alpha)$ for
$\beta\to\infty$ gives the ground-state energy per site, which is $-1$ for
$\alpha\leq 1$, and $-\alpha^{2}$ for $\alpha>1$.

The transfer-matrix solution also provides an explicit expression for the
correlation function $C(r)=\langle\sigma_{i}\sigma_{i+r}\rangle$. As shown in
Appendix D, in the $N\to\infty$ limit we obtain
$C(r)=c_{0}\exp[-r/\xi(\beta)]$ where $c_{0}$ is a prefactor that depends on
$\beta$ and $\alpha$. The correlation length $\xi(\beta)$ is given by

$\displaystyle\frac{1}{\xi(\beta)}$ $\displaystyle=$
$\displaystyle\ln\left(\cosh\beta+e^{\beta\alpha^{2}}/2\right.$
$\displaystyle+$
$\displaystyle\left.\sqrt{(\cosh\beta-e^{\beta\alpha^{2}}/2)^{2}+2}\right)-\ln
2\sinh\beta\,.$

For $\alpha<1$ the correlation length diverges exponentially at low
temperatures as

$\xi(\beta)\underset{\beta\gg 1}{\sim}\frac{1}{4}e^{2\beta}\,,$ (14)

where the exponent stems from the fact that the energy of a domain wall
between a positively polarized domain (i.e. a segment of contiguous agents
with $\sigma_{i}=1$) and a negatively polarized domain is $2$.

We now use this result to deduce the behavior of the magnetization. If
$N\gg\xi$, we can consider the system as composed of $N/\xi$ independent
domains of alternating sign, each with a length of order $\xi$. The
magnetization is thus $m\simeq\sum_{k=1}^{N/\xi}M_{k}/N$, where the domain
magnetizations $M_{k}\simeq\pm\xi$ can be treated as independent random
variables with variance of order $\xi^{2}$. From this we obtain the scaling
law

$\langle|m|\rangle^{2}N\sim\xi(\beta).$ (15)

On the other hand, when $N<\xi$, if we let the system evolve with the Glauber
or Metropolis dynamics starting from a random configuration, at infinite time
it will reach a magnetized state with $\langle|m|\rangle\simeq 1$. The
evolution towards this state takes place via a domain-growth process (known as
coarsening in statistical physics), in which the walls between domains of
different sign perform a random walk. When two domain walls collide, they
annihilate each other, fusing the domains, so $\langle|m|\rangle$ will grow in
time until approaching $\simeq 1$.

The results of our MC simulations agree very well with the above predictions.
Fig. 7(a) shows, for $N=1000$ and $\alpha=0.75$, $\langle|m|\rangle$ as a
function of temperature, for different simulation times. Indeed one can
observe that at high $T$ the data are in thermal equilibrium, while at low
$T$, $\langle|m|\rangle$ grows with time until approaching unity.

Fig. 7(b) shows the equilibrium value of the magnetization achieved after a
sufficiently long time, for different $N$: as $N$ grows, the magnetization
saturates when $T$ is below the temperatures at which $\xi$ overcomes $N$.
Fig. 8(a) shows that the data for different $N$ in the regime $N\gg\xi$ can be
rescaled very well according to Eq.(15).

Figure 7: Monte Carlo results for the average absolute magnetization of the 1D
chain as a function of temperature, for $\alpha=0.75$. For each data set, the
initial 1/5 of the time series was discarded, and data were averaged over
100-400 repetitions. Periodic boundary conditions were applied. (a) Data for
$N=1000$ and different simulations times, starting from a random
configuration. (b) Data for different $N$, for sufficiently large simulation
time so that $\langle|m\rangle$ is equilibrated ($5\times 10^{4},5\times
10^{5},5\times 10^{6}$ MCS for $N=100,500,1000$ respectively). Lines are only
a guide to the eye.

Figure 8: (a) Scaling of the average absolute magnetization of the 1D chain
for $\alpha<1$. Data are the same as in Fig. 7(b), rescaled according to
Eq.(15). The curves show the full expression of the correlation length Eq.(V)
and the asymptotic behavior Eq.(14). (b) Same as a) but for $\alpha=1$ and
scaling the data according to Eq.(17). The data for $L=5000$ were obtained
with $10^{6}$ MCS, the others with the same number of MCS as in (a). The
curves show the full expression of the correlation length Eq.(V) and the
asymptotic behavior for $\alpha=1$, Eq.(16).

At the transition point $\alpha=1$, the behavior is somewhat different. In
fact, the correlation length now diverges at low temperatures as

$\xi(\beta)\underset{\beta\gg 1}{\sim}\frac{1}{\sqrt{2}}e^{\beta}\,,$ (16)

since the energy of a domain wall between neutral domains and polarized
domains is $1$. Following the same argument given above for $\alpha<1$, the
domain magnetizations can be either $M_{k}\simeq 0$ or $M_{k}\simeq\pm\xi$
with equal probability, and thus have variance $\simeq 2\xi^{2}/3$, which
gives

$\langle|m|\rangle^{2}N\sim\frac{2}{3}\xi(\beta).$ (17)

In this case, for $N<\xi$, the equilibrium magnetization saturates to a value
less than $1$ as the system can evolve towards either one of the three
consensus states.

Finally, the case $\alpha>1$ is simpler as now there is a unique ground state,
thus the system evolves quickly towards the equilibrium configuration.

##  VI Erdös-Rényi graphs

Although the ER random graph ensemble is not representative of real-world
structures, it is often used as a null model to study complex networks Newman
_et al._ (2002). We obtain the equilibrium behavior of the model (averaged
over the ER ensemble) using the annealed approximation Bianconi (2002);
Dorogovtsev  _et al._ (2008), which consists in applying the mean-field SCEs
in Eqs.(5) and (6)) taking into account that the degree $k$ has a Poisson
distribution $P(k)=e^{-z}z^{k}/k!$ with average degree $z$ Erdös and Rényi
(1959). In this way we obtain

$\displaystyle m$
$\displaystyle=\sum_{k}P(k)\frac{2e^{\beta\alpha^{2}k(n_{w}-1)}\sinh(\beta
km_{w})}{1+2e^{\beta\alpha^{2}k(n_{w}-1)}\cosh(\beta km_{w})},$ (18)
$\displaystyle n$
$\displaystyle=\sum_{k}P(k)\frac{2e^{\beta\alpha^{2}k(n_{w}-1)}\cosh(\beta
km_{w})}{1+2e^{\beta\alpha^{2}k(n_{w}-1)}\cosh(\beta km_{w})}.$ (19)

where $m_{w}=\sum_{j}k_{j}m_{j}/zN$ and $n_{w}=\sum_{j}k_{j}n_{j}/zN$ are
weighted order parameters that satisfy two coupled SCEs (see Appendix C for
details).

###  VI.1 Phase transition

Figure 9: Magnetization vs. temperature on ER graphs of average degree $z=4$
(a) and $z=16$ (b) for $\alpha=0$ and $\alpha=0.4$. The lines corresponds to
the numerical solution in the annealed approximation, while the bars are the
results of MC simulations with $10^{4}$ MCS, averaged on 100 realizations of
different ER graphs with $N=1000$ agents.

Also in this case one expects in general a tricritical point separating a line
of continuous transitions for low $\alpha$ from a line of discontinuous
transitions at large $\alpha$, as in the Blume-Capel model on ER graphs Leone
_et al._ (2002). We determined the location of the continuous transition by
solving numerically the SCEs for $m_{w},n_{w}$ using a Broyden first Jacobian
approximation from the optimize root of SciPy. Fig. 9 shows the comparison
between the numerical solution and MC simulations for $\alpha=0$ and $0.4$.
For low connectivity (panel (a)), we found a significant discrepancy between
the annealed approximation and the simulations, while for high connectivity
(panel (b)) we obtain a better agreement.

Figure 10: Critical temperature vs. $1/z$ for different values of $\alpha$,
calculated by solving numerically the system of equations (18), (19). Mean-
field results are represented by bigger points at $z\to\infty$.

Fig. 10 shows that the critical temperature, normalized by $z$, converges
asymptotically for large $z$ to that of the FC graph. The critical temperature
decreases with $\alpha$ for any $\langle k\rangle$, as in the mean-field
solution.

The putative tricritical point and the discontinuous transition can in
principle be located from a Taylor expansion in $m$ of the annealed SCEs. We
have have not attempted this, but we note that it is more complicated than in
the Blume-Capel model Leone  _et al._ (2002), in which the SCE $m$ decouples
from $n$.

###  VI.2 Zero-temperature dynamics

We performed MC simulations with the Metropolis and Glauber dynamics at $T=0$,
starting from a random configuration as in the FC graph. We consider graphs
with only one connected component, created by generating ER graphs with the
Python $networkx$ library and then randomly adding and/or subtracting agents
and links until the desired $N$ and $z=\langle k\rangle$ are reached,
preserving the degree distribution of the original graph.

We perform many runs for each graph, letting the system relax until it reaches
a steady state in which the energy no longer changes. In this way we collect
the probability distribution $P(\epsilon)$ of the residual energy, defined as
the normalized difference between the energy of the steady state reached in a
given run and the ground-state energy of the graph,
$\epsilon=(E_{steady}-E_{0})/Nz$. Fig. 11 shows the results obtained with the
Glauber dynamics for $z\leq 4$ and different values of $\alpha$. For $z=2$, in
all cases $P(\epsilon)$ has a single peak at $\epsilon>0$, indicating that the
system never reaches the ground state as it gets stuck in a manifold of
excited isoenergetic configurations, the energy of which changes from run to
run. The same behavior persists for $\alpha>1$ (see Appendix E for $N=500$).

A similar behavior has been observed before for the Ising model Baek  _et al._
(2012), and Häggström Häggström (2002) showed rigorously that it stems from
the existence of an extensive number of subgraphs in which some nodes are
frozen (i.e. they cannot change state without increasing the energy) in an
excited configuration. In addition, some nodes are “blinkers”, namely they can
change state forever without changing the energy, thus the system gets trapped
in a manifold of configurations at constant energy above the ground state. As
an illustration, the two central nodes of Fig. 12 a) are frozen, while the
middle node in Fig. 12 b) is a blinker.

Upon increasing $z$, the number of frozen nodes and blinkers decreases, hence
so does the probability to get stuck in dynamical traps. For $z=3$, as we
increase $\alpha$ we see that the distribution is bimodal for $\alpha\leq
0.75$, with both peaks at $\epsilon>0$ (the ground state is never reached),
then the peak closer to $\epsilon=0$ disappears for $0.85\leq\alpha\leq 0.95$,
and finally the distribution becomes unimodal with a finite weigth at
$\epsilon=0$ for $\alpha>1$ (see Appendix E). Note that these dynamical traps
are not a finite-size effect: as shown in Fig. 11 d) for $z=3$ and $\alpha=0$,
the peak near $\epsilon=0$ decreases with the system size, and the one at
larger $\epsilon$ grows.

For $z=4$, we observe a large probability $P(0)$ to reach the ground state for
all values of $\alpha$. For $\alpha<1$ there is still a peak at $\epsilon>0$
but it is significantly smaller than that for $z=3$. Notice that for both
$z=3$ and $z=4$ the peak at large $\epsilon$ increases significantly in the
range $0.85\leq\alpha\leq 0.95$: it is an interesing question whether this is
related to existence of a tricritical point.

Finally, for the Metropolis dynamics we obtain very similar results to the
ones just discussed. A comparison between the two dynamics for several values
of $\alpha$, $z$ and $N$ is shown in the Appendix E.

The above results show that the relaxational $T=0$ dynamics on ER graphs is
quite different from the FC case in which, as shown in Section IV, for large
$N$ the system prepared in a random configuration always reaches the ground
state if $\alpha>1$ (namely the absorbing state $n_{0}=1$), and never reaches
it if $\alpha<1$. In ER graphs, instead, the system is able to reach the
ground state for all $\alpha$, except for small values of $z$.

Figure 11: Probability distribution function of the residual energy
$\epsilon$, for the steady state reached through Glauber dynamics on a ER
graph with $N=3000$ and $z=\langle k\rangle=2,3,4$ for (a) $\alpha=0$, (b)
$\alpha=0.5$, and (c) $\alpha=1$. (d) Comparison of the probability
distribution function for $\langle k\rangle=3$ and different values of $N$.
Data are obtained from $10^{5}$ different random initial conditions on a
single realization of the network for $N=500$ and $2\cdot 10^{4}$ different
random initial conditions for $N=1000,3000$. Figure 12: Example of subgraphs
with a stuck spin configuration, despite not being compatible with the ground
state.

##  VII Conclusions

We have studied a three-state Hamiltonian model aimed at understanding
consensus formation in the presence of neutral agents between opposite
extremes. The tendency of people to align their opinions with those of their
friends is captured by a pairwise Ising-like interaction including a neutral
orientation weighted by a neutrality parameter $\alpha$. This tendency is
opposed by individual thinking, social agitation and other factors, here
collectively represented as a temperature.

Regardless of the network type, at high temperatures the agents do not feel
the influence of their peers and evolve quickly towards a disordered state,
that can be interpreted as social unrest. At low temperature, for $\alpha>1$
the system reaches a configuration dominated by neutral agents on all
networks, while for $\alpha<1$ we observe different behaviors depending on the
network type, and on the value of $\alpha$, as we summarize below.

On the FC graph, the equilibrium phase diagram of the model exhibits a phase
boundary for $\alpha<1$ between a low-temperature phase, in which one of the
polarized opinions prevail, and the high temperature disordered phase. At
temperatures just above the transition, we observe a majority of neutral
agents and equal-sized minorities of extremists of both signs. Upon crossing
the phase boundary by decreasing the temperature, the fraction of polarized
agents increases continuously from zero when $\alpha$ is below the tricritical
value $\alpha_{tc}$, while it jumps discontinuously to a non-zero value when
$\alpha_{tc}<\alpha<1$. We also found that in the absence of social agitation,
starting from a random configuration, the population is exponentially more
likely to get stuck in neutral consensus than to reach the optimal polarized
consensus. Moderate social agitation allows the system to reach polarized
consensus for $\alpha<\alpha_{tc}$, while for $\alpha<\alpha_{tc}<1$ the
barriers to achieve equilibrium are much higher due to metastability.

It is interesting to compare our results with those of other three-state
models on the FC graph. Vazquez and Redner Vazquez and Redner (2004) studied a
kinetic model in which pairs of agents interact stochastically. Their model
exhibits the same three absorbing states, but the largest basin of attraction
is that of an absorbing boundary on the bottom line of the triangle (see
Fig.4), representing frozen mixtures of oppositely polarized agents. A large
bipolarized region was also found by Balenzuela et al.Balenzuela  _et al._
(2015), who proposed a kinetic model in which agents hold a continuous
spectrum of convictions, which is partitioned in three states according to
some thresholds. The pair interaction is such that oppositely polarized agents
tends to increase their polarization, leading to a phase transition from a
neutral to either a bipolarized or a polarized population. Svenkeson and Swami
Svenkeson and Swami (2015) considered a pair dynamics, in which polarized
agents can become neutral and viceversa, with a temperature coupled to the
instantaneous magnetization. They find a transition of the Ising type without
a tricritical point.

Our results for the one-dimensional chain, which obviously does not display a
phase transition, show that the evolution towards consensus takes place via
the growth of domains of contiguous like-minded polarized agents, until one
domain takes over the whole population.

Finally, in ER random graphs, which are somewhat closer to a real population
in which agents have a finite number of contacts, we observe a finite-
temperature transition to a polarized phase, analogous to that of the FC
graph. However, unlike in the FC graph, for large enough connectivities the
system is always able to reach the polarized state from a random start, while
for low connectivities it gets stuck in dynamical traps due to frozen nodes.

The behavior of the model on networks that more realistically represent human
connection patterns, either synthetic or extracted from real data, will be
considered elsewhere.

###### Acknowledgements.

We acknowledge financial support from MINECO via Project No.
PGC2018-094754-B-C22 (MINECO/FEDER,UE) and Generalitat de Catalunya via Grants
No. 2017SGR341 and 2017SGR1614. We are pleased to thank Conrad Pérez Vicente
for discussions.

## References

  * Abelson (1964) R. Abelson, in  _Contributions to mathematical psychology_ , edited by N. Fredericksen and H. Gullicksen (Holt, Rinehart & Winston, New York, 1964).
  * Deffuant  _et al._ (2002) G. Deffuant, F. Amblard, G. Weisbuch, and T. Faure,  _How can extremism prevail? A study based on the relative agreement interaction model_ , J. Artif. Soc. Soc. Simul. 5, 1 (2002).
  * Banisch and Olbrich (2018) S. Banisch and E. Olbrich,  _Opinion polarization by learning from social feedback_ , J. Math. Sociol. 43, 76 (2018).
  * Yang  _et al._ (2020) V. C. Yang, D. M. Abrams, G. Kernell, and A. E. Motter,  _Why are U.S. parties so polarized? A “satisficing” dynamical model_ , SIAM Rev. Soc. Ind. Appl. Math. 62, 646 (2020).
  * (5) M. Jacobs and O. Walch,  _A partial differential equations approach to defeating partisan gerrymandering_ , arXiv:1806.07725 [physics.soc-ph] .
  * Stewart  _et al._ (2019) A. J. Stewart, M. Mosleh, M. Diakonova, A. A. Arechar, D. G. Rand, and J. B. Plotkin,  _Information gerrymandering and undemocratic decisions_ , Nature 573, 117 (2019).
  * Bergstrom and Bak-Coleman (2019) C. T. Bergstrom and J. B. Bak-Coleman,  _Information gerrymandering in social networks skews collective decision-making_ , Nature 573, 40 (2019).
  * Gaisbauer  _et al._ (2020) F. Gaisbauer, E. Olbrich, and S. Banisch,  _Dynamics of opinion expression_ , Phys. Rev. E 102, 042303 (2020).
  * Kureh and Porter (2020) Y. H. Kureh and M. A. Porter,  _Fitting in and breaking up: A nonlinear version of coevolving voter models_ , Phys. Rev. E 101, 062303 (2020).
  * Baumann  _et al._ (2020) F. Baumann, P. Lorenz-Spreen, I. M. Sokolov, and M. Starnini,  _Modeling echo chambers and polarization dynamics in social networks_ , Phys. Rev. Lett. 124, 048301 (2020).
  * Holley and Liggett (1975) R. A. Holley and T. M. Liggett,  _Ergodic theorems for weakly interacting infinite systems and the voter model_ , Ann. Probab. 3, 643 (1975).
  * Suchecki  _et al._ (2005) K. Suchecki, V. M. Eguíluz, and M. S. Miguel,  _Voter model dynamics in complex networks: Role of dimensionality, disorder, and degree distribution_ , Phys. Rev. E 72, 036132 (2005).
  * Li  _et al._ (2019) L. Li, Y. Fan, A. Zeng, and Z. Di,  _Binary opinion dynamics on signed networks based on Ising model_ , Physica A 525, 433 (2019).
  * Schulze (2005) C. Schulze,  _Potts-like model for ghetto formation in multi-cultural societies_ , Int. J. Mod. Phys. C 16, 351 (2005).
  * Deffuant  _et al._ (2000) G. Deffuant, D. Neau, F. Amblard, and G. Weisbuch,  _Mixing beliefs among interacting agents_ , Adv. in Compl. Syst. 3, 87 (2000).
  * Weisbuch (2004) G. Weisbuch,  _Bounded confidence and social networks_ , Eur. Phys. J. B 38, 339 (2004).
  * Sobkowicz (2015) P. Sobkowicz,  _Extremism without extremists: Deffuant model with emotions_ , Front. Phys. 3, 17 (2015).
  * Vazquez and Redner (2004) F. Vazquez and S. Redner,  _Ultimate fate of constrained voters_ , J. Phys. A 37, 8479 (2004).
  * Svenkeson and Swami (2015) A. Svenkeson and A. Swami,  _Reaching consensus by allowing moments of indecision_ , Sci. Rep. 5, 14839 (2015).
  * Blume  _et al._ (1971) M. Blume, V. J. Emery, and R. B. Griffiths,  _Ising model for the transition and phase separation in He 3-He4 mixtures_, Phys. Rev. A 4, 1071 (1971).
  * Yang (2010) Y.-H. Yang,  _Blume-Emery-Griffiths dynamics in social networks_ , Phys. Procedia 3, 1839 (2010).
  * Fernandez  _et al._ (2016) M. A. Fernandez, E. Korutcheva, and F. J. de la Rubia,  _A 3-states magnetic model of binary decisions in sociophysics_ , Physica A 462, 603 (2016).
  * Souma  _et al._ (2003) W. Souma, Y. Fujiwara, and H. Aoyama,  _Complex networks and economics_ , Physica A 324, 396 (2003).
  * Li  _et al._ (2003) X. Li, Y. Y. Jin, and G. Chen,  _Complexity and synchronization of the world trade web_ , Physica A 328, 287 (2003).
  * Silva and Matsushita (2022) S. D. Silva and R. Matsushita,  _Editorial: Granularity in econophysics and macroeconomics_ , Front. Phys. 10 (2022).
  * Lu  _et al._ (2015) D. Lu, Y. Tian, V. Y. Liu, and Y. Zhang,  _The performance of the smart cities in China—-A comparative study by means of self-organizing maps and social networks analysis_ , Sustainability 7, 7604 (2015).
  * AlSonosy  _et al._ (2018) O. AlSonosy, S. Rady, N. Badr, and M. Hashem,  _Information innovation technology in smart cities_, edited by L. Ismail and L. Zhang (Springer Singapore, Singapore, 2018) pp. 105–122.
  * Leone  _et al._ (2002) M. Leone, A. Vázquez, A. Vespignani, and R. Zecchina,  _Ferromagnetic ordering in graphs with arbitrary degree distribution_ , Euro. Phys. J. B 28, 191 (2002).
  * De Martino  _et al._ (2012) D. De Martino, S. Bradde, L. Dall’Asta, and M. Marsili,  _Topology-induced inverse phase transitions_ , Europhysics Letters 98, 40004 (2012).
  * Balenzuela  _et al._ (2015) P. Balenzuela, J. P. Pinasco, and V. Semeshenko,  _The undecided have the key: Interaction-driven opinion dynamics in a three-state model_ , PLOS ONE 10, 1 (2015).
  * Newman  _et al._ (2002) M. E. J. Newman, D. J. Watts, and S. H. Strogatz,  _Random graph models of social networks_ , Proc. Natl. Acad. Sci. U.S.A. 99, 2566 (2002).
  * Bianconi (2002) G. Bianconi,  _Mean field solution of the Ising model on a Barabási-Albert network_ , Phys. Lett. A 303, 166 (2002).
  * Dorogovtsev  _et al._ (2008) S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes,  _Critical phenomena in complex networks_ , Rev. Mod. Phys. 80, 1275 (2008).
  * Erdös and Rényi (1959) P. Erdös and A. Rényi,  _On random graphs I_ , Publ. Math. Debr. 6, 290 (1959).
  * Baek  _et al._ (2012) Y. Baek, M. Ha, and H. Jeong,  _Absorbing states of zero-temperature Glauber dynamics in random networks_ , Phys. Rev. E 85, 031123 (2012).
  * Häggström (2002) O. Häggström,  _Zero-temperature dynamics for the ferromagnetic Ising model on random graphs_ , Physica A 310, 275 (2002).

##  Appendix A Stochastic dynamics

The following sequences decribe one elementary move of the stochastic dynamics
used for our simulations.

Metropolis dynamics

Step 1: Pick one agent $i$ at random with uniform probability among the $N$
agents.

Step 2: Propose a random opinion change with probability 1/2 for the agent $i$
and calculate the energy difference $\Delta H$ between the new and old
configuration.

Step 3: Accept the change with probability $\min\\{1,e^{-\beta\Delta H}\\}$,
otherwise remain in the current state.

Glauber dynamics

Step 1: Pick one agent $i$ at random with uniform probability among the $N$
agents.

Step 2: Calculate the energy differences $\Delta H_{1}$ and $\Delta H_{2}$
between the new and old configuration, for the two possible opinion changes
for the agent $i$.

Step 3: Accept either of the two possible changes with probability
$P_{k}=w_{k}/(1+w_{1}+w_{2})$, where $w_{k}=e^{-\beta\Delta H_{k}}$ for
$k=1,2$, or remain in the current state with probability $1-P_{1}-P_{2}$.

##  Appendix B Mean-field solution

We derive here the mean-field solution of the model, which is exact for the FC
case in the large $N$ limit. Consider first a generic graph in which node $i$
has degree $k_{i}$. The mean-field approximation consists in assuming that the
spins $\sigma_{i}$ are independent random variables that can take the values
$\sigma_{i}=\pm 1$ with probabilities $p_{i}^{\pm}$ and $\sigma_{i}=0$ with
probability $1-p_{i}^{+}-p_{i}^{-}$. The Hamiltonian in Eq.(2) takes then the
mean-field form

$H_{MF}=-\sum_{\langle i,j\rangle}[m_{i}m_{j}+\alpha^{2}(1-n_{i})(1-n_{j})]$
(B-20)

where $m_{i}=p_{i}^{+}-p_{i}^{-}$ are the local magnetizations and
$n_{i}=p_{i}^{+}+p_{i}^{-}$ is the probability that $\sigma_{i}$ is different
than zero. The mean-field entropy is given by

$\displaystyle S_{MF}$ $\displaystyle=$
$\displaystyle-\sum_{i}\left[p_{i}^{+}\ln p_{i}^{+}+p_{i}^{-}\ln
p_{i}^{-}+(1-n_{i})\ln(1-n_{i})\right]$ (B-21) $\displaystyle=$
$\displaystyle-\sum_{i}\left[\frac{n_{i}+m_{i}}{2}\ln\frac{n_{i}+m_{i}}{2}+\frac{n_{i}-m_{i}}{2}\ln\frac{n_{i}-m_{i}}{2}+(1-n_{i})\ln(1-n_{i})\right].$

The equilibrium values of $m=\\{m_{i}\\},n=\\{n_{i}\\}$ are those that
minimize the free-energy function ${\cal
L}(m,n,\beta)\equiv(H_{MF}-\beta^{-1}S_{MF})/N$ and thus satisfy the system of
$2N$ coupled equations

$\displaystyle\frac{\partial{\cal L}}{\partial m_{i}}$ $\displaystyle=$
$\displaystyle-\sum_{j\in
V(i)}m_{j}+\frac{1}{\beta}\mbox{arctanh}\frac{m_{i}}{n_{i}}=0$ (B-22)
$\displaystyle\frac{\partial{\cal L}}{\partial n_{i}}$ $\displaystyle=$
$\displaystyle\alpha^{2}\sum_{j\in
V(i)}(1-n_{j})+\frac{1}{2\beta}\ln\frac{n_{i}^{2}-m_{i}^{2}}{4(1-n_{i})^{2}}=0$
(B-23)

where $V(i)$ denotes the set of nodes connected to node $i$. The above
equations can be written as

$\displaystyle m_{i}$ $\displaystyle=\frac{2e^{\beta\alpha^{2}\sum_{j\in
V(i)}(n_{j}-1)}\sinh(\beta\sum_{j\in
V(i)}m_{j})}{1+2e^{\beta\alpha^{2}\sum_{j\in
V(i)}(n_{i}-1)}\cosh(\beta\sum_{j\in V(i)}m_{i})}\,,$ (B-24a) $\displaystyle
n_{i}$ $\displaystyle=\frac{2e^{\beta\alpha^{2}\sum_{j\in
V(i)}(n_{j}-1)}\cosh(\beta\sum_{j\in
V(i)}m_{j})}{1+2e^{\beta\alpha^{2}\sum_{j\in
V(i)}(n_{i}-1)}\cosh(\beta\sum_{j\in V(i)}m_{i})}\,,$ (B-24b)

which we will refer to as local SCEs.

###  B.1 Uniform degree

In the case in which all nodes have the same degree $k_{i}=z$, they will also
have the same $m_{i}=m$ and $n_{i}=n$. Thus the free-energy function takes the
simple form

$\displaystyle{\cal L}(m,n,\beta)$ $\displaystyle=$
$\displaystyle-\frac{z}{2}(m^{2}+\alpha^{2}n^{2})+z\alpha^{2}n$ (B-25)
$\displaystyle+$
$\displaystyle\frac{1}{\beta}\left[\frac{n+m}{2}\ln\frac{n+m}{2}+\frac{n-m}{2}\ln\frac{n-m}{2}+(1-n)\ln(1-n)\right]$

and the (global) SCEs write (in the rest of this Appendix we redefine $\beta
z$ as $\beta$ for brevity)

$\displaystyle m$ $\displaystyle=\frac{2e^{\beta\alpha^{2}(n-1)}\sinh(\beta
m)}{1+2e^{\beta\alpha^{2}(n-1)}\cosh(\beta m)}\,,$ (B-26a) $\displaystyle n$
$\displaystyle=\frac{2e^{\beta\alpha^{2}(n-1)}\cosh(\beta
m)}{1+2e^{\beta\alpha^{2}(n-1)}\cosh(\beta m)}.$ (B-26b)

If $m\neq 0$, the ratio between these two equations gives $n=m\coth(\beta m)$
and, substituting this into Eq.(B-26b), we obtain a self-consistency equation
for $m$ only:

$m=\frac{2e^{\beta\alpha^{2}[m\coth(\beta m)-1]}\sinh(\beta
m)}{1+2e^{\beta\alpha^{2}[m\coth(\beta m)-1]}\cosh(\beta m)}\,.$ (B-27)

Expanding the latter around $m=0$ gives

$a_{2}m+2a_{4}m^{3}+3a_{6}m^{4}+\dots=0$ (B-28)

where

$a_{2}=\frac{1}{2}\left(1-\frac{2\beta}{2+e^{\alpha^{2}(\beta-1)}}\right),\quad
a_{4}=\frac{\beta^{3}}{12}\frac{\left[4-(1+2\alpha^{2})e^{{\alpha^{2}(\beta-1)}}\right]}{(2+e^{\alpha^{2}(\beta-1)})^{2}},$
(B-29)

and $a_{6}>0$. We thus have a tricritical point at the values
$\alpha=\alpha_{tc}$, $\beta=\beta_{tc}$ that are solutions of the two
equations $a_{2}=a_{4}=0$ simultaneously, which gives

$\ln[2(\beta_{tc}-1)]=\frac{3-\beta_{tc}}{2},$ (B-30)

the numerical solution of which is $\beta_{tc}=1.87676\dots$, and

$\alpha_{tc}=\left(\frac{3-\beta_{tc}}{2(\beta_{tc}-1)}\right)^{1/2}=0.800354\dots\,.$
(B-31)

###  B.2 Continuous transition

For $\alpha<\alpha_{tc}$, the critical inverse temperature $\beta_{c}(\alpha)$
is determined by the condition $a_{2}=0$, $a_{4}>0$ which gives

$\alpha^{2}=\frac{1}{\beta_{c}(\alpha)-1}\ln\left[2(\beta_{c}(\alpha)-1)\right]\,.$
(B-32)

The corresponding phase boundary is represented in Fig. 2 of the main text. It
is easy to check that $a_{4}>0$ in an interval of $\beta$ from zero to a value
larger than $\beta_{c}(\alpha)$.

For $\beta$ near $\beta_{c}(\alpha)$ the spontaneous magnetization goes to
zero with the usual mean-field exponent 1/2, i.e.

$m=\left(\frac{-a_{2}}{2a_{4}}\right)^{1/2}\underset{\beta\to\beta_{c}(\alpha)^{+}}{\sim}c(\alpha)\left(\frac{\beta}{\beta_{c}(\alpha)}-1\right)^{1/2}$
(B-33)

where

$c(\alpha)=\left(\frac{-6.0}{\beta_{c}(\alpha)}\frac{(\beta_{c}(\alpha)-1)\alpha^{2}-1}{2-(2\alpha^{2}+1)[\beta_{c}(\alpha)-1]}\right)^{1/2}\,.$
(B-34)

At $\alpha=\alpha_{tc}$ the magnetization vanishes with the usual mean-field
exponent 1/4:

$m=\left(-\frac{a_{2}}{3a_{6}}\right)^{1/4}\underset{\beta\to\beta_{tc}^{+}}{\sim}c(\alpha_{tc})\left(\frac{\beta}{\beta_{c}(\alpha)}-1\right)^{1/4}\,$
(B-35)

with $c(\alpha_{tc})=1.03188$.

The fraction of non-neutral agents, $n=1-n_{0}$, is given by $n=m\coth(\beta
m)$ in the low temperature phase, where $m$ is the solution of Eq.(B-27). In
the high temperature phase, $n$ is found by solving numerically Eq.(B-26b) for
$m=0$, i.e.

$n=\frac{2e^{\beta\alpha^{2}(n-1)}}{1+2e^{\beta\alpha^{2}(n-1)}}.$ (B-36)

The two solutions merge at $\beta_{c}(\alpha)$ and we have
$\lim_{\beta\to\beta_{c}(\alpha)}n=\frac{1}{\beta_{c}(\alpha)}$ from both
sides.

###  B.3 Discontinuous transition

For $\alpha>\alpha_{tc}$ there is a discontinuous transition at an inverse
temperature $\beta_{d}(\alpha)$. The location of the transition is determined
by finding the values of $m,n,\beta_{d}$ that satisfy simultaneously the two
SCEs Eqs.(B-26a) and (B-26b) together with the condition that the free energy
of the ferromagnetic and paramagnetic phases are equal, i.e. ${\cal
L}(m,n,\beta_{d})={\cal L}(0,n(m=0,\beta_{d}),\beta_{d})$, where
$n(m=0,\beta_{d})$ is the solution of Eq.(B-36). In practice we solved
numerically Eq.(B-27) starting from large $\beta$ and reduced $\beta$ in small
increments, following the solution until the above condition was met. This
also allows us to determine the limit of metastability of the ferromagnetic
phase, namely the value of $\beta$ below which there is no longer a non-zero
solution of Eq.(B-27).

##  Appendix C Annealed mean-field approximation

We now return to the case of non-uniform degrees $k_{i}$. In this case, in
principle one could solve numerically the $2N$ local SCEs. Alternatively, one
can make the additional approximation consisting in treating the graph in an
“annealed” fashion, by replacing the adjacency matrix $\epsilon_{i,j}$ (which
is one if $i,j$ are connected and zero otherwise) by $k_{i}k_{j}/zN$, where
$z=\langle k_{i}\rangle$ Bianconi (2002); Dorogovtsev  _et al._ (2008). If we
introduce the weighted order parameters

$m_{w}=\frac{1}{zN}\sum_{j}k_{j}m_{j},\quad
n_{w}=\frac{1}{zN}\sum_{j}k_{j}n_{j},$ (B-37)

we have

$\sum_{j\in V(i)}m_{j}=\sum_{j\neq
i}\epsilon_{i,j}m_{j}=\frac{1}{zN}\sum_{j\neq i}k_{i}k_{j}m_{j}=k_{i}m_{w},$
(B-38)

and similarly $\sum_{j\in V(i)}n_{j}=k_{i}n_{w}$. Multiplying the local SCEs
by $k_{i}/(zN)$ and summing over $i$ we thus obtain two SCEs for the weighted
order parameters

$\displaystyle m_{w}$ $\displaystyle=$
$\displaystyle\frac{1}{zN}\sum_{i}k_{i}\frac{2e^{\beta\alpha^{2}k_{i}(n_{w}-1)}\sinh(\beta
k_{i}m_{w})}{1+2e^{\beta\alpha^{2}k_{i}(n_{w}-1)}\cosh(\beta k_{i}m_{w})}$
(B-39) $\displaystyle n_{w}$ $\displaystyle=$
$\displaystyle\frac{1}{zN}\sum_{i}k_{i}\frac{2e^{\beta\alpha^{2}k_{i}(n_{w}-1)}\cosh(\beta
k_{i}m_{w})}{1+2e^{\beta\alpha^{2}k_{i}(n_{w}-1)}\cosh(\beta k_{i}m_{w})}$
(B-40)

For large $N$ we can replace the sum over $i$ with a sum over all possible
degrees,

$\displaystyle m_{w}$ $\displaystyle=$
$\displaystyle\frac{1}{z}\sum_{k}kP(k)\frac{2e^{\beta\alpha^{2}k(n_{w}-1)}\sinh(\beta
km_{w})}{1+2e^{\beta\alpha^{2}k(n_{w}-1)}\cosh(\beta km_{w})}$ (B-41)
$\displaystyle n_{w}$ $\displaystyle=$
$\displaystyle\frac{1}{z}\sum_{k}kP(k)\frac{2e^{\beta\alpha^{2}k(n_{w}-1)}\cosh(\beta
km_{w})}{1+2e^{\beta\alpha^{2}k(n_{w}-1)}\cosh(\beta km_{w})}.$ (B-42)

Finally, after solving the above equations we can obtain the average order
parameters as

$\displaystyle m$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i}m_{i}=\sum_{k}P(k)\frac{2e^{\beta\alpha^{2}k(n_{w}-1)}\sinh(\beta
km_{w})}{1+2e^{\beta\alpha^{2}k(n_{w}-1)}\cosh(\beta km_{w})}$ (B-43)
$\displaystyle n$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{i}n_{i}=\sum_{k}P(k)\frac{2e^{\beta\alpha^{2}k(n_{w}-1)}\cosh(\beta
km_{w})}{1+2e^{\beta\alpha^{2}k(n_{w}-1)}\cosh(\beta km_{w})}.$ (B-44)

##  Appendix D Exact solution for the one dimensional lattice

We solve the model exactly on a 1D lattice with periodic boundary conditions
using the transfer-matrix method. By ordering the states as $1,0,-1$, the
transfer matrix
$\mathcal{T}_{\sigma_{i},\sigma_{i+1}}=\exp\\{\beta\left[\sigma_{i}\sigma_{j}+\alpha^{2}(\sigma^{2}_{i}-1)(\sigma^{2}_{j}-1)\right]\\}$
writes

$\mathcal{T}=\begin{bmatrix}w&1&w^{-1}\\\ 1&v&1\\\ w^{-1}&1&w\end{bmatrix},$
(C-45)

where we defined $w\equiv e^{\beta},v\equiv e^{\beta\alpha^{2}}$. The free
energy is given by
$F(\beta,\alpha,N)=-\frac{1}{\beta}\ln{\mbox{tr}}\,\mathcal{T}=-\frac{1}{\beta}(\lambda_{+}^{N}+\lambda_{0}^{N}+\lambda_{-}^{N})$,
where the eigenvalues of $\mathcal{T}$ are

$\displaystyle\lambda_{0}$ $\displaystyle=$ $\displaystyle w-w^{-1}$
$\displaystyle\lambda_{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[w+w^{-1}+v\right]\pm\frac{1}{2}\left[(w+w^{-1}-v)^{2}+8\right]^{1/2}$
(C-46)

In the $N\to\infty$ limit the largest eigenvalue $\lambda_{+}$ dominates the
sum, thus the free energy per spin is

$f(\beta,\alpha)=\lim_{N\to\infty}\frac{1}{N}F(\beta,N)=-\frac{1}{\beta}\ln\left[w+w^{-1}+v+\sqrt{(w+w^{-1}-v)^{2}+8}\right]\,.$
(C-47)

Substituting $w$ and $v$ in the expression above, we obtain Eq.(12) of the
main text. The internal energy per spin is given by

$u(\beta,\alpha)=\frac{\partial(\beta
f)}{\partial\beta}=-\frac{w-w^{-1}+\alpha^{2}v+B^{-1/2}(w-w^{-1}-\alpha^{2}v)(w+w^{-1}-v)}{w+w^{-1}+v+B^{-1/2}}$
(C-48)

where $B\equiv(w+w^{-1}-v)^{2}+8$. We also use the transfer matrix method to
compute the spin correlation function

$C(r)=\langle\sigma_{i}\sigma_{i+r}\rangle=\frac{{\mbox{tr}}(\Sigma\,\mathcal{T}^{r}\,\Sigma\,\mathcal{T}^{N-r})}{{\mbox{tr}}\mathcal{T}^{N}}$
(C-49)

where $\Sigma$ is a diagonal matrix with elements
$\Sigma_{\sigma,\sigma^{\prime}}=\sigma\delta_{\sigma,\sigma^{\prime}}$. We
can write the transfer matrix as
$\mathcal{T}=\sum_{\mu}\lambda_{\mu}|\mu\rangle\langle\mu|$ where $\mu=0,\pm
1$ and $|\mu\rangle$ are the eigenvectors of $\mathcal{T}$:

$|0\rangle=\frac{1}{\sqrt{2}}\begin{bmatrix}1\\\ 0\\\
-1\end{bmatrix},\quad|\pm\rangle=\frac{1}{\sqrt{1+y_{\pm}^{2}}}\begin{bmatrix}1\\\
y_{\pm}\\\ 1\end{bmatrix},$ (C-50)

where $y_{\pm}=\lambda_{\pm}-w-w^{-1}$. In this way we obtain

$C(r)=\frac{1}{\sum_{\mu}\lambda_{\mu}^{N}}\left\\{\left[\lambda_{0}^{N}\left(\frac{\lambda_{+}}{\lambda_{0}}\right)^{r}+\lambda_{+}^{N}\left(\frac{\lambda_{0}}{\lambda_{+}}\right)^{r}\right]|\langle+|\Sigma|0\rangle|^{2}+\left[\lambda_{0}^{N}\left(\frac{\lambda_{-}}{\lambda_{0}}\right)^{r}+\lambda_{-}^{N}\left(\frac{\lambda_{0}}{\lambda_{-}}\right)^{r}\right]|\langle-|\Sigma|0\rangle|^{2}\right\\}$
(C-51)

which, since $|\lambda_{+}|>|\lambda_{0}|>|\lambda_{-}|$, for $N\to\infty$
gives

$C(r)=\frac{4}{2+y_{+}^{2}}e^{-r/\xi},\quad\xi\equiv\frac{1}{\ln\frac{\lambda_{+}}{\lambda_{0}}},$
(C-52)

where we used $|\langle+|\Sigma|0\rangle|^{2}=4/(2+y_{+}^{2})$. Inserting
Eqs.(C-46) into the expression above for $\xi$, and substituting $w$ and $v$,
we obtain Eq.(V) of the main text.

##  Appendix E Residual Energy for the Erdös-Rényi graphs

In this section we present additional zero-temperature MC results for ER
graphs with low connectivity, in order to explore the effects of dynamical
traps that prevent the system from reaching optimality. We considered both the
Metropolis and Glauber dyamics, and considering that the system has reached a
steady state when the energy does not change for $100\log(N)$ MCS. with the
same energy, then calculating the difference between the measured energy and
the ground state. Fig. E.1 shows the distribution of the residual energy
$\epsilon=(E_{steady}-E_{0})/Nz$. We can observe that both algorithms give
almost identical results.

If we compare the case $\alpha=0$ with the Ising model Baek  _et al._ (2012)
we appreciate that the peaks are lower in our model, due to the existence of
the neutral opinion, which reduces the amount of energy necessary to
transition between states $+1$ and $-1$. For the case $\langle k\rangle=4$,
which presents a behavior very similar to the complete graph, deviations from
the ground state are only relevant when we are close to the mean-field
tricritical point, pointing to the existence of a first order transition also
in ER graphs.

Figure E.1: Probability distribution of the residual energy for different
values of $\alpha$ and $\langle k\rangle=2,3,4$, using Metropolis (top panel)
and Glauber dynamics (bottom panel). Results are obtained for a single
realization of a network with $N=500$, and $10^{5}$ repetitions of the
simulation.

In Fig. E.2 we can compare the probability distribution of the residual energy
for two ER networks of different sizes and $\langle k\rangle=3$ using Glauber
dynamics. The distribution is bimodal for values of $\alpha\leq 0.75$, but
becomes unimodal for larger values of $\alpha$, which is related to the
symmetry breaking caused by the neutral opinion above this threshold.

Figure E.2: Probability distribution of the residual energy for ER graphs with
$\langle k\rangle=3$ and sizes $N=500$ and $N=1000$.