A trio of heteroclinic bifurcations arising from a model of spatially-extended Rock-Paper-Scissors

TL;DR

This paper investigates how heteroclinic bifurcations in a spatial Rock-Paper-Scissors model lead to complex wave and spiral patterns, revealing three novel bifurcation types in PDEs related to population dynamics.

Contribution

It identifies and analyzes three new heteroclinic bifurcations in PDE models derived from a classical heteroclinic cycle, linking them to wave and spiral pattern formation.

Findings

Heteroclinic bifurcations create and destroy travelling waves in PDEs.

Three distinct heteroclinic bifurcations are identified and characterized.

Results enhance understanding of spiral patterns in population models.

Abstract

One of the simplest examples of a robust heteroclinic cycle involves three saddle equilibria: each one is unstable to the next in turn, and connections from one to the next occur within invariant subspaces. Such a situation can be described by a third-order ordinary differential equation (ODE), and typical trajectories approach each equilibrium point in turn, spending progressively longer to cycle around the three points but never stopping. This cycle has been invoked as a model of cyclic competition between populations adopting three strategies, characterised as Rock, Paper and Scissors. When spatial distribution and mobility of the populations is taken into account, waves of Rock can invade regions of Scissors, only to be invaded by Paper in turn. The dynamics is described by a set of partial differential equations (PDEs) that has travelling wave (in one dimension) and spiral (in two dimensions) solutions. In this paper, we explore how the robust heteroclinic cycle in the ODE manifests itself in the PDEs. Taking the wavespeed as a parameter, and moving into a travelling frame, the PDEs reduce to a sixth-order set of ODEs, in which travelling waves are created in a Hopf bifurcation and are destroyed in three different heteroclinic bifurcations, depending on parameters, as the travelling wave approaches the heteroclinic cycle. We explore the three different heteroclinic bifurcations, none of which have been observed in the context of robust heteroclinic cycles previously. These results are an important step towards a full understanding of the spiral patterns found in two dimensions, with possible application to travelling waves and spirals in other population dynamics models.

Tables

Table 1. Table 1: Classification of eigenvalues. P ⊖ L symmetric-difference 𝑃 𝐿 P\ominus L denotes the orthogonal complement in P 𝑃 P of the subspace L 𝐿 L .

Eigenvalue class	Subspace
Radial ($r$)	$L_j\equiv P_{j-1}\cap P_j$
Contracting ($c$)	$V_j\left(c\right)=P_{j-1}\ominus L_j$
Expanding ($e$)	$V_j\left(e\right)=P_j\ominus L_j$
Transverse ($s$)	$V_j\left(s\right)={\left(P_{j-1}+P_j\right)}^{⟂}$

Table 2. Table 2: Eigenvalues of the equilibrium ξ A subscript 𝜉 𝐴 \xi_{A} in ( 4 ).

Label	Eigenvalues
Radial	${\lambda }_r^{\pm }=\frac{1}{2}\left(\gamma \pm \sqrt{{\gamma }^2+4}\right)$
Contracting	${\lambda }_c^{\pm }=\frac{1}{2}\left(\gamma \pm \sqrt{{\gamma }^2+4\left(\sigma +\zeta \right)}\right)$
Expanding (${\gamma }^2-4\zeta >0$)	${\lambda }_e^{\pm }=\frac{1}{2}\left(\gamma \pm \sqrt{{\gamma }^2-4\zeta }\right)$
Expanding ()	${\lambda }_e^{\pm }={\lambda }_e^R\pm i{\lambda }_e^I=\frac{1}{2}\left(\gamma \pm i\sqrt{4\zeta -{\gamma }^2}\right)$

Table 3. Table 3: Definitions of the regions of parameter space shown in Fig.  3 and eigenvalue properties therein.

Region	Definition	Eigenvalue properties
1		${\lambda }_e^{\pm }\in ℝ$,
2	$\zeta >\frac{\sigma }{2}$,	${\lambda }_e^{\pm }\in ℝ$,
3		${\lambda }_e^{\pm }\in ℂ$,
4	$\zeta >\frac{{\gamma }^2}{4}$,	${\lambda }_e^{\pm }\in ℂ$,
5	,	${\lambda }_e^{\pm }\in ℝ$,

Equations

$$\dot{a}$$

$$\dot{b}$$

$$\dot{c}$$

$$\dot{a}$$

$$\dot{b}$$

$$\dot{c}$$

$$\dot{x}=f(x),x\in {\mathbb{R}}^N.$$

$$\frac{\partial a}{\partial t}+\gamma \frac{\partial a}{\partial z}$$

$$\frac{\partial b}{\partial t}+\gamma \frac{\partial b}{\partial z}$$

$$\frac{\partial c}{\partial t}+\gamma \frac{\partial c}{\partial z}$$

$$a_z=u,$$

$$u_z=\gamma u-a(1-(a+b+c)-(\sigma +\zeta )b+\zeta c),$$

$$b_z=v,$$

$$v_z=\gamma v-b(1-(a+b+c)-(\sigma +\zeta )c+\zeta a),$$

$$c_z=w,$$

$$w_z=\gamma w-c(1-(a+b+c)-(\sigma +\zeta )a+\zeta b).$$

$$g(a,u,b,v,c,w)=(b,v,c,w,a,u).$$

$$J_{A}=\left(\begin{array}[]{llllll}0&1&0&0&0&0\\ 1&\gamma&1+\sigma+\zeta&0&1-\zeta&0\\ 0&0&0&1&0&0\\ 0&0&-\zeta&\gamma&0&0\\ 0&0&0&0&0&1\\ 0&0&0&0&\sigma+\zeta&\gamma\\ \end{array}\right)$$

$$J_{H}=\left(\begin{array}[]{llllll}0&1&0&0&0&0\\ \frac{1}{3+\sigma}&\gamma&\frac{1+\sigma+\zeta}{3+\sigma}&0&\frac{1-\zeta}{3+\sigma}&0\\ 0&0&0&1&0&0\\ \frac{1-\zeta}{3+\sigma}&0&\frac{1}{3+\sigma}&\gamma&\frac{1+\sigma+\zeta}{3+\sigma}&0\\ 0&0&0&0&0&1\\ \frac{1+\sigma+\zeta}{3+\sigma}&0&\frac{1-\zeta}{3+\sigma}&0&\frac{1}{3+\sigma}&\gamma\\ \end{array}\right)$$

$${\gamma }_H(\sigma ,\zeta )\equiv \frac{\sqrt{3}(\sigma +2\zeta )}{\sqrt{2\sigma (\sigma +3)}},\text{and}{\omega }_H^2=\frac{\sigma }{2(\sigma +3)},$$

$${\gamma }_H(\sigma ,\zeta )<\gamma <\sqrt{\frac{2}{\sigma }}\zeta +\sqrt{\frac{\sigma }{2}},$$

$${\gamma }_H(\sigma ,\zeta )<\gamma <2\sqrt{\zeta }.$$

$$x_e^A={\lambda }_e^{-}b-v,y_e^A={\lambda }_e^{+}b-v,x_c^A={\lambda }_c^{-}c-w,y_c^A={\lambda }_c^{+}c-w.$$

$$\frac{d{x}_e^A}{dt}={\lambda }_e^{+}{x}_e^A,\frac{d{y}_e^A}{dt}={\lambda }_e^{-}{y}_e^A,\frac{d{x}_c^A}{dt}={\lambda }_c^{+}{x}_c^A,\frac{d{y}_c^A}{dt}={\lambda }_c^{-}{y}_c^A.$$

$$(r_e^A{)}^2=(x_e^A{)}^2+(y_e^A{)}^2\mathrm{and}\tan {\theta }_e^A=\frac{y_e^A}{x_e^A}.$$

$$x_e^B={\lambda }_e^{-}c-w,y_e^B={\lambda }_e^{+}c-w,x_c^B={\lambda }_c^{-}a-u,y_c^B={\lambda }_c^{+}a-u.$$

$$H_{\mathrm{in}}^A$$

$$H_{\mathrm{out}}^A$$

$$H_{\mathrm{in}}^B$$

$${\Pi }_{\mathrm{loc}}:H_{\mathrm{in}}^A\to H_{\mathrm{out}}^A$$

$$\text{bi}{x}^A(T)={\Pi }_{\mathrm{loc}}(\text{bi}{x}^A(0))$$

$$(x_e^A(T),y_e^A(T),x_c^A(T),y_c^A(T))={\Pi }_{\mathrm{loc}}(x_e^A(0),y_e^A(0),x_c^A(0),h),$$

$${\Pi }_{\mathrm{glo}}:H_{\mathrm{out}}^A\to H_{\mathrm{in}}^B$$

$$\text{bi}{x}^B={\Pi }_{\mathrm{glo}}(\text{bi}{x}^A(T))$$

$$(x_e^B,y_e^B,x_c^B,h)={\Pi }_{\mathrm{glo}}(x_e^A(T),y_e^A(T),x_c^A(T),y_c^A(T))$$

$${\gamma }_{AB}\cap H_{\mathrm{out}}^A={\widehat{\text{bi}x}}^A=({\hat{x}}_e^A,{\hat{y}}_e^A,0,0),{\gamma }_{AB}\cap H_{\mathrm{in}}^B={\widehat{\text{bi}x}}^B=(0,0,0,h)$$

$$({\hat{x}}_e^A{)}^2+({\hat{y}}_e^A{)}^2=h,\frac{{\hat{y}}_e^A}{{\hat{x}}_e^A}=\tan {\hat{\theta }}_e^A.$$

$$x_e^A(t)$$

$$y_e^A(t)$$

Full text

Nonlinearity 32 (2019) 1375–1407. https://doi.org/10.1088/1361-6544/aaf530

1 Introduction

The Rock–Paper–Scissors game, in which Rock blunts Scissors, Scissors cut Paper, and Paper wraps Rock, provides an appealing simple model of cyclic competition between different strategies or species in evolutionary game theory and biology [1, 2]. The game has been invoked as a description of three competing species of E. coli [3] and of three colour-variants of side-blotched lizards [4], but the idea of cyclic competition has arisen also in rotating convection [5] and as the simplest example of a heteroclinic cycle [6].

