Perturbative Field-Theoretical Analysis of Three-Species Cyclic Predator-Prey Models

We apply a perturbative Doi--Peliti field-theoretical analysis to the stochastic spatially extended symmetric Rock-Paper-Scissors (RPS) and May--Leonard (ML) models, in which three species compete cyclically. Compared to the two-species Lotka--Volterra predator-prey (LV) model, according to numerical simulations, these cyclical models appear to be less affected by intrinsic stochastic fluctuations. Indeed, we demonstrate that the qualitative features of the ML model are insensitive to intrinsic reaction noise. In contrast, and although not yet observed in numerical simulations, we find that the RPS model acquires significant fluctuation-induced renormalizations in the perturbative regime, similar to the LV model. We also study the formation of spatio-temporal structures in the framework of stability analysis and provide a clearcut explanation for the absence of spatial patterns in the RPS system, whereas the spontaneous emergence of spatio-temporal structures features prominently in the LV and the ML models.


Introduction
Population dynamics has been and continues to be an extremely active field of research since about forty years [1,2,3,4,5,6,7]. Steady progress in the development of mathematical and computational tools as well as the application of methods from statistical physics have allowed qualitative and quantitative insight into the behavior of interacting species. Various simplified models have been invoked to address prototypical situations in real ecosystems: The paradigmatic two-species Lotka-Volterra (LV) predator-prey model [8,9] was originally introduced to study fish population oscillations in the Adriatic sea, as well as to explain auto-catalytic chemical reaction cycles. The Rock-Paper-Scissors (RPS) model [4,10,11,12,13,14] addresses the case of three cyclically interacting species with a conserved total number of individuals, whereas the May-Leonard (ML) model [1,15,16,17,18,19,20] describes a more general, nonconserved situation. These models are obviously and necessarily rather simplified and lack many of the details of ecological neighborhoods. However, recent efforts aim at the realization and experimental implementation of such systems [21,22,23,24,25,26]. Furthermore, it is reasonable to assume that simplified constructs such as the LV, ML, and RPS systems should be useful as elementary motifs and building blocks of models for more extended ecosystems. It is therefore imperative to investigate which of their features are qualitatively and/or quantitatively robust and remain important when multiple interacting species are coupled to environments with richer structures.
Traditionally, species dynamics in ecosystems are modelled via coupled non-linear ordinary differential equations. In the case of spatially extended systems, this approach is generalized by using partial differential equations that represent species dispersion through simple diffusion, i.e., coupled reaction-diffusion equations. However, this meanfield or mass action approach fails to take into account the inherent randomness and stochastic nature of the underlying processes stemming from fluctuations in the discrete number of individuals, and neglects spatio-temporal correlations. Yet fluctuations and correlations can lead to dramatically different behavior than predicted by mean-field theory [27]. For example, the classical LV mean-field rate equations predict neutral cycles and hence non-linear oscillations around a marginal fixed point, while stochastic computer simulations of this system yield decaying oscillations towards a (quasi-)stable state [28,29,30]. This stationary state exhibits large and erratic excursions triggered by fluctuations in the species concentrations in zero-dimensional [31] as well as spatially extended systems [32]. Spatially extended stochastic LV model variants also show intriguing spatial patterns and moving activity fronts [29,33,34]. Crucially, stochastic variants of the LV model exhibit a large susceptibility to randomness in the predatorprey interaction rates [35,36].
Spatially extended cyclic models such as the RPS or ML systems are influenced by internal reaction noise and exhibit differences in species extinction times and resulting spiral pattern wavelengths compared to the mean-field approximation [13,18,37]. In one dimension, 'superdomains' may form in these cyclic models [38]. Although both models are cyclic in nature, they exhibit different sensitivity to stochastic fluctuations. The RPS model, a generalization of the LV model to three cyclically competing species, displays comparatively weak fluctuation renormalizations in the quasi-stable coexistence state and minimal modifications due to randomized reaction rates [14]. In contrast, the ML model features a stronger renormalization of the oscillation frequency in the unstable region where spiral structures form spontaneously, but appears to have an insignificant response to randomized reaction rates [17]. These observations from Monte Carlo simulations raise the intriguing question: Under what conditions will fluctuations significantly alter the system's properties and cause marked deviations from simple mean-field predictions? (c) (d) Figure 1. Snapshots of the spatial particle distribution in cyclic three-species RPS and ML models for single stochastic simulation runs (system size 100×100 lattice sites): Each lattice pixel is assigned an RGB value such that each color value is proportional to the number of individuals of a specific species. A color value 0 represents the absence of the species corresponding to that asigned color; therefore, black pixels indicate empty sites. Top: RPS model with reaction rate parameter λ = 0.5 at (a) t = 300 Monte Carlo Steps (MCS) and (b) t = 400 MCS; bottom: ML model with predation rate σ = 0.5 and reproduction rate µ = 0.5 at (c) t = 300 MCS and (d) t = 400 MCS, respectively. The red species predates on the blue species, the blue species on the green species, and the green species on the red species in both models.
To at least partially answer this question, a field-theoretical perturbation analysis was applied to the stochastic spatially extended LV model in Ref. [39]. To one-loop order, this semi-quantitative analysis confirms that i) the fluctuation-induced damping renders the system unstable against spatio-temporal structures, and ii) fluctuations significantly renormalize the oscillation frequency in the two-species co-existence phases, especially below three dimensions. Aiming to better understand the fluctuations in spatially extended RPS [ Fig. 1(a,b)] and ML models [ Fig. 1(c,d)], we utilize a similar Doi-Peliti field theory representation for their associated stochastic reaction processes. To study the impact of intrinsic fluctuations on system parameters, a one-loop calculation is carried out in the perturbative regime, where the reaction rates are small as compared with the diffusivity, and a thorough comparison between the RPS, ML, and LV systems is conducted. In contrast to earlier observations in numerical simulations, the RPS model exhibits noticeable fluctuation-induced corrections in the perturbative regime, similar to the LV model. We believe that, as the dissipation becomes non-negligible in the non-perturbative regime, the associated infra-red (IR) divergence is regularized, and thus substantial renormalizations become effectively suppressed. We note that in all investigated systems, the field-theoretic loop expansion technically only applies to the stable regions with spatially homogeneous ground states. Our results demonstrate that, at least in the stable region, the dynamical features of the ML model conversely do not receive significant modifications from fluctuations. Based on these explicit calculation results, we also provide pertinent arguments that explain the absence of spontaneous spatio-temporal patterns in the RPS model with conserved total population number, as opposed to the ML model, which for sufficiently large system sizes develops spiral oscillatory patterns, as depicted in Fig. 1(c,d).
The paper is organized as follows: Detailed perturbative field-theoretical analyses for the cyclic and symmetric RPS and ML models are performed in sections II and III, respectively, where we establish the Doi-Peliti functionals for both models and state their corresponding generalized Langevin equations. Renormalized damping coefficients, oscillation frequencies, as well as diffusivities are calculated up to one-loop order in the perturbative fluctuation expansion. In section IV, a comprehensive comparison between the LV, RPS, and ML models is provided, and pertinent distinctions between these paradigmatic systems are highlighted. Specifically, we discuss the influence of fluctuations and the stability of spatio-temporal structures, and also briefly address the effect of quenched disorder in the reaction rates. We conclude with a brief summary and outlook. Finally, Appendix A presents a succinct review of the Doi-Peliti field theory approach and also provides a brief analysis of the asymmetric RPS model, demonstrating its effective two-species limit for strong asymmetry at the mean-field level. The remaining appendices list additional technical and computational details for the symmetric ML model.

