A multi-state voter model with imperfect copying

TL;DR

This paper introduces a multi-state voter model with imperfect copying, showing how noise affects order in large populations and demonstrating that perfect order is unattainable with any imperfection.

Contribution

It develops a mean-field model incorporating copying errors, revealing how noise prevents complete ordering as population size grows.

Findings

Order decreases with population size for any non-zero noise

Complete order is impossible when copying errors are present

Monte Carlo simulations support analytical results

Abstract

The voter model with multiple states has found applications in areas as diverse as population genetics, opinion formation, species competition and language dynamics, among others. In a single step of the dynamics, an individual chosen at random copies the state of a random neighbor in the population. In this basic formulation it is assumed that the copying is perfect, and thus an exact copy of an individual is generated at each time step. Here we introduce and study a variant of the multi-state voter model in mean-field that incorporates a degree of imperfection or error in the copying process, which leaves the states of the two interacting individuals similar but not exactly equal. This dynamics can also be interpreted as a perfect copying with the addition of noise; a minimalistic model for flocking. We found that the ordering properties of this multi-state noisy voter model, measured by a parameter $\psi$ in [0, 1], depend on the amplitude $\eta$ of the copying error or noise and the population size N. In the case of perfect copying $\eta=0$ the system reaches an absorbing configuration with complete order ($\psi=1$) for all values of N. However, for any degree of imperfection $\eta>0$, we show that the average value of $\psi$ at the stationary state decreases with N as $\langle \psi \rangle \simeq 6/(\pi^2 \eta^2 N)$ for $\eta \ll 1$ and $\eta^2 N \gtrsim 1$, and thus the system becomes totally disordered in the thermodynamic limit $N \to \infty$. We also show that $\langle \psi \rangle \simeq 1-1.64 \, \eta^2 N$ in the vanishing small error limit $\eta \to 0$, which implies that complete order is never achieved for $\eta > 0$. These results are supported by Monte Carlo simulations of the model, which allow to study other scenarios as well.

Equations

$$s(t)=\left(\frac{t}{N-1}+\frac{1}{S}\right)^{-1}~{}~{}~{}\mbox{for $N\geq S\gg 1$}$$

$$\psi (t)={\left|\frac{1}{N}\sum _{m=1}^N{e}^{i{\Theta }_m(t)}\right|}^2={\left|\sum _{k=0}^{S-1}{x}_k(t)e^{i{\theta }_k}\right|}^2,$$

$${\sigma }_{\theta }^2(t)$$

$$\overline{\theta }(t)$$

$${\sigma }_u^2=\frac{4{\pi }^2}{S^3}\sum _{k=0}^{S-1}{k}^2-{\left[\frac{2\pi }{S^2}\sum _{k=0}^{S-1}k\right]}^2=\frac{{\pi }^2(S^2-1)}{3S^2},$$

$$\frac{d}{dt}P(\{x\},t)=\sum_{k=0}^{S-1}\sum_{\begin{subarray}{c}j=0\atop j\neq k\end{subarray}}^{S-1}\Big{\{}W\left(\{x^{\prime}\}_{k\;j}^{-+}\to\{x\}\right)\,P\left(\{x^{\prime}\}_{k\;j}^{-+},t\right)-W\left(\{x\}\to\{x^{\prime}\}_{k\;j}^{-+}\right)\,P\left(\{x\},t\right)\Big{\}}.$$

$$W_{kj}^{\downarrow \uparrow }\left(\{x\}\right)$$

$$W_{kj}^{\uparrow \downarrow }\left(\{x\}\right)$$

$$W_{01}^{\downarrow \uparrow }\left(\{x_0,x_1\}\right)$$

$$W_{01}^{\uparrow \downarrow }\left(\{x_0,x_1\}\right)$$

$$W_{kj}^{\uparrow \downarrow }\left(\{x^{\prime }{\}}_{kj}^{-+}\right)P\left(\{x^{\prime }{\}}_{kj}^{-+},t\right)$$

$${\partial }_{t}P(\{x\},t)$$

$$\sum _{k=0}^{S-1}\sum _{\begin{array}{c}\genfrac{}{}{0ex}{}{j=0}{j\ne k}\end{array}}^{S-1}{\partial }_j\left(W_{kj}^{\uparrow \downarrow }P\right)$$

$$\sum _{k=0}^{S-1}\sum _{\begin{array}{c}\genfrac{}{}{0ex}{}{j=0}{j\ne k}\end{array}}^{S-1}{\partial }_{kj}^2\left(W_{kj}^{\uparrow \downarrow }P\right)$$

$${\partial }_{t}P(\{x\},t)$$

$$A_k$$

$$B_{kk}$$

$$B_{kj}$$

$$\frac{\partial}{\partial t}P(x_{0},t)=\frac{1}{2}\frac{\partial}{\partial x_{0}}\left[(2x_{0}-1)P(x_{0},t)\right]+\frac{1}{4N}\frac{\partial^{2}}{\partial x_{0}^{2}}P(x_{0},t),~{}~{}\mbox{for $S=2$}.$$

$$P_{st}(x_0)=\frac{\sqrt{2N}{e}^{-2N(x_0-1/2{)}^2}}{\sqrt{\pi }\text{mbox}erf\left(\sqrt{\frac{N}{2}}\right)},$$

$$\psi (t)={\left|x_0(t)-x_1(t)\right|}^2={\left[2x_0(t)-1\right]}^2.$$

$$⟨\psi ⟩$$

$$⟨\psi ⟩\simeq \frac{1}{N},$$

$${\sigma }_{\theta }^2(t)={\pi }^2{x}_1(t)[1-x_1(t)]={\pi }^2{x}_0(t)[1-x_0(t)],$$

$$⟨{\sigma }_{\theta }^2⟩$$

$$⟨{\sigma }_{\theta }^2⟩\simeq \frac{{\pi }^2}{4}\left(1-\frac{1}{N}\right)={\sigma }_u^2\left(1-\frac{1}{N}\right),$$

$$\psi (t)={\left(\sum _{k=0}^{S-1}{x}_k(t)\cos {\theta }_k\right)}^2+{\left(\sum _{k=0}^{S-1}{x}_k(t)\sin {\theta }_k\right)}^2.$$

$$\psi (t)$$

$$⟨\psi ⟩$$

$$⟨\psi ⟩=S{z}_0+2\sum _{n=1}^{S-1}(S-n)z_n\cos (2\pi n/S).$$

$$z_0$$

$$z_n$$

$$⟨\psi ⟩=2(z_0-z_1)\simeq \frac{1}{N},$$

$$⟨\psi ⟩$$

Full text

I Introduction

Stochastic models of evolution have been successfully applied in various disciplines to study the dynamics of systems composed by many interacting entities such as genes in population genetics, animal or plant species in ecology, and people in linguistics and sociology, among others (see Blythe and McKane (2007) for a statistical physics review). The most basic –neutral– version of each of these models implements some type of copying mechanism by which an entity is removed and replaced by an exact copy of another entity in the population. For instance, in a single step of the Moran model Moran (1958) for genetic drift (similar to the Wright-Fisher model Fisher (1930); Wright (1931)) a gene is chosen at random to die and replaced by a new gene that is a replica of another gene in the population, its “parent”, also chosen at random. Similarly, neutral models for the evolution of species in ecology consider that when a tree dies is replaced by an “offspring” of a randomly chosen tree in the forest Hubbell (2001). A theory analogous to that of population genetics was presented in Baxter et al. (2006) to explore the dynamics of language change in the context of linguistic variables, such as vowel sound or grammar. The copying mechanism is also used in the voter model for opinion formation Clifford and Sudbury (1973); Holley and Liggett (1975), where each individual adopts the opinion of one of its neighbors in the population. More recently, this type of social imitation rule was introduced to study the flocking dynamics of a large group of animals Baglietto and Vazquez (2018), for instance birds, where each bird aligns its flying direction with that of a nearby random bird. In the case of all-to-all interactions, this flocking voter model is equivalent to the well known multi-state voter model (MSVM) Starnini et al. (2012); Pickering and Lim (2016) for opinion dynamics, where the moving direction of a bird is associated to its opinion or decision. The MSVM considers a population composed by a fixed number of agents (voters) subject to pairwise interactions, where each voter can hold one of $S$ possible states that represent different opinions or positions on a given issue. In a single step of the dynamics, a voter chosen at random updates its state by copying the state of another agent randomly chosen in the population. The MSVM assumes that the copying process is perfect, in the sense that once an agent copies the state of its partner these two agents are considered to be indistinguishable. However, in a real life situation one would expect some degree of inaccuracy in the copying process that translates into an imperfect copying. For instance, a person can try to adopt the exact opinion of a partner on a given opinion spectrum, but the imitation may not be perfect and the agent ends up taking an opinion very similar but not equal to that of its partner. The source of error in the copying process may also come from the fact that the perception of a person on its partner’s opinion may not be completely accurate.

The imperfect social imitation was recently modeled by adding an external noise in the original voter model, to study the outcome of electoral processes Fernández-Gracia et al. (2014). The original noisy $2$-state voter model assumes that, besides the copying dynamics, voters can randomly switch state. This variant of the model was introduced independently some years ago to study phenomena as diverse as heterogeneous catalytic chemical reactions Fichthorn et al. (1989); Considine et al. (1989), herding behavior in financial markets Kirman (1993) and species competition in probability theory Granovsky and Madras (1995). The study of the effects of noise in the voter model has lately gained attention in the physics literature. Recently, the $2$-state noisy voter model has been explored in complex networks Carro et al. (2016); Peralta et al. (2018a, b), and its dynamics has also been investigated under the presence of zealots Khalil et al. (2018) and the influence of contrarians Khalil and Toral (2019). In Diakonova et al. (2015) the authors have found that noise changes the properties of the fragmentation transition observed in a coevolving version of the voter model Vazquez et al. (2008); Demirel et al. (2014) and the MSVM on complex networks Böhme and Gross (2012).