Imagine a group of people repeatedly playing Rock–Paper–Scissors, with each person favouring one of the three choices, and let $A(t)$, $B(t)$ and $C(t)$ be the number of people playing Rock, Paper or Scissors at any moment of time $t$. Pairs of people are drawn at random and when they play, either it is a tie (if they are drawn from the same group), or one beats the other. In this case, the loser can either adopt the strategy of the winner (dominance–replacement) or the loser can withdraw from the game (dominance–removal). Once removed, players are replaced (up to a maximum number $N$) and are assigned to Rock, Paper or Scissors with probabilities proportional to the number of Rock, Paper or Scissors players. With these dynamics, if all individuals who are playing (for example) Rock are eliminated (through a random fluctuation when the number of Rock players is small), they can never return, which means that Scissors would have no competitors and would eventually wipe out Paper [3]. This process is known as fixation, and since it involves an absorbing state, is guaranteed (in a discrete stochastic model) to happen eventually [7].

In the limit of large $N$, the discrete process becomes continuous and is modelled by three ordinary differential equations (ODEs) [1, 8, 9]:

[TABLE]

where $a(t)=A/N$, $b(t)=B/N$, $c(t)=C/N$, and $\sigma$ and $\zeta$ are non-negative parameters that control the rates of dominance–removal and dominance–replacement respectively, scaled to the rate of replacement. We have assumed symmetry between Rock, Paper and Scissors. $A$, $B$ and $C$ are numbers of individuals, so $a$, $b$ and $c$ are non-negative.

The ODEs (1) have five equilibria with non-negative components: the trivial solution $(a,b,c)=(0,0,0)$, three on-axis equilibria $(1,0,0)$, $(0,1,0)$ and $(0,0,1)$, and a coexistence point with $(a,b,c)=\frac{1}{3+\sigma}(1,1,1)$. For $\sigma>0$ and $\zeta>0$, this system of ODEs has solutions that approach each of three on-axis equilibria in turn, taking progressively longer to cycle around the three points but never stopping [1] (in contrast to eventual fixation in the discrete case).

This gradual slowing down of trajectories as they spend longer and longer near a sequence of equilibria is a characteristic of asymptotically stable heteroclinic cycles. The rate of slowing down is controlled by the ratio of two of the eigenvalues of the on-axis equilibria: these are $\zeta$ and $-(\sigma+\zeta)$, and the amount of time taken for each cycle is a factor of $\frac{\sigma+\zeta}{\zeta}$ longer than the previous one [10]. In this expression it is apparent that allowing either $\zeta=0$ or $\sigma=0$ requires special attention. The situation where the eigenvalue ratio is equal to $1$ ($\sigma=0$, $\zeta>0$) is normally called a resonance bifurcation from the heteroclinic cycle, associated with the creation of a long-period periodic orbit [11, 12]. However, in the ODEs (1), letting $\sigma=0$ is degenerate, in that the coexistence equilibrium has pure imaginary eigenvalues and the ODEs have an invariant plane $a+b+c=1$ on which there is a continuous family of nested periodic orbits parameterised by $abc={}$constant.

In the last decade there has been considerable interest in the dynamics of the Rock–Paper–Scissors game where the players are distributed in space and allowed to move, for example on a two-dimensional square lattice, interacting only with their neighbours. In this case, in the limit of large $N$ and small lattice spacing, the dynamics is described by the partial differential equations (PDEs) [8, 9]:

[TABLE]

where the spatial coordinates $(x,y)$ are scaled so that the diffusion constants (assumed to be equal) are equal to $1$. Typically the PDEs are solved with periodic boundary conditions. The spatial mobility allows for persistent spiral-like or turbulent patterns of Rock, Paper and Scissors [13], in which regions dominated by Rock invade regions of Scissors, which invade regions of Paper, which in turn invade regions of Rock. In the case of spirals, these have a rotating core, with a point where $a=b=c$ at (or close to) the centre, and spiral arms that, far from the core, look like they are one-dimensional periodic travelling wave (TW) solutions of the PDEs (2) [14].

The central question we address in this paper is: what is the connection between travelling waves in the PDEs (2) and heteroclinic cycles? The TWs are periodic orbits in a moving frame of reference, and, taking the wavespeed as a parameter, these periodic orbits originate in a Hopf bifurcation and end when they collide with a heteroclinic cycle [14]. In this paper we find conditions under which TWs with arbitrarily long wavelength can exist as solutions of (2), close to a heteroclinic cycle in the sixth-order ODEs that describe the dynamics in the travelling frame. We find that there are three different ways in which this can happen:

• •

there can be a resonance bifurcation from the heteroclinic cycle in the sixth-order ODEs, at which a positive and a negative eigenvalue have equal magnitude; or

• •

there can be a bifurcation of Belyakov–Devaney type, at which the imaginary part of a pair of complex eigenvalues vanishes; or

• •

there can be a bifurcation of orbit flip type, at which there is a change in the way in which the trajectories between equilibria are oriented.

Although our analysis proceeds along reasonably standard lines, there are several unusual aspects, and the calculations are challenging, not least because the unstable manifolds of the equilibria in the heteroclinic cycles are of high dimension. It turns out that each of these three bifurcations is non-standard and, to our knowledge, has not been observed in the context of heteroclinic cycles before. We are able to find conditions under which each of these three bifurcations occurs, and, to some extent, how the transition from one type to the next occurs. Our results give a much clearer picture of the origin of the one-dimensional TW solutions of the PDEs (2), a first and necessary step in understanding their stability, which in turn is necessary for understanding the stability of the two-dimensional spiral solutions of (2).

The remainder of this paper is structured as follows. We begin in section 2 by reviewing some numerical results from [14] and showing some simulations of the PDEs (2), both in one and two spatial dimensions. We also relate properties of the travelling wave solutions of the PDEs to periodic solutions of a related set of ODEs. In section 3 we review the definitions of heteroclinic cycles and summarise what is already known about ways in which they can bifurcate. We also compare these bifurcations with those seen near homoclinic orbits, and relate these to the new bifurcations we have found. In section 4 we describe the derivation of the ODEs we will be studying for the remainder of the paper. Then in section 5 we derive a Poincaré map which describes the flow close to the heteroclinic cycle in the ODEs. This section contains a lot of calculation but the results are summarised at the start and end of the section for the reader who doesn’t wish to delve into too many of the gritty details. In section 6 we give some further numerical results from simulation of the PDEs for a range of parameter values, and finally in section 7 we look at numerically computed bifurcation diagrams as the parameter $\sigma$ is varied and discuss the limit $\sigma\rightarrow 0$. Section 8 concludes.

2 PDE simulations

We begin with the PDEs for the spatially-extended Rock–Paper–Scissors model as given in equations (2). In figure 1 we show numerical results from the integration of equations (2) in two spatial dimensions, from [14]. A variety of behaviours can be observed, but of particular interest are the spiral-type solutions. When a slice is taken radially through the centre of a spiral, the profile of the solution in the outer part of the spiral resembles a travelling wave in one spatial dimension.

Figure 2 shows the results of numerical integration of equations (2) in one spatial dimension, in a large box of size 500, for $\sigma=3.2$ and $\zeta=1.0,2.0,3.0$. Initial conditions are of small amplitude and randomly generated, and boundary conditions are periodic. The time-space plots show clearly that multiple travelling waves arise from the initial conditions after a short transient. For all three values of $\zeta$, travelling waves of different directions, wavespeeds and wavelengths are evident. In the simulation for $\zeta=1.0$, after about $t=500$, waves consistently travel to the left, and eventually (after being integrated for a longer time period than shown here), this solution has six waves of equal wavelengths (and equal wavespeed) fitting in the periodic box. For the larger values of $\zeta$, the solutions appear more complicated, in particular, faster wavespeeds and smaller wavelengths are evident. We attempt to quantify this further in section 6.

We would like, ultimately, to be able to predict the behaviour of solutions to the PDEs (2); that is, we would like to be able to say whether solutions will eventually asymptote onto a single travelling wave, and what the wavespeed and wavelength of that travelling wave will be. In order to do this, we would need to know both existence and stability criteria, as well as have information about the basins of attraction of the travelling waves. The latter two are difficult problems, and are beyond the scope of this paper, but in order to answer both of those questions, we first need to understand the existence problem, and that is what is addressed in this paper.

Specifically, by relating travelling waves of the PDEs to periodic solutions of a related set of ODEs, we are able to give existence criteria for the travelling waves, and associated with that are minimum and maximum wavespeeds. For a wide range of parameters (those associated with Belyakov–Devaney-type and resonance heteroclinic bifurcations), the wavelength increases monotonically with the wavelength, and the wavelength asymptotes to infinity as the maximum wavespeed is approached. In the parameter regime for which orbit-flip heteroclinic bifurcations are observed, the dispersion relation relating wavelength and wavespeed is non-monotonic, but we are still able to identify a minimum wavespeed and wavelength, and a wavespeed which is approached asymptotically as the wavelength goes to infinity.

3 Review of heteroclinic cycles and bifurcations

Before we begin the calculations, in this section we first include a review of heteroclinic cycles, and the definitions used by Krupa and Melbourne [15] of contracting, expanding, radial and transverse eigenvalues. In this paper, we abuse their nomenclature slightly, and give labels to eigenvalues that don’t quite fit with these definitions, but we find that this is useful nonetheless.

Consider a system of ordinary differential equations

[TABLE]

Then we have:

Definition 1

A heteroclinic cycle is a finite collection of equilibria $\{\xi_{1},\dots,\xi_{n}\}$ of (3), together with a set of heteroclinic connections $\{\phi_{1}(t),\dots,\phi_{n}(t)\}$, where $\phi_{j}(t)$ is a solution of (3) such that $\phi_{j}(t)\to\xi_{j}$ as $t\to-\infty$ and $\phi_{j}(t)\to\xi_{j+1}$ as $t\to\infty$, and where $\xi_{n+1}\equiv\xi_{1}$.

In generic systems, heteroclinic connections between saddles are of high codimension, but if a system contains invariant subspaces they can exist for open sets of parameter values, that is, they are of codimension zero, and are referred to as ‘robust’ [16, 6, 17]. In the work of Krupa and Melbourne [18, 15] and others (e.g., [19, 20, 21, 22, 23, 24, 25, 26]), robust heteroclinic cycles arise due to invariant subspaces which are a by-product of symmetry in the ODEs. In this paper, we show that for the ODEs we are studying, heteroclinic connections exist for open sets of parameter values due to a combination of invariant subspaces and the dimensions of stable and unstable manifolds of equilibria for the flow restricted to these invariant subspaces. An additional difference in our work is that the invariance of the subspaces is not forced by symmetry, but instead by the invariance of extinction in continuous-time population models.

