and the definition in Eq. (S165). We find

$\left<w(\mathbf{x})\right>=\frac{1}{2}\left[\eta-\alpha S(\mathbf{x})-\gamma
S^{(1)}(\mathbf{x})\right].$ (S168)

We also have

$\displaystyle\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}$ $\displaystyle
W^{(\ell)}(\mathbf{x},\mathbf{y})\left<s(\mathbf{x})s(\mathbf{y})\left[w(\mathbf{x})+w(\mathbf{y})\right]\right>$
$\displaystyle=\eta\nu^{(\ell)}(\mathbf{x})S^{(\ell)}(\mathbf{x})-\frac{\alpha}{2}\Bigg{[}\nu^{(\ell)}(\mathbf{x})S(\mathbf{x})+\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{x},\mathbf{y})S(\mathbf{y})\Bigg{]}$
$\displaystyle\qquad-\frac{\gamma}{2}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{x},\mathbf{y})\left[\frac{1}{\nu^{(1)}(\mathbf{y})}\sum_{\mathbf{z}\in\mathbf{y}^{(1)}}\left<s(\mathbf{x})s(\mathbf{z})\right>+\frac{1}{\nu^{(1)}(\mathbf{x})}\sum_{\mathbf{z}\in\mathbf{x}^{(1)}}\left<s(\mathbf{z})s(\mathbf{y})\right>\right].$
(S169)

Eq. (S166) then becomes

$\displaystyle\frac{\differential S^{(\ell)}(\mathbf{x})}{\differential t}$
$\displaystyle=\frac{2W^{(\ell)}(\mathbf{x},\mathbf{x})}{\nu^{(\ell)}(\mathbf{x})}\left[\eta-\alpha
S(\mathbf{x})-\gamma S^{(1)}(\mathbf{x})\right]$ $\displaystyle\qquad-2\eta
S^{(\ell)}(\mathbf{x})+\alpha\Bigg{[}S(\mathbf{x})+\frac{1}{\nu^{(\ell)}(\mathbf{x})}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{x},\mathbf{y})S(\mathbf{y})\Bigg{]}$
$\displaystyle\qquad+\frac{\gamma}{\nu^{(\ell)}(\mathbf{x})}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{x},\mathbf{y})\left[\frac{1}{\nu^{(1)}(\mathbf{y})}\sum_{\mathbf{z}\in\mathbf{y}^{(1)}}\left<s(\mathbf{x})s(\mathbf{z})\right>+\frac{1}{\nu^{(1)}(\mathbf{x})}\sum_{\mathbf{z}\in\mathbf{x}^{(1)}}\left<s(\mathbf{z})s(\mathbf{y})\right>\right].$
(S170)

In the steady-state, the derivative on the left-hand side vanishes. We can
then divide through by $\gamma$ so that the model parameters only appear in
the following combinations,

$\displaystyle\frac{\eta}{\gamma}$ $\displaystyle=1+\tau,$ (S171a)
$\displaystyle\frac{\alpha}{\gamma}$
$\displaystyle=\tau\left\langle{m}\right\rangle_{\rm st},$ (S171b)

where $\left\langle{m}\right\rangle_{\rm st}$ and $\tau$ are as in Eqs. (S24)
and (S46) respectively. As we will show below, this means that the steady-
state interface density can also be determined solely from
$\left\langle{m}\right\rangle_{\rm st}$ and $\tau$.

We now sum both sides of Eq. (S170) over all nodes $\mathbf{x}$ in the network
to derive an equation for the overall correlation $C^{(\ell)}$, defined as

$C^{(\ell)}\equiv\frac{1}{N}\sum_{\mathbf{x}}S^{(\ell)}(\mathbf{x}).$ (S172)

We find,

$\displaystyle\frac{\differential C^{(\ell)}}{\differential t}$
$\displaystyle=2\frac{1}{N}\sum_{\mathbf{x}}\left\\{\frac{W^{(\ell)}(\mathbf{x},\mathbf{x})}{\nu^{(\ell)}(\mathbf{x})}\left[\eta-\alpha
S(\mathbf{x})-\gamma S^{(1)}(\mathbf{x})\right]\right\\}-2\eta C^{(\ell)}$
$\displaystyle\qquad+\alpha\Bigg{[}\left\langle{m}\right\rangle+\frac{1}{N}\sum_{\mathbf{x}}\frac{1}{\nu^{(\ell)}(\mathbf{x})}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{x},\mathbf{y})S(\mathbf{y})\Bigg{]}$
$\displaystyle\qquad+\gamma\frac{1}{N}\sum_{\mathbf{x}}\frac{1}{\nu^{(\ell)}(\mathbf{x})}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{x},\mathbf{y})\Bigg{[}\frac{1}{\nu^{(1)}(\mathbf{y})}\sum_{\mathbf{z}\in\mathbf{y}^{(1)}}\left<s(\mathbf{x})s(\mathbf{z})\right>$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\frac{1}{\nu^{(1)}(\mathbf{x})}\sum_{\mathbf{z}\in\mathbf{x}^{(1)}}\left<s(\mathbf{z})s(\mathbf{y})\right>\Bigg{]}.$
(S173)

To proceed we introduce a number of shorthands. We write

$\displaystyle\Omega^{(\ell)}_{\mathbf{x}}$
$\displaystyle\equiv\frac{W^{(\ell)}(\mathbf{x},\mathbf{x})}{\nu^{(\ell)}(\mathbf{x})},$
(S174a) $\displaystyle\chi^{(\ell)}_{\mathbf{x}}$
$\displaystyle\equiv\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}\frac{W^{(\ell)}(\mathbf{y},\mathbf{x})}{\nu^{(\ell)}(\mathbf{y})},$
(S174b) $\displaystyle\psi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}$
$\displaystyle\equiv\sum_{\mathbf{z}\in\mathbf{x}^{(\ell)}\cap\mathbf{y}^{(1)}}\left[\frac{W^{(\ell)}(\mathbf{x},\mathbf{z})W^{(1)}(\mathbf{z},\mathbf{y})}{\nu^{(\ell)}(\mathbf{x})\nu^{(1)}(\mathbf{z})}\right],$
(S174c) $\displaystyle\phi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}$
$\displaystyle\equiv\sum_{\mathbf{z}\in\mathbf{y}^{(\ell)}\cap\,\mathbf{x}^{(1)}}\left[\frac{W^{(\ell)}(\mathbf{z},\mathbf{y})W^{(1)}(\mathbf{z},\mathbf{x})}{\nu^{(\ell)}(\mathbf{z})\nu^{(1)}(\mathbf{z})}\right].$
(S174d)

In practice these coefficients can be evaluated from the $N\times N$ adjacency
matrix $\mathbf{A}$ of the network. We here recall that the
$\mathbf{x}\mathbf{y}$ element of the $\ell$th power of the adjacency gives
the number of possible length $\ell$ walks starting at $\mathbf{x}$ and ending
at $\mathbf{y}$, denoted $W^{(\ell)}(\mathbf{x},\mathbf{y})$, see Eq. (S132).
The sum of such elements over all $\mathbf{y}$ gives the total number of walks
of length $\ell$ starting from $\mathbf{x}$, denoted
$\nu^{(\ell)}(\mathbf{x})$, see Eq. (S134). Using these, we have

$\displaystyle\Omega^{(\ell)}_{\mathbf{x}}$
$\displaystyle=\frac{(A^{\ell})_{\mathbf{x}\mathbf{x}}}{\sum_{\mathbf{y}}(A^{\ell})_{\mathbf{x}\mathbf{y}}},$
(S175a) $\displaystyle\chi^{(\ell)}_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{y}}\frac{(A^{\ell})_{\mathbf{y}\mathbf{x}}}{\sum_{\mathbf{z}}(A^{\ell})_{\mathbf{y}\mathbf{z}}},$
(S175b) $\displaystyle\psi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}$
$\displaystyle=\frac{1}{\sum_{\mathbf{z}}(A^{\ell})_{\mathbf{x}\mathbf{z}}}\left[\sum_{\mathbf{z}}(A^{\ell})_{\mathbf{x}\mathbf{z}}\frac{1}{\sum_{\mathbf{u}}A_{\mathbf{z}\mathbf{u}}}A_{\mathbf{z}\mathbf{y}}\right],$
(S175c) $\displaystyle\phi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}$
$\displaystyle=\sum_{\mathbf{z}}\left[\frac{(A^{\ell})_{\mathbf{z}\mathbf{y}}A_{\mathbf{z}\mathbf{x}}}{\sum_{\mathbf{u}}(A^{k})_{\mathbf{z}\mathbf{u}}\sum_{\mathbf{u}}A_{\mathbf{z}\mathbf{u}}}\right].$
(S175d)

We now proceed to write Eq. (S173) in more compact form. Using Eq. (S136) we
observe for the term on the second line,

$\displaystyle\frac{1}{N}\sum_{\mathbf{x}}\frac{1}{\nu^{(\ell)}(\mathbf{x})}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{x},\mathbf{y})S(\mathbf{y})$
$\displaystyle=\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}\left(\frac{W^{(\ell)}(\mathbf{x},\mathbf{y})}{\nu^{(\ell)}(\mathbf{x})}S(\mathbf{y})\right)$
$\displaystyle=\frac{1}{N}\sum_{\mathbf{x}}\left(\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}\frac{W^{(\ell)}(\mathbf{y},\mathbf{x})}{\nu^{(\ell)}(\mathbf{y})}\right)S(\mathbf{x})$
$\displaystyle=\frac{1}{N}\sum_{\mathbf{x}}\chi_{\mathbf{x}}^{(\ell)}S(\mathbf{x}).$
(S176)

For the first term on the third line,

$\displaystyle\frac{1}{N}\sum_{\mathbf{x}}\frac{1}{\nu^{(\ell)}(\mathbf{x})}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{x},\mathbf{y})\frac{1}{\nu^{(1)}(\mathbf{y})}\sum_{\mathbf{z}\in\mathbf{y}^{(1)}}\left<s(\mathbf{x})s(\mathbf{z})\right>$
$\displaystyle\qquad=\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}\sum_{\mathbf{z}\in\mathbf{y}^{(1)}}\frac{W^{(\ell)}(\mathbf{x},\mathbf{y})}{\nu^{(\ell)}(\mathbf{x})\nu^{(1)}(\mathbf{y})}\left\langle{s(\mathbf{x})s(\mathbf{z})}\right\rangle$
$\displaystyle\qquad=\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell+1)}}\left(\sum_{\mathbf{z}\in\mathbf{x}^{(\ell)}\cap\,\mathbf{y}^{(1)}}\frac{W^{(\ell)}(\mathbf{x},\mathbf{z})}{\nu^{(\ell)}(\mathbf{x})\nu^{(1)}(\mathbf{z})}\right)\left<s(\mathbf{x})s(\mathbf{y})\right>$
$\displaystyle\qquad=\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell+1)}}\left(\sum_{\mathbf{z}\in\mathbf{x}^{(\ell)}\cap\,\mathbf{y}^{(1)}}\frac{W^{(\ell)}(\mathbf{x},\mathbf{z})W^{(1)}(\mathbf{z},\mathbf{y})}{\nu^{(\ell)}(\mathbf{x})\nu^{(1)}(\mathbf{z})}\right)\left<s(\mathbf{x})s(\mathbf{y})\right>$
$\displaystyle\qquad=\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell+1)}}\psi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}\left<s(\mathbf{x})s(\mathbf{y})\right>,$
(S177)

where we have first used Eq. (S137) in going from the second to third line,
and then introduced $W^{(1)}(\mathbf{z},\mathbf{y})=1$ for
$\mathbf{z}\in\mathbf{y}^{(1)}$, before finally using the definition of
$\psi^{(\ell)}_{\mathbf{x},\mathbf{y}}$ from Eq. (S174c). In a similar manner
one shows that

$\displaystyle\frac{1}{N}\sum_{\mathbf{x}}\frac{1}{\nu^{(\ell)}(\mathbf{x})}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{x},\mathbf{y})\frac{1}{\nu^{(1)}(\mathbf{x})}\sum_{\mathbf{z}\in\mathbf{x}^{(1)}}\left<s(\mathbf{z})s(\mathbf{y})\right>$
$\displaystyle\qquad=\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}\sum_{\mathbf{z}\in\mathbf{y}^{(1)}}\frac{W^{(\ell)}(\mathbf{x},\mathbf{y})}{\nu^{(\ell)}(\mathbf{x})\nu^{(1)}(\mathbf{x})}\left<s(\mathbf{z})s(\mathbf{y})\right>$
$\displaystyle\qquad=\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell+1)}}\left(\sum_{\mathbf{z}\in\mathbf{y}^{(\ell)}\cap\,\mathbf{x}^{(1)}}\frac{W^{(\ell)}(\mathbf{z},\mathbf{y})}{\nu^{(\ell)}(\mathbf{z})\nu^{(1)}(\mathbf{z})}\right)\left<s(\mathbf{x})s(\mathbf{y})\right>$
$\displaystyle\qquad=\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell+1)}}\left(\sum_{\mathbf{z}\in\mathbf{y}^{(\ell)}\cap\,\mathbf{x}^{(1)}}\frac{W^{(1)}(\mathbf{x},\mathbf{z})W^{(\ell)}(\mathbf{z},\mathbf{y})}{\nu^{(\ell)}(\mathbf{z})\nu^{(1)}(\mathbf{z})}\right)\left<s(\mathbf{x})s(\mathbf{y})\right>$
$\displaystyle\qquad=\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell+1)}}\phi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}\left<s(\mathbf{x})s(\mathbf{y})\right>,$
(S178)

where we have first used Eq. (S138) in going from the second to third line,
then introduced $W^{(1)}(\mathbf{x},\mathbf{z})=1$ for
$\mathbf{x}\in\mathbf{z}^{(1)}$, and finally used the definition of
$\phi^{(\ell)}_{\mathbf{x},\mathbf{y}}$ from Eq. (S174d). Eq. (S173) then
becomes

$\displaystyle\frac{\differential C^{(\ell)}}{\differential t}=$
$\displaystyle
2\eta\frac{1}{N}\sum_{\mathbf{x}}\Omega^{(\ell)}_{\mathbf{x}}-\alpha\frac{1}{N}\sum_{\mathbf{x}}\left[2\Omega^{(\ell)}_{\mathbf{x}}-\chi^{(\ell)}_{\mathbf{x}}\right]S(\mathbf{x})+\alpha\left\langle{m}\right\rangle-2\eta
C^{(\ell)}$
$\displaystyle\qquad-2\gamma\frac{1}{N}\sum_{\mathbf{x}}\Omega^{(\ell)}_{\mathbf{x}}S^{(1)}(\mathbf{x})+\gamma\frac{1}{N}\sum_{\mathbf{x}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell+1)}}\left[\psi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}+\phi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}\right]\left<s(\mathbf{x})s(\mathbf{y})\right>.$
(S179)

Several of the terms in the last expression take the form of weighted
magnetisations, whose steady-state we can determine from Eq. (S164). However
the two final terms are problematic and we have not attempted to solve this
equation in full generality. Progress is possible for what we will call
‘homogeneous’ networks, as discussed in the next section.

Alternatively we can choose model parameters such that $\gamma=0$. However,
this would mean that $\tau\to\pm\infty$ and the steady-state configuration
would always be random which is not of particular interest (see also Sec. S8).

###  S6.5 Homogeneous networks

####  S6.5.1 Homogeneity assumptions

We define networks that fulfill the following two properties as homogeneous:

  1. (1)

the total number walks of length $\ell$ starting and ending at any node
$\mathbf{x}$ is the same for all $\mathbf{x}$,
$W^{(\ell)}(\mathbf{x},\mathbf{x})\equiv W^{(\ell)}$;

  2. (2)

the total number of walks of length $\ell$ starting at $\mathbf{x}$ (ending at
any point) is also the same for all $\mathbf{x}$,
$\nu^{(\ell)}(\mathbf{x})\equiv\nu^{(\ell)}$.

There are several types of network that fall under this category. For example,
infinite regular lattices in any dimension or Bethe lattices (infinite regular
trees) have these properties. Individual realisations of degree-regular
networks do not fulfill (1) and (2), but these conditions hold as averages
over the ensemble of degree-regular networks, i.e. the expected value of
$W^{(\ell)}(\mathbf{x},\mathbf{x})$ is the same for all nodes $\mathbf{x}$,
and similarly for the expected value of $\nu^{(\ell)}(\mathbf{x})$.

We now proceed making these assumptions, writing $W^{(\ell)}$ and
$\nu^{(\ell)}$ for the common values of all
$W^{(\ell)}(\mathbf{x},\mathbf{x})$ and $\nu^{(\ell)}(\mathbf{x})$
respectively. The assumptions imply in particular that all nodes in the
network have the same degree, $\nu^{(1)}$. We note that we then have
$\nu^{(\ell+1)}=\nu^{(\ell)}\nu^{(1)}$, and from this we find
$\nu^{(\ell)}=(\nu^{(1)})^{\ell}$.

Noting again that $W^{(1)}(\mathbf{z},\mathbf{y})=A_{\mathbf{z}\mathbf{y}}$ is
equal to one if and only if $\mathbf{z}$ and $\mathbf{y}$ are nearest
neighbours, and zero otherwise, Eqs. (S174a)-(S174d) can then be simplified as
follows:

