Numerical Study of Secondary Heteroclinic Bifurations near Non- Reversible Homoclinic Snaking Numerical Study of Secondary Heteroclinic Bifurations near Non-reversible Homoclinic Snaking

We discuss the emergence of isolas of secondary heteroclinic bifur-cations near a non-reversible homoclinic snaking curve in parameter space that is generated by a codimension-one equilibrium-to-periodic (EtoP) heteroclinic cycle. We use a numerical method based on Lin's method to compute and continue these secondary heteroclinic EtoP orbits for a well-known system.


Introduction
Heteroclinic and homoclinic bifurcations in vector fields typically feature rich dynamics in their neighbourhood, and often such bifurcations give rise to a cascade of related bifurcations.In recent years, heteroclinic connections that involve periodic orbits have attracted more and more attention, because solutions of this type play an important role in many applications [2,5,11].Accompanying the analytical results, advanced numerical techniques to find and to continue these orbits have been developed [7,4].The minimal configuration of a connecting cycle that involves a periodic orbit is an equilibrium-to-periodic (EtoP) cycle in R 3 , which already features a rich structure of bifurcations nearby.In this scenario, we generically find one EtoP connection to be of codimension 1 and the other EtoP connection to be of codimension 0. Figure 1 (a) shows a typical two-parameter bifurcation diagram of a specific system (introduced below) in the presence of such an EtoP cycle.The curve c 1 is the codimension-1 EtoP heteroclinic bifurcation, and the two branches t 0 are the tangencies of the codimension-0 EtoP connection.The snaking curve h b 1 is a codimension-1 homoclinic bifurcation to the equilibrium and a prominent example of bifurcations of an EtoP cycle.More precisely, h b 1 describes a 1-homoclinic bifurcation, that means a homoclinic solution that takes one excursion along the heteroclinic EtoP cycle before connecting back to the equilibrium.The curve h b 1 accumulates onto the part of c 1 where the complete EtoP cycle exists (between the two branches t 0 ).Some features of such an accumulation are known analytically [2,11,5]; in particular, the accumulation of the turning points is well-known (see panels (b)-(d), the orbit takes more and more revolutions around the periodic orbit before connecting to the equilibrium).
The goal of this paper is to apply a numerical technique to compute secondary heteroclinic bifurcations in the presence of such an EtoP cycle in a specific example system.By "secondary bifurcations" we mean N -heteroclinic connecting orbits of codimension 1, i.e. heteroclinic orbits that pass by close to the equilibrium N − 1 times and also N − 1 times close to the periodic orbit before connecting.We expect these secondary heteroclinic orbits near the original EtoP cycle, although they may leave its neighbourhood when they are continued in parameters.It turns out that the parameter values for which the secondary heteroclinic orbits exist are located on isolas, and these isolas are related to the non-reversible snaking curve h b 1 explained above.
This paper is organized as follows: Section 2 gives a general background to the class of systems where the described phenomena can be expected and the proposed methods can be applied.In section 3 we explain the numerical method that is used to compute the secondary N -heteroclinic bifurcations.The main part of the paper is section 4, in which we apply the numerical method to a specific system, that was originally derived in the context of modelling the dynamics of semiconductor lasers.Finally, section 5 concludes and gives some directions for future research.