Despite these differences, we continue in the style of Krupa and Melbourne [15]. Let $P_{j}$ be an invariant subspace which contains $\xi_{j}$ and $\xi_{j+1}$. Let $W_{u}|_{P_{j}}(\xi_{j})$ and $W_{s}|_{P_{j}}(\xi_{j+1})$ by the unstable manifold of $\xi_{j}$ and stable manifold of $\xi_{j+1}$ for the flow restricted to $P_{j}$. Then, if $\dim(W_{u}|_{P_{j}}(\xi_{j}))+\dim(W_{s}|_{P_{j}}(\xi_{j+1}))>\dim(P_{j})$, then a heteroclinic connection from $\xi_{j}$ to $\xi_{j+1}$ will be codimension zero, this is, it will persist under small changes to the ODE (so long as the changes preserve the invariant subspaces). If this is true for all $j$, then there exists a robust heteroclinic cycle between the equilibria $\xi_{1},\dots,\xi_{n}$, where robust here means codimension zero.

We further define $L_{j}\equiv P_{j-1}\cap P_{j}$ and clearly $\xi_{j}\in L_{j}$. Following [17], the eigenvalues of the linearisation of $f(x)$ about each equilibrium can be classified according to the subspaces in which the eigenspaces lie, as shown in table 1.

As we will discuss in the following, because we do not require that $P_{j}$ contains the unstable manifold of $\xi_{j}$ (unlike in the definition used by Krupa and Melbourne [15]), we are allowed to have positive radial and/or contracting eigenvalues.

Methods for determining the stability properties of an isolated heteroclinic cycle are in principal well-established [20, 18, 15, 27, 28, 11, 29, 30, 31, 32, 33, 34]: that is, one can construct a Poincaré map, by linearising the flow around the fixed points and the heteroclinic connections. Many examples have been investigated in lower dimensions ($\mathbb{R}^{3}$ and $\mathbb{R}^{4}$ in particular), but in higher dimensions, calculations can become quite intricate. A number of codimension-one bifurcations have been identified in which the stability of robust heteroclinic cycles changes, but issues of stability turn out to be more subtle than might be at first thought (for several examples, see [27, 12, 31, 29]). Two well-studied ways in which heteroclinic cycles can change stability are resonance and transverse bifurcations. A resonance bifurcation [26, 15, 28, 11, 23] occurs when an algebraic condition on the eigenvalues of the equilibria in the cycle is satisfied. Typically, resonance bifurcations are accompanied by the birth or death of a long-period periodic orbit. In a transverse bifurcation from a heteroclinic cycle [20], a local bifurcation causes a transverse eigenvalue of one of the equilibria in the cycle to change sign. This can result in a bifurcating periodic orbit or heteroclinic cycle, depending on the specific situation.

In this paper, we use the standard methods for analysing the dynamics close to a heteroclinic cycle, namely, we construct a Poincaré map which approximates the flow of the differential equations close to the heteroclinic cycle, but as mentioned in the introduction, there turn out to be several subtleties which must be carefully navigated. We do not explicitly compute the stability of the heteroclinic cycle but rather compute conditions for the existence of nearby periodic orbits. We find that long-period periodic orbits can exist close to the heteroclinic cycle, and can appear from the cycle in three different ways: at a resonance bifurcation, at a bifurcation of Belyakov–Devaney type, and at an orbit flip bifurcation. Although resonance bifurcations have been previously studied in the context of robust heteroclinic cycles, the bifurcation we find is of an unusual type, in that the eigenvalues of interest are not those that one would expect [18].

All three of these types of bifurcations have been previously studied in the context of homoclinic orbits, and in many cases are associated with complicated dynamics such as homoclinic-doubling cascades [35, 36]. Useful references for each case include the work of Chow, Deng and Fiedler for resonant homoclinic bifurcations [37], the work of Homburg, Kokubu, Krauskopf and others for orbit flip bifurcations [38, 39], and the work of Belyakov [40] and Devaney [41] for the Belyakov–Devaney bifurcation. However, homoclinic orbits cannot be robust, so each of these phenomena is at least a codimension two bifurcation (there must be another parameter associated with the existence of the homoclinic orbit). In the case of a robust heteroclinic bifurcation, then these phenomena can occur as codimension one, and as such the dynamics associated with the bifurcations may be somewhat different, and indeed, we find that this is the case.

4 Derivation of ODEs and existence of heteroclinic cycles

In this paper, we examine the behaviour of the travelling wave solutions in one dimension, and so we consider equations (2) with only one spatial dimension, so $\nabla^{2}=\frac{\partial^{2}}{\partial x^{2}}$. We move to a travelling frame with wavespeed $\gamma>0$, so define $z=x+\gamma{t}$, then $\frac{\partial}{\partial x}\rightarrow\frac{\partial}{\partial z}$ and $\frac{\partial}{\partial t}\rightarrow\gamma\frac{\partial}{\partial z}+\frac{\partial}{\partial t}$. This results in the following set of PDEs in the travelling frame:

[TABLE]

Travelling wave (TW) solutions in the moving frame have $\frac{\partial}{\partial t}=0$. We thus set $\frac{\partial}{\partial t}=0$, and add additional variables for the first derivative of $a$, $b$ and $c$ with respect to $z$. Therefore, TW solutions of (2) correspond to periodic solutions of the following set of six first-order ODEs:

[TABLE]

Since $a$, $b$ and $c$ are non-negative, we define a positive travelling wave as a periodic solution of (4) with $a,b,c>0$ for all $z$. In an abuse of notation, the independent variable $z$ will be referred to as ‘time’ (and denoted with a ‘$t$’) when we construct Poincaré maps in the next section.

Let $\mathbf{x}=(a,u,b,v,c,w)$, and note that the coexistence and on-axis equilibria of (4) correspond to $\mathbf{x}=\frac{1}{3+\sigma}(1,0,1,0,1,0)$, $\mathbf{x}=(1,0,0,0,0,0)$, $\mathbf{x}=(0,0,1,0,0,0)$ and $\mathbf{x}=(0,0,0,0,1,0)$. We label these equilibria $\xi_{H}$, $\xi_{A}$, $\xi_{B}$ and $\xi_{C}$ respectively. Also note that the ODEs (4) are invariant under the rotation symmetry $g$:

[TABLE]

The Jacobian matrix at $\xi_{A}$ is

[TABLE]

The eigenvalues of $J_{A}$ are given in table 2. Note that we frequently refer to ‘the eigenvalues of $\xi_{A}$’, by which of course we mean the eigenvalues of $J_{A}$. By the symmetry $g$, $\xi_{B}$ and $\xi_{C}$ have the same eigenvalues.

Let the four-dimensional subspace $\{c=w=0\}$ be labelled $P(\xi_{A})$. It can easily be seen that $P(\xi_{A})$ is invariant under the flow of (4). For the dynamics restricted to $P(\xi_{A})$, $\xi_{A}$ has a three-dimensional unstable manifold, and $\xi_{B}$ has a two-dimensional stable manifold. By dimension counting, it is reasonable to expect that these manifolds will intersect, and hence that there is a heteroclinic connection within $P(\xi_{A})$ between $\xi_{A}$ and $\xi_{B}$, which persists under small perturbations. We are able to numerically confirm the existence of a heteroclinic connection for a wide range of parameter values. By symmetry, there is thus a robust heteroclinic cycle between $\xi_{A}$, $\xi_{B}$ and $\xi_{C}$.

As discussed earlier, because our definition of robust heteroclinic cycle did not require the unstable manifold of $\xi_{A}$ to be contained in $P(\xi_{A})$, we can have radial or contracting eigenvalues that have positive real part, and in fact, this is what we find (see table 2). Specifically, we note that $\lambda_{c}^{-}<0<\lambda_{c}^{+}$, $\lambda_{r}^{-}<0<\lambda_{r}^{+}$, and $0<Re(\lambda_{e}^{-})\leq Re(\lambda_{e}^{+})$.

The Jacobian matrix at $\xi_{H}$ is:

[TABLE]

$J_{H}$ has pure imaginary eigenvalues $\pm i\omega_{H}$ when $\gamma=\gamma_{H}(\sigma,\zeta)$, where

[TABLE]

at which point a Hopf bifurcation creates periodic orbits of period (in the $z$ variable) $\Lambda_{H}=\frac{2\pi}{\omega_{H}}$. Numerical analysis of equations (4) with AUTO [14, 42] show that the branch of periodic orbits grows in period as $\gamma$ is increased from the Hopf bifurcation, eventually ending in a heteroclinic bifurcation. The Hopf and heteroclinic bifurcation curves can be seen in Figure 3 for $\sigma=3.2$ as the grey dashed and black solid curves respectively. Also shown in figure 3 are various curves depicting relationships between the eigenvalues (the red, yellow and blue curves), and a curve showing the location of when the heteroclinic connection undergoes an orbit flip (green curve). Recall that the heteroclinic cycle is of codimension zero, and so the orbit flip curve is of codimension one. The orbit flip curve is found by solving a boundary value problem, as explained further in section 5.1.5.

The heteroclinic bifurcation curve in figure 3 is of three different types, depending on the parameters $\zeta$ and $\sigma$. For the value of $\sigma$ used to create figure 3 ($\sigma=3.2$), we have: (a) if $\zeta>\sigma/2=1.6$, the heteroclinic bifurcation is of resonance type, and occurs when $-\lambda_{c}^{-}=\lambda_{e}^{-}$ (where the black curve coincides with the blue curve in figure 3); (b) if $\zeta^{*}<\zeta<\sigma/2=1.6$ (where $\zeta^{*}\approx 0.46$), then the heteroclinic bifurcation is of Belyakov–Devaney type, and occurs when $\lambda_{e}^{I}=0$ (where the black curve coincides with the red curve in figure 3), and (c) if $0<\zeta<\zeta^{*}$, then the heteroclinic bifurcation is of orbit flip type (where the black curve coincides with the green curve in figure 3).

