A stochastic analysis of the spatially extended May-Leonard model

TL;DR

This paper derives a coarse-grained stochastic evolution equation for the May-Leonard model, capturing spatial pattern formation and fluctuations near the Hopf bifurcation, using advanced field-theoretic methods.

Contribution

It extends deterministic models by incorporating stochastic effects, deriving a complex Ginzburg-Landau equation for the full stochastic dynamics of cyclically competing species.

Findings

Derivation of the complex Ginzburg-Landau equation from stochastic May-Leonard dynamics.

Identification of parameter constraints for time scale separation.

Reduction to two slow degrees of freedom in the model.

Abstract

Numerical studies of the May-Leonard model for cyclically competing species exhibit spontaneous spatial structures in the form of spirals. It is desirable to obtain a simple coarse-grained evolution equation describing spatio-temporal pattern formation in such spatially extended stochastic population dynamics models. Extending earlier work on the corresponding deterministic system, we derive the complex Ginzburg-Landau equation as the effective representation of the fully stochastic dynamics of this paradigmatic model for cyclic dominance near its Hopf bifurcation, and for small fluctuations in the three-species coexistence regime. The internal stochastic reaction noise is accounted for through the Doi-Peliti coherent-state path integral formalism, and subsequent mapping to three coupled non-linear Langevin equations. This analysis provides constraints on the model parameters that allow time scale separation and in consequence a further reduction to just two coarse-grained slow degrees of freedom.

Equations

$$X_{\alpha }+X_{\alpha +1}\to X_{\alpha }$$

$$X_{\alpha }\to X_{\alpha }+X_{\alpha }$$

$$X_{\alpha }+X_{\alpha }\to X_{\alpha }\mathrm{with}\mathrm{rate}{\lambda }^{\prime }.$$

$$X_{\alpha ,i}\to X_{\alpha ,i+1}\mathrm{with}\mathrm{rate}{D}^{\prime },$$

$$\frac{d{a}_{\alpha +1}(t)}{dt}=\mu a_{\alpha +1}(t)-\lambda a_{\alpha +1}(t{)}^2-\sigma a_{\alpha }(t)a_{\alpha +1}(t),$$

$$L=\frac{-\mu}{\sigma+\lambda}\ \left(\begin{array}[]{ccc}\lambda&0&\sigma\\ \sigma&\lambda&0\\ 0&\sigma&\lambda\end{array}\right)\ .$$

$${\nu }_0=-\mu \mathrm{and}\{\nu ,{\nu }^{\ast }\}=\frac{-\mu }{\sigma +\lambda }\left(\lambda -\sigma e^{\mp i\pi /3}\right);$$

$$\displaystyle\frac{\partial P(n_{\alpha,i};t)}{\partial t}=\sum_{\alpha=1,2,3}\biggl{(}\mu\,\Bigl{[}(n_{\alpha,i}-1)\,P(n_{\alpha,i}-1;t)-n_{\alpha,i}\,P(n_{\alpha,i};t)\Bigr{]}$$

$$\displaystyle\quad\!+\sigma^{\prime}\,\Bigl{[}n_{\alpha,i}(n_{\alpha+1,i}+1)\,P(n_{\alpha,i},n_{\alpha+1,i}+1;t)-n_{\alpha,i}\,n_{\alpha+1,i}\,P(n_{\alpha,i},n_{\alpha+1,i};t)\Bigr{]}$$

$$\displaystyle\qquad\qquad\qquad\qquad+\lambda^{\prime}\,\Bigl{[}(n_{\alpha,i}+1)\,P(n_{\alpha,i}+1;t)-n_{\alpha,i}\,P(n_{\alpha,i};t)\Bigr{]}\biggr{)}\ .$$

$$P(n_{\alpha ,i};0)=\prod _{\alpha =1,2,3}\frac{{\bar{n}}_{\alpha }^{n_{\alpha ,i}}}{n_{\alpha ,i}!}{e}^{-{\bar{n}}_{\alpha }}.$$

$$|\Phi (t)⟩=\sum _{\{n_{\alpha ,i}\}}P(n_{\alpha ,i};t)|n_{\alpha ,i}⟩.$$

$$\frac{\partial |\Phi (t)⟩}{\partial t}=-H|\Phi (t)⟩,$$

$$H_{\mathrm{reac}}$$

$$H_{\mathrm{diff}}=\sum _{⟨i,j⟩}\sum _{\alpha =1,2,3}{D}^{\prime }(a_{\alpha ,i}^{†}-a_{\alpha ,j}^{†})(a_{\alpha ,i}-a_{\alpha ,j}),$$

$$⟨\mathcal{O}⟩\propto \int \prod _i^N\prod _{\alpha =1,2,3}\mathcal{D}[{\psi }_{\alpha ,i}^{\ast },{\psi }_{\alpha ,i}]\mathcal{O}(\{{\psi }_{\alpha ,i}\})e^{-S[{\psi }_{\alpha ,i}^{\ast },{\psi }_{\alpha ,i};t]},$$

$$S[{\psi }_{\alpha ,i}^{\ast },{\psi }_{\alpha ,i};t]$$

$$S[\tilde{a}_{\alpha},a_{\alpha};t]=\int dt\int d^{d}x\,\Biggl{[}\ \sum_{\alpha=1,2,3}\tilde{a}_{\alpha}(\partial_{t}-D\nabla^{2})\,a_{\alpha}+H_{\rm reac}(\tilde{a}_{\alpha},a_{\alpha})\Biggr{]}\ .$$