$\displaystyle\Omega^{(\ell)}_{\mathbf{x}}$
$\displaystyle=\frac{W^{(\ell)}(\mathbf{x},\mathbf{x})}{\nu^{(\ell)}(\mathbf{x})}=\frac{W^{(\ell)}}{\nu^{(\ell)}}\equiv\Omega^{(\ell)},\hskip
28.45274pt\forall\mathbf{x}$ (S180a) $\displaystyle\chi^{(\ell)}_{\mathbf{x}}$
$\displaystyle=\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}\frac{W^{(\ell)}(\mathbf{y},\mathbf{x})}{\nu^{(\ell)}(\mathbf{y})}=\frac{1}{\nu^{(\ell)}}\sum_{\mathbf{y}\in\mathbf{x}^{(\ell)}}W^{(\ell)}(\mathbf{y},\mathbf{x})=1,\hskip
28.45274pt\forall\mathbf{x}$ (S180b)
$\displaystyle\psi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}$
$\displaystyle=\frac{1}{\nu^{(\ell)}(\mathbf{x})}\left[\sum_{\mathbf{z}\in\mathbf{x}^{(\ell)}\cap\mathbf{y}^{(1)}}W^{(\ell)}(\mathbf{x},\mathbf{z})\frac{1}{\nu^{(1)}(\mathbf{z})}W^{(1)}(\mathbf{z},\mathbf{y})\right]$
$\displaystyle=\frac{1}{\nu^{(\ell)}\nu^{(1)}}\sum_{\mathbf{z}\in\mathbf{y}^{(1)}}W^{(\ell)}(\mathbf{x},\mathbf{z})W^{(1)}(\mathbf{z},\mathbf{y})$
$\displaystyle=\frac{W^{(\ell+1)}(\mathbf{x},\mathbf{y})}{\nu^{(\ell+1)}},\hskip
28.45274pt\forall\mathbf{x},\mathbf{y}$ (S180c)
$\displaystyle\phi^{(\ell+1)}_{\mathbf{x}\mathbf{y}}$
$\displaystyle=\sum_{\mathbf{z}\in\mathbf{y}^{(\ell)}\cap\mathbf{x}^{(1)}}\left[\frac{1}{\nu^{(\ell)}(\mathbf{z})}W^{(\ell)}(\mathbf{z},\mathbf{y})\frac{1}{\nu^{(1)}(\mathbf{z})}W^{(1)}(\mathbf{z},\mathbf{x})\right]$
$\displaystyle=\frac{1}{\nu^{(\ell+1)}}\sum_{\mathbf{z}\in\mathbf{x}^{(1)}}W^{(\ell)}(\mathbf{z},\mathbf{y})$
$\displaystyle=\frac{W^{(\ell+1)}(\mathbf{x},\mathbf{y})}{\nu^{(\ell+1)}}.\hskip
28.45274pt\forall\mathbf{x},\mathbf{y}$ (S180d)

####  S6.5.2 Interface density for homogeneous networks, and resulting
simplifications

Assuming that the underlying network is homogeneous (in the sense as defined
above), Eq. (S179) reduces to

$\frac{\differential C^{(\ell)}}{\differential
t}=2\Omega^{(\ell)}\Big{[}\eta-\alpha\left\langle{m}\right\rangle-\gamma
C^{(1)}\Big{]}-2\Big{[}\eta
C^{(\ell)}-\alpha\left\langle{m}\right\rangle-\gamma C^{(\ell+1)}\Big{]}.$
(S181)

In the steady-state we then have the following relation:

$\eta C^{(\ell)}_{\rm st}-\gamma C^{(\ell+1)}_{\rm
st}=\Omega^{(\ell)}\Big{[}\eta-\alpha\left\langle{m}\right\rangle_{\rm
st}-\gamma C^{(1)}_{\rm st}\Big{]}+\alpha\left\langle{m}\right\rangle_{\rm
st}.$ (S182)

The interface density, the quantity we are aiming to calculate, is given by

$\left\langle{\sigma}\right\rangle=\frac{1-C^{(1)}}{2},$ (S183)

thus we need to determine $C^{(1)}_{\rm st}$. However setting $\ell=0$ in Eq.
(S182) is not useful, as all terms involving $C^{(1)}_{\rm st}$ then cancel
out because of $\Omega^{(\ell=0)}=1$. To see the latter we note that there is
one walk of length zero steps starting at $\mathbf{x}$, and this walk
trivially also ends at $\mathbf{x}$. Thus Eq. (S180a) evaluates to one.

Instead, to proceed we multiply both sides of Eq. (S182) by
$\gamma^{\ell}/\eta^{\ell}$, and then sum over $\ell$ from zero to infinity.
We find

$\displaystyle\sum_{\ell=0}^{\infty}\Bigg{(}\frac{\gamma^{\ell}}{\eta^{\ell-1}}C^{(\ell)}_{\rm
st}$ $\displaystyle-\frac{\gamma^{\ell+1}}{\eta^{\ell}}C^{(\ell+1)}_{\rm
st}\Bigg{)}-\alpha\left\langle{m}\right\rangle_{\rm
st}\sum_{\ell=0}^{\infty}\frac{\gamma^{\ell}}{\eta^{\ell}}=\Big{[}\eta-\alpha\left\langle{m}\right\rangle_{\rm
st}-\gamma C^{(1)}_{\rm
st}\Big{]}\sum_{\ell=0}^{\infty}\Omega^{(\ell)}\frac{\gamma^{\ell}}{\eta^{\ell}}.$
(S184)

The first term on the left-hand side is a telescopic sum and simplifies to

$\sum_{\ell=0}^{\infty}\Bigg{(}\frac{\gamma^{\ell}}{\eta^{\ell-1}}C^{(\ell)}_{\rm
st}-\frac{\gamma^{\ell+1}}{\eta^{\ell}}C^{(\ell+1)}_{\rm st}\Bigg{)}=\eta
C^{(0)}_{\rm
st}-\lim_{\ell\to\infty}\frac{\gamma^{\ell+1}}{\eta^{\ell}}C^{(\ell+1)}_{\rm
st}.$ (S185)

We can use analogous arguments to that at the end of Sec. S6.2, except
considering $C^{(\ell)}_{\rm st}$ instead of
$\left\langle{m_{n}}\right\rangle_{\rm st}$, to perform the geometric sum in
Eq. (S184) and show that the limit in Eq. (S185) vanishes. Eq. (S184) then
becomes

$\eta C^{(0)}_{\rm st}-\frac{\alpha\left\langle{m}\right\rangle_{\rm
st}}{1-\frac{\gamma}{\eta}}=\Big{[}\eta-\alpha\left\langle{m}\right\rangle_{\rm
st}-\gamma C^{(1)}_{\rm
st}\Big{]}\sum_{\ell=0}^{\infty}\Omega^{(\ell)}\left(\frac{1}{1+\tau}\right)^{\ell}.$
(S186)

Using the fact that $C^{(0)}_{\rm
st}=\left\langle{s(\mathbf{x})s(\mathbf{x})}\right\rangle_{\rm st}=1$, and the
identities from Eqs. (S171a) and (S171b), and Eq. (S183), we find

$\displaystyle\left\langle{\sigma}\right\rangle_{\rm st}$
$\displaystyle=\frac{1}{2}H_{\rm
HG}(\tau,\mathbf{A})\left(1-\left\langle{m}\right\rangle^{2}_{\rm st}\right),$
(S187)

with

$H_{\rm
HG}(\tau)=\frac{1+\tau}{\sum_{\ell=0}^{\infty}\Omega^{(\ell)}\left(\frac{1}{1+\tau}\right)^{\ell}}-\tau.$
(S188)

We note that calculation of $\Omega^{(\ell)}$ requires the adjacency matrix of
the network $\mathbf{A}$. The subscript ‘HG’ stands for ‘homogeneous graph’,
and indicates that the relation holds only when the homogeneity assumptions
apply. See Sec. S8 for detailed discussion on $\tau$ and $H(\tau)$.

Eq. (S188) is not closed form but the coefficients $\Omega^{(\ell)}$ can, in
principle, be obtained through direct enumeration of walks on a given network.
We note that this only needs to be performed once for a given network, to
obtain the function $H_{\rm HG}(\tau)$. We discuss this further in Sec.
S6.5.7.

We now move on to use Eq. (S188) for a number of special cases (finite
complete networks, hyper-cubic lattices, and Bethe lattices) where closed form
solutions can be found.

####  S6.5.3 First special case: Finite complete networks

We consider a complete network with $N$ nodes. For any fixed node, we write
$a_{N}(\ell)$ for the number of walks starting and ending at that node (this
number will be the same for all nodes, given the all-to-all connectivity).

Figure S7: Illustration of closed walks on a complete network (see text).

One then has the following recursion

$a_{N}(\ell)=(N-1)^{\ell-1}-a_{N}(\ell-1),$ (S189)

with $a_{N}(0)\equiv 1$.

This recursion can be derived as follows (see Fig. S7), keeping in mind that
the walker cannot stay put in any step of the walk. Consider a fixed node
$\mathbf{x}$, and a closed walk of length $\ell$ starting and ending at
$\mathbf{x}$. The walk consists of a sequence of $\ell+1$ nodes in total,
which we label $\mathbf{x}_{0},\mathbf{x}_{1},\dots,\mathbf{x}_{\ell}$. The
first and the last of these are given by $\mathbf{x}$
($\mathbf{x}_{0}=\mathbf{x}_{\ell}=\mathbf{x}$), shown in blue in Fig. S7. Now
focus on node $\mathbf{x}_{1}$. This node can be any node in the network other
than node $\mathbf{x}_{0}$ (the walker must make a hop in the first step).
Similarly, $\mathbf{x}_{2}$ can be any node other than $\mathbf{x}_{1}$, and
so on. If we follow this though until (and including) node
$\mathbf{x}_{\ell-1}$, we find a total of $(N-1)^{\ell-1}$ possible sequences
for the nodes $\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{\ell-1}$ such
that no node is the same as the previous one in the sequence. This produces
the term $(N-1)^{\ell-1}$ in Eq. (S189). There is however another constraint,
and not all of these sequences produce valid walks of length $\ell$ starting
and ending at $\mathbf{x}$. In a valid walk, node $\mathbf{x}_{\ell-1}$ (shown
in red in the figure) must not be node $\mathbf{x}$, as the walker must make a
hop in the $\ell$th step. We therefore need to subtract the number of all
sequences for which this is the case. This number is $a_{N}(\ell-1)$, as the
sequence
$\mathbf{x}_{0}=\mathbf{x},\mathbf{x}_{1},\dots,\mathbf{x}_{\ell-1}=\mathbf{x}$
is then a closed walk of length $\ell-1$ starting and ending at $\mathbf{x}$.
This produces the second term in Eq. (S189). Thus, we have demonstrated the
validity of the recursion.

One can directly show that (see e.g. sequences A109499, A109500, A109501 in
[30])

$a_{N}(\ell)=\frac{(N-1)^{\ell}+(N-1)(-1)^{\ell}}{N}$ (S190)

fulfills the recursion in Eq. (S189), as well as the condition $a_{N}(0)=1$.

Given that further the number of all walks of length $\ell$ starting at
$\mathbf{x}$ on a complete network of size $N$ is given by
$\nu^{(\ell)}(\mathbf{x})=(N-1)^{\ell}$, we have

$\Omega^{(\ell)}=\frac{W^{(\ell)}}{\nu^{(\ell)}(\mathbf{x})}=\frac{a_{N}(\ell)}{(N-1)^{\ell}}=\frac{(N-1)^{\ell}+(N-1)(-1)^{\ell}}{N(N-1)^{\ell}}.$
(S191)

The generating function for the sequence $\Omega^{(\ell)}$ can be directly
evaluated to be

$\sum_{\ell=0}^{\infty}\Omega^{(\ell)}z^{\ell}=\frac{1-\left(\frac{N-2}{N-1}\right)z}{1-\left(\frac{N-2}{N-1}\right)z-\left(\frac{1}{N-1}\right)z^{2}}.$
(S192)

Thus Eq. (S187) for the finite complete network is

$\displaystyle\left\langle{\sigma}\right\rangle_{\rm st}$
$\displaystyle=\frac{1}{2}\left[\frac{(N-1)\tau}{1+(N-1)\tau}\right]\left(1-\left\langle{m}\right\rangle^{2}_{\rm
st}\right)$ $\displaystyle\equiv\frac{1}{2}H_{\rm
CG}(\tau,N)\left(1-\left\langle{m}\right\rangle^{2}_{\rm st}\right),$ (S193)

(the subscript ‘CG’ stands for ‘complete graph’). We note that in the limit
$N\rightarrow\infty$ we recover Eq. (S29).

####  S6.5.4 Second special case: Infinite hyper-cubic lattice in
$d$-dimensions

We next consider the case of a hyper-cubic lattice in $d$-dimensions. We then
have $\Omega^{(\ell)}=0$ for odd values of $\ell$, as it is impossible to
return to a given point in an odd number of steps on a hyper-cubic lattice.
Thus we focus on even $\ell=2k$. So looking at Eq. (S188) we need to calculate
$\Omega^{(\ell)}\to\Omega^{(2k)}_{d}$, where the $d$ subscript denotes the
dimension of the hyper-cubic lattice. From Eq. (S180a) we have

$\Omega^{(2k)}_{d}=\frac{W^{(2k)}_{d}}{\nu^{(2k)}_{d}},$ (S194)

On a $d$-dimensional hyper-cubic lattice, the total number of walks of length
$2k$ starting at a given point is

$\nu^{(2k)}_{d}=(2d)^{2k},$ (S195)

given that the walker can choose from one of $2d$ directions at each step. The
number of walks starting and ending at the same point is

$W_{d}^{(2k)}=\sum_{n_{1}+...+n_{d}=k}\frac{(2k)!}{\prod_{i=1}^{d}(n_{i}!)^{2}}.$
(S196)

In this expression the $n_{i}$ ($i=1,\dots,d$) are integer numbers; $n_{i}$
indicates how many steps in total the walker takes in the positive $i$
direction along the walk. Given that the walker returns to the starting point,
$n_{i}$ is also the number of steps taken in the negative $i$ direction. Since
a total of $2k$ steps must be made, we have $\sum_{i}2n_{i}=2k$, i.e.
$\sum_{i}n_{i}=k$. The term inside the sum in Eq. (S196) is the number of
distinct walks that can be generated for a given choice of
$n_{1},\dots,n_{d}$. The multinomial coefficient is the number of permutations
of the $n_{1}$ forward steps in direction $1$, the $n_{1}$ backward steps in
direction $1$, the $n_{2}$ forward steps in direction $2$, the $n_{2}$
backward steps in direction $2$ and so on, keeping in mind that the total
number of steps is $2k$.

Using this we can construct the following so-called ‘exponential generating
function’ [note the factors of $1/(2k)!$ which is the difference between an
exponential generating function and an ordinary generating function],

$\sum_{k=0}^{\infty}W_{d}^{(2k)}\frac{x^{2k}}{(2k)!}=I_{0}(2x)^{d},$ (S197)

where $I_{0}(\cdot)$ is the modified Bessel function of the first kind. A
proof of this can be found in Sec. S10.1. A Laplace-Borel transform, which we
define as

$\mathcal{L}_{B}[f](x)=\int_{0}^{\infty}e^{-t}f(tx)dt,$ (S198)

can be used to convert an exponential generating function of the type in Eq.
(S197) to an ordinary generating function [31, pp. 374, 566], [32, p. 750] or
[33, p. 2]. Therefore we have

$\sum_{k=0}^{\infty}W_{d}^{(2k)}x^{2k}=\int_{0}^{\infty}e^{-t}[I_{0}(2xt)]^{d}\differential
t.$ (S199)

What we want is a summation of the form in Eq. (S188), so we need to
incorporate the $\nu^{(2k)}_{d}=(2d)^{2k}$ factor which can be done by
applying the transform $x\to\frac{x}{2d}$ to Eq. (S199),

$\sum_{k=0}^{\infty}W_{d}^{(2k)}\frac{x^{2k}}{(2d)^{2k}}=\int_{0}^{\infty}e^{-t}\left[I_{0}\left(\frac{xt}{d}\right)\right]^{d}\differential
t.$ (S200)

With $x=\frac{1}{1+\tau}$, the summation on the LHS is exactly that in Eq.
(S188) for an infinite hyper-cubic $d$-dimensional lattice. So we are left
with solving the integral on the RHS.

To proceed, it is now useful to make the substitution $xt=z$,

$\displaystyle\int_{0}^{\infty}e^{-t}\left[I_{0}\left(\frac{xt}{d}\right)\right]^{d}\differential
t$
$\displaystyle=\frac{1}{x}\int_{0}^{\infty}e^{-\frac{z}{x}}\left[I_{0}\left(\frac{z}{d}\right)\right]^{d}\differential
z$
$\displaystyle=\frac{1}{x}\mathcal{L}\left(\left[I_{0}\left(\frac{z}{d}\right)\right]^{d}\right).$
(S201)

where $\mathcal{L}$ is the standard Laplace transform defined as

$\mathcal{L}[f](s)=\int_{0}^{\infty}e^{-sz}f(z)\differential z,$ (S202)

where in this case $s=\frac{1}{x}$. We can now make use of the following
general identity [34, p. 346],

$\displaystyle\mathcal{L}\left(z^{\lambda}\prod_{i=1}^{d}I_{\nu_{i}}(a_{i}z)\right)=$
$\displaystyle\frac{\Gamma(\lambda+\nu+1)}{2^{\nu}p^{\lambda+\nu+1}}\left[\prod_{i=1}^{d}\frac{a_{i}^{\nu_{i}}}{\Gamma(\nu_{i}+1)}\right]$
$\displaystyle\times
F_{C}^{(d)}\left(\frac{\lambda+\nu+1}{2},\frac{\lambda+\nu+2}{2};\nu_{1},...,\nu_{d}+1;\frac{a_{1}^{2}}{p^{2}},...,\frac{a_{d}^{2}}{p^{2}}\right),$
(S203)

where $F_{C}^{(d)}$ is a Lauricella hypergeometric series [35, p. 113][36],

$F_{C}^{(d)}(a,b;c_{1},...,c_{d};y_{1},...,y_{d})=\sum_{i_{1},...,i_{d}=0}^{\infty}\frac{(a)_{i_{1}+...+i_{d}}(b)_{i_{1}+...+i_{d}}}{(c_{1})_{i_{1}}...(c_{d})_{i_{d}}i_{1}!...i_{d}!}y_{1}^{i_{1}}...y_{d}^{i_{d}},$
(S204)

with

$(q)_{i}\equiv q(q+1)...(q+i-1),$ (S205)

for $i=0,1,2,3,\dots$, referred to as the Pochhammer symbol or rising
factorial.

