Supplementary material for *Characterization of the Nonequilibrium Steady State of a Heteogeneous Nonlinear q-Voter Model with Zealotry*.

Authors: *Andrew Mellor, Mauro Mobilia, R.K.P. Zia*

**Items include:**

- An illustration of the model mechanism.
- Stationary state evolution as a function of $z$.
- Explicit expressions of LGA matrices, $\mathbb{F}$ and $\mathbb{D}$.
- Vorticity and Stream Function: Definitions and Figure 3.
- Comparisons between precision numerical results and simulation data.
- Comparisons between LGA predictions and simulation results.

**Abstract:**

We introduce an heterogeneous nonlinear q-voter model with two types of susceptible voters and zealots, and study its non-equilibrium properties when the population is finite and well mixed. In this two-opinion model, each individual supports one of two parties and is either a susceptible voter of type $q_1$ or $q_2$, or is an inflexible zealot. At each time step, a $q_i$-susceptible voter ($i = 1, 2$) consults a group of $q_i$ neighbors and adopts their opinion if all group members agree, while zealots are inflexible and never change their opinion. We show that this model violates the detailed balance whenever $q_1 \neq q_2$ and has surprisingly rich properties. Here, we focus on the characterization of the model’s non-equilibrium stationary state (NESS) in terms of its probability distribution and currents in the distinct regimes of low and high density of zealotry. We unveil the NESS properties in each of these phases by computing the opinion distribution and the circulation of probability currents, as well as the two-point correlation functions at unequal times (formally related to a “probability angular momentum”). Our analytical calculations obtained in the realm of a linear Gaussian approximation are compared with numerical results.

**Cartoon:** An illustration of the 2qVZ on a generic graph structure. For this letter, our analysis focuses on the case of the complete graph. In this example we have a mixture of zealots (Z), and $q_i$-susceptibles with $q_1=1$ and $q_2=2$. For each iteration, a node is first selected at random (yellow highlight). If they are not a zealot they then choose $q_i$ neighbours (blue arrows) and adopt their opinion if there is a consensus amongst the chosen neighbours.

**Video 1:** Stationary distribution, $P^*$, obtained from the exact numerical solution of the master
equation with $S = 50$, for $Z$ ranging from $15$ to $35$. Dark/light blue encodes a high/low probability.

For each fixed point, there are two matrices $\mathbb{D}$ and $\mathbb{F}$ associated with the noise correlation and linear stability. They are given explicitly by

\begin{eqnarray} \mathbb{D}^{(0)} &=&\frac{s}{8N}\left( \begin{array}{cc} 2 & 0 \\ 0 & 1 \end{array} \right) \\ \mathbb{F}^{(0)} &=&\left( \begin{array}{cc} 1-s & -s \\ -s & \frac{1}{2}-s% \end{array}% \right) \\ \mathbb{D}^{(\pm )} &=&\frac{z}{\left( 1+2z\right) N}\left( \begin{array}{cc} 1-2z & 0 \\ 0 & 2z% \end{array}% \right) \\ \mathbb{F}^{(\pm )} &=&\left( \begin{array}{cc} \frac{1+2z}{2} & \frac{-1+2z}{2} \\ -2z & \frac{1-4z-4z^{2}}{1+2z}% \end{array}% \right) \end{eqnarray}where we used different variables ($s$, $z$) as a reminder that each is small when the system is far from criticality. Note that the eigenvalues of $ \mathbb{F}$, $\lambda _{\pm }=(\mathrm{Tr}\mathbb{F}\pm \Delta )/2$ with $ \Delta =\sqrt{(\mathrm{Tr}\mathbb{F})^{2}-4\det \mathbb{F}}$, are both positive in the regions of interest.

In two-dimensions, a divergenceless vector field (in our case, $\vec{K}^{\ast }$) can be expressed as the curl of a scalar while its curl is also a scalar. Denoting the Levi-Civita symbol as $\varepsilon _{ij}$ ($\varepsilon _{12}=-\varepsilon _{21}=1$,$\varepsilon _{11}=\varepsilon _{22}=0$), these two field are known, respectively, as the stream function and vorticity (repeated indices summed)

\begin{equation}
K_{i}^{\ast }=\varepsilon _{ij}\partial _{j}\psi ^{\ast };\qquad \omega ^{\ast
}=\varepsilon _{ij}\partial _{i}K_{j}^{\ast }
\end{equation}

\begin{eqnarray}
K_{1}^{\ast }\left( n_{1},n_{2}\right) &=&\psi ^{\ast }\left( n_{1}+\frac{1}{2},n_{2}+\frac{1}{2}\right) -\psi ^{\ast }\left( n_{1}+\frac{1}{2},n_{2}-\frac{1}{2}\right) \\
K_{2}^{\ast }\left( n_{1},n_{2}\right) &=&\psi ^{\ast }\left( n_{1}+\frac{1}{2},n_{2}+\frac{1}{2}\right) -\psi ^{\ast }\left( n_{1}-\frac{1}{2},n_{2}+\frac{1}{2}\right) \\
\omega ^{\ast }\left( n_{1}+\frac{1}{2},n_{2}+\frac{1}{2}\right)
&=&K_{2}^{\ast }\left( n_{1}+1,n_{2}\right) -K_{2}^{\ast }\left(
n_{1},n_{2}\right) +K_{1}^{\ast }\left( n_{1},n_{2}\right) -K_{1}^{\ast
}\left( n_{1},n_{2}+1\right).
\end{eqnarray}

Once $P^{\ast }$ is known, it is straightforward to construct $K_{i}^{\ast }(% \vec{n})$ from $W_{i}^{+}(\vec{n})P^{\ast }\left( \vec{n}\right) -W_{i}^{-}(% \vec{n}^{\prime })P^{\ast }\left( \vec{n}^{\prime }\right) $ and the other two fields.