General background
We consider the family of ODEs where g is sufficiently smooth and λ = (λ 1 , λ 2 ) is the family parameter.Although the methods we present here also work for higher dimensions, for simplicity we restrict our considerations to systems in R 3 .We assume that the following conditions hold for λ = (0, 0): (C1) There is a hyperbolic equilibrium E and a hyperbolic periodic orbit P , where dim W u (E) = dim W u (P ) = 2, (C2) there is a robust heteroclinic orbit q 0 ⊂ W u (E) ∩ W s (P ), (C3) there is a non-robust heteroclinic orbit q 1 ⊂ W s (E) ∩ W u (P ), and for λ = (λ 1 , 0), |λ 1 | sufficiently small, we assume that (C4) W s (E) and W u (P ) split up with non-zero velocity along each curve intersecting the λ 2 -axis transversally.In other words, we assume that the system has a heteroclinic equilibrium-toperiodic-orbit cycle (EtoP cycle) for critical parameter values λ = (0, λ 2 ), |λ 2 | sufficiently small.The EtoP cycle consists of a robust connection q 0 and a connection q 1 that lies in the intersection of two manifolds that split up when λ 1 is moved away from the critical value, hence λ 1 is considered the parameter that measures the splitting of the manifolds W s (E) and W u (P ).The parameter λ 2 is considered the parameter that controls the intersection of the two two-dimensional manifolds W u (E) and W s (P ).In the following we assume that the manifolds W u (E) and W s (P ) intersect either transversally or have a quadratic tangency (and the unfolding is controlled by λ 2 ).
From theory it is known that in this scenario a homoclinic orbit to E bifurcates from the EtoP cycle [5].This homoclinic orbit stays close to the EtoP cycle in phase space and takes more and more turns around P before returning to E for parameter values approaching the critical value.If in addition to (C1)-(C4) certain conditions on the intersection of W u (E) and W s (P ) are satisfied (similar to the conditions discussed in [1]), then it can be shown that the bifurcation curve of this codimension-1 homoclinic orbit to E is indeed a snaking curve in the parameter space (λ 1 , λ 2 ), and we refer to this phenomenen as "non-reversible homoclinic snaking"; see also figure 1 (a).Note that the term "homoclinic snaking" is typically used for snaking phenomena in time-reversible systems, which is not the case in our scenario.We remark that the exponential accumulation rate of the snaking curve is given by the stable and unstable Floquet exponents of P ; see again figure 1 (a) and the explanations in sections 4 and 4.1 below.The complete discussion and analytical considerations of this phenomenon is work-in-progress and beyond the scope of this text; they will be reported elsewhere.

Computational method
We use the numerical method introduced in [7] to find and continue the secondary codimension-1 heteroclinic EtoP solutions; see [7] and references therein for a complete description of the method.
Let Σ denote a cross-section dividing the phase-space such that E and P are not in the same half-space.The main idea is to construct two orbit segments q − 1 in W u (P ) and q + 1 in W s (E), such that q − 1 (0) ∈ Σ and q + 1 (0) ∈ Σ, while both q + 1 and q − 1 are already close to the sought-after solution q 1 .According to analytical results from Lin's method [8,9,5,11], in a neighbourhood of q 1 there is one unique pair of orbit segments such that the start-respectively endpoints are in Σ and their difference is in a pre-defined given linear subspace Z of dimension 1.This is actually true for every 1-dimensional linear subspace Z that is not contained in the tangent spaces of the involved stable and unstable manifolds.In other words, for given Z the two orbit segments q − 1 and q + 1 are uniquely defined and the gap size |q − 1 (0) − q + 1 (0)| is a well-defined test function that can be used to detect real orbits by numerical continuation with respect to a family parameter and stopping when the gap size is zero.
In this paper we use exactly this approach: We construct two orbit segments q − 1 and q + 1 that both intersect Σ and have a jump in a certain direction Z.However, our goal is to detect and continue N -heteroclinic orbits, i.e. orbits that pass by N − 1 times close to E before connecting to E; hence, the constructed orbit segments q − 1 and q + 1 need to intersect Σ multiple times in a neighbourhood of Σ ∩ q 1 .More precisely, combined they both need to intersect Σ a total of N + 1 times (including the start-/endpoints) in a neighbourhood of q 1 ∩ Σ (either of them may lead from Σ directly to E/P ).In practice, we construct such orbits by numerical continuation with respect to time while keeping the startpoint of the orbit segment in the respective stable/unstable manifolds close to E or P .Note that in [10] a similar approach is used to compute N -homoclinic orbits to E, but there multiple orbit parts are constructed and the resulting gaps are closed one-by-one in order to find N -homoclinic orbits.Using only two orbit segments and consequently only one gap has a computational advantage regarding the complexity of the continuation problem, but the disadvantages that it is more difficult to construct initial orbit segments q − 1 and q + 1 , and that there is no guarantee that q − 1 and q + 1 are indeed close to the original heteroclinic EtoP cycle.
The actual computations presented in section 4 were performed with the continuation software AUTO [3].Although the described approach is straight forward, the computations are very sensitive and for N > 2 increasingly difficult.