We note that the identity in Eq. (S203) as printed in [34] contains a
typographical error, for further details see also [37]. This error has been
corrected in Eq. (S203).

We now substitute Eq. (S203) in the special case $\lambda=0$, $\nu_{i}=0$ and
$a_{i}=\frac{1}{d}$ into Eq. (S201) and find

$\displaystyle\int_{0}^{\infty}e^{-t}\left[I_{0}\left(\frac{zt}{d}\right)\right]^{d}\differential
t$
$\displaystyle=F_{C}^{(d)}\left(\frac{1}{2},1;1,...,1;\frac{1}{\left[d(1+\tau)\right]^{2}},...,\frac{1}{\left[d(1+\tau)\right]^{2}}\right),$
(S206)

where $z=\frac{1}{1+\tau}$. With this, we find the following expression for
the interface density on a $d$-dimensional hyper-cubic lattice,

$\left<\sigma\right>_{\rm
st}=\frac{1}{2}\left[\frac{1+\tau}{F_{C}^{(d)}\left(\frac{1}{2},1;1,...,1;\frac{1}{\left[d(1+\tau)\right]^{2}},...,\frac{1}{\left[d(1+\tau)\right]^{2}}\right)}-\tau\right]\left(1-\left\langle{m}\right\rangle^{2}_{\rm
st}\right).$ (S207)

####  S6.5.5 Application to infinite two-dimensional square lattice

We have already derived an analytical expression for $H(\tau)$ for infinite 2D
square lattice in Sec. S5 of this Supplement [see Eq. (S125)]. However as a
test of our results for the homogeneous network and hyper-cubic lattice in
$d$-dimensions we will re-derive it, with reference to several points of the
derivation.

First method: Closed form expression for $\Omega^{(2k)}$
We can directly evaluate the summation in Eq. (S188) by first determining a
closed form expression for $\Omega^{(\ell)}$. On an infinite square lattice
only even length walks can return us to the initial position, so it is really
$\Omega^{(2k)}$ we wish to calculate.

The number of closed walks of length $2k$ can be calculated with the following
combinatorial argument. If we take $p$ steps North, we must then also take $p$
steps South. There are $\binom{2k}{2p}$ ways to choose these $2p$ steps from
$2k$ and $\binom{2p}{p}$ ways to organise the North and South steps. From the
remaining steps, half of them must be East and the other half must be West,
there are $\binom{2k-2p}{k-p}$ ways of arranging these steps. Thus the total
number of closed walks of length $2k$ is

$\displaystyle\sum_{p=0}^{k}\binom{2k}{2p}\binom{2p}{p}\binom{2(k-p)}{k-p}$
$\displaystyle=\sum_{p=0}^{k}\frac{(2k)!}{(2p)!(2(k-p))!}\cdot\frac{(2p)!}{p!^{2}}\cdot\frac{(2(k-p))!}{(k-p)!^{2}}$
$\displaystyle=\frac{(2k)!}{k!^{2}}\sum_{p=0}^{k}\frac{k!^{2}}{p!^{2}(k-p)!^{2}}$
$\displaystyle=\binom{2k}{k}\sum_{p=0}^{k}\binom{k}{p}\binom{k}{k-p}$
$\displaystyle=\binom{2k}{k}^{2},$ (S208)

where in the last step we use Vandermonde’s identity. We also note that the
total number of length $2k$ walks on a 2D square lattice is
$\nu^{(2k)}=4^{2k}$. Thus, letting $x=\frac{1}{1+\tau}$, the summation in Eq.
(S188) becomes

$\sum_{k=0}^{\infty}\Omega^{(2k)}x^{2k}=\sum_{k=0}^{\infty}\binom{2k}{k}^{2}\frac{x^{2k}}{4^{2k}}=\frac{2}{\pi}K(x),$
(S209)

where we have used the power series for the complete elliptical integral of
the first kind [38, (1.3.6)]. With this, Eq. (S188) reduces exactly to Eq.
(S125).

Second method: Relation of elliptic integral to modified Bessel function
The second way of verifying the special case $d=2$ is by evaluation of Eq.
(S200) for $d=2$.

We are looking to prove that the Laplace-Borel transformation [see Eq. (S198)]
of the square of modified Bessel functions of the first kind is proportional
to the complete elliptical integral of the first kind i.e.

$\int_{0}^{\infty}e^{-t}I_{0}\left(\frac{xt}{2}\right)^{2}\differential
t=\frac{2}{\pi}K(x).$ (S210)

First, make the change of variables $xt=z$ to convert the Laplace-Borel
transform into a regular Laplace transform,

$\displaystyle\int_{0}^{\infty}e^{-t}I_{0}\left(\frac{xt}{2}\right)^{2}\differential
t$
$\displaystyle=\frac{1}{x}\int_{0}^{\infty}e^{-pz}I_{0}\left(\frac{z}{2}\right)^{2}\differential
z,$ (S211)

where $p=\frac{1}{x}$. From Sec. S10.1 we find the exponential generating
function for the square of modified Bessel functions of the first kind is

$\displaystyle I_{0}(2y)^{2}$
$\displaystyle=\sum_{k=0}^{\infty}W_{2}^{(2k)}\frac{y^{2k}}{(2k)!}$
$\displaystyle=\sum_{k=0}^{\infty}{\binom{2k}{k}}^{2}\frac{y^{2k}}{(2k)!},$
(S212)

where $W_{2}^{(2k)}$ can be found from Eq. (S196) with $d=2$. Substituting Eq.
(S212) into Eq. (S211) gives,

$\displaystyle\frac{1}{x}\int_{0}^{\infty}e^{-pz}I_{0}\left(\frac{z}{2}\right)^{2}\differential
z$
$\displaystyle=\frac{1}{x}\int_{0}^{\infty}e^{-pz}\sum_{k=0}^{\infty}{\binom{2k}{k}}^{2}\left(\frac{z}{4}\right)^{2k}\frac{1}{(2k)!}\differential
z$
$\displaystyle=\frac{1}{x}\sum_{k=0}^{\infty}{\binom{2k}{k}}^{2}\frac{1}{(2k)!}\frac{1}{4^{2k}}\int_{0}^{\infty}e^{-pz}{z}^{2k}\differential
z.$ (S213)

The integral here is simply the Laplace transform of ${z}^{2k}$ given by [34,
p. 15],

$\displaystyle\int_{0}^{\infty}e^{-pz}{z}^{2k}\differential z$
$\displaystyle=\frac{\Gamma(2k+1)}{p^{2k+1}}$ $\displaystyle=x^{2k+1}(2k)!,$
(S214)

where we recall $p=\frac{1}{x}$ and where we have used the definition of the
Gamma function for integer $n$ [39, p. 255],

$\Gamma(n)=(n-1)!.$ (S215)

Putting this back into Eq. (S213) we finally have

$\displaystyle\int_{0}^{\infty}e^{-t}I_{0}\left(\frac{xt}{2}\right)^{2}\differential
t$ $\displaystyle=\sum_{k=0}^{\infty}{\binom{2k}{k}}^{2}\frac{x^{2k}}{4^{2k}}$
$\displaystyle=\frac{2}{\pi}K(x),$ (S216)

where we have used the power series for the complete elliptical integral of
the first kind [38, (1.3.6)]. Thus we have shown that Eq. (S200) evaluated at
$d=2$, which is exactly the relevant summation appearing in Eq. (S188) for the
2D infinite square lattice, is proportional to $K(x)$, thus Eq. (S188) reduces
to Eq. (S125).

Third proof: Relation of elliptic integral to Lauricella hyper-geometric
series
The final method is to evaluate Eq. (S207), for the $d$-dimensional lattice,
for the case $d=2$.

For $d=2$, the following Lauricella series corresponds to the Appell
hypergeometric series [36],

$\displaystyle F_{C}^{(2)}\Bigg{(}\frac{1}{2},1;1,1;\frac{1}{4(1+\tau)^{2}}$
$\displaystyle,\frac{1}{4(1+\tau)^{2}}\Bigg{)}=F_{4}\left(\frac{1}{2},1;1,1;\frac{1}{4(1+\tau)^{2}},\frac{1}{4(1+\tau)^{2}}\right).$
(S217)

An Appell hypergeometric series of this form can be written as a generalised
hypergeometric series [40, Eq. (37)]

$F_{4}\left(\frac{1}{2},1;1,1;\frac{1}{4(1+\tau)^{2}},\frac{1}{4(1+\tau)^{2}}\right)={}_{4}F_{3}\left(\frac{1}{2},1,1,\frac{1}{2};1,1,1;\frac{1}{(1+\tau)^{2}}\right).$
(S218)

The generalised hypergeometric function ${}_{p}F_{q}$ is defined as [41, p.
439],

${}_{p}F_{q}(a_{1},...,a_{p};b_{1},...,b_{q};x)=\sum_{k=0}^{\infty}\frac{(a_{1})_{k},...,(a_{p})_{k}}{(b_{1})_{k},...,(b_{q})_{k}}\frac{x^{k}}{k!}.$
(S219)

Using this we can show that

${}_{4}F_{3}\left(\frac{1}{2},1,1,\frac{1}{2};1,1,1;x^{2}\right)={}_{2}F_{1}\left(\frac{1}{2},\frac{1}{2},1;x^{2}\right).$
(S220)

Finally it is known that [42, p. 51]

${}_{2}F_{1}\left(\frac{1}{2},\frac{1}{2},1;x^{2}\right)\equiv\frac{2}{\pi}K(x).$
(S221)

Thus we have shown that Eq. (S207) reduces to Eq. (S124) for the 2D infinite
square lattice.

####  S6.5.6 Third special case: Infinite Bethe lattice

The infinite regular tree, or Bethe lattice, is a network in which every node
has degree $n$ and which has a cascading tree-like structure such that there
are no loops or leaf nodes.

On such trees, only walks of even length can return to the starting point. So
we have $W^{(2k+1)}=0$ for all $k\in\mathbb{Z}_{\geq 0}$, and thus only need
to perform the summation in Eq. (S188) over even $k$. Additionally, it can be
seen by inspection that $\nu^{(\ell)}=n^{\ell}$. Ultimately, the summation in
Eq. (S188) reduces to

$T\left(z\right)=\sum_{k=0}^{\infty}\mathcal{T}^{(k)}z^{k},$ (S222)

for

$z=\frac{1}{n^{2}(1+\tau)^{2}},$ (S223)

where $\mathcal{T}^{(k)}\equiv W^{(2k)}$ is the number of self returning walks
of length $2k$ from any point on the network.

This means that we need the generating function for the sequence
$\mathcal{T}^{(k)}$. To calculate this we follow the approach in [43] and
first consider an infinite rooted tree, in which the root has degree $r$ and
all other nodes have degree $c+1$. We write $\mathbb{T}$ for this tree.

A closed rooted walk is defined as a walk which starts and ends at the root
node. We will refer to the nearest neighbours of the root as the ‘children’.

Any closed rooted walk of length $2k$ in $\mathbb{T}$ can be broken down as
follows:

  1. 1.

The walker starts at the root, and takes one step from the root to one of its
$r$ nearest neighbours (i.e., a child).

  2. 2.

A closed walk of length $2i$, avoiding the root, is undertaken, starting at
the child node in step 1, and bringing the walker back to the same child.

  3. 3.

The walker returns to the root node, having taken $2(i+1)$ steps in total.

  4. 4.

A closed rooted walk of length $2k-2(i+1)$ can now be undertaken.

We define ${\cal R}_{r}^{(k)}$ as the number of self-returning walks of length
$2k$ from the root node, and ${\cal C}_{c+1}^{(k)}$ as the number of self-
returning root-avoiding walks of length $2k$ from one of the children. The
subscripts indicate the degree of the root ($r$) and that of all other nodes
($c+1$) respectively. These quantities are specific to the tree $\mathbb{T}$.
However, the number of self-returning root-avoiding walks starting at a child
will not depend on the degree $r$ of the root.

From the above breakdown on rooted walks of length $2k$ we can obtain the
following recursive relation:

${\cal R}_{r}^{(n+1)}=r\sum_{i=0}^{n}{\cal R}_{r}^{(i)}{\cal
C}_{c+1}^{(n-i)}.$ (S224)

We next define

$R_{r}(z)\equiv\sum_{k=0}^{\infty}{\cal R}^{(k)}_{r}z^{k}$ (S225)

as the generating function for self-returning walks in the tree from the root
node which has degree $r$. Similarly,

$C_{c+1}(z)=\sum_{k=0}^{\infty}{\cal C}^{(k)}_{c+1}z^{k}$ (S226)

is the generating function for self-returning walks in the tree from a child
of the root node and avoiding the root.

With these definitions we have

$\displaystyle R_{r}(z)$ $\displaystyle\equiv$
$\displaystyle\sum_{k=0}^{\infty}{\cal R}^{(k)}_{r}z^{k}$ (S227)
$\displaystyle=$ $\displaystyle 1+\sum_{k=1}^{\infty}{\cal R}^{(k)}_{r}z^{k}$
$\displaystyle=$ $\displaystyle 1+r\sum_{k=1}^{\infty}\sum_{i=0}^{k-1}{\cal
R}_{r}^{(i)}{\cal C}_{c+1}^{(k-1-i)}z^{k}$ $\displaystyle=$ $\displaystyle
1+rz\sum_{k=0}^{\infty}\sum_{i=0}^{k}{\cal R}_{r}^{(i)}{\cal
C}_{c+1}^{(k-i)}z^{k}$ $\displaystyle=$ $\displaystyle
1+rzR_{r}(z)C_{c+1}(z).$

The last step can be seen by writing out the first few terms of the series for
$R_{r}(z)$ and $C_{c+1}(z)$ in Eqs. (S225) and (S226).

We can determine $R_{r}(z)$ from Eq. (S227) if we can calculate $C_{c+1}(z)$.
To do this, we note again that ${\cal C}_{c+1}^{(k)}$ does not depend on the
degree of the root. Hence $C_{c+1}(z)$ is not a function of the root-degree
either. With this in mind, we now consider a tree $\tilde{\mathbb{T}}$ with
root of degree $c$ and where, as before, all other nodes have degree $c+1$.
Setting $r=c$ in Eq. (S227) we then have

$\tilde{R}_{c}(z)=1+cz\tilde{R}_{c}(z)\tilde{C}_{c+1}(z),$ (S228)

where the tildes indicate that these are generating functions for the tree
$\tilde{\mathbb{T}}$. The only difference between the trees $\mathbb{T}$ and
$\tilde{\mathbb{T}}$ is the degree of the root. Hence
$\tilde{C}_{c+1}(z)=C_{c+1}(z)$.

Next we realise that $\tilde{R}_{c}(z)=\tilde{C}_{c+1}(z)$. This is because
$\tilde{C}_{c+1}(z)$ is the generating function of self-returning walks
starting from a child and avoiding the root on the tree $\tilde{\mathbb{T}}$.
For the purpose of counting such walks the child effectively has degree $c$
(the same degree as the root of $\tilde{\mathbb{T}}$). Hence on the tree
$\tilde{\mathbb{T}}$, the number of self-returning root-avoiding walks of a
given length and starting from a child is the same as the number of self-
returning walks of this length starting at the root.

Using $\tilde{R}_{c}(z)=\tilde{C}_{c+1}(z)$ and
$\tilde{C}_{c+1}(z)=C_{c+1}(z)$ Eq. (S228) thus becomes

$C_{c+1}(z)=1+czC^{2}_{c+1}(z).$ (S229)

We can then solve for $C_{c+1}(z)$ and find the physically relevant solution
to be

$C_{c+1}(z)=\frac{1-\sqrt{1-4zc}}{2zc}$ (S230)

(this is the solution ensuring that $\lim_{z\to 0}C_{c+1}(z)=1$, i.e. there is
exactly one walk of length zero). It is interesting to note that this
generating function is closely related to the generating function for the
Catalan numbers [44, 45].

Having now determined $C_{c+1}(z)$ we can substitute this back into Eq. (S227)
to determine $R_{r}(z)$,

$R_{r}(z)=\frac{2c}{2c-r+r\sqrt{1-4zc}}.$ (S231)

For the special case of the $n$-regular tree (i.e. the tree where all nodes
have the same degree $n$) we simply set $n=r=c+1$ in Eq. (S231). This gives
the generating function for self-returning walks on an infinite $n$-regular
tree, which we define as $T(z)$ in Eq. (S222),

$T(z)=\frac{n-1}{n-2+n\sqrt{1-4(n-1)z}}.$ (S232)

Then with $z$ as in Eq. (S223) we have a closed-form expression for the
summation in Eq. (S222):

$T(z)=\frac{2(n-1)}{n-2+n\sqrt{1-4(n-1)\left(\frac{1}{n(1+\tau)}\right)^{2}}}.$
(S233)

With the replacement $n\to\mu$, Eqs. (S187) and (S233) combine to give

$\displaystyle\left\langle{\sigma}\right\rangle_{\rm st}$
$\displaystyle=\frac{1}{2}\left[\frac{\mu-2-\mu\tau+\sqrt{(\tau+1)^{2}\mu^{2}-4(\mu-1)}}{2(\mu-1)}\right]\left(1-\left\langle{m}\right\rangle^{2}_{\rm
st}\right)$
$\displaystyle\equiv\frac{1}{2}H_{\text{BL}}(\tau)\left(1-\left\langle{m}\right\rangle^{2}_{\rm
st}\right),$ (S234)

(where ‘BL’ stands for ‘Bethe lattice’). We note that this is the same result
that we derived in Sec. S3 using the pair approximation for infinite
uncorrelated networks. In the current section we did not have to make any
approximations about spin correlations as we did for the pair approximation.
The only assumption that we made was that the dynamics was being taking place
on a particular topology, hence this result is exact for Bethe lattices. This
is not surprising as the pair approximation is known to be exact on Bethe
lattices [46, p; 346] as there are no loops.

####  S6.5.7 Calculation of $\Omega^{(\ell)}$ and validity of the method