In the LGA, $\vec{K}^{\ast }=\left[ \mathbb{L}/2\right] \vec{\nabla}P^{\ast } $ [16], so that $\psi ^{\ast }\propto P^{\ast }$ can be discerned. Fig.3 shows $\omega^*$, and $\psi^*$ for the same system as Fig.1, giving additional information about the NESS. In Fig.3(c,d), we see that $\psi ^{\ast }$ strongly resembles $P^{\ast }$, the difference is a measure of both the quality of the LGA and the finite size effects. Meanwhile, $\omega ^{\ast }\propto -\nabla ^{2}P^{\ast }$ carries the information of curvatures in $P^{\ast }$, leading us to expect regions of positive/negative vorticity near/far-from the peak(s) of $P^{\ast }$. This property is clearly seen in Fig.3(a,b). Of course, the sign of $\omega ^{\ast }$ reveals whether the microscopic net-current loop goes around the plaquette clockwise or anti-clockwise.

**Figure 3:** Vorticity and stream function obtained from the exact numerical solution of the master equation in the high/low zealotry phase (left/right panels) with $S=30$. Top row: Vorticity $\omega^*$, as function of $\vec{x}=\vec{n}/N$; red/blue encodes a positive/negative vorticity corresponding to anti-clockwise/clockwise whirls. Bottom row: Stream function $\psi^*$; darker regions are more negative. Other parameters are $Z=20$ and $N=100$ in (a) and (c), and $Z=13$ and $N=86$ in (b) and (d).

Since the numerical results from the master equation are available only for two small systems ($Z=13,20;S=30$), we quote simulation data only for these cases. For the $z < z_{c}$ case, the walker is found to visit both wells frequently. Thus, the
simulated, double peaked distribution is a reasonably good representation of
the true stationary state. The averages $\left\langle \vec{x}\right\rangle$
are exactly $s/2$ by symmetry, confirmed by simulations within statistical
error. As for the equal time correlations $\mathbb{C}$, the exact results
are $\left( C_{11},C_{12},C_{22}\right) \cong \left( 4.31,5.46,9.77\right)
\times 10^{-3}$ which are entirely consistent with the simulation values $\left( 4.30,5.43,9.69\right) \times 10^{-3}$.
Similarly, the agreement for
the $z>z_{c}$ case is equally impressive: $\left( 1.74,1.67,3.54\right) \times 10^{-3}$
*vs.* $\left( 1.74,1.66,3.53\right) \times 10^{-3}$*.* Turning to $\mathcal{L%
}$, we only examined the $Z=20$ system (for reasons stated in the main
text). Using the time series $n_{i}\left( T\right) $ to compute $%
\left\langle n_{i}n_{j}\right\rangle _{T=1}$ and dividing that by $dt=1/N$,
we obtain $\mathcal{L}=2.15\times 10^{-4}$, which is also the value
predicted by the numerical $P^{\ast },\vec{K}^{\ast }$. Finally, let us
remark that this value of $\mathcal{L}$ is completely consistent with $%
2.16\times 10^{-4}$ found by fitting the data to Eq. (7).

As true
stationary states, the exact averages $\left\langle \vec{x}\right\rangle $
are $s/2$ by symmetry, regardless of the location of the peak(s). Thus, for $%
z>z_{c}$, the MFA-LGA result for $\left\langle \vec{x}\right\rangle $
trivially agrees with simulations. The first non-trivial comparisons with
the LGA are the (truncated) correlations. For the small $Z=20,S=30$ system,
we see that the LGA predictions $\left( C_{11},C_{12},C_{22}\right) \cong
\left( 2.25,2.75,6.00\right) \times 10^{-3}$ are qualitatively correct.
Similarly, at $5.00\times 10^{-4}$, $\mathcal{L}$ is qualitatively captured
by the LGA. In the low zealotry case ($Z=13,S=30$), the exact $\left\langle
\vec{n}\right\rangle $ is still $15$, despite the peaks being located at $%
\left( 9,5\right) $ and its complement, $\left( 21,25\right) $ (Fig.1b). Comparing to the MFA-LGA fixed points, we find, say, $\vec{n}^{\left( -\right) }=N\vec{x}^{\left( -\right) }\cong \left(10.99,7.52\right) $ to be tolerably close to $\left( 9,5\right) $. Since the
LGA represents an expansion around only *one* of the two peaks, we
refrain from comparing its predictions for $\mathbb{C}$ or $\mathcal{L}$ to
simulations here.

Turning to the large system ($Z=400,S=1000$), where we
expect better agreement between simulation results (where the walker never
leaves one well) and MFA-LGA predictions (around one fixed point). Indeed,
for the averages $\left\langle \vec{x}\right\rangle $, we find quite good
agreement: ($0.120,0.0736$) *vs.* ($0.119,0.0714$). Though not as
impressive, the agreement for $\mathbb{C}$ is still quite reasonable: $%
\left( 8.79,8.01,13.4\right) \times 10^{-5}$ from simulations and $\left(
8.49,7.35,12.0\right) \times 10^{-5}$ from LGA. More remarkably, computing
the correlation $\left\langle \xi _{i}\xi _{j}\right\rangle _{\tau }$ at a
single attempt step, we find the elements of $\mathbb{FC}=\mathbb{D+L}/2$
agreeing with the LGA predictions to around 1% or better. Finally, let us
just note that the agreement between simulation results of the high zealotry
system ($Z=800,S=1000$) and the LGA predictions are much better for the
equal time correlation $\mathbb{C}$, but slightly worse for the single time
step $\mathbb{FC}$.