In the first two cases, the Hopf and heteroclinic bifurcation curves denote the existence boundaries of periodic orbits in the ODEs, and hence also of travelling waves in the PDEs. Specifically, the Hopf bifurcation curve indicates the minimum wavespeed $\gamma$ (and minimum wavelength, given by $\Lambda_{H}$ as written after equation (8)), and the heteroclinic bifurcation curve indicates the maximum wavespeed. That is, for $\zeta>\sigma/2$, the allowed wavespeeds are

[TABLE]

where $\gamma_{H}$ is given in (8). For $\zeta^{*}<\zeta<\sigma/2=1.6$, the allowed wavespeeds are

[TABLE]

For $\zeta<\zeta^{*}$, the heteroclinic bifurcation is of orbit-flip type, and there also exists a branch of saddle-node bifurcations of periodic orbits (light grey curve). Here, the right hand boundary for existence of travelling waves is the saddle-node bifurcation curve, not the heteroclinic bifurcation curve. The location of both of these curves depends on global parameters, so here we cannot give an explicit expression for the maximum wavespeed.

In figure 4 we show time-series of periodic solutions of (4) which are close to the three types of heteroclinic bifurcations. The examples are all right-travelling waves (in the PDE setup); left-travelling waves are also possible. In panel (a), we show a periodic orbit close to the heteroclinic resonance bifurcation (near the edge of region 2, the expanding eigenvalues are real). The slopes in the contracting and expanding phases can be seen to be very close to $\lambda_{c}^{-}$ and $\lambda_{e}^{-}$. In panel (b), we show a periodic orbit close to the heteroclinic Belyakov–Devaney-type bifurcation (near the edge of region 4, the expanding eigenvalues are complex). Here, $\lambda_{e}^{I}$ is very close to zero, and the slope in the expanding phase is close to $\lambda_{e}^{R}$. In the contracting phase, we see slopes equal to both the negative contracting eigenvalue, $\lambda_{c}^{-}$, and the positive contracting eigenvalues $\lambda_{c}^{+}$. In panel (c), we show a periodic orbit close to the heteroclinic orbit-flip bifurcation (in region 5, the expanding eigenvalues are real). The slope in the expanding phase is $\lambda_{e}^{+}$, because the periodic solution lies close to the heteroclinic orbit, which is close to tangent to the strong unstable manifold. Again, in the contracting phase we see both the positive and negative slopes. We refer later to periodic orbits which have both a positive and negative slope in the contracting phase as those having a kink – the kink refers to the change in growth rate at the transition from the contracting to the expanding phase.

In long-period orbits such as in Figure 4, the total amount of decay in the contracting phase must balance the growth in the expanding phase; the contracting and expanding phases must be the same length because of the symmetry between the $a$, $b$ and $c$ coordinates in the orbit. Therefore, orbits of this type cannot exist in regions 2 and 3: $\lambda_{c}^{-}$, the only negative non-radial eigenvalue, is less in absolute value than the (real part of the) smaller of the two expanding eigenvalues, and so there can’t be enough decay to balance the growth.

A further point to note is that not all periodic solutions of (4) correspond to positive travelling wave solutions of (2). In particular, because we are considering a population model, we will start with initial conditions (of (2)) which have $a,b,c\geq 0$, and, given reasonable conditions on the smoothness of the initial conditions, it can be shown that $a,b,c\geq 0$ for all $t$ (in (2)). Only periodic solutions of (4) which have $a,b,c>0$ for all $z$ correspond to positive travelling wave solutions of (2). This may be important, because the variables in (4) may change sign along the heteroclinic connections — clearly they will in the case that the expanding eigenvalues are complex.

5 Constuction of a Poincaré map and analysis of heteroclinic bifurcations

In this section we construct a Poincaré map which approximates the dynamics close to the heteroclinic cycle of equations (4) as described in section 4. We treat the cases in which the expanding eigenvalues are real and complex separately, although the computations are quite similar. Regions of real (1, 2 and 5) and complex (3 and 4) eigenvalues are divided by the red curve in figure 3. In this section, we refer to the independent variable of equations (4) as time ($t$) rather than $z$.

The Poincaré map we construct here will follow a trajectory that starts on an incoming section near $\xi_{A}$ and ends on an incoming section near $\xi_{B}$. In both real and complex cases, we define Poincaré sections close to $\xi_{A}$ and $\xi_{B}$, and derive a local map which approximates the flow close to $\xi_{A}$. We combine this with a global map linearised about the location of the heteroclinic connection from $\xi_{A}$ to $\xi_{B}$ and then use the symmetry $g$ to map the coordinates back to a Poincaré section close to $\xi_{A}$. We are able to disregard the radial directions in our computations. This is often done because eigenvalues in the radial directions are negative, with an invocation to an invariant sphere [16]. In our case, we have a positive radial eigenvalue, and the same argument may not hold. However, because of the invariance of the subspace containing the heteroclinic cycle, the radial directions decouple (to lowest order) in the construction of the Poincaré map. Since we are looking for fixed points of the map, rather than computing stability criteria, we can thus find the fixed points, and examine their properties while ignoring the radial direction.

In constructing the Poincaré map, we do not explicitly compute the amount of time $T$ that the trajectory spends close to the equilibria, but leave this as an unknown defined implicitly in terms of the other coordinates in the map: the map is defined in terms of three coordinates, and the time $T$. It then becomes possible to solve the equations for fixed points of the Poincaré map by writing each of the coordinates in term of $T$, allowing us to construct a single equation with a single unknown, $T$. Letting $T$ become large will give us the locations of the heteroclinic bifurcations.

In the case where the expanding eigenvalues are real, we are able to find two different types of solution for large $T$, depending on which terms in this equation dominate: one type of solution generates the resonant bifurcation, and the other generates the orbit-flip bifurcation. In the case where the expanding eigenvalues are complex, we find only one type of solution, corresponding to a Belyakov–Devaney-type bifurcation.

The period of the bifurcating periodic orbit scales differently with the distance from the bifurcation point, depending on the type of bifurcation. Suppose that $\mu$ is a parameter which measures the distance from the heteroclinic bifurcation in each of the three cases, then: (a) in the resonance bifurcation, $\mu\propto|\lambda_{e}^{-}+\lambda_{c}^{-}|$, and $T$ scales like $1/|\mu|$ (see equation (15qtuyabacaeam)); (b), in the Belyakov–Devaney-type bifurcation $\mu\propto|\lambda_{e}^{I}|$, and $T$ scales like $1/|\mu|$ (see equation (15qtuyabacaebt)); and (c) in the orbit flip bifurcation, $\mu\propto A_{3}$, a global constant which determines the angle at which the heteroclinic connection exits a neighbourhood of $\xi_{A}$, and $T$ scales like $\log|\mu|$ (see equation (15qtuyabacaeaq)).

In each case, once we have computed an expression for the fixed points of the Poincaré map, we also check that the corresponding periodic orbits satisfy the condition that $a,b,c>0$ for all time. The coordinates will need to be checked when they are close to $\xi_{A}$: during the transition between equilibria the coordinates will be order 1 and hence will not change sign. In a neighbourhood of $\xi_{A}$, it is clear that $a(t)$ will not change sign, as it is order 1. The heteroclinic connection leaving $\xi_{A}$ lies in an invariant subspace which has $c=0$, so $c$ cannot change sign during the transition from $\xi_{A}$ to $\xi_{B}$. Thus the coordinate which will need to be checked is $b(t)$.

Finally, we will check whether or not we expect the solution to be ‘kinked’.

5.1 Real eigenvalues

To begin, we define new coordinates which we use when the trajectory is near $\xi_{A}$

[TABLE]

Recall that we are interested in solutions which have $b(t)>0$, which in these coordinates, means we must have $y_{e}^{A}>x_{e}^{A}$. The coordinates in (9) are aligned with the eigenvectors of the Jacobian matrix, and so the linearised equations near $\xi_{A}$ can be written

[TABLE]

We will also make use of polar coordinates in the expanding directions, namely $r_{e}^{A}$ and $\theta_{e}^{A}$, defined by

[TABLE]

The constraint $y_{e}^{A}>x_{e}^{A}$ means that $\pi/4<\theta_{e}^{A}<5\pi/4$. We similarly define new coordinates for use near $\xi_{B}$:

[TABLE]

We further write $\bi{x}^{A}=(x_{e}^{A},y_{e}^{A},x_{c}^{A},y_{c}^{A})$ and $\bi{x}^{B}=(x_{e}^{B},y_{e}^{B},x_{c}^{B},y_{c}^{B})$.

We define Poincaré sections, close to $\xi_{A}$ and $\xi_{B}$:

[TABLE]

for some $h\ll 1$.

We will now construct a local map near $\xi_{A}$ and a global map from $\xi_{A}$ to $\xi_{B}$ as follows. Let the time it takes the trajectory to travel from $H_{\mathrm{in}}^{A}$ to $H_{\mathrm{out}}^{A}$ be $T$. The local map is

[TABLE]

where $x_{e}^{A}(T)^{2}+y_{e}^{A}(T)^{2}=h^{2}$, and the global map is

[TABLE]

where again $x_{e}^{A}(T)^{2}+y_{e}^{A}(T)^{2}=h^{2}$. In figure 5 we show a schematic of the expanding dynamics near $\xi_{A}$.

We label the heteroclinic connection between $\xi_{A}$ and $\xi_{B}$ as $\gamma_{AB}$. Recall that the unstable manifold of $\xi_{A}$, $W^{u}(\xi_{A})$, is four-dimensional. The heteroclinic connection is a one-dimensional sub-manifold of $W^{u}(\xi_{A})$. In addition, we also know that the connection lies in the invariant subspace $P(\xi_{A})$ (which has $c=w=0$, equivalently, $x_{c}^{A}=y_{c}^{A}=0$ near $\xi_{A}$ or $x_{e}^{B}=y_{e}^{B}=0$ near $\xi_{B}$). We consider the points at which the heteroclinic connection intersects the Poincaré sections, and write

[TABLE]

where

[TABLE]

The $x_{c}^{B}$ coordinate of $\hat{\bi{x}}^{B}$ is zero because $\gamma_{AB}$ must lie in the stable manifold of $\xi_{B}$, and $x_{c}^{B}$ is the coordinate associated with the positive contracting eigenvalue, $\lambda_{c}^{+}$. The angle $\hat{\theta}_{e}^{A}$ is marked in figure 5. Note that generically, the heteroclinic connection $\gamma_{AB}$ will be tangent to the $y_{e}^{A}$ axis at $\xi_{A}$, and so generically $\hat{\theta}_{e}^{A}$ will be order one. In the orbit flip bifurcation which we consider in section 5.1.5, the heteroclinic connection is tangent to the strong stable manifold, i.e., the $x_{e}^{A}$ axis, and then $\hat{\theta}_{e}^{A}$ will be very close to $\pi$ (so $|\tan\hat{\theta}_{e}^{A}|\ll 1$).

