# Biased-voter model: how persuasive a small group can be?

Agnieszka Czaplicka(1), Christos Charalambous(2), Raul Toral(2), Maxi San
Miguel(2) (1) Center for Humans and Machines, Max Planck Institute for Human
Development, Lentzeallee 94, Berlin 14195, Germany.
(2) Institute for Cross-Disciplinary Physics and Complex Systems IFISC (CSIC-
UIB), Campus UIB, 07122 Palma de Mallorca, Spain

###### Abstract

We study the voter model dynamics in the presence of confidence and bias. We
assume two types of voters. Unbiased voters whose confidence is indifferent to
the state of the voter and biased voters whose confidence is biased towards a
common fixed preferred state. We study the problem analytically on the
complete graph using mean field theory and on an Erdős-Rényi random network
topology using the pair approximation, where we assume that the network of
interactions topology is independent of the type of voters. We find that for
the case of a random initial setup, and for sufficiently large number of
voters $N$, the time to consensus increases proportionally to $\log(N)/\gamma
v$, with $\gamma$ the fraction of biased voters and $v$ the parameter
quantifying the bias of the voters ($v=0$ no bias). We verify our analytical
results through numerical simulations. We study this model on a biased-
dependent topology of the network of interactions and examine two distinct,
global average-degree preserving strategies (model I and model II) to obtain
such biased-dependent random topologies starting from the biased-independent
random topology case as the initial setup. Keeping all other parameters
constant, in model I, $\mu_{\text{\tiny{BU}}}$, the average number of links
among biased (B) and unbiased (U) voters is varied at the expense of
$\mu_{\text{\tiny{UU}}}$ and $\mu_{\text{\tiny{BB}}}$, i.e. the average number
of links among only unbiased and biased voters respectively. In model II,
$\mu_{\text{\tiny{BU}}}$ is kept constant, while $\mu_{\text{\tiny{BB}}}$ is
varied at the expense of $\mu_{\text{\tiny{UU}}}$. We find that if the agents
follow the strategy described by model II, they can achieve a significant
reduction in the time to reach consensus as well as an increment in the
probability to reach consensus to the preferred state. Hence, persuasiveness
of the biased group depends on how well its members are connected among each
other, compared to how well the members of the unbiased group are connected
among each other.

###### keywords:

Biased voter model, Biased-dependent complex network topologies

##  1 Introduction

The process by which people adopt an opinion about a given issue, such as
endorsing a political option or choosing a commercial product, is a complex
social phenomenon, and often the underlying mechanisms driving opinion
dynamics are not well understood. Yet, public opinion is today a key player in
most issues faced by our societies and policy makers are obliged to take into
account the evolution of public opinion. Broadly speaking, the potential
influences on an individual’s opinion can be divided in three categories:
those that are intrinsic to the individual such as, for example, personal
beliefs; global external factors such as mass media and, finally, interactions
with other members of the society. In order to capture the opinion changes
subject to this latter source of influence, and inspired by the idea of
understanding macroscopic behavior emerging from simple interactions between
particles, a plethora of models for opinion dynamics has emerged in the
statistical physics literature [1, 2, 3]. Among them, one of the simplest and
most extensively studied is the so-called voter model [4], a model of opinion
dynamics that leads to herding. It was independently proposed in various
research fields to study, among other things, neutral genetic drift in an
ideal population [5, 6], competition for territory between two countries [7],
spreading of infectious diseases [8], language competition [9, 10, 11, 12],
kinetics of catalytic reactions [13, 14, 15], coarsening phenomena [15, 16],
opinion dynamics [17], political elections [18], etc. One of the
attractiveness of the voter model is that it is one of the few known
interacting-particle models exactly solvable in regular lattices of any
spatial dimension $D$ [19, 20]. Directly, or in many variations, it can also
be related to other standard well studied models in mathematics and physics
such as coalescing random walkers [21], the zero-temperature Glauber kinetic
Ising model and the linear Glauber model [22].

In its original formulation [7, 21], the voter model was introduced as an
Ising-like model where an individual (an “agent” or “voter”) associated with a
lattice site $i$ can adopt two different values or “opinions” $s_{i}=\pm 1$.
The dynamics of the system is implemented by randomly choosing one individual
and assigning to it the value of the opinion of one of its randomly chosen
nearest neighbors. The voter model is characterized by purely noise-driven
diffusive dynamics. It exhibits two symmetric absorbing states, called
“consensus” states, from which the system cannot escape and which, for finite-
size systems, are reached almost surely. If $N$ is the total number of voters,
the mean time to reach consensus, $T_{N}$, scales in regular lattices as
$N^{2}$ in $D=1$, as $N\log N$ in $D=2$, and as $N$ in $D=3$ [23]. Many social
systems display interactions that find a better characterization as complex
networks with distinctive connectivity properties [24, 25]. For this reason,
in the last decades, an extensive effort has been devoted to studying voter-
like models on complex networks [26, 27, 28, 29, 30, 31, 32, 33]. In this
scenario, defining $\mu_{k}$ as the $k$-th moment of the degree distribution,
it is found that for uncorrelated networks, $T_{N}\sim N\mu_{1}^{2}/\mu_{2}$,
which grows sublinearly in $N$ for a sufficiently broad degree distribution
[30, 31].

Although the early versions of the model consider that all agents are
identical, it is obvious that in real applications there will be structural
differences between the agents. For instance, the number of nearest neighbors,
the rate of interactions, the preference for one or another state, etc. can
broadly vary. These inhomogeneities, modeled as quenched disorder, are known
to have an important relevance in non-equilibrium systems with absorbing
states [34, 35, 36, 37, 38]. Many variations of the voter model with quenched
disorder have been proposed [39], such as the inclusion of contrarians,
defined as agents who adopt a different opinion than that of its neighbors
[40, 41], or zealots, defined as agents who favor [42, 43] or maintain
inflexibly [44, 45] a fixed opinion. This favoring is implemented by including
in the dynamics of the zealots an spontaneous transition rate or “noise”,
independent of the copying mechanism, from the disfavored to the favored
state.

Other models include the preference or bias for one of the states in the
copying mechanism. This is the case of the partisan voter model [46] in which
the population is split into two groups, each one favoring one of the options
(Democrats or Republicans in their example), or a model for language
competition [10] in which all agents favor one of the two possible choices,
understood then as a difference of prestige between the languages. This
particular copying mechanism with preference turns out to be isomorphic to a
model of reaction-limited heterogeneous catalysis [47, 48] where the bias is
the difference in the probabilities of attempting an adsorption of one of two
reactive molecular species onto an empty substrate site.

In this work, we consider a variant of the voter model in which a fraction of
the population is biased towards one of the two options, while the rest of the
population is neutral. Our intention is to determine if the biased community
can optimize in some way its connectivity in order to have a maximum influence
in the behavior of the whole system. Previous studies on an unbiased majority-
like voter model have shown that the opinion held by a minority group can win
over that of a larger one provided that it has more internal cohesion
(stronger or more connections) than the majority group [49]. Our setup is very
different from previous models with bias in the copying mechanism, such as
that of the partisan voter model where each agent displays a bias towards one
or another option, or the language model in which all agents favor the same
option.

The structure of the paper is as follows: In Sec. 2 we present the details of
the model. In Sec 3, we study the model dynamics on a complete graph using
mean-field theory. In Sec. 4, we extend the study considering the voters’
dynamics on an Erdős-Rényi random topology, where we assume that the
probability of a connection between to sites is independent on whether the
voters are biased or not. Finally in Sec. 5, we extend our studies to the case
where voters lie on two Erdős-Rényi (ER) networks of distinct characteristics,
i.e. distinct average degrees, depending on their biased/unbiased type. For
this latter scenario, we identify a strategy by which one could change the
underlying topology of the network, while preserving the total average degree,
in such a way as to achieve a reduction of the consensus time. We show that
this is the case when the average degree of the network of interactions among
only biased individuals is increased at the cost of decreasing the average
degree of the network of interactions of only unbiased voters, while
maintaining the average degree of their in-between interactions constant.

##  2 The model

Let us consider a given lattice of $N$ nodes connected by links. The network
is single connected and can not be split in two disjoint ones. Each node
$i=1,\dots,N$ represents an agent that holds a binary state (opinion) variable
$s_{i}=\pm 1$. In the standard voter model [7, 21] those state variables
evolve by an interaction mechanism by which agents copy the opinion of a
randomly selected neighbor (those located on connected nodes). We modify these
rules by introducing a group of biased agents with preference for one of the
states. We consider that $N_{\text{\tiny{B}}}=\gamma N$ agents are biased, and
the remaining $N_{\text{\tiny{U}}}=(1-\gamma)N$ are unbiased. Bias is
introduced as a parameter $v\in[-1,1]$ which alters the probabilities of
biased agents to copy the state of a neighbor. As in the standard voter model
dynamics, we select one agent at random (node $i$). Next, one of its neighbors
(node $j$) is selected also randomly. If $s_{i}=s_{j}$, nothing happens. Else,
depending on whether $i$ is biased or unbiased and the state of agent $j$, the
following scenarios are considered:

  * 1.

If $i$ is unbiased then $i$ copies $j$’s state with probability $1/2$.

  * 2.

If $i$ is biased:

    * (a)

if $s_{i}=+1$ and $s_{j}=-1$ then $i$ copies $j$’s state with probability
$\dfrac{1-v}{2}$;

    * (b)

if $s_{i}=-1$ and $s_{j}=+1$ then $i$ copies $j$’s state with probability
$\dfrac{1+v}{2}$.

$N$ of these node selections constitute one Monte Carlo step. Under these
rules, the preferred state is $s_{i}=+1$ (resp. $s_{i}=-1$), $\forall i$ if
$v>0$ (resp. $v<0$). The case $v=\pm 1$ of extreme preference for one of the
two states will not be considered here as it bears some similarity with that
of zealot agents, those that never change their opinion, considered elsewhere
[50]. Note that the probability for a biased node to copy the state of a
neighbor is independent of whether that neighbor is biased or unbiased. The
usual voter model is recovered either for $\gamma=0$ (no presence of biased
agents) of for $v=0$ when all the nodes follow the standard voter model
dynamics with the modification that a neighbor’s state is copied only with
probability $1/2$. This modification is irrelevant as its only effect is to
rescale time (as measured by the number of Monte Carlo steps) by a factor of
$2$. As a result of the dynamical rules the system might enter an absorbing
configuration where no agent can change its state and no further evolution is
possible. With the rules considered here, the only possible absorbing
configurations are consensus situations, where all agents hold the same state,
either $+1$ or $-1$. The standard voter model also accounts in some situations
for dynamical steady states where a macroscopic fraction of agents (but not
all of them) hold a particular opinion during a long period of time until a
finite size fluctuation takes the system to one of the absorbing
configurations. In the biased voter model the symmetry between $+1$ and $-1$
sattes is broken, so that the ensemble average $\langle s_{i}\rangle$ is no
longer a conserved quantity and the system tends to reach the preferred
absorbing configuration by its intrinsic dynamics. Still, finite size
fluctuations can lead the system to the non-preferred absorbing configuration.

##  3 Mean-field approximation

We present now an analytical treatment of the biased-voter model based on a
mean-field type approximation valid for an all-to-all (or complete graph)
configuration, where each agent is connected to all other agents. Let
$\sigma=\frac{1}{2N}\sum_{i=1}^{N}(s_{i}+1)$ be the fraction of network nodes
in state $+1$. One can treat the problem as that of a random walk where the
variable $\sigma\in[0,1]$ can increase or decrease due to the dynamics. At
each time step a randomly selected node $i$ can change its state depending on
the state of the randomly selected neighbor $j$. If the change does occur,
then the fraction $\sigma$ increases or decreases by an amount
$\Delta_{\sigma}=1/N$. We denote by $R^{+}(\sigma)$ and $R^{-}(\sigma)$,
respectively, the transition probabilities
$\text{Prob}[\sigma\rightarrow\sigma+\Delta_{\sigma}]$,
$\text{Prob}[\sigma\rightarrow\sigma-\Delta_{\sigma}]$. Using the rules of the
process described in Section 2 we can write:

$\displaystyle R^{+}(\sigma)$ $\displaystyle=$
$\displaystyle\text{Prob}(s_{i}=-1,s_{j}=1)\left[\text{Prob}(i\text{ no
bias})\frac{1}{2}+\text{Prob}(i\text{ bias})\frac{1+v}{2}\right]$ (1)
$\displaystyle=$
$\displaystyle\text{Prob}(s_{j}=1|s_{i}=-1)\text{Prob}(s_{i}=-1)\left[(1-\gamma)\frac{1}{2}+\gamma\frac{1+v}{2}\right]$
$\displaystyle=$ $\displaystyle\sigma(1-\sigma)\frac{1+\gamma v}{2}.$

Here we have used $\text{Prob}(s_{i}=-1)=1-\sigma$ and the mean-field
approximation: $\text{Prob}(s_{j}=+1|s_{i}=-1)=\text{Prob}(s_{j}=+1)=\sigma$,
equivalent to assuming that the state of an agent is independent of that of
its neighbors. This implies that the density $\rho$ of active links, those
connecting nodes with different opinions, reads as
$\rho=\text{Prob}(s_{i}=-s_{j})=2\sigma(1-\sigma)$. Furthermore, we have
assumed that the label of “biased” or “unbiased” of agent $i$ is independent
on its state value $s_{i}$.