For networks that are not homogeneous in the above sense, we can still apply
the approach based on network walks [Eq. (S188)] as an approximation. We
replace $\Omega^{(\ell)}$ with its average over all nodes,

$\Omega^{(\ell)}\rightarrow\frac{1}{N}\sum_{\mathbf{x}}\Omega^{(\ell)}_{\mathbf{x}}\equiv\overline{\Omega^{(\ell)}}.$
(S235)

It is difficult to find this object in closed form. However, for many networks
we find that $\overline{\Omega^{(\ell)}}\rightarrow\frac{1}{N}$ for large
$\ell$ (where $N$ is the size of the graph). Thus we compute
$\overline{\Omega^{(\ell)}}$ directly from the adjacency matrix for $\ell$
below some cutoff $L$, and use $\overline{\Omega^{(\ell)}}\approx\frac{1}{N}$
for $\ell>L$. This approximation should be tested for each network on an
individual basis, from which an appropriate choice for $L$ should be made.
Based on this then we can write, using the geometric series,

$\sum_{\ell=0}^{\infty}\overline{\Omega^{(\ell)}}x^{\ell}\approx\frac{1}{N}\frac{x^{L+1}}{1-x}+\sum_{\ell=0}^{L}\overline{\Omega^{\ell}}x^{\ell}.$
(S236)

In Fig. S8 we show the convergence of $\overline{\Omega^{(\ell)}}$ to
$\frac{1}{N}$ for Barabási–Albert networks of differen size. These networks
are not homogeneous. In this case it would be suitable to take $L\approx 30$.

Figure S8: Demonstration of $\overline{\Omega^{(\ell)}}$ converging to
$\frac{1}{N}$ for increasing $\ell$ for difference size Barabási–Albert
networks with $4$ new links established upon adding a new node. The extent of
the vertical axis shown is from $1/N-\delta$ to $1/N+\delta$, with
$\delta=1\times 10^{-4}$.

Although the approaches based on network walks is only exact on homogeneous
networks, we find that it serves as a satisfactory approximation for networks
in which the variation of $\Omega^{\ell}_{\mathbf{x}}$ across nodes is small.
In Fig. S9 we show the frequencies at which length $\ell$ walks (top) and
length $\ell$ loops (bottom) appear as we vary $\ell$ for degree-regular,
Erdös–Rényi, and Barabási–Albert networks.

Figure S9: Histograms of the frequency of length $\ell$ walks (top) and loops
(bottom) from each node in three different networks: degree-regular (left
column) with $\mu=8$, Erdös–Rényi (middle column) with $p=8/N$ and
Barabási–Albert (right column) with with $4$ new links established upon adding
a new node. All networks have $N=2,500$.

On degree-regular networks with degree $k$, the number of walks of a given
length $\ell$ starting from any node is always $k^{\ell}$, as confirmed in
Fig. S9(a1), (a2) and (a3). The number of loops can vary slightly from node to
node (see panels d1-f1) as the network is not exactly homogeneous, but the
distributions are still relatively tight, compared to the two other types of
network shown.

For Erdös–Rényi networks there is more variance in the distributions, and for
Barabási–Albert networks there is significantly more. Thus we can anticipate
the approximation to be somewhat accurate on Erdös–Rényi networks, and less so
on Barabási–Albert networks.

In Fig. S10 we show how $H_{\rm HG}(\tau)$ from Eq. (S188) varies as a
function of $\tau$ for different degree-regular (left), Erdös–Rényi (middle)
and Barabási–Albert (right) networks. The theory works very well for degree-
regular networks and is a good approximation for Erdös-–Rényi networks as
expected because the variance of their degree distributions is small. For
Barabási–Albert there is greater deviation.

Figure S10: Plots of $H(\tau)=\left\langle{\sigma}\right\rangle_{\rm
st}/[\frac{1}{2}(1-\left\langle{m}\right\rangle_{\rm st}^{2})]$, which
characterises the scattering, as a function of the noise-to-order ratio,
$\tau$, for networks of size $N$. Degree-regular with $\mu=8$ (left),
Erdös–Rényi with $p=8/N$ (middle), Barabási–Albert with $4$ new links
established upon adding a new node (right). For all three types of network we
have $\mu\approx 8$ (deviations are due to the finite size of the graphs). The
dashed lines are the results from our theory, Eq. (S188). Each marker is the
result from averaging 1000 independent simulations in the steady-state at a
particular $\tau$ value. Model parameters for a particular $\tau$ can be
generated via the algorithm detailed in Sec. S10.2.

####  S6.5.8 Small disconnected networks

It is interesting to note that the result $H_{\rm HG}(\tau)$ applies to
homogeneous networks of any size, for example see Fig. 3 where we show the
result working for an $N=10$ degree-regular network. Further, the theory does
not assume that the network is connected, in fact many low-degree Erdös–Rényi
networks are disconnected.

Fig. S11 demonstrates that the theory works well on a network of $N=4$ with
only two connected pairs of nodes.

The pathological case of a system in which all nodes are isolated from one
another is not handled by the theory. In this case there are no links, and
thus the fraction of active links is undefined.

Figure S11: Plot of $H(\tau)=\left\langle{\sigma}\right\rangle_{\rm
st}/[\frac{1}{2}(1-\left\langle{m}\right\rangle_{\rm st}^{2})]$, which
characterises the scattering, against the noise-to-order ratio $\tau$ for a
network of size $N=4$ containing two connected pairs. Solid line is $H_{\rm
HG}(\tau)$ from Eq. (S188). Markers are the result from averaging over
$1\times 10^{5}$ independent simulations in the steady-state at a specific
value of $\tau$. Model parameters for specific values of $\tau$ are generated
via the algorithm detailed in Sec. S10.2.

##  S7 Stochastic pair approximation

The stochastic pair approximation (SPA) is an extension of the regular pair
approximation discussed in Sec. S3. The SPA targets large but finite
uncorrelated networks. This was originally proposed in [17] (pp 7-16 in
particular), which studied the noisy voter model. Our model is more general,
and reduces to the noisy voter model in the limit of no copying errors in the
horizontal process, $h_{+}^{\prime}=h_{-}^{\prime}=0$, and when the vertical
process is symmetric, $v_{+}=v_{-}$. We here construct the SPA for the more
general setup of our model. More precisely we develop what is called ‘S2PA’ in
[17] but here we refer to it as simply ‘SPA’.

Unlike in the main paper and in the remaining sections of this supplement we
do not use $-1$ and $+1$ for the node states. Instead we denote these states
by $0$ and $1$. This is to keep the notation as close to that of [17].

###  S7.1 Master equation for reduced set of degrees of freedom

To begin we consider the reduced set of variables
$\left\\{\mathbf{N}_{1},\mathcal{L}\right\\}$ where $\mathbf{N}_{1}$ is the
vector

$\mathbf{N}_{1}=\left(\mathbf{N}_{1,k_{\text{min}}},...,\mathbf{N}_{1,k_{\text{max}}}\right),$
(S237)

whose components $N_{1,k}$ are the number of spin-up nodes with degree $k$.
The notation $N_{0,k}$ will be used to represent the number of spin-down nodes
with degree $k$, and thus the total number of nodes of degree $k$ is
$N_{k}=N_{0,k}+N_{1,k}$. $\mathcal{L}$ is the number of active links. We also
define $n_{i}\in[0,k_{i}]$ as the number of neighbours of node $i$ in the up-
state.

When a node $i$ flips from the down to the up-state, we have the process
$N_{1,k_{i}}\rightarrow N_{1,k_{i}}+1$. The number of active links connected
to $i$ before flipping is $n_{i}$, after the flip it will be $k_{i}-n_{i}$,
thus a change of $k_{i}-2n_{i}$ in the number of active links connected to
$i$. Analogously, when a node flips from up to down, we have
$N_{1,k_{i}}\rightarrow N_{1,k_{i}}-1$ and the total number of active links
changes as $\mathcal{L}\rightarrow\mathcal{L}-k_{i}+2n_{i}$.

We then define the probability distribution of
$\left\\{\mathbf{N}_{1},\mathcal{L}\right\\}$ at time $t$ as
$P(\mathbf{N}_{1},\mathcal{L};t)$ and we can construct the following master
equation:

$\displaystyle\frac{\partial}{\partial
t}P(\mathbf{N}_{1},\mathcal{L};t)=\sum_{k}\sum_{n=0}^{k}\Big{\\{}$
$\displaystyle\big{(}E_{N_{1,k}}^{+}E_{\mathcal{L}}^{k-2n}-1\big{)}\big{[}W^{(-,k,n)}P(\mathbf{N}_{1},\mathcal{L};t)\big{]}$
$\displaystyle\qquad+\big{(}E_{N_{1,k}}^{-}E_{\mathcal{L}}^{-k+2n}-1\big{)}\big{[}W^{(+,k,n)}P(\mathbf{N}_{1},\mathcal{L};t)\big{]}\Big{\\}},$
(S238)

where $W^{(+,k,n)}$ and $W^{(-,k,n)}$ are the effective rates of the two
processes described above. The $E$’s are step operators which act on functions
of $N_{1,k}$ or ${\cal L}$ as follows

$\displaystyle E_{N_{1,k}}^{+}f(N_{1,k})$ $\displaystyle\rightarrow
f(N_{1,k}+1),$ (S239a) $\displaystyle E_{\mathcal{L}}^{k-2n}f(\mathcal{L})$
$\displaystyle\rightarrow f(\mathcal{L}+k-2n),$ (S239b) $\displaystyle
E_{N_{1,k}}^{-}f(N_{1,k})$ $\displaystyle\rightarrow f(N_{1,k}+1),$ (S239c)
$\displaystyle E_{\mathcal{L}}^{-k+2n}f(\mathcal{L})$
$\displaystyle\rightarrow f(\mathcal{L}-k+2n).$ (S239d)

Next we define the degree specific magnetisation $m_{k}$, i.e. the
magnetisation among the nodes of degree $k$. The interface density is denoted
by $\sigma$. These two quantities can be written in terms of the variables
$\left\\{\mathbf{N}_{1},\mathcal{L}\right\\}$ as follows:

$\displaystyle m_{k}$ $\displaystyle=2\frac{N_{1,k}}{N_{k}}-1,$ (S240a)
$\displaystyle\sigma$ $\displaystyle=\frac{\mathcal{L}}{\frac{1}{2}\mu N},$
(S240b)

where $\mu$ is the mean degree of the network, and $N$ is the total number of
nodes, thus $\frac{1}{2}\mu N$ gives the total number of links.

###  S7.2 Dynamics of first and second moments

Eq. (S238) can be used to derive differential equations for the average degree
specific magnetisation, $\left\langle{m_{k}(t)}\right\rangle$, the average
density of active links, $\left\langle{\sigma(t)}\right\rangle$, and
correlations such as $\left\langle{m_{k}(t)m_{k^{\prime}}(t)}\right\rangle$.
Angle brackets, as always, represent averages over independent realisations of
the dynamics i.e.

$\left\langle{F[\mathbf{s}(t)]}\right\rangle=\sum_{\mathbf{s}}F[\mathbf{s}(t)]P(\mathbf{s};t),$
(S241)

where the sum is over all $2^{N}$ spin configurations
$\mathbf{s}=(s_{1},\dots,s_{N})$, and where $P(\mathbf{s};t)$ is the
probability of finding the spins in the configuration $\mathbf{s}$ at time
$t$. Specifically one finds (after some algebra),

$\displaystyle\frac{\differential\left\langle{m_{k}}\right\rangle}{\differential
t}=\left\langle{F_{k}}\right\rangle,$ (S242a)
$\displaystyle\frac{\differential\left\langle{\sigma}\right\rangle}{\differential
t}=\left\langle{F_{\sigma}}\right\rangle,$ (S242b)
$\displaystyle\frac{\differential\left\langle{m_{k}m_{k^{\prime}}}\right\rangle}{\differential
t}=\left\langle{m_{k^{\prime}}F_{k}}\right\rangle+\left\langle{m_{k}F_{k^{\prime}}}\right\rangle+\delta_{k,k^{\prime}}\frac{\left\langle{G_{k}}\right\rangle}{N},$
(S242c)

where

$\displaystyle
F_{k}=\frac{2}{N_{k}}\sum_{n=0}^{k}\left[W^{(+,k,n)}-W^{(-,k,n)}\right],$
(S243a) $\displaystyle F_{\sigma}=\frac{2}{\mu
N}\sum_{k}\sum_{n=0}^{k}(k-2n)\left[W^{(+,k,n)}-W^{(-,k,n)}\right],$ (S243b)
$\displaystyle
G_{k}=\frac{4}{P_{k}N_{k}}\sum_{n=0}^{k}\left[W^{(+,k,n)}+W^{(-,k,n)}\right],$
(S243c)

and where we have written $P_{k}=\frac{N_{k}}{N}$ for the degree distribution
of the network.

We now proceed to calculate the effective rates. $W^{(+,k,n)}$ is the product
of the number of degree $k$ nodes in the down-state, the probability of a node
in the down-state with degree $k$ having $n$ neighbours in the up-state,
$P_{0}(k,n)$, and the rate at which such a down-state node flips to an up-
state node, $R^{+}_{k,n}$. The latter is determined by the dynamics of the
model. Therefore,

$W^{(+,k,n)}=(N_{k}-N_{1,k})P_{0}(k,n)R^{+}_{k,n}.$ (S244)

In a similar manner,

$W^{(-,k,n)}=N_{1,k}P_{1}(k,n)R^{-}_{k,n},$ (S245)

where $P_{1}(k,n)$ is the probability of a node in the up-state with degree
$k$ having $n$ neighbours in the up-state, and $R^{-}_{k,n}$ is the rate at
which such a node in the up-state flips to a down-state node.

The spin flip rates in our model are different from the ones in [17]. We have

$\displaystyle R^{+}_{k,n}$
$\displaystyle=v_{+}^{\prime}+h_{+}\frac{n}{k}+h_{+}^{\prime}\frac{k-n}{k},$
(S246a) $\displaystyle R^{-}_{k,n}$
$\displaystyle=v_{-}^{\prime}+h_{-}\frac{k-n}{k}+h_{-}^{\prime}\frac{n}{k}.$
(S246b)

The probabilities $P_{0}(k,n)$ and $P_{1}(k,n)$ are intricate objects,
determined by the distribution of the full spin configuration on the network.
In order to make analytical progress we use the pair approximation and
describe these probabilities in terms of
$\left\\{\mathbf{N}_{1},\mathcal{L}\right\\}$. This is the same procedure as
in Sec. S3, specifically Eq. (S34), where we assumed that the probability of a
degree $k$ node having $n$ neighbours in the opposite state is binomial. Using
this approximation we have

$\displaystyle P_{0}(k,n)$
$\displaystyle=\binom{k}{n}c_{0}^{n}(1-c_{0})^{k-n},$ (S247a) $\displaystyle
P_{1}(k,n)$ $\displaystyle=\binom{k}{n}c_{1}^{k-n}(1-c_{1})^{n}.$ (S247b)

In these expressions $c_{0/1}$ is the probability that a given neighbour of a
node in the down/up-state is in the up/down-state (i.e. in both cases, the
probability that a link, connected to the focal node, is active). Formally,

$c_{0/1}=\frac{\mathcal{L}}{\sum_{k}kN_{0/1,k}},$ (S248)

which is the ratio between the total number of active links and the number of
links connected emanating from nodes in the up/down-state respectively.
Evaluating the above expression gives

$c_{0/1}=\frac{\sigma}{1\mp m_{L}},$ (S249)

where $m_{L}$ is the degree averaged magnetisation defined by

$m_{L}=\frac{1}{\mu}\sum_{k}kP_{k}m_{k}.$ (S250)

With these ingredients we can now evaluate the objects on the right-hand sides
of Eqs. (S243a)-(S243c) further. Due to the complexity this is done using
Mathematica. An accompanying Mathematica notebook can be found at the
following GitHub repository [47].

The resulting expressions for the objects in Eqs. (S243a)-(S243c) are long, so
we do not report them here in full. We include one expression as an example,

$F_{k}=\alpha-\beta
m_{k}+\frac{\sigma}{1-m_{L}^{2}}\Big{[}2\gamma(m_{L}-m_{k})+(\gamma_{+}-\gamma_{-})(1-m_{k}m_{L})\Big{]},$
(S251)

where we have used the shorthands from Eq. (S4).

As before we now make the restriction from Eq. (S5), $\gamma_{+}=\gamma_{-}$,
which helps maintain analytic tractability. The different quantities in Eqs.
(S243a)-(S243c) can be expressed in terms of $m_{k}$, $m_{L}$, $m$, $\mu$,
$\mu_{2}$, and the shorthand parameter combinations defined in Eq. (S4). The
quantity $m$ is the global magnetisation used throughout this paper, and can
be calculated from the degree specific magnetisation $m_{k}$,

$m=\sum_{k}P_{k}m_{k}.$ (S252)

We use the notation

$\mu_{n}=\sum_{k}P_{k}\mu^{n}$ (S253)

for the $n$th moment of the degree distribution.

We also point out a typo in the second moment of $P_{1}(k,n)$ given on p. 9 of
[17]. This should in fact be
$\sum_{n=0}^{k}n^{2}P_{1}(k,n)=k\left[k-(2k-1)c_{1}+(k-1)c_{1}^{2}\right]$.

We find that Eqs. (S242a)-(S242c) are not closed. However we can use Eq.
(S242a) to obtain a differential equation for
$\left\langle{m_{L}}\right\rangle$. To do this, we evaluate
$\frac{1}{\mu}\sum_{k}kP_{k}\left\langle{F_{k}}\right\rangle$, and find

$\frac{\differential\left\langle{m_{L}}\right\rangle}{\differential
t}=\alpha-\beta\left\langle{m_{L}}\right\rangle.$ (S254)

As a consequence we have

$\left\langle{m_{L}}\right\rangle=\left\langle{m_{L}}\right\rangle_{\rm
st}-\big{[}\left\langle{m_{L}}\right\rangle_{\rm
st}-\left\langle{m_{L}(0)}\right\rangle\big{]}e^{-\beta t},$ (S255)