A mechanism of imperfect imitation was implemented in Roca, C. P. et al. (2009) within a game theory model to study the dynamics of cooperation, where the process of adopting the strategy of a neighboring player combines two different imitation dynamics, the unconditional imitation and the replicator rule Szabó and Fáth (2007). They found that cooperation is enhanced when the probability of choosing the replicator rule (the perturbation) adopts intermediate values. In the context of flocking dynamics, it is reasonable to assume that birds make an error when trying to align with a close by bird, which is modeled by adding a small perturbation (noise) to the alignment process as in Vicsek-type models Vicsek et al. (1995); Baglietto and Albano (2009). It is observed that the noise amplitude induces a transition from a –nematically– ordered phase for low noise to a disordered phase for high noise.

In this article we study a system of interacting particles subject to a multi-state voter dynamics with imperfect copying on a complete graph (all-to-all interactions). For concreteness we use the language of flocking, where the states of particles represent a finite set of angular directions equally spaced in the interval $[0,2\pi)$. In a single iteration step of the dynamics, a particle chosen at random adopts a state that is contained in an interval centered in the state of another randomly chosen particle. Thus, the level of the imperfection in the imitation process is given by the length of the error interval, which is a variable of the model. We note that the update of a particle’s state can also be thought as a two-step process where, in a first step, the particle copies the state of another particle and then, in a second step, its state is perturbed within an interval (spontaneous transitions between states). Although the MSVM with imperfect copying studied here falls in the category of the noisy $2$-state voter models mentioned above, it exhibits some crucial differences with them. That is, the multiplicity of states in the MSVM allows for different types of spontaneous transitions between states, which go beyond the stochastic transition in binary models. Specifically, we consider a system where states are ordered (a discrete set of angles ordered from [math] to $2\pi$) and noise-induced transitions are allowed only between neighboring states, and not between any two states as in most genetic models with mutations Blythe and McKane (2007).

We investigate the ordering dynamics of the system by numerical simulations and analytical techniques and found that the imperfection in the copying mechanism changes completely the ordering properties. When imitation is perfect the system reaches a state of complete order where all particles share the same state, as it happens in the original MSVM. In contrast, the addition of imperfection in the imitation rule reduces order to a level that decreases with the number of particles, leading to complete disorder in the thermodynamic limit even in the case of an infinitesimal error interval. These conclusions are supported by two complementary analytical approaches that provide accurate expressions for the order parameters in the large population limit and in the small error amplitude limit.

The article is organized as follows. We introduce the model and define its dynamics in section II. Section III presents some simulation results showing the qualitative behavior of the model for different parameter values. In sections IV and V we develop two different analytical approaches that show the scaling of macroscopic quantities in different regimes. Finally, in section VI we conclude and summarize our results.

II The Model

We consider a system of $N$ interacting particles that can take one of $S$ possible angular states $\theta_{k}=\frac{2\pi k}{S}$, with $k=0,1,2,...,S-1$, which represent their moving directions. Initially, each particle $j$ ($j=1,..,N$) adopts a state $\Theta_{j}=\theta_{k}$ at random, leading to a nearly uniform distribution of particles on a discrete angular space contained in the interval $[0,2\pi)$. In a single time step $\delta t=1/N$ of the dynamics, a particle $i$ is picked at random and its state $\Theta_{i}$ is updated according to these two steps: first, another particle $j$ with state $\Theta_{j}=\theta_{k}$ is chosen at random and, second, particle $i$ randomly adopts a state $\theta_{l}$ in the interval $[\theta_{k}-\frac{2\pi\Delta}{S},\theta_{k}+\frac{2\pi\Delta}{S}]$ centered at $\theta_{k}$, i e., with equal probability $1/(1+2\Delta)$. That is, $\Theta_{i}(t)\to\Theta_{i}(t+1/N)=\theta_{l}$ in the interval $[\theta_{k-\Delta},\theta_{k+\Delta}]$ (see Fig. 1). Here $\Delta$ is a non-negative integer parameter that defines the amplitude of the error interval ($0\leq\Delta\leq S/2$). The first step corresponds to the selection of a particle $j$ whose state is tried to be imitated by particle $i$, while the second step describes the error making in the copying process, where $i$ adopts a state similar or equal to the state of $j$. In this last step we implement periodic boundary conditions to keep the states in the $[0,2\pi)$ interval, i e., $\theta_{-k}=\theta_{S-k}$ and $\theta_{S+k}=\theta_{k}$ ($0\leq k\leq S-1$), and thus we can think the angle space as a chain ring of $S$ sites at positions $\theta_{k}$. This dynamics can also be interpreted as a perfect copying with the addition of noise, where particle $i$ first jumps to a site at position $\theta_{k}$ and then from there it jumps to any of its $2\Delta$ neighboring sites or stay in the same site with the same probability $1/(1+2\Delta)$ (see Fig. 1).

For the noiseless case $\Delta=0$ the model is equivalent to the MSVM recently studied in the literature Starnini et al. (2012); Pickering and Lim (2016) where, in the above example, particle $i$ simply jumps to the site occupied by particle $j$ and stays there. In this case, given that the system is only driven by the stochastic nature of the copying process (the so called genetic drift in population genetics), a site that becomes empty remains empty afterwards, as particles can jump to occupied sites only. Therefore, the number of sites occupied by at least one particle decreases monotonically with time until only one site becomes occupied by all particles and the system stops evolving. This configuration in which all particles share the same state –a “consensus” in the moving direction– is absorbing, and thus the system has $S$ different absorbing configurations (fixation). A magnitude of interest, which is also relevant in the analysis performed in section V, is the mean number of different states (occupies sites) in the system at time $t$, $s(t)$. It was shown in Starnini et al. (2012); Pickering and Lim (2016) that if $N$ particles are initially distributed homogeneously on the $S$ states [$s(0)=S$], then $s$ decays with time as

[TABLE]

up to a time of the order $N/2$ ($s\simeq 2$), after which $s(t)$ decays exponentially fast to $1.0$. The expected time to reach consensus can be estimated from Eq. (1) as the moment $\tau_{c}$ when $s(\tau_{c})$ becomes $1.0$, leading to the approximate mean consensus time $\tau_{c}\simeq(N-1)(S-1)/S$ Starnini et al. (2012); Pickering and Lim (2016).

Our aim is to study how the addition of imperfection ($\Delta>0$) affects the ordering properties of the system. We start by showing in section III results from Monte Carlo simulations of the dynamics, and then in sections IV and V we develop analytical approaches to gain an insight into these results.

III Simulation results

We simulated the dynamics of the model starting from a configuration in which each particle adopts one of the $S$ angular states $\theta_{k}$ at random and then evolves following the interaction rules defined in section II. The state of the system at a given time $t$ can be described by the set of $S$ variables $\{x\}(t)\equiv\{x_{0}(t),x_{1}(t),..,x_{S-1}(t)\}$, where $x_{k}(t)$ (with $k=0,..,S-1$) is the fraction of particles with state $\theta_{k}=\frac{2\pi k}{S}$ (at site $k$) at time $t$. As the total number of particles is conserved at all times, we have $\sum_{k=0}^{S-1}x_{k}(t)=1$ for all $t\geq 0$.

In order to explore how $\Delta$ affects the evolution of the system we show in Fig. 2 snapshots of the distribution of the fractions $\{x\}$ at moments $t=1,600$ and $5000$ for three distinct realizations with error amplitudes $\Delta=0,1$ and $5$, for $S=100$ states and $N=10^{3}$ particles. At the early time $t=1$, $\{x\}$ looks nearly uniform in all cases, but then evolves towards a distribution that depends on $\Delta$. In the noiseless case $\Delta=0$ (left column) the system reaches a final delta distribution corresponding to a configuration where all particles are in the same state $\theta_{k=66}=4.1448$ (bottom-left panel). This is a frozen configuration where particles’ states cannot longer evolve, and corresponds to one of the $S=100$ possible absorbing states of the MSVM Starnini et al. (2012); Pickering and Lim (2016). Instead, for $\Delta=1$ (center column) the distribution $\{x\}$ becomes narrower with time and seems to adopt a bell shape for long times, while for $\Delta=5$ (right column) $\{x\}$ looks quite uniform for any time. In the bottom row ($t=5000$) we observe that the width of $\{x\}$ increases with $\Delta$. Therefore, we can see that the imperfection in the copying process is playing the role of an external noise that allows the system to escape from an absorbing configuration.

To explore the effects of varying the number of particles $N$, we show in Fig. 3 the distribution $\{x\}$ at different times for $S=100$, $\Delta=1$, and system sizes $N=10^{2},10^{3}$ and $10^{4}$. We observe that for $N=10^{2}$ (left panels) $\{x\}$ is narrow at long times, but it becomes wider as $N$ increases (see bottom row for $t=5000$), and already looks quite uniform for large $N=10^{4}$.

In summary, the dynamics of the model can be roughly seen as a competition between two processes: the perfect copying of the voter dynamics that tries to bring all particles together around a single state, and the imperfect copying in the form of noise that spreads particles apart. When $\Delta$ and $N$ are small, the system reaches a global state of order where most particles have similar angles and thus the angles’ distribution is narrow, while increasing $\Delta$ and $N$ results in a wider distribution. One may wonder how this quasi-ordered state observed for small $\Delta$ is quantitatively affected by the system size, that is, whether it reaches a stationary value as $N$ increases. In order to investigate these issues we focus our analysis on two complementary magnitudes that characterize the system at the macroscopic level. These are the order parameter

[TABLE]

and the mean-squared deviation of the angular states

[TABLE]

Here $\left|\bullet\right|$ is the absolute value, while $\Theta_{m}(t)$ is the state of particle $m$ ($m=1,..,N$) at time $t$. The parameter $\psi$ ($0\leq\psi\leq 1$) is similar to that introduced in the context of flocking dynamics to quantify the degree of global alignment in a system of moving particles Vicsek et al. (1995); Baglietto and Vazquez (2018), while the parameter $\sigma_{\theta}$ is a measure of the width of the distribution of angular states. When all particles move in the same direction ($\theta_{m}=\theta~{}\forall m$), one can check that $\psi=1$ and $\sigma_{\theta}=0$, which corresponds to a totally ordered state. On the other extreme, when each particle moves in a random direction the distribution of angular states becomes uniform in the $[0,2\pi)$ interval, and thus $x_{k}=1/S$ for $k=1,..,S-1$. Then, defining $r\equiv e^{i2\pi/S}$ and writing $e^{i\theta_{k}}=e^{i2\pi k/S}=r^{k}$ the order parameter is $\psi=\left|\frac{1}{S}\sum_{k=0}^{S-1}r^{k}\right|^{2}=\left|\frac{1-r^{S}}{S(1-r)}\right|^{2}=0$, i e., the system is completely disordered. On its part, the mean-squared deviation takes the value