RPS model and mean-field rate equations
The RPS model consists of three particle species, subject to the cyclically coupled stochastic competition reactions In this paper, we consider the cyclic-symmetric case, such that λ 1 = λ 2 = λ 3 = λ . In this limit, the system displays a discrete S 3 symmetry among the three species. A brief analysis of the general asymmetric case is presented in Appendix A. We note that every species interacts via a standard non-linear Lotka-Volterra predation reaction with the subsequent species in the cycle, consuming a "prey" particle and reproducing at the same instant. The total number of individuals is unchanged by all reactions, hence particle number conservation holds globally and locally (except for hops to neighboring lattice sites, see below). We consider a model wherein particles from all three species perform random walks on a d-dimensional hyper-cubic lattice with L d sites and lattice constant c. We do not restrict the number of particles per lattice site, hence we do not consider finite local carrying capacities here (the total number of particles is fixed). The rate at which particles hop between sites is given by D/c 2 , where D denotes a macroscopic diffusion constant. The reactions (1) occur on-site, and only if two particles of differing species are present. Reaction products are put on the same lattice point as the reactants.
In the limit of large diffusivities (relative to the reaction rates λ ) the system can be considered well-mixed. Hence, the RPS rules can be approximated by the three coupled mean-field rate equations for the homogenized species concentrations and with the volume reactivities λ = c −d λ : This system of ordinary differential equations yields non-linear oscillations around a neutral fixed-line which is determined by the initial conditions. The fixed-line steadystate concentrations can be obtained by setting the time derivatives to zero, resulting in with the conserved total population density ρ = a 1 + a 2 + a 3 = const, parameterizing the fixed-line. Linearization about this three-species coexistence fixed-line yields the stability matrix with eigenvalues {0, −iω 0 , iω 0 }. Since the non-zero eigenvalues are purely imaginary, the mean-field RPS system performs perpetual non-linear oscillations around the coexistence fixed point with frequency (in the linearized approximation) ω 0 = ρλ/ √ 3.

Doi-Peliti field theory and generalized Langevin equations
The bulk part of the Doi-Peliti action for the stochastic spatially extended RPS model follows directly from the reactions (1) and reads ‡ For convenience, here we drop all position and time indices ( x, t) on the fields and identify a 4 = a 1 . The first term describes the random nearest-neighbor hopping of the particles in the system, while the second contribution originates from the nonlinear reactions (1). As the auxiliary fieldâ i ( x, t) corresponds to a projection dual state, with average â i ( x, t) = 1, a Doi shiftã i ( x, t) =â i ( x, t) − 1 is conveniently applied to have the new field averaged to ã i ( x, t) = 0. After the Doi shift and ignoring boundary terms, the action becomes This shifted action may now be viewed as a Janssen-De Dominicis response functional [42,43] that represents the stochastic dynamics in terms of generalized Langevin equations. Theã i fields play the role of response fields and their coupling to the particle densities, shown in the terms that are second order in these fields, entails the presence of multiplicative noise terms. This comparison leads to the formulation of equivalent Langevin stochastic differential equations encoded in the action (6), where ζ i ( x, t) are the components of multiplicative noise in the system with vanishing means and correlations The resulting action becomes A RPS = A RPS 0 + A RPS int , with the Gaussian part and the nonlinear contributions (vertices) Here, ω 0 = λρ/ √ 3 denotes the zeroth-order oscillation frequency of the φ +/− modes. A RPS 0 is the diagonalized harmonic part of the action, while A RPS int represents the "interaction" contributions for the perturbative expansion. We omit the explicit expressions for the four-point vertices as they will not contribute to the dispersion relations at one-loop order, which shall be clear in the calculation below. It is manifest that φ o = (a 1 + a 2 + a 3 − ρ)/ √ 3 represents the fluctuation of the total particle number density. Due to the conservation law (13), the φ o mode is purely diffusive, and its exact dispersion relation in the harmonic part of the action acquires no fluctuation corrections. The φ + and φ − modes may be viewed as the left-and right-rotating waves in the system. At tree level, they are purely oscillating modes without dissipation, i.e., the real part of the mass term vanishes.