Similarly, we derive

$\displaystyle R^{-}(\sigma)$ $\displaystyle=$
$\displaystyle\text{Prob}(s_{i}=1,s_{j}=-1)\left[\text{Prob}(i\text{ no
bias})\frac{1}{2}+\text{Prob}(i\text{ bias})\frac{1-v}{2}\right]$ (2)
$\displaystyle=$
$\displaystyle\text{Prob}(s_{j}=-1|s_{i}=1)\text{Prob}(s_{i}=+1)\left[(1-\gamma)\frac{1}{2}+\gamma\frac{1-v}{2}\right]$
$\displaystyle=$ $\displaystyle\sigma(1-\sigma)\frac{1-\gamma v}{2}.$

Once these transition probabilities $R^{+}(\sigma),\,R^{-}(\sigma)$ have been
derived, one can use the standard machinery of random walk theory [51, 52] to
compute several quantities of interest. We focus in this paper on the fixation
probability $P_{1}$ (probability that all agents eventually reach consensus on
the state $+1$), the average time $\tau$ to reach any absorbing state and the
average times $\tau_{1},\,\tau_{-1}$ to reach consensus on states $+1$ and
$-1$, respectively. Without lack of generality, we assume henceforth that
$v>0$ such that the preferred state is $+1$.

###  3.1 Absorbing state

The magnetization $m$ is defined in terms of the fraction $\sigma$ as
$m=2\sigma-1$. In a single time step $\Delta_{t}=1/N$ (in units of Monte Carlo
steps) $m$ can vary by an amount $\pm\Delta_{m}=\pm 2\Delta_{\sigma}=\pm 2/N$.
The evolution equation for the probability $P(m,t)$ of finding a magnetization
$m$ at time $t$ follows from the basic rules of the process as

$\displaystyle P(m,t+\Delta_{t})$ $\displaystyle=$ $\displaystyle
R^{-}(m+\Delta_{m})P(m+\Delta_{m},t)$ $\displaystyle+$ $\displaystyle
R^{+}(m-\Delta_{m})P(m-\Delta_{m},t)$ $\displaystyle+$
$\displaystyle[1-R^{+}(m)-R^{-}(m)]P(m,t).$

where $R^{\pm}(m)$ are the transition probabilities Eqs.(1,2) written in terms
of the $m$ variable. Upon Taylor expanding up to order $\Delta_{t}$ in time
and up to second order in $\Delta_{m}$ or, equivalently, taking the continuous
limit $1/N\rightarrow 0$ and keeping only terms up to order $1/N$, the time
evolution of $P(m,t)$ is given by

$\displaystyle\frac{\partial P(m,t)}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{\Delta_{m}}{\Delta_{t}}\frac{\partial[(R^{-}(m)-R^{+}(m))P(m,t)]}{\partial
m}$ (4)
$\displaystyle+\frac{(\Delta_{m})^{2}}{2\Delta_{t}}\frac{\partial^{2}[(R^{-}(m)+R^{+}(m))P(m,t)]}{\partial
m^{2}},$

which, after replacing the expressions for $R^{\pm}(m)$, $\Delta_{m}$ and
$\Delta_{t}$, becomes

$\displaystyle\frac{\partial P\left(m,t\right)}{\partial t}=-\frac{\gamma
v}{2}\frac{\partial\left[\left(1-m^{2}\right)P\left(m,t\right)\right]}{\partial
m}+\frac{1}{2N}\frac{\partial^{2}\left[\left(1-m^{2}\right)P\left(m,t\right)\right]}{\partial
m^{2}}.$ (5)

For the case $\gamma v=0$ this reproduces the results of the standard voter
model [47, 53] with the already mentioned rescaling of the time by a factor of
$2$.

Eq.(5) is a Fokker-Planck equation with state-dependent drift
$F(m)=\frac{\gamma v}{2}(1-m^{2})$ and diffusion coefficient
$D\left(m\right)=\frac{1-m^{2}}{N}$, hence the evolution of the magnetization
can be viewed as the motion of a random walk moving in a medium that is
increasingly “sticky” near the extremities of the absorbing interval. When
$m=\pm 1$ the walk stops, independently of whether there is bias or not. In
principle, the steady-state solution of the Fokker-Planck Eq.(5) would be
$P_{\text{st}}\left(m\right)=Z^{-1}\frac{e^{\gamma vNm}}{1-m^{2}}$. However,
the normalization constant is $Z=\int_{-1}^{1}\frac{e^{\gamma
vNm}}{1-m^{2}}dm=\infty$, as the integral diverges in both limits $m=\pm 1$.
This indicates that the only absorbing state for a finite size system is
consensus to either of these two states. Given that $\gamma v>0$, we conclude
that in the thermodynamic limit it is $P_{\text{st}}(m)=\delta(m-1)$, being
$\delta(\cdot)$ the Dirac-delta function, and the stationary consensus state
will be the preferred one $m=1$. To verify this, we derive from the Fokker-
Planck Eq.(5) the following equation of motion for the average magnetization
$\langle m\rangle$

$\displaystyle\frac{\partial\langle m\rangle}{\partial t}=\frac{\gamma
v}{2}(1-\langle m^{2}\rangle).$ (6)

Neglecting fluctuations, $\langle m^{2}\rangle\approx\langle m\rangle^{2}$,
the solution for an initial condition $\langle m(0)\rangle=0$ is $\langle
m(t)\rangle=\tanh(\frac{\gamma v}{2}t)$ as shown in [10] for the case
$\gamma=1$. This describes a monotonic evolution to the stationary state in a
characteristic time scale $1/(\gamma v)$. In the following we consider the
effect of finite size fluctuations that can lead the system to the non-
preferred absorbing state.

###  3.2 Fixation probability, $P_{1}$

The fixation (or exit) probability $P_{1}(\sigma)$ is defined as the
probability that a finite system with an initial fraction $\sigma$ reaches a
consensus to the preferred state $+1$ in a finite number of steps [23]. It can
be expressed as the probability of making one of the transitions
$\sigma\to\sigma-\Delta_{\sigma},\,\sigma,\,\sigma+\Delta_{\sigma}$ multiplied
by the exit probability from these intermediate points:

$\displaystyle P_{1}(\sigma)$ $\displaystyle=$ $\displaystyle
R^{+}(\sigma)P_{1}(\sigma+\Delta_{\sigma})+R^{-}(\sigma)P_{1}(\sigma-\Delta_{\sigma})$
(7) $\displaystyle+$
$\displaystyle[1-R^{+}(\sigma)-R^{-}(\sigma)]P_{1}(\sigma),$

with the boundary conditions $P_{1}(0)=0$, $P_{1}(1)=1$. Introducing the
notation $P_{1}(n)=P_{1}(\sigma=\frac{n}{N})$, and using Eqs.(1,2), Eq.(7) is
rewritten as a recurrence equation

$P_{1}(n)=\frac{1+\gamma v}{2}P_{1}(n+1)+\frac{1-\gamma v}{2}P_{1}(n-1).$ (8)

The solution of this equation satisfying the aforementioned boundary
conditions is:

$P_{1}(n)=\frac{1-a^{n}}{1-a^{N}},\quad a=\dfrac{1-\gamma v}{1+\gamma v},$ (9)

which for random initial conditions $n=\frac{N}{2}$ takes the form

$P_{1}(n=N/2)=\frac{1}{1+a^{N/2}}.$ (10)

If we approximate $\ln a=-2\gamma v+O(\gamma v)^{3}$, at this order in the
expansion of the logarithm, we obtain $a^{N/2}\approx\exp(-\gamma vN)$, which
leads to

$P_{1}(n=N/2)=\frac{1}{1+\exp(-\beta/2)},\,\beta=2\gamma vN.$ (11)

Note that the exit probability in this approximation is only a function of the
product $\gamma vN=vN_{\text{\tiny{B}}}$. It turns out that the approximation
is very good and the maximum absolute difference between Eqs.(10) and (11) is
always smaller than $0.42N^{-2}$ for all values of $\gamma v$, such that for
the system size $N=1000$ used in most numerical simulations of this paper, the
error of the approximation is less then $4.2\times 10^{-7}$.

It is worth noting that the approximate solution Eq.(11) can also be obtained
from a continuous version of Eq.(7), obtained by expanding to second order in
$\Delta_{\sigma}$:

$\frac{1}{2}\Delta_{\sigma}^{2}\left(R^{+}+R^{-}\right)\frac{d^{2}P_{1}}{d\sigma^{2}}+\Delta_{\sigma}\left(R^{+}-R^{-}\right)\frac{dP_{1}}{d\sigma}=0,$
(12)

the so-called backward Kolmogorov equation for the exit probability. After
replacing Eqs.(1,2) and $\Delta_{\sigma}=1/N$, we obtain

$\frac{d^{2}P_{1}}{d\sigma^{2}}+\beta\frac{dP_{1}}{d\sigma}=0,$ (13)

whose solution for boundary conditions $P_{1}(0)=0$ and $P_{1}(1)=1$ is 111The
general solution of  $z^{\prime\prime}(\sigma)+\beta
z^{\prime}(\sigma)=g(\sigma)$ is
$z(\sigma)=C_{2}+C_{1}e^{-\beta\sigma}+\int^{\sigma}d\sigma^{\prime}\int^{\sigma^{\prime}}d\sigma^{\prime\prime}e^{-\beta(\sigma^{\prime}-\sigma^{\prime\prime})}g(\sigma^{\prime\prime}),$
where $C_{1},\,C_{2}$ are integration constants found by fulfilling the
adequate boundary conditions.

$P_{1}(\sigma;\beta)=\frac{1-e^{-\sigma\beta}}{1-e^{-\beta}}.$ (14)

Setting the initial fraction $\sigma=\frac{1}{2}$ we recover Eq.(11).
Furthermore, when $\gamma v=0$ we recover the well known formula for the
standard voter model $P_{1}(\sigma)=\sigma$. As shown in Fig. 1 the analytical
expression Eq. (11) agrees well with the data coming from numerical
simulations of the model on a complete graph. The figure also shows results
corresponding to the random (Erdős-Rény) network for different values of the
average connectivity $\mu$ that will be analyzed in detail in Sec.4.

Figure 1: Plot of the fixation probability (probability to reach the preferred
state $s=+1$), $P_{1}(\sigma,\beta)$, as a function of the scaling variable
$\beta=2\gamma vN$ for two different values of the initial condition:
$\sigma=0.5$ (left panel) and $\sigma=0.1$ (right panel), using values
$\gamma=0.1$ and $N=1000$. The continuous (blue) lines are the analytical
result, Eq.(14), and the open square ($\square$) symbols denote the results of
computer simulations in a complete graph (CG) for $v=0.01$. We also plot the
fixation probability for two random Erdős-Rény networks with average
connectivity $\mu=5$ and $\mu=10$ and several values of $v$ as indicated in
the legend of the left panel, and $v=0.01$ in the case of the right panel. The
theoretical predictions lines, dashed (green) for $\mu=10$ and dotted (red)
for $\mu=5$, come from the same expression Eq.(14) replacing $\beta$ by
$\beta_{\mu}=2\gamma v\dfrac{\mu}{\mu+1}N$ as explained in section 4.3, while
the different symbols are the results of numerical simulations.

A consequence of Eq.(14) is that, as the system size increases, the
probability $1-P_{1}(\sigma;\beta)$ of reaching consensus in the non-preferred
opinion decreases and, eventually, tends to zero as $N$ tends to infinity for
any non-zero value of the product $\gamma v$. For a finite system size,
however, there is a finite probability to reach the non-preferred state, which
means that, in principle, it is possible to observe some realizations of the
dynamics leading to consensus for this non-preferred state. A word of caution
is relevant here: due to the smallness of this probability for large system
sizes, a large number $Q$ of realizations is needed to observe a consensus in
the non-preferred state in a numerical simulation of the process.
Alternatively, for a given number or realizations $Q$ there will be a value of
$(\gamma vN)_{0}$ above which the probability to reach the non-preferred state
is smaller that $1/Q$ and no consensus to the non-preferred state will be
likely to be observed in the numerical simulation, leading to the wrong
conclusion that order in this non-preferred state is never possible if $\gamma
vN>(\gamma vN)_{0}$. This value of $(\gamma vN)_{0}$ can be estimated by
noting that the probability that at least one of the $Q$ runs ends in the non-
preferred state ($s_{i}=-1,\,\forall i$, remember that we assume $\gamma v>0$)
is $1-P_{1}^{Q}$ and we demand this probability to be of the order the inverse
of the number of runs $1-P_{1}^{Q}\sim 1/Q$. Using Eq.(11), we arrive at the
condition

$(\gamma vN)_{0}\sim-\ln\left((1-\frac{1}{Q})^{-1/Q}-1\right)=2\ln Q+O(1/Q)$
(15)

If $\gamma vN>(\gamma vN)_{0}$ then no runs will be typically observed to
reach the non-preferred state in the $Q$ runs of the simulation. In practice,
we observe that the threshold value $(\gamma vN)_{0}$ scales roughly as
$1.3\ln Q$, see Fig. 2.

Figure 2: In this figure we plot the threshold value $(\gamma vN)_{0}$ as a
function of the number of runs $Q$ as obtained from numerical simulations of
the biased-voter model in a complete-graph (all-to-all connectivity between
the agents). As explained in the main text, $(\gamma vN)_{0}$ is defined as
the value of $\gamma vN$ above which none of the $Q$ runs is observed to end
in the non-preferred state $s=-1$. The solid line is the approximation
$2\ln(Q)$ as given by Eq.(15), while the dashed line is a fit to the form
$1.3\ln Q$.