[TABLE]

where we have used the identities Eqs. (56) and (57) in Appendix A, with $M=S-1$, to perform the summations.

In Fig. 4 we plot the time evolution of the average value of $\psi$ and $\sigma_{\theta}^{2}$ over $10^{4}$ independent realizations, denoted by $\langle\psi\rangle$ and $\langle\sigma_{\theta}^{2}\rangle$, for the same set of parameters used in Figs. 2 and 3. Both magnitudes reach a stationary value that quantifies the level of order at the stationary state that corresponds to the distributions $\{x\}$ of Figs. 2 and 3 at long times (down row). We observe that the stationary value of $\langle\psi\rangle$ decreases monotonically with $\Delta$ [Fig. 4(a)] and $N$ [Fig. 4(b)], while the stationary value of $\langle\sigma_{\theta}^{2}\rangle$ increases with $\Delta$ [Fig. 4(c)] and $N$ [Fig. 4(d)], and appears to saturate at the value of the uniform distribution $\sigma_{u}^{2}$ from Eq. (4).

These results suggest that the system reaches complete order ($\psi=1$ and $\sigma_{\theta}=0$) only for the noiseless case $\Delta=0$ [Figs. 4(a) and 4(c)], and that for any given error amplitude $\Delta>0$ the order constantly decreases with the system size $N$ and eventually vanishes in the thermodynamic limit [Figs. 4(b) and 4(d)]. This would imply that a tiny amount of error in the copying dynamics is enough to lead to complete disorder ($\psi=0$ and $\sigma_{\theta}=\sigma_{u}$) in the $N\to\infty$ limit. In order to analyze in more detail these conclusions obtained from numerical evidence, we develop in sections IV and V two analytical approaches that allow to obtain expressions for the asymptotic behavior of $\langle\psi\rangle$ and $\langle\sigma_{\theta}^{2}\rangle$ in two different limits. The approach in section IV is based on the diffusion approximation given by the Fokker-Planck equation and provides accurate results in the large $N$ limit, while the continuum approach developed in section V implements a superposition principle with open boundary conditions that works well in the limits of large $S$ and small noise.

IV The Fokker-Planck approach for large $N$

In this section we develop an analytical approach to estimate the scaling of $\psi$ and $\sigma_{\theta}^{2}$ with $S$, $N$ and $\Delta$ in the large $N$ limit. We derive a Fokker-Planck equation for the distribution of particles’ states that allows to obtain the behavior of the average value of $\psi$ and $\sigma_{\theta}^{2}$. We show that, for any $\Delta>0$ and in the $N\to\infty$ limit, the average value of the order parameter vanishes as $\langle\psi\rangle\sim 1/N$, and that the mean-squared deviation approaches the uniform value as $\sigma_{u}^{2}-\langle\sigma_{\theta}^{2}\rangle\sim 1/N$. This result implies that even the smallest error in the copying process is enough to lead the system to complete disorder in the thermodynamic limit. For the sake of simplicity, we focus on the simplest non-trivial case $\Delta=1$ where, in an iteration step, a randomly chosen particle tries to copy the state $\theta_{j}$ of another random particle, adopting either state $\theta_{j-1}$, $\theta_{j}$ or $\theta_{j+1}$ with equal probability $1/3$. We then use some heuristic arguments to extend these results to the general case $\Delta>1$.

Even though we are aware that there might be many different ways to address this problem analytically, we follow here a physics approach based on the diffusion approximation that gives the Fokker-Planck equation. This approach is particularly useful in the $N>>1$ limit because it allows to obtain rather accurate expressions for the stationary second moments $\langle x_{i}x_{j}\rangle$ that appear in the average values of both $\psi$ and $\sigma_{\theta}^{2}$ when we expand Eqs. (2) and (3a), respectively.

We start by describing the state of the system by the set of variables

$\{x\}=\{x_{0},x_{1},..,x_{S-1}\}$, where $x_{k}$ ($k=0,..,S-1$) is the fraction of particles with state $\theta_{k}$ (at site $k$) subject to the constraint $\sum_{k=0}^{S-1}x_{k}=1$ for all times. When a particle makes a transition from site $k$ to site $j\neq k$, the state of the system changes from $\{x\}$ to a new state denoted by $\{x^{\prime}\}_{k\;j}^{-+}\equiv\{x_{0},..,x_{k}-1/N,..,x_{j}+1/N,..,x_{S-1}\}$ in which, compared to $\{x\}$, site $k$ has lost a particle ($x_{k}\to x_{k}-1/N$) and site $j$ has gained a particle ($x_{j}\to x_{j}+1/N$). This is indicated in the notation $\{x^{\prime}\}_{k\;j}^{-+}$ with the $-$ and $+$ signs on top of subindices $k$ and $j$, respectively. The probability $P(\{x\},t)$ that the system is in state $\{x\}$ at time $t$ obeys the master equation

[TABLE]

The transition rate $W\left(\{x\}\to\{x^{\prime}\}_{k\;j}^{-+}\right)$ is the probability per time step $\delta t=1/N$ that a particle jumps from site $k$ to site $j$, calculated as $\frac{N}{3}x_{k}(x_{j-1}+x_{j}+x_{j+1})$. That is, a particle in site $k$ is chosen with probability $x_{k}$, then it jumps to either sites $j-1$, $j$ or $j+1$ with probability $x_{j-1}+x_{j}+x_{j+1}$, and from there jumps to site $j$ with probability $1/3$. According to the periodic character of angles, we make $x_{-1}=x_{S-1}$ and $x_{S}=x_{0}$ for the $j=0$ and $j=S-1$ cases, respectively. Similarly, we can calculate the transition rate $W\left(\{x^{\prime}\}_{k\;j}^{-+}\right)\to\{x\})$ that corresponds to a particle that jumps from site $j$ to site $k$. Then, the transition rates are given by the expressions:

[TABLE]

and

[TABLE]

For convenience, we have simplified notation using the rising and lowering operators $W_{k\,j}^{\downarrow\uparrow}$ and $W_{k\,j}^{\uparrow\downarrow}$. For instance, the down (up) arrow on top of $k$ ($j$) indicates that the operator applied on $\{x\}$ decreases (increases) $x_{k}$ ($x_{j}$) in $1/N$. The transitions for the case $S=2$ were displayed separately because they take a different form. In this particular case there are only two angular states, $\theta_{0}=0$ and $\theta_{1}=\pi$, and thus the noise step moves a particle to any of the two angles with equal probability $1/2$, instead of probability $1/3$ as explained above for any $S\geq 3$.

We can now obtain the Fokker-Planck equation by Taylor expanding the first term of Eq. (5) up to second order in $1/N\ll 1$ for large $N$:

[TABLE]

where $W_{k\,j}^{\uparrow\downarrow}$ and $P$ are short notations for $W_{k\,j}^{\uparrow\downarrow}(\{x\})$ and $P(\{x\},t)$, respectively, which are the functions $W$ and $P$ applied to the unperturbed state $\{x\}$. Inserting expression Eq. (8) into Eq. (5) leads to

[TABLE]

where $\partial_{k}\equiv\partial/\partial x_{k}$ and $\partial_{kj}^{2}\equiv\partial^{2}/\partial x_{k}\partial x_{j}$. To arrive to Eq. (9) we have made two considerations. First, we have used the following equalities to simplify the summations:

[TABLE]

Second, we have used the constraint $\sum_{k=0}^{S-1}x_{k}=1$ to write $x_{S-1}$ in terms of the other fractions, $x_{S-1}=1-\sum_{k=0}^{S-2}x_{k}$, reducing the number of independent variables to $S-1$. This makes partial derivatives $\partial/\partial x_{S-1}$ vanish, and set to $S-2$ the upper limit of the summation over $k$.

Plugging expressions from Eq. (6) and Eq. (7) for the transition rates into Eq. (9), and performing the summations inside the brackets we arrive to the Fokker-Planck equation in its final form

[TABLE]

where

[TABLE]

and

[TABLE]

Equations (10), (IV) and (12) give the time evolution of the probability distribution of angular states in a population of $N$ particles. The stationary solution of this Fokker-Planck equation, denoted by $P_{st}(\{x\})$, can be used to obtain the average value of $\psi$ and $\sigma_{\theta}^{2}$ at the stationary state, as we do in the following subsections.

IV.1 Analysis of the $S=2$ case

In order to gain an analytical insight into the behavior of the system at the stationary state we start by studying the simplest case of two angular states $S=2$ ($\theta=0,\pi$). We notice that this $2$-state model corresponds to a particular case of a surface-reaction model with noise studied in Considine et al. (1989) where, in a single step of the dynamics, one randomly chosen particle takes either state $1$ of $-1$ with probability $p_{d}/2$, or copies the state of a random neighbor with the complementary probability $1-p_{d}$. When $p_{d}=1$ the surface-reaction model turns equivalent to our model for $S=2$. Equation (12) describes the time evolution of the probability of finding a fraction $x_{0}$ of particles with angle $\theta_{0}=0$, whose stationary solution with boundary conditions $P_{st}(0)=P_{st}(1)$ and $\left.\frac{\partial P_{st}}{\partial x_{0}}\right|_{x_{0}=0}=\left.-\frac{\partial P_{st}}{\partial x_{0}}\right|_{x_{0}=1}$ is

[TABLE]

which satisfies the normalization condition $\int_{0}^{1}P_{st}(x_{0})\,dx_{0}=1$. One can check that expression Eq. (13) corresponds to the $p_{d}\to 1$ limit of the stationary solution found in Considine et al. (1989). The reason why we assumed these particular boundary conditions for $P_{st}$ is because both states $\theta_{0}=0$ and $\theta_{1}=\pi$ are equivalent, and thus we expect $P_{st}$ to be symmetric around $x_{0}=1/2$. We see that the stationary distribution of the fraction of particles with angle $\theta_{0}=0$ given by Eq. (13) is a Gaussian centered at $x_{0}=1/2$, whose width decreases as $N^{-1/2}$ with the number of particles. The order parameter from Eq. (2) becomes