One-loop fluctuation corrections
We have applied a field-theoretical perturbation theory to one-loop order and calculated the renormalized diffusion constant D r , damping constant γ r , and oscillation frequency ω r §. To all orders in the fluctuation expansion extending beyond the mean-field approximation, there should be no correction to the two-point vertex function Γ (1,1) φoφo or propagator self-energy for the φ o mode, whence it retains its tree-level purely diffusive dispersion relation as dictated by the conservation law. For the φ ± modes, the one-loop Feynman diagrams for the two-point vertex functions Γ Fig. 2. The solid lines represent the φ + and φ − propagators, while the dashed lines indicate the diffusive mode φ o . In our convention, time and hence momentum always flow from right to left in the Feynman diagrams. The analytic expression for the two-point vertex function Γ where k is short-hand for the d-dimensional wavevector integral d d k/(2π) d , and the function I is defined as The damping constant u 0 in Eq. (20) is introduced to regularize the infrared (IR) singularities that emerge in later calculations. A nonzero renormalized u r will be generated by the fluctuations, but it is of higher order in the perturbation expansion; thus, we need to take u 0 → 0 at the end. This two-point function can also be expressed with the renormalized parameters as where Z φ ± absorbs all related wave function renormalizations (ultraviolet / UV divergences) in (20). The renormalized diffusivity D r , damping u r , and oscillation frequency ω r can be inferred accordingly from the explicit one-loop result (20) and (22). We obtain the following formal expressions for D r , u r , and ω r , Hence we indeed see that a non-zero damping u r is generated at one-loop order from the fluctuations, in agreement with Monte Carlo simulation data [14]. The infra-red (IR) divergence at one-loop order that appears when in the renormalized oscillation frequency ω r in low dimensions d ≤ 2 is caused by the contribution of the massless φ o mode as u 0 is set to zero. The infra-red (IR) divergence at one-loop order that appears in the renormalized oscillation frequency ω r in low dimensions d ≤ 2 is caused by the superposition of φ + and φ − modes as u 0 is set to zero. Our analysis of the one-loop results shows that the φ ± modes are inherently massive, as they acquire non-zero damping u r . Thus, the IR divergence can be resolved simply by maintaining a finite value for u 0 . For dimensions d ≥ 2, there are also UV divergences present. Nevertheless, all systems of interest have a natural cutoff in the UV limit, which is defined by the lattice constant c. The renormalized variables in different physically accessible dimensions are presented below.
2.5.1. d = 1: In one dimension, the renormalized parameters are , some numerical results are depicted in Fig. 3. We observe that the renormalized frequency ω r diverges when u 0 → 0. This IR divergence indicates strong fluctuation corrections to the oscillation frequency in the perturbative regime where reaction rates are small, and u 0 → u r ∼ λ is of first order in the reactivity. However, numerical simulations are invariably performed outside this regime for the sake of computational efficiency, as large reaction rates are used to avoid long relaxation times. Thus, no strong  fluctuation-induced renormalization have been encountered in the simulations. To oneloop order, u r > 0, which indicates the stability of the system's spatially homogeneous ground state with respect to fluctuations.

d = 2:
In two dimensions, the renormalized variables read Note that we have explicitly introduced the cutoff Λ ∼ π/c to regularize the UV divergence for the renormalized oscillation frequency. We plot ω r and the damping parameter u r in Fig. 4. The renormalized diffusion constant D r only linearly depends on the parameter λ. As in the one-dimensional case, the oscillation frequency ω r diverges as u 0 → 0. We note that the damping constant u r is also positive at d = 2.

d = 3:
In three dimensions, we may safely set u 0 = 0 and the renormalized system parameters follow as , We notice that for d > 2 the IR divergences disappear and fluctuation effects become generally weak. u r is also positive for d = 3, as displayed in Fig. 5.
We have found that in dimensions d = 1, 2, and 3, the diffusivity D experiences an upward shift, indicating that fluctuations enhance the diffusion. Our results, depicted in Fig. 3, 4, and 5, show that as the reaction rate increases, the characteristic frequency ω r shifts to smaller values, as the reactions drive the system towards a spatially more homogeneous distribution, leading to slower oscillations. The decline 0.
2. 4. 6. 8. 10. in the characteristic frequencies is in accordance with the numerical simulation data in Ref. [14]. In contrast with the LV model, Monte Carlo simulations of the RPS system do not appear to show strong renormalization effects [14], although both models feature a logarithmic divergence in two dimensions. The IR divergence in the RPS model appears as a consequence of the fact that the corrections are built using the Gaussian theory which has zero damping, precisely as in the LV model. The positive fluctuation-induced damping µ r , in contrast to the possibly negative one in the LV model, indicates that the system remains stable against the spontaneous emergence of spatio-temporal structures.