5.1.1 Local map

We consider a trajectory which starts at time $t=0$, at a point $\bi{x}^{A}(0)\in H_{\mathrm{in}}^{A}$, and we write the solution to the equations linearised around $\xi_{A}$ as

[TABLE]

The time it takes the trajectory to travel from $H_{\mathrm{in}}^{A}$ to $H_{\mathrm{out}}^{A}$ is $T$, so $\bi{x}^{A}(T)\in H_{\mathrm{out}}^{A}$, and $T$ is defined by

[TABLE]

This gives the five equations

[TABLE]

which define $x_{e}^{A}(T),y_{e}^{A}(T),x_{c}^{A}(T),y_{c}^{A}(T)$ and (implicitly) $T$ in terms of $x_{e}^{A}(0),y_{e}^{A}(0)$ and $x_{c}^{A}(0)$, thus defining the local map from $H_{\mathrm{in}}^{A}$ to $H_{\mathrm{out}}^{A}$. Note that we do not attempt to solve for $T$ at this stage.

5.1.2 Global map

We next construct the global map from $H_{\mathrm{out}}^{A}$ to $H_{\mathrm{in}}^{B}$. We only consider trajectories which lie close to the heteroclinic connection from $\xi_{A}$ to $\xi_{B}$, so $\theta_{e}^{A}(T)$ will be close to $\hat{\theta}_{e}^{A}$ (see figure 5). Then we write

[TABLE]

and Taylor expand the right hand side around $\hat{\bi{x}}^{A}$ to get

[TABLE]

where we are assuming $(y_{e}^{A}(T)-\hat{y}_{e}^{A})$ and $(x_{e}^{A}(T)-\hat{x}_{e}^{A})$ are small and have used the fact that $({\hat{x}_{e}^{A}})^{2}+({\hat{y}_{e}^{A}})^{2}=h^{2}$.

Recall that a point on $H_{\mathrm{out}}^{A}$ can be defined by the coordinates $x_{c}^{A}(T),y_{c}^{A}(T)$ and $\theta_{e}^{A}(T)$. For a trajectory close to the heteroclinic connection, $x_{c}^{A}(T)$ and $y_{c}^{A}(T)$ are small (since the heteroclinic connection lies in $P(\xi_{A})$ which has $x_{c}^{A}=y_{c}^{A}=0$), and $(\theta_{e}^{A}(T)-\hat{\theta}_{e}^{A})$ is also small. A point on $H_{\mathrm{in}}^{B}$ is defined by the coordinates $x_{c}^{B}$, $x_{e}^{B}$ and $y_{e}^{B}$, which are also all small for a trajectory close to $\gamma_{AB}$ (see equation (13)). Thus, in the global map, to first order, $x_{c}^{B}$, $x_{e}^{B}$ and $y_{e}^{B}$ can be written as a linear combination of $x_{c}^{A}(T)$, $y_{c}^{A}(T)$ and $(\theta_{e}^{A}(T)-\hat{\theta}_{e}^{A})$. In addition, the global map must preserve the invariance of $P(\xi_{A})$. The global map can thus be written to first order as:

[TABLE]

where the $F_{j}$ are order one constants.

Using equation (15qs), we replace $\theta_{e}^{A}(T)$, and renaming the constants gives

[TABLE]

Note that $A_{3}=-F_{1}\hat{y}_{e}^{A}/h^{2}$ and $A_{4}=F_{1}\hat{x}_{e}^{A}/h^{2}$, so $\tan\hat{\theta}_{e}^{A}=-A_{3}/A_{4}$.

Usually in these sorts of calculations, it is assumed that the order one constants which arise in the global map are not functions of the eigenvalues. This is not entirely true, as they will be dependent on the global dynamics, but to leading order, if we are only considering small changes of the eigenvalues (such as near a bifurcation point), then the constants will be close enough to constant that it doesn’t matter. However, in this case, we note that the constants $B_{1}$, $B_{2}$, $C_{1}$ and $C_{2}$ do have a degeneracy near a particular degeneracy of the eigenvalues, arising because of the way we have defined our coordinates.

Specifically, consider the trajectory of $c$ and $w$ during the passage from $H_{\mathrm{out}}^{A}$ to $H_{\mathrm{in}}^{B}$. Both $c$ and $w$ are assumed small, and write $(c_{A},w_{A})$ for the coordinates on $H_{\mathrm{out}}^{A}$ and $(c_{B},w_{B})$ for the coordinates on $H_{\mathrm{in}}^{B}$. Then to lowest order, the global map can be written

[TABLE]

where the $G_{j}$ are indeed generically order one constants. When we rewrite this in terms of $x_{e}^{B}$, $y_{e}^{B}$, $x_{c}^{A}$ and $x_{c}^{B}$, we have

[TABLE]

That is, (referring to (15qtb) and (15qtc))

[TABLE]

There are thus degeneracies in $B_{1}$, $B_{2}$, $C_{1}$ and $C_{2}$ when either $\lambda_{c}^{-}=\lambda_{c}^{+}$ or $\lambda_{e}^{-}=\lambda_{e}^{+}$. The former case doesn’t occur in our system, because we assume that $\sigma,\zeta>0$ (see table 2), but the latter can occur, when $4\zeta=\gamma^{2}$: where the expanding eigenvalues change from being real to complex. In this case, when $\lambda_{e}^{-}=\lambda_{e}^{+}$, then $B_{1}=C_{1}$, and $B_{2}=C_{2}$, and the determinant of the matrix on the left hand side of (15qtux) is $\Delta_{BC}=B_{1}C_{2}-C_{1}B_{2}=0$. We assume in this section that we are away from the point where the expanding eigenvalues are equal. In section 5.2, we consider the case where the expanding eigenvalues are complex, but use a coordinate change which limits to the repeating eigenvalues case when the imaginary part of the complex pair vanishes.

5.1.3 Fixed point of the Poincaré map

Equations (15qtua) to (15qtuc) map a point on $H_{\mathrm{in}}^{A}$ to a point on $H_{\mathrm{in}}^{B}$. Due to the symmetry $g$ in equation (4), a fixed point of a full Poincaré return map will also be a fixed point of (15qtua) to (15qtuc). Fixed points of a full Poincaré return map can thus be found by dropping the $A$ and $B$ superscripts, and the dependence on [math] on the right hand side, to give the following four nonlinear equations, with four unknowns, $x_{c}$, $x_{e}$, $y_{e}$ and $T$:

[TABLE]

We substitute equations (15qtuyabacaebp) and (15qtuyc) into (15qtuya) to eliminate $x_{e}$ and $y_{e}$, which upon rearranging gives:

[TABLE]

Recall that $\lambda_{c}^{+}>0>\lambda_{c}^{-}$ and $\lambda_{e}^{+}>\lambda_{e}^{-}>0$. Since $T$ is assumed to be large, and is certainly positive, we can neglect the first ($1$) and second ($A_{1}$) terms on the left hand side, and the first ($A_{2}$) term on the right hand side, to get:

[TABLE]

We next substitute (15qtuyaa) into (15qtuyabacaebp) and (15qtuyc) and then finally into the expression for $T$ (15qtuyd), which we will then solve for $T$. This gives us:

[TABLE]

where $\Delta_{BC}=B_{1}C_{2}-C_{1}B_{2}$, and

[TABLE]

Note that when simplifying (15qtuyaba) to get (15qtuyabb) and (15qtuyabaca) to get (15qtuyabacb), terms in the numerator in $\rme^{(\lambda_{e}^{+}+\lambda_{c}^{-})T}$ and $\rme^{(\lambda_{e}^{-}+\lambda_{c}^{-})T}$, respectively, cancel out.

We substitute (15qtuyabb) and (15qtuyabacb) into (15qtuyd) to get:

[TABLE]

The final task is to solve (15qtuyabacad), which gives the period of a periodic orbit in the flow (the actual period is $3T$), close to the heteroclinic cycle, which corresponds to a fixed point in the map. For large $T$, the periodic orbit will be close to the heteroclinic cycle. We will do this in two different cases in sections 5.1.4 and 5.1.5.

Note that the left hand side of equation (15qtuyabacad) is the denominator in the equations for $x_{c}$, $x_{e}$ and $y_{e}$ (equations (15qtuyaa), (15qtuyabb) and (15qtuyabacb) respectively) so we substitute (15qtuyabacad) into these equations to simplify them, to get

[TABLE]

These three equations give the coordinates of the fixed point in terms of $T$. Note that in all three co-ordinates, the coefficient(s) of $T$ in the exponential is (are) negative, meaning the coordinates (of the fixed point) get smaller as $T$ gets larger, as would be expected.

We now check that $b(t)>0$ for all $t\in[0,T]$ for these solutions, namely that the periodic orbit corresponds to a positive travelling wave solution of (2). Note from (15a) and (15b) that $x_{e}^{A}(t)$ and $y_{e}^{A}(t)$ do not change sign, and in order to have $b(t)>0$ we require that $y_{e}^{A}(t)>x_{e}^{A}(t)$ for all $t\in[0,T]$. Writing (15a) and (15b) with initial conditions from (15qtuyabacaeb) and (15qtuyabacaec) gives us

[TABLE]

where $E_{1}=h\frac{\mathrm{sgn}(\Delta_{BC})}{\sqrt{A_{4}^{2}+A_{3}^{2}}}$, and so we require

[TABLE]

There are four cases to consider depending on the signs of $A_{4}E_{1}$ and $A_{3}E_{1}$, and the corresponding quadrant in $x_{e}^{A}$-$y_{e}^{A}$ space in which the solutions lie. Since we are only considering solutions that lie close to the heteroclinic connection, we assume in each case that the solutions $x_{e}^{A}(t)$ and $y_{e}^{A}(t)$ lie in the same quadrant as $\hat{\theta}_{e}^{A}$ .