[TABLE]

Then, the average value of $\psi$ at the stationary state can be calculated using $P_{st}(x_{0})$ from Eq. (13) as

[TABLE]

where we have made the change of variables $y=2x_{0}-1$ and integrated by parts. To first order in $1/N$, Eq. (15) is reduced to the simple expression

[TABLE]

which shows that $\langle\psi\rangle$ vanishes in the $N\to\infty$ limit. On its part, the mean-squared deviation from Eqs. (3) is

[TABLE]

and thus its stationary average value is calculated as

[TABLE]

For $N\gg 1$, Eq. (18) is reduced to the simple expression

[TABLE]

which shows that $\langle\sigma_{\theta}^{2}\rangle$ approaches the mean-squared deviation $\sigma_{u}^{2}=\pi^{2}/4$ of the uniform distribution as $N\to\infty$. Equations (16) and (19) describe the main result in the analysis of ordering in the noisy MSVM, that is, the distribution of angular states becomes uniform in the $N\to\infty$, and thus the system achieves total disorder ($\psi=0$). Even though this applies here only for the two–angle case $S=2$, we shall see in the next subsection that the same scalings with $N$ are obtained for any $S\geq 3$ as well.

To interpret these results from the dynamics of the system we resort to Eq. (13) and observe that, at the stationary state of a single realization, $x_{0}$ and $x_{1}$ fluctuate around the value $1/2$ subject to the constraint $x_{0}(t)+x_{1}(t)=1$. When $N$ increases, the amplitude of fluctuations vanishes as $N^{-1/2}$, and thus $P_{st}$ tends to the delta function $P_{st}(x_{0})=\delta(x_{0}-1/2)$. Therefore, we obtain the expected results $\psi=0$ and $\sigma_{\theta}^{2}=\sigma_{u}^{2}=\pi^{2}/4$ from Eqs. (14) and (17), respectively.

IV.2 Analysis of the general case $S\geq 2$

In the last section we obtained expressions for $\langle\psi\rangle$ and $\langle\sigma_{\theta}^{2}\rangle$ when $S=2$ from the stationary solution Eq. (13) of the Fokker-Planck equation. However, it seems hard to integrate analytically Eq. (10) and find an expression for the stationary distribution $P_{st}(\{x\})$ for the general case $S\geq 3$. Nevertheless, we shall see in the next two subsections that $\langle\psi\rangle$ and $\langle\sigma_{\theta}^{2}\rangle$ can be estimated by expressing them in terms of the second moments of $P_{st}(\{x\})$ which, in the limit of large $N$, can be obtained without knowing the explicit functional form of $P_{st}(\{x\})$.

IV.2.1 Calculation of the order parameter $\langle\psi\rangle$

We start by using the equality $e^{i\theta_{k}}=\cos\theta_{k}+i\sin\theta_{k}$ and rewriting the order parameter from Eq. (2) as

[TABLE]

Expanding the two squared terms of Eq. (20) leads to

[TABLE]

where we have used the formula for the cosine of the sum of two angles. Now, the average value of $\psi$ at the stationary state is

[TABLE]

Here we have exploited the translational symmetry of the angle’s space (a chain ring). We assumed that the second moments are invariant under translation, and thus they are a function $z_{j-k}$ of the distance $j-k\geq 0$ between $j$ and $k$, i e., $\langle x_{k}x_{j}\rangle=\langle x_{k-1}x_{j-1}\rangle=\langle x_{k-2}x_{j-2}\rangle=...=\langle x_{0}x_{j-k}\rangle\equiv z_{j-k}$. In particular, we have $\langle x_{k}^{2}\rangle=\langle x_{0}^{2}\rangle\equiv z_{0}$ for all $k=0,..S-1$. Then, expressing the double sum of the second term of Eq. (21) as a single sum over the index $n\equiv j-k\geq 0$, we arrive to

[TABLE]

To perform the summation in Eq. (22) we need to calculate the stationary value of the moments $z_{n}=\langle x_{k}x_{j}\rangle$. Given that we do not know how to obtain $P_{st}$ for $S\geq 3$ we use a different approach, developed in Appendix B. That is, starting from the Fokker-Planck equation we derive a system of coupled difference equations that relate the moments $z_{n}$ at the stationary state, whose solution gives the following approximate expressions for $z_{n}$ in the $N\gg S^{2}>1$ limit (see Appendix B for calculation details):

[TABLE]

Inserting expressions Eqs. (23a) for $z_{0}$ and $z_{1}$ into Eq. (22) we obtain for $S=2$

[TABLE]

which agrees with the expression Eq. (16) obtained in section IV.1 by direct calculation in the $N\gg 1$ limit. Now, for $S\geq 3$ we plug expression Eq. (23b) for $z_{n}$ into Eq. (22) and find, after doing some algebra,

[TABLE]

The summations above can be calculated exactly using complex variables (see Appendix C for details), obtaining

[TABLE]

Finally, inserting these expressions for the coefficients $a_{S}$ and $b_{S}$ into Eq. (24) we arrive to the following approximate expression for the order parameter

[TABLE]

This result shows that for any $S\geq 3$ the order parameter vanishes as $1/N$ in the $N\to\infty$ limit.

In Fig. 5 we compare the behavior of $\langle\psi\rangle$ from Eq. (29) (solid lines) with that obtained from Monte Carlo (MC) simulations for $\Delta=1$ and different values of $S$ and $N$ (symbols). We observe that, for a given $S$, the agreement between the analytical curve and the numerical data becomes better as $N$ increases, and it is very good when $N\gtrsim S^{2}$. This suggests plotting the data as a function of the rescaled variable $u\equiv N/S^{2}$, as we show in the inset of Fig. 5. The straight solid line is the analytical approximation

[TABLE]

obtained by expanding Eq. (29) to first order in $1/S$. We see that all data collapses into a single curve that follows the power law decay Eq. (30) when $u$ is approximately larger than $1$, i e., for $N\gtrsim S^{2}$ as mentioned above.

Even though the above analysis was performed for the particular case in which the angle perturbation is to first nearest-neighbor only ($\Delta=1$), we shall see that similar scalings hold for $\Delta>1$. In the general case $\Delta\geq 1$, it proves useful to consider each iteration of the dynamics as the two-step process (copy $+$ noise) described in section II, where the noise is represented by a uniform random variable $\xi$ that takes a value in the discrete set $-2\pi\Delta/S,-2\pi(\Delta-1)/S,..,0,..,2\pi\Delta/S$ with the same probability $1/(1+2\Delta)$. In order to generalize expression Eq. (29) for $\Delta>1$ we shall assume that $\langle\psi\rangle$ is a function of the noise variance $\sigma_{\xi}$, calculated as

[TABLE]

where we have defined

[TABLE]

Note that by letting $\Delta$ and $S$ go to infinity while keeping the ratio $\Delta/S$ fixed, $\eta$ reduces to the simple expression $\eta=2\Delta/S$ that is the noise amplitude in the case of continuum angles $\theta\in[0,2\pi)$ when $S\to\infty$ (we shall exploit this observation in section V). For $\Delta=1$, we can write $S$ in terms of $\eta$ from Eq. (31) as $S=2\sqrt{2}/\eta$. Then, replacing this expression for $S$ into Eq. (29) we obtain

[TABLE]

with $\eta=\eta(\Delta,S)$ given by Eq. (31).

To test the validity of the relation between $\langle\psi\rangle$ and $\eta$ from Eq. (32) for any $\Delta$ and $S$ we have performed simulations for $\Delta=1,2,5,10,20$, seven values of $S$ for each $\Delta$, and various system sizes. Results are shown in Fig. 6(a) (symbols) where we plot $\langle\psi\rangle$ vs $\eta$, with $\eta$ given by Eq. (31) for each pair $(\Delta,S)$. Note that, for fixed values of $N$ and $\Delta$, changing $S$ implies varying $\eta$ along the $x$-axis. Each of the solid curves corresponds to the analytical prediction Eq. (32) for a fixed value of $N$ and varying $\eta$ continuously in the range $[0.01,1.44]$. We observe that expression Eq. (32) is a good estimation of the numerical value of $\langle\psi\rangle$ within a range of $\eta$ that increases with $N$. This shows that for any $\Delta\geq 1$, $S\geq 2\Delta$ and $N$ large, $\langle\psi\rangle$ can be expressed as a function of the parameter $\eta(\Delta,S)$ given by Eq. (31). Analogously to the $\Delta=1$ case, we can expand Eq. (32) to first order in $\eta$ and obtain

[TABLE]

which suggests the scaling $\langle\psi\rangle=f(\eta^{2}N)$. Indeed, we can see in Fig. 6(b) a good data collapse when the data is plotted as a function of $\eta^{2}N$ (for $\eta<0.6$), and that obeys the power law decay from Eq. (33) (solid line) when $\eta^{2}N\gtrsim 1$.

IV.2.2 Calculation of the mean-squared deviation $\langle\sigma_{\theta}^{2}\rangle$

From the definition Eqs. (3), the mean-squared deviation can be expressed as

[TABLE]

Then, the average value $\sigma_{\theta}^{2}$ at the stationary state is

[TABLE]

where we have replaced $\langle x_{k}\rangle$ by $1/S$ [see Eq. (87) in Appendix B for calculation details], $\langle x_{k}x_{j}\rangle$ by $z_{j-k}$ and expressed the summation $\sum_{k=0}^{S-1}k^{2}$ in terms of $S$ using Eq. (57). The double summation in Eq. (34) can we rewritten in terms of the index $n=j-k$ as

[TABLE]

where we have used identities (56) and (57). Therefore, we obtain

[TABLE]

Using the approximate expression for $z_{0}$ from Eq. (23a) we obtain, for the $S=2$ case,

[TABLE]

which agrees with Eq. (19) obtained by direct integration (see section IV.1). For the case $S\geq 3$, we replace the moments $z_{n}$ in Eq. (35) by the approximate expressions from Eq. (23b) and arrive to

[TABLE]

with

[TABLE]

