A trio of heteroclinic bifurcations arising from a model of spatially-extended Rock-Paper-Scissors
AA trio of heteroclinic bifurcations arising from amodel of spatially-extended Rock–Paper–Scissors
Claire M Postlethwaite † § and Alastair M Rucklidge ‡ † Department of Mathematics, University of Auckland, Private Bag 92019, Auckland1142, New Zealand ‡ School of Mathematics, University of Leeds, Leeds LS2 9JT, UK
Abstract.
One of the simplest examples of a robust heteroclinic cycle involves threesaddle equilibria: each one is unstable to the next in turn, and connections from oneto the next occur within invariant subspaces. Such a situation can be described by athird-order ordinary differential equation (ODE), and typical trajectories approacheach equilibrium point in turn, spending progressively longer to cycle around thethree points but never stopping. This cycle has been invoked as a model of cycliccompetition between populations adopting three strategies, characterised as Rock,Paper and Scissors. When spatial distribution and mobility of the populations is takeninto account, waves of Rock can invade regions of Scissors, only to be invaded by Paperin 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. Inthis paper, we explore how the robust heteroclinic cycle in the ODE manifests itself inthe 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 createdin a Hopf bifurcation and are destroyed in three different heteroclinic bifurcations,depending on parameters, as the travelling wave approaches the heteroclinic cycle. Weexplore the three different heteroclinic bifurcations, none of which have been observedin the context of robust heteroclinic cycles previously. These results are an importantstep towards a full understanding of the spiral patterns found in two dimensions,with possible application to travelling waves and spirals in other population dynamicsmodels.
Nonlinearity (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, andPaper wraps Rock, provides an appealing simple model of cyclic competition betweendifferent strategies or species in evolutionary game theory and biology [1, 2]. The gamehas been invoked as a description of three competing species of
E. coli [3] and of threecolour-variants of side-blotched lizards [4], but the idea of cyclic competition has arisenalso in rotating convection [5] and as the simplest example of a heteroclinic cycle [6]. § Corresponding author ([email protected]) a r X i v : . [ n li n . PS ] A p r eteroclinic bifurcations in Rock–Paper–Scissors A ( t ), B ( t ) and C ( t ) be the numberof people playing Rock, Paper or Scissors at any moment of time t . Pairs of people aredrawn at random and when they play, either it is a tie (if they are drawn from the samegroup), or one beats the other. In this case, the loser can either adopt the strategy of thewinner (dominance–replacement) or the loser can withdraw from the game (dominance–removal). Once removed, players are replaced (up to a maximum number N ) and areassigned to Rock, Paper or Scissors with probabilities proportional to the number ofRock, 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 ofRock players is small), they can never return, which means that Scissors would have nocompetitors 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 modelledby three ordinary differential equations (ODEs) [1, 8, 9]:˙ a = a (1 − ( a + b + c ) − ( σ + ζ ) b + ζc ) , ˙ b = b (1 − ( a + b + c ) − ( σ + ζ ) c + ζa ) , (1)˙ c = c (1 − ( a + b + c ) − ( σ + ζ ) a + ζb ) , where a ( t ) = A/N , b ( t ) = B/N , c ( t ) = C/N , and σ and ζ are non-negative parametersthat control the rates of dominance–removal and dominance–replacement respectively,scaled to the rate of replacement. We have assumed symmetry between Rock, Paperand 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 trivialsolution ( a, b, c ) = (0 , , , , , ,
0) and (0 , , a, b, c ) = σ (1 , , σ > ζ >
0, this system of ODEshas solutions that approach each of three on-axis equilibria in turn, taking progressivelylonger to cycle around the three points but never stopping [1] (in contrast to eventualfixation in the discrete case).This gradual slowing down of trajectories as they spend longer and longer neara 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 theon-axis equilibria: these are ζ and − ( σ + ζ ), and the amount of time taken for eachcycle is a factor of σ + ζζ longer than the previous one [10]. In this expression it isapparent that allowing either ζ = 0 or σ = 0 requires special attention. The situationwhere the eigenvalue ratio is equal to 1 ( σ = 0, ζ >
0) is normally called a resonancebifurcation from the heteroclinic cycle, associated with the creation of a long-periodperiodic orbit [11, 12]. However, in the ODEs (1), letting σ = 0 is degenerate, inthat the coexistence equilibrium has pure imaginary eigenvalues and the ODEs have aninvariant plane a + b + c = 1 on which there is a continuous family of nested periodicorbits parameterised by abc = constant. eteroclinic bifurcations in Rock–Paper–Scissors N and small lattice spacing, the dynamics is describedby the partial differential equations (PDEs) [8, 9]:˙ a = a (1 − ( a + b + c ) − ( σ + ζ ) b + ζc ) + ∇ a, ˙ b = b (1 − ( a + b + c ) − ( σ + ζ ) c + ζa ) + ∇ b, (2)˙ c = c (1 − ( a + b + c ) − ( σ + ζ ) a + ζb ) + ∇ c, where the spatial coordinates ( x, y ) are scaled so that the diffusion constants (assumedto be equal) are equal to 1. Typically the PDEs are solved with periodic boundaryconditions. The spatial mobility allows for persistent spiral-like or turbulent patternsof Rock, Paper and Scissors [13], in which regions dominated by Rock invade regions ofScissors, which invade regions of Paper, which in turn invade regions of Rock. In thecase 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-dimensionalperiodic travelling wave (TW) solutions of the PDEs (2) [14].The central question we address in this paper is: what is the connection betweentravelling waves in the PDEs (2) and heteroclinic cycles? The TWs are periodic orbits ina moving frame of reference, and, taking the wavespeed as a parameter, these periodicorbits originate in a Hopf bifurcation and end when they collide with a heterocliniccycle [14]. In this paper we find conditions under which TWs with arbitrarily longwavelength can exist as solutions of (2), close to a heteroclinic cycle in the sixth-orderODEs that describe the dynamics in the travelling frame. We find that there are threedifferent ways in which this can happen: • there can be a resonance bifurcation from the heteroclinic cycle in the sixth-orderODEs, 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 partof a pair of complex eigenvalues vanishes; or • there can be a bifurcation of orbit flip type, at which there is a change in the wayin which the trajectories between equilibria are oriented.Although our analysis proceeds along reasonably standard lines, there are severalunusual aspects, and the calculations are challenging, not least because the unstablemanifolds of the equilibria in the heteroclinic cycles are of high dimension. It turns outthat each of these three bifurcations is non-standard and, to our knowledge, has not beenobserved in the context of heteroclinic cycles before. We are able to find conditions underwhich each of these three bifurcations occurs, and, to some extent, how the transitionfrom one type to the next occurs. Our results give a much clearer picture of the originof the one-dimensional TW solutions of the PDEs (2), a first and necessary step inunderstanding their stability, which in turn is necessary for understanding the stabilityof the two-dimensional spiral solutions of (2). eteroclinic bifurcations in Rock–Paper–Scissors Figure 1.
Snapshots of numerical solutions of equations (2), in two spatial dimensionswith parameters σ = 3 .
2, and (a) ζ = 1 .
0, (b) ζ = 2 .
0, (c) ζ = 3 .
0. The domain sizefor the integrations was 500 × a , b and c are dominant are shownin red, green and blue respectively. The central spiral rotates clockwise with the threecolours moving outwards. The remainder of this paper is structured as follows. We begin in section 2by reviewing some numerical results from [14] and showing some simulations of thePDEs (2), both in one and two spatial dimensions. We also relate properties of thetravelling wave solutions of the PDEs to periodic solutions of a related set of ODEs. Insection 3 we review the definitions of heteroclinic cycles and summarise what is alreadyknown about ways in which they can bifurcate. We also compare these bifurcationswith those seen near homoclinic orbits, and relate these to the new bifurcations we havefound. In section 4 we describe the derivation of the ODEs we will be studying for theremainder of the paper. Then in section 5 we derive a Poincar´e map which describesthe flow close to the heteroclinic cycle in the ODEs. This section contains a lot ofcalculation but the results are summarised at the start and end of the section for thereader who doesn’t wish to delve into too many of the gritty details. In section 6 we givesome further numerical results from simulation of the PDEs for a range of parametervalues, and finally in section 7 we look at numerically computed bifurcation diagramsas the parameter σ is varied and discuss the limit σ →
0. Section 8 concludes.
2. PDE simulations
We begin with the PDEs for the spatially-extended Rock–Paper–Scissors model asgiven in equations (2). In figure 1 we show numerical results from the integrationof equations (2) in two spatial dimensions, from [14]. A variety of behaviours can beobserved, but of particular interest are the spiral-type solutions. When a slice is takenradially through the centre of a spiral, the profile of the solution in the outer part ofthe spiral resembles a travelling wave in one spatial dimension.Figure 2 shows the results of numerical integration of equations (2) in one spatialdimension, in a large box of size 500, for σ = 3 . ζ = 1 . , . , .
0. Initial conditionsare 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 eteroclinic bifurcations in Rock–Paper–Scissors (c) (b) (a) Figure 2.
The figures show results from numerical integration of equations (2) in onespatial dimension, in a box of size 500 with periodic boundary conditions. The lefthand column shows time-space plots: time is plotted horizontally and space vertically.Areas in which a , b and c are dominant are shown in red, green and blue respectively.The right hand column shows snapshots at t = 1000. Parameters are σ = 3 .
2, and (a) ζ = 1 .
0, (b) ζ = 2 .
0, (c) ζ = 3 . conditions after a short transient. For all three values of ζ , travelling waves of differentdirections, wavespeeds and wavelengths are evident. In the simulation for ζ = 1 . t = 500, waves consistently travel to the left, and eventually (after beingintegrated for a longer time period than shown here), this solution has six waves of equalwavelengths (and equal wavespeed) fitting in the periodic box. For the larger values of ζ , the solutions appear more complicated, in particular, faster wavespeeds and smallerwavelengths 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 thePDEs (2); that is, we would like to be able to say whether solutions will eventuallyasymptote onto a single travelling wave, and what the wavespeed and wavelength ofthat travelling wave will be. In order to do this, we would need to know both existenceand stability criteria, as well as have information about the basins of attraction of thetravelling waves. The latter two are difficult problems, and are beyond the scope of this eteroclinic bifurcations in Rock–Paper–Scissors
3. Review of heteroclinic cycles and bifurcations
Before we begin the calculations, in this section we first include a review of heterocliniccycles, 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 findthat this is useful nonetheless.Consider a system of ordinary differential equations˙ x = f ( x ) , x ∈ R N . (3)Then we have: Definition 1 A heteroclinic cycle is a finite collection of equilibria { ξ , . . . , ξ n } of (3),together with a set of heteroclinic connections { φ ( t ) , . . . , φ n ( t ) } , where φ j ( t ) is a solutionof (3) such that φ j ( t ) → ξ j as t → −∞ and φ j ( t ) → ξ j +1 as t → ∞ , and where ξ n +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 parametervalues, that is, they are of codimension zero, and are referred to as ‘robust’ [6,16,17]. Inthe work of Krupa and Melbourne [15, 18] and others (e.g., [19–26]), robust heterocliniccycles 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 existfor open sets of parameter values due to a combination of invariant subspaces and thedimensions of stable and unstable manifolds of equilibria for the flow restricted to theseinvariant subspaces . An additional difference in our work is that the invariance of thesubspaces is not forced by symmetry, but instead by the invariance of extinction incontinuous-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 ξ j and ξ j +1 . Let W u | P j ( ξ j ) and W s | P j ( ξ j +1 )by the unstable manifold of ξ j and stable manifold of ξ j +1 for the flow restricted to P j . eteroclinic bifurcations in Rock–Paper–Scissors Table 1.
Classification of eigenvalues. P (cid:9) L denotes the orthogonal complement in P of the subspace L .Eigenvalue class SubspaceRadial ( r ) L j ≡ P j − ∩ P j Contracting ( c ) V j ( c ) = P j − (cid:9) L j Expanding ( e ) V j ( e ) = P j (cid:9) L j Transverse ( s ) V j ( s ) = ( P j − + P j ) ⊥ Then, if dim( W u | P j ( ξ j )) + dim( W s | P j ( ξ j +1 )) > dim( P j ), then a heteroclinic connectionfrom ξ j to ξ j +1 will be codimension zero, this is, it will persist under small changes tothe 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 ξ , . . . , ξ n , whererobust here means codimension zero.We further define L j ≡ P j − ∩ P j and clearly ξ j ∈ L j . Following [17], the eigenvaluesof the linearisation of f ( x ) about each equilibrium can be classified according to thesubspaces in which the eigenspaces lie, as shown in table 1. As we will discuss in thefollowing, because we do not require that P j contains the unstable manifold of ξ j (unlikein the definition used by Krupa and Melbourne [15]), we are allowed to have positiveradial and/or contracting eigenvalues.Methods for determining the stability properties of an isolated heteroclinic cycleare in principal well-established [11, 15, 18, 20, 27–34]: that is, one can construct aPoincar´e map, by linearising the flow around the fixed points and the heteroclinicconnections. Many examples have been investigated in lower dimensions ( R and R in particular), but in higher dimensions, calculations can become quite intricate. Anumber of codimension-one bifurcations have been identified in which the stability ofrobust heteroclinic cycles changes, but issues of stability turn out to be more subtle thanmight be at first thought (for several examples, see [12, 27, 29, 31]). Two well-studiedways in which heteroclinic cycles can change stability are resonance and transverse bifurcations. A resonance bifurcation [11, 15, 23, 26, 28] occurs when an algebraiccondition 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 periodicorbit. In a transverse bifurcation from a heteroclinic cycle [20], a local bifurcation causesa transverse eigenvalue of one of the equilibria in the cycle to change sign. This canresult in a bifurcating periodic orbit or heteroclinic cycle, depending on the specificsituation.In this paper, we use the standard methods for analysing the dynamics close to aheteroclinic cycle, namely, we construct a Poincar´e map which approximates the flowof the differential equations close to the heteroclinic cycle, but as mentioned in theintroduction, 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 computeconditions for the existence of nearby periodic orbits. We find that long-period periodic eteroclinic bifurcations in Rock–Paper–Scissors
4. Derivation of ODEs and existence of heteroclinic cycles
In this paper, we examine the behaviour of the travelling wave solutions in onedimension, and so we consider equations (2) with only one spatial dimension, so ∇ = ∂ ∂x . We move to a travelling frame with wavespeed γ >
0, so define z = x + γt ,then ∂∂x → ∂∂z and ∂∂t → γ ∂∂z + ∂∂t . This results in the following set of PDEs in thetravelling frame: ∂a∂t + γ ∂a∂z = a (1 − ( a + b + c ) − ( σ + ζ ) b + ζc ) + ∂ a∂z ,∂b∂t + γ ∂b∂z = b (1 − ( a + b + c ) − ( σ + ζ ) c + ζa ) + ∂ b∂z ,∂c∂t + γ ∂c∂z = c (1 − ( a + b + c ) − ( σ + ζ ) a + ζb ) + ∂ c∂z . Travelling wave (TW) solutions in the moving frame have ∂∂t = 0. We thus set ∂∂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 setof six first-order ODEs: a z = u,u z = γu − a (1 − ( a + b + c ) − ( σ + ζ ) b + ζc ) ,b z = v,v z = γv − b (1 − ( a + b + c ) − ( σ + ζ ) c + ζa ) , (4) c z = w,w z = γw − c (1 − ( a + b + c ) − ( σ + ζ ) a + ζb ) . eteroclinic bifurcations in Rock–Paper–Scissors a , b and c are non-negative, we define a positive travelling wave as a periodicsolution of (4) with a, b, c > z . In an abuse of notation, the independentvariable z will be referred to as ‘time’ (and denoted with a ‘ t ’) when we constructPoincar´e maps in the next section.Let x = ( a, u, b, v, c, w ), and note that the coexistence and on-axis equilibria of (4)correspond to x = σ (1 , , , , , x = (1 , , , , , x = (0 , , , , ,
0) and x = (0 , , , , , ξ H , ξ A , ξ B and ξ C respectively. Alsonote that the ODEs (4) are invariant under the rotation symmetry g : g ( a, u, b, v, c, w ) = ( b, v, c, w, a, u ) . (5)The Jacobian matrix at ξ A is J A = γ σ + ζ − ζ
00 0 0 1 0 00 0 − ζ γ σ + ζ γ (6)The eigenvalues of J A are given in table 2. Note that we frequently refer to ‘theeigenvalues of ξ A ’, by which of course we mean the eigenvalues of J A . By the symmetry g , ξ B and ξ C have the same eigenvalues.Let the four-dimensional subspace { c = w = 0 } be labelled P ( ξ A ). It can easilybe seen that P ( ξ A ) is invariant under the flow of (4). For the dynamics restricted to P ( ξ A ), ξ A has a three-dimensional unstable manifold, and ξ B has a two-dimensionalstable manifold. By dimension counting, it is reasonable to expect that these manifoldswill intersect, and hence that there is a heteroclinic connection within P ( ξ A ) between ξ A and ξ B , which persists under small perturbations. We are able to numerically confirmthe existence of a heteroclinic connection for a wide range of parameter values. Bysymmetry, there is thus a robust heteroclinic cycle between ξ A , ξ B and ξ C .As discussed earlier, because our definition of robust heteroclinic cycle did notrequire the unstable manifold of ξ A to be contained in P ( ξ A ), we can have radialor contracting eigenvalues that have positive real part, and in fact, this is what wefind (see table 2). Specifically, we note that λ − c < < λ + c , λ − r < < λ + r , and0 < Re ( λ − e ) ≤ Re ( λ + e ).The Jacobian matrix at ξ H is: J H = σ γ σ + ζ σ − ζ σ
00 0 0 1 0 0 − ζ σ σ γ σ + ζ σ
00 0 0 0 0 1 σ + ζ σ − ζ σ σ γ (7) eteroclinic bifurcations in Rock–Paper–Scissors λ ± r = (cid:16) γ ± (cid:112) γ + 4 (cid:17) Contracting λ ± c = (cid:16) γ ± (cid:112) γ + 4( σ + ζ ) (cid:17) Expanding ( γ − ζ > λ ± e = (cid:16) γ ± (cid:112) γ − ζ (cid:17) Expanding ( γ − ζ < λ ± e = λ Re ± iλ Ie = (cid:16) γ ± i (cid:112) ζ − γ (cid:17) Table 2.
Eigenvalues of the equilibrium ξ A in (4). J H has pure imaginary eigenvalues ± iω H when γ = γ H ( σ, ζ ), where γ H ( σ, ζ ) ≡ √ σ + 2 ζ ) (cid:112) σ ( σ + 3) , and ω H = σ σ + 3) , (8)at which point a Hopf bifurcation creates periodic orbits of period (in the z variable) Λ H = πω H . Numerical analysis of equations (4) with AUTO [14,42] show thatthe branch of periodic orbits grows in period as γ is increased from the Hopf bifurcation,eventually ending in a heteroclinic bifurcation. The Hopf and heteroclinic bifurcationcurves can be seen in Figure 3 for σ = 3 . ζ and σ . For the value of σ used to create figure 3 ( σ = 3 . ζ > σ/ .
6, the heteroclinic bifurcation is of resonance type, and occurswhen − λ − c = λ − e (where the black curve coincides with the blue curve in figure 3); (b) if ζ ∗ < ζ < σ/ . ζ ∗ ≈ . λ Ie = 0 (where the black curve coincides with the redcurve in figure 3), and (c) if 0 < ζ < ζ ∗ , then the heteroclinic bifurcation is of orbit fliptype (where the black curve coincides with the green curve in figure 3).In the first two cases, the Hopf and heteroclinic bifurcation curves denote theexistence boundaries of periodic orbits in the ODEs, and hence also of travelling wavesin the PDEs. Specifically, the Hopf bifurcation curve indicates the minimum wavespeed γ (and minimum wavelength, given by Λ H as written after equation (8)), and theheteroclinic bifurcation curve indicates the maximum wavespeed. That is, for ζ > σ/ γ H ( σ, ζ ) < γ < (cid:114) σ ζ + (cid:114) σ , where γ H is given in (8). For ζ ∗ < ζ < σ/ .
6, the allowed wavespeeds are γ H ( σ, ζ ) < γ < (cid:112) ζ. eteroclinic bifurcations in Rock–Paper–Scissors γζ H o p f : γ = γ H λ − e = λ + e − λ − c = λ − e − λ − c = λ + e − λ − c = λ R e ζ ∗ SN (cid:27) Figure 3.
Bifurcation diagram for the ODEs (4), in ( γ, ζ ) parameter space, with σ = 3 .
2. The blue line ( ζ = (cid:112) σ γ − σ ) and red curve (4 ζ = γ ) are tangent at( γ, ζ ) = ( √ σ, σ/ σ + ζ ) = 3 γ ). These threecurves divide the parameter space into five regions, labelled by blue numbers, anddefined in table 3. The green curve is the locus of a heteroclinic orbit flip. The dashedgrey line is a curve of Hopf bifurcations (equation (8)). Periodic orbits bifurcate tothe right of this line and disappear in a curve of heteroclinic bifurcations (black). Acurve of saddle-node bifurcations of periodic orbits (light grey) exists for smaller ζ .The inset shows a zoom near the saddle-node of periodic orbits (SN) and heteroclinicorbit flip (green) bifurcations. Table 3.
Definitions of the regions of parameter space shown in Fig. 3 and eigenvalueproperties therein.
Region Definition Eigenvalue properties1 ζ < (cid:112) σ γ − σ λ ± e ∈ R , λ − e < | λ − c | < λ + e ζ > σ , (cid:112) σ γ − σ < ζ < γ λ ± e ∈ R , | λ − c | < λ − e < λ + e γ < ζ < γ − σ λ ± e ∈ C , | λ − c | < λ Re ζ > γ , γ − σ < ζ λ ± e ∈ C , λ Re < | λ − c | ζ < σ , (cid:112) σ γ − σ < ζ < γ λ ± e ∈ R , λ − e < λ + e < | λ − c | For ζ < ζ ∗ , the heteroclinic bifurcation is of orbit-flip type, and there also exists abranch of saddle-node bifurcations of periodic orbits (light grey curve). Here, the righthand boundary for existence of travelling waves is the saddle-node bifurcation curve,not the heteroclinic bifurcation curve. The location of both of these curves dependson global parameters, so here we cannot give an explicit expression for the maximumwavespeed.In figure 4 we show time-series of periodic solutions of (4) which are close to the eteroclinic bifurcations in Rock–Paper–Scissors λ − c λ − e λ − c λ + c λ Re λ − c λ + c λ + e z l og a , l og b , l og c l og a , l og b , l og c l og a , l og b , l og c (a)(b)(c) Figure 4.
The figures show time series (in logarithmic coordinates) of periodicsolutions to the ODEs (4), computed using AUTO, near the heteroclinic bifurcation.The coordinates a , b and c are shown in red, blue and green respectively, all inlogarithmic coordinates. Parameter values are σ = 3 .
2, and (a) ζ = 3 . γ = 3 . ζ = 1 . γ = 2 . ζ = 0 . γ = 1 . three types of heteroclinic bifurcations. The examples are all right-travelling waves (inthe PDE setup); left-travelling waves are also possible. In panel (a), we show a periodicorbit close to the heteroclinic resonance bifurcation (near the edge of region 2, theexpanding eigenvalues are real). The slopes in the contracting and expanding phasescan be seen to be very close to λ − c and λ − e . In panel (b), we show a periodic orbit closeto the heteroclinic Belyakov–Devaney-type bifurcation (near the edge of region 4, theexpanding eigenvalues are complex). Here, λ Ie is very close to zero, and the slope in eteroclinic bifurcations in Rock–Paper–Scissors λ Re . In the contracting phase, we see slopes equal toboth the negative contracting eigenvalue, λ − c , and the positive contracting eigenvalues λ + 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 phaseis λ + e , because the periodic solution lies close to the heteroclinic orbit, which is closeto tangent to the strong unstable manifold. Again, in the contracting phase we seeboth the positive and negative slopes. We refer later to periodic orbits which have botha positive and negative slope in the contracting phase as those having a kink – thekink refers to the change in growth rate at the transition from the contracting to theexpanding phase.In long-period orbits such as in Figure 4, the total amount of decay in thecontracting phase must balance the growth in the expanding phase; the contractingand 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 2and 3: λ − 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 enoughdecay to balance the growth.A further point to note is that not all periodic solutions of (4) correspond to positivetravelling wave solutions of (2). In particular, because we are considering a populationmodel, we will start with initial conditions (of (2)) which have a, b, c ≥
0, and, givenreasonable conditions on the smoothness of the initial conditions, it can be shown that a, b, c ≥ t (in (2)). Only periodic solutions of (4) which have a, b, c > z correspond to positive travelling wave solutions of (2). This may be important, becausethe variables in (4) may change sign along the heteroclinic connections — clearly theywill in the case that the expanding eigenvalues are complex.
5. Constuction of a Poincar´e map and analysis of heteroclinic bifurcations
In this section we construct a Poincar´e map which approximates the dynamics closeto the heteroclinic cycle of equations (4) as described in section 4. We treat thecases in which the expanding eigenvalues are real and complex separately, althoughthe computations are quite similar. Regions of real (1, 2 and 5) and complex (3 and4) eigenvalues are divided by the red curve in figure 3. In this section, we refer to theindependent variable of equations (4) as time ( t ) rather than z .The Poincar´e map we construct here will follow a trajectory that starts on anincoming section near ξ A and ends on an incoming section near ξ B . In both realand complex cases, we define Poincar´e sections close to ξ A and ξ B , and derive a localmap which approximates the flow close to ξ A . We combine this with a global maplinearised about the location of the heteroclinic connection from ξ A to ξ B and then usethe symmetry g to map the coordinates back to a Poincar´e section close to ξ A . We areable to disregard the radial directions in our computations. This is often done becauseeigenvalues in the radial directions are negative, with an invocation to an invariant eteroclinic bifurcations in Rock–Paper–Scissors T that the trajectory spends close to the equilibria, but leave this as an unknown definedimplicitly in terms of the other coordinates in the map: the map is defined in termsof three coordinates, and the time T . It then becomes possible to solve the equationsfor fixed points of the Poincar´e 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 becomelarge will give us the locations of the heteroclinic bifurcations.In the case where the expanding eigenvalues are real, we are able to find two differenttypes of solution for large T , depending on which terms in this equation dominate: onetype of solution generates the resonant bifurcation, and the other generates the orbit-flipbifurcation. In the case where the expanding eigenvalues are complex, we find only onetype of solution, corresponding to a Belyakov–Devaney-type bifurcation.The period of the bifurcating periodic orbit scales differently with the distancefrom the bifurcation point, depending on the type of bifurcation. Suppose that µ is aparameter which measures the distance from the heteroclinic bifurcation in each of thethree cases, then: (a) in the resonance bifurcation, µ ∝ | λ − e + λ − c | , and T scales like1 / | µ | (see equation (39)); (b), in the Belyakov–Devaney-type bifurcation µ ∝ | λ Ie | , and T scales like 1 / | µ | (see equation (72)); and (c) in the orbit flip bifurcation, µ ∝ A , aglobal constant which determines the angle at which the heteroclinic connection exits aneighbourhood of ξ A , and T scales like log | µ | (see equation (43)).In each case, once we have computed an expression for the fixed points of thePoincar´e map, we also check that the corresponding periodic orbits satisfy the conditionthat a, b, c > ξ A : during the transition between equilibria the coordinates will be order 1 and hencewill not change sign. In a neighbourhood of ξ A , it is clear that a ( t ) will not change sign,as it is order 1. The heteroclinic connection leaving ξ A lies in an invariant subspacewhich has c = 0, so c cannot change sign during the transition from ξ A to ξ B . Thus thecoordinate which will need to be checked is b ( t ).Finally, we will check whether or not we expect the solution to be ‘kinked’. To begin, we define new coordinates which we use when the trajectory is near ξ A x Ae = λ − e b − v, y Ae = λ + e b − v, x Ac = λ − c c − w, y Ac = λ + c c − w. (9)Recall that we are interested in solutions which have b ( t ) >
0, which in these coordinates,means we must have y Ae > x Ae . The coordinates in (9) are aligned with the eigenvectors eteroclinic bifurcations in Rock–Paper–Scissors ξ A can be written dx Ae dt = λ + e x Ae , dy Ae dt = λ − e y Ae , dx Ac dt = λ + c x Ac , dy Ac dt = λ − c y Ac . (10)We will also make use of polar coordinates in the expanding directions, namely r Ae and θ Ae , defined by( r Ae ) = ( x Ae ) + ( y Ae ) and tan θ Ae = y Ae x Ae . (11)The constraint y Ae > x Ae means that π/ < θ Ae < π/
4. We similarly define newcoordinates for use near ξ B : x Be = λ − e c − w, y Be = λ + e c − w, x Bc = λ − c a − u, y Bc = λ + c a − u. (12)We further write x A = ( x Ae , y Ae , x Ac , y Ac ) and x B = ( x Be , y Be , x Bc , y Bc ).We define Poincar´e sections, close to ξ A and ξ B : H A in = { x | y Ac = h } H A out = { x | r Ae = h } H B in = { x | y Bc = h } for some h (cid:28) ξ A and a global map from ξ A to ξ B asfollows. Let the time it takes the trajectory to travel from H A in to H A out be T . The localmap is Π loc : H A in → H A out x A ( T ) = Π loc ( x A (0))( x Ae ( T ) , y Ae ( T ) , x Ac ( T ) , y Ac ( T )) = Π loc ( x Ae (0) , y Ae (0) , x Ac (0) , h ) , where x Ae ( T ) + y Ae ( T ) = h , and the global map isΠ glo : H A out → H B in x B = Π glo ( x A ( T ))( x Be , y Be , x Bc , h ) = Π glo ( x Ae ( T ) , y Ae ( T ) , x Ac ( T ) , y Ac ( T ))where again x Ae ( T ) + y Ae ( T ) = h . In figure 5 we show a schematic of the expandingdynamics near ξ A .We label the heteroclinic connection between ξ A and ξ B as γ AB . Recall that theunstable manifold of ξ A , W u ( ξ A ), is four-dimensional. The heteroclinic connection is aone-dimensional sub-manifold of W u ( ξ A ). In addition, we also know that the connectionlies in the invariant subspace P ( ξ A ) (which has c = w = 0, equivalently, x Ac = y Ac = 0near ξ A or x Be = y Be = 0 near ξ B ). We consider the points at which the heteroclinicconnection intersects the Poincar´e sections, and write γ AB ∩ H A out = ˆ x A = (ˆ x Ae , ˆ y Ae , , , γ AB ∩ H B in = ˆ x B = (0 , , , h ) (13)where (ˆ x Ae ) + (ˆ y Ae ) = h, ˆ y Ae ˆ x Ae = tan ˆ θ Ae . (14) eteroclinic bifurcations in Rock–Paper–Scissors y Ae x Ae ˆ θ Ae H A out x A ( T ) x A (0) γ AB ˆ x A b = 0 Figure 5.
The figure shows a schematic of the expanding subspace from ξ A , in thecase when the expanding eigenvalues are real. The bold line indicates the heteroclinicconnection γ AB , and it intersects the Poincar´e section H A out (shown by a dotted curve)at ˆ x A . A trajectory close to the the heteroclinic connection is shown, starting at apoint x A (0) and hitting H A out at x A ( T ) (both points marked with black squares). Thegrey line indicates where b = 0; b > The x Bc coordinate of ˆ x B is zero because γ AB must lie in the stable manifold of ξ B , and x Bc is the coordinate associated with the positive contracting eigenvalue, λ + c . The angleˆ θ Ae is marked in figure 5. Note that generically, the heteroclinic connection γ AB will betangent to the y Ae axis at ξ A , and so generically ˆ θ Ae will be order one. In the orbit flipbifurcation which we consider in section 5.1.5, the heteroclinic connection is tangent tothe strong stable manifold, i.e., the x Ae axis, and then ˆ θ Ae will be very close to π (so | tan ˆ θ Ae | (cid:28) We consider a trajectory which starts at time t = 0, at a point x A (0) ∈ H A in , and we write the solution to the equations linearised around ξ A as x Ae ( t ) = x Ae (0)e λ + e t (15 a ) y Ae ( t ) = y Ae (0)e λ − e t (15 b ) x Ac ( t ) = x Ac (0)e λ + c t (15 c ) y Ac ( t ) = h e λ − c t (15 d )The time it takes the trajectory to travel from H A in to H A out is T , so x A ( T ) ∈ H A out , and T is defined by r Ae ( T ) = x Ae ( T ) + y Ae ( T ) = h . (16)This gives the five equations x Ae ( T ) = x Ae (0)e λ + e T (17 a ) eteroclinic bifurcations in Rock–Paper–Scissors y Ae ( T ) = y Ae (0)e λ − e T (17 b ) x Ac ( T ) = x Ac (0)e λ + c T (17 c ) y Ac ( T ) = h e λ − c T (17 d ) h = x Ae (0) e λ + e T + y Ae (0) e λ − e T (17 e )which define x Ae ( T ) , y Ae ( T ) , x Ac ( T ) , y Ac ( T ) and (implicitly) T in terms of x Ae (0) , y Ae (0) and x Ac (0), thus defining the local map from H A in to H A out . Note that we do not attempt tosolve for T at this stage. We next construct the global map from H A out to H B in . We onlyconsider trajectories which lie close to the heteroclinic connection from ξ A to ξ B , so θ Ae ( T ) will be close to ˆ θ Ae (see figure 5). Then we write θ Ae ( T ) = arctan (cid:18) y Ae ( T ) x Ae ( T ) (cid:19) (18)and Taylor expand the right hand side around ˆ x A to get θ Ae ( T ) = arctan (cid:18) ˆ y Ae + ( y Ae ( T ) − ˆ y Ae )ˆ x Ae + ( x Ae ( T ) − ˆ x Ae ) (cid:19) = ˆ θ Ae − ˆ y Ae (ˆ x Ae ) + (ˆ y Ae ) ( x e ( T ) − ˆ x Ae ) + ˆ x Ae (ˆ x Ae ) + (ˆ y Ae ) ( y e ( T ) − ˆ y Ae )= ˆ θ Ae − ˆ y Ae h x e ( T ) + ˆ x Ae h y e ( T ) (19)where we are assuming ( y Ae ( T ) − ˆ y Ae ) and ( x Ae ( T ) − ˆ x Ae ) are small and have used the factthat (ˆ x Ae ) + (ˆ y Ae ) = h .Recall that a point on H A out can be defined by the coordinates x Ac ( T ) , y Ac ( T ) and θ Ae ( T ). For a trajectory close to the heteroclinic connection, x Ac ( T ) and y Ac ( T ) are small(since the heteroclinic connection lies in P ( ξ A ) which has x Ac = y Ac = 0), and ( θ Ae ( T ) − ˆ θ Ae )is also small. A point on H B in is defined by the coordinates x Bc , x Be and y Be , which arealso all small for a trajectory close to γ AB (see equation (13)). Thus, in the global map,to first order, x Bc , x Be and y Be can be written as a linear combination of x Ac ( T ), y Ac ( T )and ( θ Ae ( T ) − ˆ θ Ae ). In addition, the global map must preserve the invariance of P ( ξ A ).The global map can thus be written to first order as: x Bc = F ( θ Ae ( T ) − ˆ θ Ae ) + F x Ac ( T ) + F y Ac ( T ) (20 a ) x Be = F x Ac ( T ) + F y Ac ( T ) (20 b ) y Be = F x Ac ( T ) + F y Ac ( T ) (20 c )where the F j are order one constants.Using equation (19), we replace θ Ae ( T ), and renaming the constants gives x Bc = A x Ac (0)e λ + c T + A h e λ − c T + A x Ae (0)e λ + e T + A y Ae (0)e λ − e T (21 a ) x Be = B x Ac (0)e λ + c T + B h e λ − c T (21 b ) y Be = C x Ac (0)e λ + c T + C h e λ − c T (21 c ) eteroclinic bifurcations in Rock–Paper–Scissors A = − F ˆ y Ae /h and A = F ˆ x Ae /h , so tan ˆ θ Ae = − A /A .Usually in these sorts of calculations, it is assumed that the order one constantswhich arise in the global map are not functions of the eigenvalues. This is not entirelytrue, as they will be dependent on the global dynamics, but to leading order, if we areonly 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 , B , C and C do have a degeneracy neara particular degeneracy of the eigenvalues, arising because of the way we have definedour coordinates.Specifically, consider the trajectory of c and w during the passage from H A out to H B in . Both c and w are assumed small, and write ( c A , w A ) for the coordinates on H A out and ( c B , w B ) for the coordinates on H B in . Then to lowest order, the global map can bewritten (cid:32) c B w B (cid:33) = (cid:32) G G G G (cid:33) (cid:32) c A w A (cid:33) (22)where the G j are indeed generically order one constants. When we rewrite this in termsof x Be , y Be , x Ac and x Bc , we have (cid:32) x Be y Be (cid:33) = (cid:32) λ − e − λ + e − (cid:33) (cid:32) G G G G (cid:33) (cid:32) λ − c − λ + c − (cid:33) − (cid:32) x Ac ( T ) y Ac ( T ) (cid:33) (23)That is, (referring to (20 b ) and (20 c )) (cid:32) B B C C (cid:33) = (cid:32) λ − e − λ + e − (cid:33) (cid:32) G G G G (cid:33) (cid:32) λ − c − λ + c − (cid:33) − (24)There are thus degeneracies in B , B , C and C when either λ − c = λ + c or λ − e = λ + e . Theformer case doesn’t occur in our system, because we assume that σ, ζ > ζ = γ : where the expanding eigenvalues change frombeing real to complex. In this case, when λ − e = λ + e , then B = C , and B = C , and thedeterminant of the matrix on the left hand side of (24) is ∆ BC = B C − C B = 0. Weassume in this section that we are away from the point where the expanding eigenvaluesare equal. In section 5.2, we consider the case where the expanding eigenvalues arecomplex, but use a coordinate change which limits to the repeating eigenvalues casewhen the imaginary part of the complex pair vanishes. Equations (21 a ) to (21 c ) map a point on H A in to a point on H B in . Due to the symmetry g in equation (4), a fixed point of a fullPoincar´e return map will also be a fixed point of (21 a ) to (21 c ). Fixed points of a fullPoincar´e return map can thus be found by dropping the A and B superscripts, and thedependence on 0 on the right hand side, to give the following four nonlinear equations,with four unknowns, x c , x e , y e and T : x c = A x c e λ + c T + A h e λ − c T + A x e e λ + e T + A y e e λ − e T (25 a ) eteroclinic bifurcations in Rock–Paper–Scissors x e = B x c e λ + c T + B h e λ − c T (25 b ) y e = C x c e λ + c T + C h e λ − c T (25 c ) h = x e e λ − e T + y e e λ + e T (25 d )We substitute equations (68) and (25 c ) into (25 a ) to eliminate x e and y e , which uponrearranging gives: x c (1 − e λ + c T ( A + A B e λ + e T + A C e λ − e T )) = h e λ − c T ( A + A B e λ + e T + A C e λ − e T ) (26)Recall that λ + c > > λ − c and λ + e > λ − e >
0. Since T is assumed to be large, and iscertainly positive, we can neglect the first (1) and second ( A ) terms on the left handside, and the first ( A ) term on the right hand side, to get: x c = − h e ( λ − c − λ + c ) T A B e λ + e T + A C e λ − e T A B e λ + e T + A C e λ − e T (27)We next substitute (27) into (68) and (25 c ) and then finally into the expression for T (25 d ), which we will then solve for T . This gives us: x e = B x c e λ + c T + B h e λ − c T = − B h e λ − c T A B e λ + e T + A C e λ − e T A B e λ + e T + A C e λ − e T + B h e λ − c T (28 a )= − hA ∆ BC (cid:32) e ( λ − e + λ − c ) T A B e λ + e T + A C e λ − e T (cid:33) (28 b )where ∆ BC = B C − C B , and y e = C x c e λ + c T + C h e λ − c T = − C h e λ − c T A B e λ + e T + A C e λ − e T A B e λ + e T + A C e λ − e T + C h e λ − c T (29 a )= hA ∆ BC (cid:32) e ( λ + e + λ − c ) T A B e λ + e T + A C e λ − e T (cid:33) (29 b )Note that when simplifying (28 a ) to get (28 b ) and (29 a ) to get (29 b ), terms in thenumerator in e ( λ + e + λ − c ) T and e ( λ − e + λ − c ) T , respectively, cancel out.We substitute (28 b ) and (29 b ) into (25 d ) to get: h = x e e λ + e T + y e e λ − e T | ∆ BC | (cid:113) A + A e ( λ − e + λ + e + λ − c ) T A B e λ + e T + A C e λ − e T A B e λ + e T + A C e λ − e T = | ∆ BC | (cid:113) A + A e ( λ − e + λ + e + λ − c ) T (30)The final task is to solve (30), which gives the period of a periodic orbit in the flow (theactual period is 3 T ), close to the heteroclinic cycle, which corresponds to a fixed pointin the map. For large T , the periodic orbit will be close to the heteroclinic cycle. Wewill do this in two different cases in sections 5.1.4 and 5.1.5. eteroclinic bifurcations in Rock–Paper–Scissors x c , x e and y e (equations (27), (28 b ) and (29 b ) respectively) so we substitute (30)into these equations to simplify them, to get x c = − h e − λ + c T A B e − λ − e T + A C e − λ + e T | ∆ BC | (cid:112) A + A (31 a ) x e = − h A sgn(∆ BC ) (cid:112) A + A e − λ + e T (31 b ) y e = h A sgn(∆ BC ) (cid:112) A + A e − λ − e T (31 c )These three equations give the coordinates of the fixed point in terms of T . Note thatin 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 beexpected.We now check that b ( t ) > t ∈ [0 , T ] for these solutions, namely that theperiodic orbit corresponds to a positive travelling wave solution of (2). Note from (15 a )and (15 b ) that x Ae ( t ) and y Ae ( t ) do not change sign, and in order to have b ( t ) > y Ae ( t ) > x Ae ( t ) for all t ∈ [0 , T ]. Writing (15 a ) and (15 b ) with initialconditions from (31 b ) and (31 c ) gives us x Ae ( t ) = − A E e − λ + e ( T − t ) (32) y Ae ( t ) = A E e − λ − e ( T − t ) (33)where E = h sgn(∆ BC ) √ A + A , and so we require A E e − λ − e ( T − t ) > − A E e − λ + e ( T − t ) (34)There are four cases to consider depending on the signs of A E and A E , andthe corresponding quadrant in x Ae - y Ae space in which the solutions lie. Since we are onlyconsidering solutions that lie close to the heteroclinic connection, we assume in eachcase that the solutions x Ae ( t ) and y Ae ( t ) lie in the same quadrant as ˆ θ Ae .If A E , A E >
0, then y Ae ( t ) > x Ae ( t ) < A E , A E <
0, then y Ae ( t ) > x Ae ( t ) < b ( t ) < t and this is not a positive travelling wave. If A E < < A E , then x Ae , y Ae < − A A = tan ˆ θ Ae < e ( λ − e − λ + e )( T − t ) (35)Note that the left-hand side of the inequality is positive and the right-hand side isbetween 0 and 1, and so we require 0 < tan ˆ θ Ae < x Ae , y Ae <
0, putting these together means that π < ˆ θ Ae < π/
4. Furthermore, solutionsmust satisfy
T < − λ + e − λ − e log(tan ˆ θ Ae ) (36) eteroclinic bifurcations in Rock–Paper–Scissors λ + e − λ − e must decrease to 0 as T → ∞ , as the heteroclinic bifurcationis approached.Finally, suppose A E < < A E , then x Ae , y Ae > − A A = tan ˆ θ Ae > e ( λ − e − λ + e )( T − t ) (37)Again, the left-hand side of the inequality is positive and the right-hand side is between0 and 1, so tan ˆ θ Ae >
1. With x Ae , y Ae >
0, this means that π/ < ˆ θ Ae < π/ b ( t ) > t if the heteroclinic connection issuch that π/ < ˆ θ Ae < π . If π < ˆ θ Ae < π/
4, then we can also find solutions with large T with b ( t ) >
0, so long as λ + e − λ − e decreases to zero as T tends to infinity. For othervalues of ˆ θ Ae , periodic solutions close to the heteroclinic cycle will not correspond topositive travelling wave solutions of the PDEs (2).In the following two sections, we consider two different cases depending on therelative size of the the two terms on the left-hand side of equation (30). λ − c + λ − e = 0 In this section, we will show that a resonant-type heteroclinic bifurcation occurs when λ − c + λ − e = 0.Suppose that A B e λ + e T (cid:29) A C e λ − e T . This will be the case if A , B , A and C are order 1, since T is large and λ + e > λ − e . Then equation (30) simplifies to1 = | ∆ BC | (cid:112) A + A A B e ( λ − c + λ − e ) T (38)or T = 1 λ − c + λ − e log( D ) (39)for D = A B | ∆ BC | √ A + A . If D <
1, then we see a branch of long-period periodic orbitsemerging from the curve λ − c + λ − e = 0 into the region where λ − c + λ − e <
0. If D > λ − c + λ − 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 − λ − c = λ − e . At this fixed point, taking the leading order term for x c in (31 a ) gives x c = − h A B | ∆ BC | (cid:112) A + A e − ( λ − e + λ + c ) T (40)This resonant bifurcation is unusual: usually you expect to see a resonantbifurcation when the contracting eigenvalue is equal to the leading expanding eigenvalue,that is, when − λ − c = λ + e [18], but here it is − λ − c = λ − e .Numerical simulations of periodic orbits close to the resonance bifurcation indicatethat x e and y e are both positive (ˆ θ Ae ≈ π/ θ Ae > π/ b ( t ) > t . Indeed, this is what we see in the numericalsimulations.We next assess whether we expect to see a ‘kink’ in the shape of the profile ofthe long-period solutions as the bifurcation point is approached. As can be seen in thetime-series plots in figure 4, a kink is observed when there is a period of time during eteroclinic bifurcations in Rock–Paper–Scissors λ + c . When the trajectory is near ξ A ,the contracting components are c and w , which are linear combinations of x Ac and y Ac ,which grow/decay exponentially at rates λ + c and λ − c respectively. Observing a kinkcorresponds to having x Ac ( t ) > y Ac ( t ) for some range of time t . Since y Ac is decaying and x Ac is growing, we will observe a kink if | y Ac ( T ) | (cid:28) | x Ac ( T ) | . We have that x Ac ( T ) = x c e λ + c T = − h A B | ∆ BC | (cid:112) A + A e − λ − e T , (41) y Ac ( T ) = h e λ − c T . (42)At the resonant bifurcation, λ − c = − λ − e and so x Ac ( T ) and y Ac ( T ) are the same order andso 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 − λ − c = λ − e ,that is, on the blue line in figure 3 with ζ > σ/
2, at the boundary between regions 1and 2. A = 0 In this section we show that a branch oflong-period periodic orbits can emerge when the heteroclinic cycle undergoes an orbitflip: that is, in the case when the heteroclinic connection is tangent to the strongunstable manifold. Recall that ˆ θ Ae gives the position at which the heteroclinic connectionintersects H A out . We have that tan(ˆ θ Ae ) = − A /A , and so as A goes to zero, ˆ θ Ae goes to π ,which corresponds to the heteroclinic connection being tangent to the strong unstablemanifold (the x Ae axis; see figure 5), that is, a point of heteroclinic orbit flip.We suppose that A is small enough that the two terms on the right hand side of (30)are of the same order, that is, neither can be discarded. We then rewrite equation (30)as A = − A C B e ( λ − e − λ + e ) T + | ∆ BC | A B e ( λ − e + λ − c ) T (43)where we have assumed A (cid:28) A and so it can be dropped from the square root. Notethat A will only be small if the expressions in both exponentials are negative, namelyif λ − c + λ − e <
0, and then as T goes to infinity, A goes to zero. This holds in regions 1and 5 of figure 3.For fixed points in this case, we find the leading order term in x c to be x c = − h e − λ + c T A B e − λ − e T + A C e − λ + e T | ∆ BC | (cid:112) A + A (44)= h e − λ + c T B (cid:16) sgn(∆ BC )e − λ + e T + B e λ − c T (cid:17) (45)Numerical simulations of periodic orbits close to the heteroclinic orbit flipbifurcation indicate that x e < < y e , ˆ θ Ae is very close to (but just less than) π , and sowe automatically satisfy the condition that b ( t ) remains positive for all time. eteroclinic bifurcations in Rock–Paper–Scissors | y Ac ( T ) | (cid:28) | x Ac ( T ) | . For solutionswhich start at the fixed point, this gives x Ac ( T ) = x c e λ + c T = hB (cid:16) sgn(∆ BC )e − λ + e T + B e λ − c T (cid:17) (46) y Ac ( T ) = h e λ − c T (47)If λ − c < − λ + e , then | y Ac ( T ) | < | x Ac ( T ) | , and we will see a kinked solution. If λ − c > − λ + e ,then | y Ac ( T ) | and | x Ac ( 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 theycontract, we would instead observe a kink in the expanding phase, that is, a change ingrowth rate from λ − e to λ + 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 thelocation of the orbit flip numerically, by solving a boundary value problem to locatethe heteroclinic orbit between ξ A and ξ B , and insisting that the heteroclinic orbit istangent to the strong unstable manifold at ξ A . The location of the orbit flip is shownby a green curve in region 5 of figure 3. This green curve coincides with the black curveof heteroclinic bifurcations. In region 5, λ − c < − λ + e , and the periodic orbits close to thisheteroclinic 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 λ − e = λ + e (the red curvein figure 3), which is to be expected, as equation (43) clearly does not generate large T solutions at this point. Equation (43) gives the possibility ofa saddle-node bifurcation between periodic orbits near the orbit-flip bifurcation. Wecompute: dA dT = − ( λ − e − λ + e ) A C B e ( λ − e − λ + e ) T + ( λ − e + λ − c ) | ∆ BC | A B e ( λ − e + λ − c ) T (48)and set dA dT = 0 to find C ( λ − e − λ + e ) | ∆ BC | ( λ − e + λ − c ) = e ( λ + e + λ − c ) T (49)giving a branch of saddle-node bifurcations of periodic orbits at A = − A C B (cid:18) C ( λ − e − λ + e ) | ∆ BC | ( λ − e + λ − c ) (cid:19) ( λ − e − λ + e )( λ + e + λ − c ) + | ∆ BC | A B (cid:18) C ( λ − e − λ + e ) | ∆ BC | ( λ − e + λ − c ) (cid:19) ( λ − e + λ − c )( λ − c + λ + e ) (50)Recall that the orbit flip bifurcations only occur if λ − c + λ − e <
0, so for the lefthand side of equation (49) to be positive, we require C >
0. The branch of saddle-nodebifurcations can terminate in the branch of orbit flip bifurcations if the right hand sideof (50) becomes equal to zero. This can happen in a number of different ways, forinstance, by the eigenvalue condition − λ − c = λ + e , or if one of the constants A or C become equal to zero. In figure 3 it appears that the first of these does not happen, andsince A and C do not depend on the eigenvalues in an obvious way, we cannot say forsure what happens at the end of the branch of saddle-node bifurcations. eteroclinic bifurcations in Rock–Paper–Scissors x Ae y Ae ˆ θ Ae H A out x A (0) x A ( T ) γ AB ˆ x A Figure 6.
The figure shows a schematic of the expanding subspace from ξ A , in thecase where the expanding eigenvalues are complex. The bold line is the heteroclinicconnection γ AB , and it intersects the Poincar´e section H A out (shown by a dotted curve)at ˆ x A . A trajectory close to the the heteroclinic connection is shown, starting at apoint x A (0) and hitting H A out at x A ( T ) (both start and end points are marked bysquares). Note that we have positive travelling wave solutions to (2) ( b >
0) when y Ae > We now repeat the Poincar´e map calculations in the region where the expandingeigenvalues are complex (regions 3 and 4). We make a different change of coordinatesnear ξ A , and instead write x Ae = λ Re b − v, y Ae = b, x Ac = λ − c c − w, y Ac = λ + c c − w. (51)In the new x Ae , y Ae coordinates, the local part of the flow becomes ddt (cid:32) x Ae y Ae (cid:33) = (cid:32) λ Re ( λ Ie ) − λ Re (cid:33) (cid:32) x Ae y Ae (cid:33) . Note that in the limit as λ Ie →
0, the Jordan form of the linear part here becomes whatone would use in the case of repeated eigenvalues. The solution to the local flow is x Ae ( t ) = e λ Re t (cid:0) x Ae (0) cos( λ Ie t ) + y Ae (0) λ Ie sin( λ Ie t ) (cid:1) y Ae ( t ) = e λ Re t (cid:18) − x Ae (0) sin( λ Ie t ) λ Ie + y Ae (0) cos( λ Ie t ) (cid:19) We note again that in the limit λ Ie →
0, these solutions are exactly those that one wouldexpect for the case with two repeated eigenvalues (in particular, the term sin( λ Ie t ) /λ Ie limits to t ).We also define new coordinates for use near ξ B : x Be = λ Re c − w, y Be = c, x Bc = λ − c a − u, y Bc = λ + c a − u. (52) eteroclinic bifurcations in Rock–Paper–Scissors r Ae and θ Ae as before, as in equation (11).We use the same Poincar´e sections, close to ξ A and ξ B : H A in = { x | y Ac = h } H A out = { x | r Ae = h } H B in = { x | y Bc = h } for some h (cid:28) ξ A is now x Ae ( t ) = e λ Re t (cid:0) x Ae (0) cos( λ Ie t ) + y Ae (0) λ Ie sin( λ Ie t ) (cid:1) y Ae ( t ) = e λ Re t (cid:18) − x Ae (0) sin( λ Ie t ) λ Ie + y Ae (0) cos( λ Ie t ) (cid:19) x Ac ( t ) = x Ac (0)e λ + c t y Ac ( t ) = y Ac (0)e λ − c t the local map then gives us x Ae ( T ) = e λ Re T (cid:0) x Ae (0) cos( λ Ie T ) + y Ae (0) λ Ie sin( λ Ie T ) (cid:1) (53) y Ae ( T ) = e λ Re T (cid:18) − x Ae (0) sin( λ Ie T ) λ Ie + y Ae (0) cos( λ Ie T ) (cid:19) (54) x Ac ( T ) = x Ac (0)e λ + c T (55) y Ac ( T ) = h e λ − c T (56)where T is again defined by x Ae ( T ) + y Ae ( T ) = h . We now note that in these coordinates, to ensure that b ( t ) > t , we willrequire that y Ae > t ∈ (0 , T ). We can write b ( t ) = e λ Re t (cid:18) − x Ae (0) sin( λ Ie t ) λ Ie + y Ae (0) cos( λ Ie t ) (cid:19) (57)= K e λ Re t sin( λ Ie t + φ ) (58)where K = x Ae (0) λ Ie + y Ae (0) , and tan φ = − λ Ie y Ae (0) x Ae (0)Thus, in order for b ( t ) to remain positive for all t ∈ [0 , T ], a clear upper bound on λ Ie T is π (because the sin changes sign with frequency π ). Thus, λ Ie < π/T , and sincewe are interested in solutions for which T is large, λ Ie will be small.The global part of the map doesn’t change, that is, we still have x Bc = F ( θ Ae ( T ) − ˆ θ Ae ) + F x Ac ( T ) + F y Ac ( T ) x Be = F x Ac ( T ) + F y Ac ( T ) (59) y Be = F x Ac ( T ) + F y Ac ( T )for some order one constants F j . eteroclinic bifurcations in Rock–Paper–Scissors F ( θ Ae ( T ) − ˆ θ Ae ) = A x e ( T ) + A y e ( T )where A and A have the same expression as in the real part, namely tan ˆ θ Ae = − A /A (although note that the values are different because the angle ˆ θ Ae is defined differentlybecause of the different coordinate transformations made).Substituting in for the right-hand side, we get F ( θ Ae ( T ) − ˆ θ Ae ) = e λ Re T (cid:0) A (cid:0) x Ae (0) cos( λ Ie T ) + y Ae (0) λ Ie sin( λ Ie T ) (cid:1) + A (cid:18) − x Ae (0) sin( λ Ie T ) λ Ie + y Ae (0) cos( λ Ie T ) (cid:19)(cid:19) (60)Putting the global map (59) together with the local map (53) to (56), using (60),renaming the constants, and finallydropping the superscripts and the dependence on 0,gives us the following equations for the fixed points: x c = A x c e λ + c T + A h e λ − c T + e λ Re T (cid:18) x e (cid:18) A cos( λ Ie T ) − A λ Ie sin( λ Ie T ) (cid:19) + y e (cid:0) A λ Ie sin( λ Ie T ) + A cos( λ Ie T ) (cid:1)(cid:1) (61) x e = B x c e λ + c T + B h e λ − c T (62) y e = C x c e λ + c T + C h e λ − c T (63) h = ( x e + y e )e λ Re T (64)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 ofthe constants B , B , C , and C . Here, we find (cid:32) B B C C (cid:33) = (cid:32) λ Re −
11 0 (cid:33) (cid:32) G G G G (cid:33) (cid:32) λ − c − λ + c − (cid:33) − (65)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, weare interested in the limit when T is large and hence λ Ie is small. To make the notationclear, we write (cid:15) = λ Ie . Recall that 0 < (cid:15)T < π , and in particular, we make the ansatz (cid:15)T = π − K(cid:15) + O ( (cid:15) )for some order-one unknown K . We demonstrate below that this ansatz is correct. Wecan then writesin (cid:15)T = K(cid:15) + O ( (cid:15) ) , cos (cid:15)T = − O ( (cid:15) )Again, we substitute the expressions (62) and (63) for x e and y e into theexpression (61) for x c . This gives x c = A x c e λ + c T + A h e λ − c T − e λ Re T ( B x c e λ + c T + B h e λ − c T )( A + A K )e λ Re T ( C x c e λ + c T + C h e λ − c T )( (cid:15) A K − A ) . (66) eteroclinic bifurcations in Rock–Paper–Scissors x c (cid:16) − A e λ + c T − e ( λ + c + λ Re ) T (cid:0) − B ( A + A K ) + C ( (cid:15) A K − A ) (cid:1)(cid:17) = h (cid:16) A e λ − c T + e ( λ + c + λ Re ) T (cid:0) − B ( A + A K ) + C ( (cid:15) A K − A ) (cid:1)(cid:17) (67)The first term in the parentheses on the left hand side of equation (67) is clearlysmaller than the others. Dropping this term, and the terms of O ( (cid:15) ) gives x c = − h e ( λ − c − λ + c ) T A − e λ Re T ( B ( A + A K ) + A C ) A − e λ Re T ( B ( A + A K ) + A C ) (68)Substituting this into the expressions (62) and (63) for x e and y e gives, after somecancellation, x e = B x c e λ + c T + B h e λ − c T = h e λ − c T ∆ AB + e λ Re T ∆ BC A A − e λ Re T ( B ( A + A K ) + A C ) (69)where ∆ AB = A B − A B . Similarly, y e = C x c e λ + c T + C h e λ − c T = h e λ − c T ∆ AC − e λ Re T ∆ BC ( A + A K ) A − e λ Re T ( B ( A + A K ) + A C ) (70)where ∆ AC = A C − A C . At this point, we note that the numerators anddenominators of the expressions in (68), (69) and (70) all contain one term which ismultiplied by e λ Re T , and one which is not. Since λ Re > T is large, we might thinkthat 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 e λ Re T consists of O (1) globalconstants which generically are non-zero. However, the expression multiplying e λ Re T inthe denominators of these fractions contains the unknown constant K . It turns out thatthis expression is very small, and in fact, in the calculations below, we approximate itto lowest order by zero when finding K . Thus, both terms in the denominators must bekept.Following this observation, we use equation (64) to compute an expression for T ,where we will ignore the terms not multiplied by e λ Re T in the numerators of both (69)and (70). First, we use (69) and (70) to compute x e + y e = h e λ − c + λ Re ) T (∆ BC ) ( A + ( A + A K ) ) (cid:0) A − e λ Re T ( B ( A + A K ) + A C ) (cid:1) Substituting into (64) and rearranging, we get B ( A + A K ) + A C = A e − λ Re T − | ∆ BC | e ( λ Re + λ − c ) T (cid:113) A + ( A + A K ) . (71)Both terms on the right-hand side are non-zero, but exponentially small (as T is large),if λ Re + λ − c <
0. If λ Re + λ − c >
0, there are no solutions to this equation as there is nothingto balance the second term on the right-hand side, which would be exponentially large. eteroclinic bifurcations in Rock–Paper–Scissors λ Re + λ − c < B ( A + A K ) + A C = 0or K = − A A − C B , confirming that K is order 1. Thus our solution for T is given by T = πλ Ie + (cid:18) A A + C B (cid:19) + O ( λ Ie ) (72)As noted above, in order for b ( t ) to remain positive for all t ∈ (0 , T ), we must have T < π/λ Ie . For large T solutions then, we require A A + C B <
0. Since A , A , B and C are functions of the global dynamics (that is, the are not solely dependent on theeigenvalues of the equilibria), we cannot say where in parameter space this conditionholds (apart from being within region 4, close to the boundary between regions 4 and 5).We note that K can change sign when A passes through zero: this occurs whenthe heteroclinic connection is tangent to the positive y e -axis. Since the coordinatechanges we have used in the real and complex cases are different, this corresponds tothe heteroclinic connection in the real case being tangent to the negative x e -axis, whichis 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 bifurcationto occur, at a location determined by the global constants. This is consistent with whatis observed in figure 3.Again, we check to see whether we expect to see kinked solutions. Recall thatto get a kinked solution, we require that | y Ac ( T ) | (cid:28) | x Ac ( T ) | . Using equation (71) inthe denominator in the x c equation (68) we can see that the denominator scales likee (2 λ Re + λ − c ) T . The numerator will be order e λ Re T . Thus x Ac ( T ) = E h e ( − λ + c − λ Re ) T e λ + c T = E h e − λ Re T for some O (1) constant E and y Ac ( T ) = h e λ − c T As noted above, λ − c < − λ Re , so | y Ac ( T ) | (cid:28) | x Ac ( T ) | , and we expect to see a kinkedsolution, as observed (see panel (b) of figure 4).In summary, we expect to find a Belyakov–Devaney-type heteroclinic bifurcationwith λ − e = λ + e , that is, on the red curve in figure 3 with ζ < σ/
2, at the boundarybetween regions 4 and 5.
In summary, we conclude that heteroclinic bifurcations can only occur on the boundarybetween regions 1 and 2 (see figure 3), or on the boundary between regions 4 and 5, orwithin regions 1 or 5. All our numerical results (detailed below) point to all periodic eteroclinic bifurcations in Rock–Paper–Scissors ζ > σ/ − λ − c = λ − e ), on the boundary betweenregions 2 and 1. If ζ < σ/
2, there are two possibilities: the heteroclinic bifurcationcan either be of Belyakov–Devaney type (expanding eigenvalues changing from real toa complex-conjugate pair), on the boundary between regions 4 and 5, or of orbit fliptype (when a constant in the global part of the map vanishes as the way in which thetrajectories between equilibria change their orientation), within region 5 or region 1.The transition between the resonance and Belyakov–Devaney-type bifurcations occursat ζ = σ/
2. The transition from Belyakov–Devaney-type to orbit flip occurs when aglobal coefficient changes sign and so cannot be deduced only from considerations ofeigenvalues.
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 heteroclinicbifurcations. We begin by showing dispersion relations computed from the ODEs (4),which relate the period of the orbit to the parameter γ (in the terminology of theODEs), or equivalently, the wavelength of the travelling wave Λ, to the wavespeed γ (in the terminology of the PDEs). The dispersion relations are computed in AUTO, byfollowing periodic orbits created the Hopf bifurcation given in (8), as done in [14]. Theresulting curves are shown in figure 7.We ran numerical simulations of the PDEs (2) for values of ζ ∈ { . , . , . , . } ,on a periodic domain of size 500. For each simulation, we started with small, randomlygenerated initial conditions, and integrated for a time period of 10 ,
000 to remove anytransient behaviour. We then sample the solution at timepoints t = 10 ,
000 + 100 k , for k = 1 , . . . ,
40. At each sample point, we compute the wavelengths and wavespeeds ofthe current solution profile. The wavelengths are computed by calculating the distances(in x ) between points which have both log( a ) = −
1, and dadx >
0. Wavespeeds at each ofthese points are computed by locally calculating dadx and dadt and using γ = dadt / dadx . Wavesare only included in the analysis if each of the three variables log a , log b and log c goesboth above and below −
1, over the wavelength of the wave. The results are plotted infigure 7(a).The first thing to notice about these results is that there is a lot of scatter. Thisis for two main reasons. Firstly, although we can compute a ‘local’ wavespeed (i.e. awavespeed for some specific point ( x, t )), we cannot reliably compute a ‘local’ value of thewavelength. Secondly, in the simulations there are many different waves travelling bothleft and right (see figure 2); whenever the waves collide there is a region of time and spacefor which the wavespeed and wavelength are not well-defined, and our computations donot take account of this. However, for each of the values of ζ shown, it can be seen thatthere is a concentration of points along the AUTO-computed dispersion relation curve.Further observations of numerical experiments indicate that for ζ = 1 .
0, solutions eteroclinic bifurcations in Rock–Paper–Scissors ζ = 0 . γ Λ (b) ζ = 1 . γ Λ(c) ζ = 2 . γ Λ (d) ζ = 3 . γ Λ Figure 7.
In each panel, the solid curve shows the wavelength (period in ξ ) Λ, as γ isvaried, of periodic orbits in the ODEs (4), computed using AUTO, with σ = 3 . ζ as indicated. Each curve of periodic orbits arises in a Hopf bifurcation onthe left (black dot), and ends in a heteroclinic (long-period) bifurcation on the right.Effectively these curves are nonlinear dispersion relations for travelling wave solutionsin the PDEs (2). We additionally show estimated wavespeeds and wavelengths fromPDE simulations as red points; the points are transparent, so darker areas indicatean accumulation of points. In (b), in addition, different coloured dots correspondto estimated wavespeeds and wavelengths in long-time behaviour for different initialconditions. Note that the scale on the y -axes is different in (b) so that these pointscan be seen. Further details can be found in the text. will often converge to a single travelling wave after sufficient time has passed (sometimesin excess of t = 50 , ζ used in these experiments, we donot observe this convergence. In figure 7(b) we show the results of further similarcomputations for ζ = 1 .
0, but now we run multiple simulations from randomly choseninitial conditions, and sample the wavelengths and wavespeeds of the solution at a singletimepoint, after the solutions have become close to a single travelling wave. Differentinitial conditions converge to travelling waves with different numbers of waves fittinginto the box, but all of these lie very close to the AUTO-computed dispersion relationcurves.The values of ζ used above, together with σ = 3 .
2, correspond to observing each of eteroclinic bifurcations in Rock–Paper–Scissors ζ = 0 . ζ = 1 . ζ = 2 . ζ = 3 .
7. Bifurcation diagrams for varying σ In this section we give some numerical results showing different bifurcation diagramsin the ( γ, ζ ) plane as σ is varied. Most of the bifurcation curves were computed usingAUTO [42]. Maintaining computational accuracy for periodic orbits close to heterocliniccycles can be difficult for two reasons. Firstly, because the periodic orbits are of verylong period, it is necessary to have a large number of mesh points defining the periodicorbit. Secondly, the heteroclinic connections lie in invariant planes where some of thecoordinates are zero. The nearby periodic orbits thus will have coordinates which arevery close to zero. In order to overcome the numerical issues associated with smallnumbers we make the following change of coordinates: A = log( a ) , U = ua , B = log( b ) , V = vb , C = log( c ) , W = wc , (73)and use differential equations for A, B, C, U, V and W in our numerical computationsinstead of the original equations. Since the periodic orbits which we are interestedin exist entirely in the positive orthant (they correspond to positive travelling wavesof (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 inthe following calculations), and say that this curve well approximates the curve ofheteroclinic bifurcations.Figure 8 shows bifurcation diagrams of system (4) for various values of σ . For easeof comparison, we rescale ζ and γ by writing ˆ ζ = ζ/σ and ˆ γ = γ/ √ σ . Note that withthis rescaling, all of the coloured lines (given by equations involving eigenvalues) in thebifurcation diagrams do not depend on σ . However, the location of the Hopf curvechanges. The grey curve shows the Hopf bifurcation curve, as given by (8), and theblack 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-nodebifurcations 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. eteroclinic bifurcations in Rock–Paper–Scissors (a) σ = 0 .
32 (b) σ = 1 . σ = 3 . σ = 10 . γ ˆ γ ˆ γ ˆ γ ˆ ζ ˆ ζ Figure 8.
Bifurcation diagrams for the ODEs (4), in (ˆ γ, ˆ ζ ) parameter space, with σ as indicated for each column. The blue line ( ζ = (cid:112) σ γ − σ ) and red curve (4 ζ = γ )are tangent at ( γ, ζ ) = ( √ σ, σ/ σ + ζ ) = 3 γ ).The purple curve ( σ + ζ = 2 γ ). The green curve is the locus of a heteroclinic orbitflip. The dark grey line is a curve of Hopf bifurcations. Periodic orbits bifurcate tothe right of this line and disappear in a curve of heteroclinic bifurcations (black). Acurve of saddle-node bifurcations of periodic orbits (light grey) exists for smaller ζ .The lower panels show zooms of the upper panels near the orbit-flip bifurcation. For all four values of σ shown, the heteroclinic curve coincides with the curve λ − c = λ − e (the light blue curve) for values of ζ greater than σ/
2. For values of ζ below σ/
2, there is a range of ζ = [ ζ ∗ , σ/
2) for which the heteroclinic curve coincides with thered curve, where the expanding eigenvalues are equal: the expanding eigenvalues arereal to the right of this curve and complex to the left of this curve. Then for ζ < ζ ∗ ,the heteroclinic curve coincides with the green curve: the curve of orbit flip heteroclinicorbits. We note that the transition point ζ ∗ is dependent on the global dynamics,and varies as σ is varied. The curve of saddle-node of periodic orbits also appears toterminates at ζ = ζ ∗ .In the lower panels of figure 8, we show zooms of each set of curves near to ˆ ζ = 0,showing the orbit flip and saddle-node curves more clearly. We note that the numericalcalculations become more difficult as σ decreases, and for this reason we do not show theorbit flip or saddle-node curves on the panel for σ = 0 .
32. In particular, we note that forthe original Rock–Paper–Scissor equations with no diffusion (1), there is a degeneracywhen σ = 0, namely that the Hopf and heteroclinic resonant bifurcations are degenerate(the Jacobian matrix at the coexistence point has imaginary eigenvalues for all valuesof ζ , and the heteroclinic orbit is at resonance for all values of ζ ). Something similarhappens in this six-dimensional system: it is simple to shown that the Hopf bifurcationcurve and the resonance bifurcation curve collapse onto one another as σ is reduced to eteroclinic bifurcations in Rock–Paper–Scissors
8. Discussion
We have summarised the results of our Poincar´e map construction in section 5.3: thePDEs (2) have been shifted into a travelling frame of reference moving at speed γ (4).In these sixth-order ODEs, travelling waves correspond to periodic solutions that arecreated in a Hopf bifurcation and, with increasing wavespeed, are destroyed in one ofthree types of heteroclinic bifurcation. Although the construction of the Poincar´e mapfollows reasonably standard lines, there are technical issues: the unstable manifoldsof the equilibria are four-dimensional, and we restated some standard definitionsof heteroclinic cycles in order to accommodate (for example) positive contractingeigenvalues. We find it advantageous to delay solving for the period T of the orbit untilthe very end, since due to cancellation of some exponential terms in the calculations, itisn’t obvious which terms can be safely neglected.The periodic orbits we find can be kinked because one of the contracting eigenvaluesis positive and the growth rate changes in magnitude but not sign at the transition fromthe contracting phase to the expanding phase. In addition, each of the three heteroclinicbifurcations is non-standard or new in some way. The resonance bifurcation, with − λ − c = λ − e , involves the leading expanding (that is, smallest positive) eigenvalue; usually itwould be the non-leading (largest positive) expanding eigenvalue, that is, − λ − c = λ + e [18].The Belyakov–Devaney-type bifurcation, with the expanding eigenvalues changing fromreal to a complex-conjugate pair, and the orbit flip heteroclinic bifurcation, where thetrajectories between equilibria change their orientation, are both new because theyinvolve a robust (codimension zero) heteroclinic cycle, rather than a higher codimensionhomoclinic orbit [39–41].It seems to be the case that stability conditions of heteroclinic cycles can be muchmore complicated than perhaps was thought several decades ago when the study ofrobust heteroclinic cycles was in its infancy. Much of this complexity perhaps arises incases where unstable manifolds have dimensions greater than one. In the case in thispaper, we have a positive contracting eigenvalue, and other types of stability are oftenseen when cycles have positive transverse eigenvalues (see, e.g. [12, 20, 31, 34]). Recently,the study of heteroclinic networks is receiving increasing attention in the literature: bydefinition, such networks must have at least one equilibrium with an unstable manifoldof dimension greater than one. The stability of heteroclinic networks is almost certainlyvery subtle [33, 43], and we expect many interesting results in this area in the future.Viewing the ODEs (4) as an Initial Value Problem, the heteroclinic cycles wedescribe are hopelessly unstable. However, the fixed points of the map correspond totravelling waves in the PDEs (2), and these may (viewed as a Boundary Value Problemon an appropriate periodic domain) be stable. We plan in future to use the results inthis paper to address the stability of the travelling waves within the PDEs: intriguing eteroclinic bifurcations in Rock–Paper–Scissors Acknowledgements
We are grateful to Graham Donovan, Cris Hasan, Mauro Mobilia and Hinke Osingafor helpful conversations, to Beverley White for reminding us to “think about theeigenvalues”, and to two anonymous referees for some very helpful reports. Weacknowledge the London Mathematical Society for financial support through a Researchin Pairs (Scheme 4) grant, and the hospitality of and financial support from theDepartment of Mathematics, University of Auckland.
References [1] R. M. May and W. J. Leonard. Nonlinear aspects of competition between 3 species.
SIAM J.Appl. Math , 29(2):243–253, 1975.[2] A. Szolnoki, M. Mobilia, L. Jiang, B. Szczesny, A. M. Rucklidge, and M. Perc. Cyclic dominancein evolutionary games: a review.
J. Roy. Soc. Interface , 11(100):20140735, 2014.[3] B. Kerr, M. A. Riley, M. W. Feldman, and B. J. M. Bohannan. Local dispersal promotesbiodiversity in a real-life game of rock-paper-scissors.
Nature , 418(6894):171–174, 2002.[4] B. Sinervo and C. M. Lively. The rock-paper-scissors game and the evolution of alternative malestrategies.
Nature , 380(6571):240–243, 1996.[5] F. H. Busse and K. E. Heikes. Convection in a rotating layer: A simple case of turbulence.
Science ,208(4440):173–175, 1980.[6] J. Guckenheimer and P. Holmes. Structurally stable heteroclinic cycles.
Math. Proc. Camb. Phil.Soc. , 103:189–192, 1988.[7] M. Kimura and T. Ohta. The average number of generations until fixation of a mutant gene in afinite population.
Genetics , 61:763–771, 1969.[8] E. Frey. Evolutionary game theory: Theoretical concepts and applications to microbialcommunities.
Physica A , 389:4265–4298, 2010.[9] B. Szczesny, M. Mobilia, and A. M. Rucklidge. Characterization of spiraling patterns in spatialrock-paper-scissors games.
Phys. Rev. E , 90(3):032704, 2014.[10] J. Hofbauer and K. Sigmund.
Evolutionary Games and Population Dynamics . CambridgeUniversity Press, Cambridge, 1998.[11] A. Scheel and P. Chossat. Bifurcation d’orbites p´eriodiques `a partir d’un cycle homoclinesym´etrique.
C. R. Acad. Sci. Paris Ser. I , 314:49–54, 1992.[12] C. M. Postlethwaite. A new mechanism for stability loss from a heteroclinic cycle.
Dyn. Syst.Int. J. , 25:305–322, 2010.[13] T. Reichenbach, M. Mobilia, and E. Frey. Self-organization of mobile populations in cycliccompetition.
J. Theor. Biol. , 254(2):368–383, 2008.[14] C. M. Postlethwaite and A. M. Rucklidge. Spirals and heteroclinic cycles in a spatially extendedRock-Paper-Scissors model of cyclic dominance.
EPL , 117(4):48006, 2017.[15] M. Krupa and I. Melbourne. Asymptotic stability of heteroclinic cycles in systems with symmetry.II.
Proc. Roy. Soc. Edin. A , 134:1177–1197, 2004. eteroclinic bifurcations in Rock–Paper–Scissors [16] M. Field. Lectures on Bifurcations, Dynamics and Symmetry , volume 356 of
Pitman ResearchNotes in Mathematics . Longman, 1996.[17] M. Krupa. Robust heteroclinic cycles.
J. Nonlin. Sci. , 7(2):129–176, 1997.[18] M. Krupa and I. Melbourne. Asymptotic stability of heteroclinic cycles in systems with symmetry.
Ergod. Theory Dyn. Syst. , 15:121–147, 1995.[19] M. Field and J. W. Swift. Stationary bifurcation to limit cycles and heteroclinic cycles.
Nonlinearity , 4(4):1001–1043, 1991.[20] P. Chossat, M. Krupa, I. Melbourne, and A. Scheel. Transverse bifurcations of homoclinic cycles.
Physica D , 100(1-2):85–100, 1997.[21] C. A. Jones and M. R. E. Proctor. Strong spatial resonance and traveling waves in Benardconvection.
Phys. Lett. A , 121(5):224–228, 1987.[22] I. Melbourne, P. Chossat, and M. Golubitsky. Heteroclinic cycles involving periodic-solutions inmode interactions with O (2) symmetry. Proc. Roy. Soc. Edin. A , 113:315–345, 1989.[23] C. M. Postlethwaite and J. H. P. Dawes. A codimension-two resonant bifurcation from aheteroclinic cycle with complex eigenvalues.
Dyn. Syst. Int. J. , 21(3):313–336, 2006.[24] V. Kirk and M. Silber. A competition between heteroclinic cycles.
Nonlinearity , 7(6):1605–1621,1994.[25] P. Ashwin and M. Field. Heteroclinic networks in coupled cell systems.
Arch. Rat. Mech. Anal. ,148(2):107–143, 1999.[26] R. Driesse and A. J. Homburg. Essentially asymptotically stable homoclinic networks.
Dyn. Syst.Int. J. , 24:459–471, 2009.[27] I. Melbourne. An example of a non-asymptotically stable attractor.
Nonlinearity , 4(3):835–844,1991.[28] C. M. Postlethwaite and J. H. P. Dawes. Resonance bifurcations from robust homoclinic cycles.
Nonlinearity , 23(3):621–642, 2010.[29] O. Podvigina. Stability and bifurcations of heteroclinic cycles of type Z.
Nonlinearity , 25(6):1887–1917, 2012.[30] O. Podvigina. Classification and stability of simple homoclinic cycles in R . Nonlinearity ,26(5):1501–1528, 2013.[31] O. Podvigina and P. Ashwin. On local attraction properties and a stability index for heteroclinicconnections.
Nonlinearity , 24(3):887–929, 2011.[32] O. Podvigina and P. Chossat. Asymptotic Stability of Pseudo-simple Heteroclinic Cycles in.
Journal of Nonlinear Science , 27(1):343–375, 2017.[33] S.B.S.D. Castro and A. Lohse. Stability in simple heteroclinic networks in R . Dynamical Systems ,29(4):451–481, 2014.[34] A. Lohse. Stability of heteroclinic cycles in transverse bifurcations.
Physica D: NonlinearPhenomena , 310:95–103, 08 2015.[35] B. E. Oldeman, B. Krauskopf, and A. R. Champneys. Death of period-doublings: locating thehomoclinic-doubling cascade.
Physica D , 146(1-4):100–120, 2000.[36] B. E. Oldeman, B. Krauskopf, and A. R. Champneys. Numerical unfoldings of codimension-threeresonant homoclinic flip bifurcations.
Nonlinearity , 14(3):597–621, 2001.[37] S.-N. Chow, B. Deng, and B. Fiedler. Homoclinic bifurcation at resonant eigenvalues.
Journal ofDynamics and Differential Equations , 2(2):177–244, 1990.[38] H. Kokubu, M. Komuro, and H. Oka. Multiple homoclinic bifurcations from orbit-flip I: Successivehomoclinic doublings.
International Journal of Bifurcation and Chaos , 06(05):833–850, 1996.[39] A. J. Homburg and B. Krauskopf. Resonant homoclinic flip bifurcations.
Journal of Dynamicsand Differential Equations , 12(4):807–850, 2000.[40] L. A. Belyakov. Bifurcation of systems with homoclinic curve of a saddle-focus with saddle quantityzero.
Mat. Zametki , 36:681–689, 1984.[41] R. L. Devaney. Homoclinic orbits in hamiltonian systems.