Demonstration
In this section we use the method introduced in section 3 to find and continue N -heteroclinic orbits in a dynamical system given by the following ODE: (2) The given vector field describes the dynamics near a saddle-node Hopf bifurcation in the presence of a global reinjection mechanism and was introduced in [6].The parameters ν 1 and ν 2 are the unfolding parameters of the saddle-node Hopf bifurcation, while the parameters ω, α, β, s, c, d and f are fixed at the values This choice corresponds to the unfolding of type A that has been studied in [6].
Note that the variable ϕ is 2π-periodic and global reinjection is realized by trajectories that connect a neighbourhood of a saddle-node Hopf point with one of its symmetric copies.Hence, a global reinjection corresponds to a large excursion near the circle S 1 = {x = y = 0}.Note that this circle is not invariant because d = 0 and f = 0 (where rational ratios are avoided).
In order to plot some of the computed orbits in phase-space, we use the following transformation, which maps orbits that are close to S 1 into a torus in (u, v, w)-space: To avoid self-intersections, we fix R = 2.
As shown numerically in [7] and [5], system (2) satisfies conditions (C1)-(C4) for a certain parameter range (implicitely assuming an appropriate parameter transformation such that (ν 1 , ν 2 ) satisfy the conditions).Hence, there is a codimension-1 heteroclinic EtoP connection q 1 which exists on a curve c 1 in (ν 1 , ν 2 ) parameter space, and a codimension-0 heteroclinic EtoP connection q 0 which is robust, its tangencies are the branches of the bifurcation curve t 0 .Moreover, from this EtoP cycle bifurcates a homoclinic orbit to E, and the corresponding bifurcation curve h b 1 shows "non-reversible homoclinic snaking", i.e. h b 1 snakes and accumulates onto the curve segment of c 1 where the full EtoP cycle exists; see figure 1 (a).The turning points of the snaking curve h b 1 seem to lie on the two branches of the tangency curve t 0 .However, note that they are not exactly on t 0 , but exponentially close: they are accumulating onto t 0 with an exponential rate given by the unstable Floquet exponent of P , and they are accumulating onto c 1 with an exponential rate given by the stable Floquet exponent of P .
During the accumulation process, the homoclinic orbit stays longer and longer near P , meaning that it performs more and more revolutions along P before returning to E; cf. the time-vs-w plots in panels (b)-(d) of figure 1. 4.1.Rescaling.In this section we explain the rescaling that is applied to the bifurcation diagrams in panels (b) of figures 2, 4 and 6.According to the analytical results [5], the accumulation rate of h b 1 is related to the stable/unstable Floquet multipliers of the periodic orbit.More precisely, the overall exponential accumulation rate of the snaking curve to the curve segment of c 1 on which the EtoP cycle exists is given by the stable Floquet exponent of P , while the exponential rate of the displacement of the turning points of the snaking curve relative to t 0 is given by the unstable Floquet exponent of P .
In order to stretch and turn the snaking curve, we use the following nonlinear transformation: We approximate the curve c 1 by a parabola, whose vertex is then moved in parallel to the upper branch of t 0 .Likewise, we approximate the upper branch of t 0 by a parabola whose vertex is moved parallel to c 1 .(We choose parabolas to account for the fact that c 1 and t 0 indeed are not straight lines.)Each point on h b 1 can then be expressed as an intersection point of two such parabolas, and its coordinates are simply the arclength of the respective parabola parts, meaning the arclength of the segment from the intersection point to t 0 (c 1 , respectively).In a final step, the transformed ν 1 -coordinate is logarithmized to equalize the exponentially fast accumulation and then scaled by the stable Floquet multiplier.The resulting curve h b 1 in the new coordinates (ν 1 , ν2 ) is sinusoidal with approximately constant amplitude and period, giving a much clearer picture of the snaking phenomenon.
Note that the rescaling is only applied to the bifurcation diagram, all computations are performed in the original parameter space.4.2.Two-heteroclinic isolas.Using the method described in section 3, it is possible to find and continue 2-heteroclininc EtoP connections of codimension one near the homoclinic snaking curve h b 1 .Figure 2 shows the bifurcation diagram, panel (a) in the original parameter space (ν 1 , ν 2 ) and panel (b) in the rescaled parameter space (ν 1 , ν2 ) (see section 4.1).
We observe that the 2-heteroclinic EtoP connections are on closed curves in parameter space, we denote these 2-isolas by h i 2 , i = 1, 2, 3.They all lie in vales of the homoclinic snaking curve h b 1 .We observe that the shapes of the isolas h i 2 change: the left fold point enters the neighbouring vale of h b 1 and reaches further down for increasing i.Further we observe that the isolas do not reach all the way down the vales, and, what is more, they rise beyond the upper branch of t 0 , which suggests that it is not the tangency curve t 0 that bounds the isolas from below and from above.At present it is not known how the isolas are bounded.Some representative orbits on the 2-isolas are plotted in figure 3. The parameter values of all these orbits are set to the leftmost lower turning points of the isolas.We notice that the orbit on h 1 2 , shown in panel (a1) as phase-space and in panel (a2) as time-vs-w plots, starts off close to P , passes by E rather quickly, and then takes three turns close to P before connecting to E. The orbit on h 2 2 (panel (b)) behaves almost the same, but it takes four turns close to P before connecting to E, and the orbit on h 3 2 (panel (c)) takes five turns close to P .This observation suggests that the isolas h i 2 are defined by the number of revolutions the 2-heteroclinic orbit takes before connecting to E, which increases by one for each consecutive 2-isola.We remark that the number of revolutions also varies along the same isola, the orbits on h 1 2 can have three or four revolutions (not shown in the figure), this value seems to depend on the vale of h b 1 that the orbit is located in.4.3.Three-and four-heteroclinic isolas.In the same manner as for the 2heteroclinic EtoP connections of codimension 1, it is possible to compute 3-and 4-heteroclinic EtoP connecing orbits.Figure 4 shows the bifurcation diagram for 3-heteroclinic EtoP connections, panel (a) in the original parameter space and panel  1 and appear to be bounded from below and from above similar to h 1 2 , h 2 2 , h 3 2 .Even in the rescaled figure it is hard to see that the isolas are indeed closed curves; see also figure 8 below for a sketch of the isolas.
The corresponding 3-heteroclinic orbits are shown in figure 5.After starting close to P , they pass by E rather quickly, take one additional turn around P , then pass by E for the second time, then take four, five or six additional turns near P before connecting to E. The number of revolutions around P prior to connecting to E is again the factor that defines the different isolas h i 3 ; this number increases by one for increasing i.Note that this number is constant along one isola.
Figure 6 shows the only 4-isola h 2 4 we were able to compute; constructing the orbit segments to initialize the method explained in section 3 becomes increasingly difficult.Although it is not quite clear from the figure, the 4-isola h 2 4 is indeed a closed curve.The sketch in figure 8 explains the relative position of the isola h 2 4 compared to the isolas h 2 2 and h 2 3 , note that h 2 4 borders h 2 3 closely.Figure 7 shows a representative orbit on h 2 4 in phase space (panel (a)) and as a time-vs-w plot (panel (b)).We observe that the orbit starts close to P , passes by E for the first time, then takes four revolutions around P , passes by E for the second time, takes one additional loop close to P , passes by E for the third time, finally takes five revolutions around P and connecting to E. Note that the number of revolutions that the orbit takes before connecting to E is the same as for the orbit on h 2 3 , so it seems again to be this number that defines the different isolas.Also note that the first two numbers of revolutions around P indeed vary for different orbits along h 2 4 (not shown), but they always sum up to five.