To calculate the coefficients $c_{S}$ and $d_{S}$ above we expand the terms of each summation in powers of $n$ and use the identities Eqs. (56-60) to obtain, after some algebra,

[TABLE]

Replacing the above expressions for $c_{S}$ and $d_{S}$ in Eq. (36) and simplifying the resulting expression we finally arrive to

[TABLE]

Equation (37) tells that the average width of the angular states distribution is smaller than that of the uniform distribution $\sigma_{u}$, and that approaches $\sigma_{u}^{2}=\pi^{2}(S^{2}-1)/(3S^{2})$ as $1/N$ when $N$ increases. As previously suggested, this result shows that the distribution of angular states becomes uniform in the $N\to\infty$ limit, where disorder is total ($\psi=0$).

Figure 7 shows MC simulation results for the behavior of $\langle\sigma_{\theta}^{2}\rangle$ with $N$, for $\Delta=1$ and various values of $S$ (symbols). We see that the agreement with Eq. (37) (solid lines) is good for $N\gtrsim S^{2}$, as it happens with $\langle\psi\rangle$. This approximate lower limit for the validity of Eq. (37) can be better checked in the data collapse shown in the inset.

To generalize Eq. (37) for any noise amplitude $\eta$, we follow an analysis similar to the one done in section IV.2.1 for $\langle\psi\rangle$. Replacing the expression $S=2\sqrt{2}/\eta$ obtained for $\Delta=1$ into Eq. (37) we obtain

[TABLE]

with $\eta=\eta(\Delta,S)$ given by Eq. (31). To test Eq. (38) we performed MC simulations for various system sizes and various different values of the set $(\Delta,S)$, and calculated the average mean-squared deviation. Results are shown by symbols in Fig. 8. In panel (a) we plot $\langle\sigma_{\theta}^{2}\rangle$ as a function of $\eta(\Delta,S)$. We see that $\langle\sigma_{\theta}^{2}\rangle$ increases with $\eta$ and saturates at the value $\sigma_{u}^{2}\simeq\pi^{2}/3$ corresponding to the large number of states ($S>10$) used in simulations, and then decreases for larger $\eta$. We also observe that, for a given $N$, Eq. (38) (solid lines) gives a good estimation of the numerical data for the largest values of $\eta$. The collapse of the data in panel (b) shows that $\langle\sigma_{\theta}^{2}\rangle$ is a function of $\eta^{2}N$, as long as $\eta$ is small enough. This behavior is in agreement with Eq. (38), from where we see that for small $\eta$ is $\langle\sigma_{\theta}^{2}\rangle\simeq\frac{\pi^{2}}{3}\left(1-\frac{2}{5\,\eta^{2}N}\right)$. This theoretical approximation is plotted by a solid line in the inset, showing the power-law approach $2/(5\,\eta^{2}N)$ of $\langle\sigma_{\theta}^{2}\rangle$ to $\pi^{2}/3$ that is valid when $\eta^{2}N\gtrsim 1$. For the sake of clarity, data points for large values of $\eta$ that fall off the straight line were removed.

We can also see in the main plot of Fig. 8(b) that for $\eta^{2}N\lesssim 10^{-2}$ the data seems to follow a power-law increase with an exponent similar to $1/2$, suggesting the scaling $\langle\sigma_{\theta}^{2}\rangle\sim\eta N^{1/2}$ (solid line). In the next section we give a theoretical insight into this particular behavior of the system in the low noise limit.

V The continuum approach for large $S$

In this section we consider the limiting case $S\gg\Delta\gg 1$, that is, the limit of a large number of particles’ states and a very small amplitude of the relative error, $\eta\simeq 2\Delta/S\ll 1$ [see Eq. (31)]. Before entering into the definition of the dynamics, we describe bellow a series of assumptions that we make to simplify our analysis.

In the limit of very large $S$ the angular space becomes continuous, and thus the state of a given particle $i$ at time $t$ can be taken as a real variable $\Theta_{i}(t)$ in the continuous space. The noise $\xi_{i}(t)$ introduced in the state of particle $i$ after the copying process at time $t$ can also be considered as a continuous variable uniformly distributed in the interval $[-\eta\pi,\eta\pi]$ (uniform white noise), with first moment $\langle\xi_{i}(t)\rangle=0$ and second moment $\langle\xi_{i}^{2}(t)\rangle=\pi^{2}\eta^{2}/3$ for all $i$ and $t$. Moreover, when the noise amplitude $\eta$ is very small we expect that the distribution of particles’ states will be very narrow at the stationary state, as compared to the length $2\pi$ of the angular space (see for instance down-left panel of Fig. 3). Thus, all particles will be far from the borders at [math] and $2\pi$ during the time we consider in this analysis. Therefore, we assume that particles can freely diffuse in the entire real axis with no boundary conditions. Besides, we consider that the updates of particles’ states take place in parallel (all at the same time) at discrete integer times $t=0,1,2...$, which resembles the update of the Wright-Fisher model in contrast to the sequential update of the Moran or voter dynamics implemented in our model. As each particle interacts once in average per unit time, we expect that both the sequential and the parallel updates have similar mesoscopic behavior. In fact, it was shown in Blythe and McKane (2007) that the Wright-Fisher and Moran models give the same mesoscopic Fokker Planck equation.

Let us consider that at a given time $t-1$ the states of particles are described by the angles $\Theta_{i}(t-1)$, with $i=1,..,N$. Then, in a single step of the parallel dynamics, for each particle $i$ we select a random particle $j\neq i$ and update its state according to $\Theta_{i}(t-1)\to\Theta_{i}(t)=\Theta_{j}(t-1)+\xi_{i}(t)$, where $\xi_{i}(t)\in[-\eta\pi,\eta\pi]$ is a small perturbation in the form of noise. We remark that the state of every particle at time $t$ depends on the state of another particle at the previous time $t-1$, and that time is increased by $\Delta t=1$ after all particles have updated their states. Suppose that at time $t-1$ particle $i$ copies the state of particle $j$, who has copied the state of particle $k$ at the previous time $t-2$ with a perturbation $\xi_{j}(t-1)$. Then, the state of $i$ at time $t$ can be expressed as $\Theta_{i}(t)=\Theta_{k}(t-2)+\xi_{j}(t-1)+\xi_{i}(t)$ and, iterating back in time until $t=0$, as

[TABLE]

where $\delta_{i}^{0}(t)=\Theta_{m}(0)$ is the initial state of some particle $m$, $\delta_{i}^{\tau}(t)=\xi_{n}(\tau)$ (for $1\leq\tau\leq t-1$) corresponds to the noise added to some particle $n$ at time $\tau$, and $\delta_{i}^{t}(t)=\xi_{i}(t)$. The reason why we use the notation $\delta_{i}^{\tau}(t)$ is because each term in the summation of Eq. (39) depends on the present time $t$ and the past time $\tau$, as we shall see bellow. Given that the set of numbers $\delta_{i}^{\tau}(\tau)=\xi_{i}(\tau)$ introduced in the system at $t=\tau$ will play a special role in the rest of this section, we will refer to them with the name of generation $\tau$.

To help better understand the dynamics of the system we show in Fig. 9 a simple example for the evolution of a four-particle system in a single realization. Initially, particles have states $\Theta_{i}(0)$ ($i=1,2,3,4$). In the first time step, particle $1$ copies the state of particle $2$, $2$ copies the state of $3$, while $3$ and $4$ copy the state of $1$, and then a different noise $\xi_{i}(1)$ is added to each particle. Thus, particles’ states at $t=1$ can be written as

[TABLE]

In terms of generations and by looking at the right-hand side of these equations, we can say that at time $t=1$ only $3$ different states are alive in generation $\tau=0$ [$\Theta_{1}(0)$, $\Theta_{2}(0)$ and $\Theta_{3}(0)$], while generation $\tau=1$ has just been born with its $4$ different states [$\xi_{1}(1)$, $\xi_{2}(1)$, $\xi_{3}(1)$ and $\xi_{4}(1)$]. The terms in the summation Eq. (39) for particle $1$ are $\delta_{1}^{0}(1)=\Theta_{2}(0)$ and $\delta_{1}^{1}(1)=\xi_{1}(1)$, and similarly for the states of the other three particles. In the second time step $t=2$, $1$ and $3$ copy $2$, $2$ copies $4$, and $4$ copies $3$. Then, $\Theta_{1}(2)=\Theta_{2}(1)+\xi_{1}(2)=\Theta_{3}(0)+\xi_{2}(1)+\xi_{1}(2)$ and thus $\delta_{1}^{0}(2)=\Theta_{3}(0),\delta_{1}^{1}(2)=\xi_{2}(1)$ and $\delta_{1}^{2}(2)=\xi_{1}(2)$. We note that the dependence in $t$ of each term $\delta_{1}^{\tau}(t)$ is due to the fact that when particle $1$ copies a particle $j$ it replaces all $t-1$ terms in the summation by those of $j$ at $t-1$. This process can be seen as particle $1$ copying the “state history” of particle $j$ from time [math] to $t-1$, formed by the list of $t-1$ terms in the summation $\Theta_{j}(t-1)=\sum_{\tau=0}^{t-1}\delta_{j}^{\tau}(t-1)$. In a similar way, we can iterate the copying processes at each time step as depicted in Fig. 9, and find the following states at $t=4$:

[TABLE]