ML model and mean-field rate equations
The following discussion of the mean-field theory, Doi-Peliti action, and Langevin representation for the spatially extended stochastic ML model was laid out in detail in Ref. [20]. Here we summarize the pertinent points needed for our comparison with the RPS model and the computation of the fluctuation corrections to one-loop order. Following the conventions in Ref. [20], the reactions in the ML model read where i = 1, 2, 3 denotes the three competing species, and we identify B 4 = B 1 as before.
In contrast to the RPS system, the reactions in the ML model do not conserve the total particle number. The first two reactions implement predation and reproduction independently, while the third reaction implements "soft" site occupation constraint to effectively represent a finite carrying capacity. As in the RPS model, we consider a model wherein particles from all three species perform random walks on a d-dimensional hyper-cubic lattice with L d sites and lattice constant c. In the large diffusivity limit, the dynamics is governed by the mean-field rate equations where σ = c d σ and κ = c d κ are the volume reaction rates. Instead of a fixed line defined by the initial condition, the ML system displays a unique fixed point at mean-field level. By setting the time derivatives to be 0, the steady-state concentrations are found to be and the associated stability matrix reads Its eigenvalues at the coexistence fixed point are {−µ, −µ(2κ−σ±i √ 3σ)/2(σ+κ)}. The first eigenvalue −µ is always negative which implies the stability of the corresponding eigenmode, namely the exponential decay of the total particle number, see below. The imaginary part of the two complex conjugate eigenvalues, ± √ 3µσ/2(σ + κ), represents the frequency of temporal oscillations for the associate modes, whose amplitudes are either exponentially damped or growing. When 2κ > σ, the real part of the complex eigenvalues is negative and the limit circles are stable. Otherwise, for 2κ < σ, the limit circles are unstable and one observes the spontaneous formation of spiral structures in the system. In the vicinity of the Hopf bifurcation at 2κ = σ, the time evolution of the two modes corresponding to the complex conjugate eigenvalues becomes much slower than the fast relaxing mode, which introduces a natural time scale separation. As a consequence of the critical slowing down near the Hopf bifurcation, the fast relaxing mode can be integrated out and the system is effectively governed by the complex timedependent Ginzburg-Landau equation [20].

Doi-Peliti field theory and generalized Langevin equations
The Doi-Peliti action follows from the reactions of the ML model and reads This action does not obey the U (1) global symmetry present in the RPS model; indeed, the total particle number is not conserved under the ML reactions (27). Following the Doi shift to the fluctuating auxiliary fieldsb i ( x, t) =b i ( x, t) − 1, the action becomes As in the RPS model above, we may now view this shifted action as a Janssen-De Dominicis functional which is equivalent to the corresponding generalized Langevin equations where ξ i ( x, t) represent the multiplicative noise components with correlators

Diagonalization of the harmonic action
Before diagonalizing the quadratic action, we first shift to fluctuating fields, Here C is a counterterm which encodes the fluctuation corrections to the average concentrations. Owing to the cyclic symmetry among the three different species, we only need to introduce a single counterterm. The harmonic part of the action in terms of the new fieldsd i and d i reads (36) Since the counterterm C is of first order in the perturbative expansion parameters, up to zeroth order the harmonic part of the action can be diagonalized by the following linear transformation and    After applying this linear transformation, the action can be expressed as representing, respectively, the harmonic, source, and non-linear interaction terms. Again, we omit the four-point vertices, since they will not contribute to the dispersion relation renormalizations at one-loop order. The other terms are The ψ o mode corresponds to the fluctuation of the total particle density and decays exponentially at tree level. In contrast to the RPS model, the ML ψ ± modes display non-vanishing dissipation γ 0 already on the mean-field level. As mentioned above, a Hopf bifurcation occurs at 2κ = σ; when 2κ < σ, the system is unstable and spiral structures are spontaneously generated. In the perturbative regime, the assumed tiny fluctuation corrections should not change the overall stability features, but will only shift the Hopf bifurcation point by a small amount.