$$H_{\mathrm{reac}}({\tilde{a}}_{\alpha },a_{\alpha })$$

$$\frac{\partial \mathbf{a}(\vec{x},t)}{\partial t}=D{\nabla }^2\mathbf{a}(\vec{x},t)+\mathbf{F}[\mathbf{a}(\vec{x},t)]+\mathbf{\zeta }(\vec{x},t),$$

$$⟨{\zeta }_{\alpha }(\vec{x},t){\zeta }_{\beta }({\vec{x}}^{\prime },t^{\prime })⟩=2L_{\alpha \beta }[\mathbf{a}(\vec{x},t)]\delta (\vec{x}-{\vec{x}}^{\prime })\delta (t-t^{\prime }).$$

$$\displaystyle S[{\bf a}]=\int dt\int d^{d}x\sum_{\alpha}\tilde{a}_{\alpha}\biggr{[}(\partial_{t}-D\nabla^{2})\,a_{\alpha}-F_{\alpha}[{\bf a}]-\sum_{\beta}L_{\alpha\beta}[{\bf a}]\,\tilde{a}_{\beta}\biggr{]}\ .$$

$$\frac{\partial a_{\alpha +1}(\vec{x},t)}{\partial t}$$

$$\displaystyle L_{\alpha\beta}[{\bf a}]=\left(\begin{array}[]{ccc}\mu\,a_{1}-\lambda\,a_{1}^{2}&-\sigma\,a_{1}a_{2}/2&-\sigma\,a_{1}a_{3}/2\\ -\sigma\,a_{1}a_{2}/2&\mu\,a_{2}-\lambda\,a_{2}^{2}&-\sigma\,a_{2}a_{3}/2\\ -\sigma\,a_{1}a_{3}/2&-\sigma\,a_{2}a_{3}/2&\mu\,a_{3}-\lambda\,a_{3}^{2}\end{array}\right).$$

$$a_{\alpha }(\vec{x},t)=\frac{\mu }{\sigma +\lambda }+b_{\alpha }(\vec{x},t),{\tilde{a}}_{\alpha }(\vec{x},t)={\tilde{b}}_{\alpha }(\vec{x},t).$$

$$\frac{\partial \mathbf{b}(\vec{x},t)}{\partial t}=D{\nabla }^2\mathbf{b}(\vec{x},t)+\mathbf{f}[\mathbf{b}(\vec{x},t)]+{\mathbf{\zeta }}^{\prime }(\vec{x},t),$$

$$\mathcal{R}=\left(\begin{array}[]{rcr}1/\sqrt{2}&0&-1/\sqrt{2}\\ -1/\sqrt{6}&\sqrt{2/3}&-1/\sqrt{6}\\ 1/\sqrt{3}&1/\sqrt{3}&1/\sqrt{3}\end{array}\right).$$

$${\partial }_t\mathbf{c}=D{\nabla }^2\mathbf{c}+\mathcal{R}f[{\mathcal{R}}^{-1}\mathbf{c}]+\mathbf{\eta }$$

$${\partial }_t{c}_1$$

$${\partial }_t{c}_2$$

$${\partial }_t{c}_3$$

$$c_3=K\left(c_1^2+c_2^2\right),$$

$$K=\frac{(\sigma -2\lambda )(\sigma +\lambda )}{2\sqrt{3}\mu (2\sigma -\lambda )}{\left(1-\frac{D(\sigma +\lambda )}{\mu (2\sigma -\lambda )}{p}^2\right)}^{-1}.$$

Full text

1 Introduction

In the recent past, ecologists and applied mathematicians have sought to quantitatively delineate dynamic emerging phenomena such as biodiversity and population extinction [1, 2, 3, 4, 5]. The payoff for these theoretical endeavors potentially includes stabilizing and protecting endangered ecosystems in addition to establishing a fundamental understanding of the myriad forms of pattern formation observed in nature. Populations with a cyclic competition motif have been studied in several contexts such as ecology, epidemiology and opinion poll research. For example, cyclic dominance is observed in certain Californian lizard subspecies [6] and in petri dish experiments that involve E.-coli bacteria with three distinct strains [7]. Stochastic spatially extended population dynamics with cyclic motifs such as the May–Leonard model [8] have been investigated numerically, mostly on two-dimensional lattices [9, 10, 11]. In contrast to simpler cyclic rock-paper-scissors models for which the total particle number is conserved [12, 13, 14], May–Leonard systems display striking spiral patterns in certain parameter ranges. Linking the characteristic length and time scales in these spontaneously emerging spatio-temporal structures to the basic rates of the underlying stochastic processes is an important fundamental problem. In the framework of the deterministic rate equation time evolution, Reichenbach, Mobilia, and Frey [15] demonstrated that the formation of spirals in the May–Leonard model is effectively governed by the paradigmatic complex Ginzburg–Landau equation (CGLE).