It is interesting to note in Eqs. (40) that, even though all four particles started from different states $\Theta_{i}(0)$ at $t=0$ [$\delta_{i}^{0}(0)=\Theta_{i}(0)$], the state history of all particles at $t=4$ have the same “parent” particle $3$ (root) at $\tau=0$ [$\delta_{i}^{0}(4)=\Theta_{3}(0)~{}\forall i$]. That is, generation $\tau=0$ has already converged to $\Theta_{3}(0)$ at time $t=4$. They also have the same parent particle $2$ at $\tau=1$ [$\delta_{i}^{1}(4)=\xi_{2}(1)~{}\forall i$], i e., generation $\tau=1$ converged to $\xi_{2}(1)$ at $t=4$. We can also see that in Fig. 9, where the states of all particles at $t=4$ can be traced back in time to the state of particle $2$ at $t=1$, and to the state of particle $3$ at $t=0$. That is, when we look backwards in time, this “tree structure” that emerges can be seen as a system of coalescing random walks, which is the dual process of the voter model, as it is known in the mathematics literature on voter models Cox (1989). In general, if we look at a given time $\tau$, all four different numbers $\delta_{i}^{\tau}(\tau)=\xi_{i}(\tau)$ introduced at $t=\tau$ (generation $\tau$) change at each time step $t>\tau$ to the new values $\delta_{i}^{\tau}(t)$ following the copying dynamics between the numbers $\delta_{i}^{\tau}$ of the same generation $\tau$. Therefore, we can see $\delta_{i}^{\tau}$ as the state of the particle $i$ in generation $\tau$ that evolve under the rules of the original MSVM (without noise). In general terms, each generation $\tau$ behaves as a MSVM without noise that starts with $N$ different random states $\delta_{i}^{\tau}(\tau)=\xi_{i}(\tau)$ at time $\tau>0$, or from $\delta_{i}^{0}(0)=\Theta_{i}(0)$ at $\tau=0$. Given that convergence into a single state is eventually reached in the voter dynamics, we expect that the number of different states $s_{\tau}$ of a given generation $\tau$ decreases with time until consensus ($s_{\tau}=1$) is achieved when the generation has evolved for a time $t-\tau\sim N$, which happens in the example of Fig. 9 for $t=4$ and $\tau=0$ and $\tau=1$. These observations will be relevant for the calculations ahead. Contrary to the order of section IV.2, we first derive a scaling behavior for $\sigma_{\theta}^{2}$ that we use then to estimate $\psi$.

V.1 Calculation of the mean-squared deviation $\langle\sigma_{\theta}^{2}\rangle$

Now that we have described basic properties of the dynamics, let us start the analysis of the system by calculating the mean-squared deviation (sample variance) $\sigma_{\theta}^{2}(t)=\frac{1}{N}\sum_{i=1}^{N}\left[\Theta_{i}(t)-\overline{\theta}(t)\right]^{2}$, which can be expressed as

[TABLE]

where we have used Eq. (39) and introduced the new variable

[TABLE]

Then, the expectation value of $\sigma_{\theta}^{2}(t)$ over different runs of the dynamics can be written from Eq. (41) as

[TABLE]

Given that the perturbations introduced in a given generation $\tau$ are independent of the perturbations introduced in a different generation $\tau^{\prime}\neq\tau$, we have $\langle\mathcal{\delta}_{i}^{\tau}(t)\,\delta_{j}^{\tau^{\prime}}(t)\rangle=\langle\xi_{n}(\tau)\,\xi_{m}(\tau^{\prime})\rangle=\langle\xi_{n}(\tau)\rangle\langle\xi_{m}(\tau^{\prime})\rangle=0$ for all $i$ and $j$. Using this last relationship and Eq. (42) one can check that $\langle\Delta\delta_{i}^{\tau}(t)\,\Delta\delta_{i}^{\tau^{\prime}}(t)\rangle=0$ for $\tau^{\prime}\neq\tau$, and thus the second term of Eq. (43) vanishes, leading to the simple expression

[TABLE]

Here $\sigma_{\tau}^{2}(t)\equiv\frac{1}{N}\sum_{i=1}^{N}\left[\Delta\delta_{i}^{\tau}(t)\right]^{2}$ is the variance at time $t$ of the perturbations $\delta_{i}^{\tau}(\tau)$ introduced at time $\tau$ (we have used $\sigma_{\tau}=\sigma_{\delta^{\tau}}$ to simplify notation). We see that the expectation value of the dispersion in the particles’ states at time $t$ can be expressed as the superposition of the corresponding expectation values in the generations introduced in the previous times $0\leq\tau\leq t$, which have evolved under the multi-state voter dynamics during a time $t-\tau$.

We can now use the basic known results of the MSVM described in section II to study the evolution of the variance in a given generation $\tau$. Starting from $N$ states $\delta_{i}^{\tau}(\tau)$ at time $\tau$, the number of different states occupied by particles in generation $\tau$, $s_{\tau}$, decreases with time due to the copying process until all particles condensate into one state. Therefore, we expect that the dispersion will decrease with time and reach the value $\sigma_{\tau}^{2}=0$ when consensus in that generation is reached. From Eq. (1), the mean number of states in generation $\tau$ behaves approximately as

[TABLE]

where we have assumed an initial value $s_{\tau}(\tau)=S=N\gg 1$ and introduced the prefactor $2$ that appears in the parallel update (see Baglietto and Vazquez (2018)). Note also that $s_{\tau}$ is $1$ after consensus is reached in a time of order $N$. Then, if at time $t$ there are $s_{\tau}(t)\leq N$ surviving states that we call $\tilde{\delta}_{k}^{\tau}(t)$ [$k=1,..,s_{\tau}(t)$], we can express the states’ variance as $\sigma_{\tau}^{2}(t)\simeq\frac{1}{N}\sum_{k=1}^{s_{\tau}(t)}n_{k}(t)\left[\tilde{\delta}_{k}^{\tau}(t)-\overline{\tilde{\delta}^{\tau}}(t)\right]^{2}$, where $n_{k}$ is the number of particles with state $\tilde{\delta}_{k}^{\tau}$ and $\overline{\tilde{\delta}^{\tau}}=\frac{1}{s_{\tau}}\sum_{k=1}^{s_{\tau}}\tilde{\delta}_{k}^{\tau}$ is their mean value (note that $k$ is the state label rather than the particle label). Assuming that the $N$ particles are distributed uniformly among the $s_{\tau}(t)$ states, we can write $n_{k}(t)\simeq N/s_{\tau}(t)$ $\forall k$ and thus

[TABLE]

Expanding the right-hand side of Eq. (46) we find the products $\tilde{\delta}_{k}^{\tau}\,\tilde{\delta}_{k^{\prime}}^{\tau}$, whose expectation value is

[TABLE]

where now $\delta_{k,k^{\prime}}$ is the Kronecker delta. Applying brackets to both sides of Eq. (46) and using relations Eq. (47) we obtain

[TABLE]

which finally becomes

[TABLE]

after replacing expression Eq. (45) for $s_{\tau}$. Plugging Eq. (49) into Eq. (44) and performing the summation

[TABLE]

we finally arrive to

[TABLE]

where $c$ is a constant. It is interesting to note from Eq. (50) that for any time $t\geq N$, on average, only the $N$ generations that are less than a distance $N$ from $t$ contribute to $\langle\sigma_{\theta}^{2}(t)\rangle$, given that generations earlier than $t-N$ have reached consensus and thus they have zero variance. This makes $\langle\sigma_{\theta}^{2}\rangle$ reach a stationary value when $t\geq N$.

In order to test the previous result Eq. (51), we have run simulations of the model in continuous state space (states $\theta$ and noise $\xi$ are real numbers in $[-\infty,\infty]$ and $[-\eta,\eta]$, respectively) under the parallel update. Fig. 10(a) shows that, for various system sizes, $\sigma_{\theta}^{2}$ grows with $\eta$ as a power law with exponent $2$ (dashed lines). Also, the data collapse in the inset confirms the scaling given by Eq. (51), where the straight line is the function $c\,\eta^{2}N$, with $c=3.24$ corresponding to the best fit (solid line).

Eq. (51) predicts that the distribution of angular states at the stationary state has a width that increases linearly with $\eta^{2}$ and $N$ when the system is unbounded, given that particles’ states can freely spread on the real axis. However, if the system has periodic boundaries at $\theta=0$ and $\theta=2\pi$, particles are bounded in $[0,2\pi)$ and thus the width saturates to the value $\sigma_{u}\simeq\pi/\sqrt{3}$ when $N$ and $\eta$ increase, as explained in section IV.2.2. Besides that this saturation is obviously not captured by the open boundary approach developed in this section, neither the scaling $\langle\sigma_{\theta}^{2}\rangle\sim\eta N^{1/2}$ observed in Fig. 8(b) for small $\eta^{2}N$ (solid line) agrees with that of Eq. (51). The reason for that discrepancy lies in the fact that, for a single realization in the periodic system, all particles concentrate in a narrow interval whose distribution has a mean-squared deviation that scales as $\langle\sigma_{\theta}^{2}\rangle\sim\eta^{2}N$ as long as particles are far from the boundaries (not shown), but when particles reach a boundary the distribution splits into two separate sharp distributions peaked at [math] and $2\pi$, increasing its mean-squared deviation from the small value $\sigma_{\theta}^{2}\sim\eta^{2}N$ to a much larger value $\sigma_{\theta}^{2}\simeq 2\pi^{2}$. This last value of $\sigma_{\theta}^{2}$ is roughly estimated assuming that particles are concentrated around the boundaries and thus the distribution looks like two Dirac delta functions at [math] and $2\pi$. Then, when $\sigma_{\theta}^{2}$ is measured at a given time of a single realization, we estimate that the probability $\omega$ that the distribution is split into two parts is proportional to the distribution’s width. That is, with probability $\omega\propto\sigma_{\theta}/2\pi\propto\eta N^{1/2}$ is $\sigma_{\theta}^{2}\simeq 2\pi^{2}$ and with the complementary probability $1-\omega$ is $\sigma_{\theta}^{2}\propto\eta^{2}N$. This leads to an average mean-squared deviation over many runs that scales as

[TABLE]

Equation (52) gives the correct scaling $\langle\sigma_{\theta}^{2}\rangle\sim\eta N^{1/2}$ obtained in the simulations for the discrete system under the sequential update [Fig. 8(b)] when $\eta^{2}N\lesssim 10^{-2}$ (solid line with slope $1/2$).

V.2 Calculation of the order parameter $\langle\psi\rangle$

For the low noise case, we shall see bellow that there is a simple relationship between the order parameter $\psi(t)$ and the variance $\sigma_{\theta}^{2}(t)$ of the particle system, which allows to estimate $\langle\psi\rangle$ from the behavior of $\langle\sigma_{\theta}^{2}\rangle$ found in the last subsection.

As we showed before, when the noise is very low the states of all particles are within a very narrow angular window that, for the sake of simplicity, it is assumed to be centered at $\theta=0$ given that $\psi$ is invariant under angular translations. Indeed, one can check from the definition Eq. (2) of the order parameter that any translation $\Theta_{m}\to\Theta_{m}+\alpha$ ($m=1,..,N$) returns the same value of $\psi$. Therefore, we can approximate the exponential functions in Eq. (2) as $e^{i\Theta_{m}}\simeq 1+i\Theta_{m}-\Theta_{m}^{2}/2$ to second order in $|\Theta_{m}|\ll 1$, and write the order parameter as