Conclusion
We showed that the numerical method introduced in [7] can be used to find and continue secondary heteroclinic bifurcations of EtoP cycles.In an example system, we computed the 2-, 3 and 4-heteroclinic bifurcations that lie on isolas in parameter space and find that these isolas are closely related to a snaking curve of a homoclinic bifurcation.
An analytical verification of this observation as well as finding the mechanism and bifurcation curves that bound the isolas from below and from above is left for future research.

Figure 1 .
Figure 1.Non-reversible homoclinic snaking.The bifurcation diagram (a) shows the bifurcation curve of the codimension-one EtoP connection c 1 , the tangencies of the codimension-zero EtoP connection t 0 and the homoclinic bifurcation h b 1 , which accumulates onto the curve segment of c 1 where the complete EtoP cycle exists.Selected homoclinic orbits to E along h b 1 are shown as integration time-vs-w plots in panels (b)-(d).T denotes the total integration time.

Figure 7 .
Figure 7. Representative orbit on the 4-isola.Panel (a1) shows a representative 4-heteroclinic orbit on h 2 4 in phase space, panel (a2) shows the same orbit as a time-vs-w plot.

Figure 8 .
Figure 8. Sketch of the heteroclinic isolas.All of the bifurcation curves are closed curves, panel (a) shows a sketch of the 2-isola h 2 2 , panel (b) shows h 2 2 in grey and the 3-isola h 2 3 in black, and panel (c) shows h 2 2 and h 2 3 in grey and the 4-isola h 2 4 in black.The sketch is not accurate in scale, but gives a good impression of the relative positions of the isolas.