###  3.3 Time to reach consensus, $\tau$

Let $\tau(\sigma)$ be the average time to reach any absorbing state (all
agents taking the same value, $s_{i}=s,\,\forall i$ with either $s=+1$ or
$s=-1$) starting from an initial fraction $\sigma$. As before, we consider the
transitions $\sigma\to\sigma-\Delta_{\sigma},\sigma,\sigma+\Delta_{\sigma}$
and relate $\tau(\sigma)$ to the average times from these intermediate points,

$\displaystyle\tau(\sigma)$ $\displaystyle=$ $\displaystyle
R^{+}(\sigma)\left[\tau(\sigma+\Delta_{\sigma})+\Delta_{t}\right]+R^{-}(\sigma)\left[\tau(\sigma-\Delta_{\sigma})+\Delta_{t}\right]$
(16)
$\displaystyle+\left[1-R^{+}(\sigma)-R^{-}(\sigma)\right]\left[\tau(\sigma)+\Delta_{t}\right],$

and boundary conditions $\tau(0)=\tau(1)=0$. Replacing $\sigma=n/N$,
introducing the notation $\tau(n)=\tau(\sigma=n/N)$, and using Eqs.(1,2) we
arrive at the difference equation:

$\frac{1+\gamma v}{2}\tau(n+1)+\frac{1-\gamma
v}{2}\tau(n-1)=\tau(n)-\frac{1}{n}-\frac{1}{N-n},$ (17)

to be solved with the aforementioned boundary conditions. Due to the linearity
of this difference equation and the symmetry of the process under the change
$n\to N-n$ and $v\to-v$, the solution can be written as

$\tau(n)=T_{1}(n;\gamma v)+T_{1}(N-n;-\gamma v),$ (18)

where $T_{1}(n;\gamma)$ is the solution satisfying the boundary conditions
$T_{1}(0)=T_{1}(N)=0$ of the following difference equation:

$\displaystyle\frac{1+\gamma v}{2}T_{1}(n+1)+\frac{1-\gamma v}{2}T_{1}(n-1)$
$\displaystyle=$ $\displaystyle T_{1}(n)-\frac{1}{n}.$ (19)

As explained in the Appendix, the solution can be written in terms of the
harmonic function $H_{n}$ and the function $f(n,a)$ defined in Eqs.(A6,A7) as:

$\displaystyle
T_{1}(n;a)=\frac{1+a}{1-a}\biggl{[}\dfrac{1-a^{n}}{1-a^{N}}\left[H_{N-1}-f(N,a)\right]-\left[H_{n-1}-f(n,a)\right]\biggr{]}.$
(20)

According to the definition Eq.(9), when $\gamma v$ is replaced by $-\gamma
v$, $a$ becomes $\frac{1}{a}$ and Eq.(18) can be written as

$\tau(n;a)=T_{1}(n;a)+T_{1}(N-n;\frac{1}{a}).$ (21)

It is also possible to obtain an approximation to this expression starting
directly from the differential equation that follows from the expansion of
Eq.(16) to second order in $\Delta_{\sigma}$:

$\frac{1}{2}\frac{\Delta_{\sigma}^{2}}{\Delta_{t}}\left(R^{+}+R^{-}\right)\frac{d^{2}\tau}{d\sigma^{2}}+\frac{\Delta_{\sigma}}{\Delta_{t}}\left(R^{+}-R^{-}\right)\frac{d\tau}{d\sigma}=-1,$
(22)

or, replacing the expressions for $R^{+},R^{-},\Delta_{\sigma},\Delta_{t}$,

$\frac{d^{2}\tau}{d\sigma^{2}}+\beta\frac{d\tau}{d\sigma}=-2N\left(\frac{1}{\sigma}+\frac{1}{1-\sigma}\right).$
(23)

The solution of this equation with boundary conditions $\tau(0)=\tau(1)=0$ is

$\tau(\sigma;\beta)=2N\left[T(\sigma;\beta)+T(1-\sigma;-\beta)\right],$ (24)

with $T(\sigma;\beta)$ the solution of

$\displaystyle\frac{d^{2}T}{d\sigma^{2}}+\beta\frac{dT}{d\sigma}=-\frac{1}{\sigma},$
$\displaystyle T(0)=T(1)=0.$ (25)

given explicitly by (see footnote 1)

$\displaystyle
T(\sigma;\beta)=\frac{1}{\beta}\biggl{[}\frac{e^{-\beta\sigma}-1}{e^{\beta}-1}\mathrm{Ei}(\beta)+\frac{1-e^{\beta(1-\sigma)}}{e^{\beta}-1}(\ln(|\beta|)+\gamma_{\text{e}})+e^{-\beta\sigma}\mathrm{Ei}(\beta\sigma)-\ln(\sigma)\biggr{]},$
(26)

where $\gamma_{e}\approx 0.577$ is the Euler-Mascheroni constant and
$\mathrm{Ei}(x)=-\int_{-x}^{\infty}\exp(-z)/zdz$ the exponential integral.
Again, the continuous approximation Eq.(24) and the discrete counterpart
Eq.(21) are almost indistinguishable for all system sizes and parameter values
used in the figures. An example is given in Fig. 3, where we compare the
results of this theoretical analysis with those of computer simulations on a
complete-graph and a random initial condition $\sigma=0.5$. In the same plot
we see results corresponding to the random (Erdös-Rény) network for different
values of the average connectivity $\mu$ that we will address in a more
detailed theory developed in Sec. 4.

Figure 3: Average time $\tau$ to reach an absorbing state, in units of MCS and
rescaled by the system size $N$, for an initial condition $\sigma=0.5$ as a
function of $\beta=\gamma vN$. At the scale of the figure the analytical
predictions of the continuous approach Eq.(24), continuous (blue) line, and
the discrete time approach Eq.(21), long-dashed (light blue) line, are
indistinguishable. The open square ($\square$) symbols denote the results of
computer simulations in a complete graph (CG) for $v=0.01$, $\gamma=0.1$ and
$N=1000$. We also plot the results for the random Erdős-Rény networks with
average degree $\mu$. The theoretical prediction lines, dashed (green) for
$\mu=5$ and dotted (red) for $\mu=10$ are given by Eq.(56), and the different
symbols are the results of numerical simulations.

Using the known asymptotic expansions of the exponential integral
$\lim_{x\to\infty}\mathrm{Ei}(x)=e^{x}(x^{-1}+O(x^{-2}))$ we find that:

$\lim_{\beta\rightarrow\infty}\tau(\sigma=1/2;\beta)\rightarrow
2N\frac{\ln(\beta)+\gamma_{\text{e}}}{\beta}\to\frac{\ln(N)}{\gamma v},$ (27)

which means that the average time, in units of MCS, to reach consensus
starting from $\sigma=1/2$ in the presence of a group of biased agents,
$\gamma v\neq 0$, scales with the number of agents as $\ln(N)$. A result that
is confirmed by the numerical simulations, see Fig. 4. This is to be compared
with the limit of no bias $\gamma v\to 0$ which can be obtained directly from
Eq.(23) setting $\beta=0$, or from Eq.(24) using the expansion $\lim_{x\to
0}\mathrm{Ei}(x)=\ln(|x|)+\gamma_{\text{e}}+x+O(x^{2})$

$\lim_{\beta\rightarrow 0}\tau(\sigma=1/2;\beta)\rightarrow N\ln(4),$ (28)

a much slower and well known approach to consensus than in the biased case.

Figure 4: We plot the average time $\tau$ to reach an absorbing state, in
units of MCS, in order to check the logarithmic dependence with system size
$N$ as it follows from the asymptotic expression Eq. (27), written as
$\frac{\tau}{2N}\to\frac{\log\beta}{\beta}$, solid line. The symbols
correspond to numerical simulations on a complete graph for different values
of $\gamma v$, as indicated in the legend.

###  3.4 Time to reach preferred state, $\tau_{1}$

The average time $\tau_{1}(\sigma)$ to reach the preferred absorbing state
(for $\gamma v>0$ the preferred state is $+1$) starting from an initial
fraction $\sigma$ satisfies a recurrence relation:

$\displaystyle P_{1}(\sigma)\tau_{1}(\sigma)$ $\displaystyle=$ $\displaystyle
R^{+}(\sigma)\left[P_{1}(\sigma+\Delta_{\sigma})\tau_{1}(\sigma+\Delta_{\sigma})+P_{1}(\sigma)\Delta_{t}\right]$
(29) $\displaystyle+$ $\displaystyle
R^{-}(\sigma)\left[P_{1}(\sigma-\Delta_{\sigma})\tau_{1}(\sigma-\Delta_{\sigma})+P_{1}(\sigma)\Delta_{t}\right]$
$\displaystyle+$
$\displaystyle\left[1-R^{+}(\sigma)-R^{-}(\sigma)\right]\left[P_{1}(\sigma)\tau_{1}(\sigma)+P_{1}(\sigma)\Delta_{t}\right].$

We do not solve this recurrence relation, but proceed directly to the
continuous limit approach, in view of its accuracy. Expanding $P_{1}\tau_{1}$
to second order in $\Delta_{\sigma}=\Delta_{t}=1/N$, Eq.(29) becomes

$\frac{1}{2N}\left(R^{+}+R^{-}\right)\frac{d^{2}(P_{1}\tau_{1})}{d\sigma^{2}}+\left(R^{+}-R^{-}\right)\frac{d(P_{1}\tau_{1})}{d\sigma}=-P_{1}.$
(30)

Replacing $P_{1}(\sigma)$ from Eq.(14) in the right-hand-side, we obtain

$\frac{d^{2}(P_{1}\tau_{1})}{d\sigma^{2}}+\beta\frac{d(P_{1}\tau_{1})}{d\sigma}=-\frac{2N}{1-e^{-\beta}}\left(\frac{1}{\sigma}+\frac{1}{1-\sigma}-\frac{e^{-\beta\sigma}}{\sigma}-\frac{e^{-\beta\sigma}}{1-\sigma}\right)$
(31)

with boundary conditions $\tau_{1}(1)P_{1}(1)=\tau_{1}(0)P_{1}(0)=0$. To solve
this equation we note that the solution of

$\displaystyle\frac{d^{2}\hat{T}}{d\sigma^{2}}+\beta\frac{d\hat{T}}{d\sigma}=-\frac{e^{-\beta\sigma}}{\sigma},$
$\displaystyle\hat{T}(0)=\hat{T}(1)=0,$ (32)

is $\hat{T}(\sigma;\beta)=e^{-\beta\sigma}T(\sigma,-\beta)$, with
$T(\sigma;\beta)$ as given by Eq.(26). The solution of (31) is hence

$\displaystyle\tau_{1}(\sigma;\beta)=\frac{2N}{1-e^{-\beta\sigma}}\bigl{[}T(\sigma;\beta)+T(1-\sigma,-\beta)$
$\displaystyle-e^{-\beta\sigma}T(\sigma,-\beta)-e^{-\beta\sigma}T(1-\sigma;\beta)\bigr{]}.$
(33)

In Fig. 5 we plot the times $\tau(\sigma;\beta)$ and $\tau_{1}(\sigma;\beta)$
as given by Eqs.(24, 33), as a function of the initial value $\sigma$ for two
different values of $\beta$. For the particular case of random initial
conditions $\sigma=1/2$, it turns out that
$\tau(\sigma=1/2;\beta)=\tau_{1}(\sigma=1/2;\beta)$. When the initial fraction
$\sigma$ of $+1$ agents is smaller than $1/2$, then even in the presence of
bias it is $\tau_{1}>\tau$. In the opposite case, when we run the dynamics
starting with more than half of the agents in the preferred state and small
values of $\beta$ (see Fig. 5(a)), we observe that $\tau_{1}<\tau$. A similar
relation between $\tau$ and $\tau_{1}$ holds for the standard voter model
[30]. As showed in Fig. 5(b), for a large bias parameter $\beta=2\gamma vN$
both times converge to the same value, a result which is a consequence of the
very small probability to reach the non-preferred state.

(a) $\beta=1$

(b) $\beta=10$

Figure 5: Plot of the average time to reach consensus in any state starting
from an initial fraction $\sigma$ of agents in the preferred state,
$\tau(\sigma)$, as given by Eq.(24) continuous (black) line, and the average
time to reach consensus in the preferred state, $\tau_{1}(\sigma)$, given by
Eq.(33) dotted (red) line for $\beta=1$ (left panel), and $\beta=10$ (right
panel). Note that both times, and hence also the time $\tau_{-1}$ to reach
consensus in the non-preferred state, see Eq.(34), coincide for $\sigma=1/2$.
Parameters values: $\gamma=0.1,v=0.01,N=1000$.

The average time to reach the non-preferred state, $\tau_{-1}(\sigma;\beta)$
can be found by noting the relation :

$\tau(\sigma;\beta)=P_{1}(\sigma;\beta)\tau_{1}(\sigma;\beta)+(1-P_{1}(\sigma;\beta))\tau_{-1}(\sigma;\beta)$
(34)

and using Eqs.(14, 24, 33). We show in Figs 6(a), 6(b) that the analytical
expression for the times to reach the absorbing states $\tau_{1}$ or
$\tau_{-1}$ agree well with computer simulations of the system dynamics both
in the complete graph and in the random networks to be discussed in section 4.
The presence of a bias decreases both times $\tau_{1}$ and $\tau_{-1}$.