and the steady-state value of $\left\langle{m_{L}}\right\rangle$ is

$\left\langle{m_{L}}\right\rangle_{\rm
st}=\frac{\alpha}{\beta}=\frac{(h_{+}^{\prime}-h_{-}^{\prime})+(v_{+}^{\prime}-v_{-}^{\prime})}{(h_{+}^{\prime}+h_{-}^{\prime})+(v_{+}^{\prime}+v_{-}^{\prime})}.$
(S256)

This is the same solution that we have found in all other cases, for instance
see Sec. S3.2 where we analysed the standard pair approximation. This is not
surprising, as all we have done in this section is perform the pair
approximation again but in a style better suited to the extension we will make
in Sec. S7.3. In fact, we have already shown in Sec. S6.3 that Eq. (S256) is a
very general solution, a fact we will use in the next section.

###  S7.3 System-size expansion

We now use what is called ‘S2PA’ (stochastic pair approximation 2) in [17] to
address Eqs. (S242a)-(S242c). This method focuses on expanding around the
deterministic solution for the variable $m_{L}$.

Motivated by the fact that $\left\langle{m_{L}}\right\rangle$ is a self
governed variable, as we saw from Eq. (S254), and that the steady-state
interface density under the standard pair approximation is as in Eq. (S43), we
purpose the following expansion:

$\displaystyle m_{k}$
$\displaystyle=m_{L}+\frac{\lambda_{k}}{\sqrt{N}}+\frac{\nu_{k}}{N}+...,$
(S257a) $\displaystyle\sigma$
$\displaystyle=\frac{1}{2}H_{\text{PA}}(\tau)(1-m_{L}^{2})+\frac{\lambda_{\sigma}}{\sqrt{N}}+...,$
(S257b)

where $\lambda_{k}$, $\nu_{k}$ and $\lambda_{\sigma}$ are stochastic
variables. All objects in these equations, excluding $H_{\rm PA}(\tau)$ which
is given in Eq. (S45), are understood to be time-dependent. Ultimately
however, we are only interested in the stationary state values. We neglect
higher-order terms in the expansion in powers of $N^{-1/2}$. The goal is to
calculate $\left\langle{\sigma}\right\rangle_{\rm st}$ to leading and sub-
leading order. Any expansions we perform in the remainder of this section can
be found in a Mathematica notebook at the following GitHub repository [47].

This expansion is further justified by the fact that individual simulated
trajectories in $(m,\sigma)$-space, performed on finite size uncorrelated
networks, fluctuate around a parabola with height $H_{\rm PA}(\tau)$. This is
demonstrated in Fig. 2 in the main paper. Simulated trajectories on a truly
infinite uncorrelated network would move exactly along this parabola, in fact
they would follow the path described by Eqs. (S25) and (S42) which is also
shown in Fig. 2. Time averaging such analytical trajectories on infinite
networks in the steady-state (equivalently we can average over many such
independent trajectories at a single time in the steady state) will result in
the point given by Eqs. (S256) and (S43). However, stochasticity means that
the steady-state averages in finite networks sit on a different parabola. This
is what we wish to determine and what motivates the system-size expansion in
Eqs. (S257a) and (S257b).

To continue, we note that based on the results from Sec. S6.3, that any
normalised weighted magnetisation in the steady-state is equal to the steady-
state global magnetisation, Eq. (S164). The degree-specific average
magnetisation $\left\langle{m_{k}}\right\rangle_{\rm st}$ is effectively a
weighted magnetisation with weights $\delta_{k,k_{i}}$, i.e. we have

$\left\langle{m_{k}}\right\rangle=\frac{\sum_{i}\delta_{k,k_{i}}S_{i}}{\sum_{i}\delta_{k,k_{i}}},$
(S258)

where $S_{i}$ is the average spin state at site $i$. Thus we have
$\left\langle{m_{k}}\right\rangle_{\rm st}=\left\langle{m}\right\rangle_{\rm
st}$. Given the definition of $m_{L}$ in Eq. (S250) this then also means that
$\left\langle{m_{k}}\right\rangle_{\rm
st}=\left\langle{m_{L}}\right\rangle_{\rm st}$ for all $k$. With this we see
immediately from Eq. (S257a) that $\left\langle{\lambda_{k}}\right\rangle_{\rm
st}=0$ and $\left\langle{\nu_{k}}\right\rangle_{\rm st}=0$.

Next we substitute Eqs. (S257a) and (S257b) into Eqs. (S242b) and (S242c), and
equate powers of $N^{-1/2}$. Many model parameters that appear, i.e. $\alpha$,
$\beta$, $\gamma$ etc. from Eq. (S4), can be expressed in terms of the noise-
to-order ratio $\tau=\frac{\beta}{\gamma}$ from Eq. (S46) and the steady-state
global magnetisation $\left\langle{m}\right\rangle_{\rm
st}=\frac{\alpha}{\beta}$ from Eq. (S24). We find

$\displaystyle\frac{\differential\left\langle{\lambda_{\sigma}}\right\rangle}{\differential
t}$ $\displaystyle=-2\left\\{\tau-1+2H_{\rm
PA}(\tau)+\frac{2}{\mu}\big{[}1-H_{\rm
PA}(\tau)\big{]}\right\\}\left\langle{\lambda_{\sigma}}\right\rangle$
$\displaystyle\qquad+\frac{2\gamma}{\mu}H_{\rm PA}(\tau)\big{[}1-H_{\rm
PA}(\tau)\big{]}\sum_{k}P_{k}\left\langle{\lambda_{k}m_{L}}\right\rangle$
(S259)

and

$\frac{\differential\left\langle{\lambda_{k}m_{L}}\right\rangle}{\differential
t}=-\gamma\big{[}H_{\text{PA}}(\tau)+2\tau\big{]}\left\langle{\lambda_{k}m_{L}}\right\rangle+\gamma\tau\left\langle{m}\right\rangle_{\rm
st}\left\langle{\lambda_{k}}\right\rangle$ (S260)

respectively. We note that to get Eq. (S260) we set $k=k^{\prime}$ in Eq.
(S242c).

Taking both of these equations to the steady-state, and given that we already
know that $\left\langle{\lambda_{k}}\right\rangle_{\rm st}=0$, we also find
that $\left\langle{\lambda_{\sigma}}\right\rangle_{\rm st}=0$. With this
$\left\langle{\sigma}\right\rangle_{\rm st}$ reduces to

$\left\langle{\sigma}\right\rangle_{\rm st}=\frac{1}{2}H_{\rm
PA}(\tau)\Big{(}1-\left\langle{m_{L}^{2}}\right\rangle_{\rm st}\Big{)},$
(S261)

and we are left with determining $\left\langle{m_{L}^{2}}\right\rangle_{\rm
st}$. Comparing with the result from the conventional pair approximation in
Eq. (S43), we note that the final result contains the object
$\left\langle{m_{L}^{2}}\right\rangle_{\rm st}$ instead of
$\left\langle{m}\right\rangle_{\rm st}^{2}$. In particular the square is now
inside the average.

To calculate $\left\langle{m_{L}^{2}}\right\rangle_{\rm st}$, we start from
Eq. (S242c),

$\frac{\differential\left\langle{m_{k}m_{k^{\prime}}}\right\rangle}{\differential
t}=\left\langle{m_{k^{\prime}}F_{k}}\right\rangle+\left\langle{m_{k}F_{k^{\prime}}}\right\rangle+\delta_{k,k^{\prime}}\frac{\left\langle{G_{k}}\right\rangle}{N}.$
(S262)

We expand this equation in powers of $\frac{1}{\sqrt{N}}$ using Eqs. (S257a).
The LHS of Eq. (S262) then becomes

$\displaystyle\frac{\differential\left\langle{m_{k}m_{k^{\prime}}}\right\rangle}{\differential
t}$
$\displaystyle=\frac{\differential\left\langle{m_{L}^{2}}\right\rangle}{\differential
t}+\frac{1}{\sqrt{N}}\left(\frac{\differential\left\langle{\lambda_{k}m_{L}}\right\rangle}{\differential
t}+\frac{\differential\left\langle{\lambda_{k^{\prime}}m_{L}}\right\rangle}{\differential
t}\right)$
$\displaystyle\qquad+\frac{1}{N}\left(\frac{\differential\left\langle{\lambda_{k}\lambda_{k^{\prime}}}\right\rangle}{\differential
t}+\frac{\left\langle{\nu_{k}m_{L}}\right\rangle}{\differential
t}+\frac{\left\langle{\nu_{k^{\prime}}m_{L}}\right\rangle}{\differential
t}\right)+{\cal O}(N^{-3/2}).$ (S263)

Next we consider the term $\left\langle{m_{k^{\prime}}F_{k}}\right\rangle$ on
the RHS of Eq. (S262). This expands as [47]

$\displaystyle\left\langle{m_{k^{\prime}}F_{k}}\right\rangle$
$\displaystyle=\gamma\tau\left[\left\langle{m}\right\rangle_{\rm
st}\left\langle{m_{L}}\right\rangle-\left\langle{m_{L}^{2}}\right\rangle\right]$
$\displaystyle\qquad+\frac{\gamma}{\sqrt{N}}\Big{[}-\Big{(}H_{\rm
PA}(\tau)+\tau\Big{)}\left\langle{\lambda_{k}m_{L}}\right\rangle-\tau\left\langle{\lambda_{k^{\prime}}m_{L}}\right\rangle+\tau\left\langle{\lambda_{k^{\prime}}}\right\rangle\left\langle{m}\right\rangle_{\rm
st}\Big{]}$ $\displaystyle\qquad+\frac{\gamma}{N}\Bigg{[}-\Big{(}H_{\rm
PA}(\tau)+\tau\Big{)}\left\langle{\lambda_{k}\lambda_{k^{\prime}}}\right\rangle+\tau\left\langle{\nu_{k}}\right\rangle\left\langle{m}\right\rangle_{\rm
st}-2\left\langle{\frac{m_{L}\lambda_{\sigma}\lambda_{k^{\prime}}}{1-m_{L}^{2}}}\right\rangle$
$\displaystyle\qquad\qquad-\tau\left\langle{\nu_{k}m_{L}}\right\rangle-\Big{(}H_{\rm
PA}(\tau)+\tau\Big{)}\left\langle{m_{L}\nu_{k^{\prime}}}\right\rangle\Bigg{]}$
$\displaystyle\qquad+{\cal O}(N^{-3/2}).$ (S264)

The $\left\langle{m_{k}F_{k^{\prime}}}\right\rangle$ term gives a similar
expression, but with $k$ and $k^{\prime}$ interchanged.

Since $\left\langle{G_{k}}\right\rangle$ carries a pre-factor $1/N$ in Eq.
(S262) we only consider the leading order term of
$\left\langle{G_{k}}\right\rangle$. We have [47]

$\left\langle{G_{k}}\right\rangle=\frac{2\gamma}{P_{k}}\Big{[}H_{\rm
PA}(\tau)+\tau-H_{\rm
PA}(\tau)\left\langle{m_{L}^{2}}\right\rangle-\tau\left\langle{m}\right\rangle_{\rm
st}\left\langle{m_{L}}\right\rangle\Big{]}+{\cal O}(N^{-1/2}).$ (S265)

Collecting terms in the expansions of the left-hand and right-hand sides of
Eq. (S262), and keeping only expressions up to and including order $1/N$, we
have

$\displaystyle\frac{\differential\left\langle{m_{L}^{2}}\right\rangle}{\differential
t}$
$\displaystyle+\frac{1}{\sqrt{N}}\left(\frac{\differential\left\langle{\lambda_{k}m_{L}}\right\rangle}{\differential
t}+\frac{\differential\left\langle{\lambda_{k^{\prime}}m_{L}}\right\rangle}{\differential
t}\right)+\frac{1}{N}\left(\frac{\differential\left\langle{\lambda_{k}\lambda_{k^{\prime}}}\right\rangle}{\differential
t}+\frac{\left\langle{\nu_{k}m_{L}}\right\rangle}{\differential
t}+\frac{\left\langle{\nu_{k^{\prime}}m_{L}}\right\rangle}{\differential
t}\right)=$ $\displaystyle\qquad
2\gamma\tau\left[\left\langle{m}\right\rangle_{\rm
st}\left\langle{m_{L}}\right\rangle-\left\langle{m_{L}^{2}}\right\rangle\right]$
$\displaystyle\qquad+\frac{\gamma}{\sqrt{N}}\Bigg{\\{}-\Big{(}H_{\rm
PA}(\tau)+2\tau\Big{)}\Big{[}\left\langle{\lambda_{k}m_{L}}\right\rangle+\left\langle{\lambda_{k^{\prime}}m_{L}}\right\rangle\Big{]}+\tau\left\langle{m}\right\rangle_{\rm
st}\Big{[}\left\langle{\lambda_{k}}\right\rangle+\left\langle{\lambda_{k^{\prime}}}\right\rangle\Big{]}\Bigg{\\}}$
$\displaystyle\qquad+\frac{\gamma}{N}\Bigg{\\{}-2\Big{(}H_{\rm
PA}(\tau)+\tau\Big{)}\left\langle{\lambda_{k}\lambda_{k^{\prime}}}\right\rangle+\tau\left\langle{m}\right\rangle_{\rm
st}\Big{[}\left\langle{\nu_{k}}\right\rangle+\left\langle{\nu_{k^{\prime}}}\right\rangle\Big{]}-2\left\langle{\frac{(\lambda_{k}+\lambda_{k^{\prime}})m_{L}\lambda_{\sigma}}{1-m_{L}^{2}}}\right\rangle$
$\displaystyle\qquad\qquad\qquad-\Big{(}H_{\rm
PA}(\tau)+2\tau\Big{)}\Big{[}\left\langle{\nu_{k}m_{L}}\right\rangle+\left\langle{\nu_{k^{\prime}}m_{L}}\right\rangle\Big{]}$
$\displaystyle\qquad\qquad\qquad+\frac{2}{P_{k}}\Big{[}H_{\rm PA}(\tau)+\tau-
H_{\rm
PA}(\tau)\left\langle{m_{L}^{2}}\right\rangle-\tau\left\langle{m}\right\rangle_{\rm
st}\left\langle{m_{L}}\right\rangle\Big{]}\delta_{k,k^{\prime}}\Bigg{\\}}.$
(S266)

Performing the operation $\sum_{kk^{\prime}}kk^{\prime}P_{k}P_{k^{\prime}}$ on
both sides, several of the contributions can be eliminated. To see this, we
perform the sum $\sum_{k}kP_{k}$ on the expansion of $m_{k}$ from Eq. (S257a),
and find

$\sum_{k}kP_{k}m_{k}=\sum_{k}kP_{k}m_{L}+\frac{1}{\sqrt{N}}\sum_{k}kP_{k}\lambda_{k}+\frac{1}{N}\sum_{k}kP_{k}\nu_{k},$
(S267)

resulting in

$\mu m_{L}=\mu
m_{L}+\frac{1}{\sqrt{N}}\sum_{k}kP_{k}\lambda_{k}+\frac{1}{N}\sum_{k}kP_{k}\nu_{k}.$
(S268)

Equating powers of $N^{-1/2}$ we therefore conclude

$\displaystyle\sum_{k}kP_{k}\lambda_{k}$ $\displaystyle=0,$ (S269a)
$\displaystyle\sum_{k}kP_{k}\nu_{k}$ $\displaystyle=0.$ (S269b)

We now use the identities in Eqs. (S269a) and (S269b) when performing the
operation
$\sum_{kk^{\prime}}kk^{\prime}P_{k}P_{k^{\prime}}$ in Eq. (S266). This results
in all stochastic terms vanishing and we are left with

$\displaystyle\frac{\differential\left\langle{m_{L}^{2}}\right\rangle}{\differential
t}$
$\displaystyle=-2\gamma\tau\Big{[}\left\langle{m_{L}^{2}}\right\rangle-\left\langle{m}\right\rangle_{\rm
st}\left\langle{m_{L}}\right\rangle\Big{]}$
$\displaystyle\qquad+\frac{2\gamma\mu_{2}}{N\mu^{2}}\Big{[}H_{\rm
PA}(\tau)+\tau-H_{\rm
PA}(\tau)\left\langle{m_{L}^{2}}\right\rangle-\tau\left\langle{m}\right\rangle_{\rm
st}\left\langle{m_{L}}\right\rangle\Big{]},$ (S270)

where we note the introduction of $\mu_{2}$ and $\mu^{2}$ arising from terms
of the form $\sum_{k}k^{2}P_{k}$ and $(\sum_{k}kP_{k})^{2}$. It is possible to
solve this equation analytically and obtain a full time dependent solution for
$\left\langle{m_{L}^{2}(t)}\right\rangle$ and thus
$\left\langle{\sigma(t)}\right\rangle$ using Eq. (S255).

However, we are mostly interested in the steady-state, so we move to that
regime. Eq. (S270) then becomes

$0=-\tau N\Big{[}\left\langle{m_{L}^{2}}\right\rangle_{\rm
st}-\left\langle{m}\right\rangle^{2}_{\rm
st}\Big{]}+\frac{\mu_{2}}{\mu^{2}}\Big{[}H_{\rm PA}(\tau)+\tau-H_{\rm
PA}(\tau)\left\langle{m_{L}^{2}}\right\rangle_{\rm
st}-\tau\left\langle{m}\right\rangle^{2}_{\rm st}\Big{]}$ (S271)

and we can solve for $\left\langle{m_{L}^{2}}\right\rangle_{\rm st}$:

$\left\langle{m^{2}_{L}}\right\rangle_{\rm st}=\frac{(H_{\rm
PA}(\tau)+\tau)\mu_{2}-\tau(\mu_{2}-N\mu^{2})\left\langle{m}\right\rangle_{\rm
st}^{2}}{N\mu^{2}\tau+H_{\rm PA}(\tau)\mu_{2}}.$ (S272)

Finally, we substitute this back into Eq. (S261) and obtain