The CGLE appears extensively in various contexts in physics [16], ranging from second-order phase transitions in condensed matter systems to string theory and ubiquitous non-equilibrium phenomena (see, for example, Refs. [17, 18, 19, 20, 21, 22, 23]). This complex partial differential equation typically describes a slowly varying continuous order parameter field in the presence of weak non-linearities near a bifurcation point that governs the instability of a spatially homogeneous state. It exhibits gauge invariance of the modulating variable under a global phase change, usually as a consequence of periodicity in the extended space-time. Intriguing more recent applications of the CGLE include the synchronization of coupled non-linear noisy oscillators [24] and driven-dissipative Bose–Einstein condensation (through an equivalent Gross–Pitaevskii equation with complex parameters) [25, 26, 27], underscoring the remarkable universality of the CGLE.

More than thirty years ago, Kuramoto demonstrated that a lattice of diffusively coupled oscillators is governed by a generic coarse-grained evolution equation near the Hopf bifurcation, namely the CGLE [17]. In 2007, Reichenbach, Mobilia, and Frey constructed the CGLE as a convenient effective description of the May–Leonard model near the Hopf bifurcation of the three-species coexistence fixed point in parameter space [15]. They showed that the emerging spiral wavelength and wavefront velocity are encoded in the coefficients of the CGLE near this fixed point. In this present paper, we establish a full derivation of the stochastic CGLE in this context, which properly accounts for intrinsic reaction noise, and hence extends the deterministic analysis of Ref. [15], and also the perturbative multi-scale expansion around the bifurcation performed by Szczesny, Mobilia, and Rucklidge [28]. We remark that the incorporation of intrinsic stochasticity is crucial, as in some prominent situations, spatio-temporal patterns cannot be adequately characterized by a mere deterministic treatment. This is true, for example, in stochastic spatially extended lattice Lotka–Volterra models for predator-prey competition and coexistence [29, 30, 31]. Under more general settings, Butler and Goldenfeld demonstrated that stochastic fluctuations may cause significant alterations to otherwise simpler deterministic patterns [32, 33].

In our derivation of the stochastic CGLE for the May–Leonard model, valid near its Hopf bifurcation and for small deviations from the stationary population densities in the three-species coexistence phase, we account for the systematic treatment of fluctuations due to internal reaction noise through a bosonic field-theoretic formalism. Thus, by including the inherent stochasticity of this system with cyclic species competition as encoded in the microscopic master equation for its defining reaction processes, then studying (small) non-linear fluctuations about the mean-field stationary densities in the three-species coexistence region, and finally exploiting time scale separation in the vicinity of the Hopf bifurcation that allows us to eliminate one fast relaxing mode, we arrive at a Langevin-type extension of the complex Ginzburg–Landau equation, for which random effects are superimposed on the non-linear deterministic behavior ascribed to the CGLE through (to leading order) additive noise terms. Based on this resulting effective dynamical theory, one could now, e.g., in a perturbative analysis akin to Ref. [30], evaluate fluctuation-induced renormalizations of the characteristic oscillation frequencies and attenuation, as well as typical spiral pattern wavelengths, and thereby quantitatively relate the stochastic CGLE components directly to numerical or actual observations in pattern formation. Yet our formalism more also yields a set of three coupled Langevin equations that faithfully describes the stochastic spatially extended May–Leonard system under quite general circumstances, not subject to the additional constraints required for the applicability of the CGLE.

In the following, we first describe and define our spatially extended stochastic version of the May–Leonard model, and provide the analysis of its relevant fixed points in section 2. Next in section 3 we derive the action describing the stochastic master equation evolution of the system through the Doi–Peliti coherent-state path integral formalism. We proceed by reducing the system to a coupled set of stochastic non-linear partial differential equations relevant to the three-species coexistence fixed point in section 4. In section 5, we make crucial use of time scale separation near the Hopf bifurcation to obtain the invariant two-dimensional reactive manifold. On dimensional reduction, we recast the ensuing dynamical problem using normal forms, and derive the CGLE in section 6 for small fluctuations within the species coexistence phase. We explicitly articulate the noise contributions to the stochastic CGLE in section 7. In the concluding section 8, we summarize our assumptions pertinent to this derivation, and comment on the range of applicability of the CGLE mapping.

2 Stochastic May-Leonard model

2.1 Model description

The May–Leonard model for cyclic competition comprises the following independent stochastic processes

[TABLE]

Here, the subscript $\alpha=1,2,3$ denotes the three competing populations, and the identification $X_{4}=X_{1}$ is implicit. Note that for simplicity we study the symmetric situation, for which identical reaction rates are implemented for all three species.

In addition to these cyclic predation and reproduction reactions, we prescribe a population-limiting intra-species competition reaction,

[TABLE]

We justify the addition of this reaction (which is not explicitly listed, e.g., in Ref. [15]) as follows: The May–Leonard model is usually simulated with population restrictions on the lattice (often, at most a single particle is allowed to occupy any site), representing finite local carrying capacities $\rho$ for each species. For our subsequent theoretical analysis, it turns out that the pair coagulation (3) represents a simpler ‘soft’ implementation of local population density suppression than enforcing ‘hard’ site restrictions; upon coarse-graining, both description become essentially equivalent (with an effective rate $\lambda^{\prime}$ that can be expressed in terms of the reproduction rate $\mu$ and the local carrying capacity $\rho$) [30]. Indeed, allowing multiple particles of either species to occupy each lattice site will enable us to utilize bosonic field operators in section 3 below.

In addition to the above on-site reactions, we allow for populations to migrate across the lattice through nearest-neighbor hopping (diffusion in the continuum limit),