(a) $\tau_{1}$

(b) $\tau_{-1}$

Figure 6: Average times $\tau_{\pm 1}(\sigma)$ to reach consensus in $+1$ or
$-1$ states, in units of MCS and rescaled by the system size $N$, for an
initial condition $\sigma=0.5$ as a function of $\beta=\gamma vN$. Numerical
results (symbols) are compared with analytical predictions (solid lines) both
for the complete-graph (CG) and for the random Erdős-Rény networks with
different average degree $\mu$, as given by Eq.(33) and Eq.(4.5) for
$\tau_{1}$ and Eq.(34) for $\tau_{-1}$. Same line, symbol meanings and
parameter values than in Fig. 3.

##  4 Random networks

So far, we have only discussed the situation of a complete graph. In this
section we want to consider more general networks of interactions. A network
is fully defined through its adjacency matrix ${\cal A}$ whose elements are
${\cal A}_{i,j}=1$ if nodes $i$ and $j$ are connected and ${\cal A}_{i,j}=0$
otherwise. This detailed information is most times simplified to the knowledge
of the degree distribution $P_{k}=N_{k}/N$, being $N_{k}$ the number of nodes
with degree $k$ and $N$ the total number of nodes. The average degree is
$\mu=\sum_{k}kP_{k}$, and the second moment $\mu_{2}=\sum_{k}k^{2}P_{k}$. A
random or connected Erdős-Rény network is constructed by linking each possible
pair of nodes with a given probability $p$. In the large $N$ limit such a
network follows a Poisson distribution for $P_{k}$, $k>0$, with an average
value $\mu=pN$ and second moment $\mu_{2}=\mu^{2}+\mu$. As usual, in the
numerical simulations we disregard those networks that can be split in two
disconnected parts.

###  4.1 Pair approximation

At the level of the mean-field approximation, the only relevant variable is
the fraction $\sigma(t)$ of sites in the state $s=+1$ as a function of time.
Within this mean-field approximation, the probability that a randomly selected
pair of neighbors is “active”, i.e. both sites are in different states, is
$2\sigma(1-\sigma)$, a result coming from the statistical independence assumed
in the approximation. At a more detailed level, the pair approximation
considers correlations between the states of different connected sites by
introducing the density of active links, $\rho(t)$, as a new dynamical
variable [31, 54, 10, 55, 56]. This approach is reasonably successful to treat
random network configurations without degree-degree correlations such as an
Erdős-Rény network. If $i,j$ are connected nodes, we define
$\rho^{i,j}=\text{Prob}(s_{i}=-s_{j})$ as the probability that the link $i,j$
is active. The global density of active links is then
$\rho=\frac{1}{L}\sum_{\langle i,j\rangle}\rho^{i,j}$, being $L=\mu N/2$ the
total number of links. The pair approximation assumes that the probability of
a link being active is independent of the state of the other links, hence
$\rho^{i,j}\approx\rho$. Consistent with this approximation it is further
assumed that $\text{Prob}(s_{i}=-s,s_{j}=s)\approx\rho/2$ independently of the
value of $s=\pm 1$.

Beyond the mean-field approach developed in the previous sections, the pair
approximation uses $\rho(t)$ and $\sigma(t)$ as an independent pair of
dynamical variables to describe the state of the system. Note, however, that
$\rho(t)=0$ is only consistent with $\sigma(t)=0,1$, coming from the fact that
a consensus state, one in which all nodes hold the same value of their state
variable, has no active links. It is of course possible to include further
variables in the analysis. For instance, the set of degree-dependent fractions
defined as the ratio $\sigma_{k}=n_{k}/N_{k}$ between the number $n_{k}$ of
nodes which are in state $+1$ and have degree $k$ and the total number $N_{k}$
of nodes with degree $k$. It is obviously $\sigma=\sum_{k}P_{k}\sigma_{k}$. A
better description of the state of the network replaces the fraction of nodes
in the state $+1$ by the degree-weighted fraction
$\sigma_{L}=\frac{1}{\mu}\sum_{k}P_{k}k\sigma_{k}$. For a regular or all-to-
all connected network where $P_{k}=\delta_{k,\mu}$, the degree-weighted
fraction $\sigma_{L}(t)$ coincides with $\sigma(t)$.

A complete and closed description of the evolution of the dynamical variables
$\rho(t),\,\sigma(t),\,\sigma_{L}(t)$ is possible within the context of the
pair approximation. This description is, however, rather cumbersome [57] and
we have decided to present here a simplified treatment based on [31]. The idea
is to consider a random walk not for the variable $\sigma(t)$ but for the
variable $\sigma_{L}(t)$. A microscopic update $s_{i}\leftarrow s_{j}$ where
node $i$ with degree $k_{i}=k$ copies the state of node $j$ modifies
$\sigma_{L}$ by an amount $\pm\dfrac{k}{\mu N}\equiv\pm\Delta_{k}$. We now
compute the probability $R_{k}^{+}$ that, given that node $i$ with degree
$k_{i}=k$ has been chosen for updating, the change of the degree-weighted
fraction $\sigma_{L}$ is $+\Delta_{k}$:

$\displaystyle R_{k}^{+}$ $\displaystyle=$
$\displaystyle\text{Prob}(s_{i}=-1,s_{j}=1|k_{i}=k)\left[\text{Prob}(i\text{
no bias})\frac{1}{2}+\text{Prob}(i\text{ bias})\frac{1+v}{2}\right]$ (35)
$\displaystyle=$ $\displaystyle\frac{\rho}{2}\frac{1+\gamma v}{2},$

where we have used the approximation
$\text{Prob}(s_{i}=-1,s_{j}=+1|k_{i}=k)=\text{Prob}(s_{i}=-1,s_{j}=+1)=\dfrac{\rho}{2}$
and that the probability that a node is biased or unbiased is independent of
its degree. Similarly for the probability $R_{k}^{-}$ that, given that node
$i$ with degree $k_{i}=k$ has been chosen for updating, the change of the
degree-weighted fraction $\sigma_{L}$ is $-\Delta_{k}$ we obtain:

$\displaystyle R_{k}^{-}$ $\displaystyle=$
$\displaystyle\frac{\rho}{2}\frac{1-\gamma v}{2}.$ (36)

Note that we still can define the rates $R^{+}(\sigma),\,R^{-}(\sigma)$ that
the fraction $\sigma$ decreases or increases, respectively, by an amount
$\Delta_{\sigma}=1/N$. Within this context, they are equal to $R_{k}^{+}$ and
$R_{k}^{-}$, respectively, as these are independent of $k$ due to the
approximations considered.

To proceed, we need an equation for the time evolution of $\rho$. We follow
closely the derivation of [31] and note that every time a node with degree $k$
and $\ell$ active links is updated, the density of active links varies in an
amount $\Delta\rho=\dfrac{2(k-2\ell)}{\mu N}$. As time increases by
$\Delta_{t}=1/N$ after every node update, we write:

$\displaystyle\frac{d\rho}{dt}$ $\displaystyle=$
$\displaystyle\sum_{k}P_{k}\left.\frac{\Delta\rho}{\Delta_{t}}\right|_{k}$
$\displaystyle=$ $\displaystyle\sum_{k}\frac{P_{k}}{1/N}\sum_{s=\pm
1}P(s\to-s)\sigma_{s}\sum_{\ell=0}^{k}B(\ell,k|s)\frac{\ell}{k}\frac{2(k-2\ell)}{\mu
N},$

where $\left.\dfrac{\Delta\rho}{\Delta_{t}}\right|_{k}$ denotes the average
change in $\rho$ when a node of degree $k$ is chosen, $P(s\to-s)$ is the
probability that the proposed change $s\to-s$ is accepted, and $B(\ell,k|s)$
is the conditional probability that $\ell$ of the $k$ links connected to a
node are active, given that the node is in the state $s$. We have introduced
the notation $\sigma_{1}=\sigma$, $\sigma_{-1}=1-\sigma$. This expression is
equivalent to

$\frac{d\rho}{dt}=\frac{2}{\mu}\sum_{k}P_{k}\sum_{s=\pm
1}P(s\rightarrow-s)\sigma_{s}\left(\langle\ell\rangle_{k,s}-\frac{2}{k}\langle\ell^{2}\rangle_{k,s}\right)$
(38)

where $\langle\ell\rangle_{k,s}$, is the average number of active neighbors of
a node in state $s$ and degree $k$. Using the pair approximation and
neglecting correlation of second and higher neighbors, it turns out that
$B(\ell,k|s)$ becomes a binomial distribution, whose first and second moments
are

$\displaystyle\langle\ell\rangle_{k,s}=\frac{\rho
k}{2\sigma_{s}},\,\langle\ell^{2}\rangle_{k,s}=\frac{\rho
k}{2\sigma_{s}}+\frac{k(k-1)\rho^{2}}{4\sigma_{s}^{2}}.$ (39)

Replacing in Eq.(38), using $P(-1\rightarrow 1)=\frac{1+\gamma v}{2}$,
$P(1\rightarrow-1)=\frac{1-\gamma v}{2}$, we arrive at

$\frac{d\rho}{dt}=\frac{2\rho}{\mu}\left[\left(\frac{\mu}{2}-1\right)-\rho\left(\mu-1\right)\frac{1+\gamma
v(2\sigma-1)}{4\sigma\left(1-\sigma\right)}\right].$ (40)

This equation has to be combined with the evolution equation for the fraction
$\sigma$:

$\frac{d\sigma}{dt}=R^{+}(\sigma)-R^{-}(\sigma)=\frac{\gamma v}{2}\rho.$ (41)

The set of coupled equations (40) and (41) are the basis of our subsequent
analysis. They are the result of the pair approximation which neglects finite
size fluctuations and it is therefore valid in the thermodynamic limit. In
this limit and in the absence of bias ($v=0$), $\sigma$ is a conserved
quantity and there is a stationary solution with a finite value of $\rho$.
However, when bias is present, the stationary solution fulfills $\rho=0$,
indicating that the absorbing state is reached by the intrinsic dynamics of
the system in the absence of finite size fluctuations. The dynamical equations
(40) and (41) reproduce the ones obtained in the preferred language study
[10], when setting $\gamma=1$, i.e. when all agents are biased, although our
analysis is different. Instead of finding the general solution
$\sigma(t),\,\rho(t)$ with given boundary conditions, we note that for $\gamma
v$ small the time scale of Eq.(41) indicates that $\sigma(t)$ is a slow
variable, and we assume that the dependence of $\rho(t)$ in time is through
the relation $\rho(t)=\rho(\sigma(t))$. Dividing Eq.(40) by Eq.(41) we get a
closed differential equation to find the dependence $\rho(\sigma)$:

$\frac{d\rho}{d\sigma}=\frac{4}{\mu\gamma
v}\left[\left(\frac{\mu}{2}-1\right)-\rho\left(\mu-1\right)\frac{1+\gamma
v(2\sigma-1)}{4\sigma\left(1-\sigma\right)}\right].$ (42)

The solution satisfying the boundary conditions
$\rho(\sigma=1)=\rho(\sigma=0)=0$ is

$\displaystyle\rho(\sigma)$ $\displaystyle=$
$\displaystyle\frac{\mu-2}{\mu-1+\gamma
v}2\sigma(1-\sigma)^{a_{1}}{{}_{2}F_{1}}\left(a_{1},a_{2};1+a_{2};\sigma\right),$
(43) $\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle\frac{1+\gamma
v}{\gamma v}\frac{\mu-1}{\mu},$ $\displaystyle a_{2}$ $\displaystyle=$
$\displaystyle\frac{\gamma v+\mu-1}{\gamma v\mu},$

where ${{}_{2}F_{1}}(\cdot)$ is the hypergeometric function. It is possible to
check the limits

$\displaystyle\lim_{v\to 0}\rho(\sigma)$ $\displaystyle=$
$\displaystyle\frac{\mu-2}{\mu-1}2\sigma(1-\sigma),\quad\forall\mu,$ (44)
$\displaystyle\lim_{\mu\to\infty}\rho(\sigma)$ $\displaystyle=$ $\displaystyle
2\sigma(1-\sigma),\quad\forall\gamma v.$ (45)

While the first limit coincides with the one obtained in [31], the last limit
is an important check of the consistency of the calculation. When $\mu=N-1$
every two nodes are connected and it follows the exact relation
$\rho=\dfrac{2\sigma(1-\sigma)}{1-1/N}$ and, for $N\to\infty$,
$\rho=2\sigma(1-\sigma)$ is mandatory independently on the value of the bias
parameter $v$ or the fraction of biased agents $\gamma$.

We simplify the complicated functional relation $\rho(\sigma)$ as given in
Eq.(43) in order to use it in further calculations and get full analytical
expressions for $P_{1}$ and $\tau$. To this end we use the previous asymptotic
limits and expand around $\sigma=1/2$:

$\frac{4\sigma(1-\sigma)\rho(1/2)}{\rho(\sigma)}=1-c_{1}(\mu,v)(2\sigma-1)+O((2\sigma-1)^{2}).$
(46)

While $\rho(1/2)$ and $c_{1}(\mu,v)$ can be fully expressed in terms of the
hypergeometric function, it is possible to use approximate expressions valid
for small $v$, namely $\rho(1/2)\approx\dfrac{\mu-2}{2(\mu-1)}$ and
$c_{1}\approx\dfrac{\gamma v}{\mu-1}$. This leads to the approximation

