{'fraction_non_alphanumeric': 0.03598257685752162, 'fraction_numerical': 0.010542263745975633, 'mean_word_length': 4.3553752535496955, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ordering in voter models on networks: Exact reduction to a single-coordinate diffusion

We study the voter model and related random-copying processes on arbitrarily complex network structures. Through a representation of the dynamics as a particle reaction process, we show that a quantity measuring the degree of order in a finite system is, under certain conditions, exactly governed by a universal diffusion equation. Whenever this reduction occurs, the details of the network structure and random-copying process affect only a single parameter in the diffusion equation. The validity of the reduction can be established with considerably less information than one might expect: it suffices to know just two characteristic timescales within the dynamics of a single pair of reacting particles. We develop methods to identify these timescales, and apply them to deterministic and random network structures. We focus in particular on how the ordering time is affected by degree correlations, since such effects are hard to access by existing theoretical approaches.


Introduction
The voter model [1][2][3][4] is arguably the simplest nonequilibrium ordering process. Its dynamics are perhaps most easily described within its original formulation as a model for population dynamics [1]. Each site of a lattice (or network) is occupied by an individual belonging to one of a number of species. In each timestep, a site is chosen at random, and the individual at that site is replaced by an offspring (copy) of an individual selected at random from a neighbouring site. In a finite system, one of the possible absorbing states-in which only a single species remains-is eventually reached. From a physical perspective, the interest in this process is that order is arrived at purely through fluctuations, rather than driven by surface tension as in the Ising model (see e.g., [5,6]; also [7] for a recent review).
Many variants of the basic voter model have been investigated, some quite recently: examples include those where multiple individuals occupy a site [8], where the site chosen in each timestep contains the individual that is copied rather than replaced [9,10], or where different sites have different probabilities of being chosen for update [11]. Specific instances drawn from this wider class of random-copying models have applications in population genetics, ecology and cultural evolution, as reviewed in [12]. Despite their simplicity, these models reproduce observed patterns of species diversity in ecosystems [13] and of the cultural dynamics of features like baby names [14]. They also play an important role as null models of evolution in the absence of selection [15].
A stochastic equation of motion for the frequency ξ of a chosen species is straightforwardly obtained in the mean-field limit (i.e., when the system lacks spatial structure). More precisely, when each site is visited for update at a uniform rate m, and the individual landing on the target site is sampled uniformly from the entire population, the probability distribution of ξ is-for sufficiently large systems-governed by Derivations of this equation can be found in many places, e.g., [12,16,17]. These stochastic dynamics have three key properties. First, there is no deterministic (drift) term in the above equation, confirming that order is attained through fluctuations alone. Second, all timescales of the ordering process are proportional to the parameter N e : in the mean field, this parameter is simply equal to the number of sites, N. Finally, the full, time-dependent solution for P (ξ, t) from arbitrary initial condition is known [16]. Then, for example, the probability that complete order has been attained by time t is given by the combination P (0, t) + P (1, t).
Many recent works have switched focus to models defined on nontrivial network structures, notably heterogeneous networks in which the degree can vary widely between nodes [9,10,[18][19][20][21][22][23][24][25][26][27]. Such networks are particularly relevant to cultural applications, where the interaction frequency between agents need not be purely a function of spatial proximity. Remarkably, it has been found that even in the presence of significant heterogeneity, the ordering dynamics is well-described by the mean-field equation (1). For example, the mean time to reach complete order as predicted by (1) has been found to apply in models with complex spatial structure [10,[23][24][25][26][27]. The only differences are: (i) that ξ is a weighted average of a species frequency over the network; and (ii) the effective size, N e , need not equal the actual network size N (or even scale linearly with it). Thus, a good understanding of ordering behaviour within random-copying models on networks can be gained if one can identify when a reduction to (1) is possible, and if so, what the value of N e is. For example, the entire distribution of ordering times can be read off from the exact solution of (1).
In this work we will demonstrate that the mean-field diffusion equation (1) is in fact exact for general random-copying processes on arbitrary networks at all but the very earliest times when there is a separation of emergent timescales within a two-particle coalescence reaction dynamics that is dual to the random-copying process [28,29]. In several earlier works [9,10,[23][24][25][26][27], the equation (1) was typically obtained by making specific assumptions and approximations, and the fundamental physical reason why the simple diffusion (1) should apply was not clearly identified.
The central object in our analysis is the quasi-stationary distribution of a pair of particles that have survived for a long time without reacting. This state decays at some characteristic rate, and if its lifetime greatly exceeds the relaxation time to the quasi-stationary state from a given initial condition, the model exhibits the separation of timescales that is required for the validity of (1). The effective size N e is then proportional to the single, dominant timescale in the reaction dynamics, that is, the lifetime of the quasi-stationary state.
This paper is divided into two parts. In the first part (Sections 2-5) we define the general class of random-copying models that includes well-studied instances, like the voter model, as special cases, and derive our main results. These are: Eq. (1), the condition for its validity and the expression for the fundamental parameter N e . As we will need to switch repeatedly between the forward-time formulation of the dynamics in terms of a diffusion equation and an equivalent backward-time particle reaction process whose details are not widely discussed in the physics literature, we must unfortunately contend with lengthy, but necessary, preliminaries in Sections 2-4. Once this is done, the derivations of our main results follow in Section 5. We emphasise that the validity of Eq. (1) can be determined by considering only a single pair of coalescing random walkers on the network. In principle, ascertaining whether the system is completely ordered requires knowledge of the state of all N sites simultaneously: hence being able to determine whether the considerably simpler single-coordinate diffusion applies without having to solve the full dynamics is of great practical benefit.
In the second part of this paper, we go beyond purely formal statements about random-copying dynamics and focus on strategies to estimate the emergent timescales of the two-particle reaction process. These methods are developed in Section 6 and used to demonstrate the validity of the diffusion (1) for specific combinations of network structures and copying dynamics in Section 7. Although these methods are restricted to random-copying processes that map onto time-reversal symmetric random walks, it turns out that the most widely-studied variations on voter-model dynamics [9,10,20,26,27] fall into this category. On the other hand, we do not (in contrast to some existing treatments [9,10,23,[25][26][27]) invoke a coarse-graining over multiple sites of the network to obtain a diffusion equation. Thus, we can use our analytical methods to investigate the effects of degree correlations on the effective size N e that appears in Equation (1), which until now has been accessible only by computer simulation.
We remark that the idea that random-copying dynamics in the presence of spatial structure may be adequately represented by a mean-field model with some effective size N e has a very long history in population genetics, dating back to the work of Wright from the 1930s [30]. Furthermore, various separation of timescales arguments have been employed as a means to identify an effective size [31][32][33][34][35][36]. In this work, our main contribution to this line of enquiry is to highlight the two particular emergent timescales that are relevant in determining the applicability (or otherwise) of a single effective size: their importance does not appear to have been fully recognised before. We furthermore show that this separation exists even in cases where it has been argued that there will be more than one relevant timescale in the dynamics, and that a description in terms of a single effective size will not apply [35]. We discuss these points in much more detail in the conclusion, Section 8, along with a number of remaining open questions.

Definition of a class of spatially-structured random-copying models
We consider a network of N nodes (or sites) labelled by the integers i = 1, 2, 3, . . . , N. Site i hosts a population of M individuals, each of which is a member of one of a number of species. The dynamics is specified in terms of a set of copying rates µ ij . Each corresponds to an independent Poisson process: at rate µ ij one of the M individuals on site j is copied, and the offspring individual replaces a randomly-chosen individual on site i. See Figure 1. Note that the case i = j is legitimate, and indeed the rates µ ii will typically be nonzero in any particular instance of this model. We could also consider models where M varies from site to site within the framework presented below: however, this serves only to complicate certain definitions and expressions so we shall restrict the discussion to a spatially-uniform M.
The rates µ ij determine the model's network structure. Elements of the adjacency matrix A ij , which indicate if sites i and j are neighbours, are zero if µ ij = 0, and unity otherwise (except the diagonal elements A ii which are always taken to be zero). The degree z i of node i can then be defined as Although we no not require the adjacency matrix to be symmetric in the first part of this paper, we will for simplicity assume this is the case. (For example, we do not then have to distinguish between in-degree and out-degree). Note, however, that the symmetry A ij = A ji does not imply that µ ij = µ ji . In a finite system, there is some probability that a species eventually goes extinctthis occurs if all its representatives have been replaced by individuals of other species. Over time, the number of species decreases until only a single one remains. This stationary absorbing state of complete order is referred to variously as consensus or fixation (in the sociodynamics and genetics literature respectively).

Voter, invasion and link dynamics
The original voter model dynamics specifies that, in each timestep, a site is first chosen at random, and then a copy taken from one of its neighbours. This constrains the copying rates µ ij to have the form where µ is the rate at which a site is chosen for update. This is, however, only one of a vast range of possible ways to construct a set of copying rates. Also discussed in the literature is the invasion process, whereby the parent site j is chosen uniformly at random in each timestep, and the offspring 'invades' a randomly chosen neighbour (see, e.g., [9,10]). This corresponds to the choice Another possibility is link dynamics [9,10], in which an edge between two nodes is chosen at random, and the copy made in one of the two possible directions along that edge. One way to specify this process is through the rates wherez is the mean degree within the network. In networks where all nodes have the same degree, z i = z j =z for all i, j, and these three choices are all equivalent. On heterogeneous networks, these different rules have been found to lead to different scalings of the effective size N e with the actual network size N. For concreteness, we will mostly consider the above copying schemes. However, we stress that many of our results are quite general, and are not restricted to these specific choices.

Discrete-and continuous-frequency models
In the original formulation of the voter model, there is one individual per site: M = 1. Then, the frequency of a given species on a site can take only one of two values: 0 or 1. To obtain a diffusion equation for the copying dynamics without the need to average over multiple sites (for example, all sites of the same degree, as has been done elsewhere [9,10,23,[25][26][27]), we will take the limit M → ∞, so that the frequency x i of a chosen species on a site becomes a continuous variable in the range [0, 1]. We will use the terminology discrete-frequency and continuous-frequency models to refer to the extreme cases M = 1 and M → ∞ respectively. It is not a priori obvious that results obtained for continuous-frequency models will apply to their discrete-frequency counterparts. However, it turns out that for any fixed, finite M, discrete-frequency dynamics are recovered in a limit where the copying rates between sites, µ ij where i = j, vanish. This fact can be deduced from a convergence theorem for Markov chains due to Möhle [37]. The basic physical idea behind this result is that in the time it takes for an individual to arrive from a neighbouring site, sufficiently many replacements have been made by sampling from the local population of M individuals that each site is in one of its locally absorbing states-that is, each site hosts only a single species. As a consequence, the frequency of each species on any chosen site is always close to 0 or 1 (although globally it may be somewhere in between) and thus the discrete-frequency dynamics is recovered.
The way in which the M = 1 discrete-frequency dynamics is obtained mathematically is summarised in Appendix A. We illustrate here with the concrete examples of the previous subsection. Discrete-frequency versions of voter, invasion or link dynamics can be obtained by taking the parameter µ that appears in (3)-(5) to zero, whilst holding all µ ii at fixed, finite (but otherwise arbitrary) values. To recover the M = 1 version of the voter, invasion or link dynamics, with the unit of time defined such that each site is visited on average once per timestep, one must multiply all times by µ/M (see Appendix A).
A key point is that for any finite M, no matter how large, the discrete-frequency dynamics are recovered for sufficiently small µ. That is, any difference between a set of frequencies {x i } in the M > 1 model and their values in corresponding realisations of the M = 1 model can be made arbitrarily small by reducing µ. In the next section, and the remainder of this work, we will be working in the M → ∞ limit, such that all frequencies x i can be regarded as continuous variables. We will assume that results for discrete frequency models can then be obtained by subsequently taking µ → 0. In other words, we assume that the limits M → ∞ and µ → 0 commute. We do not know a priori that this assumption is valid. In practice, we find that we recover known results for the voter model, or good agreement with simulation, when we perform this procedure (see Section 7).

Forward-time formulation as a Fokker-Planck equation
We now derive the Fokker-Planck equation for probability distribution over the Ndimensional space of local frequencies x i of a chosen species. Although the method used to derive this equation, based on jump moments, is standard [38], we will recapitulate it here as it central to the derivation Eq. (1) that follows in Section 5 below. The key point is that the Fokker-Planck equation has the expansion where τ is the time variable and the jump moments are defined as a (n) In calculating the jump moments, the initial coordinates x i (τ ) are assumed known and fixed, and the angle brackets denote an average over all possible realisations of the dynamics in the time interval [τ, τ + δτ ]. These averages can be computed from the underlying stochastic dynamical rules. Within the present class of random-copying models, a local frequency x i may increase or decrease by an amount 1 M in any single update. The continuous-time transition rates for these processes are In the time interval [τ, δτ ], only events involving a single change to one of the x i contributes to order δτ in the average in (7), for any value of n. Thus, for the first jump moment we have, and for the second moments while a (2) ij = 0 for i = j, since two different x i never change at the same time. If all µ ij remain fixed as M → ∞, the second jump moments vanish relative to the first. The resulting Fokker-Planck equation then has only a drift (deterministic) term, and no diffusive component. Both terms are contribute in this limit if one takes µ ij to vanish as 1/M for all j = i, as is standard in population genetics [12,16]. We thus define rescaled rates m ij and c i through By introducing a rescaled time t = τ /M 2 we find as M → ∞ the Fokker-Planck equation The third and higher jump moments are all of order 1/M 3 and higher, and thus do not contribute in the M → ∞ limit. Although this Fokker-Planck equation is strictly valid in this limit, it is understood that it serves as an excellent approximation for large but finite M [16]. Our experience suggests that M ∼ 100 yields results consistent with infinite M (see Section 7.2 for a comparison between finite M simulations and the infinite-M theory).
In principle, to study the ordering dynamics of random-copying processes one would need to solve the full time-dependence of this N-dimensional diffusion. It is therefore i j μ ij Figure 2. Backward-time formulation of the random-copying dynamics as a random walk process. Shown is the configuration of the system of Figure 1 after the copying event has occurred. Looking backward in time, the tagged individual on site i, shown with a solid fill, has a probability µ ij per unit time of having been copied from site j. The tagged individual thus hops to site j at rate µ ij . apparent that the reduction to the single coordinate diffusion (1), when it is possible, represents a dramatic reduction in the complexity of the problem.

Backward-time formulation as a particle reaction process
As mentioned in the introduction, the copying process defined in Section 2 has an equivalent formulation as a particle reaction dynamics (see e.g. [28,29] for mathematical formulations and [12] for a more detailed account from a statistical physics perspective than the brief outline we give here). This process is identified by looking backward in time at the history of an individual in the current configuration. At some point in the past, this individual was created by the copying process. The rate at which an individual currently on site i was created by copying an individual on site j is µ ij /M, where the factor M appears because µ ij is the total rate at which one of the M individuals on site i is replaced by an offspring arriving from site j. This individual we call the ancestor of that in the current configuration. As we look further back in time, the ancestor performs a random walk, hopping from site i to site j at rate µ ij . See Figure 2. We will also refer to these ancestors as particles.
One can also consider the ancestry of multiple individuals. These independently perform random walks, as just described. However there is also the possibility that at some point in the history, one ancestor is the direct offspring of one of the others. The rate at which one of a pair of ancestors occupying sites i and j was copied from the other is (µ ij + µ ji )/M 2 . This time, the factor M 2 appears because both parent and offspring have a probability of 1/M of being involved in the copying process. After this copying event, the two lines of ancestry merge: one individual is the common ancestor of the original pair. This interaction can be viewed as a particle coalescence reaction.
To convert these rates into the system of rescaled time units defined in the previous section, we multiply them by M 2 . We then find that an ancestor hops from site i to j at rate m ij , and a pair of particles on site i react at rate c i . The rate for particles on different sites i = j to react vanishes under this rescaling.
Within this backward-time picture, it is important to realise that if one starts with an arbitrary pair of individuals in a finite system and looks far enough back in time the ancestor of one individual will at some point be the direct descendant of the other. That is, any two particles eventually react, leading to an overall reduction in the number of particles over time (i.e., fewer ancestors who contribute descendants to the present-day population). The steady state of the particle reaction process thus comprises a single common ancestor of any set of individuals chosen from the present-time population. If there is a single ancestor common to the whole population, then it follows that all individuals at the present time are of the same species-i.e., that complete order has been reached at the present time. Thus the stationary absorbing state of full order in the forward-time dynamics corresponds directly to the stationary state of a single particle in the backward-time picture.
In the remainder of this work we need only to consider the backward-time dynamics starting from a single pair of particles. Two key properties of this reaction-diffusion process are the probability distributions Q k (t) that a single particle on site i at time t = 0 occupies site k at time t and Q kℓ (t) that a pair of particles occupy sites k and ℓ. Note that in the reaction dynamics, we take increasing time to mean further back into the past, with t = 0 indicating the present time. When talking about a particle pair, we will consider the particles to be distinguishable, so that the configuration kℓ is not the same as ℓk if k = ℓ. These conventions make writing down the master equations for these probability distributions slightly more straightforward.
The master equation governing the single-particle distribution function is d dt Here we have introduced the matrix where the sign of the entries has been chosen so that the eigenvalues of M (1) all have nonnegative real parts. Assuming that the network is such that the stationary distribution Q k = lim t→∞ Q k (t) is unique, we have at late times where λ (1) r is the eigenvalue of M (1) with the smallest, nonzero real part. Meanwhile, the joint probability distribution for a pair of unreacted particles is governed by the master equation The first two terms in this equation give the contributions to this probability from particles hopping into, or out of, sites k and ℓ. The final term arises from the decrease in the probability that the two particles have not yet reacted due to the coalescence reaction. We can write this equation more compactly as d dt where the matrix M Here we use δ ijkℓ as a shorthand for δ ij δ jk δ kℓ , i.e., the object that is zero unless all four indices are the same, in which is it equals unity. On a finite network where the stationary distribution for a single particle is unique, the two particles will eventually react. Therefore Q kℓ (t) vanishes as t → ∞. The matrix M (2) has eigenvalues whose real parts are all positive. The smallest eigenvalue λ (2) q is strictly positive, and is equal to the asymptotic rate at which a pair of unreacted particles coalesce. The associated eigenvector we will denote Q kℓ and normalise so that kℓ Q kℓ = 1. This is the quasi-stationary distribution of the two-particle system. That is, if we condition on the two particles never having reacted, the probability they are found on sites k and ℓ approach the stationary value Q kℓ [39]. The rate of relaxation to the quasi-stationary state is given by the real part of the second-smallest eigenvalue λ (2) r . That is, at late times As stated in the introduction, and to be demonstrated below, it is a separation of timescales between the relaxation to and decay of the quasi-stationary state, controlled by λ (2) r and λ (2) q respectively, that allows for the reduction of the full Fokker-Planck equation (13) to the single-coordinate diffusion of Eq. (1).

Exact reduction to a single-coordinate diffusion equation
Now that the forward-and backward-time formulations of the model defined in Section 2 have been set up, we may derive the single-coordinate diffusion (1) and establish the condition for its validity. This amounts to identifying the quantity ξ whose ensembleaverage is conserved by the dynamics, and calculating its second jump moment. The procedure is not hard, but it does involve several manipulations that require explanation.
The key quantity ξ has been identified a number of times before (see e.g., [9,10,20,24]). It is where we recall that Q i is the stationary distribution of the single-particle random walk, and x i is the frequency of the species of interest on site i. To show that the diffusion equation (1) has no deterministic term, we calculate the first jump moment By multiplying (13) by x i and integrating over all the x coordinates (see [24,40]), we find d dt where M (1) is as defined in Eq. (15). Hence, where the second equality follows because Q i is the stationary distribution of the singleparticle random walk dynamics. In other words, it is the left eigenvector of the matrix M (1) with eigenvalue zero. Thus the first jump moment of ξ vanishes, as required.
We now turn to the second jump moment of ξ. For brevity, we define Notice that we take all the x i variables to have known values at time t. There is a slight complication here: only the value of ξ is known at time t, not the individual x i . However, we will shortly find that at late times a (2) (ξ) depends only on ξ, and hence the resulting expression will apply no matter what the actual individual values of x i are in any given realisation of the dynamics. As we have shown above, ξ(t) has a zero time derivative, so only the first term in the previous expression contributes to the second jump moment. To make progress, it is useful to introduce the quantity This is the probability that two individuals randomly selected from sites i and j are of different species. On a finite network, this quantity asymptotically decays to zero, and in fact, its time evolution is governed by the equation where M (2) is the transition matrix for the two-particle reaction process defined by Eq. (19). One way to see this is to calculate the time derivative of x i x j by multiplying (13) by x i x j and integrating, as described previously. A more direct way is to consider events that cause the probability that two individuals are of different species to change in a time interval [t, t + δt] (see, e.g., [41]). The rate at which P ij (t) changes due to a parent on site k leaving an offspring on site i is m ik [P kj (t) − P ij (t)]. Meanwhile, if both individuals are on the same site i, there is a probability c i per unit time that one of them replaces the other (see Section 4), an event that causes the probability that the two individuals are distinct to decrease at rate c i P ii (t). Summing over all events that change P ij we find which is equivalent to (30).
Since the time evolution of P ij (t) and the particle pair distribution function Q kℓ (t) discussed in Section 4 are both determined by the same matrix M (2) , P ij (t) has a latetime expansion of the same general form as that given by Eq. (20) for Q kℓ (t). That is, Now let us suppose there is some limit, for example N → ∞, in which Then, under a further rescaling of time t ′ = λ (2) q t, we have in this limit at any t ′ > 0 In other words, at times t ′ > 0, the system is assumed to be in the quasi-stationary state. To avoid too many different time units, we shall remain in the unscaled time t and assume that the quasi-stationarity condition (33) is satisfied, and thus that d dt P ij (t) = −λ (2) q P ij (t) when t > 0. To complete the derivation of the second jump moment a (2) (ξ), we differentiate both sides of (29) with respect to time t, rearrange and substitute into (28). We find To arrive at the second line we used the fact that the average of ξ has zero time derivative, that (33) holds and t > 0. The derivation is completed using the definition of P ij (t), Eq. (29), and of ξ, Eq. (21). As advertised, this second jump moment is a function only of the collective coordinate ξ. Using the expansion (6) we find that the motion of the collective coordinate ξ is governed by the Fokker-Planck equation at any time t > 0 in a limit where the condition (33) is satisfied. This is Eq. (1) with an effective size where we recall m is the baseline rate of the copying process. This concludes the derivation of the main result of this work. In the remaining sections we focus on strategies for determining when the key condition (33) is satisfied, and for estimating the single parameter λ (2) q (or equivalently, N e ) that depends on the network structure and choice of copying dynamics when the diffusion (1) or (38) is valid.
We conclude this section with a couple of remarks on the derivation of (38).
Behaviour at time t = 0 At very early times, and technically only at time t = 0 in the limit where (33) holds, there are additional contributions to the second jump moment.
Whether looking forwards or backwards in time, their net effect is to project the initial condition onto the quasi-stationary state.
To be more precise, if the leading eigenvectors of the matrix M (2) are normalised so that the probability that two particles initially occupying sites i and j have not reacted by time t and occupy sites k and ℓ is if all but the quasi-stationary modes are neglected at times t > 0. By taking the limit t → 0 from above, and summing over all k and ℓ, we find that P ij is the probability that this particle pair does not react before the quasi-stationary state is reached. This immediate collapse onto the quasi-stationary state has implications on the initial condition for the diffusion equation (1), given an initial set of frequencies x i (0). It is possible to calculate the variance in ξ at the onset of quasi-stationarity by appealing to (29). In the simple, solvable examples discussed in Section 7 below where we find (1) to be valid (Wright's island model and the bipartite network), we find that for any given initial combination of is appropriate. Whether this is generally the case is a question we leave open for future work.
Diffusion of other collective coordinates It turns out that any collective coordinate of the form where the ω i are an arbitrary set of positive weights that sum to unity, satisfies the diffusion equation (1) as long as, in addition to (33) the condition is also satisfied. This can be seen by repeating the above derivation of the jump moments with this more general set of weights. Ultimately, one finds The leading decay modes of the first jump moment is governed by the eigenvalue λ (1) r , while the modes corresponding to λ (2) q , λ (2) r and λ (1) r contribute to the second. Thus we require both (33) and (43) in order to obtain a vanishing first jump moment, and a time-dependent second jump moment of the form at times t > 0. Although all collective coordinates of the form (42) have the same diffusion equation-and crucially, the same effective size N e -when these conditions hold, we prefer the definition (21) because its first jump moment vanishes at all times, not just in the quasi-stationary state, and only the two-particle decay modes contribute to the second, as opposed to a mixture of the one-and two-particle modes. As a consequence, the condition for validity of the diffusion (21) is slightly weaker and cleaner, as it makes reference only to the two leading decay modes within the twoparticle sector. Furthermore, the effect of the dynamics at time t = 0 before the onset of the quasi-stationary state is harder to predict when the first jump moment does not vanish identically.

Methods of analysis for models with time-reversal symmetry
So far in this work, we have shown that if there is a separation of timescales in randomcopying models on networks between the rate of relaxation to the quasi-stationary state and the lifetime of that quasi-stationary state, the exactly solvable diffusion equation (1) correctly describes the ordering dynamics at all but the earliest times. To demonstrate the validity of (1), and compute the crucial parameter N e , one must in principle find the two smallest eigenvalues of the two-particle reaction matrix M (2) and verify the condition (33). Whilst this represents a major simplification of the full dynamics, there is of course no simple, general expression for these eigenvalues for an arbitrary set of transition rates m ij and c i .
In this section, we describe simplifications that arise in models where the singleparticle hopping process exhibits time-reversal symmetry (detailed balance). In particular, this allows one to relate the relaxation rate to the two-particle quasistationary state, λ (2) r , to the relaxation rate to the single-particle stationary state λ (1) r . This leads to a weaker, but more easily verified, condition for the validity of the diffusion (1). Furthermore, the symmetry allows one to develop a systematic variational method to bound the asymptotic reaction rate λ (2) q from above, or equivalently, the effective size N e from below. We develop these methods in this section, and apply them to specific combinations of network structure and update rules in the next.

Time-reversal symmetry and its consequences
The stationary state of the single-particle random walk is time-reversal symmetric if the detailed balance condition is satisfied. Although this may seem an overly restrictive constraint on the rates m ij , it applies, for example, if they have the general form m ij = mA ij α i β j , where A ij is the adjacency matrix and α i and β j are arbitrary positive site-dependent quantities. The three sets of rules discussed in Section 2.1 fall into this category, and so the methods we derive here will be applicable to these cases of interest. The relation (47) implies a corresponding symmetry in the matrix M (2) : As a consequence of this symmetry, all eigenvalues λ of M (2) are real, and for any right [42]. The eigenvectors form an orthonormal set, i.e., In particular, this means that the key quantities P ij and Q ij characterising the quasistationary state are related by Q ij = AQ i Q j P ij , where the constant A is fixed by the normalisation conditions ij Q ij = ij P ij Q ij = 1. These properties are crucial in the analysis of the quasi-stationary state that follows.

Relationship between one-and two-particle relaxation rates
We now show that, as a consequence of (48), λ (2) r ≥ λ (1) r , and hence that the condition (43) implies the condition (33) on the validity of the diffusion equation (1). The condition (43) may be easier to demonstrate than (33), since the relaxation rate of the single-particle dynamics is in principle more accessible than that for the two-particle dynamics. Indeed there are methods to bound λ (1) r from above [43] which would provide a means to demonstrate that (43) is satisfied, although we will not use them here.
One way to show that λ (2) r ≥ λ (1) r is to use standard perturbation theory [42]. We assume we have the full set of eigenvalues and eigenvectors of the matrix M (2) for some combination of m ij and c i , and then increase the coalescence rate on site x by an amount δc. The perturbation of M (2) is then δM (2) ij,kℓ = δc if i = j = k = ℓ = x, and zero otherwise. Then, for an eigenvalue λ of M (2) we have from first order perturbation theory [42] δλ = δc ijkℓ ψ λ ij δM This is a nonnegative quantity. Hence, in particular, λ (2) r must increase (or remain constant) as reaction rates c i are increased from zero. When all c i = 0, the reactions have no effect, and the two particles independently evolve according to the single-particle process governed by M (1) . Hence, the second smallest eigenvalue of M (2) equals that of M (1) when all c i = 0. Since λ (2) r cannot decrease as the reactions are turned on, we have λ (2) r ≥ λ (1) r , as claimed.

Variational approach for bounding the asymptotic coalescence rate
Since the eigenvalue λ (2) q that sets the characteristic timescale of the diffusion (1) is the smallest of a set of real eigenvalues, we can use a variational principle [42] to bound it from above. Let us denote such a bound λ + . The utility of this is twofold. First, if one can show that in some limit then the condition (33) for the validity of (1) follows. Second, the expression for λ + can be used as an estimate of the asymptotic reaction rate λ (2) q . One can construct a bound λ + from an arbitrary vector φ ij through the expression which is obtained by expanding the vector φ ij as a sum over the eigenvectors of the matrix M (2) (see, e.g., [42] for a derivation). There are two key points here. First, this result does not rely on being able to solve the two-particle reaction process. Any choice of φ ij leads to an expression for λ + . Even if this choice departs dramatically from the exact quasi-stationary eigenvector P ij , it may still be possible to show that one of the ratios (51) vanishes, in turn assuring validity of (1). The second point is that, in the usual way, one can include variational parameters in the expression for φ ij , and minimise λ + with respect to these parameters to obtain the best estimate of λ (2) q available within with the space of vectors spanned by them. We will encounter explicit examples in the next section, details of which are in the Appendix.

Applications
In order to make concrete the ideas and methods developed in the previous two sections, we now discuss their application to specific combinations of copying dynamics and network structures. In addition to checking certain known limits, we will use the fact that our approach does not invoke a coarse-graining over multiple lattice sites to obtain new theoretical results for the effect of correlations on the effective size N e .

Wright's island model
A spatially-structured random-copying model that is prominent in the population genetics literature is Wright's island model [44]. In this model, each site receives an offspring from one of the other sites (chosen uniformly) at a constant rate m. That is, m ij = m N −1 for all i = j. The coalescence rate is also taken to be uniform: c i = c for all i. The limited spatial structure in this model admits exact calculation of the eigenvalues λ (1) r and λ (2) q , and therewith demonstration that (1) is exact in the N → ∞ limit. Despite the venerability of Wright's island model, we believe that the validity of the diffusion equation (1) has not been shown explicitly before, although the resulting expression for the effective population size N e is well known (see e.g., [32]; also [45] for a review).
We begin by examining the single-particle dynamics. Since all sites are equivalent, the stationary distribution is uniform: Q i = 1 N . The rate of relaxation to this steady state, λ (1) r , is the smallest nonzero eigenvalue of the eigenvector M ij = m N −1 (Nδ ij − 1). One can verify by inspection that any N-dimensional vector will a single entry equal to 1, another equal to −1 and the rest zero is an eigenvector of M (1) . Since there are N − 1 linearly-independent vectors of this form, all N − 1 nonzero eigenvalues are degenerate, and one finds We now turn to the two-particle dynamics. It is possible to calculate λ (2) q exactly. One way to do this is with the variational method described above. We note that in the quasi-stationary state, all knowledge of the initial condition is lost. Thus there are only two distinct two-particle configurations in the quasi-stationary state: one in which both particles are on the same site, and one in which they are on distinct sites. To construct the ansatz, we set all diagonal entries φ ii = α and all off-diagonal entries equal to unity. The overall normalisation of the vector is unimportant, since this cancels in the numerator and denominator of (52). Thus we can without loss of generality fix one of the variational parameters and optimise the bound (52) with respect to those that remain. Since the exact quasi-stationary eigenvector is obtained with a particular value of α, it follows that by minimising the right-hand of (52) with respect to α will give us the exact decay rate of the quasi-stationary state.
In the Appendix, we show that in the quasi-stationary state, α and λ (2) q are given by the solution of the pair of equations Although these equations admit two solutions for λ (2) q , only one is positive. This is We can now use condition (43) to identify limits in which the diffusion (1) with effective size N e = 2m/λ (2) q is valid. One is the fast-migration limit, c m → 0, for which exact results are available on arbitrary network structures [46]. For Wright's island model we have indicating that in this limit, (1) is valid for any N. In this regime, λ (2) q ∼ c/N. This is consistent with the general result obtained by other means [46], and thus serves as a basic check of the methods derived here.
For general c and m, we have in the infinite system size limit, N → ∞, that Applying the condition (43) we find that the diffusion (1) is valid for any combination of m and c with an effective size N e = N(1 + 2 m c ). This is in agreement with the large-N expression for the effective population size quoted in the literature (see, e.g., [32,45]). To achieve this agreement, however, one needs to match up the choices of time units and copying rates that are used in standard population genetics models: the correct choice has the literature effective size related to our own by N (lit) e = N e /4m and the reaction rate c = 1/2M where M is the subpopulation size. In these units, m is the fraction of individuals replaced by immigrants in each generation. What is perhaps more significant is that, as far as we are aware, the applicability of this single effective population size over extended timescales has not been demonstrated explicitly even for this extremely simple model, even though it is routinely assumed. We remark that in the limit m → 0 where discrete-frequency models, like the original voter model, are recovered (see section 2.2), we find that N e = N, as expected. Thus at least in this case, the limits of infinite subpopulation size M and vanishing migration rates m indeed commute.

Bipartite network
The bipartite network, previously examined in [10,23,26], is insightful because it is perhaps the simplest network structure for which the three rules of Section 2.1 are not equivalent. Furthermore, as we will see below, degree correlations can be relevant, in the sense that they can lead to departures from results previously obtained by assuming an absence of such correlations [10,19,23,27].
The N nodes of the bipartite network are divided into two groups, A and B, with a number n A of the nodes in group A, and n B in group B. The rate m ij at which an individual i from group A is replaced from a copy of an individual j from group B is m ij = m AB . At rate m BA , copying takes place in the reverse direction. As in Wright's island model, the rate c is taken to be the same on all nodes. Note that individuals are never copied between two nodes belonging to the same group. See Figure 3.

Eigenvalue
Degeneracy Left eigenvector Table 1. Eigenvalues of the single-particle random walk on the bipartite network. (x; y) denotes a N -dimensional vector whose first n A elements are given by x and the remaining n B elements by y. The notation [1, −1] indicates a vector of the required dimensionality that has one element equal to 1, another equal to −1 and the rest equal to 0.
To analyse the copying dynamics on this network, we follow the strategy used above for Wright's island model, starting with the single-particle dynamics. One can verify by inspection that the eigenvalues and eigenvectors of the matrix M (1) have the structure presented in Table 1. Note in particular that the elements of the zero left eigenvector give the stationary distribution Q i (up to a normalisation constant). Hence, Meanwhile, the second-smallest eigenvalue is To bound the decay rate of the quasi-stationary state, λ (2) q , we make use of an extension to the variational ansatz described in the previous section (see Appendix for details). This extended ansatz still takes the quasi-stationary eigenvector φ ij = 1 when two particles occupy different sites i = j. However, it allows for site-dependent values φ ii = α i when both particles occupy the same site. The physical idea behind this ansatz is that the reaction process most strongly affects the probability of finding two particles on the same site, and has a weaker effect elsewhere.
A key quantity in this analysis is the inverse participation ratio, which is defined as Roughly speaking, I −1 measures the number of sites on which the particle may be found in the stationary state (see [47] for further discussion). If I vanishes as N → ∞, the stationary distribution is delocalised, and λ (2) q can be shown to be asymptotically bounded from above by where m i = j =i m ij is the total rate at which a particle departs from site i, and c i is the coalescence rate on site i. Furthermore, if λ * + also vanishes with N, the expression given above is the optimal bound on λ (2) q within the space of quasi-stationary eigenvectors admitted by the ansatz. These result are derived the Appendix. We note also that (62) coincides with an estimate for the effective size that we obtained previously by more approximate means [24].
In the case of the bipartite network, m i takes on two values: if the inverse participation ratio vanishes with N. We can use this bound, along with (60), in (51) to establish cases where the diffusion equation (1) is valid for the bipartite network.
In fact, one can obtain λ (2) q exactly for the bipartite network by following the same procedure as for Wright's island model-that is, by recognising that there are five distinct two-particle configurations in the quasi-stationary state, and identifying each with a variational parameter. If one follows this path, one finds that when I vanishes, so does the ratio (43), and that λ (2) q = λ * + asymptotically. That is, we reach exactly the same conclusions with this more complicated approach, so we present only the results obtained using the simpler analysis just described.
Of the three rules discussed in Section 2.1, the voter dynamics and invasion process are the most interesting. We discuss these two cases separately.

Voter dynamics Within voter dynamics, we have
The rate of relaxation to the single-particle steady state is simply λ (1) r = m. Meanwhile the inverse participation ratio I = N 4n A n B , which vanishes as long as both n A and n B = N − n A grow with N as N → ∞. When this is the case (51) is satisfied, and the asymptotic rate of decay of the quasi-stationary state is  . Mean reaction rate per particle pair as a function of time for coalescing random walks on bipartite graphs with hop rates chosen to correspond to the voter model and infinite on-site reaction rate. Symbols show the mean reaction rate within an n-particle sector plotted at the mean time t n between entering and exiting that sector. These averages were taken over 100000 realisations of the reaction dynamics on any given network structure. Time has been rescaled according to t 2 to allow direct comparison of two different bipartite network structures: one has equal group sizes whilst in the other they differ greatly.
In fact, since in the standard voter dynamics, the total rate of copying into a site, m i , is the same for every site, λ * + is proportional to the inverse participation ratio on any network. We remark that in the limit m c → 0 in which the traditional discretefrequency voter model (see Section 2.2) is expected to be recovered, the result (65) corresponds with that obtained for the voter model on bipartite networks obtained by other means [23,26].
If, say, n A remains finite as N → ∞, the inverse participation ratio no longer vanishes. The stationary distribution is localised within the smaller group. As a consequence, condition (51) is no longer satisfied, suggesting (but not necessarily implying) that there is no separation of timescales in the two-particle dynamics. We test this explicitly with Monte Carlo simulations of the backward-time particle reaction dynamics. Initially, we place one particle on each site, and use the Gillespie algorithm [48] to determine the next particle to hop and the time that it does so. If two particles occupy a site simultaneously, they immediately coalesce: thus we are directly simulating the limit m c → 0. In Figure 4, we show an empirical measure of the pairwise reaction rate as a function of time. This is achieved by taking the reciprocal of the mean time spent in the n-particle state (averaged over multiple runs), and dividing by n 2 to convert the total reaction rate to a rate per particle pair. As can be seen from its asymptotic value.
Behaviour in the limit N → ∞ is revealed by Figure 5. Here, a fit to the asymptotic reaction rate as a function of N is compared for the case n A = const versus n A ∝ N. This is done by fitting a constant function to data series like those shown in Figure 4. The error bar indicates the standard deviation of the set of rates: this gives an indication of how much transients are contributing to the average. As can be seen, when n A remains constant as N increases, there is no evidence for a decrease in the contribution of the transient, suggesting that the early-time relaxational dynamics cannot be neglected in the limit N → ∞, and thus that the diffusion (1) is indeed not valid in this case.
We also examine the effect of intermediate subpopulation size M (see Section 2 for the precise definition) on ordering times. We recall that the limit M → ∞ is used to obtain the Fokker-Planck equation (13). At finite M it is possible that additional terms may contribute. To investigate this, we performed direct Monte Carlo simulations of the forward-time random-copying dynamics specified in Section 2 with copying rates given by (3). We took the overall rate of copying between sites µ = 1 M , so that then m ij = limit M → ∞, the effective size appearing in (1) is found from (65) to be From a solution of a backward equation corresponding to (1), as described for example in [23,49], one finds that the mean time T o for complete order to be reached (i.e, for ξ to reach 0 or 1) from the initial condition ξ = 1 2 is given by T o = N e ln 2. The result of performing many repetitions of the forward-time copying dynamics until order is reached, and averaging the ordering time, is shown in Figure 6 for various M. We see that even with M = 10, agreement with the M → ∞ result is reasonable. Closer inspection of the data reveals that T o is systematically overestimated in this small M regime; however, already for M = 80 the data are statistically indistinguishable from the infinite M prediction. Thus we find that even relatively small subpopulations are reasonably well-described by the continuous frequencies x i .

Invasion process
In the invasion process, the hop rates of Eq. (64) are exchanged. The relaxation rate of the single-particle random walk is Meanwhile, the inverse participation ratio is given by which vanishes as N → ∞ for any combination of n A and n B , as can be verified by taking n A ∼ aN γ with 0 ≤ γ ≤ 1 and n B = N − n A . Then, (63) is asymptotically optimal. Using (51), we then find, again for any combination of n A and n B , that the single-coordinate diffusion (1) applies for the invasion process on the bipartite graph.
For γ < 1 in the above group-size definitions, the estimate for the decay rate of the quasi-stationary state takes the particularly simple large-N form We compare the predictions of this formula in the limit m c → 0 against simulations of the type discussed above in the context of voter dynamics. As can be seen from Figure 7, numerical measurements of the asymptotic decay rate of the quasi-stationary state agree well with the prediction from Eq. (69). Again this provides evidence that the dynamics of discrete-frequency models is well-represented by the m → 0 limit of their continuous-frequency counterparts.
We remark that the result (69) differs from that previously obtained within an approximation that assumes an absence of degree correlations in the network [10]. Thus the invasion process on the bipartite network is a simple, illustrative example of a case where degree correlations must be taken into account to obtain the correct reduction onto a single-coordinate diffusion of the form (1). Although the correct result was obtained in [26], the separation of timescales required for its validity was assumed, rather than explicitly demonstrated.

Correlations in random network structures
We now explore this last point in a little more detail by investigating a transformation to a network structure that induces degree anti-correlations. We achieve this by beginning with some template network and then 'exploding' it by inserting a node in the middle of each edge of the template. See Figure 8. This way, every node of the template network is then next to a node of degree two: thus high-degree nodes end up next to low-degree nodes.
Clearly, one consequence of this transformation is to increase the size of the network, with a likely corresponding increase in the effective size N e . Our aim in this final section is to demonstrate that degree anti-correlations lead to a further increase in N e beyond what is predicted by (62).
As template networks we use Molloy-Reed networks [50] with an approximate power-law degree distribution. These networks are designed to have weak degree correlations. They are constructed by assigning to each node i its desired degree by attaching z i half-edges to it. Then, these half-edges are formed into complete edges by picking a random pair and joining them together if the two nodes in question are not already connected. Once all half-edges are joined, the probability that nodes i and j are connected is z i z j /K 1 where For this study, we use template networks with the degree distribution where A γ is a normalisation constant. For large z, this approaches the power-law distribution p z ∼ z −γ . We shall however mostly be looking at the behaviour on relatively small networks: the purpose of using the above form for p z is to allow the possibility of highly heterogeneous networks (i.e., a large variance in the degree z). We shall also, for simplicity, restrict ourselves to voter dynamics, i.e., the case where m i = m on every node, and the limit m c → 0. Then, asymptotically, we have from (65)  bound λ * + ∼ 2mI = 2m on λ (2) q , since Q i = z i K 1 under voter dynamics [20,23]. As previously, this bound is valid in the large N limit when I vanishes with N.
In Figure 9 we show empirical measurements of the asymptotic decay rate of the quasi-stationary state, divided by (72). We used two degree distributions, one with γ = 2.5, and one with γ = 3.5. In both cases the minimum degree of any node was 4. For the original template we find deviations from this prediction of 10-20%. Since (72) was obtained by assuming an absence of correlations between particles occupying different sites in the quasi-stationary two-particle state, it would appear that (at least at the network sizes simulated) such correlations exist, despite the fact that degrees of neighbouring nodes are designed to be weakly correlated in Molloy-Reed networks. We also see that the deviations from (72) increase dramatically when anti-correlations are introduced by the explosion procedure discussed above. Although we have not been able to demonstrate analytically that there is a separation of timescales in the dynamics on these networks, numerical data (not shown) suggest that the transient is irrelevant except possibly on the uncorrelated template network with γ = 2.5.
To understand more quantitatively the decrease of λ (2) q (and hence increase of N e ) under the explosion operation, we employ a refinement of the ansatz used to obtain (72) that allows for correlations between unreacted particles on neighbouring sites in the quasi-stationary state. Details are presented in the Appendix. The bound on λ (2) q that is ultimately obtained in the m c → 0 limit is, for large N, Here, S is a measure of how much high-degree nodes tend to be connected to other highdegree nodes: it is related to the assortativity coefficient of a network [51]. Meanwhile, T comprises a sum over triangles in the network, weighted by the total degree of the vertices. It is thus similar in spirit to a clustering coefficient [52] in that it vanishes if the network contains no triangles, that is, if locally it has a tree-like structure. The way these quantities interact in (73) is unfortunately not straightforward to interpret. However, the deviation from the baseline (72) will typically be large when degrees of neighbouring nodes are anticorrelated and/or triangles tend to be formed by high-degree nodes.
If nodes are uncorrelated, A ij can be replaced with its ensemble average value in the limit N → ∞, that is, with the probability z i z j /K 1 that nodes i and nodes j are connected [23]. In this case, one finds asymptotically that λ # + and λ * + converge. On the other hand, the bounds are distinct after the explosion operation. Suppose the template network has given values of K 1 and K 2 . Then, after the explosion, there are an additional K/2 additional nodes of degree two. This implies that for the exploded network we have K ′ = 2K and L ′ = L + 2K (where we use a prime to indicate the exploded network). To compute the statistic S ′ , we note that in the exploded graph, each original node i is connected to z i nodes of degree two. Taking into account that the sum over all neighbouring pairs i, j can be written as twice the sum over all edges, and that each edge is connected to one of the original template nodes, we find that By design, the exploded network contains no triangles, so T ′ = 0. Using these values in (73) we find, again in the m c → 0 limit, which we can contrast with the baseline λ * + ∼ 2mI ′ . Relative to this baseline we see that the asymptotic decay rate of the quasi-stationary state is reduced by a factor of at least 3 4 as a consequence of the degree anti-correlations in the exploded networks. In Figure 10 the same data as in Figure 9 are shown, but this time relative to the improved bound (73). We see that in all cases, the improved bound is closer to the  Figure 10. The asymptotic decay rates of Figure 9 rescaled relative to the improved estimate λ # of Eq. (73). Error bars are to be interpreted as in Figure 9. measured asymptotic decay rate, although significant discrepancies remain. The reason for this is that the ansatz for the quasi-stationary eigenvector φ ij does not contain sufficient structure to capture all the relevant correlations between a pair of unreacted particles. However, the fact that the correction factors shown in Figures 10 and 9 are only weakly dependent on system size N suggests that the correct scaling of the quasistationary reaction rate with N is captured by the expressions (72) and (73). This is consistent with previous works in which the scaling of N e with N has been established through fairly simple approximation schemes [9,10,19,23], further work having been necessary to improve on the prefactors [24,27]. Further improvements demand a more detailed understanding of the structure of the quasi-stationary two-particle state.

Discussion and outlook
In this work, we have shown that the onset of order within a large class of random copying processes (including the voter model) on spatially-structured networks can be described exactly in terms of the diffusion of a single collective coordinate ξ. This quantity is a weighted average of a species frequency over nodes of the network. When the reduction onto the single-coordinate diffusion process is possible, the only significant effect of the spatial structure relative to the exactly-solvable 'mean-field' diffusion of Equation (1) is to replace the actual network size N with an effective size N e .
As we have shown in Section 5, a remarkable feature of the analysis is that in order to establish the validity of the single-coordinate diffusion (1), one only needs to establish two properties of a two-particle reaction-diffusion process. The first is the lifetime of the quasi-stationary probability distribution for two unreacted particles. The second is a separation of timescales between this quasi-stationary lifetime and that of the relaxation to quasi-stationarity. In principle, to establish the probability that complete order has set in by a time t, one would need to compute an N-point correlation function, starting, for example, with the full N-dimensional diffusion equation (13). The fact that one can show that the one-dimensional diffusion (1) is valid by appealing only to a twoparticle problem represents a considerable simplification in the analysis of voter and other random-copying models.
In Section 7 we demonstrated the separation of timescales within certain exactlysolvable models, namely, Wright's island model and on the bipartite graph. Although this property has previously been assumed with varying degrees of explicitness within these models [10,24,26], we are not aware of this property actually being confirmed for these models before. Furthermore, many existing treatments of randomcopying dynamics on heterogeneous networks [9,10,23,[25][26][27] have derived a diffusion approximation by aggregating multiple sites in order to construct continuous spatial variables. Such approaches have the drawback that the network structure gets 'washed out' in the coarse-graining. Here we have described an alternative approach whereby a continuous variable is associated with each site, and the network structure is retained. As we have shown in Section 7, anti-correlations between the degree of neighbouring nodes leads in the example of the invasion process on a bipartite network to a scaling of the effective size N e with network size N that differs from that predicted by an approach that neglects these correlations [10]. For voter dynamics on random networks, the effect of degree correlations appears to be less severe, affecting the prefactor rather than the scaling exponent. Since the inverse participation ratio I seems to be the quantity that dominates the scaling of N e with N in these models, the physical origin of the sublinear scaling of N e with N that can be observed [9-11, 23, 26] lies in the localisation of the quasi-stationary single-particle distribution Q i in the vicinity of high-degree nodes.
In existing treatments of random-copying dynamics on heterogeneous networks, the methods used to arrive at the diffusion equation (1) have to some extent been ad-hoc.
Here we have formulated a more systematic treatment based on timescales of the twoparticle reaction-diffusion dynamics. The remaining task is to find good methods to estimate these timescales-and thereby demonstrate their separation. In cases where the single-particle reaction dynamics is time-reversal symmetric, we have shown that one can use the single-particle relaxation rate as a proxy for the corresponding quantity in the two-particle sector; we have also demonstrated the utility of a variational approach to bound the effective size N e from below. While this has allowed the separation of timescales to be demonstrated within very simple models, it was still necessary to resort to simulations to establish whether it was present on random network structures. Development of these methods-perhaps using bounds on the single-particle relaxation dynamics that have been derived elsewhere [43]-would be needed to determine the class of networks and dynamical processes for which (1) is valid.
We remark that (1) is known to apply to the voter model on square lattices in two dimensions and above [28]. In two dimensions, it is found that the effective population size grows logarithmically with system size (relative to the rate at which a particle relaxes to the stationary uniform distribution). This can be understood as a statement that domains coarsen logarithmically in the two-dimensional voter model [7]. Although the separation of timescales is weak in this case, it is nonetheless present. Thus one would expect the separation of timescales to be observed on random networks with small diameters, on the grounds that a single particle would be expected to explore the system on a faster timescale than would be possible on a two-dimensional square lattice of the same size. It would be interesting to pursue this line of enquiry further.
As mentioned in the introduction, the use of a single effective size N e to characterise purely stochastic population dynamics in subdivided populations is widespread in population genetics (see e.g., [45] for a review). A fair amount of discussion has centred around the fact that different definitions of N e are possible, depending on which property of the stochastic dynamics is measured and calibrated against an unstructured (mean-field) population [41]. Whilst it is widely understood that asymptotically many definitions converge onto the lifetime of the quasi-stationary state (also known as the eigenvalue effective size [53]), the importance of a separation of timescales between this lifetime and others seems to be less well established. For example, Sjödin et al [35], identify the asymptotic rate of particle coalescence as defining 'the' effective population size. In these analyses, they examine a separation of timescales in the definition of migration and coalescence rates, m ij and c i . That is, migration rates that vanish with local population size M with differing exponents. This approach potentially misses the possibility-demonstrated here in the examples of Section 7-that one can have a separation of emergent timescales that is sufficient to allow a single effective size that characterises almost the entire history of the dynamics. The importance of emergent timescales has perhaps been most clearly recognised by Wakeley and coworkers [33,34,36]. In particular, the rapid relaxation of the single-particle dynamics as a condition for a single, well-defined effective size was noted in [36]. However, that study was restricted to a rather special class of vertex-transitive networks and dynamical update rules. Furthermore, the irrelevance of all timescales other than the relaxation of the quasi-stationary state was implicitly assumed rather than explicitly demonstrated. We believe that the present analysis offers considerable additional insight into the applicability of a one-dimensional diffusion characterised by a single effective population size for a general class of random-copying dynamics on complex networks.
Nevertheless, a number of open questions remain. As we have already noted, the task of calculating the relevant decay rates λ (2) q and λ (2) r is in general challenging (despite being simpler than solving the full N-body dynamics). Although we have shown that certain simplifications arise when the single-particle random walk dynamics are timereversal symmetric, there is of course a much larger class of processes that do not respect this symmetry: these have nonequilibrium steady states. Here we have not been able to show a result analogous to that derived in Section 6.2, that is, to show that the relaxation rate to two-particle quasi-stationarity exceeds the relaxation rate to a single-particle nonequilibrium steady state, although intuitively one might expect this to be the case. Furthermore, we do not know of any convenient means to bound the asymptotic reaction rate in these models.
Although we have found that we recover known results for discrete-frequency models, or obtain good agreement with simulations, in the limit m c → 0, as we noted in Section 2.2, this correspondence relies on the limits of infinite local subpopulation size and vanishing migration rate to commute. Again, we are not aware of a proof that this should be the case. At the end of Section 5 we commented that the initial condition on the diffusion equation may be modified as the system evolves into the quasi-stationary state: further clarity here would be beneficial. Additionally, although the coordinate ξ is sufficient to track the degree of order present in a system, other aspects of the evolution cannot be determined from this quantity alone-for example, the typical state of the system or correlation length given some value of ξ. It is possible that some of these quantities can be obtained with knowledge of the quasi-stationary eigenvector φ ij , this by definition encapsulates two-point correlations. It is unclear whether higher-point correlation functions could be inferred from the two-point functions, or whether these would need to be determined separately.
Finally, we revisit the fundamental property shared by the entire set of randomcopying processes discussed here, namely, that all species are equivalent. There is considerable interest in those cases where one species is fitter than another. Although in some sense this is only a minor modification to the dynamics, the mathematics become much more involved. In particular, the backward-time dynamics are more complicated: the reaction rates then depend on the frequencies of different species in the population at the point in history at which the coalescence reaction takes place. It may, however, be possible to apply the results obtained here for the neutral theory in some kind of perturbation theory to study the effects of a small selective advantage. We leave this as one of many possibilities for future work. point in time, the state of the system, σ, is specified by the set of numbers n i which count the number of individuals of the species of interest on site i. In any one timestep, one of these numbers may change by ±1. The idea now is to split the matrix A of transition probabilities that appears in the master equation P σ ′ (t + δt) = σ A σ ′ ,σ P σ (t) into two pieces. The first, L, accounts for local interactions: that is, when an individual is replaced by sampling from the same site. The second, I, accounts for interactions between sites. Both are stochastic matrices. If the latter happen a fraction ǫ of the time, then The theorem in [37] then states that the vector of probabilities P (t) has the limiting form as ǫ → 0, subject to the existence of the projector The way to read (A.2) is as follows. The projector P performs an infinite number of local copying events on every site. If on site i there are initially n i individuals of the species of interest, then at the end of this process n i → M with probability n i M , and n i → 0 with probability 1 − n i M . Thus the initial condition is projected onto a state in which every site is occupied by individuals of a single species. The exponential term in (A.2) is the generator of a continuous-time process. If we take the limit δt → 0, a copy between sites is taken at rate ǫ per timestep. A single copy may leave the system unchanged; or it may change one of the n i values by ±1. In the latter case, all but one of the individuals are of the same species. Then, the projector acts, either to restore the existing value of n i (with probability 1 − 1 M ) or to 'flip' it (i.e., from M → 0 or vice versa) with probability 1 M . Putting this together, one finds that the frequency x i = n i M flips from 0 to 1 (or vice versa) due to a copy taken from a neighbouring site j at a rate ǫ M multiplied by the probability of copying from j to i in the original discrete-time process. In Section 2, the random-copying process was set up as a continuous-time process in which copies from j arriving at i at a rate µ ij . This is obtained here by putting the transition probabilities encoded in I equal to µ ij δt. Then, in the original time units, a copy between sites is taken at rate µ ′ ij = ǫ M µ ij . Rescaling time such that t ′ = ǫ M t, we obtain in the limit ǫ → 0 the M = 1 (discrete frequency) dynamics from any model with M > 1. Note that in the case of the dynamical rules (3)-(5) we can identify ǫ with µ.

Appendix B. Variational bounds on the quasi-stationary decay rate
As stated in the main text, when the backward-time hopping process is time-reversal symmetric, all eigenvalues of the matrix M (2) are real, and hence the smallest eigenvalue satisfies the bound for any vector φ ij . In order to construct bounds on λ (2) q one can use a variational approach. That is, one can include various parameters in φ ij , and minimise the righthand side of (B.1), thereby obtaining the most stringent bound on λ (2) q available within the space of vectors spanned by these parameters.
The key to designing a variational ansatz is to recall that the quasi-stationary distribution of two unreacted particles Q ij ∝ Q i Q j φ ij . Thus φ ij measures the extent to which the probability of finding a pair of particles on sites i and j is affected by the interaction between them. Alternatively, one can infer from the interpretation of P ij ∝ φ ij given at the end of Section 5 that φ ij is the probability that two particles starting on sites i and j do not react before the onset of the quasi-stationary state. Either way anticipate that φ ij will be close to unity for particles that are far apart, and decrease as their separation decreases.
We make use of two variational ansätze in the main text. In the first it is assumed that φ ij differs from unity only when both particles are on the same site. In the second we extend to include the effects of interactions when particles occupy neighbouring sites (but at the expense of restricting to voter-type update rules only).
Ansatz invoking on-site interactions only Restricting the effect of interactions to particles occupying the same site, corresponds to an ansatz for φ ij of the form φ ii = α i and φ ij = 1 if i = j. Differentiating (B.1) with respect to α m , and setting it to zero, we find that ij,kℓ φ kℓ ij Q i Q j φ 2 ij = λ + . where m i = j =i m ij , the total exit rate from site i in the random-walk picture. Substituting into (B.2), we find that This is not sufficient to determine the optimal α i , since the optimal value of λ + is not known. To compute this, we need an additional relation between them. This can be obtained as follows. First, we note that (B.2) can be written as ii,kℓ φ kℓ = λ + Q 2 i φ ii . (B.5) Meanwhile, if we temporarily set φ ij = α N +1 for i = j, differentiate with respect to this additional α parameter and set to zero, we find analogously that i =j ij,kℓ φ kℓ = λ + i =j Summing (B.5) over all i, and adding to this expression, we find that where the second equality holds because Q i is the zero left eigenvector of the matrix M (1) . Using the definition (19) and the form of the ansatz for φ ij , we ultimately obtain . (B.8) The optimal bound on λ (2) q is found by solving the simultaneous set of N + 1 equations given by (B .4) and (B.8).
For Wright's island model, discussed in Section 7.1, m i , c i and α i are all siteindependent quantities, and Q i = 1 N . Making these substitutions in (B.4) and (B.8), we arrive at the pair of equations (54) and (55) studied in the main text.
For more general models, solving (B.4) and (B.8) exactly is difficult. Nevertheless, one can still make some progress here. An alternative way to formulate this problem is to define Then, the optimal bound λ + is given by the value of x for which y(x) = x. Differentiating, we find dy (B.10) Hence, the root of y(x) = x must lie at a value of y(x) ≤ y(0); hence λ + ≤ y(0). Thus, one has that Note that the right-hand side of this expression is not necessarily the optimal bound. However, if the inverse participation ratio I = i Q 2 i vanishes as N → ∞, then asymptotically we obtain the simpler expression which was used in the main text. This bound is optimal within the present ansatz for φ ij if λ + also vanishes as N → ∞. Then, since (B.14) If both I = i Q 2 i and λ * + vanish as N → ∞, this lower bound on the optimal λ + asymptotically approaches the upper bound (B.12), and hence the right-hand side of (B.12) is optimal in this case.
Ansatz invoking on-site and nearest-neighbour interactions As noted in the main text, the traditional voter model dynamics is obtained by taking m i = m and c i = c followed by the limit m c → 0. In the above ansatz, this leads to all α i = 2m 2m+c and λ * + ∝ I. To examine the effects of degree correlations between neighbouring nodes, it is natural to introduce an ansatz for φ ij that takes one value when both particles are on the same site, a different value when they are on neighbouring sites, and unity elsewhere. That is, where A ij is the adjacency matrix for the network. Recalling that A ij = 0 if i = j, we see that this expression takes the value α if i = j, β if i and j are neighbours, and unity otherwise. For voter dynamics, we have from (3) and (19) that where z i is the degree of node i. Differentiating the right-hand side of (B.1) with respect to α and setting to zero, we find We notice the similarity with (B.4): the on-site correlation φ ii is multiplied by the nearest-neighbour correlation β. To obtain an expression for this correlation, we differentiate the right-hand side of (B.1) with respect to β, and set to zero. This yields This expression can be evaluated by using the fact that, for the voter model, the stationary probabilities are Q i = z i /K 1 . The sum on the right-hand side becomes 1 where K n and S are defined by Equations (70) and (74) respectively. The sum on the left requires a bit more work. One finds it can be written as Substituting (B. 19), and rearranging, we finally find Similar manipulations allow one to determine that When the inverse participation ratio I vanishes with N, we know that λ + does too, since (B.12) bounds it from above, and is proportional to I under voter dynamics. The upshot of this is that in this case, we can neglect the terms in λ + in ( is assumed to vanish as N → ∞, we ultimately obtain the expression (73) quoted in the main text.