[TABLE]

where $\overline{\theta}$ and $\overline{\theta^{2}}$ are the first and second moments of the states’ distribution defined in Eqs. (3b), and $\sigma_{\theta}^{2}=\overline{\theta^{2}}-\overline{\theta}^{2}$. From Eq. (53), the average value of the order parameter at the stationary state is related to the average mean-squared deviation of the population by the simple expression

[TABLE]

As expected, the order increases as the width of the distribution of particles decreases, and reaches its maximum value $\psi=1$ at consensus ($\sigma_{\theta}^{2}=0$). Finally, replacing Eq. (51) for $\langle\sigma_{\theta}^{2}\rangle$ into Eq. (54) we obtain

[TABLE]

Figure 10(b) shows the behavior of the order parameter $\langle\psi\rangle$ with the rescaled variable $\eta^{2}N$, obtained from simulations of the continuous model under the parallel update for various system sizes. We observe that for $\eta^{2}N\lesssim 10^{-2}$ the data is in good agreement with Eq. (55) (solid line) using the value $c=3.24$ obtained from the best fit of the $\langle\sigma_{\theta}^{2}\rangle$ vs $\eta^{2}N$ data shown in the inset of Fig. 10(a). In the inset of Fig. 6(b) we show that Eq. (55) with a constant $c_{s}=1.64$ (solid line) describes the behavior of $\langle\psi\rangle$ at low noise for the discrete sequential version of the model as well. The prefactor similar to $1/2$ between the constants $c_{s}$ and $c$ of the two versions of the model also appears in the expression $3/(\pi^{2}\eta^{2}N)$, as compared to Eq. (33), which reproduces very well the behavior of $\langle\psi\rangle$ for large $N$ of the continuous parallel model shown in the inset of Fig. 10(b). We do not know how to explain this prefactor.

Equation (55) shows that for a system of size $N$, total order $\langle\psi\rangle=1$ is eventually achieved as the noise vanishes. This result completes the picture of the behavior of the model in the two limits. That is, for fixed $\eta>0$, complete disorder is reached in the $N\to\infty$ limit, while for fixed $N>0$, complete order is achieved as $\eta\to 0$.

VI Summary and Conclusions

We studied the dynamics of a multi-state voter model in mean-field with a degree of error or imperfection in the copying process. Starting from a uniform distribution of particles over a discrete angular space, we investigated the dynamics of ordering and its stationary state. When the copying is perfect the number of different angular states occupied by particles decreases monotonically with time until eventually the system reaches an absorbing state of complete order where all particles share the same state. However, when we add a source of imperfection in the copying process that leaves the states of two interacting particles similar but not exactly equal (an imperfect copying) a new scenario appears. The system evolves towards a stationary state characterized by an ordering level $\psi$ that depends on the number of particles $N$, the number of possible particle states $S$ and the copying error amplitude $\Delta$. We analyzed two different limits. In the large $N$ limit we proved by means of a Fokker-Planck equation approach that the average order decreases with $N$ and the relative error amplitude $\eta=2\Delta/S$ as $\langle\psi\rangle\simeq 6/(\pi^{2}\eta^{2}N)$ for $0<\eta\ll 1$ and $\eta^{2}N\gtrsim 1$. Besides, when $N$ and $\eta$ increase, the average mean-squared deviation of particles’ states approaches the value $\sigma_{u}^{2}$ corresponding to the uniform distribution as $\langle\sigma_{\theta}^{2}\rangle\simeq\sigma_{u}^{2}[1-2/(5\,\eta^{2}N)]$. These results imply that for any degree of error $\eta>0$ the system gets completely disordered in the thermodynamic limit $N\to\infty$, where the distribution of particles over the angular space is perfectly uniform. In the large $S$ limit we developed an analytical approach that assumes a continuum angular space and showed that when $\eta\to 0$ the system approaches total order as $\langle\psi\rangle\simeq 1-1.64\,\eta^{2}N$, while $\langle\sigma_{\theta}^{2}\rangle$ vanishes as $\langle\sigma_{\theta}^{2}\rangle\sim\eta N^{1/2}$. This result also shows that complete order $\langle\psi\rangle=1$ is only achieved for perfect copying $\eta=0$.

As mentioned in section IV.1, the $2$-state case of our MSVM is equivalent to the noisy voter model studied in Considine et al. (1989). This suggests that it might be possible to map the MSVM with $S>2$ states to the $2$-state noisy voter model studied in Considine et al. (1989); Peralta et al. (2018a) by finding appropriate copying and noise rates that depend on $\eta$, and then use known results on these studied models to derive the scaling relations obtained in this article. We have become aware of a recent unpublished article Herrerías-Azcué and Galla (2019) that investigates multi-state noisy voter models where imitation and mutation (noise) events occur at respective rates $r_{ji}$ and $\epsilon_{ji}$ that may depend on states $j$ and $i$. It seems that our model corresponds to an homogeneous $r_{ji}$ and a very particular choice of $\epsilon_{ji}$ that depends on the fraction of particles $x_{i+1}$ and $x_{i-1}$. This particular case is not explored by the authors who rather focus their study on the multistability properties of the system considering rates that are independent on the particles’ fractions $\{x\}$. We also need to mention that the version of the MSVM with continuous states studied in section V is related to a family of processes with $N$ branching particles, initially proposed in Brunet and Derrida (1997) to investigate selection mechanisms in biological systems and later extended to a continuous time version known as “N-branching Brownian motions” recently explored in Maillard (2016); A. De Masi and Soprano-Loto (2017a, b). However, these models introduce a type of asymmetric copying process that gives rise to a traveling wave of particles that moves to the right, which is absent in our model due to the symmetry of interactions.

It is also worth mentioning some possible implications that the studied model could have on some related problems. Within the context of flocking dynamics, the appearance of complete disorder for $\eta>0$ in the thermodynamic limit suggests an order-disorder transition at zero noise $\eta=0$, something unseen in related Vicsek-type models where a transition occurs at a finite critical value $\eta_{c}>0$. Within the context of population genetics, the addition of imperfection in the process of gene replication would lead to a population characterized by a diversity of gene types that would increase with the population size. Given that the results in this article are of mean-field type (all to all interactions), it should be worthwhile to study the imperfect MSVM in two-dimensional systems to investigate the effects of spatial interactions on the ordering dynamics. Finally, it might be interesting to study the imperfect copying mechanism on a constrained version of the MSVM Vazquez et al. (2003); Lanchier and Scarlatos (2017), where interactions are only allowed between agents whose opinion distance is smaller than a fixed threshold. These are possible topics for further investigation.

Acknowledgments

We acknowledge financial support from CONICET (PIP 11220150100039CO) and (PIP 0443/2014). We also acknowledge support from Agencia Nacional de Promoción Científica y Tecnológica (PICT-2015-3628) and (PICT 2016 Nro 201-0215).

Appendix A Power sums

Below we write closed expressions for five different power sums that are useful in several calculations along the article.

[TABLE]

Appendix B Calculation of the moments $z_{n}$

In this section we calculate approximate expressions for the moments $z_{n}=\langle x_{k}x_{j}\rangle$ at the stationary state, for any $S\geq 2$ and $N\gg 1$. For that, we derive a set of coupled equations that relate the first and second moments using the Fokker-Planck equation derived in section IV. We now illustrate this procedure for the simplest case $S=2$. From Eq. (12), the time evolution of the first moment $\langle x_{0}\rangle$ obeys the equation

[TABLE]

and thus at the stationary state we have

[TABLE]

where $P_{st}(x_{0})$ is the stationary distribution given by Eq. (13). The above integrals can be exactly calculated using $P_{st}(x_{0})$ but, instead, we can perform the integrals by parts and assume that $P_{st}(x_{0})$ and its derivatives are zero at the boundaries

[TABLE]

One can check directly from Eq. (13) that these boundary conditions are satisfied in the $N\to\infty$ limit. The reason is that $P_{st}(x_{0})$ is a Gaussian of width $\frac{1}{2}N^{-1/2}$ centered at $x_{0}=1/2$. Therefore, when $N$ is very large $P_{st}(x_{0})$ quickly drops to zero outside the interval $[1/2-N^{-1/2},1/2+N^{-1/2}]$, and thus $P_{st}(x_{0})$ and $\partial P_{st}(x_{0})/\partial x_{0}$ are expected to be similar to zero at $x_{0}=0$ and $x_{0}=1$. Then, integrating by parts Eq. (62) we obtain

[TABLE]

Using the boundary conditions Eqs. (63–64) we see that only the second term of Eq. (65) is not zero, leading to the simple relation

[TABLE]

from where the first moments read

[TABLE]

Following the same approach for $d\langle x_{0}^{2}\rangle/dt$ we obtain the relation

[TABLE]

Then, combining Eqs. (67) and (68) we obtain the following second moments for the $S=2$ case:

[TABLE]

quoted in Eq. (23a) of the main text. We can now apply the same procedure to calculate the moments $z_{n}$ ($0\leq n\leq S-1$) for the general case scenario $S\geq 3$. Even though the calculations are analogous to the ones for the $S=2$ case described above, the generalization is not straight forward because new types of integrals appear due the existence of crossed derivatives for $S\geq 3$. Same as before, the idea is to write a differential equation for the time evolution of each of the moments, $\langle x_{k}\rangle$, $\langle x_{k}x_{j}\rangle$ and $\langle x_{k}^{2}\rangle$, using the Fokker-Planck equation (10). We illustrate here this procedure for the second moment $\langle x_{l}x_{m}\rangle$ ($l\neq m$), and leave for the interested reader the corresponding calculations for the other two moments. From the Fokker-Planck equation (10) for $S\geq 3$, the time evolution of $\langle x_{l}x_{m}\rangle$ obeys the equation

[TABLE]

where

[TABLE]

Then, at the stationary state we have

[TABLE]