One-loop fluctuation corrections
Prior to calculating the renormalized quantities, we need to determine the counterterm C by requiring the average fluctuation of the total density to be zero, ψ o = 0. Up to one-loop order, the contributions to ψ o are shown in Fig. 6. We note that the second diagram in Fig. 6 is of second order and thus can be dropped. The corresponding analytic expression results in We may now proceed to the fluctuation renormalization of the two-point vertex functions Γ (1,1) , which encode the self-energies entering the dispersion relation of the ψ o and ψ ± modes. As the total particle number is not conserved in the ML model, the dispersion relation acquires non-trivial corrections from the perturbation expansion. The one-loop diagrams that contribute to the vertex function Γ (1,1) ψoψo are pictured in Fig. 7. The last diagram is of higher order and thus can be omitted; this Provided γ 0 > 0, all corrections at one-loop level are real, and no imaginary part appears in the mass term of the ψ o mode, hence there are no total particle number oscillations. However, for γ 0 < 0 the system exhibits emergent oscillations of the total particle number, indicating the spontaneous formation of spatio-temporal structures. As the renormalized two-point vertex function can also be written as the renormalized diffusivity D o r and mass parameter µ r can be calculated accordingly. Here, Z φo absorbs all wave function renormalizations. Since the rotating wave modes ψ ± acquires a different diffusivity renormalization from the o mode, we carefully distinguish the renormalized quantities D ± r and D o r . The explicit expressions for D o r and µ r read After evaluating the integrals one arrives at (46) For d ≥ 2, UV divergences in µ r are manifest; but since the lattice constant c serves as a natural UV cutoff in lattice models, we will not discuss these UV divergences further.
The renormalized vertex function Γ can also be calculated according to the one-loop diagrams in Fig. 8, resulting in Upon invoking the relation with renormalized quantities the expressions for the renormalized parameters γ r and ν r are readily inferred, where the coefficients M (±) i are defined in Eq. (B.2) in the appendix, and the angles θ and η are given by tan θ = ν 0 /γ 0 and tan η = ν 0 /(γ 0 + µ). We note that at odd dimensions d, the first term in the bracket in Eq. (49) switches from real to imaginary as γ 0 changes its sign; however, at even dimensions, this term is always real. Finally, the renormalized diffusivity reads In the appendices, we provide additional details and the definitions of the various coefficients, as well as explicit evaluations for d = 1, 2, 3. Here, we focus on the behavior of the damping parameter µ r across different dimensions. It is important to note that µ r is generally a complex number: Its imaginary part conveys information about spatial oscillations, while its real part represents either exponential decay or growth of the average particle density.

d = 2:
At two dimensions, we need to introduce the UV cutoff Λ ∼ π/c; the damping parameter becomes As in the one-dimensional case, when γ 0 > 0, we have pure damping, whereas the system displays population oscillations if γ 0 < 0. The damping coefficients µ r at d = 2 for different bare parameters are plotted in Fig. 10.

d = 3:
In three dimensions, the damping parameter reads Again, oscillations emerge in the region where γ 0 < 0. The real and imaginary parts of the damping parameter µ are depicted in Fig. 11. For ν 0 = 1, the damping µ r decreases in one and two dimensions, but increases in three dimensions in the stable region where γ 0 > 0, as shown in Figs. 9, 10, and 11. This indicates that the reactions effectively slow down the relaxation processes in one and two dimensions, but speed them up in three dimensions. However, due to the highly nonlinear dependence on the parameters µ 0 , γ 0 , and ν 0 , a more general conclusion cannot be made. It is worth noting that in the unstable region, an imaginary part of µ r is generated, formally resulting from the subtraction of an inadequate homogeneous steady state. Yet these emergent oscillations manifestly indicate the instability with respect to spontaneous formation of spatio-temporal patterns in this regime. We emphasize again that the one-loop fluctuation corrections should merely induce small quantitative corrections, and cannot induce qualitative changes in the region where perturbation theory is applicable. The ML system thus maintains a bifurcation point, below which the homogeneous ground state is rendered unstable and spiral structures emerge.

Comparison between the RPS, ML, and LV models
In this section, we present a thorough comparison between the spatically extended stochastic LV, RPS, and ML models. We specifically discuss the spontaneous formation of spatio-temporal structures and the stability of the homogeneous state up to perturbative one-loop order in the fluctuation corrections. We also briefly address the influence of quenched spatial disorder in the reaction rates in the perturbative regime.

Spiral formation from a single lattice site point of view
It is well-established that in sufficiently large spatial systems, the ML model displays spontaneously emerging dynamic spiral structures in individual simulation runs (these are of course averaged out in ensemble averages). Here, we propose that one crucial necessary condition for the formation of such persistent spatio-temporal patterns is the existence of a stable and uniform oscillation frequency at the local lattice site level, which then allows spatially extended coherent oscillatory features.
Each site in a lattice subject to stochastic reactions and spreading processes can be regarded as a separate system that is coupled to a particle reservoir (in the thermodynamic limit). In the ML model, at the (linearized) mean-field level, the local oscillation frequency is uniquely determined by the reaction rates, as is also apparent in the diagonalized ML Doi-Peliti action (39) at tree level. A similar definition of the oscillation frequency at linearized mean field level can also be obtained in the LV model [39]. However, in the RPS model at tree level, the oscillation frequency is set by the global conserved particle number ρ. As the particle number at each site is changing all the time owing to its coupling to the environment that serves as a nonlocal reservoir, there does not exist a unique characteristic oscillation frequency for each site during any single run stochastic realization. Ultimately, these nonlocal effects originate from the long-range correlations introduced by the global particle number conservation law, whose relevance in the context of pattern formation was demonstrated in Ref. [11]. The average of the oscillation frequency with a fixed total particle number is given by the expression in the diagonalized RPS Doi-Peliti action (18).
These straightforward tree-level arguments are readily generalized to all (loop) orders in the perturbation expansion. In the ML model, no nonlinear terms that depend on the total particle density (or any other global quantity) are present in the action, and thus the renormalized frequency will also be independent of the particle number density. This is not the case for the RPS model, where ρ manifestly enters the vertices. A perhaps more straightforward way to understand this distinction invokes the fixed points (stationary species densities) of the systems. In the RPS system, there exists a fixed line that is parameterized by the conserved total particle number. For each lattice site in the RPS model, the population oscillations wander along this fixed line with a mean that corresponds to the averaged local total particle number at this specific site. In contrast, the ML and LV models have unique fixed points, independent of any global constraints that originate from conservation laws, and consequently all lattice sites are governed by identical oscillation frequencies determined by these fixed points.
To illustrate these points, we plot the time evolution of the population density of a single species at two randomly chosen lattice sites in a single simulation run for the RPS and ML models in Figs. 12 and 13, respectively. For the RPS species density, the time intervals separating different peaks are not regularly spaced for the two lattice sites. Correspondingly, its Fourier transform displays multiple peaks with almost equal intensities, and no well-defined characteristic frequency is discernible. In contrast, the peaks in the ML species density time evolution are evenly separated, and hence a dominant Fourier peak emerges. These numerical observations support our analysis that the frequency in a single run in the ML model is well-defined, while that is clearly not the case for the RPS system. Yet a well-defined oscillation frequency clearly constitutes a necessary condition to form spatio-temporal patterns.
Although the arguments in this section are formulated in the framework of individual sites in a regular lattice, they are readily extended to an effective unit cell that is similar in size to the characteristic diffusion length scale, or to models defined on a continuum which require a finite reaction range. As long as the system size is much larger than any of these small length scales (that provide suitable ultraviolet cutoffs for the continuum field theory), the specific microscopic details implemented in numerical simulations are not expected to have a significant impact on our conclusions.