$\rho(\sigma)\approx\frac{2(\mu-2)\sigma(1-\sigma)}{\mu-1+\gamma
v(1-2\sigma)}.$ (47)

The essence of our approximation, in contrast to other approaches in the
literature [31, 10], is that $\rho(t)$ follows adiabatically $\sigma(t)$ and
we can use Eq.(47) using the time-dependent values $\sigma(t)$ and $\rho(t)$.
The comparison with computer simulations shown in Fig. 7(a) proves the
goodness of this approximation for Erdős-Rény networks and two different
values of the average degree $\mu$. In Fig. 7(b) we plot the time dependence
of the interface density in a single realization of the dynamics of the system
and compare it against the value of the plateau
$\rho(1/2)\approx\dfrac{\mu-2}{2(\mu-1)}$ that follows from the stationary
solution of Eq. (40) [31]. We now introduce this approximation to analyze the
behavior of the fixation probability and the times to reach the different
consensus states.

(a) $\rho(\sigma)$

(b) $\rho(t)$

Figure 7: Density of active links $\rho$ for a single realization of the
system dynamics of $N=10000$ agents on an Erdős-Rény network and with initial
conditions $\sigma(0)=1/2$. In panel $(a)$ we display with symbols
$\rho(\sigma)$ as obtained from a parametric plot of the numerical simulation
values $(\sigma(t),\rho(t))$ for average connectivity $\mu=5$ (circles) and
$\mu=10$ (triangles) . The solid lines correspond to the adiabatic
approximation Eq.(43). Panel $(b)$ shows with symbols the time evolution of
$\rho(t)$ coming from the numerical simulations. The corresponding solid lines
are the result of the numerical integration of Eqs.(40,41). The horizontal
solid lines correspond to the plateau values
$\rho\approx\dfrac{\mu-2}{2(\mu-1)}$ as given by Eq.(47) for $\sigma=1/2$.

###  4.2 Absorbing state

In this subsection we go beyond the pair approximation, introducing finite
size effects which are neglected in the previous treatment. In order to
identify the absorbing state of the system, we resort to the relevant master
equation for the probability $P(m,t)$ that the system has magnetization $m$ at
time $t$. In a time step, a node with opinion $s$ flips with probability
$\sigma_{s}P(-s|s)P(s\rightarrow-s)$ after which the magnetization
$m=2\sigma-1$ changes by $\Delta_{m}=2s\Delta_{k}$ where
$\Delta_{k}=\frac{k}{\mu N}$ with $k$ the number of neighbors of the selected
node.

Following a similar approach as before, we can arrive to the following Fokker-
Planck equation

$\displaystyle\frac{\partial P\left(m,t\right)}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{\mu-2}{\mu-1}\frac{\gamma
v}{2}\frac{\partial\left[\frac{1-m^{2}}{1-\frac{\gamma
v}{\mu-1}m}P\left(m,t\right)\right]}{\partial m}$ $\displaystyle+$
$\displaystyle\frac{\mu-2}{\mu-1}\frac{1}{2N_{\mu}}\frac{\partial^{2}\left[\frac{1-m^{2}}{1-\frac{\gamma
v}{\mu-1}m}P\left(m,t\right)\right]}{\partial m^{2}}$

where $N_{\mu}=\dfrac{\mu^{2}}{\mu_{2}}N$ or $N_{\mu}=\dfrac{\mu}{\mu+1}N$ for
the Erdős-Rény network. Here we see that if $\gamma$ or $v$ are set to 0 we
obtain the same result as in [31]. We can now perform a similar analysis to
that of Eq.(5) to conclude that for $\gamma v>0$ it is
$P_{\text{st}}=\delta(m-1)$, and that the approach to the stationary state
$m=1$ occurs in a characteristic time scale $\frac{\mu-1}{(\mu-2)\gamma v}$.

###  4.3 Fixation probability $P_{1}$

To compute the fixation probability $P_{1}$ we consider the random walk in the
degree-weighted fraction $\sigma_{L}$ which takes the value $\sigma_{L}=1$ in
the consensus state $s_{i}=s,\,\forall i$. It can be computed by a reasoning
similar to the one that led to Eq.(7), but considering the different
contributions according to the degree of the selected node for updating:

$\displaystyle
P_{1}(\sigma_{L})=\sum_{k}P_{k}\biggl{[}R_{k}^{+}P_{1}(\sigma_{L}+\Delta_{k})+R_{k}^{-}P_{1}(\sigma_{L}-\Delta_{k})+[1-R_{k}^{+}-R_{k}^{-}]P_{1}(\sigma_{L})\biggr{]}.$
(49)

Expanding up to second order in $\Delta_{k}=\dfrac{k}{\mu N}$ and replacing
the rates (35,36) we arrive at:

$\frac{d^{2}P_{1}}{d\sigma_{L}^{2}}+\beta_{\mu}\frac{dP_{1}}{d\sigma_{L}}=0,$
(50)

here $\beta_{\mu}=2\gamma vN_{\mu}$. Hence the solution Eq.(14) is still valid
if we replace $\beta$ by $\beta_{\mu}$ and $\sigma$ by $\sigma_{L}$.
Nevertheless, as proven by a more detailed analysis [55], the variable
$\sigma(t)$ follows $\sigma_{L}(t)$ and we can replace one variable by the
other. However, it is essential to do the random-walk analysis in terms of the
variable $\sigma_{L}$, otherwise the dependence on $N_{\mu}$ is lost. As shown
in Fig. 1, where we plot the fixation probability for the complete-graph and
two Erdős-Rény networks with main degree $\mu=5$ and $\mu=10$, the agreement
of Eq.(14) with the numerical results is very good if we include the system-
size dependence in $N_{\mu}=\dfrac{\mu}{\mu+1}N$.

###  4.4 Time to reach consensus, $\tau$

We modify the approach used for the complete-graph by noticing that in the
case of a heterogeneous network, the change in the fraction of nodes in state
$+1$ depends now on the degree $k$ of the node selected for update as
$\Delta_{k}$. Hence we add all contributions weighted by its probability and
modify Eq.(16) as

$\displaystyle\tau(\sigma_{L})=$
$\displaystyle\sum_{k}P_{k}\biggl{[}R_{k}^{+}\left[\tau(\sigma_{L}+\Delta_{k})+\Delta_{t}\right]+R_{k}^{-}\left[\tau(\sigma_{L}-\Delta_{k})+\Delta_{t}\right]$
(51)
$\displaystyle+\left[1-R_{k}^{+}-R_{k}^{-}\right]\left[\tau(\sigma_{L})+\Delta_{t}\right]\biggr{]}.$

Expanding to second order in $\Delta_{k}=\frac{k}{\mu N}$ and replacing
$\Delta_{t}=\frac{1}{N}$, this equation becomes

$\sum_{k}P_{k}\left[\frac{k^{2}}{2\mu^{2}N}\left(R_{k}^{+}+R_{k}^{-}\right)\frac{d^{2}\tau}{d\sigma^{2}}+\frac{k}{\mu}\left(R_{k}^{+}-R_{k}^{-}\right)\frac{d\tau}{d\sigma}\right]=-1,$
(52)

and after replacing the rates as given by Eqs. (35,36,47)

$\frac{d^{2}\tau}{d\sigma^{2}}+\beta_{\mu}\frac{d\tau}{d\sigma}=-2N_{\mu}\left(\frac{A}{\sigma}+\frac{B}{1-\sigma}\right),$
(53)

where

$\displaystyle A\equiv\frac{\mu-1+\gamma v}{\mu-2},$ (54) $\displaystyle
B\equiv\frac{\mu-1-\gamma v}{\mu-2},$ (55)

and, in view of the aforementioned equivalence, we have replaced $\sigma_{L}$
by $\sigma$. The solution of Eq.(53) with the boundary conditions
$\tau(0)=\tau(1)=0$ can be written in terms of the function $T(\sigma;\beta)$
defined in Eq.(26) as

$\tau(\sigma;\beta_{\mu})=2N_{\mu}\left[AT(\sigma;\beta_{\mu})+BT(1-\sigma;-\beta_{\mu})\right].$
(56)

In Fig. 3 we compare this analytical solution with the results of computer
simulations for Erdős-Rény networks with average degree $\mu=5$, $\mu=10$. We
observe that the time to reach consensus increases for decreasing average
degree $\mu$ and that it is larger in an Erdős-Rény network than in the all-
to-all configuration. Note, again, that analyzing the random walk in terms of
$\sigma$ instead of $\sigma_{L}$ we would have missed the dependence on
$N_{\mu}$ which provides a much better fit to the numerical data.

###  4.5 Time to reach preferred state, $\tau_{1}$

We start by writing an analog expression of Eq.(29) for heterogeneous networks
(see Sec. 4.4) but considering a random walk in the $\sigma_{L}$ variable with
contributions depending on the degree $k$

$\displaystyle P_{1}(\sigma_{L})\tau_{1}(\sigma_{L})$ $\displaystyle=$
$\displaystyle\sum_{k}P_{k}\left\\{R_{k}^{+}\left[P_{1}(\sigma_{L}+\Delta_{k})\tau_{1}(\sigma_{L}+\Delta_{k})+P_{1}(\sigma_{L})\Delta_{t}\right]\right.$
(57)
$\displaystyle+R_{k}^{-}\left[P_{1}(\sigma_{L}-\Delta_{k})\tau_{1}(\sigma_{L}-\Delta_{k})+P_{1}(\sigma_{L})\Delta_{t}\right]$
$\displaystyle+\left.\left[1-R_{k}^{+}(\sigma_{L})-R_{k}^{-}(\sigma_{L})\right]\left[P_{1}(\sigma_{L})\tau_{1}(\sigma_{L})+P_{1}(\sigma_{L})\Delta_{t}\right]\right\\}.$

Expanding to second order in $\Delta_{k}=\frac{k}{\mu N}$ and replacing
$\Delta_{t}=\frac{1}{N}$, we obtain

$\displaystyle\sum_{k}P_{k}\biggl{[}\frac{k^{2}}{2\mu^{2}N}\left(R_{k}^{+}+R_{k}^{-}\right)\frac{d^{2}(P_{1}\tau_{1})}{d\sigma_{L}^{2}}+\frac{k}{\mu}\left(R_{k}^{+}-R_{k}^{-}\right)\frac{d(P_{1}\tau_{1})}{d\sigma_{L}}\biggr{]}=-P_{1}.$
(58)

Using in the right-hand-side the expression of $P_{1}(\sigma)$ from Eq.(14)
with $\beta$ replaced by $\beta_{\mu}$, and substituting again $\sigma_{L}$ by
$\sigma$, we obtain

$\displaystyle\frac{d^{2}(P_{1}\tau_{1})}{d\sigma^{2}}+\beta_{\mu}\frac{d(P_{1}\tau_{1})}{d\sigma}=-\frac{2N_{\mu}}{1-e^{-\beta_{\mu}}}\left(\frac{A}{\sigma}+\frac{B}{1-\sigma}-\frac{Ae^{-\beta_{\mu}\,\sigma}}{\sigma}-\frac{Be^{-\beta_{\mu}\,\sigma}}{1-\sigma}\right),$
(59)

where $A$, $B$ are given by Eqs. (54,55). The solution with boundary
conditions $\tau_{1}(1)P_{1}(1)=\tau_{1}(0)P_{1}(0)=0$ can be expressed in
terms of the function defined in Eq. (26) as

$\displaystyle\tau_{1}(\sigma;\beta_{\mu})=$
$\displaystyle\frac{2N_{\mu}}{1-e^{-\beta_{\mu}\,\sigma}}\biggl{[}AT(\sigma;\beta_{\mu})+BT(1-\sigma,-\beta_{\mu})-Ae^{-\beta_{\mu}\,\sigma}T(\sigma;-\beta_{\mu})-Be^{-\beta_{\mu}\,\sigma}T(1-\sigma;\beta_{\mu})\biggr{]}$

In this case, it is no longer true that
$\tau_{1}(\sigma=1/2;\beta_{\mu})=\tau(\sigma=1/2;\beta_{\mu})$, as it was in
the mean-field approximation, although both times scale in the same way
$\tau_{1}(\sigma=1/2;\beta_{\mu}),\tau(\sigma=1/2;\beta_{\mu})\sim\frac{B}{\gamma
v}\ln(N_{\mu})$ as $N\to\infty$. In Fig. 6 we compare these analytical results
for $\tau_{1}(\sigma)$ and $\tau_{-1}(\sigma)$, derived from Eq. (34) with the
results of numerical simulations in Erdős-Rény networks.

##  5 Biased-dependent topology

So far we have assumed that links amongst agents, and hence the possibility of
interaction, form, not just randomly, but also independently of whether they
are biased or unbiased. In this section we will assume that the bias
influences not only the interactions between individuals but, more
importantly, the way they are connected, their network topology. Our goal is
to determine whether the biased community is able to influence the whole
system more efficiently by establishing their links in a more organized
fashion. By more efficiently we mean that consensus to the preferred state
occurs with a higher probability and in a shorter average time. To this end we
will consider different network structures in which the connections between
nodes will depend on their biased/unbiased label.