If $A_{4}E_{1},A_{3}E_{1}>0$, then $y_{e}^{A}(t)>0$ and $x_{e}^{A}(t)<0$ and we are done. If $A_{3}E_{1},A_{4}E_{1}<0$, then $y_{e}^{A}(t)>0$ and $x_{e}^{A}(t)<0$ so solutions have $b(t)<0$ for all $t$ and this is not a positive travelling wave. If $A_{3}E_{1}<0<A_{4}E_{1}$, then $x_{e}^{A},y_{e}^{A}<0$ and (15qtuyabacaeah) gives us:

[TABLE]

Note that the left-hand side of the inequality is positive and the right-hand side is between [math] and $1$, and so we require $0<\tan\hat{\theta}_{e}^{A}<1$ for the solution to be positive. Since $x_{e}^{A},y_{e}^{A}<0$, putting these together means that $\pi<\hat{\theta}_{e}^{A}<5\pi/4$. Furthermore, solutions must satisfy

[TABLE]

which implies that $\lambda_{e}^{+}-\lambda_{e}^{-}$ must decrease to [math] as $T\rightarrow\infty$, as the heteroclinic bifurcation is approached.

Finally, suppose $A_{4}E_{1}<0<A_{3}E_{1}$, then $x_{e}^{A},y_{e}^{A}>0$ and (15qtuyabacaeah) gives:

[TABLE]

Again, the left-hand side of the inequality is positive and the right-hand side is between [math] and $1$, so $\tan\hat{\theta}_{e}^{A}>1$. With $x_{e}^{A},y_{e}^{A}>0$, this means that $\pi/4<\hat{\theta}_{e}^{A}<\pi/2$.

In summary, solutions will have $b(t)>0$ for all $t$ if the heteroclinic connection is such that $\pi/4<\hat{\theta}_{e}^{A}<\pi$. If $\pi<\hat{\theta}_{e}^{A}<5\pi/4$, then we can also find solutions with large $T$ with $b(t)>0$, so long as $\lambda_{e}^{+}-\lambda_{e}^{-}$ decreases to zero as $T$ tends to infinity. For other values of $\hat{\theta}_{e}^{A}$, periodic solutions close to the heteroclinic cycle will not correspond to positive travelling wave solutions of the PDEs (2).

In the following two sections, we consider two different cases depending on the relative size of the the two terms on the left-hand side of equation (15qtuyabacad).

5.1.4 Resonant bifurcation at $\lambda_{c}^{-}+\lambda_{e}^{-}=0$

In this section, we will show that a resonant-type heteroclinic bifurcation occurs when $\lambda_{c}^{-}+\lambda_{e}^{-}=0$.

Suppose that $A_{3}B_{1}\rme^{\lambda_{e}^{+}T}\gg A_{4}C_{1}\rme^{\lambda_{e}^{-}T}$. This will be the case if $A_{3},B_{1},A_{4}$ and $C_{1}$ are order 1, since $T$ is large and $\lambda_{e}^{+}>\lambda_{e}^{-}$. Then equation (15qtuyabacad) simplifies to

[TABLE]

or

[TABLE]

for $D_{1}=\frac{A_{3}B_{1}}{|\Delta_{BC}|\sqrt{A_{4}^{2}+A_{3}^{2}}}$. If $D_{1}<1$, then we see a branch of long-period periodic orbits emerging from the curve $\lambda_{c}^{-}+\lambda_{e}^{-}=0$ into the region where $\lambda_{c}^{-}+\lambda_{e}^{-}<0$. If $D_{1}>1$ then the solutions branch into $\lambda_{c}^{-}+\lambda_{e}^{-}>0$. This bifurcation curve can be seen in figure 3, where the black curve of heteroclinic bifurcations coincides with the light blue curve at $-\lambda_{c}^{-}=\lambda_{e}^{-}$. At this fixed point, taking the leading order term for $x_{c}$ in (15qtuyabacaea) gives

[TABLE]

This resonant bifurcation is unusual: usually you expect to see a resonant bifurcation when the contracting eigenvalue is equal to the leading expanding eigenvalue, that is, when $-\lambda_{c}^{-}=\lambda_{e}^{+}$ [18], but here it is $-\lambda_{c}^{-}=\lambda_{e}^{-}$.

Numerical simulations of periodic orbits close to the resonance bifurcation indicate that $x_{e}$ and $y_{e}$ are both positive ($\hat{\theta}_{e}^{A}\approx\pi/2$), and so from (15qtuyabacaeak), we must have $\hat{\theta}_{e}^{A}>\pi/4$ in order for solutions to have $b(t)>0$ for all $t$. Indeed, this is what we see in the numerical simulations.

We next assess whether we expect to see a ‘kink’ in the shape of the profile of the long-period solutions as the bifurcation point is approached. As can be seen in the time-series plots in figure 4, a kink is observed when there is a period of time during which the solution grows exponentially with rate $\lambda_{c}^{+}$. When the trajectory is near $\xi_{A}$, the contracting components are $c$ and $w$, which are linear combinations of $x_{c}^{A}$ and $y_{c}^{A}$, which grow/decay exponentially at rates $\lambda_{c}^{+}$ and $\lambda_{c}^{-}$ respectively. Observing a kink corresponds to having $x_{c}^{A}(t)>y_{c}^{A}(t)$ for some range of time $t$. Since $y_{c}^{A}$ is decaying and $x_{c}^{A}$ is growing, we will observe a kink if $|y_{c}^{A}(T)|\ll|x_{c}^{A}(T)|$. We have that

[TABLE]

At the resonant bifurcation, $\lambda_{c}^{-}=-\lambda_{e}^{-}$ and so $x_{c}^{A}(T)$ and $y_{c}^{A}(T)$ are the same order and so a kink won’t be observed in solutions. This is indeed what is observed, see panel (a) of figure 4.

In summary, we expect to find a resonant heteroclinic bifurcation with $-\lambda_{c}^{-}=\lambda_{e}^{-}$, that is, on the blue line in figure 3 with $\zeta>\sigma/2$, at the boundary between regions 1 and 2.

5.1.5 Orbit flip bifurcation at $A_{3}=0$

In this section we show that a branch of long-period periodic orbits can emerge when the heteroclinic cycle undergoes an orbit flip: that is, in the case when the heteroclinic connection is tangent to the strong unstable manifold. Recall that $\hat{\theta}_{e}^{A}$ gives the position at which the heteroclinic connection intersects $H_{\mathrm{out}}^{A}$. We have that $\tan(\hat{\theta}_{e}^{A})=-A_{3}/A_{4}$, and so as $A_{3}$ goes to zero, $\hat{\theta}_{e}^{A}$ goes to $\pi$, which corresponds to the heteroclinic connection being tangent to the strong unstable manifold (the $x_{e}^{A}$ axis; see figure 5), that is, a point of heteroclinic orbit flip.

We suppose that $A_{3}$ is small enough that the two terms on the right hand side of (15qtuyabacad) are of the same order, that is, neither can be discarded. We then rewrite equation (15qtuyabacad) as

[TABLE]

where we have assumed $A_{3}\ll A_{4}$ and so it can be dropped from the square root. Note that $A_{3}$ will only be small if the expressions in both exponentials are negative, namely if $\lambda_{c}^{-}+\lambda_{e}^{-}<0$, and then as $T$ goes to infinity, $A_{3}$ goes to zero. This holds in regions 1 and 5 of figure 3.

For fixed points in this case, we find the leading order term in $x_{c}$ to be

[TABLE]

Numerical simulations of periodic orbits close to the heteroclinic orbit flip bifurcation indicate that $x_{e}<0<y_{e}$, $\hat{\theta}_{e}^{A}$ is very close to (but just less than) $\pi$, and so we automatically satisfy the condition that $b(t)$ remains positive for all time.

To get a kinked solution, we again require that $|y_{c}^{A}(T)|\ll|x_{c}^{A}(T)|$. For solutions which start at the fixed point, this gives

[TABLE]

If $\lambda_{c}^{-}<-\lambda_{e}^{+}$, then $|y_{c}^{A}(T)|<|x_{c}^{A}(T)|$, and we will see a kinked solution. If $\lambda_{c}^{-}>-\lambda_{e}^{+}$, then $|y_{c}^{A}(T)|$ and $|x_{c}^{A}(T)|$ will be of the same order, and we will not observe a kink. However, we note that in order for solutions in this region to expand as much as they contract, we would instead observe a kink in the expanding phase, that is, a change in growth rate from $\lambda_{e}^{-}$ to $\lambda_{e}^{+}$.

The location of the orbit flip (if it exists at all) is determined by the global dynamics (that is, it can not be predicted by the eigenvalues). For equations (4), we find the location of the orbit flip numerically, by solving a boundary value problem to locate the heteroclinic orbit between $\xi_{A}$ and $\xi_{B}$, and insisting that the heteroclinic orbit is tangent to the strong unstable manifold at $\xi_{A}$. The location of the orbit flip is shown by a green curve in region 5 of figure 3. This green curve coincides with the black curve of heteroclinic bifurcations. In region 5, $\lambda_{c}^{-}<-\lambda_{e}^{+}$, and the periodic orbits close to this heteroclinic bifurcation do indeed show a kinked solution — see panel (c) of figure 4. We note that the orbit flip curve terminates on the curve where $\lambda_{e}^{-}=\lambda_{e}^{+}$ (the red curve in figure 3), which is to be expected, as equation (15qtuyabacaeaq) clearly does not generate large $T$ solutions at this point.

5.1.6 Saddle-node bifurcation of periodic orbits

Equation (15qtuyabacaeaq) gives the possibility of a saddle-node bifurcation between periodic orbits near the orbit-flip bifurcation. We compute:

[TABLE]

and set $\frac{dA_{3}}{dT}=0$ to find

[TABLE]

giving a branch of saddle-node bifurcations of periodic orbits at

[TABLE]

Recall that the orbit flip bifurcations only occur if $\lambda_{c}^{-}+\lambda_{e}^{-}<0$, so for the left hand side of equation (15qtuyabacaeaw) to be positive, we require $C_{1}>0$. The branch of saddle-node bifurcations can terminate in the branch of orbit flip bifurcations if the right hand side of (15qtuyabacaeax) becomes equal to zero. This can happen in a number of different ways, for instance, by the eigenvalue condition $-\lambda_{c}^{-}=\lambda_{e}^{+}$, or if one of the constants $A_{4}$ or $C_{1}$ become equal to zero. In figure 3 it appears that the first of these does not happen, and since $A_{4}$ and $C_{1}$ do not depend on the eigenvalues in an obvious way, we cannot say for sure what happens at the end of the branch of saddle-node bifurcations.