[TABLE]

where the subscript $i$ denotes a lattice site (vector) index. The subsequent analysis of this stochastic spatially extended May–Leonard model variant will be carried out in the thermodynamic limit on an infinite lattice.

2.2 Mean-field analysis

We remark that our model represents a specialized case of generalized Lotka–Volterra systems, for which analyses of Hopf bifurcations and global Lyapunov functions are well-established, see, e.g., Ref. [34]. In the (much simplified) case of a well-mixed system or for very fast diffusivity $D^{\prime}\gg\mu,\sigma^{\prime},\lambda^{\prime}$, spatial correlations are washed out and mean-field mass action factorization is applicable. The corresponding coupled rate equations for the three spatially uniform particle densities $a_{\alpha}(t)$ read

[TABLE]

where we define the continuum reaction rates $\lambda=c^{d}\lambda^{\prime}$ and $\sigma=c^{d}\sigma^{\prime}$, with $c$ denoting the lattice spacing. Mean-field steady states are stationary solutions of the rate equations (5), $da_{\alpha}(t)/dt=0$. There exist two sets of extinction fixed points, namely (i) $a_{1}=\mu/\lambda$, $a_{2}=0=a_{3}$ (and cyclic permutations thereof); and, provided $\lambda>\sigma$, (ii) $a_{1}=\mu(\lambda-\sigma)/\lambda^{2}$, $a_{2}=0$, $a_{3}=\mu/\lambda$ (and cyclic permutations thereof). Both sets of extinction fixed points describe absorbing states in the stochastic system. In addition, we have a symmetric three-species coexistence fixed point $a_{1}=a_{2}=a_{3}=\bar{a}=\mu/(\sigma+\lambda)$.

Linearizing about this coexistence fixed point, and collecting the species densities in a three-component vector $\boldsymbol{a}=(a_{1},a_{2},a_{3})^{T}$, we have $\delta\dot{\boldsymbol{a}}=L\,\delta\boldsymbol{a}$, with the linear stability matrix

[TABLE]

Its eigenvalues at the coexistence fixed point are

[TABLE]

the associated eigenvectors are depicted in Fig. 1. The first eigenvalue $\nu_{0}=-\mu$ is always negative, implying stability against small perturbations along its corresponding eigenvector, which will relax exponentially in time $\sim e^{-\mu t}$ with decay rate $\mu$. The other two complex conjugate eigenvalues describe either exponentially damped or growing temporal oscillations with linear frequency $\omega_{0}=\mu\sigma\sqrt{3}/2(\sigma+\lambda)$. For $\lambda>\sigma/2$, the real part of $\nu$ is negative, and the ensuing limit cycle stable, contracting exponentially in time with decay rate $\mu(2\lambda-\sigma)/2(\sigma+\lambda)$. Conversely, for $\sigma>2\lambda$ we obtain an unstable limit cycle, with an exponentially growing amplitude. We note that in the associated spatially extended system, these unstable limit cycles generate spiral structures, whose amplitudes are ultimately constrained by the non-linear terms.

In the language of bifurcation theory, when a pair of complex conjugate stability eigenvalues cross the imaginary axis of the complex plane, the associated dynamical system displays a Hopf bifurcation; in our May–Leonard model variant it is located at $\epsilon=(\sigma-2\lambda)/2(\sigma+\lambda)\to 0$. In the vicinity of the Hopf bifurcation, i.e., as the dimensionless parameter $|\epsilon|\ll 1$, the temporal evolution in the ‘reactive’ plane spanned by the two eigenvectors corresponding to the eigenvalues $\nu$ and $\nu^{*}$ is very slow compared to the fast relaxing mode along the orthogonal eigenvector corresponding to $\nu_{0}$. The presence of the Hopf bifurcation thus provides us with a natural time scale separation for the dynamical eigenmodes of the May–Leonard system (see also Ref. [35]). The provision of a small expansion parameter $\epsilon$ distinguishes our model variant from the simpler one studied in Ref. [15]. The deterministic derivation of the CGLE carried out in Ref. [28] too utilizes time scale separation afforded through the small value of $\epsilon$ near the Hopf bifurcation.

3 Doi–Peliti coherent-state path integral

In order to systematically account for the intrinsic fluctuations in the system, we begin with the stochastic master equation. By applying the Doi–Peliti formalism [36, 37] (for more detailed pedagogical expositions, see Refs. [38, 39]), we derive an effective field theory action that captures the evolution of the system in the continuum limit, while faithfully incorporating its non-linearities and stochasticity due to the on-site reactions. We then construct an equivalent system of coupled Langevin equations which will serve as starting point to the subsequent derivation of the CGLE.

3.1 Doi–Peliti operator representation of the master equation

A specific configuration in this context entails enumerating the integer occupation numbers $n_{\alpha,i}\geq 0$ for each species $\alpha$ at every lattice site $i$. The state of the system at time $t$ is then given as a sum over all possible such configurations, weighted with their probabilities $P(n_{\alpha,i};t)$ which change over time through transitions with rates associated with the possible reactions allowed in the model. The continuous-time stochastic master equation describes the dynamical evolution of the system through balancing gain and loss terms for the configurational probabilities. For the on-site stochastic reactions (1)–(3), excluding for now nearest-neighbor hopping, the master equation reads explicitly:

[TABLE]

For the initial configuration, we assume the particle numbers on each site $i$ to be drawn from independent Poisson distributions with mean initial population densities $\bar{n}_{\alpha}=N_{\alpha}/N$, with the total number of lattice sites $N$ and $N_{\alpha}=\sum_{i}n_{\alpha,i}$, i.e.:

[TABLE]

We then adopt the ladder operator approach first associated with quantum harmonic oscillators to build up a many-particle Fock space on each site with the basic bosonic commutation relations $[a_{\alpha,i},a_{\beta,j}]=0=[a^{\dagger}_{\alpha,i},a^{\dagger}_{\beta,j}]$, $[a_{\alpha,i},a^{\dagger}_{\beta,j}]=\delta_{ij}\delta_{\alpha\beta}$ and particle number eigenstates $|n_{\alpha i}\rangle$ satisfying $a_{\alpha,i}\,|n_{\alpha,i}\rangle=n_{\alpha,i}\,|n_{\alpha,i}-1\rangle$ and $a^{\dagger}_{\alpha,i}\,|n_{\alpha,i}\rangle=|n_{\alpha,i}+1\rangle$. Any arbitrary state can then be written as a series of creation operators acting on an empty vacuum state, $|n_{\alpha,i}\rangle=\prod_{i}\prod_{\alpha}(a^{\dagger}_{\alpha,i})^{n_{\alpha,i}}\,|0\rangle$. A general state vector of the system $|\Phi(t)\rangle$ is then conveniently defined as the linear superposition of the various Fock space configurations of the system weighted by their associated probabilities,

[TABLE]

The stochastic master equation can then be rewritten in this operator formalism as a time evolution operator $H$ acting on the state vector,

[TABLE]

where the ‘pseudo-Hamiltonian’ is a sum of products of (normal-ordered) creation and annihilation operators.

For our stochastic May–Leonard problem, the evolution operator $H_{\rm reac}$ corresponding to the local on-site reactions becomes

[TABLE]

We may also construct the non-local pseudo-Hamiltonian $H_{\rm diff}$ describing hopping transport on the lattice or unbiased diffusion with continuum diffusivity $D=c^{2}D^{\prime}$:

[TABLE]

where we sum over nearest-neighbor pairs $\langle i,j\rangle$. We then express the total pseudo-Hamiltonian of the system as a sum of both these contributions, $H=H_{\rm reac}+H_{\rm diff}$.

3.2 Field theory action in the coherent-state basis

Following the procedures detailed in Refs. [38, 39], one can compute the expectation values of any observable $\mathcal{O}(\{n_{\alpha,i}\})$ as a path integral over a coherent-state basis,

[TABLE]

the $\psi_{\alpha,i}^{*}$ and $\psi_{\alpha,i}$ respectively denoting the complex-valued left and right eigenvalues of $a^{\dagger}_{\alpha,i}$ and $a_{\alpha,i}$. The associated Doi–Peliti action for local on-site reactions is

[TABLE]

In this expression, the last term originates from the initial Poissonian product distribution, while the penultimate term corresponds to the field computed at the final time. The terms in square brackets, referred to as the ‘bulk’ part of the action, are relevant to our subsequent analysis and derivation. The $H(\psi_{\alpha,i}^{*}(t),\psi_{\alpha,i}(t))$ term in the bulk action is simply the evolution operator $H$, obtained by replacing the creation and annihilation operators with the corresponding coherent-state eigenvalue fields.

Finally, we proceed to the continuum limit of our problem (lattice constant $c\to 0$) by substituting $\sum_{i=1}^{N}\to c^{-d}\int d^{d}x$, $\psi_{\alpha,i}(t)\to c^{d}a_{\alpha}(\vec{x},t)$, and $\psi_{\alpha,i}^{*}(t)\to 1+\tilde{a}_{\alpha}(\vec{x},t)$, and thus obtain the coarse-grained bulk action for the May–Leonard model,

[TABLE]

Here, the continuum pseudo-Hamiltonian for the on-site reactions reads explicitly

[TABLE]

with all contributions written in terms of the continuum reaction rates $\sigma=c^{d}\sigma^{\prime}$, $\lambda=c^{d}\lambda^{\prime}$, and $D=c^{2}D^{\prime}$.

4 Langevin description

In order to derive the CGLE, we seek a set of coupled stochastic partial differential equations

[TABLE]

with $\langle\boldsymbol{\zeta}\rangle=0$ and associated noise correlations

[TABLE]

A set of coupled stochastic partial differential equations of this form can be cast in terms of an equivalent dynamical Janssen de-Dominicis response functional [40, 41, 39]

[TABLE]

Hence, upon identifying the response functional (20) with the bulk Doi–Peliti action in eqs. (16) and (17), one arrives at coupled Langevin equations

[TABLE]

for the three complex fields $a_{\alpha}(\vec{x},t)$.111The identification of the ‘shifted’ Doi–Peliti action with a dynamical response functional is associated with certain mathematical subtleties; for an up-to-date exposition and analysis, see Ref. [42]. On comparison with the mean-field rate equations (5), we note the presence of additional diffusion and multiplicative noise contributions that are governed by the (symmetric) stochastic correlation matrix

[TABLE]