For the sake of simplicity, we assume that the total number of links in the
system, $L$, is fixed and given. Consequently, the total average degree of
nodes $\mu=2L/N$ is also fixed. We denote the total number of links between
biased-biased, unbiased-unbiased and biased-unbiased pairs of nodes as
$L_{\text{\tiny{BB}}}$, $L_{\text{\tiny{UU}}}$, and $L_{\text{\tiny{UB}}}$
respectively. Thus the total number of links is
$L=L_{\text{\tiny{BB}}}+L_{\text{\tiny{UU}}}+L_{\text{\tiny{UB}}}=\frac{1}{2}\mu
N$. Let us denote by $\mu_{\text{\tiny{B}}}$, $\mu_{\text{\tiny{U}}}$ the
average degrees of biased and unbiased nodes, respectively. They are related
to the global average degree by
$\mu=\gamma\mu_{\text{\tiny{B}}}+(1-\gamma)\mu_{\text{\tiny{U}}}$, if one
assumes that the degree distribution $P_{k,U}$ for unbiased nodes is
independent of the degree distribution $P_{k,B}$ for biased nodes. We can
write the average degree of the given node as a sum of the connections to
biased and unbiased neighbors, i.e.,
$\mu_{\text{\tiny{B}}}=\mu_{\text{\tiny{BB}}}+\mu_{\text{\tiny{BU}}}$ and
$\mu_{\text{\tiny{U}}}=\mu_{\text{\tiny{UB}}}+\mu_{\text{\tiny{UU}}}$, where
$\mu_{{}_{\text{XY}}}$ is the average number of links from an X-type node to a
Y-type neighbor. Note that in general $\mu_{\text{\tiny{UB}}}$ and
$\mu_{\text{\tiny{BU}}}$ are not equal but they are related as
$\mu_{\text{\tiny{BU}}}N_{\text{\tiny{B}}}=\mu_{\text{\tiny{UB}}}N_{\text{\tiny{U}}}$
or $\gamma\mu_{\text{\tiny{BU}}}=(1-\gamma)\mu_{\text{\tiny{UB}}}$. For the
generation of the biased-dependent topology networks we will use as a control
parameter the ratio $\delta$ of links a biased node has to biased neighbors
with respect to the number of links an unbiased node shares with its unbiased
neighbors, i.e.
$\delta=\dfrac{\mu_{\text{\tiny{BB}}}}{\mu_{\text{\tiny{UU}}}}$. In the case
of a biased-independent random network this control parameter takes the value
$\delta^{\text{rand}}=\dfrac{\gamma}{1-\gamma}$.

There are many ways in which one can modify the links in order to vary the
parameter $\delta$ above or below its random value $\delta^{\text{rand}}$.
Amongst all possibilities we have chosen to compare one case (so-called model
I) in which the average connectivity of each agent remains always equal to
$\mu$ and another case (so-called model II) in which the number of UB links is
kept equal to that of the random network. In summary, in model I, we set
$\mu_{\text{\tiny{B}}}=\mu_{\text{\tiny{U}}}=\mu$, while in model II we set
$\mu_{\text{\tiny{UB}}}=\gamma\mu$. The fulfillment of these conditions, given
the constraints listed in Table 1, leads, after a simple but lengthy algebra,
to the values of the parameters listed in the same table.

quantity |  | Erdős-Rény |  | model I |  | model II
---|---|---|---|---|---|---
$\mu_{\text{\tiny{B}}}$ |  | $\mu$ |  | $\mu$ |  | $\dfrac{1+\delta-\gamma[2+\delta-\gamma(1+\delta)]}{1-\gamma(1-\delta)}\mu$
$\mu_{\text{\tiny{U}}}$ |  | $\mu$ |  | $\mu$ |  | $\dfrac{1-\gamma[1-\gamma(1+\delta)]}{1-\gamma(1-\delta)}\mu$
$\mu_{\text{\tiny{UU}}}$ |  | $(1-\gamma)\mu$ |  | $\dfrac{1-2\gamma}{1-\gamma(1+\delta)}\mu$ |  | $\dfrac{(1-\gamma)^{2}+\gamma^{2}}{1-\gamma(1-\delta)}\mu$
$\mu_{\text{\tiny{UB}}}$ |  | $\gamma\mu$ |  | $\dfrac{\gamma(1-\delta)}{1-\gamma(1+\delta)}\mu$ |  | $\gamma\mu$

$\mu_{\text{\tiny{BU}}}=\dfrac{1-\gamma}{\gamma}\mu_{\text{\tiny{UB}}}$
---
$\mu_{\text{\tiny{BB}}}=\delta\mu_{\text{\tiny{UU}}}$
$\mu=\gamma\mu_{\text{\tiny{B}}}+(1-\gamma)\mu_{\text{\tiny{B}}}$
$\mu_{\text{\tiny{B}}}=\mu_{\text{\tiny{BB}}}+\mu_{\text{\tiny{BU}}}$
$\mu_{\text{\tiny{U}}}=\mu_{\text{\tiny{UB}}}+\mu_{\text{\tiny{UU}}}$
$N_{\text{\tiny{B}}}=\gamma N$
$N_{\text{\tiny{U}}}=(1-\gamma)N$
$L_{\text{\tiny{BB}}}=\frac{1}{2}\mu_{\text{\tiny{BB}}}N_{B}$
$L_{\text{\tiny{UU}}}=\frac{1}{2}\mu_{\text{\tiny{UU}}}N_{U}$
$L_{\text{\tiny{UB}}}=\frac{1}{2}(\mu_{\text{\tiny{BU}}}N_{B}+\mu_{\text{\tiny{UB}}}N_{U})=\mu_{\text{\tiny{BU}}}N_{B}=\mu_{\text{\tiny{UB}}}N_{U}$
$L=L_{\text{\tiny{UU}}}+L_{\text{\tiny{BB}}}+L_{\text{\tiny{UB}}}=\frac{1}{2}\mu
N$
$p_{\text{\tiny{UU}}}=\mu_{\text{\tiny{UU}}}/N_{U}$
$p_{\text{\tiny{BB}}}=\mu_{\text{\tiny{BB}}}/N_{B}$
$p_{\text{\tiny{BU}}}=\mu_{\text{\tiny{BU}}}/N_{U}=\mu_{\text{\tiny{UB}}}/N_{B}$

Table 1: Comparison of relevant networks quantities, and definitions and
relations among them, for an Erdős-Rény network and the biased-dependent
topologies. Setting the values of $\mu$, $\gamma$, $\delta$ and $N$ determines
all other possible quantities. For the Erdős-Rény network it is
$\delta^{\text{rand}}=\frac{\gamma}{1-\gamma}$, while for model I and II,
$\delta$ is a free parameter.

As, obviously, the quantities
$\mu_{\text{\tiny{UU}}},\,\mu_{\text{\tiny{BB}}},\,\mu_{\text{\tiny{UB}}},\,\mu_{\text{\tiny{BU}}}$
must be all non-negative, it follows from Table 1 that when constructing model
I it should be
$\text{sign}(1-2\gamma)=\text{sign}(1-\delta)=\text{sign}(1-\gamma(1+\delta))$.
A simple manipulation of these conditions allows us to conclude that there are
two regions of allowed parameters $(\gamma$, $\delta)$ for model I, namely,
$(\gamma<1/2,\delta<1)$, and $(\gamma>1/2,\delta>1)$. Although in the case of
Model II we do not have any of those limitations, we restrict our posterior
analysis to $\gamma<1/2$, where the biased community is a minority.

When analyzing the behavior of the average number of links
$\mu_{\text{\tiny{BB}}},\mu_{\text{\tiny{UU}}},\mu_{\text{\tiny{BU}}},\mu_{\text{\tiny{UB}}}$
as given in Table. 1, for model I it turns out that if $\delta$ increases
above the value $\delta^{\text{rand}}=\frac{\gamma}{1-\gamma}$, then
$\mu_{\text{\tiny{BB}}}$ and $\mu_{\text{\tiny{UU}}}$ increase and
$\mu_{\text{\tiny{UB}}}$ and $\mu_{\text{\tiny{BU}}}$ decrease with respect to
the respective values $\mu_{\text{\tiny{BB}}}^{\text{rand}}=\gamma\mu$,
$\mu_{\text{\tiny{UU}}}^{\text{rand}}=(1-\gamma)\mu$,
$\mu_{\text{\tiny{BU}}}^{\text{rand}}=\gamma\mu$,
$\mu_{\text{\tiny{UB}}}^{\text{rand}}=(1-\gamma)\mu$ they adopt in a random
Erdős-Rény network with average degree $\mu$. The opposite behavior, namely
$\mu_{\text{\tiny{BB}}}<\mu_{\text{\tiny{BB}}}^{\text{rand}}$,
$\mu_{\text{\tiny{UU}}}<\mu_{\text{\tiny{UU}}}^{\text{rand}}$,
$\mu_{\text{\tiny{BU}}}>\mu_{\text{\tiny{BU}}}^{\text{rand}}$,
$\mu_{\text{\tiny{UB}}}>\mu_{\text{\tiny{UB}}}^{\text{rand}}$, occurs for
$\delta<\delta^{\text{rand}}$. In the case of Model II,
$\mu_{\text{\tiny{UB}}}$ and $\mu_{\text{\tiny{BU}}}$ do not vary with
$\delta$, but $\mu_{\text{\tiny{BB}}}$ increases and $\mu_{\text{\tiny{UU}}}$
decreases with respect to their random values
$\mu_{\text{\tiny{BB}}}^{\text{rand}}$, $\mu_{\text{\tiny{UU}}}^{\text{rand}}$
when $\delta>\delta^{\text{rand}}$, and the opposite behavior
$\mu_{\text{\tiny{BB}}}<\mu_{\text{\tiny{BB}}}^{\text{rand}}$,
$\mu_{\text{\tiny{UU}}}>\mu_{\text{\tiny{UU}}}^{\text{rand}}$ if
$\delta<\delta^{\text{rand}}$. This allows us to interpret Model I as follows:
In order to increase $\delta>\delta^{\text{rand}}$, start from a random Erdős-
Rény network and rewire the necessary number of UB links (with equal
probability) either as BB or as UU links. Analogously, in order to decrease
$\delta<\delta^{\text{rand}}$, rewire an equal number of BB and UU links as UB
links. Similarly, we can interpret model II as follows: In order to increase
$\delta>\delta^{\text{rand}}$, start from a random Erdős-Rény network and move
UU links to BB links. Analogously, in order to decrease
$\delta<\delta^{\text{rand}}$, move BB links to UU links. See Fig. 8 for a
schematic representation of these two strategies. In Fig. 9 we show some
characteristic networks for values of $\delta$ equal, smaller and larger than
the $\delta^{\text{rand}}$ value of the biased-independent Erdős-Rény case in
the case of model II.

In practice, those biased-dependent networks are constructed starting from $N$
disconnected nodes, a fraction $\gamma$ of which are biased, and linking each
possible pair of nodes with probabilities $p_{\text{\tiny{BB}}}$,
$p_{\text{\tiny{UU}}}$,or $p_{\text{\tiny{UB}}}$ if both nodes are biased,
both nodes are unbiased, or one node is biased and the other unbiased,
respectively. The random Erdős-Rény network uses the same probability
$p=\mu/N$ for the three cases. In order to achieve the correct network
characteristics as before one must use the values of
$p_{\text{\tiny{BB}}},\,p_{\text{\tiny{UU}}},\,p_{\text{\tiny{UB}}}$ listed in
Table 1.

(a) model I

(b) model II

Figure 8: Schematic representation of the strategy to form biased-dependent
topologies. Starting from a bias-independent Erdős-Rény network, model I moves
UB links to either UU or BB in order to increase
$\delta>\delta^{\text{rand}}$, while the opposite moves of UU or BB to UB
links decrease $\delta<\delta^{\text{rand}}$. This keeps the average
connectivity of each agent unchanged. In model II the movement of a UU link to
a BB type increases $\delta>\delta^{\text{rand}}$, while the opposite move of
BB to UU decreases $\delta<\delta^{\text{rand}}$. This leaves constant the
number of UB links joining both communities.

(a) $\gamma=0.2$ ($\delta^{\text{rand}}=0.25$)

(b) $\gamma=0.2$, $\delta=0.1<\delta^{\text{rand}}$

(c) $\gamma=0.2$, $\delta=10>\delta^{\text{rand}}$

Figure 9: Here we plot for Model II representative Erdős-Rény networks of
$N=50$ nodes with average degree $\mu=5$ with a $\gamma$ fraction of biased
nodes. $\delta^{\text{rand}}=\frac{\gamma}{1-\gamma}$. On the left (red
circles) we draw biased nodes and on the right (green circles) unbiased nodes.
Black links connect nodes of the same type (i.e. UU and BB pairs), whereas
orange links represent connections between pairs of nodes of different types
(i.e. BU).