Stability of the homogeneous ground states against fluctuations (to one-loop order)
While we posit above that a stable oscillation frequency at each lattice site is a necessary requirement to form oscillatory spatio-temporal patterns such as spirals, this is not a sufficient condition. Another crucial ingredient is of course that the spatially homogeneous stationary state is rendered unstable; fluctuations of a certain wavevector range will then spontaneously generate spatial patterns [45,46].
In the stochastic spatially extended LV model [39], in the absence of site restrictions, the population oscillation damping vanishes in the Gaussian approximation. To oneloop order, a negative damping frequency is generated which indicates an instability of the spatially uniform system, inducing nontrivial particle transport that in turn drives the formation of expanding evasion-pursuit activity waves. In the ML model, the damping coefficient is already nonzero in the Gaussian approximation. As stated above, when 2κ < σ, the limit cycles of the ML model become unstable and spiral structures emerge. Yet within the realm of a perturbative approach, any fluctuation corrections will only shift the Hopf bifurcation point, but cannot qualitatively change the system's overall features. In contrast, the RPS system behaves similarly to the LV system, as we encounter a vanishing damping coefficient in the Gaussian approximation. However, at one-loop level, the generated mass term is positive indicating merely a generated finite damping constant, as opposed to its negative counterpart in the LV model. Hence the RPS model remains inert against spontaneous pattern formation. Higher-order terms in the fluctuation expansion should not overturn the sign of the one-loop results within the perturbative regime, which renders this argument perturbatively robust.

Influence of quenched spatially disordered reaction rates
In this part, we briefly investigate the effect of spatially disordered, uniformly distributed quenched reaction probabilities on the fluctuation corrections. We start with the Doi-Peliti action where L 0 and L are (usually polynomial) functions of the fields {â i } and {a i }. We introduce spatial disorder to L , which we take to represent stochastic reactions while L 0 encodes particle diffusions. The action with quenched spatial disorder is then given by here, we assume η( x) to be uniformly distributed in the interval [0, 2] with mean η = 1.
Note that the reaction rates are only spatially disordered and remain fixed in time. The average of any observables follows from Eq. (A.9) through where the overline denotes the quenched disorder average. We can readily integrate out the disorder η( x) first, since the observable O is independent of the random variables η, and arrive at with Incorporating disorder in the original Doi-Peliti action (54) thus effectively leads to an additional term that is nonlocal in the time domain, owing to the temporally fixed reaction rates. However, in the perturbative regime where L is small, this extra term can be expanded near L = 0: where the labels t and t distinguish the different time dependences. It is evident from this expansion that the extra temporally nonlocal term "entangles" different replicas of the system. However, as the first term is already second-order in L , it is of higher order than the original action A, and should not markedly affect the system's fluctuation corrections in the naive perturbation regime, at least to low loop orders. This observation may account for the insensitivity of the RPS and ML models to quenched randomness in the reaction rates, as reported in Refs. [14,17]. We remark that the stronger effects of varying predation rates λ in the LV system can be largely traced to the sensitivity ∼ 1/λ of the stationary species densities already on the mean-field approximation level [35].