$\left\langle{\sigma}\right\rangle_{\rm
st}=\frac{1}{2}H_{\text{SPA}}(\tau)\Big{(}1-\left\langle{m}\right\rangle_{\rm
st}^{2}\Big{)},$ (S273)

where

$H_{\text{SPA}}(\tau)=\frac{\tau
H_{\text{PA}}(\tau)(N\mu^{2}-\mu_{2})}{N\mu^{2}\tau+H_{\text{PA}}(\tau)\mu_{2}},$
(S274)

(the subscript ‘SPA’ of course stands for ‘stochastic pair approximation’).
Expanding in $1/N$ Eq. (S274) can be written as

$H_{\text{SPA}}(\tau)=H_{\text{PA}}(\tau)\left[1-\frac{\mu_{2}}{N\mu^{2}}\frac{\tau+H_{\text{PA}}(\tau)}{\tau}\right]+{\cal
O}(N^{-2}),$ (S275)

making it transparent that the result from the stochastic pair approximation
is a finite-size correction to the expression obtained from the conventional
pair approximation.

Due to the terms involving $H_{\rm PA}(\tau)$ this solution applies in the
same regions as that of the conventional pair approximation, see Sec. S3.3 for
a detailed discussion of this, as well as Secs. S8 and S10.2 for discussions
on the $\tau$ limits in general.

Focusing on positive $\tau$, we substitute the solution for $H_{\rm PA}(\tau)$
from Eq. (S45) into Eq. (S274) and then expand in powers of $\tau$. To lowest
order we find

$H_{\rm
SPA}(\tau)=\left(\frac{N\mu^{2}}{\mu_{2}}-1\right)\tau+\order{\tau^{2}}.$
(S276)

Thus it is clear that in the limit $\tau\to 0$, $H_{\rm SPA}(\tau)\to 0$. This
matches what we see in simulations on finite size networks.

In the limit $\tau\to\infty$, Eq. (S275) becomes

$\lim_{\tau\to\infty}H_{\rm SPA}(\tau)=1-\frac{\mu_{2}}{N\mu^{2}}.$ (S277)

This is of similar form to Eq. (S111) when we studied the annealed
approximation. We know the large $\tau$ regime is dominated by random state
changes so there are no correlations between neighbours, hence we expect
$H(\tau)\to 1$, as this corresponds to random scattering. Eq. (S277)
reproduces this, up to corrections in $1/N$.

In Fig. S12 we demonstrate how Eq. (S274) compares to simulation for different
size degree-regular (left), Erdös-–Rényi (middle), and Barabási–Albert (right)
networks. We find convincing agreement across all ranges of $\tau$. We note
that the approximation becomes inaccurate for very small networks [see Fig.
3(b) in the main paper].

Figure S12: Plots of $H(\tau)=\left\langle{\sigma}\right\rangle_{\rm
st}/[\frac{1}{2}(1-\left\langle{m}\right\rangle_{\rm st}^{2})]$, which
characterises the scattering, as a function of the noise-to-order ratio,
$\tau$, for networks of size $N$. Degree-regular with $\mu=8$ (left),
Erdös–Rényi with $p=8/N$ (middle), Barabási–Albert with $4$ new links
established upon adding a new node (right). For all three types of network we
have $\mu\approx 8$ (deviations are due to the finite size of the graphs). The
dashed lines are the results from the SPA approximation, Eq. (S274). Each
marker is the result from averaging 1000 independent simulations in the
steady-state at a particular value of $\tau$. Model parameters for a
particular $\tau$ value can be generated via the algorithm detailed in Sec.
S10.2.

##  S8 Scattering and the noise-to-order ratio

###  S8.1 Relation between $\mathbb{H}$ and $H(\tau)$

####  S8.1.1 Degree of scatter for individual realisations

For individual configurations of the system, the parameter $\mathbb{H}$ is
defined as [Eq.(2)]

$\mathbb{H}\equiv\frac{\sigma}{\frac{1}{2}(1-m^{2})}.$ (S278)

$\mathbb{H}$ quantifies the ‘scatter’ of a spin configuration, relative to a
random configuration with the same magnetisation (same fraction of up-spins).
More specifically, the denominator in Eq. (S278) is the (expected) fraction of
links that are active if spins of total magnetisation $m$ are placed randomly
on the network. The quantity $\mathbb{H}$ is not defined if $m=\pm 1$, i.e.
all spins are in the same state. From an ensemble of independent simulation
runs we can measure the quantity $\left\langle{\mathbb{H}(t)}\right\rangle$.
In computing this average, we discard any configurations where $\mathbb{H}$ is
not defined. Such configurations are very rare provided $\tau$ is not too
small. The exclusion is always implied when we write
$\left\langle{\mathbb{H}}\right\rangle$.

In Fig. S13 we show $\mathbb{H}(t)$ (green line) for a single simulation on a
Barabási-Albert network of size $N=2500$, and
$\left\langle{\mathbb{H}(t)}\right\rangle$ obtained from many independent
simulations (blue line). The value of $\mathbb{H}(t)$ fluctuates around
$\left\langle{\mathbb{H}}\right\rangle_{\rm st}$. The amplitude of the
fluctuations is about one tenth of the mean in this example. Therefore, the
value of $\mathbb{H}$ measured from a single configuration in the stationary
state will be well approximated by $\left\langle{\mathbb{H}}\right\rangle_{\rm
st}$.

Figure S13: Time evolution of various order parameters for a Barabási-Albert
network with $\mu\approx 8$ and $N=2500$. The red and purple dashed lines are
theory lines from $H_{\rm PA}(\tau)$ [Eq. (S45)] and $H_{\rm SPA}(\tau)$ [Eq.
(S274)] respectively. The solid light green line is the measured value of
$\mathbb{H}(t)$ in a single simulation. The solid blue line is
$\left\langle{\mathbb{H}(t)}\right\rangle$ over many such simulations. The
solid orange line is $H(\tau,t)$, as defined in the text and in the legend,
and obtained as an average over many simulations. Data in this example is for
$\tau=10^{-2}$.

####  S8.1.2 Relation between average density of active interfaces and average
magnetisation

We define a second order parameter

$H(\tau)\equiv\frac{\left\langle{\sigma}\right\rangle_{\rm
st}}{\frac{1}{2}(1-\left\langle{m}\right\rangle_{\rm st}^{2})}.$ (S279)

Even if individual configurations are in ‘consensus’ states ($m=\pm 1$), time
averages will lead to $\left\langle{m}\right\rangle_{\rm st}\neq\pm 1$ (except
for pathological choices of the model parameters). We can also define
$H(\tau,t)\equiv\frac{\left\langle{\sigma(t)}\right\rangle}{\frac{1}{2}(1-\left\langle{m(t)}\right\rangle^{2})}$.
This quantity is shown as an orange line in Fig. S13 (for $\tau=10^{-2}$). The
data demonstrates that $H(\tau)$ will in general not take the same value as
$\left\langle{\mathbb{H}}\right\rangle_{\rm st}$ (blue line).

We note in Fig. S13 that, in the stationary state, $\mathbb{H}$ for a single
realisation (green line) fluctuates around $H_{\rm PA}(\tau)$ while
$H(\tau,t)$ (orange line) approaches $H_{\rm SPA}(\tau)$. This is analogous to
Fig. 2 in the main paper, where individual trajectories fluctuate around a
parabola whose height is set by $H_{\rm PA}(\tau)$ but where the points
defined by $\left\langle{\sigma}\right\rangle_{\rm st}$ and
$\left\langle{m}\right\rangle_{\rm st}$ sit on a parabola with a different
height $H_{\rm SPA}(\tau)$.

####  S8.1.3 $H(\tau)$ can be inferred from measurement $\mathbb{H}$ in
individual configurations, for large $\tau$

In Fig. S14 we compare $\left\langle{\mathbb{H}}\right\rangle_{\rm
st}=\left\langle{\frac{\sigma}{\frac{1}{2}(1-m^{2})}}\right\rangle_{\rm st}$
and $H(\tau)=\frac{\left\langle{\sigma}\right\rangle_{\rm
st}}{\frac{1}{2}(1-\left\langle{m}\right\rangle_{\rm st}^{2})}$ over a range
of values of $\tau$ for square lattices and Barabási–Albert networks of
different sizes. We find that $H(\tau)$ measured in simulations follows the
analytical predictions $H_{\rm HG}(\tau)$ in the case of finite lattices, and
$H_{\rm SPA}(\tau)$ for the Barabási–Albert network.

We note in Fig. S14 that $\left\langle{\mathbb{H}}\right\rangle_{\rm st}$
measured in simulations agrees with the analytical predictions for $H(\tau)$
(and with the measurements of this quantity in simulation) above some
$\tau_{C}$. The value of $\tau_{C}$ depends on the type and size of the
network. Deviations are seen for $\tau<\tau_{C}$. The value of $\tau_{C}$
decreases as $N$ becomes larger, see e.g. the triangle markers in Fig. S14(a).

Figure S14: Plots comparing $\left\langle{\mathbb{H}}\right\rangle_{\rm st}$
(triangles) and $H(\tau)$ (circles), from Eqs. (S278) and (S279), in
simulations on square lattices with periodic boundary conditions (left) and
Barabási-Albert networks with $\mu\approx 8$ (right), of varying size $N$.
Left panel: dotted line is the $H_{\rm SL}(\tau)$ theory line for an infinite
square lattice [Eq. (S125)], dashed lines are the $H_{\rm HG}(\tau)$ theory
lines for homogeneous networks [Eq. (S188)]. Right panel: dotted lines are
$H_{\rm PA}(\tau)$ theory lines for infinite uncorrelated networks [Eq.
(S45)], dashed lines are the $H_{\rm SPA}(\tau)$ theory lines for large but
finite uncorrelated networks [Eq. (S274)]. In both panels we find that there
is much more stochasticity in $\left\langle{\mathbb{H}}\right\rangle_{\rm st}$
(triangles) than in $H(\tau)$ (circles). This is because if the system reaches
consensus (all spins up or down), $\mathbb{H}_{\rm st}$ cannot be measured and
the sample is discarded. This happens the majority of the time at low $\tau$
values meaning many more simulations are needed to measure
$\left\langle{\mathbb{H}}\right\rangle_{\rm st}$ than $H(\tau)$

Analytical functions for the order parameter $H(\tau)$ can be derived for many
networks. These are the functions $H_{X}(\tau)$ [Eqs. (S45), (S110), (S188)
and (S274)]. These functions depend only on properties of the network (e.g.
the average degree $\mu$), and their argument is $\tau$ [see Sec. S8.2]. This,
together with the above, provides a method to infer $\tau$ from snap shots of
spin configurations in the steady-state.

####  S8.1.4 Range of possible $H(\tau)$ values

Since $m$ is bound on $[-1,1]$, and $\sigma$ on $[0,1]$, $H(\tau)$ can, in
principle, take any real value greater than or equal to 0. However, not all
combinations of $\left\langle{m}\right\rangle_{\rm st}$ and
$\left\langle{\sigma}\right\rangle_{\rm st}$ are realised by the model
dynamics. In practice the maximum value of $H(\tau)$ that can be achieved
depends on the topology of the network. For general networks, finding the
assignment of states $\pm 1$ for each node leading to the largest possible
number of active links is closely related to the maximum-cut problem, which is
known to NP-hard [48].

###  S8.2 Noise-to-order ratio, $\tau$

####  S8.2.1 Definition and interpretation

The noise-to-order ratio, $\tau$, is given by [see Eq. (S46)]

$\tau=2\frac{h_{+}^{\prime}+h_{-}^{\prime}+v_{+}^{\prime}+v_{-}^{\prime}}{(h_{+}+h_{-})-(h_{+}^{\prime}+h_{-}^{\prime})}.$
(S280)

Interpretation of the numerator
To interpret the numerator we imagine a fully ordered spin configuration, that
is, all spins up or all spins down. For example, if all spins are up, then the
only possible events are spontaneous down-flips or down-flips due to incorrect
copying. As seen from Eqs. (S2) and (S3) the total event rate for the former
is $v_{-}^{\prime}N$ and for the latter $h_{-}^{\prime}N$. So the total event
rate, starting from the all-up-state is $(v_{-}^{\prime}+h_{-}^{\prime})N$.
Similarly, the total event rate starting from the all down configuration is
$(v_{+}^{\prime}+h_{+}^{\prime})N$. Up to the factor of $N$, we can thus think
of the numerator in Eq. (S280) as the average event rate in a fully ordered
state. Any event in such a state moves the system towards disorder.

Interpretation of the denominator
We now show that the denominator describes the initial rate of the build-up of
nearest-neighbour correlations if the system is started from a fully
uncorrelated state with $m=0$. In other words, we assume that, in the initial
state, each spin is up or down with equal probability, and independently of
its neighbours. Then we calculate $\frac{\differential}{\differential
t}\left[\left\langle{ss^{\prime}}\right\rangle-\left\langle{s}\right\rangle\left\langle{s^{\prime}}\right\rangle\right]$.

If a give node is the state $+1$ and if $n$ of its $k$ nearest neighbours are
in state $-1$, then the rate at which it flips to $-1$ is

$\displaystyle T_{n,k,t}^{-}$
$\displaystyle=h_{-}\frac{n}{k}+h_{-}^{\prime}\frac{k-n}{k}+v_{-}^{\prime}.$
(S281)

Spontaneous down flips in one node (described by $v_{-}^{\prime}$) cannot
affect correlations between neighbours. In fact, in the fully random state
with $m=0$ such a spontaneous flip will, on average, convert as many active
interfaces into inactive interfaces as the other way round. Therefore, we do
not consider this term any further for the purposes of calculating the build-
up of correlations.

The first two terms in Eq. (S281) describe the following: A random neighbour
of the node is chosen. If that neighbour is in the down-state (this happens
with probability $1/2$ in the fully uncorrelated state) the focal node flips
with rate $h_{-}$. The pair of focal node and neighbour change from $+-$ to
$--$. The change in $ss^{\prime}$ is thus $+2$. If instead the neighbour is in
the up-state (which also happens with probability $1/2$), the focal node flips
with rate $h_{-}^{\prime}$. The pair of focal node and neighbour change from
$++$ to $-+$, and the change in $ss^{\prime}$ is $-2$.

We note that interfaces of the focal node with any other nearest neighbours
(other than the designated interaction partner) are equally likely to change
from active to inactive or from inactive to active, given that we are starting
from the fully uncorrelated state. These interfaces therefore do not
contribute to any build-up of correlations.

The total expected rate of change of $\left\langle{ss^{\prime}}\right\rangle$
due to events described by Eq. (S281) is thus
$\frac{1}{2}h_{-}\cdot(+2)+\frac{1}{2}h_{-}^{\prime}\cdot(-2)=h_{-}-h_{-}^{\prime}$.

Similarly, if the focal node is in state $-1$ and $n$ is the number of
neighbours in the up-state, then the rate at which it flips to state $+1$ is

$\displaystyle T_{n,k,t}^{+}$
$\displaystyle=h_{+}\frac{n}{k}+h_{+}^{\prime}\frac{k-n}{k}+v_{+}^{\prime}.$
(S282)

We can make an analogous argument as above to show that the total expected
rate of change of $\left\langle{ss^{\prime}}\right\rangle$ due to events
described by Eq. (S282) is thus $h_{+}-h_{+}^{\prime}$.

Finally, we note that

$\frac{\differential}{\differential
t}\big{(}\left\langle{s}\right\rangle\left\langle{s^{\prime}}\right\rangle\big{)}=\left\langle{s}\right\rangle\frac{\differential}{\differential
t}\left\langle{s^{\prime}}\right\rangle+\left\langle{s^{\prime}}\right\rangle\frac{\differential}{\differential
t}\left\langle{s}\right\rangle=0$ (S283)

in the initial state in which $m=0$, i.e. $\left\langle{s}\right\rangle=0$.

Putting everything together we therefore find that

$\frac{\differential}{\differential
t}\left\langle{ss^{\prime}}\right\rangle\propto
h_{-}-h_{-}^{\prime}+h_{+}-h_{+}^{\prime}.$ (S284)

Thus, we have shown that the denominator of the expression in Eq. (S280)
describes (up to constant pre-factors) the initial rate of the build-up of
nearest neighbour correlations if the dynamics is started from a fully
disordered state in which each spin is up or down with probability $1/2$, and
independently of its neighbours.

####  S8.2.2 Range of possible values for $\tau$

The definition in Eq. (S280) is of the form

$\tau=2\frac{a+b}{c-a},$ (S285)

with $a,b,c,d\geq 0$.

It is clear from this definition that all values $\tau\geq 0$ can be realised
through a suitable choice of $a,b,c,d$.

The quantity $\tau$ takes negative values, if $a>c$. In this case we have

$-\tau^{-1}=\frac{1}{2}\frac{a-c}{a+b},$ (S286)

and $0<a-c\leq a\leq a+b$. Thus -$\tau^{-1}\leq 1/2$, and hence $\tau\leq-2$.

We conclude that $\tau$ can only take values in the ranges $\tau\geq 0$ or
$\tau\leq-2$.

###  S8.3 Behaviour in different regimes

In Fig. S15 we show snapshots of the spin configuration in the stationary
state, obtained from simulations of the model on a $4\times 4$ square grid.

Figure S15: Snapshots of the spin distribution of a single realisation on a
grid network after running the model dynamics for a long time, for different
values of $\tau$. Nodes in state $+1$ are blue. Nodes in state $-1$ are
yellow. Links connecting different state nodes are active (solid). Links
connecting same state nodes are inactive (dashed).

####  S8.3.1 The limit $\tau\to-\infty$

The limit $\tau\to-\infty$ emerges in two different ways:

(i) One or both of $v_{\pm}^{\prime}$ tend to infinity, while
$h_{+}^{\prime}+h_{-}^{\prime}>h_{+}+h_{-}$; or

(ii) The denominator of Eq. (S280) approaches zero from below.

Typical configurations of the system will be mostly random, but there will be
some small anti-correlation between nodes so the number of active links is
slightly greater than random. We find $H(\tau)\to 1^{-}$ in this case. This
case is illustrated by the configuration in Fig. S15(a), for which we find
$\mathbb{H}=\sigma/[\frac{1}{2}(1-m^{2})]=1.02$.