We now discuss the type of communities that biased and unbiased agents form in
each model for different values of the parameters $(\gamma,\delta)$. To this
end we introduce, as a measure of how strongly united a community is, the
ratio of the number of links that this community holds inside to the number of
links it holds outside. For the biased community the measure is defined as
$\delta_{\text{\tiny{B}}}\equiv
2\dfrac{L_{\text{\tiny{BB}}}}{L_{\text{\tiny{UB}}}}=\dfrac{\mu_{\text{\tiny{BB}}}}{\mu_{\text{\tiny{BU}}}}$,
(the factor of $2$ in the definition is arbitrary) and for the unbiased
community we use $\delta_{\text{\tiny{U}}}\equiv
2\dfrac{L_{\text{\tiny{UU}}}}{L_{\text{\tiny{UB}}}}=\dfrac{\mu_{\text{\tiny{UU}}}}{\mu_{\text{\tiny{UB}}}}$.
For the biased-independent topology they adopt the values
$\delta_{\text{\tiny{B}}}^{\text{rand}}=\frac{\gamma}{1-\gamma}$ and
$\delta_{\text{\tiny{U}}}^{\text{rand}}=\frac{1-\gamma}{\gamma}$. Therefore,
whenever $\delta_{\text{\tiny{B}}}>\delta_{\text{\tiny{B}}}^{\text{rand}}$,
the biased community is more strongly linked internally than in the case that
links are formed randomly without taking into account the preference of the
agents, and we talk about a closed biased community. Similarly, for
$\delta_{\text{\tiny{B}}}<\delta_{\text{\tiny{B}}}^{\text{rand}}$, the biased
community has less internal links that those corresponding to a complete
random assignment and we speak of an open biased community. A similar
classification of closed or open applies to the community of unbiased agents
for $\delta_{\text{\tiny{U}}}>\delta_{\text{\tiny{U}}}^{\text{rand}}$ or
$\delta_{\text{\tiny{U}}}<\delta_{\text{\tiny{U}}}^{\text{rand}}$,
respectively. As shown in Fig. 10, it turns out that, for fixed $\gamma<1/2$,
the biased community is closed for $\delta>\delta^{\text{rand}}$ and open for
$\delta<\delta^{\text{rand}}$, independently of the model I or II considered.
However, for model I, the unbiased community is open for
$\delta<\delta^{\text{rand}}$ and closed for $\delta>\delta^{\text{rand}}$ and
the opposite behavior for model II: open for $\delta>\delta^{\text{rand}}$ and
closed for $\delta<\delta^{\text{rand}}$. This allows us to plot the phase
diagram of Fig. 11, where, for completeness, we also include the
characteristics of the communities for $\gamma>1/2$, a case not considered
here.

Figure 10:
$\delta_{\text{\tiny{U}}}\equiv\frac{\mu_{\text{\tiny{UU}}}}{\mu_{\text{\tiny{UB}}}}$
and
$\delta_{\text{\tiny{B}}}\equiv\frac{\mu_{\text{\tiny{BB}}}}{\mu_{\text{\tiny{BU}}}}$,
defined in terms of the average number of links between biased and unbiased
communities, rescaled by the values they adopt in a random network where links
are not determined by the preference of the agents,
$\delta_{\text{\tiny{B}}}^{\text{rand}}=\frac{\gamma}{1-\gamma}$,
$\delta_{\text{\tiny{U}}}^{\text{rand}}=\frac{1-\gamma}{\gamma}$, as a
function of $\delta=\frac{\mu_{\text{\tiny{BB}}}}{\mu_{\text{\tiny{UU}}}}$ for
both model I and II, defined in the main text. Values of
$\delta_{\text{\tiny{B}}}$ larger (resp. smaller) than
$\delta_{B}^{\text{rand}}$ indicate a closed (resp. open) biased community,
and similarly for the unbiased community. We consider a minority fraction
$\gamma=0.3$ of biased agents.

Figure 11: Here we plot the phase diagram according to the type, closed
$\delta^{\text{I,II}}_{\text{B,U}}>\delta^{\text{rand}}_{\text{B,U}}$ or open
$\delta^{\text{I,II}}_{\text{B,U}}<\delta^{\text{rand}}_{\text{B,U}}$, of each
of the two communities, biased and unbiased, as a function of the parameters
$\gamma$, the fraction of biased voters, and $\delta$, the ratio of the
average number of connections among only biased voters over the average number
of connections among only unbiased voters. On the left, we plot the phase
diagram for model I. We see that in all of the phase space, the type of
community is the same for both communities. We have a transition of the type
when the line $\delta^{\text{rand}}$ is crossed, and also when the line
$\gamma=0.5$ is crossed. On the contrary, for model II, the type of the two
communities are always opposite, i.e. when the biased community is closed the
unbiased is open and vice versa.

Having established the strategies with which we obtain the biased-dependent
topology of interactions of the network, we proceed to examine the dependence
of the main observables, such as the fixation probability $P_{1}$, the
consensus time $\tau$ and the time to reach consensus to the preferred state
$\tau_{1}$, on the different points of the phase diagram and its dependence
with the closed/open property of the communities. As the results are rather
representative, we focus on a fixed value of $\gamma=0.1$ and vary the
parameter $\delta$ taking values smaller and larger than the purely random
case $\delta^{\text{rand}}=\frac{\gamma}{1-\gamma}=\frac{0.1}{0.9}=0.11$.

The results can be found in Fig 12, where we plot the aforementioned
observables as a function of $\delta$. In this figure we see that increasing
$\delta$ beyond the random value $\delta^{\text{rand}}$ and using the strategy
proposed in Model II to change the biased-dependent topology of the network,
results in a significant reduction in the time to reach consensus, while it
also results in a significant increase in the probability to reach consensus
to the preferred state, when compared to a homogeneous Erdős-Rényi random
network topology of interactions where the establishment of a link is not
influenced by the preference of the nodes. Similarly, decreasing $\delta$
below the random value $\delta^{\text{rand}}$ results in a decrease of the
probability to reach consensus on the preferred state and a slight increase in
the time to consensus (although the results in this latter case are not
conclusive due to the statistical errors). In what concerns the use of the
strategy proposed by model I to change the biased-dependent topology of the
network, the results indicate that neither the fixation probability $P_{1}$,
nor the consensus time $\tau$ show any significant deviation with respect to
the homogeneous Erdős-Rényi case both for $\delta$ larger or smaller than the
random value $\delta^{\text{rand}}$.

The conclusion we draw from these results, is that what matters the most for
the biased group in being more efficient in convincing the rest of the
community to reach consensus on their preferred state, is that these agents,
on average, have more connections with each other, compared to the random
topology of interactions scenario, while at the same time the unbiased group
has, on average, less number of connections with each other, compared to the
random topology of interactions scenario. To be more specific, from the
results concerning the topology variation following the strategy corresponding
to model I, we conclude that when both the biased and unbiased agents decide
to “clash”, i.e. to interact more among each other at the expense of losing
interactions with their peers, then this offers no significant advantage to
the biased agents (blue lines in Fig. 12 for $\delta>\delta^{\text{rand}}$).
The same conclusion arises when both groups decide to interact less among each
other at the expense of winning interactions with their peers (blue lines in
Fig. 12 for $\delta<\delta^{\text{rand}}$). On the contrary, when the biased
agents decide to interact less with the external group, i.e. the unbiased
agents, and more among themselves, and at the same time the unbiased agents
interact less among themselves, then at the long run this gives the biased
agents an advantage compared to the unbiased agents and results in the
preferred opinion being reached faster and more often (green lines in Fig. 12
for $\delta>\delta^{\text{rand}}$). This is because if the group of biased
agents is more compact, i.e. when they interact much more among themselves
compared to their number of interactions with unbiased agents, and at the same
time the unbiased agents are less internally connected, this allows the biased
group to have on average a constant bias towards the preferred state for large
periods of time, bias that can not be reversed by a weak unbiased group. In
this way the biased agents drag the unbiased ones towards their preferred
state.

(a) $\tau(\delta)$

(b) $P_{1}(\delta)$

Figure 12: (a) Time to reach consensus $\tau$ and (b) probability of reaching
consensus in the preferred state $P_{1}$, as a function of the parameter
$\delta=\frac{\mu_{\text{\tiny{BB}}}}{\mu_{\text{\tiny{UU}}}}$, starting from
an initial condition $\sigma=0.5$. The red dashed horizontal lines corresponds
to the theoretical result for an Erdős-Rény network and a biased-independent
topology, Eq.(56) and Eq.(14), respectively, replacing $\beta$ with
$\beta_{\mu}$, as explained in the main text. The vertical line is the value
$\delta^{\text{rand}}=\frac{\gamma}{1-\gamma}$, that separates closed from
open communities (see Fig.11). Blue continuous lines represent numerical
results for model I (note that model I by construction is only defined up to
$\delta=1$), while green dotted lines correspond to model II (see the text for
a detailed explanation of these two models). We can see that model II results
in a significant reduction of the consensus time as $\delta$ is increased,
while it also results in a significant increase in the probability to reach
consensus to the preferred state, when compared to both model I and the case
of a homogeneous Erdős-Rényi random network topology of interactions. The plot
for the consensus time $\tau_{1}$ to the preferred state resembles (a) so it
is omitted. The parameters are $v=0.01$, $\gamma=0.1$, $\mu=5$ and $N=1000$,
and the numerical results have been averaged over $1000$ realizations.

##  6 Conclusions

In this work, we have studied a variation of the standard classical voter
model, where voters have a constant confidence. This means that with a
constant probability $p$ they keep their opinion upon an interaction with
their neighbor, instead of copying. We assumed that the default confidence is
$p=1/2$. However for a fraction $\gamma$ of these voters we assumed that they
are biased towards a fixed opinion, in the sense that their confidence when
changing from a fixed state, say $s_{i}=-1$ to $s_{i}=+1$ is given by
$p=(1+v)/2$ with a bias parameter $v$, while the reverse switch $s_{i}=+1$ to
$s_{i}=-1$ occurs with a confidence equal to $p=(1-v)/2$. We have considered
two distinct scenarios in our studies. First, we assumed that there was no
dependence of the topology of the network on which the dynamics took place and
the type of voters, in which case we studied the model on the complete graph
as well as on an Erdős-Rényi (ER) network. Then we assumed that the topology
of interactions of the two distinct type of voters was indeed dependent on
their type and we examined strategies that the biased voters could follow to
convince faster the rest of the voters to adopt their opinion.

For the biased-independent topology of interactions, we showed that the
fraction of biased agents $\gamma$ and the bias parameter $v$ are not
independent parameters, but rather $\gamma v$ is the only relevant parameter.
Bias breaks the symmetry of the problem. As a consequence, in the
thermodynamic limit the system reaches the absorbing preferred state
$m=1,\rho=0$ in a characteristic time $\frac{1}{\gamma v}$. This is at
variance with the unbiased voter model ($\gamma v\rightarrow 0$) in which
magnetization is conserved and the system remains in a dynamically disordered
state with a finite value of $\rho$. Moving beyond this, we considered finite
size effects on the system, for the same biased-independent topology of
interactions setup. We focused on three observables, the fixation probability
$P_{1}(\sigma)$, or probability to reach the preferred state, as a function of
the initial density $\sigma$ of nodes in state $+1$, the average consensus
time $\tau(\sigma)$ and the average consensus time to the preferred state
$\tau_{1}(\sigma)$ for which we derive analytical expressions. We show that
local effects introduced by the Erdős-Rény network, as compared with the
complete graph case, are accounted by replacing $N$ by
$N_{\mu}=\frac{\mu}{\mu+1}N$. The main effect of the bias is to reduce the
consensus times so that $\tau$ scales as $\log(N)$ while it scales linearly
with $N$ in the limit of no bias ($\gamma v\rightarrow 0$).

Finally, we also have studied the case where the voters lie on Erdős-Rényi
networks of distinct characteristics, i.e. distinct average degrees, depending
on the type of the voter. We refer to this as a biased-dependent topology of
interactions. In this case, we defined as the parameter that quantifies our
deviation from the homogeneous random network, the ratio
$\delta=\dfrac{\mu_{\text{\tiny{BB}}}}{\mu_{\text{\tiny{UU}}}}$, where
$\mu_{\text{\tiny{XY}}}$ represents the average degree of connections between
voters of type X to voters of type Y. With this in mind we identified two
candidate rewiring strategies that keep the total average degree $\mu$
constant and we studied how they affect the dynamics to consensus as a
function of the parameter $\delta$. In strategy I we considered the case of
varying $\mu_{\text{\tiny{BU}}}$ at the expense of $\mu_{\text{\tiny{BB}}}$
and $\mu_{\text{\tiny{UU}}}$. We found that the consensus time and the
probability to reach consensus in the preferred option is indifferent to this
variation. On the contrary in strategy II we considered the scenario of
varying $\mu_{\text{\tiny{BB}}}$ (and accordingly $\mu_{\text{\tiny{UU}}}$)
while keeping $\mu_{\text{\tiny{BU}}}$ constant. We found that increasing
$\mu_{\text{\tiny{BB}}}$ resulted in a significant reduction of the consensus
time $\tau$, as well as to a significant increment to the probability $P_{1}$
of reaching consensus to the preferred state. This evidences that what matters
the most for the members of the biased group in being more efficient in
convincing the rest of the community to reach consensus on their preferred
state, is to increase their internal connections within the group and, at the
same time, decrease the interactions between the members of the unbiased
group. Or, in other words, that a closed community, defined as one where its
members have a higher proportion of inside to outside links that it
corresponds to a completely random Erdős-Rényi network, when put into contact
with an open one, can lead the rest of the system faster and with higher
probability to an agreement to its preferred state.

Acknowledgments. Partial financial support has been received from
MCIN/AEI/10.13039/501100011033 and the Fondo Europeo de Desarrollo Regional
(FEDER, UE) through project PACSS (RTI2018-093732-B-C21) and the María de
Maeztu Program for units of Excellence in R&D, grant MDM-2017-0711.

## Appendix