Summary and conclusions
In this paper, we have investigated the dynamics of the RPS and ML models up to one-loop order in the fluctuation corrections by means of a perturbative field-theoretical analysis. We utilized the Doi-Peliti formalism to obtain the dynamical probability functional for the stochastic Markovian dynamics, and also extracted the equivalent generalized Langevin equations. In the Gaussian theory, as expected, the RPS model displays only purely oscillatory modes in addition to the strictly diffusive conserved total particle density. In the ML model, a Hopf bifurcation point appears at 2κ − σ = 0 that separates the parameter space into stable and unstable regions; in the latter regime, spiral structures are spontaneously generated. The one-loop fluctuation corrections in the RPS model, which are of first order in the effective nonlinear coupling λ/D 1, generate dissipation. We have found that in the physically accessible dimensions d = 1, 2, 3, the damping coefficient is always positive. This indicates the stability of the spatially uniform stationary state in the RPS model, at least in the perturbative regime. Thus, our analysis sheds light on the absence of spatio-temporal structures in the RPS model. In addition, the oneloop correction to the oscillation frequency is IR divergent due to the dissipation-free nature of the mean-field modes, very similar to the LV model. However, outside the range of validity of the perturbation expansion, the damping terms become sizeable; hence we argue that this IR divergence becomes naturally regularized. This explains why no significant fluctuation corrections to the oscillation frequency have as yet been numerically observed in computer simulations, which have invariably been situated far away from the perturbative regime.
In the ML model, as both the dissipation and oscillation frequencies are already finite at mean-field level, the one-loop fluctuation corrections should not qualitatively modify the mean-field conclusions. Since both propagating modes are massive due to the finite damping coefficients, the ML system does not display any IR singularities. Hence the ML model is insensitive to fluctuations at least perturbatively, and in the region where the homogeneous steady state is stable. Moreover, we have argued that uniformly distributed quenched random disorder in the reaction rates only weakly influences fluctuation corrections in either system, which is in agreement with earlier Monte Carlo simulation data.
Finally, we provided two decisive criteria that determine the possibility of spatiotemporal structures in the LV, RPS, and ML models. The first argument considers a single lattice site point of view, while the second is based on studying the global stability of the spatially homogeneous stationary state of the system. From a single lattice site perspective, a necessary (but not sufficient) condition for the emergence of spatially extended coherent oscillatory behavior is that the oscillation frequency is constant over space and time in each run. Different oscillation frequencies on distinct sites would not allow the formation of stable coherent patterns. From a global point of view, only if the spatially uniform quasi-steady state is unstable against finite-wavelength fluctuations can non-trivial spatio-temporal structures be generated, as is evident in the ML model. Both these criteria explain the absence of spatio-temporal patterns in the RPS system, as a consequence of the relevant conservation law for the total particle number, in contrast with the otherwise apparently similar LV and ML models. We remark that adding some external noise to the RPS model that explicitly invalidates total particle number conservation might induce the formation of spiral patterns at sufficiently large length and long time scales.

Appendix A. Doi-Peliti formalism and the asymmetric RPS model
In this appendix, we provide a brief overview of the construction of the time evolution operator and the resulting field-theoretic action via the Doi-Peliti mapping of stochastic reactions to a non-Hermitean many-body quantum action [47,48,49,50] (for recent reviews and additional details, see Refs. [40,41,51,52]). We then proceed to construct coupled Langevin equations describing the dynamics of the system similar to previous work in the context of the LV system [32,39,40], plankton-based predator-prey models [53], and Turing patterns [54]. An interesting corner limit is analyzed via the generalized Langevin equations; we will show that the strongly asymmetric RPS system reduces to an effective LV model. Finally, we construct the diagonalized action for the RPS model and briefly compare the general situation with the symmetric version discussed in the bulk of this paper. This appendix is written in a self-consistent way and we hope it will be of use to readers who would like to delve into the Doi-Peliti formalism.