5.2 Complex eigenvalues

We now repeat the Poincaré map calculations in the region where the expanding eigenvalues are complex (regions 3 and 4). We make a different change of coordinates near $\xi_{A}$, and instead write

[TABLE]

In the new $x_{e}^{A},y_{e}^{A}$ coordinates, the local part of the flow becomes

[TABLE]

Note that in the limit as $\lambda_{e}^{I}\rightarrow 0$, the Jordan form of the linear part here becomes what one would use in the case of repeated eigenvalues. The solution to the local flow is

[TABLE]

We note again that in the limit $\lambda_{e}^{I}\rightarrow 0$, these solutions are exactly those that one would expect for the case with two repeated eigenvalues (in particular, the term $\sin(\lambda_{e}^{I}t)/\lambda_{e}^{I}$ limits to $t$).

We also define new coordinates for use near $\xi_{B}$:

[TABLE]

We define $r_{e}^{A}$ and $\theta_{e}^{A}$ as before, as in equation (11).

We use the same Poincaré sections, close to $\xi_{A}$ and $\xi_{B}$:

[TABLE]

for some $h\ll 1$.

The solution to equations (4) linearised around $\xi_{A}$ is now

[TABLE]

the local map then gives us

[TABLE]

where $T$ is again defined by

[TABLE]

We now note that in these coordinates, to ensure that $b(t)>0$ for all $t$, we will require that $y_{e}^{A}>0$ for all $t\in(0,T)$. We can write

[TABLE]

where

[TABLE]

Thus, in order for $b(t)$ to remain positive for all $t\in[0,T]$, a clear upper bound on $\lambda_{e}^{I}T$ is $\pi$ (because the $\sin$ changes sign with frequency $\pi$). Thus, $\lambda_{e}^{I}<\pi/T$, and since we are interested in solutions for which $T$ is large, $\lambda_{e}^{I}$ will be small.

The global part of the map doesn’t change, that is, we still have

[TABLE]

for some order one constants $F_{j}$.

Using (15qs), we can again write

[TABLE]

where $A_{3}$ and $A_{4}$ have the same expression as in the real part, namely $\tan\hat{\theta}_{e}^{A}=-A_{3}/A_{4}$ (although note that the values are different because the angle $\hat{\theta}_{e}^{A}$ is defined differently because of the different coordinate transformations made).

Substituting in for the right-hand side, we get

[TABLE]

Putting the global map (15qtuyabacaebg) together with the local map (15qtuyabacaeba) to (15qtuyabacaebd), using (15qtuyabacaebh), renaming the constants, and finallydropping the superscripts and the dependence on [math], gives us the following equations for the fixed points:

[TABLE]

There are again four unknowns, $x_{c}$, $x_{e}$, $y_{e}$ and $T$.

As in the case for real expanding eigenvalues, we again consider the magnitudes of the constants $B_{1}$, $B_{2}$, $C_{1}$, and $C_{2}$. Here, we find

[TABLE]

So, in this case, there are no degeneracies in these constants.

We now continue to find the fixed points of the return map. As noted earlier, we are interested in the limit when $T$ is large and hence $\lambda_{e}^{I}$ is small. To make the notation clear, we write $\epsilon=\lambda_{e}^{I}$. Recall that $0<\epsilon T<\pi$, and in particular, we make the ansatz

[TABLE]

for some order-one unknown $K$. We demonstrate below that this ansatz is correct. We can then write

[TABLE]

Again, we substitute the expressions (15qtuyabacaebj) and (15qtuyabacaebk) for $x_{e}$ and $y_{e}$ into the expression (15qtuyabacaebi) for $x_{c}$. This gives

[TABLE]

Rearranging gives

[TABLE]

The first term in the parentheses on the left hand side of equation (15qtuyabacaebo) is clearly smaller than the others. Dropping this term, and the terms of $O(\epsilon^{2})$ gives

[TABLE]

Substituting this into the expressions (15qtuyabacaebj) and (15qtuyabacaebk) for $x_{e}$ and $y_{e}$ gives, after some cancellation,

[TABLE]

where $\Delta_{AB}=A_{1}B_{2}-A_{2}B_{1}$. Similarly,

[TABLE]

where $\Delta_{AC}=A_{1}C_{2}-A_{2}C_{1}$. At this point, we note that the numerators and denominators of the expressions in (15qtuyabacaebp), (15qtuyabacaebq) and (15qtuyabacaebr) all contain one term which is multiplied by $\rme^{\lambda_{e}^{R}T}$, and one which is not. Since $\lambda_{e}^{R}>0$ and $T$ is large, we might think that the latter term is much smaller than the former and at lowest order, can be ignored. This is true for the numerators, since the term multiplying $\rme^{\lambda_{e}^{R}T}$ consists of $O(1)$ global constants which generically are non-zero. However, the expression multiplying $\rme^{\lambda_{e}^{R}T}$ in the denominators of these fractions contains the unknown constant $K$. It turns out that this expression is very small, and in fact, in the calculations below, we approximate it to lowest order by zero when finding $K$. Thus, both terms in the denominators must be kept.

Following this observation, we use equation (15qtuyabacaebl) to compute an expression for $T$, where we will ignore the terms not multiplied by $\rme^{\lambda_{e}^{R}T}$ in the numerators of both (15qtuyabacaebq) and (15qtuyabacaebr). First, we use (15qtuyabacaebq) and (15qtuyabacaebr) to compute

[TABLE]

Substituting into (15qtuyabacaebl) and rearranging, we get

[TABLE]

Both terms on the right-hand side are non-zero, but exponentially small (as $T$ is large), if $\lambda_{e}^{R}+\lambda_{c}^{-}<0$. If $\lambda_{e}^{R}+\lambda_{c}^{-}>0$, there are no solutions to this equation as there is nothing to balance the second term on the right-hand side, which would be exponentially large. Thus, we require $\lambda_{e}^{R}+\lambda_{c}^{-}<0$ (that is, we can’t be in region 3 of figure 3 and so must be in region 4), and so to lowest order, solutions will have

[TABLE]

or

[TABLE]

confirming that $K$ is order $1$. Thus our solution for $T$ is given by

[TABLE]

As noted above, in order for $b(t)$ to remain positive for all $t\in(0,T)$, we must have $T<\pi/\lambda_{e}^{I}$. For large $T$ solutions then, we require $\frac{A_{3}}{A_{4}}+\frac{C_{1}}{B_{1}}<0$. Since $A_{3},A_{4},B_{1}$ and $C_{1}$ are functions of the global dynamics (that is, the are not solely dependent on the eigenvalues of the equilibria), we cannot say where in parameter space this condition holds (apart from being within region 4, close to the boundary between regions 4 and 5).

We note that $K$ can change sign when $A_{4}$ passes through zero: this occurs when the heteroclinic connection is tangent to the positive $y_{e}$-axis. Since the coordinate changes we have used in the real and complex cases are different, this corresponds to the heteroclinic connection in the real case being tangent to the negative $x_{e}$-axis, which is exactly the point where the orbit-flip bifurcation curve terminates (see section 5.1.5). We thus expect a transition between orbit-flip and Belyakov–Devany-type bifurcation to occur, at a location determined by the global constants. This is consistent with what is observed in figure 3.

Again, we check to see whether we expect to see kinked solutions. Recall that to get a kinked solution, we require that $|y_{c}^{A}(T)|\ll|x_{c}^{A}(T)|$. Using equation (15qtuyabacaebs) in the denominator in the $x_{c}$ equation (15qtuyabacaebp) we can see that the denominator scales like $\rme^{(2\lambda_{e}^{R}+\lambda_{c}^{-})T}$. The numerator will be order $\rme^{\lambda_{e}^{R}T}$. Thus

[TABLE]

for some $O(1)$ constant $E_{1}$ and

[TABLE]

As noted above, $\lambda_{c}^{-}<-\lambda_{e}^{R}$, so $|y_{c}^{A}(T)|\ll|x_{c}^{A}(T)|$, and we expect to see a kinked solution, as observed (see panel (b) of figure 4).

In summary, we expect to find a Belyakov–Devaney-type heteroclinic bifurcation with $\lambda_{e}^{-}=\lambda_{e}^{+}$, that is, on the red curve in figure 3 with $\zeta<\sigma/2$, at the boundary between regions 4 and 5.

5.3 Summary

In summary, we conclude that heteroclinic bifurcations can only occur on the boundary between regions 1 and 2 (see figure 3), or on the boundary between regions 4 and 5, or within regions 1 or 5. All our numerical results (detailed below) point to all periodic orbits created in the Hopf bifurcation ending in heteroclinic bifurcations. If $\zeta>\sigma/2$, the heteroclinic bifurcation is of resonance type ($-\lambda_{c}^{-}=\lambda_{e}^{-}$), on the boundary between regions 2 and 1. If $\zeta<\sigma/2$, there are two possibilities: the heteroclinic bifurcation can either be of Belyakov–Devaney type (expanding eigenvalues changing from real to a complex-conjugate pair), on the boundary between regions 4 and 5, or of orbit flip type (when a constant in the global part of the map vanishes as the way in which the trajectories between equilibria change their orientation), within region 5 or region 1. The transition between the resonance and Belyakov–Devaney-type bifurcations occurs at $\zeta=\sigma/2$. The transition from Belyakov–Devaney-type to orbit flip occurs when a global coefficient changes sign and so cannot be deduced only from considerations of eigenvalues.

6 Further PDE simulations

In this section, we continue the numerical PDE simulations first discussed in section 2, and relate the results of these to the results of our calculations of the heteroclinic bifurcations. We begin by showing dispersion relations computed from the ODEs (4), which relate the period of the orbit to the parameter $\gamma$ (in the terminology of the ODEs), or equivalently, the wavelength of the travelling wave $\Lambda$, to the wavespeed $\gamma$ (in the terminology of the PDEs). The dispersion relations are computed in AUTO, by following periodic orbits created the Hopf bifurcation given in (8), as done in [14]. The resulting curves are shown in figure 7.