To find $T_{1}(n)$, the solution of the recurrence equation Eq.(19) satisfying
the boundary conditions $T_{1}(0)=T_{1}(N)=0$ we first define $x(n)$ as

$x(n)=T_{1}(n+1)-T_{1}(n)\Longleftrightarrow T_{1}(n)=\sum_{j=0}^{n-1}x(j).$
(A1)

After substitution in Eq.(19), $x(n)$ is found to satisfy the recurrence
relation

$x(n)=ax(n-1)-\frac{1+a}{n},\,a=\frac{1-\gamma v}{1+\gamma v},$ (A2)

whose solution is

$x(n)=x(0)a^{n}-(1+a)L_{n}(a),\,L_{n}(a)=\sum_{j=1}^{n}\frac{a^{n-j}}{j},\,L_{0}(a)=0.$
(A3)

Replacing in Eq.(A1) and setting the constant $x(0)$ by imposing the boundary
condition $T_{1}(N)=0$, we get

$T_{1}(n,a)=(1+a)\left[\frac{1-a^{n}}{1-a^{N}}\sum_{j=1}^{N-1}L_{j}(a)-\sum_{j=1}^{n-1}L_{j}(a)\right].$
(A4)

The sums of the function $L_{j}(a)$ can be written as

$\displaystyle\sum_{j=1}^{n-1}L_{j}(a)=\sum_{j=1}^{n-1}\sum_{q=1}^{j}\frac{a^{j-q}}{q}=\sum_{q=1}^{n-1}\sum_{j=q}^{n-1}\frac{a^{j-q}}{q}$
$\displaystyle=\sum_{q=1}^{n-1}\frac{1}{q}\frac{1-a^{n-q}}{1-a}=\frac{1}{1-a}\left[H_{n-1}-f(n,a)\right]$
(A5)

where $H_{n}$ (the harmonic function) and $f(n,a)$ are defined by:

$\displaystyle H_{n}$ $\displaystyle=$
$\displaystyle\sum_{q=1}^{n}\frac{1}{q},$ (A6) $\displaystyle f(n,a)$
$\displaystyle=$ $\displaystyle\sum_{q=1}^{n-1}\frac{a^{n-q}}{q},$ (A7)

which replaced in Eq.(A4) yields Eq.(20) in the main text.

It turns out that the function $f(n,a)$ can be written as

$f(n,a)=-\Phi(\frac{1}{a},1,n)-a^{n}\ln\left(\frac{a-1}{a}\right),$ (A8)

in terms of the Lerch transcendent function

$\Phi(x,s,n)=\sum_{j=0}^{\infty}\frac{x^{j}}{\left(j+n\right)^{s}},$ (A9)

a relation valid for all values of $a$ and $n$ as the imaginary parts of the
logarithm and the Lerch transcendent function cancel out for $a<1$.

## References

  * [1] Claudio Castellano, Santo Fortunato, and Vittorio Loreto.  Statistical physics of social dynamics.  Rev. Mod. Phys., 81(2):591–646, may 2009.
  * [2] Jean-Philippe Bouchaud and Jean-Pierre Nadal.  Foreword.  Comptes Rendus Physique, 20(4):241–243, 2019.
  * [3] Maxi San Miguel and Raúl Toral.  Introduction to the Chaos focus issue on the dynamics of social systems.  Chaos, 30(12):120401, dec 2020.
  * [4] Sidney Redner.  Reality-inspired voter models: A mini-review.  Comptes Rendus Physique, 20(4):275–292, 2019.
  * [5] Motoo Kimura and Tomoko Ohta.  The average number of generations until fixation of a mutant gene in a finite population.  Genetics, 61(3):763, 1969.
  * [6] James Franklin Crow, Motoo Kimura, et al.  An introduction to population genetics theory. New York, Evanston and London: Harper & Row, Publishers, 1970.
  * [7] Peter Clifford and Aidan Sudbury.  A model for spatial conflict.  Biometrika, 60(3):581–588, dec 1973.
  * [8] Oscar A Pinto and Miguel A Muñoz.  Quasi-neutral theory of epidemic outbreaks.  PloS ONE, 6(7):e21946, 2011.
  * [9] Xavier Castelló, Víctor M. Eguíluz, and Maxi San Miguel.  Ordering dynamics with two non-excluding options: bilingualism in language competition.  New Journal of Physics, 8:308, 2006.
  * [10] Federico Vazquez, Xavier Castelló, and Maxi San Miguel.  Agent based models of language competition: macroscopic descriptions and order-disorder transitions.  Journal of Statistical Mechanics: Theory and Experiment, 2010(04):P04007, 2010.
  * [11] Laetitia Chapel, Xavier Castelló, Claire Bernard, Guillaume Deffuant, Víctor M Eguíluz, Sophie Martin, and Maxi San Miguel.  Viability and resilience of languages in competition.  Plos ONE, 5(1):e8681, 2010.
  * [12] Marco Patriarca, Xavier Castelló, J. R. Uriarte, Víctor M. Eguíluz, and Maxi. San Miguel.  Modeling two-language competition dynamics.  Advances in Complex Systems, 15(03n04):1250048, 2012.
  * [13] Pavel L. Krapivsky.  Kinetics of a monomer-monomer model of heterogeneous catalysis.  Journal of Physics A: Mathematical and General, 25(22):5831, 1992\.
  * [14] Laurent Frachebourg and Pavel L Krapivsky.  Exact results for kinetics of catalytic reactions.  Physical Review E, 53(4):R3009, 1996.
  * [15] Eli Ben-Naim, Laurent Frachebourg, and Paul L Krapivsky.  Coarsening and persistence in the voter model.  Physical Review E, 53(4):3078, 1996.
  * [16] Ivan Dornic, Hugues Chaté, Jérôme Chave, and Haye Hinrichsen.  Critical coarsening without surface tension: The universality class of the voter model.  Physical Review Letters, 87(4):045701, 2001.
  * [17] Federico Vazquez, Paul L Krapivsky, and Sidney Redner.  Constrained opinion dynamics: Freezing and slow evolution.  Journal of Physics A: Mathematical and General, 36(3):L61, 2003\.
  * [18] Juan Fernández-Gracia, Krzysztof Suchecki, José J. Ramasco, Maxi San Miguel, and Víctor M. Eguíluz.  Is the voter model a model for voters?  Phys. Rev. Lett., 112:158701, Apr 2014.
  * [19] Thomas M Liggett.  Stochastic interacting systems: contact, voter and exclusion processes, volume 324.  Springer Science & Business Media, 2013.
  * [20] Thomas M Liggett.  Continuous time Markov processes: an introduction, volume 113.  American Mathematical Soc., 2010.
  * [21] Richard A Holley and Thomas M Liggett.  Ergodic theorems for weakly interacting infinite systems and the voter model.  The Annals of Probability, 3:643–663, 1975.
  * [22] Claudio Castellano, Vittorio Loreto, Alain Barrat, Federico Cecconi, and Domenico Parisi.  Comparison of voter and glauber ordering dynamics on networks.  Phys. Rev. E, 71:066107, Jun 2005.
  * [23] Pavel L. Krapivsky, Sidney Redner, and Eli (Eli) Ben-Naim.  A kinetic view of Statistical Physics.  Cambridge University Press, Cambridge, 2010.
  * [24] S H Strogatz.  Exploring complex networks.  Nature, 410:268, 2001.
  * [25] Réka Albert and Albert-László Barabási.  Statistical mechanics of complex networks.  Reviews of Modern Physics, 74(1):47, 2002.
  * [26] Claudio Castellano, Daniele Vilone, and Alessandro Vespignani.  Incomplete ordering of the voter model on small-world networks.  Europhysics Letters (EPL), 63(1):153–158, jul 2003.
  * [27] Krzysztof Suchecki, Victor M Eguíluz, and Maxi San Miguel.  Conservation laws for the voter model in complex networks.  Europhysics Letters (EPL), 69(2):228–234, jan 2005.
  * [28] Krzysztof Suchecki, Victor M Eguiluz, and Maxi San Miguel.  Voter model dynamics in complex networks: Role of dimensionality, disorder, and degree distribution.  Phys. Rev. E, 72:36132, 2005.
  * [29] Xavier Castelló, Riitta Toivonen, Víctor M Eguíluz, Jari Saramäki, Kimmo Kaski, and Maxi San Miguel.  Anomalous lifetime distributions and topological traps in ordering dynamics.  EPL (Europhysics Letters), 79(6):66006, 2007.
  * [30] Vishal Sood, Tibor Antal, and Sidney Redner.  Voter models on heterogeneous networks.  Physical Review E, 77(041121):41121, 2008.
  * [31] Federico Vazquez and Víctor M. Eguíluz.  Analytical solution of the voter model on uncorrelated networks.  New Journal of Physics, 10(6):63011, jun 2008.
  * [32] Yuni Iwamasa and Naoki Masuda.  Networks maximizing the consensus time of voter models.  Physical Review E, 90(1):012816, 2014.
  * [33] Gareth J Baxter, Richard A Blythe, and Alan J McKane.  Fixation and consensus times on a network: A unified approach.  Physical Review Letters, 101(25):258701, 2008.
  * [34] Adriana G Moreira and Ronald Dickman.  Critical dynamics of the contact process with quenched disorder.  Physical Review E, 54(4):R3090, 1996.
  * [35] Raffaele Cafiero, Andrea Gabrielli, and Miguel A Muñoz.  Disordered one-dimensional contact process.  Physical Review E, 57(5):5060, 1998.
  * [36] Jef Hooyberghs, Ferenc Iglói, and Carlo Vanderzande.  Strong disorder fixed point in absorbing-state phase transitions.  Physical Review Letters, 90(10):100601, 2003.
  * [37] Géza Ódor and Nóra Menyhárd.  Critical behavior of an even-offspringed branching and annihilating random-walk cellular automaton with spatial disorder.  Physical Review E, 73(3):036130, 2006.
  * [38] Claudio Borile, Amos Maritan, and Miguel A Muñoz.  The effect of quenched disorder in neutral theories.  Journal of Statistical Mechanics: Theory and Experiment, 2013(04):P04032, 2013.
  * [39] Arkadiusz Jȩdrzejewski and Katarzyna Sznajd-Weron.  Person-situation debate revisited: Phase transitions with quenched and annealed disorders.  Entropy, 19(8):415, aug 2017.
  * [40] Serge Galam.  Contrarian deterministic effects on opinion dynamics: “the hung elections scenario”.  Physica A: Statistical Mechanics and its Applications, 333:453–460, 2004.
  * [41] N Masuda.  Voter models with contrarian agents.  Physical Review E, 88:52803, 2013.
  * [42] Mauro Mobilia.  Does a Single Zealot Affect an Infinite Group of Voters?  Physical Review Letters, 91(2):28701, jul 2003.
  * [43] Mauro Mobilia and Ivan T Georgiev.  Voting and catalytic processes with inhomogeneities.  Phys. Rev. E, 71:46102, 2005.
  * [44] Mauro Mobilia, Alexander M Petersen, and Sidney Redner.  On the Role of Zealotry in the Voter Model.  Journal Statistical Mechanics: Theory and Experiment, P08029:1–17, 2007.
  * [45] Serge Galam and Frans Jacobs.  The role of inflexible minorities in the breaking of democratic opinion dynamics.  Physica A: Statistical Mechanics and its Applications, 381:366–376, 2007.
  * [46] Naoki Masuda, Nicolas Gibert, and Sidney Redner.  Heterogeneous voter models.  Physical Review E, 82(1):010103, 2010.
  * [47] D ben Avraham, D Considine, P Meakin, S Redner, and H Takayasu.  Saturation transition in a monomer-monomer model of heterogeneous catalysis.  Journal of Physics A: Mathematical and General, 23(19):4297, 1990\.
  * [48] James W Evans and TR Ray.  Kinetics of the monomer-monomer surface reaction model.  Physical Review E, 47(2):1018, 1993.
  * [49] Krzysztof Suchecki, Andrea Scharnhorst, and Janusz A Holyst.  Highly connected-a recipe for success.  arXiv preprint arXiv:1007.0671, 2010.
  * [50] Nagi Khalil, Maxi San Miguel, and Raul Toral.  Zealots in the mean-field noisy voter model.  Physical Review E, 97(1):012310, 2018.
  * [51] Sidney Redner.  A guide to first-passage processes.  Cambridge university press, 2001.
  * [52] Kurt Jacobs.  Stochastic processes for physicists: understanding noisy systems.  Cambridge University Press, 2010.
  * [53] Frantisek Slanina and Hynek Lavicka.  Analytical results for the Sznajd model of opinion formation.  The European Physical Journal B-Condensed Matter and Complex Systems, 35(2):279–288, 2003.
  * [54] E Pugliese and C Castellano.  Heterogeneous pair approximation for voter models on networks.  EPL (Europhysics Letters), 88(5):58004, dec 2009.
  * [55] Antonio F Peralta, Adrián Carro, Maxi San Miguel, and Raúl Toral.  Stochastic pair approximation treatment of the noisy voter model.  New Journal of Physics, 20(10):103045, 2018.
  * [56] Antonio F Peralta and Raúl Toral.  Binary-state dynamics on complex networks: Stochastic pair approximation and beyond.  Physical Review Research, 2(4):043370, dec 2020.
  * [57] Antonio F Peralta, Adrián Carro, Maxi San Miguel, and Raúl Toral.  Analytical and numerical study of the non-linear noisy voter model on complex networks.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(7):075516, 2018.