We proceed with a linear variable transformation to fluctuating dynamical fields relative to the three-species (mean-field) coexistence fixed-point densities $\bar{a}=\mu/(\sigma+\lambda)$,

[TABLE]

The $b$ fields thus represent deviations about the mean-field stationary concentrations, governed by stochastic partial differential equations of the form

[TABLE]

with ${\bf f}[b_{\alpha}]={\bf F}[a_{\alpha}\to\bar{a}+b_{\alpha}]$ and $\langle\zeta^{\prime}_{\alpha}\,\zeta^{\prime}_{\beta}\rangle=2B_{\alpha\beta}\,\delta(\vec{x}-\vec{x\,}^{\prime})\delta(t-t^{\prime})$, $B_{\alpha\beta}=L_{\alpha\beta}[a_{\alpha}\to\bar{a}+b_{\alpha}]$.

5 Invariant manifold

Our main goal is to explore if the stochastic May–Leonard problem can possibly be further simplified by reducing the dynamical degrees of freedom from three to two, at least for sufficiently small fluctuations in the species coexistence regime. To this end, we essentially follow the procedure for the deterministic model outlined in Ref. [15], and utilize its invariant manifold, i.e., the dynamical subspace of the system invariant to perturbations. Within the linear approximation, the deterministic dynamics will stabilize on the reactive plane spanned by the eigenvectors associated with the complex eigenvalues $\nu$ and $\nu^{*}$ normal to the eigendirection with negative eigenvalue $\nu_{0}=-\mu$. Near the Hopf bifurcation, this stable mode relaxes very fast in comparison with the oscillations in the invariant manifold and can hence be eliminated. We shall first apply a suitable linear variable transformation that allows us to orient the dynamical degrees of freedom along the stable (fast) and reactive (slow) directions. Note that the dynamics of the (to linear order) fast relaxing mode will be affected by non-linear couplings as well as noise cross-correlations to the two oscillating modes, and in turn feed back into the slow dynamics on the invariant manifold. Consequently, as a second step, we shall exploit time scale separation afforded by the vicinity to the Hopf bifurcation and slave the fast mode to the two slow degrees of freedom and thereby account for these non-linear effects. We remark that while we shall carry out the analysis in the Langevin representation here, we could have equivalently performed all required transformations within the associated dynamical functionals.

5.1 Dynamical variable transformation

We proceed with a non-orthogonal dynamical variable transformation but otherwise akin to a rotation aligning the dynamical variables within and normal to the invariant manifold. Introducing the column vector ${\bf c}=(c_{1},c_{2},c_{3})^{T}$ for the new fields, we apply a transformation $\bf{c}=\mathcal{R}\bf{b}$, where

[TABLE]

Following this transformation, we obtain a set of stochastic partial differential equations in the new $c$ fields [43],

[TABLE]

where $\langle\eta_{\alpha}\,\eta_{\beta}\rangle=2\tilde{B}_{\alpha\beta}\,\delta(\vec{x}-\vec{x\,}^{\prime})\delta(t-t^{\prime})$, $\tilde{B}_{\alpha\beta}=[\mathcal{R}B\mathcal{R}^{T}]_{\alpha\beta}$ ($\alpha,\beta=1,2,3$). Explicitly, the resulting coupled Langevin equations become

[TABLE]

Inspection of eq. (32) confirms that the linear couplings of $c_{3}$ to the other modes have indeed been eliminated. Note that the ‘mass’ terms, i.e., constant coefficients of the $c_{1}$ and $c_{2}$ terms in eqs. (30) and (31), respectively, are $\mu\epsilon$. When compared to the corresponding relaxation rate $\mu$ in the term linear in $c_{3}$ in eq. (32), we see again that as $\epsilon\to 0$, $c_{3}$ relaxes much faster than $c_{1}$ and $c_{2}$ (to linear order). This time scale separation provides the rationale for the subsequent elimination of the $c_{3}$ degree of freedom. Yet if the dimensionless parameter $\epsilon$ is not small, no such reduction to just two dynamical degrees of freedom can be justified, and hence a mapping to the CGLE is inapplicable. One then needs to retain all three dynamical modes and their non-linear couplings and noise (cross-)correlations to faithfully describe the dynamics of the system. The effective description of the May–Leonard model in terms of the CGLE consequently holds only near the Hopf bifurcation.

Let us further consider the linear terms in the equations of motion for the ‘slow’ modes $c_{1}$ and $c_{2}$. We see that for $\epsilon<0$, i.e., $\sigma<2\lambda$, deviations from the coexistence fixed point densities will exponentially relax to zero. Comparison with the diffusive spreading term yields a characteristic length scale $\xi_{c}=\sqrt{\frac{2D\,(\lambda+\sigma)}{\mu\,(2\lambda-\sigma)}}$ which describes the typical extent of spatial patches for each species. One would not expect to encounter other more interesting spatio-temporal structures in this regime. On the other hand, for $\epsilon>0$ or $\sigma>2\lambda$, the coexistence fixed point becomes unstable, and correspondingly, spatially homogeneous species distributions develop instabilities at wavelengths larger than $\lambda_{c}=\sqrt{\frac{2D\,(\sigma+\lambda)}{\mu\,(\sigma-2\lambda)}}$. Along with the associated periodic temporal oscillations, this generates spiral structures of typical size $\lambda_{c}$. Note that these are stabilized due to saturating non-linearities in the fluctuating fields $c_{1}$ and $c_{2}$, which we capture in the next section.