We ran numerical simulations of the PDEs (2) for values of $\zeta\in\{0.2,1.0,2.0,3.0\}$, on a periodic domain of size $500$. For each simulation, we started with small, randomly generated initial conditions, and integrated for a time period of $10,000$ to remove any transient behaviour. We then sample the solution at timepoints $t=10,000+100k$, for $k=1,\dots,40$. At each sample point, we compute the wavelengths and wavespeeds of the current solution profile. The wavelengths are computed by calculating the distances (in $x$) between points which have both $\log(a)=-1$, and $\frac{da}{dx}>0$. Wavespeeds at each of these points are computed by locally calculating $\frac{da}{dx}$ and $\frac{da}{dt}$ and using $\gamma=\frac{da}{dt}/\frac{da}{dx}$. Waves are only included in the analysis if each of the three variables $\log a$, $\log b$ and $\log c$ goes both above and below $-1$, over the wavelength of the wave. The results are plotted in figure 7(a).

The first thing to notice about these results is that there is a lot of scatter. This is for two main reasons. Firstly, although we can compute a ‘local’ wavespeed (i.e. a wavespeed for some specific point $(x,t)$), we cannot reliably compute a ‘local’ value of the wavelength. Secondly, in the simulations there are many different waves travelling both left and right (see figure 2); whenever the waves collide there is a region of time and space for which the wavespeed and wavelength are not well-defined, and our computations do not take account of this. However, for each of the values of $\zeta$ shown, it can be seen that there is a concentration of points along the AUTO-computed dispersion relation curve.

Further observations of numerical experiments indicate that for $\zeta=1.0$, solutions will often converge to a single travelling wave after sufficient time has passed (sometimes in excess of $t=50,000$). For the other values of $\zeta$ used in these experiments, we do not observe this convergence. In figure 7(b) we show the results of further similar computations for $\zeta=1.0$, but now we run multiple simulations from randomly chosen initial conditions, and sample the wavelengths and wavespeeds of the solution at a single timepoint, after the solutions have become close to a single travelling wave. Different initial conditions converge to travelling waves with different numbers of waves fitting into the box, but all of these lie very close to the AUTO-computed dispersion relation curves.

The values of $\zeta$ used above, together with $\sigma=3.2$, correspond to observing each of the three types of heteroclinic bifurcation discussed in section 5: the curve of heteroclinic bifurcations at $\zeta=0.2$ is of orbit flip type, at $\zeta=1.0$ is of Belyakov–Devaney type, and at $\zeta=2.0$ is of resonance type (see figure 3). We also include $\zeta=3.0$ to match with the data shown in figures 1 and 2. The heteroclinic bifurcation calculations we have done tell us about existence criteria for periodic orbits in the ODEs, which correspond to existence criteria for travelling waves in the PDEs. For instance, we can say that travelling wave solutions exist to the left of the heteroclinic bifurcation curve, and to the right of the Hopf bifurcation curve in figure 3, and this in turn gives us a maximum and minimum wavespeed for observed travelling waves. In order to be able to give firm predictions about whether these travelling waves would be observed in simulations, we would also need to understand the stability of the travelling waves, which is beyond the scope of this paper (but the subject of future work).

7 Bifurcation diagrams for varying $\sigma$

In this section we give some numerical results showing different bifurcation diagrams in the $(\gamma,\zeta)$ plane as $\sigma$ is varied. Most of the bifurcation curves were computed using AUTO [42]. Maintaining computational accuracy for periodic orbits close to heteroclinic cycles can be difficult for two reasons. Firstly, because the periodic orbits are of very long period, it is necessary to have a large number of mesh points defining the periodic orbit. Secondly, the heteroclinic connections lie in invariant planes where some of the coordinates are zero. The nearby periodic orbits thus will have coordinates which are very close to zero. In order to overcome the numerical issues associated with small numbers we make the following change of coordinates:

[TABLE]

and use differential equations for $A,B,C,U,V$ and $W$ in our numerical computations instead of the original equations. Since the periodic orbits which we are interested in exist entirely in the positive orthant (they correspond to positive travelling waves of (2)), we have no issues with taking the logarithm of a negative number. In AUTO, we compute a curve of periodic orbits which has a large, fixed, period ($T=300$ in the following calculations), and say that this curve well approximates the curve of heteroclinic bifurcations.

Figure 8 shows bifurcation diagrams of system (4) for various values of $\sigma$. For ease of comparison, we rescale $\zeta$ and $\gamma$ by writing $\hat{\zeta}=\zeta/\sigma$ and $\hat{\gamma}=\gamma/\sqrt{\sigma}$. Note that with this rescaling, all of the coloured lines (given by equations involving eigenvalues) in the bifurcation diagrams do not depend on $\sigma$. However, the location of the Hopf curve changes. The grey curve shows the Hopf bifurcation curve, as given by (8), and the black curve shows the heteroclinic bifurcation curve, as computed by AUTO. In (b), (c) and (d), a light grey curve, also computed by AUTO, shows a curve of saddle-node bifurcations of periodic orbits. The green curve is a curve of orbit-flip heteroclinic orbits, computed by solving a boundary value problem, as explained in section 5.1.5.

For all four values of $\sigma$ shown, the heteroclinic curve coincides with the curve $\lambda_{c}^{-}=\lambda_{e}^{-}$ (the light blue curve) for values of $\zeta$ greater than $\sigma/2$. For values of $\zeta$ below $\sigma/2$, there is a range of $\zeta=[\zeta^{*},\sigma/2)$ for which the heteroclinic curve coincides with the red curve, where the expanding eigenvalues are equal: the expanding eigenvalues are real to the right of this curve and complex to the left of this curve. Then for $\zeta<\zeta^{*}$, the heteroclinic curve coincides with the green curve: the curve of orbit flip heteroclinic orbits. We note that the transition point $\zeta^{*}$ is dependent on the global dynamics, and varies as $\sigma$ is varied. The curve of saddle-node of periodic orbits also appears to terminates at $\zeta=\zeta^{*}$.

In the lower panels of figure 8, we show zooms of each set of curves near to $\hat{\zeta}=0$, showing the orbit flip and saddle-node curves more clearly. We note that the numerical calculations become more difficult as $\sigma$ decreases, and for this reason we do not show the orbit flip or saddle-node curves on the panel for $\sigma=0.32$. In particular, we note that for the original Rock–Paper–Scissor equations with no diffusion (1), there is a degeneracy when $\sigma=0$, namely that the Hopf and heteroclinic resonant bifurcations are degenerate (the Jacobian matrix at the coexistence point has imaginary eigenvalues for all values of $\zeta$, and the heteroclinic orbit is at resonance for all values of $\zeta$). Something similar happens in this six-dimensional system: it is simple to shown that the Hopf bifurcation curve and the resonance bifurcation curve collapse onto one another as $\sigma$ is reduced to zero, and it appears numerically that in fact the whole heteroclinic bifurcation curve approaches the Hopf curve in this limit.

8 Discussion

We have summarised the results of our Poincaré map construction in section 5.3: the PDEs (2) have been shifted into a travelling frame of reference moving at speed $\gamma$ (4). In these sixth-order ODEs, travelling waves correspond to periodic solutions that are created in a Hopf bifurcation and, with increasing wavespeed, are destroyed in one of three types of heteroclinic bifurcation. Although the construction of the Poincaré map follows reasonably standard lines, there are technical issues: the unstable manifolds of the equilibria are four-dimensional, and we restated some standard definitions of heteroclinic cycles in order to accommodate (for example) positive contracting eigenvalues. We find it advantageous to delay solving for the period $T$ of the orbit until the very end, since due to cancellation of some exponential terms in the calculations, it isn’t obvious which terms can be safely neglected.

The periodic orbits we find can be kinked because one of the contracting eigenvalues is positive and the growth rate changes in magnitude but not sign at the transition from the contracting phase to the expanding phase. In addition, each of the three heteroclinic bifurcations is non-standard or new in some way. The resonance bifurcation, with $-\lambda_{c}^{-}=\lambda_{e}^{-}$, involves the leading expanding (that is, smallest positive) eigenvalue; usually it would be the non-leading (largest positive) expanding eigenvalue, that is, $-\lambda_{c}^{-}=\lambda_{e}^{+}$ [18]. The Belyakov–Devaney-type bifurcation, with the expanding eigenvalues changing from real to a complex-conjugate pair, and the orbit flip heteroclinic bifurcation, where the trajectories between equilibria change their orientation, are both new because they involve a robust (codimension zero) heteroclinic cycle, rather than a higher codimension homoclinic orbit [40, 41, 39].

It seems to be the case that stability conditions of heteroclinic cycles can be much more complicated than perhaps was thought several decades ago when the study of robust heteroclinic cycles was in its infancy. Much of this complexity perhaps arises in cases where unstable manifolds have dimensions greater than one. In the case in this paper, we have a positive contracting eigenvalue, and other types of stability are often seen when cycles have positive transverse eigenvalues (see, e.g. [20, 12, 31, 34]). Recently, the study of heteroclinic networks is receiving increasing attention in the literature: by definition, such networks must have at least one equilibrium with an unstable manifold of dimension greater than one. The stability of heteroclinic networks is almost certainly very subtle [43, 33], and we expect many interesting results in this area in the future.

Viewing the ODEs (4) as an Initial Value Problem, the heteroclinic cycles we describe are hopelessly unstable. However, the fixed points of the map correspond to travelling waves in the PDEs (2), and these may (viewed as a Boundary Value Problem on an appropriate periodic domain) be stable. We plan in future to use the results in this paper to address the stability of the travelling waves within the PDEs: intriguing preliminary results have shown that long-wavelength travelling wave solutions are stable when they bifurcate from the Belyakov–Devaney-type bifurcation curve, but unstable otherwise. In further work we will address the problem of existence and stability of spiral waves in the two-dimensional problem. Speculating further, it seems plausible that these ideas can be used to examine spiral waves in other systems, such as reaction–diffusion systems, or other spatially extended population models (e.g., Lokta–Volterra).

Acknowledgements

We are grateful to Graham Donovan, Cris Hasan, Mauro Mobilia and Hinke Osinga for helpful conversations, to Beverley White for reminding us to “think about the eigenvalues”, and to two anonymous referees for some very helpful reports. We acknowledge the London Mathematical Society for financial support through a Research in Pairs (Scheme 4) grant, and the hospitality of and financial support from the Department of Mathematics, University of Auckland.