####  S8.3.2 Finite negative values of $\tau<-2$

The copying process is mostly unfaithful, but there can also be faithful
copying and/or spontaneous spin flips. There will be some anti-correlation
between neighbouring spins, and thus, for a given trait frequency $m$, one
expects more active interfaces than a configuration in which these spins are
placed randomly on the network. Hence, we expect values of $H(\tau)>1$.

A case in which $\tau$ is negative and finite is illustrated by the
configuration in Fig. S15(b), for which
$\mathbb{H}=\sigma/[\frac{1}{2}(1-m^{2})]=1.27$.

####  S8.3.3 The limit $\tau\to-2^{-}$

From the definition of $\tau$ in Eq. (S280) we have

$\tau=-2\frac{h_{+}^{\prime}+h_{-}^{\prime}+v_{+}^{\prime}+v_{-}^{\prime}}{(h_{+}^{\prime}+h_{-}^{\prime})-(h_{+}+h_{-})}.$
(S287)

Hence $\tau=-2$ when either $h_{+}^{\prime}+h_{-}^{\prime}$ diverges and the
$v_{\pm}^{\prime},h_{\pm}$ remain finite, or if the $v^{\prime}_{\pm}$ and
$h_{\pm}$ all vanish. In either case, the dynamics is dominated by unfaithful
copying from neighbours, i.e. the model is type of ‘contrarian voter model’,
in which a spin copies the state opposite to that of a neighbour. This creates
strong anti-correlations between the states of neighbouring spins, and
$H(\tau)$ will take its maximum value which is specific to the network.

This scenario is illustrated in Fig. S15(c), where we show a perfect
‘checkerboard’ pattern seen in simulations on a small square grid. For this
pattern we have $\mathbb{H}=\sigma/[\frac{1}{2}(1-m^{2})]=2$.

####  S8.3.4 The limit $\tau\to 0^{+}$

The limit $\tau\to 0^{+}$ can occur in one of two ways:

(i) The numerator in Eq. (S280) goes to zero (from above) and, at the same
time, $h_{+}+h_{-}>h_{+}^{\prime}+h_{-}^{\prime}$, so that denominator is
positive; or

(ii) One has $h_{+}+h_{-}\to\infty$.

In the former case we have $h_{\pm}^{\prime}\to 0$ and $v_{\pm}^{\prime}\to
0$. There is no transmission noise. In the latter case faithful copying from
neighbouring spins dominates the dynamics.

In either case the number of active links is low and we expect configurations
in which either all spins are up or all spins are down, as illustrated in Fig.
S15(d). Due to the convergence to a fully ordered state, the ratio
$\mathbb{H}=\sigma/[\frac{1}{2}(1-m^{2})]$ cannot be defined for single
realisations. However all realisations have $\sigma=0$ eventually, and we
measure $H(\tau)=\left\langle{\sigma}\right\rangle_{\rm
st}/[\frac{1}{2}(1-\left\langle{m}\right\rangle_{\rm st}^{2})]=0$.

####  S8.3.5 Finite $\tau>0$

For $\tau>0$ and finite, there are some unfaithful events occurring, and the
rate of faithful horizontal events is greater than the total rate of
unfaithful horizontal events, i.e.
$h_{+}+h_{-}>h_{+}^{\prime}+h_{-}^{\prime}$. This leads to steady-state
situations where $0<H(\tau)<1$. The case is illustrated in Fig. S15(e), where
we measure $\mathbb{H}=\sigma/[\frac{1}{2}(1-m)^{2}]=0.33$.

####  S8.3.6 The limit $\tau\to\infty$

This limit emerges in two different ways:

(i) The denominator in Eq. (S280) remains positive, and either
$v_{+}^{\prime}$ or $v_{-}^{\prime}$ (or both) tend to infinity; or

(ii) The $v_{\pm}^{\prime}$ remain finite, but the denominator in Eq. (S280)
diverges, i.e., $h_{+}+h_{-}=h_{+}^{\prime}+h_{-}^{\prime}$.

In the former case, the system is dominated by random spin flips. In the
second scenario, correct and incorrect copying from neighbouring nodes occur
at the same rates. This means that the copying spin takes a state that is
effectively unrelated to that of the state it is trying to copy.

In neither of these two situations correlations between neighbouring states
can build up, and so the resulting configurations are effectively random. In
general $H(\tau)\to 1^{-}$. This case is illustrated in Fig. S15(f), where we
find $\mathbb{H}=\sigma/[\frac{1}{2}(1-m^{2})]=0.98$.

####  S8.3.7 Illustration and test against simulations

Fig. S16 shows $H_{\rm SL}(\tau)$, from Eq. (S125), as an example, for both
positive and negative values of $\tau$.

Figure S16: Plot of $H_{\rm SL}(\tau)$ from Eq. (S125) over a range of
positive and negative values of $\tau$. The horizontal orange dashed line
shows $H=1$, and the vertical green dotted line is at $\tau=-2$ for reference.

In Fig. S17 we show simulation results to illustrate the effects of negative
values of $\tau$, and to demonstrate the validity of our analytical
description. Specifically, simulations are on a square lattice and we show
plots of $\left\langle{\sigma}\right\rangle_{\rm st}$ as a function of
$\left\langle{m}\right\rangle_{\rm st}$. As before these curves are obtained
by varying the detailed model parameters, while keeping $\tau$ fixed [see Sec.
S10.2]. We also show the analytical predictions for the square lattice, i.e.
$\left\langle{\sigma}\right\rangle_{\rm st}=\frac{1}{2}H_{\rm
SL}(\tau)(1-\left\langle{m}\right\rangle_{\rm st}^{2})$ with $H_{\rm
SL}(\tau)$ from Eq. (S125).

For $\tau=1$, the model dynamics contains some transmission noise, but copying
from neighbours is more likely to be faithful than unfaithful. We have
$H(\tau)<1$.

For very large positive $\tau$, noise dominates, and the resulting
configurations at a given $m$ are effectively random, so $H(\tau)=1$.

For $\tau=-3$, copying events are mostly unfaithful (that is, spins take the
state opposite to the one they are ‘copying’). There is also spontaneous
flipping, and some faithful copying. Overall, some degree of anti-correlation
between neighbours ensues, and so $H(\tau)>1$.

The stationary average density of active interfaces at the apexes of the three
parabolas (located at $\left\langle{m}\right\rangle_{\rm st}=0$) are given by
$\left\langle{\sigma}\right\rangle_{\rm st}=\frac{H(\tau)}{2}$, and are
smaller than $0.5$ for $\tau=1$, equal to $0.5$ for $\tau\to\infty$, and
larger than $0.5$ for $\tau=-3$.

Figure S17: Parabolas $\left\langle{\sigma}\right\rangle_{\rm
st}=\frac{1}{2}H_{\rm SL}(\tau)(1-\left\langle{m}\right\rangle_{\rm st}^{2})$
for the model on an infinite square lattice, and for different values of
$\tau$. Dashed lines are from the theory [Eq. (S125)], markers are the results
from averaging the steady-state values of $10$ independent simulations on a
lattice with $N=10,000$ nodes (periodic boundary conditions). For $\tau=-3$
the algorithm we use to generate model parameters [see Sec. S10.2] limits the
values of the stationary magnetisation that we can probe. More specifically,
for parameters generated by this algorithm, $\left\langle{m}\right\rangle_{\rm
st}$ can only take values in the range given in Eq. (S300). For the other
cases $\left\langle{m}\right\rangle_{\rm st}$ can take any value in the
interval $[-1,1]$. The dashed line labelled $\tau\to\pm\infty$ corresponds to
$H(\tau)=1$, the orange markers are from simulations with $\tau=1000$.

##  S9 Analysis of linguistic data

###  S9.1 Dataset

In this section we provide further information on the analysis of data from
WALS (The World Atlas of Language Structures) [22] in the main paper. WALS
contains information on $2,662$ human languages, most importantly the co-
ordinates of where that language is most densely spoken. We will refer to
languages by their WALS-ID, which is a combination of numbers and letters such
as ‘11A’, ‘143F’ etc.

Language possess a number of language traits (known as ‘typological features
of language’). WALS details 192 such traits. Not all traits are recorded for
all languages in WALS. We focus on a curated set of $35$ language traits which
can be quantified in a sufficient set of languages, and which are
‘binarizable’, that is, if the trait can be quantified in the language it can
meaningfully be classified as present or absent. The set of traits and
languages is that of [13], where further details can be found. Ultimately, we
consider $2,248$ languages from WALS, these are those for which at least one
trait is quantifiable in binary form. This means that not all nodes
(languages) have a well defined binary state for all $35$ traits.

Binary traits can all be naturally mapped to our theoretical set-up. Nodes
represent languages, the spin states indicate whether that language posses or
lacks a specific language trait, links between nodes indicate interaction.

The relevant data to perform this analysis including details about traits,
distances between languages, and information on which languages contain which
features can be found at the following GitHub repository [47].

###  S9.2 Empirical noise-to-order ratio

To determine the density of active interfaces for each of the $35$ features
from the WALS data we need to decide how the languages are connected. We do
not know the true interaction network. We therefore make the following choice:
two languages A and B are connected if B is among the $10$ geographically
closest neighbours of A, or if A is among the $10$ geographically closest
neighbours of B. This creates an undirected network.

With this method, languages can end up having more than $10$ neighbours in the
network. For example, start by considering A. We want to determine which other
languages it connects to. We find the $10$ geographically nearest languages of
A and connect it to those languages. Assume that B is one of these languages.
We then consider B. It already has one connection to A. We then find its $10$
geographically nearest languages. It is possible that A is not within this
set. We connect B to these $10$ languages. B now has $11$ neighbours in the
network.

The choice of $10$ nearest neighbours is somewhat arbitrary, however, the
results of [13] for example have been shown to be robust against variation in
neighbourhood size.

The degree distribution for the above network is shown in Fig. S18. The
network is not regular, neither is it uncorrelated as it is likely that two
neighbours of a given node are connected to one another. However the degree
distribution is fairly tight, so we can still apply the various analytical
results to good effect.

Figure S18: Degree distribution of our approximation to the real-world network
of language interactions (see text).

To determine an empirical noise-to-order ratio for a particular binarisable
trait, we assign states $\pm 1$ to those nodes (languages), for which the
binary state of the trait is defined. The state $\pm 1$ indicates whether the
language represented by the node posses or lacks the trait. If the trait is
not quantifiable on a given node, then no spin state is assigned to the node.
Taking into account only nodes with states $\pm 1$, we then determine the
feature frequency $\rho$, and, using the putative interaction network, the
density of active links, $\sigma$. From this we obtain
$\mathbb{H}=\sigma/[2\rho(1-\rho)]$ as a quantifier of the spatial scattering.
Finally, the different analytical functions, $H_{\rm AA}(\tau)$ [Eq. (S110)],
$H_{\rm SL}(\tau)$ [Eq. (S125)], $H_{\rm HG}(\tau)$ [Eq. (S188)] and $H_{\rm
SPA}(\tau)$ [Eq. (S274)], can be used to infer values for the noise-to-order
ratio, $\tau_{X}$. The $X$ subscript stands for one of AA, SL, HG, or SPA
depending on the specific $H_{X}(\tau)$ function that was used.

Using this procedure it is sometimes possible to obtain
$\mathbb{H}=\sigma/[2\rho(1-\rho)]>1$. As discussed in Sec. S8 this would lead
to negative values of $\tau$ values. Such occurrences are rare in our dataset
(less than $1\%$), so we discard these cases.

###  S9.3 Comparison against longitudinal stability estimate

To validate our results we compare to a phylogenetic method which does not
take into account any spatial data. This was done in [49] using a principal
component analysis, and in particular the first component PC1. Higher PC1
values corresponding to greater instability. PC1 values were only measured for
a subset of the 35 traits that we consider.

In Fig. S19 (left) we show regression plots of PC1 against $\tau_{X}$. The
different rows correspond to different $X$ and thus use different analytical
forms $H_{X}(\tau)$, as described above and indicated in the figure.

A similar but restricted analysis was carried out in [13]. This was solely
based on $H_{\rm SL}(\tau)$ from the theory for infinite square lattices. That
earlier analysis identified a small set of outlier traits for which the
estimates of $\tau$ and Dediu’s PC1 gave conflicting indications regarding
relative stability. Some of these outliers could be justified based on known
properties of some of the WALS features, see the discussion in [13] for
details.

We have carried out a similar identification of outliers to test robustness
against the use of different analytical forms $H_{X}(\tau)$. Outliers are
determined by recursively pruning the dataset and removing the trait that
would provide the largest increase in Pearson correlation coefficient $r$
(measuring linear correlation [50]). Results are shown in the right-hand
panels of Fig. S19. Independent of the analytical form used for the function
$H_{X}(\tau)$ we find the same set of outliers, namely WALS features 8A, 11A,
107A and 57A. This is the same set as found in [13], confirming the robustness
against variation of the assumed network structure [reflected by the choice of
analytical form of $H_{X}(\tau)$].

After outlier removal we find strong correlation between the inferred values
of $\tau_{X}$ and Dediu’s PC1, with an $r$ consistently above $0.8$, and
$p$-value $<0.01$, for all theories $X$.

Figure S19: Panels on the left show regression plots of Dediu’s PC1 [49]
against the noise-to-order ratio, $\tau$. Red points are deemed outliers (see
text). The red and blue dashed lines are the best linear fit before and after
outliers are removed (we draw attention to the logarithmic scale on the
horizontal axis, making the fits look curved). The different rows correspond
to different choices of the analytical form of $H(\tau)$, as indicated on the
right-hand side. From top to bottom we have used the annealed approximation
[Sec. S4], the infinite 2D square lattice [Sec. S5], the network walk method
[Sec. S6] and the stochastic pair approximation [Sec. S7] After outlier
removal the $r$ values are $0.813$, $0.857$, $0.848$, $0.821$ from top to
bottom, all with $p$-value$<0.01$. The right column illustrates the outlier
removal process. At each pruning iteration the point that gives the largest
increase in $r$ is removed. For clarity, we only show the first $15$
iterations.

###  S9.4 Bootstrapping and rankograms

####  S9.4.1 Construction of rankograms

The method in Sec. S9.2 can be used to rank the $35$ traits in terms of noise-
to-order ratio, which in turn can be regarded as a proxy for stability.

Our analysis is based on the WALS dataset, involving $N=2,248$ languages, each
with at least one binarisable feature among our selected set of $35$ traits.
In reality, there are other languages for which we do not have the binarisable
data. This can be due to lack of observation, or WALS simply is missing the
data. We can therefore think of the WALS dataset as an effective random sub-
sample of all human languages. To quantify the uncertainty associated with
this, we use bootstrap re-sampling [51].

In our context, bootstrapping involves sampling $N=2,248$ languages with
replacement from the original WALS dataset. This gives rise to an (assumed)
interaction network as described in Sec. S9.2. This network is specific to the
bootstrap sample. We note that is possible to select the same language
multiple times. This leads to multiple identical languages with geographical
distance zero, thus they are essentially guaranteed to be connected. This
represents potentially very close neighbouring languages, not included in
WALS, which have the same trait values.

For each sample, we can perform the procedure in Sec. S9.2 to calculate a
noise-to-order ratio $\tau_{X}$ for each trait (discarding nodes in the
network for which the status of a given trait cannot be quantified in binary
form). Repeating this many times, we generate distributions of $\tau_{X}$ for
each trait.

To produce rankograms for a chosen $X$ we can then draw a value of $\tau_{X}$
for each trait from the respective distributions. The traits can then be
ranked based off of these values. This is then repeated many times and thus a
distribution of ranks is constructed for each trait. These are the rankograms
shown in Figs. 4 in the main paper, and in S20 in this supplement. These are
histograms where the $i^{\rm th}$ bar indicates the probability that the trait
represented by that rankogram ranks in the $i^{\rm th}$ position amongst all
traits.

Figure S20: Rankograms for all 35 language traits for 10,000 bootstrap
samples. We only show Gaussian kernel density estimates of the distributions
for clarity. Each rankogram has five distributions. Four of which are based on
inferring values of $\tau_{X}$ from the theories $H_{\rm AA}(\tau)$ [Eq.
(S110)], $H_{\rm SL}(\tau)$ [Eq. (S125)], $H_{\rm HG}(\tau)$ [Eq. (S188)] and
$H_{\rm SPA}(\tau)$ [Eq. (S274)]. The fifth is from ranking raw $\mathbb{H}$
values calculated using Eq. (2). The height at a given rank is the probability
a trait is ranked in that position amongst all traits.

We observe that the traits on the extreme ends, such as 130A or 53A, have very
tight rankograms. This means that we can be confident that those traits will
rank in those positions amongst all features. We are less confident about the
traits with broad distributions, such as 100A or 48A.

####  S9.4.2 Distributions of $\tau_{X}$

From the rankograms in Fig. S20 we see that the rankings produced from
$\mathbb{H}$ and those produced from $\tau_{X}$ using each of the theories $X$
are very similar. This is a non-trivial result and is not a consequence of the
fact that all the $H_{X}(\tau)$ functions are monotonic increasing functions
of $\tau$ as explained below.

In a single bootstrap sample, we have a set of languages and construct one
interaction network. We then measure 35 values of $\mathbb{H}$, one for each
trait, using the procedure in Sec. S9.2. These $\mathbb{H}$ values will have a
certain ranking. With the single interaction network we construct the
functions $H_{X}(\tau)$, one for each theory $X$. For a chosen $X$ the 35
$\mathbb{H}$ values can be mapped to 35 $\tau_{X}$ values. Since all
$H_{X}(\tau)$ are monotonic increasing functions of $\tau$, ordering the
traits based on $\tau_{X}$ will be the same as ordering based on $\mathbb{H}$,
for all theories $X$.