Appendix A.1. Stochastic time evolution operator
The RPS rules (1) mandate that particle numbers are discrete, hence the occupation numbers of lattice sites can be written as positive integers n iα , where the index i accounts for different particle species, while α enumerates the sites on a d-dimensional hyper-cubic lattice. The master equation for the local, on-site RPS reactions then reads: λ i (n iα − 1)(n i+1α + 1)P (n iα − 1, n i+1α + 1; t) −n iα n i+1α P (n iα , n i+1α ; t) , where the index i wraps around (i.e., i = 4 is to be identified with i = 1). Note that for brevity we have not included hopping to adjacent lattice sites here. As an initial state, we assume a uniform distribution of particles with an average initial number of particles N i per lattice site of species i. This corresponds to a Poisson distribution for the occupation number of each species i at all lattice sites α, P (n iα ; t) = i=1,2,3 N n iα i e −N i /(n iα !). The discrete nature of the possible states of the RPS systems suggests the introduction of a product Fock space state vector where the |n iα represent the occupation states of species i on lattice site i. In analogy with the quantum-mechanical harmonic oscillator, the single-site states (and thereby the full state vector) can be acted upon by bosonic ladder operators obeying the commutation relations [a iα , a jβ ] = 0 and [a iα , a † jβ ] = δ αβ δ ij . The occupation number eigenstates are constructed via a iα |n iα = n iα |n iα − 1 , a † iα |n iα = |n iα + 1 , and the empty state |0 is defined by a iα |0 = 0. where H denotes the (time-independent) Liouville operator which can be split into a diffusion and a reaction term, H = H diff + H reac , where the on-site reaction contribution is a sum of local terms H reac = L d i=1 H α , and specifically for the RPS model Similarly, since on-lattice diffusion is implemented by particles performing simple jumps between nearest-neighbor lattice sites, the diffusion part of H reads where αβ indicates a sum over all possible nearest-neighbor lattice site pairs in the system.
Appendix A.2. Coherent-state path integral and equivalent Langevin partial differential equations Following the steps of Refs. [39,40,41], we write averages for observables O = O({n iα }) as a multi-dimensional integral over coherent states where the ψ iα and ψ * iα are complex eigenvalues describing the coherent right and left eigenstates of the ladder operators a iα and a † iα , respectively. The coherent-state "action" is given by where we have to replace the ladder operators by their eigenvalues in the Liouville operator H.
In the spatial continuum limit (lattice constant c → 0) we may replace the sum over lattice sites with a d-dimensional volume integral t). Hence, the "bulk" part of action (not considering the terms from the initial conditions at t = 0 and the projection states at t = t ) of the RPS system is given by (A.8) In the continuum limit we can thus write averages in the following coherent-state path integral form Our aim here is to derive stochastic partial differential (Langevin) equations for the species concentrations that accurately capture the intrinsic reaction noise. To this end, we note that the Janssen-de Dominicis response functional [40,42,43,55] is equivalent to the set of SPDEs with the associated noise (cross-)correlations This correspondence allows the immediate derivation of a coupled Langevin equation formulation of any system that exhibits an action functional of the form (A.10). Hence, via a direct comparison with the action of the RPS system (A.7), we can extract the deterministic part of the SPDEs describing the RPS system which equal the right-hand side of the mean-field equations, as they should. Furthermore, the effective noise correlations are given by the matrix L ij : (A.14) Hence, the SPDEs (A.11) can be constructed from the mean-field equations by including a term that accounts for diffusion and multiplicative noise terms obeying the given (cross-)correlations. Note that the noise auto-correlations L ii are always determined by the concentration of the predator species A i and its respective prey A i+1 , and the scale is set by the associated predation rate λ i . Thus, the noise directly associated with a given species is solely determined by its role as predator.
Appendix A.3. Strongly asymmetric RPS model: mapping to the LV system In order to investigate the asymmetric "corner" limit of the RPS system, we re-define the interaction rates as λ 1 = λ/x, λ 2 = λ and λ 3 = κλ. The dimensionless variable x varies in the interval (0, 1] and describes the asymmetry of the rates, while the equally dimensionless parameter κ is of order unity and describes the difference between the predation rates of species A 2 and A 3 . We are interested in the limit x → 0 in which the predation reactions between species A 1 and A 2 dominate. The concentrations at the coexistence fixed point become Hence, the densities of species A 1 and A 2 become small as x → 0, while species A 3 makes up most of the overall species abundance. This is the "corner" limit in which RPS can be approximated by a two-species Lotka-Volterra system, with the third, most abundant species serving as a mean-field like reservoir to feed the first species, and to provide the effective spontaneous death reaction for the second species, as explained above [56]. Species A 1 thus effectively turns into prey, while A 2 becomes the sole predator species. The noise correlation matrix L in this corner case reads The noise strength of species A 1 , as well as the noise cross-correlations between species A 1 and A 2 , are inversely proportional to the large rate scaling factor x, indicating that fluctuations of species A 1 and A 2 (the LV predator and prey, respectively) become strong in the limit x → 0. Indeed, writing the resulting effective SPDEs in the limit of large λ 1 and assuming a homogeneous and stationary distribution of species A 3 yields with the noise correlations This set of Langevin equations precisely matches those derived directly for the LV model [39].

Appendix A.4. Fluctuation corrections
In order to gain more insight into the role of fluctuations in the RPS system, we study the non-linear vertices arising from the Doi-Peliti action (A.7). To this end, we first need to diagonalize the action by transforming to appropriate field combinations. We then list the resulting vertices that capture fluctuation corrections beyond the Gaussian mean-field approximation.
Our goal is to find a transformation that diagonalizes the mass matrix A, and thus the harmonic part of the action (A.7), if we set all diffusivities equal, D i = D. The matrix A is asymmetric, hence we make use of its orthogonal left and right eigenvectors u i A = e i u i and A v i = e i v i , respectively. The resulting eigenvector matrices Q = ( u 1 , u 2 , u 3 ) and where Λ = λ 1 λ 2 λ 3 (λ 1 + λ 2 + λ 3 ), and (A.26) The right and left eigenvector matrices then transform the mass matrix to the diagonal form Q T AP (Q T P ) −1 = diag(e i ). Defining new fieldsφ i and φ i according to the transformationc i = j Q T jiφ j and c i = j P ij φ j , we arrive at It is already obvious from this structure that the fieldsψ and ψ describe the fluctuation of the total population, while the other fields are oscillatory in nature. Employing these Note that we have employed a different diagonalization convention in this appendix as compared to the main text. This difference is reflected in the constant factors in the propagators.

transformations, one arrives at the diagonalized harmonic action
The field ψ = λ 2 (c 1 + c 2 + c 3 )/(λ 1 + λ 2 + λ 3 ) is massless, encodes no reactions, and is purely diffusive, as it represents the total local concentration of all three species. The corresponding harmonic propagator is given in Fourier space by while the oscillating field propagators display poles at finite eigenfrequencies ∓ω 0 , similar to the LV case [39]. The action transformed to the new fields becomes quite cumbersome to write out in full, hence we merely provide the coefficients of the possible field combinations in the vertices in Table A1. Note that we omit the coefficients of any four-point vertices as these do not contribute to one-loop order corrections.

Appendix B. Detailed calculations for the ML model
Here we provide intermediate steps for the one-loop calculation of the ML model. The renormalized frequencies are 3ν 0 )(µ 2 + 3ν 2 0 ) − ν 2 0 (µ 2 + ν 2 0 ) , The renormalized diffusivity is and can be further simplified to D ± r =D + κ + σ d P where P (±) i = ReP i ± iImP i and Q (±) i = ReQ i ± iImQ i , with , .

(B.5)
Appendix C. Renormalized variables in the ML model at physically accessible dimensions Finally, we list the expressions of the renormalized variables at physically accessible dimensions d = 1, 2, and 3.