5.2 Dimensional reduction

We seek to describe the invariant manifold through a function $c_{3}=\mathcal{G}(c_{1},c_{2})$ that expresses the fast variable in our problem through the two slow modes. It is a difficult problem to obtain $\mathcal{G}$ exactly to all orders. However, one can try the simplest non-trivial ansatz compatible with the rotational symmetry in the reactive plane spanned by $c_{1}$ and $c_{2}$; to first approximation, we set

[TABLE]

where $K$ is a constant, to be determined next. Differentiation with respect to time gives $\partial_{t}c_{3}=2K\,(c_{1}\partial_{t}c_{1}+c_{2}\partial_{t}c_{2})$. We then just substitute the deterministic part of eqs. (30)-(32) to identify $K$ to second order in $c_{1}$ and $c_{2}$. To simplify the calculation involving the Laplacian operators, we operate in spatial Fourier space, replacing the $\nabla^{2}$ operators with $-p^{2}$. This yields

[TABLE]

Near the Hopf bifurcation, where $\sigma\approx 2\lambda$ and in the long-wavelength limit $p\ll\sqrt{\mu/D}$, the inverse bracket in (34) approximately becomes $1+Dp^{2}/\mu$, whence we may set

[TABLE]

and subsequently restore all $p^{2}$ terms with the differential operators $-\nabla^{2}$. Inserting this result (35) into eqs. (30), (31) we obtain

[TABLE]

The wavevector-dependent contributions from (34) would induce additional subleading non-linearities containing spatial derivatives. In the noise covariance matrix we also simply apply the substitution (35), as detailed in section 7 below.

6 Derivation of the CGLE from the normal form

The normal form of a dynamical system encapsulates its essential behavior. Normal forms facilitate the description of non-linear dynamics near bifurcations in a natural way, and thus enable classification schemes [44, 45, 46, 47]. Reichenbach, Mobilia, and Frey showed that the normal form of the deterministic May–Leonard model allows its characterization in terms of the CGLE as its effective dynamical description [15]. We therefore proceed to obtain a non-linear variable transformation $c_{i}\to z_{i}$ with the goal to eliminate the quadratic terms in eqs. (36) and (37):

[TABLE]

This non-linear transformation enables us to write the effective dynamical system in the following form, up to quartic terms in the complex fields ${\bf z}=(z_{1},z_{2})$:

[TABLE]

with coefficients that depend on the original reaction rates according to

[TABLE]

These coefficients encode information about the spatio-temporal pattern formation present in this system. Its linear instability is apparent for $k_{1}>0$. The resulting oscillatory or spiral instability is saturated by the coefficient of the nonlinear term $k_{3}>0$. The associated stochastic noise terms convey information about the intrinsic fluctuations in the system and are described in the subsequent section 7.

In general, the fields $z_{1}$ and $z_{2}$ are complex-valued, and the two independent Langevin equations (40) hence contain twice as many degrees of freedom as a single dynamical equation for one complex field. For small fluctuations near the species coexistence fixed point, one may assume the deviations from $\bar{a}$ to be constrained to the real axis; eqs. 40 then precisely match the partial differential equations for the real and imaginary parts of the CGLE complex order parameter field, respectively. The connection to the noisy complex Ginzburg–Landau equation (CGLE) is thus borne out upon constructing the two linear combinations $\phi=z_{1}+iz_{2}$ and $\chi=z_{1}-iz_{2}$, or $z_{1}=(\phi+\chi)/2$, $z_{2}=(\phi-\chi)/2i$, with $z_{1}^{2}+z_{2}^{2}=\phi\chi$. These obey the Langevin equations

[TABLE]

For small and real fluctuations $z_{1}$, $z_{2}$, obviously $\chi=\phi^{*}$, and eq. (42) turns into the desired CGLE,

[TABLE]

while (43) is merely its complex conjugate.

7 Noise covariance matrix calculation

An analysis of the fully stochastic system enables us to systematically account for internal reaction noise in the system. The noise correlation matrix obtained from the stochastic partial differential equations through this path integral approach is modified during the course of our derivation of the CGLE. Here we describe the steps that lead to the additive noise contributions in the final CGLE (44). 222Tracking the full noise correlation matrix in the course of all intermediate steps is rather cumbersome. We have employed Mathematica to aid us with detailed book-keeping and algebraic simplifications.

• •

The noise correlation matrix $L$ written for the fields $a$ is modified into the matrix $B$ expressed in terms of the $b$ field variables which are just deviations from the mean density $\bar{a}$,

[TABLE]

• •

The rotation-like dynamical variable transformation outlined in section 5.1 modifies the correlation matrix as follows,

[TABLE]

• •

In section 5.2, we simply use the substitution (35) in $\tilde{B}$. We note that the resulting noise contributions for the fast field are purely multiplicative and of the order $c_{1}^{2}$, $c_{2}^{2}$; there are no constant terms, independent of the fluctuating fields. Dimensional reduction of the matrix is withheld until the last step, see below.

• •

Our final transformation is the non-linear one outlined in section 6, where we use eqs. (38), (39) in the matrix obtained in the previous step. To zeroth order in the fluctuating fields its entries are constants:

[TABLE]