However, the procedure for constructing rankograms is different in that we do
not perform any rankings until distributions of $\tau_{X}$ have been
constructed for each trait. To construct distributions of $\tau_{X}$ for a
trait we repeat the procedure in Sec. S9.2 many times. Each bootstrap sample
consists of a different set of languages, and thus a different interaction
network. Focusing on two traits $A$ and $B$, we measure the corresponding
values $\mathbb{H}^{A}=\sigma^{A}/[2\rho^{A}(1-\rho^{A})]$ and
$\mathbb{H}^{B}=\sigma^{B}/[2\rho^{B}(1-\rho^{B})]$ in each sample. This gives
distributions for $\mathbb{H}^{A}$ and $\mathbb{H}^{B}$. For each bootstrap
sample we also map the empirical values of $\mathbb{H}^{A}$ and
$\mathbb{H}^{B}$ to values $\tau^{A}_{X}$ and $\tau^{B}_{X}$, using one of the
analytical functions $H_{X}(\tau)$. This gives us distributions of
$\tau^{A}_{X}$ and $\tau^{B}_{X}$.

Since the interaction network is different in each bootstrap sample, the
function $H_{X}(\tau)$ will also change in each sample (for a fixed choice of
$X$). This is because the functional form of $H_{X}(\tau)$ makes use for
example of the mean degree or other characteristics of the network. This
means, that even if the distributions of $\mathbb{H}^{A}$ and $\mathbb{H}^{B}$
do not overlap, for a given $X$ it is possible that the distributions of
$\tau^{A}_{X}$ and $\tau^{B}_{X}$ do overlap. Fig. S21 gives a visual
representation of this idea.

Figure S21: Pictorial demonstration of how non-overlapping $\mathbb{H}$
distributions from two traits A and B can lead to $\tau_{X}$ distributions
which do overlap through the bootstrapping process (see text). The yellow,
green and red curves are examples of $H_{X}(\tau)$ from different bootstrap
samples.

###  S9.5 SUCRA scores and comparison with Dediu

To quantify the rankograms for each of the 35 linguistic traits shown in Fig.
S20 further we use a statistic called SUCRA (surface under the cumulative
ranking line) first introduced in [52] to rank medical treatments. This is a
transform of the expected rank,

$\text{SUCRA}_{a}=\frac{N-\mathbb{E}(r)_{a}}{N-1},$ (S288)

where $\mathbb{E}(r)_{a}$ is the expected rank of item $a$ and where we are
ranking $N$ items total. SUCRA values go from 0 to 1 but are typically written
as a percentage. If item $a$ ranks first with probability 1 then it will have
$\text{SUCRA}_{a}=1$ (or 100%). Similarly if item $a$ ranks last with
probability 1 then it will have $\text{SUCRA}_{a}=0$ (or 0%).

In Tab. S1 we list the SUCRA scores for all $35$ traits for the different
choices of the function $H_{X}(\tau)$ as well as the SUCRA derived from a
ranking based on $\mathbb{H}=\sigma/[2\rho(1-\rho)]$ itself. We also list
Dediu’s PC1 values where possible. We can calculate the Spearman rank
correlation coefficient, $r_{S}$, between the traits sorted by SUCRA value and
PC1. All give $r_{S}>0.8$ with $p$-value$<0.01$ after removing the outliers.

WALS ID | SUCRA values (%) |  Dediu’s PC1
---|---|---
|  $\mathbb{H}$ |  $H_{\rm AA}(\tau)$ |  $H_{\rm SL}(\tau)$ |  $H_{\rm HG}(\tau)$ |  $H_{\rm SPA}(\tau)$ |
130A |  96.97 |  96.97 |  97.02 |  96.90 |  97.03 |  97.03
92A |  94.52 |  94.36 |  94.33 |  94.60 |  94.44 |  87.10
11A |  88.42 |  88.07 |  88.43 |  87.67 |  88.12 |  -101.68
143F |  88.02 |  87.99 |  87.87 |  87.38 |  88.01 |
143E |  85.30 |  85.61 |  85.42 |  84.97 |  85.54 |
33A |  79.73 |  79.70 |  79.66 |  79.57 |  79.66 |
38A |  79.49 |  79.12 |  79.23 |  79.56 |  79.34 |  76.40
37A |  78.77 |  78.77 |  78.64 |  78.77 |  78.77 |  87.38
8A |  77.89 |  77.25 |  77.42 |  77.93 |  77.04 |  -57.03
77A |  73.38 |  73.26 |  73.10 |  73.01 |  73.01 |  50.18
69A |  72.26 |  72.28 |  72.26 |  72.26 |  72.26 |
4A |  66.43 |  66.44 |  66.52 |  66.43 |  66.40 |  -16.94
26A |  61.30 |  61.49 |  61.35 |  62.21 |  61.60 |
100A |  58.05 |  57.77 |  57.63 |  57.21 |  57.47 |
57A |  55.69 |  55.90 |  55.76 |  55.51 |  55.83 |  40.60
6A |  52.71 |  52.75 |  52.96 |  52.70 |  52.77 |  -68.03
7A |  50.79 |  51.04 |  50.97 |  51.34 |  51.14 |  -42.87
129A |  50.69 |  50.65 |  50.68 |  50.82 |  50.62 |
120A |  45.45 |  45.55 |  45.79 |  45.46 |  45.80 |
48A |  44.79 |  45.16 |  44.99 |  45.16 |  45.32 |  18.44
44A |  44.47 |  44.48 |  44.68 |  44.08 |  44.27 |  -8.95
107A |  41.41 |  41.64 |  41.23 |  41.44 |  41.50 |  79.84
70A |  40.29 |  40.27 |  40.03 |  40.48 |  40.17 |  -23.10
87A |  33.27 |  33.33 |  33.30 |  33.65 |  33.23 |  -11.23
73A |  31.48 |  31.71 |  31.86 |  31.67 |  31.86 |  -64.56
119A |  31.42 |  31.37 |  31.35 |  31.98 |  31.21 |
82A |  25.51 |  25.54 |  25.47 |  25.85 |  25.51 |  -69.15
13A |  17.81 |  18.01 |  18.04 |  17.98 |  18.04 |  -35.89
9A |  17.76 |  17.80 |  17.80 |  18.05 |  18.05 |  -20.70
91A |  17.15 |  17.28 |  17.28 |  17.69 |  17.28 |  -23.96
27A |  16.65 |  16.94 |  16.85 |  16.69 |  16.85 |
83A |  12.20 |  12.27 |  12.15 |  12.07 |  12.20 |  -28.65
89A |  11.62 |  11.59 |  11.65 |  11.81 |  11.54 |  -30.62
86A |  6.82 |  6.99 |  6.87 |  7.07 |  6.99 |  -60.98
53A |  1.23 |  1.19 |  1.10 |  1.29 |  1.13 |  -43.21
Table S1: Table showing the SUCRA values for each of the 35 rankograms in Fig.
S20 for all theories. Dediu’s PC1 values are also shown where possible. The
Spearman rank correlation coefficients, $r_{S}$, between each column and
Dediu’s PC1 is 0.824, except $H_{\rm SPA}(\tau)$ which is 0.832, all with
$p$-value$<0.01$.

To determine the outliers we recursively prune the sample. At each step we
remove the trait that would give the largest increase in $r_{s}$. We show this
in Fig. S22. We find that the four traits 11A, 8A, 6A and 7A give significant
improvement when removed. This is a different ordering than in Fig. S19, here
it seems that 6A and 7A are outliers whereas 107A and 57A are not. This is
likely simply due to using different pruning methods. It is encouraging
however that 11A and 8A are still outliers which is the same as before. In
[13] there are strong linguistics arguments as to why Dediu incorrectly
classifies 11A.

Figure S22: Visual representation of outlier removal when determining the
Spearman rank correlation coefficient, $r_{S}$, between the linguistic traits
ranked by SUCRA value and Dediu’s PC1 [see Tab. S1]. We recursively prune the
data and remove the feature which would give the largest improvement in
$r_{S}$.

##  S10 Further mathematical background

###  S10.1 Exponential generating function for modified Bessel functions of
the first kind

The exponential generating function for the modified Bessel function of the
first kind is [53, 49:6:1]

$I_{0}(2x)=\sum_{k=0}^{\infty}\frac{x^{2k}}{(k!)^{2}}=\frac{x^{0}}{(0!)^{2}}+\frac{x^{2}}{(1!)^{2}}+\frac{x^{4}}{(2!)^{2}}+....$
(S289)

Squaring this equation provides us with the exponential generating function
for $I_{0}(2x)^{2}$,

$\displaystyle I_{0}(2x)^{2}$
$\displaystyle=\left(\frac{1}{(0!)^{2}}\frac{1}{(0!)^{2}}\right)x^{0}+\left(2\frac{1}{(0!)^{2}}\frac{1}{(1!)^{2}}\right)x^{2}$
$\displaystyle\qquad+\left(2\frac{1}{(0!)^{2}}\frac{1}{(2!)^{2}}+\frac{1}{(1!)^{2}}\frac{1}{(1!)^{2}}\right)x^{4}+...$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{a_{1}+a_{2}=k}\frac{x^{2k}}{(a_{1}!)^{2}(a_{2}!)^{2}}.$
(S290)

To generalise, first define the combinatoric sum

$\sum_{a_{1}+...+a_{d}=k}\equiv\sum_{A_{d}^{k}}.$ (S291)

In general then we can form the exponential generating function for
$I_{0}(2x)^{d}$,

$\displaystyle I_{0}(2x)^{d}$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{A_{d}^{k}}\frac{x^{2k}}{\prod_{i}^{d}(a_{i}!)^{2}}$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{A_{k}^{d}}\frac{(2k)!}{\prod_{i}^{d}(a_{i}!)^{2}}\frac{x^{2k}}{(2k)!}.$
(S292)

Now recall the expression for $W_{d}^{(2k)}$ from equation (S196),

$W_{d}^{(2k)}=\sum_{A_{k}^{d}}\frac{(2k)!}{\prod_{i=1}^{d}(a_{i}!)^{2}},$
(S293)

which gives us the final expression

$I_{0}(2x)^{d}=\sum_{k=0}^{\infty}W_{d}^{(2k)}\frac{x^{2k}}{(2k)!}.$ (S294)

###  S10.2 Determining standard model parameters for a given $\tau$ and
$\left\langle{m}\right\rangle_{\rm st}$

In this section we describe the algorithm we used to generate model parameters
in the simulation. The challenge is to generate sets of non-negative
$h_{\pm},h^{\prime}_{\pm},v^{\prime}_{\pm}$ leading to different stationary
trait frequencies $\left\langle{\rho}\right\rangle_{\rm st}$ but at the same
time keeping
$\tau=2(h_{+}^{\prime}+h_{-}^{\prime}+v_{+}^{\prime}+v_{-}^{\prime})/[(h_{+}+h_{-})-(h_{+}^{\prime}+h_{-}^{\prime})]$
fixed. Fixing $\left\langle{\rho}\right\rangle_{\rm st}$ is equivalent to
fixing the stationary magnetisation $\left\langle{m}\right\rangle_{\rm
st}=2\left\langle{\rho}\right\rangle_{\rm st}-1$.

We apply the restriction in Eq. (S5), which we write in the form

$h_{+}=h_{+}^{\prime}-h_{-}^{\prime}+h_{-}.$ (S295)

Equations (S24) and (S46) can be solved for $v_{+}^{\prime}$ and
$v_{-}^{\prime}$, assuming $\tau\neq 0$

$\displaystyle
v_{+}^{\prime}=\frac{\tau}{2}(1+\left\langle{m}\right\rangle_{\rm
st})(h_{-}-h_{-}^{\prime})-h_{+}^{\prime},$ (S296a) $\displaystyle
v_{-}^{\prime}=\frac{\tau}{2}(1-\left\langle{m}\right\rangle_{\rm
st})(h_{-}-h_{-}^{\prime})-h_{-}^{\prime}.$ (S296b)

####  S10.2.1 The case $\tau>0$

For $\tau>0$ and $-1<\left\langle{m}\right\rangle_{\rm st}<1$, we set
$h_{+}^{\prime}$ and $h_{-}^{\prime}$ to small random values and then choose a
value of $h_{-}$ large enough so that $v_{+}^{\prime}$, $v_{-}^{\prime}$ and
$h_{+}$ are non-negative. This amounts to choosing $h_{-}$ such that

$h_{-}\geq\text{max}\left(0,h_{-}^{\prime}-h_{+}^{\prime},h_{-}^{\prime}+\frac{2}{\tau(1+\left\langle{m}\right\rangle_{\rm
st})}h_{+}^{\prime},\left[1+\frac{2}{\tau(1-\left\langle{m}\right\rangle_{\rm
st})}\right]h_{-}^{\prime}\right).$ (S297)

For the case $\tau>0$ and $\left\langle{m}\right\rangle_{\rm st}=-1$, Eq.
(S296a) means that we need to set $v_{+}^{\prime}=h_{+}^{\prime}=0$. Eq.
(S296b) then gives $v_{-}^{\prime}=\tau h_{-}-(1+\tau)h_{-}^{\prime}$. So to
ensure $v_{-}^{\prime}\geq 0$ we draw random values of $h_{-}$ and
$h_{-}^{\prime}$ such that
$h_{-}\geq\left(1+\frac{1}{\tau}\right)h_{-}^{\prime}$. This allows us to
calculate $v_{-}^{\prime}$. Finally $h_{+}$ is determined from Eq. (S295).

For the case $\tau>0$ and $\left\langle{m}\right\rangle_{\rm st}=+1$, Eq.
(S296b) means that we should set $v_{-}^{\prime}=h_{-}^{\prime}=0$. Eq.
(S296a) then gives $v_{+}^{\prime}=\tau h_{-}-h_{+}^{\prime}$. So to ensure
$v_{-}^{\prime}\geq 0$ we draw random values of $h_{-}$ and $h_{+}^{\prime}$
such that $h_{-}\geq\tau h_{+}^{\prime}$. This allows us to calculate
$v_{+}^{\prime}$. The coefficient $h_{+}$ is again determined from Eq. (S295).

####  S10.2.2 The case $\tau=0$ (standard voter model)

For $\tau=0$ Eqs. (S296a) and (S296b) dictate that we need to set the rates of
all unfaithful events to zero
($h_{+}^{\prime}=h_{-}^{\prime}=v_{+}^{\prime}=v_{-}^{\prime}=0$). We then
have $h_{+}=h_{-}$ from Eq. (S295). The model therefore reduces to the
standard voter model. The common value of the rates $h_{\pm}$ can be drawn at
random (the only significance of the value it to set a time scale). We note
however that in this case the expression for
$\left\langle{m}\right\rangle_{\rm st}$, Eq. (S24), is not defined. In the
voter model average magnetisation is conserved in time, and every initial
magnetisation is therefore a steady-state value
($\left\langle{m}\right\rangle_{\rm st}=\left\langle{m}\right\rangle_{0}$).

####  S10.2.3 Negative values of $\tau$ ($\tau\leq-2$)

As discussed in Sec. S8 the quantity $\tau$ cannot take values in the interval
$-2<\tau<0$, that is, the only negative values of $\tau$ permitted by the
model setup are those in the range $\tau\leq-2$. We can confirm this through
numerical search for model parameters. More precisely, we conducted searches
of parameter space, enforcing the constrain in Eq. (S5) and fixing
$\left\langle{m}\right\rangle_{\rm st}$ and $\tau$, in Eqs. (S24) and (S46)
respectively. For $-2<\tau<0$ the search returns no solutions (see the
Mathematica notebook at the following GitHub repository [47] for a
demonstration).

Focusing on $\tau\leq-2$, a simple method to generate model parameters is to
set $h_{+}^{\prime}=h_{-}^{\prime}$ and then draw this number at random.
Consequently, from Eq. (S295) we have $h_{+}=h_{-}$. Then, to ensure
$v_{+}^{\prime}\geq$ and $v_{-}\geq 0$, Eqs. (S296a) and (S296b) give the
following bounds on $h_{-}$:

$0\leq
h_{-}\leq\text{min}\left(1+\frac{2}{\tau(1+\left\langle{m}\right\rangle_{\rm
st})},1+\frac{2}{\tau(1-\left\langle{m}\right\rangle_{\rm
st})}\right)h_{-}^{\prime}.$ (S298)

To ensure the upper bound here is at least zero we require

$\tau\leq\text{min}\left(\frac{-2}{1+\left\langle{m}\right\rangle_{\rm
st}},\frac{-2}{1-\left\langle{m}\right\rangle_{\rm st}}\right).$ (S299)

This means that for a given $\tau\leq-2$, $\left\langle{m}\right\rangle_{\rm
st}$ must be within the following bounds.

$-1-\frac{2}{\tau}\leq\left\langle{m}\right\rangle_{\rm st}\leq
1+\frac{2}{\tau}.$ (S300)

A visual representation of these bounds is shown in Fig. S23.

The inequalities in Eq. (S300) bound the possible values of
$\left\langle{m}\right\rangle_{\rm st}$ that can be reached for a given $\tau$
by model parameters generated using our algorithm. In particular we used the
restriction $h_{+}=h_{-}$ to derive these bounds. We can therefore not exclude
that models with $h_{+}\neq h_{-}$ can cover a wider range of
$\left\langle{m}\right\rangle_{\rm st}$ for a given $\tau$. However, numerical
searches (allowing for $h_{+}\neq h_{-}$) find that the property from Eq.
(S300) seems to continue to hold [47].

We here provide further intuitive support for these bounds. From Sec. S8 we
know that $\tau\to-2^{-}$ corresponds to increasing amounts of anti-
correlation between the spin states of nearest-neighbour nodes. This requires
an increasingly equal proportion of up and down spin nodes [see Fig. S15(e)].
At $\tau=-2$ this requires $\left\langle{m}\right\rangle_{\rm st}=0$, and as
$\tau$ becomes more negative from there the range of possible
$\left\langle{m}\right\rangle_{\rm st}$ values increases.

Figure S23: Diagram showing the values of $\tau\leq-2$ and
$\left\langle{m}\right\rangle_{\rm st}$ for which it is possible to generate
model parameters. The boundaries are given by Eq. (S299).