where the integrals $\tilde{I}_{lmk}$, $\tilde{I}_{lmkk}$ and $\tilde{I}_{lmkj}$ have the same form as those from Eqs. (72), (73) and (74), respectively, but integrating over the stationary distribution $P_{st}(\{x\})$ instead of $P(\{x\},t)$. To calculate these integrals we are going to assume that $P_{st}(\{x\})$ and its first derivatives are zero at the boundaries

[TABLE]

This is because, in analogy to the $S=2$ case, we expect for $S\geq 3$ a bell-shaped $P_{st}(\{x\})$ peaked at the point ($x_{k}=1/S\,\,\forall k$) of the $(S-1)$–dimensional space $\{x_{0},..,x_{S-1}\}\in[0,1]^{S-1}$, and that the width of $P_{st}(\{x\})$ vanishes as $N^{-1/2}$ with the system size. Therefore, if $N^{-1/2}$ is much smaller than the distance $1/S$ that separates the location of the peak and the closest boundary $x_{k}=0$, then $P_{st}(\{x\})$ and its first derivatives should be similar to zero at both boundaries $x_{k}=0$ and $x_{k}=1$ for all $k$. This allows to give the rough estimation $N\gg S^{2}$ that relates the system size $N$ and the number of angular states $S$ for which the approximations we make in this section are valid. In the next three subsections we calculate the integrals $\tilde{I}_{lmk}$, $\tilde{I}_{lmkk}$ and $\tilde{I}_{lmkj}$.

B.1 Calculation of $\tilde{I}_{lmk}$

Case $k\neq l\neq m$:

[TABLE]

where we have used the simplified notation $P_{st}=P_{st}(\{x\})$, and the boundary condition Eqs. (76) to set to the term inside the curly brackets to zero.

Case $k=l$:

[TABLE]

Similarly, we can show that for $k=m$ is $\tilde{I}_{lmm}=-\langle x_{l}A_{m}\rangle$. Then, combining all cases we can write

[TABLE]

B.2 Calculation of $\tilde{I}_{lmkk}$

Case $k\neq l\neq m$:

[TABLE]

where we have used the boundary conditions Eqs. (76–77) to set to the term inside the brackets to zero.

Case $k=l$:

[TABLE]

Similarly, for $k=m$ we obtain $\tilde{I}_{lmmm}=0$. Then,

[TABLE]

B.3 Calculation of $\tilde{I}_{lmkj}$

Case $k\neq l\neq m\neq j$:

[TABLE]

where we have used the boundary conditions Eqs. (76–77) to set to the term inside the brackets to zero.

Case $k=l,j=m\neq k$:

[TABLE]

Following the same type of calculations we find $\tilde{I}_{lmlj}=0$ for $k=l$ and $\tilde{I}_{lmkm}=0$ for $j=m$. Then,

[TABLE]

Finally, plugging expressions (79), (80) and (81) for $\tilde{I}_{lmk}$, $\tilde{I}_{lmkk}$ and $\tilde{I}_{lmkj}$, respectively, into Eq. (75) we obtain the following equation for the moments:

[TABLE]

Now that we have calculated, starting from Eq. (71) for $d\langle x_{m}x_{l}\rangle/dt$, the first equation that relates the moments, it is possible to obtain two more equations by following the same procedure for $d\langle x_{l}\rangle/dt$ and for $d\langle x_{l}^{2}\rangle/dt$ (calculations not shown). This results in the following system of equations

[TABLE]

where indexes $l$ and $m$ were renamed as $k$ and $j$. Plugging expressions for $A_{k}$, $A_{j}$, $B_{kk}$ and $B_{kj}$ from Eqs. (IV) into Eqs. (83) leads to

[TABLE]

The solution to Eq. (84) with periodic boundary conditions $\langle x_{-1}\rangle=\langle x_{S-1}\rangle$ and $\langle x_{S}\rangle=\langle x_{0}\rangle$ that satisfies the constraint $\sum_{k=0}^{S-1}\langle x_{k}\rangle=1$ is

[TABLE]

If now we use the definition $z_{n}\equiv z_{j-k}\equiv\langle x_{k}x_{j}\rangle$ and $\langle x_{k}\rangle=1/S$, we can rewrite Eq. (85) in terms of $z_{0}$ and $z_{1}$, and Eq. (86) in terms of $z_{n-1}$, $z_{n}$ and $z_{n+1}$. For that, we need to take into account the identity $z_{-n}=\langle x_{j}x_{k}\rangle=\langle x_{k}x_{j}\rangle=z_{n}$. Then, after some algebra and regrouping terms we arrive to the following system of equations for the moments

[TABLE]

where we defined

[TABLE]

The system of Eqs. (88–89) must also satisfy the constraint

[TABLE]

which is derived from the normalization condition $\sum_{k=0}^{S-1}x_{k}=1$ by multiplying both sides of this equality by $x_{j}$, then taking the average $\langle\cdot\rangle$ at both sides, and setting $\langle x_{j}\rangle=1/S$. The number of independent equations in the system of Eqs. (88–89) can be reduced by half by implementing the periodic property

[TABLE]

obtained by direct calculation: $z_{S-n}=z_{n-S}=z_{(j-k)-S}=\langle x_{S}x_{j-k}\rangle=\langle x_{0}x_{j-k}\rangle=z_{j-k}=z_{n}$. We now show how the system of equations is reduced for the case of even $S$. Equation (89) for $n=S/2$ reads $z_{S/2+1}-rz_{S/2}+z_{S/2-1}$ which, after replacing $z_{S/2+1}$ by $z_{S/2-1}$ from the periodic relation Eq. (93) becomes $2z_{S/2-1}-rz_{S/2}=0$. In general, one can prove that the equation for an index $n$ in the interval $S/2+1\leq n\leq S-1$ becomes the same equation as that for index $S-n$. Therefore, the system of Eqs. (88–89) for $S>3$ even is

[TABLE]

where Eq. (97) comes from the constraint Eq. (92). The same analysis applied to $S\geq 3$ odd leads to the following system

[TABLE]

Even though the set of equations for $S$ odd looks different from that of $S$ even, the solutions turn out to be the same, and thus we focus now on $S$ even. The solution of Eqs. (94–97) for any $S$ and $N$ is rather complicated, but because we are interested in the limit of $N\gg 1$ we look for solutions of the form $z_{n}=C\left(1-\alpha_{n}\frac{1}{N}\right)+\mathcal{O}\left(1/N^{2}\right)$, where $C$ is a constant and $\alpha_{n}$ are functions that depend on $n$ and $S$. We note that this proposed ansatz agrees with the corresponding expressions found for the $S=2$ case [Eqs.(69) and (70)], with $C=1/4$, $\alpha_{0}=-1/N$ and $\alpha_{1}=1/N$. Then, to first order in $\epsilon\equiv 1/N\ll 1$ the solutions take the approximate form

[TABLE]

Inserting the above expressions for $z_{n}$ into Eqs. (94–97) and neglecting $\epsilon^{2}$ terms, we arrive to the following closed system of $S/2+2$ equations for $C$ and $\alpha_{n}$ with $S/2+2$ unknowns

[TABLE]

To solve the system of Eqs. (103–106) we define $\beta_{n}\equiv\alpha_{n}-\alpha_{n-1}$. Then, from Eqs. (103) and (105) we get $\beta_{1}=\frac{3}{2}\left(\frac{1}{SC}-1\right)$ and $\beta_{S/2}=3/2$, respectively, while Eq. (104) becomes $\beta_{n+1}=\beta_{n}-3$, whose solution is $\beta_{n}=\beta_{1}-3(n-1)$. This last equation for $n=S/2$ leads to a simple relation between $S$ and $C$, from where we obtain $C=1/S^{2}$ and thus $\beta_{1}=\frac{3}{2}(S-1)$ and $\beta_{n}=\frac{3}{2}(S+1-2n)$. Therefore, we get the relation $\alpha_{n}=\alpha_{n-1}+\frac{3}{2}(S+1-2n)$ which can be solved by simple iteration, leading to $\alpha_{n}=\alpha_{0}+\frac{3}{2}n(S-n)$. Using this last expression for $\alpha_{n}$ in Eq. (106) and setting $C=1/S^{2}$ we arrive to a closed equation for $\alpha_{0}$, with solution $\alpha_{0}=(1-S^{2})/4$. Thus, the final expression for $\alpha_{n}$ becomes

[TABLE]

Finally, using expression (107) for $\alpha_{n}$ and $C=1/S^{2}$ in Eq. (102) we obtain the following approximate expression for $z_{n}$ to first order in $1/N\ll 1$:

[TABLE]

quoted in Eq. (23b) of the main text. We can check that expression (108) is a solution of the system of Eqs. (98–101) for $S$ odd as well.

Appendix C Calculation of the coefficients $a_{S}$ and $b_{S}$

In this section we derive the expressions Eqs. (27) and (28) for the coefficients $a_{S}$ and $b_{S}$, respectively. We start by rewriting these coefficients given by the summations in Eqs. (25) and (26), as the real part of complex numbers $\mathcal{A}_{S}$ and $\mathcal{B}_{S}$, respectively

[TABLE]

where

[TABLE]

with $r\equiv e^{i2\pi/S}$. To perform the summations in Eqs. (109) and (110) we first extend the upper limit to $n=S$ and the lower limit of Eq. (110) to $n=0$, then expand the terms in brackets and define

[TABLE]

Thus, Eqs. (109) and (110) can be written as

[TABLE]

In order to find $y_{m}$ ($m=0,1,2,3$) we start from the well known geometric series

[TABLE]

and differentiate this formula with respect to $r$ to obtain

[TABLE]

These formulas can be greatly simplified by noting that $r^{S}=e^{i2\pi}=1$, and thus $r^{S+1}=r$ and $r^{S+2}=r^{2}$, which leads to

[TABLE]

Replacing these expressions for $y_{m}$ in Eqs. (111) and (112) we obtain

[TABLE]

Substituting in the above expressions $r$ by $e^{i2\pi/S}$ and using for convenience the identities $1-r=-2\,i\,e^{i\pi/S}\sin(\pi/S)$ and $1+r=2\,e^{i\pi/S}\cos(\pi/S)$ we finally arrive to

[TABLE]

whose real parts correspond to Eqs. (27) and (28), respectively, of the main text.