As a consequence of rotational symmetry, this matrix is diagonal. Note that to this order the dynamical variable transformation (46) does not ‘rotate’ the noise correlation matrix, and hence generates no constant additive term for the fast field, rendering it kinetics deterministic.

• •

We note that in general, one would need to eliminate the fast degree of freedom in the correlation matrix to obtain the conditioned noise correlator [48], $\tilde{B^{\prime}}=\tilde{B}_{22}-\tilde{B}_{23}\tilde{B}_{33}\tilde{B}_{32}$; here, \tilde{B}_{22}=\left(\begin{array}[]{cc}3\mu^{2}\sigma/2(\sigma+\lambda)^{2}&0\\ 0&3\mu^{2}\sigma/2(\sigma+\lambda)^{2}\end{array}\right), $\tilde{B}_{23}^{T}=\tilde{B}_{32}=(0\;\ 0)$, $\tilde{B}_{33}=0$.

• •

To lowest order in the field fluctuations, the final noise correlator $\tilde{B}^{\prime}$ is simply diagonal and constant

[TABLE]

and hence describe mere additive noise in the stochastic CGLE (44).

We emphasize that the assumption of small amplitude fluctuations near the Hopf bifurcation enables us to justify keeping only the zeroth-order constant terms in the noise correlators. We note that for larger fluctuations, additional multiplicative noise terms would come into play and indeed become dominant near the absorbing extinction fixed points. As stochastic trajectories reach the heteroclinic cycles in the system, the effective description in terms of the CGLE is thus rendered invalid. As seen in (47), for small field fluctuations the noise in the fast relaxational eigendirection is decoupled from the stochastic dynamics on the slow reactive manifold; there are no cross-correlations between the fast relaxing mode and the slow damped oscillatory modes.

8 Results and conclusions

In the spatially extended stochastic May–Leonard model variant under consideration here, a Hopf bifurcation separates two regimes: (i) For $\sigma<2\lambda$, there exists a stable limit cycle, and all dynamical degrees of freedom relax towards the stationary three-species coexistence fixed point. The lattice system is correspondingly characterized by finite species patches with typical size $\xi_{c}$. (ii) For $\sigma>2\lambda$, in contrast the deterministic limit cycles become unstable, inducing an instability of a spatially homogeneous state against the spontaneous formation of spiral structures. Analyzing the associated coefficients of the CGLE (44), we note that this parameter region corresponds to the real part $k_{1}>0$ of the coefficient of the linear term. The imaginary part $k_{2}$ is always positive. The growth of the spiral spatio-temporal patterns is ultimately inhibited by the non-linear terms in the deterministic dynamics. For stable spirals, the real part $k_{3}$ of the coefficient of the non-linear term must be greater than zero, which again implies $\sigma>2\lambda$, consistent with the instability condition for the limit cycle. We note that for small fluctuations in the dynamical fields, the noise components of the CGLE in the three-species coexistence region are merely additive to the lowest order. Hence, at least in this regime we may draw a correspondence between the original microscopic reaction rates of the model and the coefficients of the effective continuum evolution equation, thus quantitatively describing spiral pattern formation in the coexistence phase.

In forthcoming work, we plan to verify our results through detailed Monte Carlo simulations for stochastic May–Leonard models on a two-dimensional lattice. We hope to map the continuum theoretical parameters to corresponding rates on a discrete lattice. We furthermore intend to corroborate the change in pattern formation from patches to spirals as we move across the Hopf bifurcation in parameter space, i.e., from $\sigma<2\lambda$ to $\sigma>2\lambda$.

To summarize, we demonstrate that a stochastic treatment of the fluctuations due to internal reaction noise is possible for the May–Leonard model through the field-theoretic formalism. We derive a fully stochastic set of partial differential equations (30)-(32) that incorporates the intrinsic stochasticity of the system in the continuum limit. Specifically near the Hopf bifurcation, i.e., for $0<\epsilon\ll 1$, we may exploit the emerging time scale separation and eliminate one fast relaxing mode, which at least for small amplitude fluctuations leads us to a stochastic complex Ginzburg–Landau equation as a coarse-grained dynamical description. We derive the most relevant non-trivial noise effect terms for this effective CGLE and see that to the lowest order in the fluctuating fields, additive noise dominates, and cross-correlations between from the fast and slow degrees of freedom are absent.

We emphasize that for generic parameter values away from the Hopf bifurcation, one cannot achieve a similar dimensional reduction as there is no adequate separation of time scales. Indeed, the relaxing degree of freedom will couple to the two oscillating modes through non-linear feedback. A mapping of the stochastic May–Leonard model to the CGLE is thus not in general possible. Also, the validity of the non-linear transformation employed in our derivation is constrained to small fluctuations of the dynamical variables. Higher-order terms neglected in this procedure become non-negligible when deviations from the stationary coexistence fixed point become appreciable. For example, phase space trajectories could then traverse the heteroclinic orbit not captured by the simple CGLE. It is important to note that the field-theoretic Doi–Peliti formalism and the equivalent Langevin description in terms of three dynamical fields remain applicable for arbitrary values of the intrinsic rates, and could be utilized for further detailed mathematical analysis of the spatially extended stochastic May–Leonard model.

We would like to acknowledge and thank Darka Labavić, Hildegard Meyer-Ortmanns, and Mauro Mobilia for fruitful discussions and helpful suggestions associated with this work.