Les Canards de Turing
Theodore Vo Arjen Doelman Tasso J. Kaper

abbrv
Heisenberg-limit spin squeezing with spin Bogoliubov Hamiltonian
    Wenxian Zhang
    September 9, 2024
================================================================




§ ABSTRACT 
In this article, we study a prototypical system of reaction-diffusion equations in which the diffusivities are widely separated. We report on the discovery of families of spatially periodic canard solutions that emerge from singular Turing bifurcations.  
We show that the small-amplitude, spatially periodic solutions that emerge from the Turing bifurcations form families of spatially periodic canards that oscillate about the homogeneous equilibrium, with wavenumbers near the critical value obtained from the Turing analysis. 
The emergence of these spatially periodic canards asymptotically close to the Turing bifurcations, which are reversible 1:1 resonant Hopf bifurcations in the spatial ODE system, is an analog in spatial dynamics of the emergence of limit cycle canards in the canard explosions that occur asymptotically close to Hopf bifurcations in time-dependent ODEs. 
We also find families of large-amplitude, spatially periodic canards, including some with 𝒪(1) wavenumber and some with small wavenumbers. 
These lie further from the homogeneous state and have a “fast-slow" spatial structure, with segments of steep gradients and segments of gradual variation.

In the full PDE system, we show that for most parameter values under study the Turing bifurcation is sub-critical, and we present the results of some direct numerical simulations showing that several of the different types of spatial canard patterns are attractors in the prototypical PDE.
    
To support the main numerical discoveries, we use the method of geometric desingularization and geometric singular perturbation theory on the spatial ODE system to demonstrate the existence of these families of spatially periodic canards. 
Crucially, in the singular limit, 
we study a novel class of reversible folded singularities of the spatial ODE system. 
In particular, there are two reversible folded saddle-node bifurcations of type II (RFSN-II), each occurring asymptotically close to a Turing bifurcation. We derive analytical formulas for these singularities and show that their canards play key roles in the observed families of small-amplitude and large-amplitude spatially periodic canard solutions. 
Then, for an interval of values of the bifurcation parameter further below the Turing bifurcation and RFSN-II point, the spatial ODE system also has spatially periodic canard patterns, however these are created by a reversible folded saddle (instead of the RFSN-II).   
It also turns out that there is an interesting scale invariance, so that some components of some spatial canards exhibit nearly self-similar dynamics. 



Key words. 
folded singularities,
spatial canards,
singular Turing bifurcation,
periodic solutions,
Turing instability,
spatial dynamics,
reversible systems,
nearly self-similar dynamics,
subcritical Ginzburg-Landau

MSC codes. 
35B36, 34E17, 34E15, 35B25




§ INTRODUCTION
		


In mathematical models of pattern formation in biology, chemistry, ecology, engineering, material science, physics, and many other fields, the Turing bifurcation <cit.> is one of the key mechanisms that generates spatially periodic patterns. 
It was a pioneering discovery of Alan Turing in 1952 that diffusion can destabilize spatially homogeneous steady states, which are stable states of the associated reaction kinetics, and that this instability to plane wave perturbations results in the formation of spatially periodic patterns as the attractors. See for example <cit.>.

In this article, we report on the discovery and analysis of reaction-diffusion systems in which the spatially periodic solutions created in Turing bifurcations are spatially periodic canards. 
Like the known spatially periodic patterns that emerge from Turing bifurcations, 
these spatially periodic canards oscillate about the homogeneous equilibrium state. 
However, unlike the known periodic patterns, they consist of canard segments generated by folded singularities.

In the systems of spatial ordinary differential equations (ODEs) that govern the time-independent states of reaction-diffusion models,  Turing bifurcations correspond to reversible 1:1 resonant Hopf bifurcation points, see for example  <cit.>.
We show that, asymptotically close to these reversible 1:1 Hopf bifurcations, the spatial ODE systems can have reversible folded saddle-node singularities of type II (RFSN-II) and that, together with the true and faux canards attached to them, these singularities can serve as the mechanisms responsible for the creation of the observed spatially periodic canard patterns.
We focus first on the van der Pol partial differential equation (PDE),


    u_t    = v - f(u) + d u_xx, 
     
        v_t    =  (a-u) + v_xx,


and later generalise to a class of activator-inhibitor systems.
The van der Pol PDE is a prototypical reaction diffusion system of activator-inhibitor type.
Here, t ≥ 0, x ∈ℝ, u represents a voltage (or more genrally an activator), v represents a recovery variable (or more generally an inhibitor), the nonlinear reaction function is f(u) = 13u^3-u, a is a threshold parameter,  measures the timescale separation for the underlying kinetics, and 0< d ≪ 1 is the ratio of the diffusivities. 
With this small parameter, the activator diffuses more slowly than the inhibitor.
This spatial scale separation, i.e., the difference between the diffusivities of the interacting species, is an important –though not necessary– feature of the emergence of periodic patterns in spatially extended systems.

For the van der Pol PDE (<ref>), we consider all positive 𝒪(1) values of the  parameter . 
This will allow us to consider the general kinetics of the van der Pol system. 
We recall that, in the regime of small , the van der Pol ordinary differential equation (ODE), which corresponds to the ODE satisfied by spatially-homogeneous solutions of (<ref>), is in the strongly nonlinear limit, and the limit cycles of the temporal ODE are relaxation oscillations, created in an explosion of temporal limit cycle canards. 
By contrast, in the regime of large , the associated oscillations are in the weakly nonlinear limit, and the limit cycles of the kinetics ODE are small perturbations of circular orbits. 
(We refer the reader to <cit.> for analysis of the classical van der Pol ODE, i.e. of the ODE satisfied by the spatially-independent solutions of (<ref>).)

The main outcomes of this article are the numerical and analytical demonstrations that the van der Pol  PDE  (<ref>) possesses several main types of spatially periodic canards.
We find small-amplitude canards with 𝒪(1) wavenumbers (and hence 𝒪(1) spatial periods), and small-amplitude canards with small wavenumbers (and hence large spatial periods). 
These small-amplitude solutions emerge along branches emanating from the Turing point, and they oscillate spatially near the homogeneous state (see Fig. <ref>). 
In addition, we find large-amplitude canards with 𝒪(1) wavenumbers, as well as large-amplitude canards with small wavenumbers (and hence large periods, which are also referred to as “near-homoclinic" periodic solutions). 
These lie further from the homogeneous state in norm and have a distinct “fast-slow" spatial structure, with segments of steep gradients and segments of gradual variation. 
The small-amplitude canards transition continuously into those with large-amplitudes. 
Moreover, at each value of a in the main interval studied, one can find spatial canards of different types.



We use the method of geometric desingularization (also known as the blow up method) and analysis of folded singularities,  especially of the new RFSN-II singularities, to establish the existence of the different types of spatially periodic canards.
The time-independent solutions of (<ref>) satisfy the following fourth-order system of ODEs, in which the spatial variable x is the independent variable: 


    δ u_x    = p, 
    δ p_x    = f(u) - v, 
     
        v_x    = q, 
    
        q_x    = (u-a),


where δ = √(d).

We show that this spatial ODE system has a two-dimensional critical manifold in the four-dimensional (u,p,v,q) phase space.
This manifold is induced by the cubic function v=f(u), and hence it has three branches. 
We will see that two branches consist of saddle equilibria and one of center equilibria.
Crucially, there are fold sets that separate these  branches, and the RFSN-II singularities that are responsible for creating the spatial canards for values of a near the Turing value 

    a_T = √(1 -2δ√())

lie on these fold sets.

Additionally, we study the desingularized reduced vector field on the critical manifold.
We will show that the desingularized reduced system has a cusp point precisely at a=1, where the RFSN-II singularity occurs in the limit δ=0. 
Then, for all sufficiently small values of δ>0, we show that key center-unstable and center-stable manifolds coincide
at locally unique critical values of a.
The first of these is 

    a_c (δ) = 1 - 5/48δ^2 +𝒪(δ^3).
 
It corresponds to the locally unique parameter value at which key center-unstable and center-stable manifolds coincide.
Of all of the canards that pass through the neighborhood of the cusp point, the maximal canard has the longest segments near the stable and unstable manifolds of the cusp point (Fig. <ref>(a)). 
There are also maximal canards that make a small loop about the equilibrium. as seen in the projection on to the (u,q) coordinates  in phase space (Fig. <ref>(b)). 
Then, under further variations in a, additional small loops develop around the equilibrium state, and the solutions exhibit nearly self-similar dynamics (Fig. <ref>(c)). 
The analysis near the RFSN-II singularity will be valid for all values of a and δ such that a = 1 + 𝒪(δ^3/2). 

 


In the spatial ODE, we also study the spatially periodic solutions for values of a further from the Turing value.
We will show that there are reversible folded saddles (RFS) on the fold sets for all 𝒪(1) values of a ∈ [0,1). 
Together with their true and faux canards, these RFS points will be shown numerically to be responsible for creating the spatial canards in this interval away from the Turing bifurcation a_T. 

To complement the analysis of the periodic canard solutions of the spatial ODE system, we also present the results of some direct numerical simulations of the PDE (<ref>). 
These reveal that several of the large-amplitude spatial canards are stationary attractors. 
Moreover, we will find that there are also parameter regimes in which the attractors are small-amplitude solutions that are periodic in both space and time, with profiles given by spatial canards of the form discovered here. 
Overall, it turns out that the main Turing bifurcations we study (in the parameter regime in which δ^2 < 64/625) are sub-critical bifurcations (also referred to as the focusing case), in which the coefficient on the cubic term in the Ginzburg-Landau equation is positive.
We will discuss how the spatial canard patterns here constitute a new class of nonlinear attractors in the sub-critical case, and how the canards provide a new selection mechanism.

Having outlined the main ODE and PDE results in this article, we briefly highlight the motivation for this study, and we provide a short comparison of the new spatial canards to classical temporal canards in fast-slow systems of ODEs.
To highlight the motivation, we recall that Turing bifurcations in two-component reaction-diffusion systems with widely-separated spatial scales, such as (<ref>), are reversible, 1:1 resonant Hopf bifurcations in the associated spatial ODE systems 
<cit.>.
In particular, the critical points that undergo Turing/1:1 Hopf bifurcations cannot lie on normally hyperbolic slow manifolds, but must instead be located on the boundaries of such manifolds. 
This is precisely the case with the RFSN-II singularities that lie on the fold sets between the two-dimensional, saddle and center slow manifolds here in the spatial ODE. 
Therefore, in some respect, from the point of view of spatial dynamics, there should be a natural connection between Turing bifurcations and the formation of canards. 

As to the comparison, the spatial canards introduced in this article may be viewed as analogs in spatial dynamics of the temporal limit cycle canards first discovered in the van der Pol ODE <cit.>, as well as of the canards of folded saddles and folded saddle-nodes in many other (temporal) systems, see for example <cit.>. 
As just recalled, the parameter values at which Turing bifurcations occur are reversible 1:1 resonant Hopf bifurcations in the spatial ODE system <cit.>. 
These are the analogs in spatial dynamics of the (singular) Hopf bifurcations that occur in fast-slow ODEs.
Then, just after the Turing bifurcations, families of spatially periodic canards are created, with small-amplitudes for a close to a_T,  wavenumbers close to the critical wavenumber k_T (determined by the point of marginal stability), and spatial profiles close to the plane wave e^ik_T x. 
Hence, these are analogs in spatial dynamics of the small-amplitude, temporally oscillating solutions that exist close to the Hopf bifurcation in fast-slow systems of ODEs and that oscillate essentially as e^iω t, where ω is the imaginary part of the eigenvalues at the Hopf bifurcation. 
Moreover, for parameter values further from a_T,  the spatial canards appear to be created by folded saddles, and they appear numerically to exist over a broad range of parameter values, as is also observed for folded saddles and their temporal canards in fast-slow systems of ODEs 
(see for example <cit.>).
Also, in the spatial dynamics, the maximal spatial canards act as separatrices that locally partition the phase and paramater spaces into regions of distinct spatial behavior, which is also analagous to the roles played by maximal canards in the phase spaces and parameter spaces of temporal fast-slow ODEs. Finally, in comparing, we will see that there are also key new features of the spatial canards that arise due to the dimension and geometry of the critical manifolds.

Zooming out more broadly, we suggest that the results presented here for spatially periodic canards also contribute to the growing literature about canards and bifurcation delay in spatially-extended systems.
Canards in the kinetics of a two-component model for the Belousov-Zhabotinsky reaction were demonstrated to play a crucial role in the nucleation and annihilation of trigger waves in a 1-D medium of phase waves <cit.>. 
Canards arise in the traveling wave ODEs of some reaction-diffusion systems <cit.>. 
Slow passage through a saddle-node bifurcation in linear and semi-linear heat equations in 1-D can lead robustly to solutions that spend long times near unstable states as shown analytically in <cit.>. 
In astrophysics, folded saddles and their canards play a central role in a model of solar wind  when there is a steady, spherically-symmetric outflow from the surface of a star <cit.>. 
In an Amari-type neural field integral model <cit.>, temporal canards were observed in the spatial patterns of coherent structures, as were  some more complex spatio-temporal patterns containing canard segments.

More recently for PDEs, canard solutions have been studied in an ODE model derived from a sub-critical, infinite-dimensional, pattern-forming system with nonlinear advection on a bounded domain
<cit.>, where they play a role in the nonlinear transitions between two primary states of the system describing the locations of stationary fronts.
Spatio-temporal canards serve as boundaries in multi-mode attractors of reaction-diffusion systems, in which different regions of the domains exhibit different modes of stable oscillation and the canards mediate the transition intervals, keeping  the regions separated, see 
<cit.>.
Delayed loss of stability occurs in nonlinear PDEs that undergo slow passage through Hopf bifurcations, with examples including the CGL equation, the Brusselator model, the FitzHugh-Nagumo PDE, and the Hodgkin-Huxley PDE (see <cit.>).
The solutions remain for long times near unstable states in a rich manner governed by space-time buffer curves. 
A rigorous framework for the local analysis of canard solutions and other forms of bifurcation delay was developed in <cit.> for systems in which the fast variables are governed by a PDE (i.e., infinite-dimensional dynamical system) and the slow variables are governed by ODEs (i.e., by a finite dimensional dynamical system). 
Slow passage through fold bifurcations has been studied in fast-slow systems of reaction-diffusion equations, using a Galerkin approach <cit.>. 
Slow passage through Turing bifurcations has been studied in reaction-diffusion systems (see  <cit.>), as has slow passage through pitchfork bifurcations in Allen-Cahn type equations, in the presence of quenching fronts with small spatial gradients (see <cit.>).




The article is organized as follows.
Section <ref> contains the application of classical analysis to identify the Turing bifurcations and Turing-Hopf bifurcations in (<ref>) and the application of the classical normal form analysis to show that the Turing bifurcation  to spatially periodic canards is sub-critical in the main parameter regimes we study, and also to identify where it is super-critical. 
In Section <ref>, we introduce the four main types of spatial canards, and present the results obtained from numerical continuation to identify regimes in the (a,k) parameter plane in which spatial canards exist. 
In Section <ref>, we begin the geometric singular perturbation analysis by analyzing the fast system, also known as the layer problem. We identify the jump conditions for the fast homoclinics that generate spikes and for the sharp-interfaces (Proposition <ref>). Also, we identify the cusp of the fast system singularity that will be the organizing center for the spatial canard dynamics. 
We continue the geometric singular perturbation analysis in Section <ref> by deriving the desingularized reduced vector field on the critical manifold and by studying the slow flow.
We show that there are folded saddle singularities with reversibility symmetry. Moreover, these undergo reversible folded saddle-node bifurcations of type II (RFSN-II) under variation of the system parameters. 
Then, in Section <ref>, we rigorously analyze the dynamics around the RFSN-II using the blow-up technique. We determine the key parameter values for which there is an explosion of spatial canards with reversibility symmetry. 
Next in Section <ref>, we analyze in detail the geometry of the small-amplitude and large-amplitude spatially periodic canard solutions where the analysis is informed by the rigorous results about the fast system, the desingularized reduced vector field, and the folded singularity with its canards 
in Sections <ref>–<ref>.  
In Section <ref>, we use that same information to deconstruct the bifurcation sequences along isolas of spatially periodic canards.
Then, in Section <ref>, we analyze the aspects of the spatial canard dynamics that are nearly self-similar. 
Section <ref> describes how the spatial canards are analogs in spatial dynamics of the temporal limit cycle canards found in many fast-slow ODEs, such as the van der Pol ODE, FitzHugh-Nagumo ODE, the Lengyel-Epstein model, the Kaldor model, among others, and we identify important differences.
The final new results are in Section <ref>, where we present results from some direct PDE simulations that complement the analysis. Several of the different types of spatial canards discovered here are    observed to be attractors in the PDE (<ref>).
We conclude the article in Section <ref> with a summary of the main results, as well as a generalization from the prototype (<ref>) to a class of activator-inhibitor systems, and discussion of future work and open problems.
The appendices contain further information about the main numerical methods we employed, as well as the proofs of some of the propositions and lemmas. 





§ TURING BIFURCATION TO SPATIALLY PERIODIC SOLUTIONS
		

In this section, we apply the classical Turing analysis <cit.> (see also <cit.>) to derive the neutral stability curve of the homogeneous steady state (u,v)=(a,f(a)) of (<ref>).
There are Turing and Hopf bifurcations that occur asymptotically close to each other in the parameter a. 
We will obtain the critical wavenumbers and system parameters for the onset of spatially periodic patterns. 
In addition, we will apply the classical normal form theory <cit.> for reversible 1:1 resonant Hopf points to show that, in the main parameter regimes studied here with 0<δ≪ 1,  the bifurcation to spatially periodic solutions is sub-critical; and, we also identify the conditions under which  it is super-critical.

We linearize (<ref>) about the homogeneous state  and then Fourier transform the system in space, so that the governing equations become


    [ U̇; V̇ ] = DF [ U; V ] = [ -f^'(a)-dk^2            1;            -         -k^2 ][ U; V ].


The overdot denotes the time derivative, k is the wavenumber, and U(k,t) and V(k,t) denote the Fourier transforms of u(x,t) and v(x,t). 
For each k ∈ℝ, the Jacobian has two eigenvalues λ_±(k), and the state (a,f(a)) is spectrally stable as a solution of the PDE (<ref>) if Re(λ_±(k)) < 0 for all k∈ℝ.  
The trace of the Jacobian, tr DF= -f^'(a)-(d+1)k^2, is negative whenever 
f'(a) = a^2 - 1 > 0, and positive for -1<a<1 (where f'(a)<0). 
Hence, we recover the classical (non-spatial) Hopf bifurcation of the van der Pol equation at a=1. 
See the left panel in Fig. <ref>.



Now, for a<1, there are two key dynamical features. First, a narrow interval of unstable wavenumbers exists for a<1, centered around k=0, so that the homogeneous state becomes unstable to long wavelength perturbations that oscillate in time as e^i ω_H t, with ω_H =Im(λ_+(0)).  However, this instability is only a weak instability, since the rate constant (Re(λ_+(0))=√( d)) is asymptotically small for 0<d ≪ 1 and =𝒪(1).

Second, a Turing instability occurs slightly below a=1.
It appears at the critical wavenumber and critical parameter determined by the conditions


    DF = 0    and   ∂/∂ k^2 DF = 0


(for which λ_± (k) ∈ℝ).
For the van der Pol PDE (<ref>), we find 


    k_T^2 = - f^'(a_T)/2d = √(/d)   and    a_T = ±√(1 - 2√( d)).


The focus in this article will be on the positive value of a_T; the results for the negative value may be obtained by using the symmetry (u,v,a) → (-u,-v,-a) of (<ref>).

Precisely at a=a_T, the real part of the dominant eigenvalue of the Jacobian is zero for wavenumbers k = ± k_T = ±(/d)^1/4. 
See the right panel in Fig. <ref>. 
The Turing point a_T marks the boundary between two different (linear) regimes.
On one side, for a ≲ a_T, there are intervals of wavenumbers k, one about each of ± k_T, over which the homogeneous steady state is linearly unstable to plane waves e^ikx since the determinant is strictly negative there. 
On the other side, for a ≳ a_T, the homogeneous state is linearly stable to plane waves with k near k_T (though not to those with k=0  if a<1 by the above).  

Critically, the modes with k near k_T (which is ∼ k_m for small d) are only weakly stable for a ≳ a_T, as may be seen at a=1 in the left panel in Fig. <ref>.
There, λ_+(k_m)=-2√( d), so that the rate constant is asymptotically small, and nonlinear terms also play a central role.
In comparison, the magnitude of this weak stability of the Turing modes (k∼ k_T) is of the same size asymptotically as the magnitude of the weak instability of the Hopf modes (k∼ 0).  Therefore, it is important and useful first to study the two competing mechanisms individually and then to study their interactions. 

In this article, we focus primarily on the spatial dynamics of stationary solutions of (<ref>), which turn out to capture important features of the overall dynamics of (<ref>) for a near a_T. 
In addition, toward the end of the article, we identify some of the rich dynamics created by the interactions of the Turing and Hopf modes.

We study the Turing bifurcation by using the equivalent formulation obtained through the spatial ODE (<ref>),  which we recall is


    δ u_x    = p, 
    δ p_x    = f(u) - v, 
     
        v_x    = q, 
    
        q_x    = (u-a).


and we recall the small parameter is δ = √(d). 
(We refer to <cit.> for the general theory about the equivalence of this spatial ODE formulation.) 

The equilibrium at (u,p,v,q) = (a,0,f(a),0) corresponds to the homogeneous steady state. Linear stability analysis shows that the Jacobian of (<ref>) has a quartet of eigenvalues

    μ = ±1/√(2)δ[
    f'(a) ±√( (f'(a))^2 - 4 δ^2)]^1/2,

which are symmetric about the real- and imaginary-axes in the spectral plane. See Fig. <ref>.



Exactly at the critical value a_T=√(1 - 2 δ√()) given by (<ref>), where f'(a)<0 and (f'(a))^2 = 4 δ^2, the quartet of eigenvalues consists of two coincident pairs of pure imaginary eigenvalues. 
Hence, at a_T, the equilibrium is a reversible, 1:1 resonant Hopf bifurcation point. 
Moreover, this point is non-degenerate, because the two pairs pass through this configuration transversely as a passes through a_T. 
See Fig. <ref>. 
Furthermore, this Turing bifurcation is singular because one pair of the eigenvalues becomes singular in the limit δ→ 0. 
For -a_T < a < a_T, the quartet consists of two pairs of pure imaginary eigenvalues with different imaginary parts, while for a_T < a < √(1+2δ√()) it consists of two pairs of complex conjugate eigenvalues, one with negative real parts and the other with positive real parts of equal magnitude, where the upper bound of this interval corresponds to the parameter value at which the eigenvalues become real.
(See <cit.> for other examples of the relation between a Turing bifurcation in a PDE and a reversible 1:1 resonant Hopf bifurcation in the associated spatial ODE system.)

It is also useful to study the spatial ODE system (<ref>) in the stretched spatial variable y = x/δ, 

    u_y    = p 
    
    p_y    = f(u) - v 
    
    v_y    = δ q 
    
    q_y    = δ (u-a).

This system will be used to analyze  the steep gradients in the spatially periodic solutions that we study, while system (<ref>) will be used for the slowly-varying segments of the spatial periodic solutions.
For all δ>0, the systems of ODEs (<ref>) and (<ref>) are equivalent.




For a spatially periodic solution of (<ref>) with period T, the wavenumber with respect to the stretched y variable is k_y = 2πT. For comparison with the wavenumbers in (<ref>), we work with the wavenumber given by k = 1δ k_y = 2πδ T. Throughout the remainder of the article, we will use the wavenumber k = 2πδ T. 


System (<ref>) has a reversibility symmetry, which it inherits from the x→ -x symmetry of the system of two coupled second-order ODEs that govern time-indepednent solutions of (<ref>). Let 

    ℛ
         = [
        [  1  0  0  0;  0 -1  0  0;  0  0  1  0;  0  0  0 -1 ]],    u = [
        [ u; p; v; q ]],    and   F = [ [              p; 1/3u^3 - u - v;            δ q;        δ (u-a) ]].

Then, ℛ is a reversibility symmetry of (<ref>), because it anti-commutes with the vector field:

    ℛ F ( u) = -  F( ℛ u).

(See for example <cit.> for general results about reversibility in the spatial ODEs governing stationary solutions of reaction-diffusion systems.)  The presence of this reversibility symmetry is also what guarantees the symmetries about the real- and imaginary-axes of the resolvent and spectrum of the operator L obtained by linearizing the vector field F at a stationary state, since u is real-valued and since the reversibility ℛL = -Lℛ implies ( μ𝕀 + L )^-1ℛ = ℛ(μ𝕀 - L )^-1. 

For all >0, system (<ref>) (equivalently (<ref>)) has the following conserved quantity/first integral: 

    𝒢(u,p,v,q,a)=
        1/2 p^2 - 1/2 q^2 -  [f̃(u) + (a-u)v],

where f̃(u) = ∫ f(u)  du = u^4/12 - u^2/2.
(It may be derived for example by setting u_t=0 in the first equation of (<ref>) and multiplying it by u_x, setting v_t=0 in the second equation in (<ref>) and multiplying it by v_x, then subtracting the second equation from the first so that all terms are perfect derivatives, and recalling the definitions of p and q given in (<ref>).)

At the stationary state,  𝒢(a,0,f(a),0,a)=1/12 a^2 ( 6 - a^2 ). 

The existence of small-amplitude, spatially periodic solutions of the ODE system near the critical Turing bifurcation and with wavenumber k near k_T is established by the following proposition. 
The proof is presented in Appendix <ref>, where we use the standard procedure to derive the normal form for the reversible 1:1 resonant Hopf bifurcation in the spatial ODE and then directly apply Theorem 3.21 in Chapter 4.3.3 of <cit.> to derive the conclusions about the system dynamics. 





 For a-a_T sufficiently small, where a_T=√(1-2δ√()), we have:
 

  (i)
For δ^2 < 64/625, the spatial ODE system (<ref>) has a one-parameter family of spatially periodic solutions that surround the symmetric equilibrium (a,0,f(a),0), and it has a two-parameter family of quasi-periodic solutions located on KAM tori. 
The same conclusion also holds for δ^2 > 64/625 but then only when a-a_T<0.

  (ii) 
For a-a_T<0 and δ^2 > 64/625, there exists a one-parameter family of pairs of reversible homoclinics to periodic solutions.

  (iii)
For a-a_T>0 and δ^2 < 64/625, there is a pair of reversible homoclinic orbits to the symmetric equilibrium.

  (iv)
For a-a_T>0 and δ^2 > 64/625, there is just the symmetric equilibrium and no other bounded solutions.



We are most interested in the parameter values for which δ^2 < 64/625, which corresponds to the main case in part (i) and to part (iii). 
These parts of the proposition  establish the existence of one-parameter families of small-amplitude, spatially periodic solutions of (<ref>) near a_T, and the limiting homoclinics.
The Turing bifurcation is sub-critical, and the coefficient on the cubic term in the Ginzburg-Landau PDE is positive <cit.>.
(In contrast, for δ^2 > 64/625, the Turing bifurcation is super-critical, and the coefficient on the cubic term in the Ginzburg-Landau equation is negative.) 



Some representative spatially periodic solutions are shown (in blue and red) in Fig. <ref> for parameter values a close to a_T. 
These solutions exhibit small-amplitude oscillations about the homogeneous state, and their wavenumber k is close to k_T. 
Moreover, the closer  a is to a_T, the smaller the amplitude of the solution, the closer k is to k_T, and the closer the spatial profile is to that of the plane wave e^ik_T x.  
Then, for a slightly further away, nonlinear terms have a dominant effect (black orbit), and a corner begins to form in the orbit, which will develop into a cusp.




From the point of view of Ginzburg-Landau theory for the weakly nonlinear stability of spatially periodic patterns created in Turing bifurcations, the sub-critical case that we focus on is the least studied and understood. In fact, in the sub-critical case, the Ginzburg-Landau formalism typically predicts that the amplitude of the solutions under consideration will grow beyond that of their assumed near-onset magnitude, and thus beyond the region of validity of the Ginzburg-Landau set-up. 
As the analysis below will show, it turns out to be possible to go beyond this magnitude, and our results may thus be expected to shed a fundamental light on the behavior of patterns near a sub-critical Turing bifurcation (see also Section <ref>).





§ SPATIALLY PERIODIC CANARD SOLUTIONS
		


In this section, we present representative numerical simulations of the plethora of spatially periodic canard solutions of (<ref>) that are created asymptotically close to the Turing bifurcation. There are at least four different types of spatial canards, depending on amplitudes and wavenumbers.




 §.§ Small-amplitude spatial canards with Lg wavenumbers
		  


Immediately beyond the parameter regime of essentially regular plane wave solutions, where the nonlinear terms become dominant, we find small-amplitude spatially periodic canard solutions with 𝒪(1) wavenumbers. 
A representative spatial canard is shown in Fig. <ref>.
This solution has period T ≈ 353.366 (k ≈ 1.778), and 
it lies on the level set 𝒢 = ( 23a-14).
The spatial profile (see Fig. <ref>(a)) shows that it lies in a small neighborhood of the u=1 state over the entire period. 
Moreover, on a large portion of the period, the solution slowly increases or slowly decreases. In between these slow regimes, there is a spatial interval on which the solution makes three small-amplitude oscillations around the state u=a. 




During the long portion of slow increase,  the orbit gradually moves up the right branch of the cubic nullcline v=f(u) (black curve), with u slowly increasing. 
See the projection onto the (u,v) plane shown in Fig. <ref>(b).
The orbit then transitions away and oscillates rapidly about the middle branch of the cubic nullcline. 
In the (u,v) projection, the oscillations are manifested by the zigzag upward to the upper left extremum, and then back along the same path (due to the reversibility symmetry of the solution) to the right branch, completing the oscillatory segment on the central spatial interval. 
Along the remainder of the period, the solution slowly moves down along the right branch of the v=f(u) curve, to its minimum value. Along this segment, u slowly decreases. 
The projections onto the (u,p) and (u,q) planes (see panels (c) and (d), respectively) also show that the solution oscillates about the equilibrium state. 

We will show in Section <ref> that system (<ref>) has a folded singularity –in particular a reversible folded saddle-node of type II– which lies on the fold set L^+ and that portions of this  orbit (and of the other small-amplitude solutions of this type) lie along the canards of this singularity. Hence, these are spatial canard solutions. 




 §.§ Small-amplitude spatial canards with small wavenumbers
		


The spatial ODE system (<ref>) also possesses small-amplitude spatial canards with small (𝒪(δ))  wavenumbers. 
A representative solution is shown in Fig. <ref>. 
The solution stays close to the equilibrium state u=a over almost the entire period, as may be seen from the spatial profile of u shown in Fig. <ref>(a).
(Here, the value of a is 𝒪(δ^2) below a=1.)
Then, on an extremely short interval in space, it exhibits a 
small-amplitude excursion, with magnitude of at most 𝒪(δ) (see also the inset in Fig. <ref>(a)). 



During the long spatial interval over which u is close to a, the solution slowly evolves in a neighborhood of the local minimum (see the projection of the solution onto the (u,v) plane shown in Fig. <ref>(b)).
In particular, the solution slowly moves down  the right branch of the nullcline with u decreasing, 
and then it passes underneath  the local minimum 
(right black dot in Fig. <ref>(b)), slowly spiraling in toward the equilibrium at (u,v)=(a,f(a)) (left black dot).
Subsequently, the solution slowly spirals away from the equilibrium, retracing its path (in the (u,v) projection, due to the reversibility symmetry) back underneath the local minimum and toward the right branch of the nullcline, completing the long segment of slow evolution. 
On the short spatial interval that interrupts the long segment of slow evolution, the fast excursion corresponds to the jump away from the right branch to just beyond the middle branch, followed by the return along the same segment (in the projection), due to the reversibility symmetry.
We emphasize that these small-amplitude spatially periodic solutions are also canard solutions, because they also have long segments near the canards of a reversible folded saddle-node singularity, as we will show in Sections <ref> and <ref>.
For completeness, we note that, in frames (c) and (d),  the fast portion of the solution on the short spatial interval corresponds to the segments located to the left of the box marking the first inset, i.e. to the left of approximately u=0.998. 

The dynamics of the solution in the neighborhood of the equilibrium also exhibit a degree of self-similarity. This may be seen in the projections onto the (u,p) and (u,q) planes; see frames (c) and (d), respectively.  
As u and v slowly evolve, their spatial derivatives p and q alternate in sign, and the solution exhibits nested spirals as it slowly approaches –and then slowly recedes from– the equilibrium. 
Three levels of the nested, nearly self-similar spiraling are shown, and the solution exhibits additional, more-deeply nested levels (not shown). 
The nearly self-similar dynamics of this representative solution (and of other small-amplitude solutions with small wavenumbers) will be shown to follow from the self-similarity in a level set of the Hamiltonian that arises naturally in the desingularized reduced system (see Section <ref>).




 §.§ Large-amplitude spatial canards with moderately small wavenumbers
		


In this section, we turn our attention to large-amplitude spatially periodic canard solutions, and we present a class of spatial canards  
that have moderately small (𝒪(√(δ))) wavenumber. 
A representative spatial canard solution of this type is shown in Fig. <ref>. 
Over a large portion of the spatial period, u slowly decreases monotonically from √(3) to just above 1. 
Near its minimum, the u-component exhibits a small-amplitude spike (see the inset in Fig. <ref>(a)).
Subsequently, u slowly increases monotonically  near 1 back up to near √(3). 
Finally, there is a fast downward jump
to u=-√(3) on a narrow spatial interval, followed by a slow segment along which u slowly increases to a local maximum just above -√(3), and then a fast upward jump back to √(3), to complete the period. 



The dynamics along each of the slow and fast segments of the solution observed in the spatial profile of u may be understood from the projections onto the (u,v), (u,p), and (u,q) planes, which are shown in Figs. <ref>(b)–(d), respectively.
The long portion on which u slowly decreases from √(3) (near v=0) down toward 1 corresponds  to the slow evolution near the right branch of the cubic nullcline v = f(u), see  Fig. <ref>(b).
The solution comes close to the equilibrium (right black dot). 
Then, it exhibits a small, localized oscillation about the middle branch of the nullcline, before returning to the neighborhood of the right branch of the v = f(u) nullcline, where it slowly moves upward until it reaches the neighborhood of the v=0 state. 

Subsequently, at the point where the solution reaches the { v=0 } state, a fast jump takes place, along which the solution transitions to the left branch of the v=f(u) nullcline. 
It slowly flows up the left branch of the cubic nullcline until it reaches a maximum, turns around, and then heads back toward the v=0 state. 
Once it reaches the v=0 state, another fast transition is initiated, and the solution returns to the right branch of the cubic. 
The projections of the large-amplitude solution onto the (u,p) and (u,q) phase spaces are shown in Figs. <ref>(c) and (d), respectively. 

These solutions are called spatial canards, 
because they have long segments near the canards of a reversible folded saddle-node singularity, as we will show in Section <ref>. 
Moreover, there are also large-amplitude spatial canards with 𝒪(1) wavenumbers, as we will see in Section <ref>.
These have no loops and only relatively short segments near the true and faux canards of the RFSN-II point.





 §.§ Large-amplitude spatial canards with small  wavenumbers
		


System (<ref>) has a fourth main family of spatially periodic canards, which are large-amplitude, small wavenumber solutions.
In Fig. <ref>, we show a  representative solution with wavenumber k ≈ 0.000628 (T = 1 × 10^6). 
Over most of the spatial period, the solution slowly varies near u=a 
(see Fig. <ref>(a)).
There is a narrow interval in which the solution has five, symmetrically-disposed, large-amplitude spikes, with peaks near u=√(3) and troughs near u=-√(3). 
This narrow interval is magnified in the main inset, so that the individual spikes are visible, including the  pairs of spikes at the edges and the central spike. 
Also, within the main inset, two further insets zoom onto the sharp, small spikes at the outer edges.


 
In the phase space, the solution is near the equilibrium during the large portion on which u near a, see the left black dot in Fig. <ref>(b).  
Then, outside of this large portion of slow evolution, the solution has a narrow interval (see the main inset in Fig. <ref>(a)) in which it makes the two small and five large-scale oscillations. In particular, the small-amplitude, moderately fast oscillation about u=1 (secondary inset) corresponds in phase space to the small-amplitude excursion around the manifold of center points (near the local minimum of the nullcline, as shown in the inset oanel (b)).  
Then, the orbit travels fairly rapidly up the right branch of the cubic nullcline until it reaches a neighborhood of v=0, where the profile has the first of the five large-amplitude peaks. Subsequently, it makes a fast jump toward the left branch of the cubic, oscillating once before it reaches the left branch. 
This oscillation corresponds to the second spike at the left edge of the interval shown in the spatial profile (recall Fig. <ref>(a)). 

Further along, the orbit travels up the left branch of the cubic until it reaches some maximum, where it exhibits a fast jump corresponding to the middle large-amplitude spike of the five in the inset in Fig. <ref>(a). 
Finally, the solution travels back down near the left branch of the cubic nullcline until it reaches the neighborhood of v=0, where another fast jump is initiated, and it makes another large-amplitude oscillation on its way back to the right branch of the cubic. 
These final large-scale oscillations correspond  to the spikes at the right edge in Fig. <ref>(a), thereby completing the periodic solution. 

This large-amplitude, small wavenumber spatial canard exhibits dynamics that are close to being self-similar in the neighborhood of the u=1 state. 
This may be seen in the projections of this solution onto the (u,p) and (u,q) planes (see  Figs. <ref>(c) and (d), respectively). This near self-similarity will be discussed in Section <ref>.






 §.§ Numerical continuation of the families of spatial canards
		



Having presented four main types of spatially periodic canard solutions, we now identify regions in the (a,k) parameter space containing spatially periodic canards, as obtained using numerical continuation.  
(We refer to Appendix <ref> for details about the numerical continuation.) 
We began by finding a stable periodic orbit for a=0, continued this with respect to the period to find the orbits with minimal and maximal periods for a=0, and then performed two-parameter continuation on these in a and the period.
The continuation results are shown in Fig. <ref>. 
For each of four values of δ,  namely for 10^-2, 10^-1, 0.5, and 1, spatially periodic canard solutions exist in the region between two curves.
The upper curve corresponds to minimal period solutions (with maximal k) and the lower curve to the maximal period solutions (with minimal k). 



The solutions with small k have long segments along which they are near the equilibrium,  making successively smaller nested loops, which are close to being self-similar.
The solutions with 𝒪(1) or higher k have shorter or no segments near the equilibrium.
This trend is illustrated for values of a on both sides of a=1 by the representative canards shown in Figs. <ref>–<ref>.
Consider first a=0.999433, so that the equilibrium lies to the left of the local minimum.
The small-amplitude canard with small k shown in Fig. <ref> has a finite sequence of successively smaller, nearly self-similar loops, and hence it has a long segment near the equilibrium (and we note that similar canards exist also for a>1, e.g. along the rightmost blue stalk in Fig. <ref>(c)). 
By contrast, small-amplitude canards with 𝒪(1) k
(one of which is shown for example below in Section <ref>) have few or no loops, and hence have only short or no segments near the equilibrium.
At this same value of a, the large-amplitude canards exhibit the same trend with respect to k, as may be seen for example by comparing the canard with small k shown in Fig. <ref> with a canard with 𝒪(1)  k.
Similar trends are seen for a=1.008, as well as for values of a further from a=1 on both sides.

The lower curves in Figs. <ref> (a), (b), and (c) exhibit several sharp drops in k near the critical Turing value a_T = √(1-2δ√()). 
These sharp drops occur because of the birth of infinite-period (i.e., zero wavenumber),  homoclinic solutions. 
The existence of these homoclinics follows from part (iii) of  Proposition <ref>.
Finally, as illustrated in Fig. <ref> (d), the rightmost edge of the region of spatial periodics converges to the critical Turing bifurcation as δ→8^-25√(), which is in accord with part (i) of Proposition <ref>, since spatially periodic solutions only exist for a<a_T for δ on the other side.
(For reference, we also show the linear stability boundary, δ^2 k^4 + (a^2-1) k^2+ = 0,
of the equilibrium state, which is obtained from DF = 0 (<ref>), as the black dashed curve.)



A more detailed view of the bifurcation structure of the spatially periodic solutions of (<ref>) is  presented in Fig. <ref>, with frames (a)–(c) showing the wavenumbers k of sizes k ≫ 1, k=𝒪(1),
and k ≪ 1, respectively. 
First, the same diagram from Fig. <ref>(a) is shown in Fig. <ref>(a), now also with  curves of constant 𝒢 (red curves), fixed at the special value 𝒢 = ( 23a-14). These particular spatially periodic solutions only exist at small and 𝒪(1) wavenumbers. At the 𝒪(1) scale of wavenumbers (Fig. <ref>(b)), the fine structure of the spatially periodic solutions becomes more apparent. Each of the different colored curves (red, orange and blue) correspond to different branches of spatially periodic solutions on the level set 𝒢 = ( 23a-14). The red branch emerges from the Turing bifurcation (black marker), continues into large a, decreases to small k, and remains at small wavenumbers (horizontal segments).
The detailed structure of the nearly horizontal segments for small wavenumbers is shown in Fig. <ref>(c). It can be seen that there are even more branches of spatially periodic solutions with 𝒢 = ( 23a-14). Some of these branches exhibit sharp drops in k and continue in toward the homoclinic solution along k=0 for a>a_T. 
Each of these branches in Fig. <ref>(c) contains spatially periodic canard solutions. 



Overall, the numerical continuation results show that system (<ref>) has many different types of spatially periodic canard solutions in the (a,k,d) space for a wide range of values of >0. 
We have found these types of spatially periodic canards for other values of , including =0.8.
We add that, in Section <ref>, 
the solutions and bifurcations along these branches will be studied in detail by using the results obtained in the analysis of the fast system, the desingularized reduced vector field, and the RFSN-II singularity (Sections <ref>–<ref>).

To conclude this section, we note that the regions shown in Fig. <ref> correspond to what are called `existence balloons' in <cit.>, i.e., regions in (parameter,wave number)-space in which spatially periodic patterns exist. 
Naturally, the eventual goal would be to determine the Busse balloon, the sub-balloon of stable spatially periodic patterns <cit.>. 
The Busse balloon plays a central role in pattern forming systems in fluid mechanics, reaction-diffusion systems, ecological models, and many other scientific and engineering problems <cit.> and is also expected to play a similar role in understanding the impact of the canard-type solutions of this work on the dynamics exhibited by PDE (<ref>).




§ FAST SYSTEM (LAYER PROBLEM)




In this section, we study the fast system (layer problem) of the spatial ODE, which is obtained by taking the limit δ→ 0 in (<ref>): 

    u_y    = p 
    
      p_y    = f(u) - v 
    
      v_y    = 0 
    
      q_y    = 0.

Hence, v and q are constants, and the fast system is independent of a.
Moreover, for each fixed v, the fast system is Hamiltonian, 

    u_y = ∂ H_ fast/∂ p,    
      p_y = -∂ H_ fast/∂ u,
         with   
      H_ fast(u,p;v)= 1/2p^2 - f̃(u) + uv.

The functional form of H_ fast is equivalent to that obtained by reducing 𝒢 to the fast system and scaling by a factor of . We recall from (<ref>) that f̃= 1/12u^4 - 1/2 u^2.





 §.§ The critical manifold of the fast system




In the four-dimensional (u,p,v,q) space, the fast system (<ref>) has a two-dimensional critical manifold, 

    S := { (u,p,v,q)
            : u∈ℝ,   p=0,   v=f(u),    q ∈ℝ},

which corresponds to the union of all of the fixed points of (<ref>).
Moreover, given that the cubic function f(u)=1/3u^3 - u has three segments separated by the local extrema at (u,v)=(± 1, ∓2/3), the critical manifold S has three branches:

    S_s^-
               = { (u,p,v,q) ∈ S   :  u=u_- < -1,     p = 0,    v=f(u_-),    q ∈ℝ}, 
    
         S_c 
                = { (u,p,v,q) ∈ S :  u=u_0,    -1 < u_0 < 1,    p=0,     v=f(u_0),     q ∈ℝ}, 
    
         S_s^+ 
                = { (u,p,v,q) ∈ S : u=u_+>1,    p=0,    v=f(u_+),     q ∈ℝ}.

S_s^± are hyperbolic saddle invariant manifolds, because the fixed points (u_±,0) are saddle fixed points of the (u,p) subsystem for each v>-2/3 and each v < 2/3, respectively. Similarly, S_c is a non-hyperbolic, elliptic invariant manifold, because the fixed point (u_0,0) of the (u,p) subsystem is a nonlinear center for each -2/3 < v < 2/3.
The three invariant manifolds  are illustrated in Fig. <ref>.



The fold curves are

    L^±=
    {
    (u,p,v,q) | u = ± 1, p=0,
    v = ∓2/3, q ∈ℝ}.

They form the boundaries of the invariant manifolds, where S_s^± lose normal hyperbolicity.
Together with the three invariant manifolds, they provide the following natural decomposition of the critical manifold, 

    S  = S_s^- ∪ L^- ∪ S_c ∪ L^+ ∪ S_s^+.









 §.§ Homoclinics of the fast system



We prove the existence of homoclinic orbits in the fast system, which will be useful in the deconstruction of spatially periodic canard solutions.

  
For each fixed v ∈ (-23,0), there exists a homoclinic solution of (<ref>) corresponding to the intersection W^s(u_+) ∩ W^u(u_+), where u_+ ∈ S_s^+ is the root of v = f(u) for which u_+ >1. 
Similarly, for each fixed v ∈ (0,23), there exists a homoclinic solution of (<ref>) corresponding to the intersection W^s(u_-) ∩ W^u(u_-), where u_- ∈ S_s^- is the root of v = f(u) for which u_- <-1. 



We examine the fast system Hamiltonian along the cross-section { p =0 }:


    H_ fast(u,0; v) = u v - ( 112u^4-12u^2 ).


For a fixed v in the interval (0,23), H_ fast(u,0; v) is a quartic polynomial (Fig. <ref>) in u with local maxima at u_± and a local minimum at u_0.
Recall that u_± are saddle equilibria of (<ref>) and u_0 is a center. 



The Hamiltonian H_ fast(u,0; v) is monotonically increasing for all u<u_-, decreasing for u_-<u<u_0, increasing for u_0<u<u_+, and decreasing for u>u_+.
Moreover,


    H_ fast(u_+,0; v) > H_ fast(u_-,0; v) > H_ fast(u_0,0; v),


so that, by the Intermediate Value Theorem, there exists a u_* ∈ (u_0,u_+) such that H_ fast(u_-,0; v) = H_ fast(u_*,0; v).
Since the only equilibrium between u_- and u_* is the nonlinear center u_0, the orbits with u ∈ [u_-,u_*] are closed. 
Hence, there exists a closed orbit of (<ref>) that connects the fixed point (u_-,0) to a point (u_*,0) and then back to itself. 

The proof for v ∈ (-23,0) is similar.





 §.§ Heteroclinics in the intersections Lg




In this section, we establish

 
In the (u,p,v,q) space, W^u(S_s^-) ∩ W^s(S_s^+) and W^u(S_s^+) ∩ W^s(S_s^-). Moreover, these intersections lie in the plane { v=0 }, and they are transverse.


The stable and unstable manifolds are shown in Fig. <ref>,
along with the heteroclinic orbits that lie in their transverse intersections.




For v=0,
the saddles (-√(3),0)
and (√(3),0)
lie on the level set
H_ fast (±√(3),0;0)=3/4.
They are connected by heteroclinic orbits
in the (u,p) plane,
as may be verified directly
from the functional form
of the Hamiltonian. 


Next, for values of v near –but not equal to– zero,
the saddles (u^±,0) lie on different level sets of the Hamiltonian.
In particular, for each fixed v>0 and close to zero,
H_ fast < 3/4 at the saddle (u^-,0) on S_s^-,
and H_ fast > 3/4 at the saddle (u^+,0) on S_s^+.
Similarly, for each fixed v≲  0 ,
H_ fast > 3/4 at the saddle (u^-,0) on S_s^-,
and H_ fast < 3/4 at the saddle (u^+,0) on S_s^+.
Therefore,
W^u(S_s^-) and W^s(S_s^+)
intersect transversely
in the plane { v=0 },
as do the invariant manifolds
W^u(S_s^+) and W^s(S_s^-).
 




 §.§ The cusp point at Lg




For v=-2/3,
the fast system (<ref>) has a saddle-node point at (u,p)=(1,0). The Jacobian at that point has a double zero eigenvalue, with
Jordan block given by 

    [  
    [ 0 1; 0 0 ]].

In the phase plane, it is a cusp point with a single-branched stable manifold
and a single-branched unstable manifold. In parametrized form, the functions whose graphs give these manifolds have the form 


    p = ∓ C (u-1)^3/2,
 

for some C>0, respectively, and they lie on the level set H_ fast(u,p)=- 1/4.

This cusp point is the point we will desingularize.  The stable manifold will correspond to one fixed point on the equator of the blown-up hemisphere, and the unstable manifold to another fixed point. The dynamics of the orbits near the cusp point will then be determined by studying the dynamics over the hemisphere.





§ SLOW SYSTEM ON LG, WITH THE FOLDED AND ORDINARY SINGULARITIES
 



In this section, we study the slow system induced by the spatial ODE system on the two-dimensional critical manifold S = { u ∈ℝ,   p=0,   v=f(u),   q ∈ℝ}, recall (<ref>). We derive the locations and stability types of the folded and ordinary singularities. 

Recall that the slow manifold S consists of the fixed points of the fast system in the limit δ→ 0 (i.e., of (<ref>) with δ=0).
One obtains the slow dynamics on S in the slow variable by differentiating the constraints p=0 and v=f(u) that define S:

    [
            [     0     1; f'(u)     0 ]]
            [
            [ u_x; p_x ]]
            =
            [
            [   0; v_x ]].

The adjoint operator is

    adj[
    [     0     1; f'(u)     0 ]]
    =
    [
    [      0     -1; -f'(u)      0 ]].

Hence, left multiplying (<ref>)
by the adjoint and recalling the third and fourth components of
(<ref>), we find

    -f'(u)
            [
            [ u_x; p_x ]]
               = [
            [ -q;  0 ]] 
    [
            [ v_x; q_x ]]
               = [
            [     q; (u-a) ]].


Next, we desingularize the reduced system on S
by rescaling time with a factor of f'(u).
Let dx = f'(u) dx_d.
In this new, dynamic time variable,
the reduced system (<ref>) is equivalent to

    [
            [ u_x_d; p_x_d; v_x_d; q_x_d ]]
            = [
            [           q;           0;     f'(u) q; f'(u) (u-a) ]].

Finally, we observe that, on S, the constant solution for p is p=0.
Also, the equation for v decouples (and may be solved by quadrature).
Hence, we arrive at the desingularized reduced system

    u_x_d   = q 
     
            q_x_d   =  f'(u) (u-a).


This is the main system studied in this section; it defines the folded and ordinary singularities. 
We note that on S_s^± the direction of the flow of (<ref>) is the same as that of (<ref>), whereas on S_c the flow direction is reversed since f'(u)<0 on S_c.
Also, (<ref>) is Hamiltonian with

    H_ desing(u,q;a)= -1/2q^2 + ( 1/4u^4
               - 1/3a u^3 - 1/2u^2 + au ),

where H_ desing is induced by (<ref>), i.e., H_ desing = .𝒢|_S, and u_x_d=-∂ H_ desing∂ q and q_x_d=∂ H_ desing∂ u.

The ordinary and folded singularities of the desingularized reduced system (<ref>) are found as follows.
Off the fold set, both components of the vector field vanish at (u,q)=(a,0). Hence, (<ref>) has an ordinary singularity (equilibrium) at

    E = { u=a, q=0 }.

Then, on the fold set where f'(u)=0, there are two  folded singularities 

    M^± = { u = ± 1, q=0 },

where the vector field in (<ref>) vanishes.

We classify the linear stability types of the ordinary and folded singularities using the eigenvalues of the Jacobian of (<ref>),

    [
            [                    0                    1; f”(u) (u-a) +  f'(u)                    0 ]]



























The classifications are listed in Table <ref>.






The phase plane of (<ref>)
is sketched for different values of a>0
in Fig. <ref>. (Those for a<0 may be obtained using the symmetry (u,q,a) → (-u,-q,-a).)
The orbits lie on the level sets
of the Hamiltonian H_ desing
given by (<ref>). 
The equilibrium E lies on the level set H_ desing= a^2/12 ( 6 - a^2 ), and the folded singularities M^± lie on the level sets H_ desing=( ±2/3a - 1/4). Furthermore, we observe that ∂ H_ desing/∂ u (u,0) =  (u^2-1)(u-a)> 0
for u ∈ (-1,a) in the case 0<a<1 and for u ∈ (-1,1)
in the case a≥ 1. 


By Fenichel theory <cit.>, the normally hyperbolic saddle sheets, S_s^±, of the critical manifold perturb to nearby invariant slow manifolds, S_s,δ^±, of (<ref>). Moreover, the flow of (<ref>) restricted to S_s,δ^± is a small 𝒪(δ) perturbation of the reduced flow on S_s^±. Canard theory can then be used to establish that the singular true and faux canards of the folded saddle persist as maximal canard solutions for 0<δ≪ 1 provided the eigenvalues of the folded saddle remain bounded away from zero <cit.>. 



The saddle slow manifolds, S_s,δ^+, and associated maximal canards are shown in Fig. <ref>. The numerical method used to compute these is outlined in Appendix <ref>.
The saddle slow manifolds consist of two types of solutions. The first type are the solutions that approach the neighbourhood of the fold curve L^+ and subsequently escape via the fast dynamics. The other type are the solutions that turn away from the fold region and remain on S_s,δ^+. The separatrices that divide between the two types of solutions are the maximal canards of the folded saddle. In the RFSN-II limit, the saddle slow manifolds and the canard dynamics become complicated (see Fig. <ref>(b).)







§ DESINGULARIZATION OF THE REVERSIBLE FSN-II POINT LG




In this section, we focus on the dynamics of the system for a values near a=1 and desingularize the reversible FSN-II point M^+ located at (u, p, v, q, a, δ) = (1, 0, -2/3, 0, 1, 0).
First, we translate coordinates so that the RFSN-II point is shifted to the origin.
Let


    (u,p,v,q,a,δ)
    =(1+ũ,p̃,-2/3 + ṽ, q̃,
    1+ã,δ̃).

In these new coordinates,
the spatial system (<ref>) becomes

    ũ_y    = p̃
    p̃_y    = ũ^3/3 + ũ^2 - ṽ
    ṽ_y    = δ̃q̃
    q̃_y    = δ̃ (ũ-ã) 
    ã_y    = 0 
    δ̃_y    = 0.


Also, the conserved quantity (<ref>) may be rewritten as

    𝒢(u,p,v,q,a) 
        =  𝒢̃(ũ,p̃,ṽ,q̃,ã) + ( 2/3ã + 5/12)  
        =  𝒢̃(ũ,p̃,ṽ,q̃,ã) + ( 2/3a - 1/4),

where 𝒢̃(ũ,p̃,ṽ,q̃,ã)= 1/2 (p̃^2 - q̃^2 ) 
- ( 1/12ũ^4 + 1/3ũ^3 + (ã-ũ)ṽ). 


The point 
(ũ,p̃,ṽ,q̃,ã,δ̃)
=(0,0,0,0,0,0), which corresponds to the RFSN-II point at a=1

(recall Table <ref>), 
is a nilpotent point
of (<ref>).
Hence, in order to unfold the dynamics, we desingularize this point 
by introducing the additional dependent variable r
and by dynamically rescaling the variables:

    ũ = √() r^2 u̅,  p̃ =  r^3 p̅,  ṽ =  r^4 v̅,   q̃ = ^3/2  r^3 q̅,   ã = √() r^3 a̅,   δ̃ = r^2 δ̅.

Here, (u̅,p̅,v̅,q̅,a̅, δ̅) ∈𝕊^5.
The factors of  in (<ref>) were chosen 
for symmetry in the equations for ũ and ṽ.

We will examine system
(<ref>)
in two coordinate charts:
the entry/exit chart

    K_1 = {u̅=1 }

and the rescaling (or central) chart

    K_2 = {δ̅=1 }.

As is customary in articles on geometric desingularization (a.k.a. blow up), in chart K_i, i=1,2,
the variable □̅
will be denoted by □_i.
In chart K_1,
the rescaling (<ref>)
is given by

    ũ   = √() r_1^2,  p̃ =  r_1^3 p_1,  ṽ =  r_1^4 v_1,  q̃ = ^3/2  r_1^3 q_1,  ã = √()r_1^3 a_1,  δ̃ = r_1^2 δ_1.

Then, in chart K_2,
the rescaling (<ref>)
is given by

    ũ   = √() r_2^2 u_2,  p̃ =  r_2^3 p_2,  ṽ =  r_2^4 v_2,  q̃ = ^3/2  r_2^3 q_2,  ã = √() r_2^3 a_2, δ̃ = r_2^2.



This section is organized as follows.
In Subsection <ref>, we present the analysis in the rescaling chart K_2, identifying the key algebraic solution in this chart and the self-similarity of the level set on which it lies. 
These results also justify the range of validity of the asymptotics, as stated in the Introduction, since  we have ã = ^3/2δ̃^3/2 q_2 by the rescalings (<ref>) and δ̃=δ by (<ref>).
Subsection <ref> contains the analysis in the entry/exit chart K_1.
Then, the transition maps between the two charts are presented in Subsection <ref>, and these are used to track the key solution from chart K_2 into chart K_1. 
Next, in Subsection <ref>, the transverse intersection of the main center-unstable and center-stable manifolds is shown, thereby establishing the formula (<ref>) for the critical value a_c(δ).
Finally, Subsection <ref> presents sp,e analysis of the smoothness of key solutions. 





 §.§ Rescaling (or central) chart Lg




To obtain the governing equations
in chart K_2,
we substitute (<ref>)
into (<ref>):

    du_2/dy   = r_2 √() p_2 
    dp_2/dy   = r_2 (u_2^2 - v_2) + 1/3√() r_2^3 u_2^3 
    dv_2/dy   = r_2 √() q_2
    dq_2/dy   = r_2 (u_2 - r_2 a_2) 
    da_2/dy   = 0.

Here, we used dr_2/dy=0 in K_2,
which follows since δ̅≡ 1 implies
dδ̅/dy=0.
Then, introducing the rescaled independent variable
y_2 = r_2 y
and letting prime denote
d/dy_2,
we obtain the desingularized vector field in chart K_2,

    u_2'    = √() p_2 
    p_2'    = u_2^2 - v_2 + 1/3√() r_2^2 u_2^3 
    v_2'    = √() q_2
    q_2'    = u_2 - r_2 a_2 
    
      a_2'    = 0.


For each a_2, this system is Hamiltonian

    [ u_2^'; q_2^'; p_2^'; v_2^' ]   = [  ∂ H_2/∂ p_2;  ∂ H_2/∂ v_2; -∂ H_2/∂ u_2; -∂ H_2/∂ q_2 ] = J  ∇_(u_2,q_2,p_2,v_2) H_2

Here, J is the 4× 4 anti-symmetric matrix [  0  𝕀; -𝕀  0 ], 𝕀 is the 2× 2 identity matrix, and H_2 is the Hamiltonian:

    H_2(u_2,p_2,v_2,q_2;r_2)
       = 1/2√()( p_2^2 - q_2^2 )
       + (u_2 - r_2 a_2) v_2
         - ( 1/3 u_2^3 + 1/12√() r_2^2 u_2^4 ),

in which r_2 is a constant parameter.
Furthermore, this Hamiltonian H_2 is directly related through the transformation  (<ref>) to the conserved quantity via:
𝒢̃  (ũ,p̃,ṽ,q̃,ã) 
= ^5/2 r_2^6 H_2 (u_2, p_2, v_2, q_2; r_2),  
recall (<ref>).



 
The hyperplane { r_2=0 }
is invariant in chart K_2.
On this hyperplane,
the system reduces to

    u_2^'   = √() p_2 
    p_2^'   = u_2^2 - v_2 
    v_2^'   = √() q_2
    q_2^'   = u_2.





  §.§.§ The key algebraic solution in chart Lg




In this section, we establish the existence of a key algebraic solution that acts as a separatrix in phase space.


The system (<ref>) possesses an explicit algebraic solution given by


    Γ_0 = (u_2(y_2), p_2(y_2), v_2(y_2), q_2(y_2))
             = (
               1/12√() y_2^2,
               1/6 y_2,
               1/144 y_2^4 - 1/6,
               1/36√() y_2^3
             ).


Moreover, the orbit Γ_0 has the following two halves:


    Γ_0^- = . Γ_0 |_y_2≤ 0   and   Γ_0^+ = . Γ_0 |_y_2 ≥ 0.


The segment Γ_0^- corresponds to the true canard of the RFSN-II singularity and the segment Γ_0^+ corresponds to the faux canard of the RFSN-II. 


As we will see, both halves of the orbit Γ_0  play central roles in the dynamics induced by the RFSN-II singularity and also in the singular limit of the spatially periodic canards. 
It lies in the zero level set of the Hamiltonian H_2(u_2,p_2,v_2,q_2;0), as may be verified by direct substitution and observing that all y_2-dependent terms sum to zero.  
For reference, it may also be parametrized by p_2: u_2 = 3 √() p_2^2, v_2 = 9 p_2^4 - 1/6, and q_2 = 6√() p_2^3.






  §.§.§ Geometric role of the algebraic solution Lg


Here, we study the unperturbed problem (<ref>) and provide a geometric interpretation of the algebraic solution Γ_0. 
Recall from (<ref>) that the system (<ref>) has the Hamiltonian 


    H_2(u_2,p_2,v_2,q_2;0) = 12√() (p_2^2-q_2^2)+u_2 v_2 - 13 u_2^3,
 

and the canard solutions, Γ_0^±, lie in the zero level contour { H_2(u_2,p_2,v_2,q_2;0) = 0 }. 
By using this constraint to eliminate the v_2 variable, the dynamics restricted to the { H_2(u_2,p_2,v_2,q_2;0) = 0 } level set are given by 


    u_2^'   = √() p_2 
     
    p_2^'   = 23u_2^2 + 12u_2√() (p_2^2-q_2^2) 
     
    q_2^'   = u_2.


To remove the singularity at u_2 = 0, we desingularize the vector field via the transformation dy_2 = u_2 dη_2, which gives


    du_2/dη_2   = √() u_2 p_2 
    dp_2/dη_2   = 23u_2^3 + 12√() (p_2^2-q_2^2) 
    dq_2/dη_2   = u_2^2.


The system (<ref>) is topologically equivalent to (<ref>) for u_2 >0 and has opposite orientation for u_2<0. Moreover, the system (<ref>) has the symmetry 


    (u_2,p_2,q_2,η_2) → (u_2,-p_2,-q_2,-η_2),


which it inherits from the reversibility symmetry (<ref>) of the original problem. 

The system (<ref>) possesses two lines of equilibria, 


    ℒ_± = { u_2 = 0, p_2 = ± q_2: q_2 ∈ℝ}.


The line ℒ_- has hyperbolic and center spectra, σ_ s/u = { -√()q_2,-√()q_2} (where s/u depends on the sign of q_2) and σ_c = { 0 }, respectively, with corresponding subspaces 


    𝔼^ s/u( ℒ_- ) = span{[ 1; 0; 0 ], [ 0; 1; 0 ]}   and   𝔼^c( ℒ_- ) = span[  0; -1;  1 ].


Thus, the half-lines 


    ℒ_-^s := { (u_2,p_2,q_2) ∈ℒ_- : q_2 >0 }   and   ℒ_-^u := { (u_2,p_2,q_2) ∈ℒ_- : q_2 < 0 }
 

are center-stable and center-unstable, respectively.
Similarly, the line ℒ_+ has hyperbolic and center spectra, σ_ s/u = {√() q_2, √() q_2 } and σ_c = { 0 }, respectively, with corresponding subspaces 


    𝔼^ s/u( ℒ_+ ) = span{[ 1; 0; 0 ], [ 0; 1; 0 ]}   and   𝔼^c ( ℒ_+ ) = span[ 0; 1; 1 ].


Thus, the half-lines 


    ℒ_+^s := { (u_2,p_2,q_2) ∈ℒ_+ : q_2 <0 }   and ℒ_+^u := { (u_2,p_2,q_2) ∈ℒ_+ : q_2 >0 }


are center-stable and center-unstable, respectively.

In Appendix <ref>,
we prove the following proposition that establishes the key geometric properties of the algebraic solution, Γ_0, and its stable and unstable manifolds in phase space:


The restriction of the system (<ref>) to the half-space { u_2 ≥ 0} possesses three classes of heteroclinic solutions. 


  * A class 1 heteroclinic is a solution that connects a point (0,-q̃_2,q̃_2) ∈ℒ_-^u for some q̃_2<0 to a point (0,-q̂_2,q̂_2) ∈ℒ_-^s for some q̂_2>0.

  * A class 2 heteroclinic is a solution that connects a point (0,-q̃_2,q̃_2) ∈ℒ_-^u for some q̃_2<0 to a point (0,q̂_2,q̂_2) ∈ℒ_+^s for some q̂_2<0.

  * A class 3 heteroclinic is a solution that connects a point (0,q̃_2,q̃_2) ∈ℒ_+^u for some q̃_2>0 to a point (0,-q̂_2,q̂_2) ∈ℒ_-^s for some q̂_2>0.

The stable manifold, W^s(Γ_0^-), of the algebraic solution Γ_0^- is the phase-space boundary between class 1 and class 2 heteroclinics. 
The unstable manifold, W^u(Γ_0^+), of the algebraic solution Γ_0^+ is the phase-space boundary between class 1 and class 3 heteroclinics. 



The results of this proposition are illustrated in Fig. <ref>. 


  





 §.§ Entry/exit chart Lg




We substitute (<ref>) into (<ref>) and desingularize the vector field by transforming to a dynamic independent variable y_1= r_1 y. Hence, the system in chart K_1 is

    dr_1/dy_1   = 1/2√() p_1 r_1 
    d δ_1/dy_1   = -√()p_1 δ_1 
    dp_1/dy_1   = 1-v_1-3/2√()p_1^2+1/3√()r_1^2 
    dv_1/dy_1   = √()(-2p_1 v_1 + δ_1 q_1 ) 
    dq_1/dy_1   = -3/2√() p_1 q_1+δ_1 (1-r_1 a_1 ) 
    da_1/dy_1   = -3/2√() p_1 a_1.


This system has a series of invariant hyperplanes defined by { r_1=0}, {δ_1=0}, { a_1=0}, and their intersections. In this section, we present the analysis of system (<ref>) in the hyperplane {r_1=0}∩{δ_1=0}∩{ a_1=0}, where the dynamics are the simplest.  Then, we go directly to the analysis of system (<ref>) in the full (r_1,δ_1,p_1,v_1,q_1,a_1) space. 
Intermediate results that are needed for this full analysis are presented in Appendix <ref>, where we analyze system (<ref>) in the following sequence of invariant hyperplanes: 
{δ_1=0}∩{ a_1=0}, 
{ r_1=0}∩{ a_1 = 0 }, 
{ r_1 = 0 }∩{δ_1 = 0 }, 
{δ_1=0 }, 
and { r_1=0},

respectively.
For the remainder of this section and in Appendix <ref>, we will use an overdot to denote the derivative with respect to y_1: = d/dy_1.

In the invariant hyperplane
{r_1=0}∩{δ_1=0}∩{ a_1 = 0 },
the equations (<ref>) reduce to

    ṗ_1    = 1 - v_1 - 3/2√() p_1^2 
    v̇_1    = -2 √() p_1 v_1 
    q̇_1    = -3/2√() p_1 q_1.

Hence, the q_1 equation decouples from the equations for (p_1,v_1). 
This system (<ref>) has two invariant lines. The first is a line of equilibria

    ℓ = { r_1=0,  δ_1=0,   p_1=0,   v_1=1,   q_1 ∈ℝ,   a_1=0 }.

It is a line of saddle points,
since the eigenvalues are

    λ_s = -√(2)^1/4,   λ_u = √(2)^1/4,    and   λ_c = 0.

The associated eigenvectors are

    w_s = [ [      2√(2);      4^1/4; 3 ^1/4 q_1 ]],
       w_u = [ [     -2√(2);      4^1/4; 3 ^1/4 q_1 ]],
       and    w_c = [ [ 0; 0; 1 ]].

See Fig. <ref>. 
Therefore, by standard center manifold theory (see for example <cit.>), system (<ref>) has a one-dimensional center manifold W^c(ℓ) in the invariant hyperplane
{ r_1 = 0 }∩{δ_1 = 0 }∩{ a_1=0 }, and it has one-dimensional fast stable and unstable fibers in the transverse directions.
This manifold, which we will show is embedded in a larger center manifold in the full system, and the fast stable and unstable fibers,
will play central roles in the dynamics of the full system (<ref>).



The second invariant line in the phase space of equations (<ref>) is

    I = {   r_1=0,  δ_1=0,  
                          p_1 ∈ℝ,  
                          v_1=0,   q_1=0,   a_1=0
                }.

On I, the vector field reduces to ṗ_1 = 1 - 3/2√()p_1^2, and there are two equilibria:

    E_±= { r_1 =0,  δ_1=0,   
                          p_1 = ±√(23)^-1/4,  
                          v_1=0,   q_1=0,   a_1=0
                }.


The equilibrium E_+ is stable with eigenvalues -√(6)^1/4, -2√(2/3)^1/4, -√(3/2)^1/4,
and eigenvectors

    w_1^s = [ [ 1; 0; 0 ]],
       w_2^s = [ [    -√(3); √(2)^1/4;        0 ]],
       w_3^s = [ [ 0; 0; 1 ]].

The equilibrium E_- is unstable with eigenvalues
√(3/2)^1/4,
2√(2/3)^1/4,
and √(6)^1/4,
and eigenvectors

    w_1^u = [ [ 0; 0; 1 ]],
       w_2^u = [ [     √(3); √(2)^1/4;        0 ]],
       w_3^u = [ [ 1; 0; 0 ]].


The equilibria E_± are shown in the (p_1,v_1) subsystem in 
Fig. <ref>(a) 
and in the (p_1,v_1,q_1) system in Fig. <ref>(b).
The number of arrows on the principal components of the stable and unstable manifolds indicates  the relative magnitudes of the eigenvalues. 
The orbits in this phase plane are defined implicitly by the one-parameter family  √() p_1^2 - ( 2/3 - 2 v_1 ) -c v_1^3/2 = 0, where c denotes the parameter.
Overall, the dynamics of (<ref>) are organized by the line of equilibria ℓ and the invariant line I. 

Now, we turn to the invariant sets and dynamics of the full system (<ref>), recalling that some of the intermediate results are in Appendix <ref>.
The line I given by (<ref>) is invariant for the full system (<ref>), and E_± are again fixed points on I, recall (<ref>).

The equilibrium E_+ is a hyperbolic saddle with stable and unstable spectra 

    σ^s_+ = { -√(6)^1/4,
    -2√(2/3)^1/4,
    -√(3/2)^1/4,
    -√(3/2)^1/4,
    -√(2/3)^1/4}   and  σ_+^u = {1/√(6)^1/4}.

At E_+, the stable and unstable subspaces are given by 

    𝔼^s =  span{[ [ 0; 0; 1; 0; 0; 0 ]],
    [ [        0;        0;    -√(3); √(2)^1/4;        0;        0 ]],
    [ [ 0; 0; 0; 0; 1; 0 ]],
    [ [ 0; 0; 0; 0; 0; 1 ]],
    [ [    0; ^1/4;    0;    0; √(6);    0 ]]
    }   and  𝔼^u =  span{[ [ 1; 0; 0; 0; 0; 0 ]]
                 }.


The other equilibrium, E_-, is also a hyperbolic saddle with stable and unstable spectra given by

    σ^s_- = 
        {-1/√(6)^1/4}   and  σ_-^u = {√(2/3)^1/4,
        √(3/2)^1/4,
        √(3/2)^1/4,
    2√(2/3)^1/4, √(6)^1/4}.

At E_-, the stable and unstable subspaces are

    𝔼^s =  span{[ [ 1; 0; 0; 0; 0; 0 ]]
                 }   and  𝔼^u =  span{[ [     0; -^1/4;     0;     0;  √(6);     0 ]],
    [ [ 0; 0; 0; 0; 1; 0 ]],
    [ [ 0; 0; 0; 0; 0; 1 ]],
    [ [        0;        0;     √(3); √(2)^1/4;        0;        0 ]],
    [ [ 0; 0; 1; 0; 0; 0 ]]
    }.

Thus, the blow-up transformation (<ref>) splits the non-hyperbolic equilibrium state (a,0,f(a),0) into the pair of hyperbolic fixed-points E_±. Moreover, by (<ref>) E_+ has a five-dimensional stable manifold and a one-dimensional unstable manifold, with the latter being in the r_1 direction,
and by (<ref>) E_- has a one-dimensional stable manifold (in the r_1 direction) and a five-dimensional unstable manifold. 

Now, in addition to the equilibria E_± and invariant line I, the full system (<ref>) also has an important, three-dimensional manifold of equilibria:

    𝒮_1 = { r_1 ∈ℝ,  δ_1=0,    p_1=0,   v_1=1 + 1/3√() r_1^2,   q_1 ∈ℝ,   a_1 ∈ℝ}.

It contains the line ℓ, and it is a manifold of saddle fixed points, since the eigenvalues are

    λ_s = -√(2 √() +   r_1^2),   λ_u = √(2 √() +   r_1^2), 
              λ_c = 0, 0, 0, 0.

The associated stable and unstable eigenspaces are


    𝔼^s =  span{[ [             -√() r_1;                    0;   2 √(2√() +  r_1^2); 4/3√()(3 + √()r_1^2);            3 √() q_1;              3√()a_1 ]]
                 },
                   𝔼^u =  span{[ [             -√() r_1;                    0; - 2 √(2√() +  r_1^2); 4/3√()(3 + √()r_1^2);            3 √() q_1;              3√()a_1 ]] },


and the associated center subspace is 


    𝔼^c =  span{[ [ 0; 0; 0; 0; 1; 0 ]],
    [ [ 0; 0; 0; 0; 0; 1 ]],
    [ [        3;        0;        0; 2√() r_1;        0;        0 ]],
    [ [             0; 2 + √() r_1^2;           q_1;             0;             0;             0 ]]
    }.

Therefore, application of  center manifold theory
(see for example <cit.>) establishes that (<ref>)
has a four-dimensional
center manifold

    M = W^c(ℓ).

This manifold contains the surface 𝒮_1 of equilibria.
It also contains the three-dimensional center manifold N 
in the invariant hyperplane { r_1=0 }, see (<ref>) in Appendix <ref>.
Moreover, the branch of N in p_1<0 and δ_1>0 is unique.

The above analysis establishes the following lemma:

The following properties hold for the system (<ref>).
  
  
  * There exists a stable invariant foliation, ℱ^s, with base M and one-dimensional fibers. For any constant c>-√(2), the contraction along ℱ^s during an interval [0,X] is stronger than e^cX. 
  
  * There exists an unstable invariant foliation, ℱ^u, with base M and one-dimensional fibers. For any constant c< √(2), the expansion along ℱ^u during an interval [0,X] is stronger than e^cX. 
  






 §.§ Forward and backward asymptotics of Lg in chart Lg




In this section, we study the forward and backward asymptotics (as y_2 →±∞) of the “Vo connection" Γ_0, given by (<ref>), along with its tangent vectors. 

We start with the coordinate transformation maps between the two charts. 
For δ_1>0, the coordinate transformation from K_1 to K_2 is

    κ_12 (r_1,δ_1,p_1,v_1,q_1,a_1)
                 = (r_2,u_2,p_2,v_2,q_2,a_2) 
               = (
                   r_1 √(δ_1), 1/δ_1,
                   p_1/δ_1^3/2, v_1/δ_1^2,
                   q_1/δ_1^3/2,
                   a_1/δ_1^3/2).

For u_2>0,
the coordinate transformation from K_2 to K_1 is

    κ_21 (r_2,u_2,p_2,v_2,q_2,a_2)
               = (r_1,δ_1,p_1,v_1,q_1,a_1) 
             = (
                 r_2 √(u_2), 1/u_2, p_2/u_2^3/2,
                 v_2/u_2^2, q_2/u_2^3/2,
                 a_2/u_2^3/2).


Now, for the forward asymptotics of Γ_0, we use κ_21 to find 

    lim_y_2 →∞κ_21(Γ_0)
            = lim_y_2 →∞(
              0, 12/√() y_2^2, 4√(3)/^3/4 y_2^2,
              1-24/ y_2^4, 2√(3)/3^1/4, 0
              )
            = (
              0, 0, 0, 1,
              2√(3)/3^1/4, 0
              ).

Hence,
on the invariant hyperplane { r_1=0 },
orbits on Γ_0
are forward asymptotic
to the fixed point
(0,0,0,1,q^+,0)
with q^+=2√(3)/3^1/4,
on the line of equilibria ℓ.

Similarly, the backward asymptotics of Γ_0
are given by:

    lim_y_2 → -∞κ_21(Γ_0)
            = lim_y_2 → -∞(
              0, 12/√() y_2^2, -4√(3)/^3/4 y_2^2,
              1-24/ y_2^4, -2√(3)/3^1/4, 0 
              )
            = (
              0, 0, 0, 1,
              -2√(3)/3^1/4, 0 
              ).

Hence, on the invariant hyperplane
{r_1=0}, solutions on Γ_0
are backward asymptotic to the fixed point
(0,0,0,1,q^-,0) with
q^-=-2√(3)/3^1/4,
on the line of equilibria ℓ.



The tangent vectors of Γ_0 in chart K_1 at the equilibria (0,0,0,1,q^±,0) are given by


    lim_y_2 →±∞d/dy_2κ_21(Γ_0)/‖d/dy_2κ_21(Γ_0) ‖ 
    = lim_y_2 →±∞(0,∓ 24 √() y_2^2,-8√(3)^1/4 y_2^2,± 96,0 ,0)/8√(3)√(48+√()y_2^4(1+3√())) 
    = ( 0, ∓√(3)^1/4/√(1+3√()),-1/√(1+3√()), 0, 0, 0).


Hence, the tangent vectors of Γ_0 in chart K_1 at the equilibria (0,0,0,1,q^±,0) are tangent to the center manifolds M. In fact, at these equilibrium points, the tangent vectors are parallel to the generalized eigenvector of the center subspace.





 §.§ Intersection of center-unstable and center-stable manifolds at Lg




In this section, we present the derivation of (<ref>): 



   For 0 < δ≪ 1, the center-unstable manifold W^cu(0,0,0,1,q^-,0) coincides with the center-stable manifold W^cs(0,0,0,1,q^+,0) for 
   
    a_c(δ) = 1 - 5 /48δ^2 + 𝒪(δ^3).



This lemma is proven in Appendix <ref>. 
It involves a Melnikov type calculation.

An immediate consequence of the lemma is that, for all δ>0 sufficiently small, there is a maximal canard orbit exactly at a_c(δ)=1 - 5 /48δ^2 + 𝒪(δ^3).
Moreover, the orbit is unique up to translations since  the ODEs are autonomous. 
Of all the canard orbits in the spatial ODE system (<ref>) that pass through a neighborhood of the cusp point, the maximal canard has the longest segments near the stable and unstable manifolds of the cusp point, all the way out to where u=+√(3).
This critical value is the analog in spatial dynamics of the critical parameter value at which the slow attracting and repelling manifolds coincide to all orders, and hence the explosion of temporal limit cycle canards occurs, in planar fast-slow ODEs.





 §.§ Calculation of the center-unstable and center-stable manifolds for Lg



In this section, we establish the following proposition: 

 
The branches of the critical manifold corresponding to the level set 𝒢 = (23a-14) perturb to invariant slow manifolds. These slow manifolds are smooth to at least 𝒪(δ^2) in p and smooth to at least 𝒪(δ^3) in q, for 

    a_c(δ) = 1 + a_2 δ^2 + 𝒪(δ^3).

for any real a_2.


The proof uses the iterative scheme devised in <cit.> to calculate the invariant slow manifolds of (<ref>), and determine the key parameter values for which the slow manifolds intersect. 
To do this, we follow the approach in <cit.> in which the parameters are carefully chosen to remove poles in the expansions of the slow manifolds. 

Proof of Proposition <ref>. 
We first translate the reversible FSN-II to the origin via the linear shift


    ũ = u-1,   p̃ = p,   ṽ=v+23,   q̃ = q,   ã = a-1.


Dropping the tildes, the dynamics are governed by 


    u̇   = p 
    ṗ   = f(u) - v 
    v̇   = δ q 
    q̇   = δ (u-a)


where f(u) = u^2+13u^3. In this formulation, the conserved quantity is 


    𝒢 = (u-a)v - f(u) + 12 p^2 - 12q^2,


where f(u) = 13u^3+112u^4 is an antiderivative of f, and the level set 𝒢 = (23a-14) becomes 𝒢 = 0. By eliminating the q-variable, we reduce the study of the 2-fast/2-slow system (<ref>) with constraint 𝒢=0 to the study of the 2-fast/1-slow system 


    u̇   = p 
    ṗ   = f(u) - v 
    v̇   = √(2) δ √((u-a)v-f(u)+p^2).


Here, we focus on the positive root for q, and then obtain results for the negative root for q by using the symmetry (<ref>). 

For the 2-fast/1-slow system (<ref>), the critical manifold is given by 


    S_0 = { p = 0 , v = f(u) }


with fold points at (u,p,v)=(0,0,0) and (u,p,v)=(-2,0,43), and canard points at (u,p,v,a)=(0,0,0,0) and (u,p,v,a)=(-2,0,43,-2). The canard point at the origin corresponds to the RFSN-II point of interest. 

We assume that the invariant slow manifold, S_δ, has a graph representation in u, and that it can be expanded as a power series in δ,


    S_δ = { (u,p,v,q) : p = ∑_k=0^∞ p_k(u) δ^k,    v = ∑_k=0^∞ v_k(u) δ^k,    q = ∑_k=0^∞ q_k(u) δ^k }


We also expand the parameter a as a power series


    a = ∑_k=0^∞ a_k δ^k.


Substituting these expansions into (<ref>) and equating coefficients of like powers of δ, we find that the leading 𝒪(1) terms are given by 


    p_0 = 0,    v_0 = f(u),    and    q_0 = √(2)√((u-a_0)f(u)-f(u)),


corresponding to the critical manifold. Proceeding to the next terms in the expansion, we find that the 𝒪(δ) terms are given by 


    p_1 = √(2)√((u-a_0)f(u)-f(u))/f^'(u),    v_1 = 0,    and    q_1 = 0.


To remove the singularity at u=0 in p_1, we set a_0 = 0 so that every term in the argument of the square root has a factor of u^2. In that case, the coefficients become


    p_1 = √(2)√(u ( 14u + 23))/2+u,    v_1 = 0,    and    q_1 = 0.



Continuing to higher order, we find that the 𝒪(δ^2) terms are given by 


    p_2 = -a_1 √(23) (3+u)/√(u) (2+u) √(8+3u),    v_2 = -(4+u)/3(2+u)^3,    and    q_2 = -√( u) 4a_2(2+u)^3(3+u)-(10+3u)/2√(6) (2+u)^3 √(8+3u).


To eliminate the singularity at u=0 in p_2, we set a_1 = 0, so that 


    p_2 = 0,    v_2 = -(4+u)/3(2+u)^3,    and    q_2 = -√( u) 4a_2(2+u)^3(3+u)-(10+3u)/2√(6) (2+u)^3 √(8+3u).



The 𝒪(δ^3) terms are given by 


    p_3 = √()( -12a_2(2+u)^5(3+u)+(9u^3+54u^2+64u-40) )/6√(6) (2+u)^6 √(u)√(8+3u),    v_3 = 0,    and    q_3 = 0.


Here, there is a singularity at u=0 in p_3, however, we choose not to enforce smoothness of the p-component of the slow manifold and leave a_2 free. 

Therefore, the proposition is proven, since p is smooth up to and including  𝒪(δ^2) and q up to and including 𝒪(δ^3) for any choice of a_2. 

We add that the 𝒪(δ^4) terms are given by 


    p_4 = -√(2/3)a_3 (3+u)/√(u)(2+u)√(8+3u),    v_4 = P_6(u)/3(2+u)^8,    and    q_4 = -√() P_12(u)/24 √(6)
       (u+2)^8 u^3/2 (8+3u)^3/2.


where P_6(u) and P_12(u) are, respectively, the 6^ th and 12^ th order polynomials


    P_6(u)    = -a_2(2+u)^5(4+u)+(3u^3+18u^2+14u-34) 
     
    P_12(u)    = 24  a_2 (9 u^3+61 u^2+110 u+32) (u+2)^5 + 48 a_2^2 u (u+3)^2 (u+2)^8    + 
            48 a_4
       u^2 (3 u^2+17 u+24) (u+2)^8- ^2 (513 u^5+4608 u^4+12168 u^3+5120 u^2-12400 u-2560)


To cancel one factor of u in the denominator of q_4, we must set a_2 = -5/48. 
Thus, upon reverting to the original coordinates and parameters, Fraser-Roussel iteration shows that in order for the slow manifolds to be smooth to at least 𝒪(δ^2) in p and smooth to at least 𝒪(δ^3) in q, the parameters must satisfy


    a = 1 - 5/48δ^2 + 𝒪(δ^3).


consistent with the values obtained from the Melnikov analysis, recall (<ref>). □



On the level set 𝒢 = (23a-14), one can alternatively get smoothness to at least 𝒪(δ^3) in p and to at least 𝒪(δ^3) in q for 


    a_c2(δ) = 1 - 5/144δ^2 + 𝒪( δ^3 ).

This follows from the calculations in the proof of Proposition <ref>. At 𝒪(δ^3), the pole at u=0 in p_3 is removed by seeking a value for a_2 such that the numerator in the expression for p_3 has u as a factor, namely a_2 = -5144. The lack of smoothness in the singular limit is visible in the (u,p) and (u,q) projections for small δ), with corners and cusps (in the projections) developing. There appears to be a discrete set of critical values of a marking the transitions between solutions with different numbers of loops, and different numbers of points of nonsmoothness in the limit of small δ. 





§ ANALYZING THE GEOMETRY OF THE SPATIAL CANARDS
 


In this section, we use the analytical results from  Sections <ref>–<ref>,  for the analyses of the fast system, the desingularized reduced system, and the desingularization of the RFSN-II singularity, respectively, to deconstruct and understand the spatially periodic canard solutions.
This analysis reveals how the RFSN-II point at a=a_T  and its canards, as identified in Section <ref>, are responsible for the creation of these spatial canard solutions.  
The geometry of small-amplitude spatial canards is analyzed in Section <ref>, and that of  large-amplitude spatial canards in Section <ref>.




 §.§ Geometry of the Small-Amplitude Spatial Canards
 

We present the geometry of small-amplitude spatially periodic canards.
For illustration, we use the solution shown in Fig. <ref>, and we recall also the representative small-amplitude canards in Figs. <ref> and <ref>. 



In the singular limit, the small-amplitude orbit consists of the following three segments, all confined to the level set {𝒢 = ( 23a-14) }: 



  * A short slow segment from the folded saddle at (u,p,v,q) = (1,0,-23,0)  (blue ducky) that closely follows the faux canard, Γ_0^+, of the folded saddle up to some nearby v-value (near the red ducky).


  * A fast segment following a homoclinic solution (see Proposition <ref>) of the layer problem until it returns to a neighborhood of S_s^+ (near the green ducky). 


  * A short slow segment that closely follows the true canard, Γ_0^-, of the folded saddle until it returns to the folded saddle. 


In this manner, the singular solutions  of small-amplitude, spatially periodic canards with 𝒪(1) wavenumbers are completely understood through the dynamics of the RFSN-II singularity, its true and faux canards, and some invariant manifold theory, following the results of the analysis in Sections <ref>–<ref> .
The geometry of other small-amplitude spatially periodic canard solutions with 𝒪(1) wavenumbers, such as that shown in Fig. <ref>, is similar. 

The geometry of small-amplitude canards with asymptotically small wavenumbers may also be understood in a similar manner. 
These are very close in profile to the small-amplitude canards with 𝒪(1) wavenumbers over most of the period, but they also exhibit nested, successively smaller, nearly  self-similar twists (loops in certain projections) that are confined to the neighborhood of the equilibrium and the folded saddle
(recall for example the solution shown in Fig. <ref>).
These nearly self-similar loops arise for small δ due to scale invariance of the zero level of the Hamiltonian H_2 in the rescaling  chart K_2. 






 §.§ Geometry of the Large-Amplitude Spatial Canards




In this section, we use the analytical results 
derived in Sections <ref>–<ref>, 
about the fast system, the desingularized reduced system, and the blowup analysis, to fully deconstruct the large-amplitude, spatially periodic canard solutions of (<ref>). 
We choose the solution from Fig. <ref> for the deconstruction, and we review some of the main properties of that solution again in  Fig. <ref>.




The singular limit solution is obtained by concatenating orbit segments from the fast and slow subsystems, (<ref>) and (<ref>). 
The singular orbit consists of seven segments, all of which are constrained to the level set {𝒢 = ( 23a-14) }. 
For convenience, we denote the level set {𝒢 = ( 23a-14) } restricted to orbits of the layer problem by 𝒢_ fast and we denote the level set {𝒢 = ( 23a-14) } restricted to the slow flow on S by 𝒢_ slow.



  * Slow segment on S_s^+ along the faux canard Γ_0^+ from the folded saddle at (u,p,v,q)=(1,0,-23,0) (grey ducky) to the upper right corner point (red ducky) at


    (u,p,v,q)=( √(3),0,0,√(23(3-2a))).


(For a = 0.998512, the upper right corner point is located at (u,p,v,q) ≈ (1.7321,0,0,0.2586).) 


  * Fast jump along the heteroclinic orbit in the transverse intersection W^u(S_s^+) ∩ W^s(S_s^-) in { v = 0} (see Proposition <ref>) from the take off point at the upp[er right corner (red ducky) to the touch down point at the upper left corner (green ducky) at 


    (u,p,v,q)=( -√(3),0,0,√(23(3-2a))).
 This heteroclinic consists of three segments that occur in rapid succession, one from a neighborhood of S_s^+ to near S_s^-, followed by another back to the neighborhoof of S_s^+, and then the third is forward asymptotic to S_s^-; it lies in a secondary intersection of the two invariant manifolds.  

(For a = 0.998512, this corner point is located at (u,p,v,q) ≈ (-1.7321,0,0,0.2586).)

  * Slow segment on S_s^- from the touch down point (green ducky) to the local maximum (in u) of the 𝒢_ slow contour (gold ducky) at 


    (u,p,v,q) = ( u_-,0,f(u_-),0 ),
 

where u_- = 13( 2a-3-2√(a^2+3a)). (For a = 0.998512, this gold ducky is located at (u,p,v,q) ≈ (-1.6664,0,0.1239,0).)  

  * Fast jump consisting of the homoclinic orbit of the layer problem (see Proposition <ref>) from the take off point (gold ducky) to the local maximum of the homoclinic (cyan ducky) 
and then to the touch down point on S_s^- at (u,p,v,q) = ( u_-,0,f(u_-),0 ) (magenta ducky).  

  * Slow segment on S_s^- from the touch down point (magenta ducky) until the solution reaches the { v = 0 } hyperplane at 
(u,p,v,q)=(-√(3),0,0,-√(23(3-2a)))
(blue ducky).

  * Fast jump in the transverse intersection W^u(S_s^-) ∩ W^s(S_s^+) in { v = 0 } (see Proposition <ref>) from the take off point (blue ducky) to the touch down point at the lower right corner (yellow ducky) on S_s^+ at 


    (u,p,v,q) = (√(3),0,0,√(23(3-2a)))

This heteroclinic consists of three segments in rapid succession and lies in a secondary intersection of the manifolds.

(For a = 0.998512, the lower right corner point is located at (u,p,v,q) ≈ (1.7321,0,0,-0.2586).)

  * Slow segment on S_s^+ along the true canard, Γ_0^-, of the RFSN-II point
from the touch down point (yellow ducky) up to the cusp  (grey ducky), thus completing the singular cycle. 






For the solution shown in Fig. <ref>, the connection from S_s^+ to S_s^- (and vice versa) is a secondary heteroclinic and features additional large-amplitude spikes not present in  the singular limit. 
We observed other large-amplitude, small wavenumber solutions without additional spikes that  correspond to primary heteroclinics. 
The deconstruction of such solutions is similar to that presented above. 
In addition, we observed solutions that have more of these spikes along the outer edges of the orbit.   






§ ISOLAS OF CANARD-INDUCED SPATIAL PERIODICS

 


In this section, we employ many of the results derived in Sections <ref>-<ref>, about the fast/layer problem, the desingularized reduced vector field, and the geometry of the RFSN-II singularity, to study the bifurcations of the spatially periodic canard solutions along a representative isola in the (a,k) parameter plane (see Fig. <ref>). 
The branches of the isola have been color coded (green, orange, blue, and red) according to the nature of the spatially periodic solutions that lie on each segment.


 




 §.§ Growth of canard segments
 

We analyze the growth of the canard segments for solutions on the green branch of the isola (see Fig. <ref>(a)). 
These solutions are four-stroke, spatial relaxation cycles. 



As illustrated in Fig. <ref>(c) and (d), in the singular limit, the relaxation cycles consist of:  


  * A slow segment on S_s^+ that starts at the touchdown point (u,p,v,q)=(√(3),0,0,-√(2(1-23a))), flows down the critical manifold toward the fold set L^+ to the turning point (local minimum) at (u,p,v,q) = (u_ turn,0,13u_ turn^3-u_ turn,0) where u_ turn = 23a-1+23√(a^2+3a), and then flows back up the critical manifold to the takeoff point at (u,p,v,q)=(√(3),0,0,√(2(1-23a))).


  * A fast jump at v = 0 corresponding to the heteroclinic of the layer problem that takes the orbit from the takeoff point at (u,p,v,q)=(√(3),0,0,√(2(1-23a))) to the touchdown point at (u,p,v,q)=(-√(3),0,0,√(2(1-23a))). 

  * A slow segment that flows along S_s^- from the touchdown point at (u,p,v,q)=(-√(3),0,0,√(2(1-23a))) up to the turning point (local maximum) at (u,p,v,q) = (u_ turn,0,13u_ turn^3-u_ turn,0) where u_ turn = 23a-1-23√(a^2+3a), and then down to the takeoff point at (u,p,v,q)=(-√(3),0,0,-√(2(1-23a))). 

  * Another fast jump at v = 0 corresponding to the heteroclinic of the layer problem that takes the orbit from the takeoff point at (u,p,v,q)=(-√(3),0,0,-√(2(1-23a))) to the touchdown point at (u,p,v,q)=(√(3),0,0,-√(2(1-23a))).


At the right edge of this branch, there is a    
saddle-node bifurcation near a ≈ 1.5 (Fig. <ref>(a)), which marks the transition to the blue branch.
As the parameter a is decreased toward a=1, the slow segment on S_s^+ with q<0 converges to the true canard of the folded saddle-node, and the slow segment on S_s^+ with q>0 converges to the faux canard of the folded saddle-node. 
Then, as a is further decreased along this branch to a=0, the slow segment on S_s^- with q>0 grows toward the true canard of the folded saddle at (u,p,v,q)=(-1,0,23,0), and the slow segment on S_s^- with q<0 approaches the faux canard of the folded saddle.




 §.§ Spike formation
 

In this section, we analyze the sequence of canards that exhibit spike formation along the orange branch of the isola (see  Fig. <ref> (a)). 
The solutions here consist of large-amplitude cycles with spikes initiated on S_s^-. The geometry of these solutions is similar to that of the solution presented in Section <ref>.  



At a=0 (left green triangle), the solution is a four-stroke, spatial relaxation cycle, with slow segments given by the true and faux canards of the FS points on L^±, and fast jumps in { v=0 }.
Then, moving left to right along the orange branch, as the value of a is increased from a=0, an extra spike forms at the local maximum of the slow segment on S_s^-. 
It is a given to leading order by a homoclinic of the fast subsystem. 
The size (in the (u,p) projection) of the extra spike increases with a, whilst the v-height at which it occurs decreases away from v = 23. 
At the saddle-node bifurcation at a ≈ 1.5, the v-level of the fast jump is at v = 0, the fast jump corresponds to the heteroclinic of the layer problem  
(recall Section <ref>), and the slow segment on S_s^- has vanished. 
We also note that at the other end, there is a saddle-node bifurcation near a ≈ -0.0002, where the orange branch transitions to the green branch 





 §.§ Two-stroke relaxation cycles
 

In this section, we analyze the solutions along the blue branch of the isola (Fig. <ref>(a)).
This branch emanates from one of the homoclinic orbits (k=0) near a=1. 
The homoclinic solution has small amplitude and is  centered on the steady state. Continuation of the solution away from the homoclinic limit causes the wavenumber to increase to a value close to the critical Turing value k_T. As this occurs, the corresponding spatially periodic solution remains small amplitude and confined to a small neighborhood of the fold L^+.



As the value of a is increased along this branch, the solutions develop into two-stroke relaxation oscillations that consist of a slow segment on S_s^+ and a fast (near homoclinic) jump. 
The slow segment on S_s^+ starts at a `touchdown' point at (u,p,v,q) = (u_ touchdown(v),0,v,q_ touchdown(v)), where 


    u_ touchdown(v)    = 2^1/3/(3v+√(9v^2-4))^1/3 + (3v+√(9v^2-4))^1/3/2^1/3, 
     
        q_ touchdown(v)    = -√()√(23a-14-au_ touchdown+12u_ touchdown^2+13au_ touchdown^3-14u_ touchdown^4),


for v ∈ (-23,0). It flows down the 𝒢= (23a-14) contour to the turning point (u,p,v,q) = (u_ turn,0,v_ turn,0), where 


    13 u_ turn^3 - u_ turn = v_ turn = 2/81(2a-3)(2a-3+2√(a^2+3a))(2a+3+2√(a^2+3a)).


Then, it flows back up the 𝒢= (23a-14) contour to the `takeoff' point at (u,p,v,q) = (u_ takeoff(v),0,v,q_ takeoff(v)), where u_ takeoff = u_ touchdown and q_ takeoff = - q_ touchdown.
At the takeoff point, a fast jump is initiated which connects the takeoff point to the touchdown point via the homoclinic of the layer problem that exists at that value of v (see Proposition <ref>). 

When the blue branch of the isola reaches the saddle-node bifurcation at a ≈ 1.49881, the v-level at which the fast jump occurs is v = 0,  and the homoclinic jump mechanism switches to a heteroclinic jump mechanism.  





 §.§ Two-stroke double-loop cycles
 

Finally, we analyze solutions along the red branch of the isola, which emerges from the orange branch  at the saddle-node bifurcation near a ≈ 1.5. 
The solutions on the red branch of the isola are two-stroke relaxation oscillations, like those in Section <ref>. However, the fast segments of these orbits consist of two near-homoclinic cycles. 



As the parameter a is decreased along this branch  away from the saddle-node bifurcation near a≈ 1.5, the v-level at which the fast jump occurs decreases,  and the amplitudes of the double-loop homoclinics shrink. 
Then, near a=1, the red branch of the isola becomes nearly vertical. Along this nearly vertical  part of the isola, the amplitudes of the solutions are very small, and the solution becomes homoclinic to the equilibrium near the fold L^+ (still with small amplitude). 




§ THE NEARLY SELF-SIMILAR DYNAMICS OF SOME SPATIAL CANARDS




In this section, we discuss some of the spatial canard solutions that exhibit nearly self-similar dynamics. Recall for example the canards shown in Figs. <ref> and <ref>. First, we study self-similarity in the equations in the rescaling chart, (<ref>), and then we study it in the full fourth order system (<ref>).

In the coordinate chart K_2, the nearly self-similar dynamics for 0<δ≪ 1 may be understood as follows. 
For δ=0, the equations are Hamiltonian with

    H_2(u_2,p_2,v_2,q_2;0)
       = 1/2√()( p_2^2 - q_2^2 )
       + u_2 v_2
         - 1/3 u_2^3,

where we recall (<ref>).
Examination of the level set { H_2 = 0 } reveals that it is scale-invariant. 
In particular, on the hyperplane {r_2=0}, which corresponds to the singular limit of δ=0, the zero level set of the Hamiltonian H_2(u_2,p_2,v_2,q_2;0) is invariant under the scaling
(u_2,p_2,v_2,q_2) → (ζũ_2, ζ^3/2p̃_2, ζ^2 ṽ_2, ζ^3/2q̃_2) for any real number ζ.
In fact, H_2(u_2,p_2,v_2,q_2;0)=ζ^3 H_2(ũ_2,p̃_2,ṽ_2,q̃_2;0).
The projection of this level set onto the (u,q) plane has infinitely many self-crossings and nested loops. It is infinitely self-similar.

Then, for 0<δ≪ 1, this scale invariance is broken. Nevertheless, for sufficiently small values of δ, the solutions exhibit nearly self-similar dynamics.
This nearly self-similar dynamics manifests itself in the maximal canards of the full equations with 0 < r_2 ≪ 1 in chart K_2. 
A boundary value problem was set up using these equations, together with the boundary conditions
that enforce the symmetry (u_2,p_2,v_2,q_2) → (u_2,-p_2,v_2,-q_2) and an integral constraint to preserve the Hamiltonian.
(See Appendix <ref> for details.) 
The continuation method starts with the algebraic solution Γ_0 at r_2=0 and a_2=0 given by (<ref>). 
That solution is continued into r_2>0 until r_2=√(δ). 
Then, it is continued in a_2 to yield the maximal canards for the system in chart K_2.
Throughout, the constraint is imposed to conserve the value of the Hamiltonian.
The key observation is that these maximal canards  approach the infinitely self-similar structure in the limit.
The continuation of the maximal canards is shown in the plane of a_2 and the L^2 norm in Fig. <ref> (a), and orbits of the maximal canards are shown in Figs. <ref> (b)–(g) for a sequence of a_2 values. 



The bifurcation curve undulates about the line a_2=0, as shown in Fig. <ref>(a).
The magnitude of the undulations decreases as the L^2 norm decreases. 
Also, there appears to be a self-similarity in the curve. 
For example, if one takes the lower part of the bifurcation curve shown in Fig. <ref>(a) below the second local extremum on the left, stretches that vertically and horizontally so that it has the same height as the curve shown, and overlays that stretched version onto the curve shown, then they look almost identical.

The nearly self-similar structure of the bifurcation diagram persists into the full fourth-order spatial system (<ref>), see Fig. <ref>. 
The undulations decay and the branch eventually converges to a single value of a. For = 1 and δ = 0.05, we find that the a value to which the branch converges is a_c,num≈ 0.99972881. Comparing this with the critical value, a_c(δ), from the Melnikov analysis (see formula (<ref>)), we have that a_c(0.05) ≈ 0.99973958. Hence, 


    a_c(0.05)-a_c,num = 𝒪(δ^4).








§ THE SPATIAL CANARDS OF LG ARE  ANALOGS IN SPATIAL DYNAMICS OF TEMPORAL LIMIT CYCLE CANARDS IN FAST-SLOW ODES




The spatially periodic canards in the PDE system (<ref>) are spatial analogs of the time-periodic canards known to exist in fast-slow systems of ODEs <cit.>. 
(For reference, we recall that for example in the van der Pol ODE u̇=v-f(u), v̇=(a-u), the Hopf bifurcation occurs at a=1 and the canard explosion is centered on a_ ODE∼ 1-1/8- 3/32^2, asymptotically as → 0, recall <cit.>.)


First, the Turing points a_T are reversible 1:1 resonant Hopf bifurcations in the spatial ODE system, as shown in Section <ref>. 
The reversible Hopf points are the analogs in spatial dynamics of the Hopf bifurcation points that occur in (temporal) fast-slow ODEs.
Second, families of small-amplitude, spatially periodic solutions are created in the Turing bifurcation.
These solutions have wavenumbers close to k_T and profiles close to the plane wave e^ik_Tx (recall the red and blue orbits in Fig. <ref>).
These spatially periodic solutions are the analogs of the small-amplitude, temporally oscillating solutions that are created in (singular) Hopf bifurcations in fast-slow systems of ODEs. 
Third, the critical value a_c(δ), recall (<ref>), asymptotically close to a_T is the analog in spatial dynamics of the canard explosion value (a_ ODE for the van der Pol ODE, for example) in fast-slow oscillators, which is asymptotically close to the Hopf point.
Fourth, on each side of the maximal spatial canards, there are families of (non-maximal) spatial canards, just as there are (non-maximal) temporal canards (headless ducks and ducks with heads) on each side of maximal limit cycle canards.
In this manner, the maximal spatial canards act as separatrices that locally partition the phase space into regions of distinct behaviour, and the bifurcations of maximal canards delimit distinct modes of activity in the parameter space, just as is the case for maximal temporal canards. 
Fifth, the singularity responsible for the spatial canards is an RFSN-II point. 
This is the analog of the canard point in the temporal fast-slow systems, which may be viewed as an FSN-II point in the extended (u,v,a) phase space of the kinetics ODE.
Sixth, in carrying out the geometric desingularization analysis of the RFSN-II point in the spatial ODE system (<ref>), we identified a key algebraic separatrix solution Γ_0, recall (<ref>), in the rescaling chart that is the key segment of the singular limit of the spatial canards. 
It is the spatial analog of the parabolic separatrix solution in the rescaling chart (see <cit.>) in the analysis of the explosion of temporal limit cycle 
canards in fast-slow ODEs.

There are also interesting differences between the spatial canards and the temporal limit cycle canards, which arise due to differences between the spatial ODEs and the temporal (kinetics) ODEs, as well as due to the increased dimension of the phase space. 
In the classical explosion of temporal canards, the left and right branches are one dimensional attracting slow manifolds and the middle branch a one dimensional repelling slow manifold. 
Moreover, the maximal canards exist for the parameter values when these two slow manifolds coincide to all orders. 
In contrast, in the spatial ODE, the left and right branches are two dimensional saddle slow manifolds, and the middle branch a two dimensional elliptic slow manifold, recall Fig. <ref>. 
Further, the maximal canards exist when the single-branched stable and unstable manifolds of the cusp point, which are the true and faux canards of the RFSN-II point, 
continue into each other to all orders.





§ PDE DYNAMICS
   

In Sections <ref>–<ref>, we focused primarily on studying the rich family of stationary, spatially periodic canard patterns created by Turing bifurcations in the singularly perturbed van der Pol PDE (<ref>) for values of a near the  Turing bifurcation a_T = √(1-2δ√(ε)), as well as numerically for a∈ [0,a_T). 
A natural –and necessary– next question is: which of these patterns are observable, i.e., which are  stable as solutions of the PDE?  
This is in general a hard and deep question.
In this section, we will scratch the surface of this challenge, presenting a number of PDE simulations that give some hints about what kind of (stability) results may be expected.

We fix (,δ) at (0.1,0.01), which implies that a_T = 0.996832....
Moreover, we consider two primary values of a: a=1 > a_T, so that the background state (a,f(a)) is stable as solution of (<ref>) and a=0.99 < a_T for which (a,f(a)) has undergone the Turing bifurcation. 
We consider the dynamics generated by (<ref>) from the point of view of the onset of pattern formation, i.e., we consider initial conditions (u_0(x),v_0(x)) that are close (in L^1 norm) to (u,v) ≡ (a,f(a)), and we look for the a priori small (∼ close to (a,f(a)))  patterns generated by (<ref>). 
For a=1, |a-a_T|= 0.00316...
This implies, from the point of view of weakly nonlinear stability theory/the Ginzburg-Landau approach (see for example <cit.>), that the magnitudes of the patterns expected near the Turing bifurcation are of the order of  √(|a-a_T|) = 0.0562.... 
Therefore, we only work with initial conditions (u_0(x),v_0(x)) that satisfy |u_0(x)-a|, |v_0(x)-f(a)| ≤ 0.10 < 2√(|a-a_T|) uniformly on the x-interval under consideration. 

Now, since δ^2 < 64/625 for 𝒪(1) values of >0 and 0<δ≪ 1, the Ginzburg-Landau equation associated to the Turing bifurcation has a positive Landau coefficient i.e., a positive coefficient on the cubic term  (Section <ref>).  
Thus, one does not expect small stable patterns for a < a_T.
In fact, one expects the amplitude of the patterns generated by the Turing bifurcation to grow beyond the validity of the Ginzburg-Landau approach <cit.>. 
Then, for a > a_T, only the trivial state (corresponding to (a,f(a))) will be stable. 
Moreover, its domain of attraction will be small, in the sense that solutions with initial data that is  𝒪(√(|a-a_T|)) close to (a,f(a)), but not sufficiently close, are also expected to outgrow the Ginzburg-Landau domain of validity.   



Fig. <ref> shows the final attracting patterns of two numerical simulations with a=1.00 that both start with localized initial conditions at most 0.10 away from the (1,-2/3) background state, ı.e., the initial perturbations of (a,f(a)) = (1,-2/3) can be considered to be sufficiently small so that the dynamics can be expected to be within the range of validity of the Ginzburg-Landau approach. 
In fact, the only difference between the two  simulations is the widths (in x) of the `plateaus' where u_0(x) differs from a=1. 
In both cases, the initial conditions are outside the domain of attraction of (u(x,t),v(x,t)) ≡ (1,-2/3) (by construction: further reducing the amplitude of the initial conditions brings (u_0(x),v_0(x)) into the domain of attraction of (1,-2/3)). 
Indeed, the dynamics generated by (<ref>) depart from the Ginzburg-Landau region near (1,-2/3), and the magnitude of the pattern grows to 𝒪(1). 
In both cases, a (stationary) localized homoclinic multi-front pattern appears as attractor.
Depending on the width of the initial plateau, it is either a 1-pulse/2 front or a 2-pulse/4-front pattern, see Fig. <ref>. 
In the second row of Fig. <ref>, the orbits  of the final attracting states are plotted in the (u,v) phase space (parameterized by x). 
This further confirms that the attractors are of the large-amplitude spatial canard types studied in Sections <ref>, <ref>, and <ref>. 



Next, we consider exactly the same setting as in the simulations of Fig. <ref> but decrease a to a=0.99 < a_T. 
As a consequence, the pattern dynamics of (<ref>) must be significantly different, since the background state (u(x,t),v(x,t)) ≡ (a,f(a)) has become unstable. 
In Fig. <ref>, the attracting pattern is plotted that emerges from the simulations with a=0.99 from initial conditions that are similar to those in the second column of Fig. <ref> (in the sense that the deviations from (a,f(a)) are identical in both cases). 
Here also, the attractor is a stationary large-amplitude homoclinic pattern, but now not homoclinic to the homogeneous state (a,f(a)), but to a spatially periodic pattern – with a well-selected wavenumber. 

This periodic pattern appears step-by-step, on a long time scale. 
Up to t ≈ 380, the pattern is  `transitionally stable', in the sense that it remains very similar to the attractor observed for a=1.00 in Fig. <ref>. 
Then, a fast growing spike appears that evolves into a new localized 2-front pulse – or better: two of these, one each for x > 0 and for x <0 – between t = 390 and t=400 (see the first row of Fig. <ref>). 
In the second row of Fig. <ref>, the same patterns are plotted, now again in (u,v)-space. 
The (spatial) canard dynamics near the point (a,f(a)) ≈ (1,-2/3) seem to play a dominant role in the formation and shape of the spike. 
Each `new period' of the periodic pattern to which the final state is asymptotic is formed in the same manner – with long periods of `transient stability' in which the pattern only evolves marginally and slowly. 
Here, we leave it as a subject for future research to go deeper into the relation between the spatial dynamics near (1,-2/3) and the formation of the periodic pattern and especially the mechanism by which the periodic pattern –that must be an element in a one-parameter family of spatially periodic patterns– is selected.    



Overall, the simulations of the van der Pol PDE (<ref>) that we have conducted so far indicate that there are substantial classes of initial conditions that are in the basins of attraction of (large-amplitude) stationary patterns of the types  considered in the analysis here. 
Moreover, this behavior appears to be related to the width of plateaus in both components u_0(x) and v_0(x) of the initial data.
If one keeps all other aspects as in the simulations shown, then there appears to be a critical threshold for the width of the plateaus.
In particular, data for which the plateau width is below the threshold evolve to these attractors. Whereas, for data in which the plateau width exceeds the threshold, the system no longer exhibits stationary large-amplitude fronts and pulses. Instead, the PDE exhibits attracting patterns that vary periodically in time with a small amplitude, while their spatial variation may either be uniformly small, or consist of distinct small- and large-amplitude parts (see also Fig. <ref> below).

To gain an intuitive understanding of this threshold, we recall from Section <ref> that, near the (spatial) Turing bifurcation at a=a_T, the trivial state (a,f(a)) is already unstable with respect to the long (spatial) wavelength perturbations associated to the (temporal) Hopf bifurcation at a=1, since a_T<1 by (<ref>) –see Fig. <ref>. 
Hence, if the initial data consist of perturbations of (a,f(a)) for which the support is sufficiently large, then one expects that the unstable Hopf modes will be triggered.
On the other hand, if the initial data consists of sufficiently localized perturbations, i.e., has  plateaus that are too narrow, then the dynamics of the PDE are predominantly driven by the Turing mode, as seen above. 
Moreover, these observations also suggest that to understand the dynamics of patterns in (<ref>) that do not originate from such localized perturbations, one needs to consider both instability mechanisms, i.e. one needs to study the co-dimension 2 bifurcation in which the spatial Turing mode and the temporal Hopf mode interact. 

In that case, the natural Ginzburg-Landau set-up involves two (complex) amplitudes: A(ξ,τ) that modulates the (linear) Turing mode e^i k_T x and B(ξ,τ) that modulates the Hopf mode e^i ω_T t (cf. Section <ref>). In the generic setting, the dynamics of small-amplitude solutions near this bifurcation has been studied in <cit.> (in the setting of reaction-diffusion systems).
In fact, the validity of a coupled (cubic) system of equations for A and B has been established in <cit.> for the generic situation in which
the shapes of the eigenvalue curves λ_±(k) are independent of the small parameter δ.  
Here, significant additional analysis appears to be required to go beyond the generic situation, because in (<ref>) the distance between the Hopf and Turing bifurcations is 2 √( d) = 2 √()δ. Thus, it is natural to take  δ > 0 as the main small parameter in the Ginzburg-Landau analysis (recall also that δ = 0.01 (and = 0.1) in most  simulations here). 
Furthermore, in order to derive the (A,B)-system that governs the dynamics of small patterns 𝒪(δ) close to a_T, one needs to adapt the generic scalings of <cit.>. We also leave this analysis for future work. 

Finally, regarding the PDE dynamics, there is also the question whether (some of) the stationary small-amplitude canard patterns presented in Sections <ref>, <ref>, and <ref> may also be stable as solutions of the PDE (<ref>). 
In our (limited) numerical simulations to date, we have not yet encountered such patterns. As observed above, by choosing the extent of the support of the initial conditions sufficiently large, i.e., above threshold, we see that system (<ref>) may exhibit small-amplitude patterns. 
In the simulation shown in the top row of Fig. <ref>, we consider localized initial data for which the spatial extent is somewhat larger than in Figs. <ref>, <ref>, and <ref> and take (again)  a=1. 
We observe the formation of a small-amplitude pattern that varies periodically in time and in space. Its spatial variation has a long wavelength character, which indicates that the amplitude A(ξ,τ) of the Turing mode has vanished. 
Thus, the PDE dynamics seem to be dominated by the Hopf bifurcation associated to the van der Pol equation and the associated (complex, uncoupled) Ginzburg-Landau equation for B(ξ, τ). 

This changes drastically if we decrease a to 0.999 (see the bottom row of Fig. <ref>). 
Here, as also in the simulations shown in Figs. <ref> (second column), <ref>, and <ref>, a localized, four-front, large-amplitude spatial structure appears near x=0 (see Fig. <ref>). 
However, unlike in Figs. <ref> and <ref>, the background state is now unstable with respect to the Hopf bifurcation, and the tails of the large-amplitude pattern vary periodically in time – and on the long spatial scale. (In fact, the localized structure itself also oscillates in time, but with an amplitude that is much smaller than that of the oscillations beyond the large-amplitude part of the pattern.) 
Note that the attractor of Fig. <ref> has the (expected) nature of a stationary localized structure as in Fig. <ref> that is destabilized by a Hopf bifurcation stemming from its essential spectrum (cf. <cit.>).

Of course, this evidence is far from sufficient to conclude that the stationary small-amplitude canard patterns cannot be stable.
So far, our findings only indicate that in the study of potentially stable small-amplitude patterns, one necessarily must also take the spatial aspects of the Hopf bifurcation into account. 
In other words, one needs to merge the (non-generic) coupled Ginzburg-Landau dynamics of the Turing mode A(ξ,τ) e^i k_T x and the Hopf mode B(ξ,τ) e^i ω_H t with the canard analysis in the present work. Moreover, a much more detailed and extensive (numerical) investigation of the PDE dynamics is crucial: our limited experiments show that the nature of the final attractor depends in a subtle way on small variations in the initial conditions and (consequentially) on the accuracy of the numerical procedure. 
These questions would be interesting and relevant subjects for future work. 






The types of stationary, multi-front patterns shown in Fig. <ref> have been constructed for various problems in the literature on singularly perturbed reaction-diffusion equations (see for example <cit.> and references therein). 
However, in the known constructions, the asymptotic homogeneous states of the homoclinic patterns correspond to critical points on normally hyperbolic manifolds, which is of course not the case here.






§ CONCLUSIONS AND OUTLOOK 




 §.§ Summary 



In this article, we reported on the discovery of classes of spatially periodic canard solutions that emerge from Turing bifurcations in the van der Pol PDE (<ref>) in one dimension, a phenomenon that we have dubbed “Turing's canards".
The canards that we studied analytically and numerically include classes of small-amplitude and large-amplitude spatially periodic canard solutions, with large, 𝒪(1), and small wavenumbers. (See Figs. <ref>–<ref>, as well as  Figs. <ref> and <ref> for representative canards, and see  Figs. <ref>,  <ref>, and <ref>–<ref> for bifurcation diagrams.) 
Furthermore, we observed numerically that several of these classes of spatially-periodic canards are attractors in the PDE.

The spatial ODE system, recall (<ref>), governs time-independent solutions of the PDE.
It has reversible, 1:1 resonant Hopf bifurcation points exactly at the parameter values where the PDE undergoes Turing bifurcations.
Quartets of eigenvalues merge there into two identical purely imaginary pairs, and hence hyperbolicity of the equilibrium/homogeneous state is lost,
recall Fig. <ref> and Proposition <ref> in Sec. <ref>. 
We performed a complete analysis of the Turing/reversible 1:1 Hopf point at a_T=√(1-2δ√()), recall (<ref>), in the full four-dimensional phase space (and the results for the other RFSN-II points follow by symmetry).
We studied the two-dimensional fast/layer system (Sec. <ref>) and the two-dimensional slow system (Sec. <ref>).
Both are one-degree-of-freedom Hamiltonian systems, due to the reversibility symmetry of the full spatial ODE system.  
Our analysis showed that the critical manifolds, which govern the slow dynamics to leading order, are two-dimensional, cubic-shaped manifolds consisting of saddle points of the fast system on the left and right branches and of  center points on the middle branch. 

We identified the key folded singularities, namely the reversible folded saddle-node points of type II (RFSN-II points), that lie on the fold sets between the saddle and center branches of the cubic crirical manifold. 
These RFSN-II exist asymptotically close to the Turing bifurcations. 
Using the method of geometric desingularization (see Section <ref>), we showed that the true and faux canards of the RFSN-II points are responsible for creating the spatially periodic canard patterns that we discovered, with the spatial canards having long segments near the true and faux canards. 

The analysis of the dynamics in the coordinate charts led to the discovery that there is a special algebraic solution Γ_0 (see Section <ref>) in the rescaling chart that constitutes the core component of the maximal spatial canards in the full system. 
The orbit Γ_0 consists of two branches, one corresponding to each of the true and faux canards of the RFSN-II point, that approach the cusp point from above and below. 
It is the geometrically unique orbit in the rescaling chart that asymptotes to the two critical points in the extry-exit chart on the equator of the hemisphere, thereby serving as a separatrix (or “river" type solution) over the hemisphere (in analogy to the parabola in the rescaling chart of the canard explosion in the fast-slow van der Pol ODE, recall <cit.>). 
Perturbation analysis of this algebraic solution then led to the calculation of a critical value  a_c(δ) 
(recall (<ref>)) at which this orbit persists to leading order in δ (where the asympotic expansion was calculated using the dynamic, small-amplitude, perturbation parameter r_2 in the rescaling chart). 
This is the value at which the true and faux canards continue into each other.
Further analysis of key solutions focused on their smoothness in δ, and critical values of a were identified, the next one of which is a_c2(δ)  (recall (<ref>)). 
These calculations were performed on a fixed level set of the conserved quantity 𝒢̃, and can be generalized.  
All of these canards are maximal canards that have the longest segments near the true and faux canards of the RFSN-II points, and they serve as boundaries in phase and parameter space separating spatially periodic canards of different profiles.


The dynamics of the canards change as one moves along the isolas of periodic solutions in the (a,k) parameter plane. 
One elementary change along branches of isolas occurs when the length of the canard segments grows as a changes (recall the green branch in Fig. <ref>). 
Next, for branches of isolas of spatially periodic solutions with canard segments along the true and faux canards of the RFSN-II point on the right slow manifold S_s^+, we studied the nucleation of interior spikes at the RFSN-II point on the left critical manifold S_s^- (recall the orange branch in Fig. <ref>). 
In addition, we followed these new spikes through the transition in parameter space from small to large spikes (recall  Section <ref>).
A number of further bifurcations of spatially periodic canards were also observed, including a bifurcation to canards with double loops (recall Fig. <ref>).

Self-similarity plays a central role in the spatially periodic canards. 
We showed that, in the singular limit δ=0, the zero level set of the Hamiltonian of the slow system is scale invariant, and hence that it has an infinite self-similarity. 
This self-similarity manifests through a sequence of crossing points and small, nested loops. 
Then, we observed that for 0<δ≪ 1, the self-similarity is broken. 
The spatial patterns exhibit nearly self-similar dynamics, recall Fig. <ref>.
Solutions with different numbers of nested, successively-smaller, nearly self-similar loops are shown in Figs. <ref> and <ref>. 
Moreover, the wavenumber decreases as the number of successive loops increases.

The spatial canards are analogs in spatial dynamics of the situation in fast-slow ODEs with time-periodic limit cycle canards, where the explosion of temporal canards occurs near –and asymptotically close to– the singular Hopf bifurcation. 
The new spatially periodic canards were found  (recall Sec. <ref>) to have many features in common with the classical temporal limit cycle canards, as well as several  important new features. The most interesting of these new features is that the critical manifold in the spatial ODE system is two-dimensional and consists of branches of saddle equilibria and center equilibria of the fast system, which contrasts with the one-dimensional critical manifolds of attracting and repelling equilibria that give rise to canard explosions in fast-slow ODEs.

To complement the above results for the spatial ODE system, we also studied some basic dynamical properties of the PDE (<ref>), recall Section <ref>.  
We showed that the Turing bifurcation from which the spatially periodic canards is sub-critical for the PDE.
The standard Ginzburg-Landau theory shows that small-amplitude perturbations should grow because the coefficients of both the linear and cubic terms are positive, and it cannot determine what the nonlinear saturation mechanisms might be. 
Hence, this study also sheds new light on what attractors can exist in the sub-critical case.

The PDE simulations that we have performed so far showed that the large-amplitude, spatially periodic canards are attractors in the full PDE (<ref>). With localized initial data, which are at the homogeneous state over most of the interval and have localized tanh-shaped perturbations (`plateaus'), we observed that, when the plateaus are not too wide, the attractors are stationary, large-amplitude canard patterns that are homoclinic in space to the (stable) homogeneous state (a,f(a)) for a ≳ a_T.
Then, for a ≲ a_T, with the same initial data (plateaus not too wide), the attractors are large-amplitude canard patterns homoclinic to spatially periodic patterns, since the homogeneous state is linearly unstable here. Recall Figs. <ref>-<ref>. 
By contrast, when the plateaus in the initial data are too wide, then the attractors exhibit time periodic dynamics for a<1, due to the Hopf bifurcation.
In addition, a class of small-amplitude attractors is observed in the PDE. 
For 1 > a ≳ a_T, these are periodic in both time and space, and they are observed for initial data in which the support is sufficiently wide, greater than a threshold. 
By contrast, for a ≲  a_T, the same iniital data develops into a large-amplitude attractor whose spatial tails vary periodically in time.   
Recall Fig. <ref>. 
Furthermore, the simuations suggest that there is a rich interplay between the Turing and Hopf modes, especially since the real parts of the dominant eigenvalue are both of size 𝒪(δ) at the Hopf value a=1 and the Turing bifurcation a_T=√(1-2δ√()) (recall Fig. <ref>), and it is expected that the PDE dynamics are governed by coupled Ginzburg-Landau equations (one for each mode), as described in Section <ref>.




 §.§ RFSN-II points and spatial canards in general reaction-diffusion systems 



In this brief section, we generalize some of the main results about the classes of small-amplitude, spatially periodic canards asymptotically close to the Turing bifurcation value a_T established here for the van der Pol PDE (<ref>) to other reaction-diffusion systems with separated diffusivities, which undergo Turing bifurcations and which have reaction kinetics consisting of two or more branches separated by non-degenerate fold points.
In particular, we consider 
systems of reaction-diffusion PDEs


    u_t    = f̂(u,v) + d u_xx, 
     
        v_t    = ĝ(u,v) + v_xx,


where 0<d≪ 1.
To generalize the results about the RFSN-II points and their canards, system (<ref>) needs to have a sufficiently large open set in the (u,v) plane on which the critical set f̂(u,v)=0 has a locally unique (non-degenerate) quadratic solution v=h_0(u)  i.e., a point (u_0,v_0) in the open set such that f̂(u_0,h_0(u_0))=0, ∂f̂/∂ u (u_0,h_0(u_0))=0, and  ∂^2 f̂/∂ u^2(u_0,h_0(u_0))  0.   
With this assumption, there will be a critical manifold consisting of a saddle branch and a center branch that meet along a fold curve in the (u,p,v,q) space
(recall (<ref>) and Fig. <ref> for the van der Pol system).


Like all reaction-diffusion systems in one dimension, these multi-scale reaction-diffusion systems are invariant under the transformation x → -x. 
Hence, the systems of spatial ODEs governing the shapes of time-independent patterns have a reversibility symmetry of the form ℛ, recall (<ref>), and the Turing bifurcation points in these systems correspond to reversible 1:1 resonant Hopf bifurcations in the spatial ODE systems. 
In addition, the fast system (or layer problem) 


    u_y    = p 
    
        p_y    = -f̂(u,v)


will be Hamiltonian, with v as a parameter, recall (<ref>).
Furthermore, the desingularized reduced vector fields on the critical manifold S={  p = 0, v = h_0(u) } will be


    u_x_d   = q, 
     
        q_x_d   = f̂_u/f̂_vĝ(u,h_0(u)).


These systems are also Hamiltonian, just as the desingularized reduced system (<ref>) is for the van der Pol system.  
Due to the reversibility symmetry, the ordinary singularities (fixed points) of the desingularized vector fields must be saddles, centers, or saddle-nodes, and the folded singularities must be reversible folded saddle (RFS), reversible folded center, or reversible folded saddle-node (RFSN) points, since the eigenvalues of the Jacobians must be symmetric with respect to the real and imaginary axes.
Finally, the true and faux canards of RFSN-II and RFS points will give rise to the canard segments of the spatially periodic solutions.




 §.§ Discussion 



Open questions about the full PDE (<ref>) have been listed at the end of Section <ref>. 
Here, we list some open questions about the spatial ODE system.

We are presently performing the geometric desingularization of the reversible folded saddle (RFS) points on the fold sets L^± in (<ref>), which were shown to exist for parameter values further from the Turing bifurcation (and hence from the RFSN-II singularity) in Section <ref>.  
Simulations (see for example  Figs. <ref>–<ref>) show that these RFS singularities are responsible for the creation of the spatially periodic canard solutions observed for a values sufficiently far away from a_T. The further a is from a_T, the true and faux canards of the RFS play a similar role in the geometric deconstructions of the spatially periodic solutions as those of the RFSN-II for a near a_T. Also, the ordinary singularity E, which is a center for a<1, is located further from the fold set, and hence it would be of interest to determine how the number and structure of the loops changes as a decreases.
Furthermore, in the opposite limit as a → 1, preliminary calculations indicate  that the canards of the RFS points converge to the canards created by the RFSN-II.


The results here motivate a more general rigorous analysis of folded singularities in singularly perturbed systems of ODEs with reversibility symmetry, especially into the geometry of the invariant manifolds associated to RFSN-II and RFS points. 
We refer to <cit.> for the theory of general (not necessarily reversible) folded saddle nodes, and to <cit.> for the theory of general folded saddles.  
 
The three-component reaction-diffusion model studied in <cit.> consists of one activator and two inhibitors. 
The kinetics of the activator and first inhibitor are essentially those of the FitzHugh-Nagumo ODE. 
The kinetics of the second inhibitor are also linear.
All three species diffuse, with the diffusivity of the activator being asymptotically smaller than those of the inhibitors.
Overall, the reaction-diffusion subsystem for the activator and first inhibitor is of the general form (<ref>).
Hence, it would be of interest to study the roles of the folded singularities and their canards in the onset of the spikes in the two-component FHN system in which the inhibitor also diffuses, as well as in the full three-component reaction-diffusion system.
The formation of the spike on the left branch, from S_s^- due to the nearby RFSN-II point and its canards observed here for <ref> (recall Sec. <ref>) may explain the nucleation of some spikes in the three-component model of <cit.>.
In the other direction, analysis similar to that in <cit.> of how a small spike grows into a full-fledged spike should apply here for system <ref>. 

Another question concerns the double asymptotic limit in which  and δ are small. It is well-known that the asymptotic limit of small  gives rise to temporal canards in the kinetics, and we have established here that the asymptotic limit of small δ gives rise to canards in the spatial dynamics. 
How do the temporal limit cycle canards of the kinetics problem for small  interact with spatial canards?

Acknowledgments.
The authors gratefully acknowledge Irv Epstein, Guido Schneider, and Gene Wayne,
for useful comments. We also thank Irv Epstein for bringing reference <cit.> to our attention, for the possible evidence of spatial canards. 
The results in this article were presented by T.V. at the conference Multiscale Systems: Theory and Applications, held July 8-12, 2024, at the Lorentz Center, Leiden University, Leiden, NL.
NSF-DMS 1616064 and NSF-DMS 1853342 provided partial support to T.K. and T.V., respectively. 
The research of A.D. is supported by the ERC-Synergy project RESILIENCE (101071417).







§ NUMERICAL METHODS







 §.§ Continuation of periodic solutions
 


Periodic solutions of the system (<ref>) with non-trivial first integral 𝒢 given by (<ref>) were computed and numerically continued using the method developed in <cit.>. More specifically, families of periodic solutions were computed by calculating solutions of the auxiliary system 


    u̇   = T (F( u) + η∇𝒢),


subject to periodic boundary conditions u(0) =  u(1). Here,
u = (u,p,v,q), F is the vector field in (<ref>), T is the spatial period, η is a new auxiliary parameter (fixed at zero), and the overdot denotes the derivative with respect to the scaled spatial variable y = yT.

The branches of 𝒢( u) = g, where g is a constant, were computed by appending the integral constraint


    ∫_0^1 𝒢( u)   dy = g,


to the system (<ref>) subject to the periodic boundary conditions u(0) =  u(1).

The numerical continuation was implemented using the continuation software AUTO <cit.>. 
The bifurcation diagrams for the families of spatially periodic canard solutions of (<ref>) shown in Figs. <ref>,   <ref>, and later figures were also computed using AUTO.  We note that  additional branches of periodic solutions may exist.




 §.§ Computation of saddle slow manifolds and maximal canards
  

The saddle slow manifolds shown in Fig. <ref> were computed following the method developed in <cit.>. More specifically, the saddle slow manifolds were computed in two parts: one for the solutions enclosed by the true and faux canards of the folded saddle and one for the solutions outside that lie outside the region enclosed by the true and faux canards of the folded saddle. 

The subsets of S_s,δ^+ outside the region enclosed by the true and faux canards of the folded saddle were computed by solving the system (<ref>), where u = (u,p,v,q) and F is the vector field in (<ref>), and subject to the boundary conditions


    u(0) ∈{ v = 13u^3-u,    q = q_0 : q_0 < 0 }   and    u(1) ∈{ p = 0,    u = 1 }.


The left-end condition enforces the constraint that solutions enter the saddle slow manifold along the p-nullcline at a fixed q-distance from the folded singularity. The right-end condition is a statement that solutions leave the neighbourhood of the slow manifold along the u-nullcline and terminate at the fold. 

The subsets of S_s,δ^+ enclosed by the true and faux canards of the folded saddle were computed by solving the system (<ref>), where u = (u,p,v,q) and F is the vector field in (<ref>), and subject to the boundary conditions


    u(0) ∈{ v = 13u^3-u,    u = u_0 : u_0 > 1 },
       u(1) ∈{ v = 13u^3-u }, 
       and   { q(0) + q(1) = 0 }.


The left-end condition specifies that solutions enter the saddle slow manifold along the p-nullcline at a fixed u-distance from the fold. The remaining boundary conditions ensure that the solutions turn away from the fold and stay on the slow manifold. 

The maximal canards were then computed as saddle-node bifurcations of the above two-point boundary value problems. 





 §.§ Continuation of canard orbits in chart Lg
  

For the results shown in Section <ref>, solutions of the blown-up system (<ref>) on the zero level set H_2(u_2,p_2,v_2,q_2,r_2) = 0 were computed and numerically continued by solving the auxiliary problem (<ref>), where u = (u_2,p_2,v_2,q_2), F is the vector field in (<ref>), and the conserved quantity is 𝒢 = H_2. These equations were solved subject to the boundary conditions 


    p_2(0) + p_2(1) = 0,    and    q_2(0) + q_2(1) = 0,


which enforces the symmetry (u_2,p_2,v_2,q_2) → (u_2,-p_2,v_2,-q_2), together with the integral constraint


    ∫_0^1 H_2(u_2,p_2,v_2,q_2,r_2) dy = 0


which constrains the solution to the zero level contour of the Hamiltonian. 




 §.§ ODE and PDE Simulations
 

Finally, direct numerical simulations of solutions of the spatial dynamics of (<ref>) were carried out using Mathematica's in-built ODE solvers.
Direct numerical solutions of the full PDE (<ref>) were performed using
Mathematica's NDSolveValue-package with the “MethodOfLines” and “SpatialDiscretization” prescibed by {“TensorProductGrid”, “MinPoints” → 10000}.




§ THE PROOF OF PROPOSITION  LG 



In this appendix, we prove Proposition <ref>. by applying Theorem 3.21 from Chapter 4.3.3  of <cit.>.  
We work with the spatial ODE system in a general form, 

    0    = v - f(u) + d_u u_xx
    
            0    = (a-u) + d_v v_xx,

where the diffusivities, d_u and d_v are positive. 

To derive the normal form, we rewrite these spatial ODEs as the following fourth-order system:

    u_x    = p 
    
            p_x    = 1/d_u 
            (f(u)-v) 
    
            v_x    = q 
    
            q_x    = /d_v (u-a).

We translate the variables so that the equilibrium (a,0,f(a),0) is at the origin and use the notation u=(u_1,u_2,u_3,u_4) of Chapter 4.3.3  for the dependent variables: 
u=a+u_1, p = u_2, v = f(a) + u_3, and q=u_4. 
Hence, the system is

    u_1_x    = u_2 
    u_2_x    = 1/d_u[ (a^2-1) u_1 + a u_1^2 + 1/3 u_1^3 - u_3  ]
    u_3_x    = u_4 
    u_4_x    = /d_v u_1.

The quartet of eigenvalues is given by

    ±1/√(2)d_u d_v√( (a^2-1)d_u d_v^2 ± i d_v^3/2√( 4  d_u - (a^2-1)^2 d_v)).

At a=a_T=√(1 - 2 √( d_u/d_v)), the quartet degenerates into two coincident pairs of pure imaginary eigenvalues 

    ± i ω = ± i (/d_u d_v)^1/4.

The eigenvalues ± i ω are algebraically double and geometrically simple.
This is the non-degenerate, reversible, 1:1 resonant Hopf bifurcation identified in Section <ref>, recall also Fig. <ref> (where we note that d_u=d=δ^2 and d_v=1 in (<ref>)).

In order to unfold this point, we set 

    a = a_T + μ,

where the parameter μ here is different from the spatial eigenvalue μ used in Section <ref>,
and we write the system (<ref>) as

    u_x = ℱ (  u, μ) =  L u 
          + R_20( u, u) 
        +  R_30( u, u, u) 
        + μ R_11 ( u)
        + μ R_21( u, u)
        + μ^2  R_12( u),

where the operators are defined as

    L = [
        [      0      1      0      0; -2 ω^2      0 -1/d_u      0;      0      0      0      1;   /d_v      0      0      0 ]],    R_20( u, v) = [
        [               0; a_T/d_u u_1 v_1;               0;               0 ]],    R_30( u, v,  w) = [
        [                  0; 1/3d_u u_1 v_1 w_1;                  0;                  0 ]],


    R_11( u) = [
        [            0; 2a_T/d_u u_1;            0;            0 ]],    R_21( u, v) = [
        [             0; 1/d_u u_1 v_1;             0;             0 ]],    R_12( u) = [
        [         0; 1/d_u u_1;         0;         0 ]].

The terms represent, respectively  L: the Jacobian at the origin at a_T; 
R_20: 
the quadratic terms in u;
R_30: 
the cubic terms in u;  R_11 and R_21: the terms  that are linearly proportional to the unfolding parameter μ, and R_12: the term proportional to μ^2. (We  follow the general notation in <cit.> for the unfolding of this bifurcation, labeled as (iω)^2 in Chapter 4.3.3. The subscripts on R indicate the powers of u and μ, respectively, in each term in (<ref>); for example, the term involving R_21 is quadratic in u and linear in μ.)

For L, the Jacobian at a_T, we use the following eigenvector and generalized eigenvector associated to iω: 

    ζ_0 = [
        [         1;       i ω; -/ω^2 d_v; -i /ω d_v ]], 
          ζ_1 = [
        [       -i/ω;          2; -i/ω^3 d_v;          0 ]].

These satisfy ( L-iω)ζ_0=0 and ( L-iω)ζ_1 = ζ_0. 
Moreover, since L is real, one also has 
( L + i ω) ζ̅_0 = 0
and ( L + i ω) ζ̅_1 = ζ̅_0,
where the overbar denotes the complex conjugate.
The set {ζ_0, ζ_1, ζ̅_0,ζ̅_1 } is used as a basis for ℝ^4, and we represent the vector u = A ζ_0 + B ζ_1 + A̅ζ̅_0 + B̅ζ̅_1 by u=(A,B,A̅,B̅), where A, B ∈ℂ.

Now, we derive the normal form of (<ref>). 
It may be obtained directly by applying Lemma 3.17 in Chapter 4.3.3 of <cit.>, as follows. 
The linear part, L, is conjugate to the block Jordan matrix,

    J = [
        [  i ω    1    0    0;    0  i ω    0    0;    0    0 -i ω    1;    0    0    0 -i ω ]].

For the nonlinear part, which also includes the μ-dependent terms, we observe that Hypotheses 3.1, 3.2, and 3.14 from <cit.> are satisfied
for every integer k ≥ 3, i.e., the vector is C^k for every K≥ 3 since it is polynomial.  Hence, Lemma 3.17 establishes that, for any integer p with 2 ≤ p ≤ k,  there exist neighborhoods 𝒱_1 of the origin in ℝ^4 and 𝒱_2 of the origin in ℝ and a real-valued polynomial Φ of degree p such that the near-identity coordinate change,
u = A ζ_0 + B ζ_1 + A̅ζ̅_0 + B̅ζ̅_1 + Φ(A,B,A̅,B̅,μ) defined on 𝒱_1 ×𝒱_2, transforms (<ref>) into the following normal form:

    A_x    = i ω A + B + i A P(μ, A A̅, i/2(A B̅ - A̅ B) ) + ρ_A 
    
            B_x    = i ω B 
            + i B P(μ,AA̅, I/2(AB̅-A̅B)) + A Q(μ, A A̅, i/2(A B̅-A̅ B)) + ρ_B.

Moreover, Φ(0,0,0,0,0)=0, ∂_(A,B,A̅,B̅)Φ(0,0,0,0,0)=0, the coefficients of the monomials of degree q in Φ(·,μ) are C^k-q in μ, and 
Φ satisfies the reversibility symmetry
ℛΦ(A,B,A̅,B̅,μ) = Φ(A̅,-B̅, A, -B, μ).
Here, P and Q are real-valued polynomials of degree p-1.
The remainder terms ρ_A(A,B,A̅,B̅,μ) and ρ_B(A,B,A̅,B̅,μ) are C^k functions, satisfy the estimate 
|ρ_A| + |ρ_B | = o( (|A| + |B|)^p), and have the following symmetries:
ℛρ_A(A,B,A̅,B̅,μ) = -ρ̅_A(A̅,B̅, A, B, μ) and
ℛρ_B(A,B,A̅,B̅,μ) = ρ̅_B(A̅,B̅, A, B, μ) for each μ.
(See the general normal form for reversible 1:1 resonant Hopf bifurcations given by equations (3.25) in Chapter 4.3.3 of  <cit.>.)

As shown in <cit.>, the existence of equilibria, spatially periodic orbits, quasi-periodic orbits, and homoclinic orbits of (<ref>) may be determined by working with the principal parts of the polynomials P and Q. That is, it suffices to work with p=2. 
Hence, throughout the rest of this appendix, we set

    P= α̂μ + β̂ A A̅ + i/2γ̂ (A B̅ - A̅ B)      and     Q = âμ + b̂ A A̅ + i/2ĉ(AB̅-A̅B).



After lengthy calculations using the invariance equations, one finds:

    â   = a_T/2 d_u,  b̂=25-8√(d_v/ d_u)/36d_u,   ĉ= 109 d_u (d_v/)^1/4 - 32(d_v/)^3/4/72 d_u^5/4, 
    α̂   = -a_T/4 ω d_u,   β̂= -1/16( d_v/ d_u^3)^1/4,  γ̂=  -32 d_v +37√( d_u d_v)/216  d_u.

(Note that these six parameters are the same as those in Hypothesis 3.18 in <cit.>, except that we have added the hats.)

Next, we turn to the signs of these six parameters in the principal parts of P and Q. Since a_T>0, we know â>0 and α̂<0.
Also, β̂<0. 
The sign of b̂ depends on the ratios of the diffusivities:

    b̂ < 0    for   d_u/d_v < 64/625,
               and  b̂ > 0    for   d_u/d_v > 64/625.

This sign analysis of b̂ determines the parameter conditions in Proposition <ref>. The former is referred to as the focusing/sub-critical case, and the latter as the defocusing/super-critical case.
The sign of ĉ is given by:

    ĉ > 0    for  d_u >  32/109√(d_v/),    and  ĉ < 0    for  d_u < 32/109√(d_v/).

Finally, the sign of γ̂ is given by

    γ̂ > 0    for  d_u >  (32/37)^2d_v/,    and  γ̂ < 0    for  d_u < (32/37)^2 d_v/.





§ THE PROOF OF PROPOSITION <REF>
  


In this appendix, we prove Proposition <ref>.

We analyze the dynamics around the nilpotent equilibrium (0,0,0) by performing the blow-up transformation


    u_2 = r u_2,    p_2 = r p_2,    and    q_2 = r q_2,


with (u_2,p_2,q_2)∈𝕊^2 and r≥ 0, which inflates the nilpotent equilibrium to the unit sphere. For our analysis, we restrict attention to the half-space {u_2 ≥ 0 } which completely contains the algebraic solutions Γ_0^±. We will examine the dynamics in three coordinate charts: 


    K_21: {u_2=1 },    K_22: {p_2=1 },    and    K_23: {q_2 = 1},


and then appeal to the symmetry (<ref>) to obtain the dynamics on the remainder of the blown-up hemisphere. We denote an object Φ of (<ref>) in the blown-up coordinates by Φ and in the charts K_2j by Φ_2j for j=1,2,3. 

Dynamics in chart K_21: {u_2=1 }: 
The blow-up transformation in chart K_21 is 


    u_2 = r_21,    p_2 = r_21  p_21,    and    q_2 = r_21 q_21.


Transformation and desingularization (dζ_21 = r_21 dη_2) gives


    ṙ_21   = √() r_21 p_21
    ṗ_21   = 23r_21-12√()( p_21^2+q_21^2 ) 
    q̇_21   = 1-√() p_21q_21


where the overdot denotes the derivative with respect to the rescaled spatial coordinate ζ_21. 

The set { r_21 = 0 } is invariant. The dynamics in the { r_21 = 0 } subspace are given by 


    ṗ_21   = -12√()( p_21^2+q_21^2 ) 
    q̇_21   = 1-√() p_21q_21.


The coordinate change (z_21,w_21) = 12( p_21-q_21, p_21+q_21) transforms the system (<ref>) to the decoupled system of Riccati equations


    ż_21   = -√() z_21^2-12
    ẇ_21   = -√() w_21^2+12.


From this, we find that solutions of (<ref>) are given by 


    p_21(ζ_21)    = 1√(2)^1/4[ tanh( ^1/4√(2)ζ_21 + tanh^-1( √(2)^1/4w̃_21) ) - tan( ^1/4√(2)ζ_21 - tan^-1( √(2)^1/4z̃_21) ) ] 
     
        q_21(ζ_21)    = 1√(2)^1/4[ tanh( ^1/4√(2)ζ_21 + tanh^-1( √(2)^1/4w̃_21) ) + tan( ^1/4√(2)ζ_21 - tan^-1( √(2)^1/4z̃_21) ) ]


where z̃_21 and w̃_21 are the values of z_21 and w_21 at ζ_21=0. Among these solutions, we distinguish two particular solutions:


    ℓ_21^± := { p_21+q_21 = ±√(2)^-1/4}.


The solutions with initial conditions above ℓ_21^- are forward asymptotic (i.e., ζ_21→ +∞) to the line ℓ_21^+





and the solutions with initial conditions below ℓ_21^+ are backward asymptotic (i.e., ζ_21→ -∞) to the line ℓ_21^-. 




The dynamics in the invariant subspace {r_21=0} are shown in Fig. <ref>(a).



Dynamics in chart K_22: {p_2 = 1 }: The blow-up transformation in chart K_22 is


    u_2 = r_22  u_22,    p_2 = r_22,    and    q_2 = r_22  q_22.


Transformation and desingularization (dζ_22 = r_22 dη_2) gives 


    ṙ_22   = r_22( 12√()(1-q_22^2)+23 r_22 u_22^3 ) 
    u̇_22   = 12√()u_22(1+q_22^2)-23r_22 u_22^4 
    q̇_22   = -12√()q_22(1-q_22^2)+u_22^2-23r_22u_22^3 q_22,


where the overdot has been recycled to denote the derivative with respect to ζ_22. 

The line { r_22 = 0, u_22 = 0 } is invariant with dynamics governed by 


    q̇_22 = -12√() q_22(1-q_22^2).


The equilibrium at q_22=0 is stable, and the equilibria at q_22 = ± 1 are unstable. 

The plane { u_22 = 0 } is also invariant. The dynamics restricted to this plane are given by 


    ṙ_22   = 12√()r_22(1-q_22^2) 
    q̇_22   = -12√()q_22(1-q_22^2).


The system possesses a pair of lines of equilibria, ℒ_22,+^u = { q_22 = 1 } and ℒ_22,-^u = { q_22 = -1 }, which are both center-unstable with eigenvalues λ_u = √() and λ_c = 0. The associated eigenspaces are given by 


    𝔼^u( ℒ_22,±^u ) = [ ∓ r_22;      1 ]   and   𝔼^c( ℒ_22,±^u ) = [ 1; 0 ].


The lines ℒ_22,±^u correspond to the center-unstable lines of equilibria ℒ_±^u. 

The plane { r_22 = 0 } is invariant. The dynamics restricted to this subspace are given by 


    u̇_22   = 12√()u_22(1+q_22^2) 
    q̇_22   = -12√()q_22(1-q_22^2)+u_22^2.


As shown in Fig. <ref>(b), there is a saddle equilibrium at the origin, with stable eigenvalue λ_s = -12√() and stable eigendirection aligned with the q_22-axis, and unstable eigenvalue λ_u = 12√() and unstable eigendirection aligned with the u_22-axis. 
Moreover, there is a pair of unstable degenerate nodes at (u_22,q_22)=(0,± 1) with spectra given by σ_u = {√(), √()} and corresponding eigenspaces


    𝔼^u( 0,± 1 ) = span{[ 1; 0 ], [ 0; 1 ]}.



Thus, the equilibria (r_22,u_22,q_22)=(0,0,± 1) correspond to the intersections of the center-unstable lines, ℒ_±^u, of equilibria (see (<ref>) and (<ref>)) with the blown-up hemisphere. 

We concl;ude with the following two key observations. First, there are two straight line solutions, 


    ℓ_22^± := { r_22=0, u_22 = ±1√(2)^1/4( q_22+1 ) }.


Second, the unstable manifold, W^u(0,0)=:Γ_22^+, of the saddle equilibrium at the origin splits the (u_22,q_22) phase space into two regions. Solutions with initial conditions to the left of Γ_22^+ are backward asymptotic to the equilibrium at (r_22,u_22,q_22)=(0,0,-1), and solutions with initial conditions to the right of Γ_22^+ are backward asymptotic to the equilibrium at (r_22,u_22,q_22)=(0,0,1). The dynamics in the { r_22=0 } subspace are shown in Fig. <ref>(b).



Dynamics in chart K_23: {q_2 = 1 }: The blow-up transformation in chart K_23 is given by 


    u_2 = r_23 u_23,    p_2 = r_23 p_23,    and    q_2 = r_23.


Transformation and desingularization (dζ_23=r_23dη_2) gives


    ṙ_23   = r_23u_23^2 
    u̇_23   = u_23( √()p_23-u_23^2 ) 
    ṗ_23   = 12√()( p_23^2-1 )-u_23^2p_23+23 r_23 u_23^3,


where the overdot now denotes derivatives with respect to ζ_23. 

The line { r_23=0, u_23=0 } is invariant. The dynamics on this line are governed by 


    ṗ_23 = 12√()( p_23^2-1 ).


There is a stable equilibrium at p_23=-1 and an unstable equilibrium at p_23=1. 

The plane { u_23=0 } is invariant. The dynamics in the (r_23,p_23) plane are given by 


    ṙ_23   = 0 
    ṗ_23   = 12√()( p_23^2-1 ).


The line ℒ_23,-^s = { p_23=-1 } of equilibria is center-stable and corresponds to the line ℒ_-^s. 
The line ℒ_23,+^u = { p_23=1 } of equilibria is center-unstable and corresponds to the line ℒ_+^u. 

The plane { r_23=0 } is invariant. The dynamics on this subspace are governed by 


    u̇_23   = u_23( √()p_23-u_23^2 ) 
    ṗ_23   = 12√()( p_23^2-1 )-u_23^2p_23.


As shown in Fig. <ref>(c), there is a stable degenerate node at (u_23,p_23)=(0,-1) with spectrum σ_s={ -√(),-√()} and stable subspace 


    𝔼^s( 0,-1 ) = span{[ 1; 0 ], [ 0; 1 ]}.


There is also an unstable degenerate node at (u_23,p_23)=(0,1) with spectrum σ_u={√(),√()} and unstable subspace 


    𝔼^u( 0,1 ) = span{[ 1; 0 ], [ 0; 1 ]}.


Solutions converge to the equilibium at (u_23,p_23)=(0,-1) along the lines 


    ℓ_23^± := { u_23=±1√(2)^1/4(p_23+1) }.


The dynamics in the {r_23=0} subspace are shown in Fig. <ref>(c). 

The equilibrium (r_23,u_23,p_23)=(0,0,1) corresponds to the intersection of the center-unstable line, ℒ_+^u, of equilibria with the blown-up hemisphere. Similarly, the equilibrium (r_23,u_23,p_23)=(0,0,-1) corresponds to the intersection of the center-stable line, ℒ_-^s, of equilibria with the blown-up hemisphere.


Transition between charts K_21:{u_2=1 } and K_22:{p_2=1 }: The coordinate change, φ_12(r_21,p_21,q_21), from chart K_21 to chart K_22 is given by 


    φ_12(r_21,p_21,q_21) = (r_22, u_22, q_22) = ( r_21p_21, 1p_21, q_21p_21),    for  p_21 >0.


The inverse map, φ_21(r_22,u_22,q_22), which transports coordinates from chart K_22 to K_21, is given by 


    φ_21(r_22,u_22,q_22) = (r_21,p_21,q_21) = ( r_22u_22, 1u_22,q_22u_22),    for   u_22>0.


Under this coordinate change, we find that 


    lim_ζ_21→ -∞φ_12( ℓ_21^±) = lim_ζ_22→ -∞ℓ_22^± = { (r_22,u_22,q_22) = (0,0,-1) }.


Thus, solutions emanate from the point ℒ_-^u ∩𝕊^2 on the equator of the blown-up hemisphere. 
See Fig. <ref>.


Transition between K_21:{u_2=1 } and K_23:{q_2=1 }: The coordinate change, φ_13(r_21,p_21,q_21), from chart K_21 to chart K_23 is given by 


    φ_13(r_21,p_21,q_21) = (r_23,u_23,p_23) = ( r_21q_21, 1q_21, p_21q_21),    for   q_21 >0.


The inverse map, φ_31(r_23,u_23,p_23), which transports coordinates from chart K_23 to K_21, is given by


    φ_31(r_23,u_23,p_23) = (r_21,p_21,q_21) = ( r_23u_23, p_23u_23, 1u_23),    for   u_23 >0.


From these transition maps, we find that the image of the attracting line ℓ_21^+ is the line ℓ_23^+ and is forward asymptotic to the equilibrium (r_23,u_23,p_23)=(0,0,-1). 
That is, 


    lim_ζ_21→∞φ_13( ℓ_21^+ ) = lim_ζ_23→∞ℓ_23^+ = { (r_23,u_23,p_23)=(0,0,-1) }.


Thus, solutions terminate at the point ℒ_-^a ∩𝕊^2 on the blown-up hemisphere. 
See Fig. <ref>.


Transition between K_22:{p_2=1 } and K_23:{q_2=1 }:
The coordinate change, φ_23(r_22,u_22,q_22), from chart K_22 to chart K_23 is given by 


    φ_23(r_22,u_22,q_22) = (r_23,u_23,p_23) = ( r_22q_22, u_22q_22, 1q_22),    for   q_22 > 0.


The inverse map, φ_32(r_23,u_23,p_23), which transports coordinates from chart K_23 to K_22, is given by


    φ_32(r_23,u_23,p_23) = (r_22,u_22,q_22) = ( r_23p_23, u_23p_23, 1p_23),    for   p_23 >0.


The images of the lines ℓ_22^± are the lines ℓ_23^± and they are forward asymptotic to the equilibrium at (r_23,u_23,p_23)=(0,0,-1). That is,


    lim_ζ_22→∞φ_23( ℓ_22^±) = lim_ζ_23→∞ℓ_23^± = { (r_23,u_23,p_23)=(0,0,-1) }.


Thus, solutions terminate at the point ℒ_-^a ∩𝕊^2 on the equator of the blown-up hemisphere. 
See Fig. <ref>.


Symmetry: With the above analysis in hand, we can obtain the dynamics in the coordinate charts K_24:={p_2 = -1 } and K_25:={q_2 = -1 } via the symmetry (<ref>). In particular, the image of the unstable manifold, Γ_22^+, of the saddle equilibrium at the origin in chart K_22: {p_2=1} under the symmetry transformation (<ref>) is the stable manifold, W^s(0,0,0) =: Γ_24^-, of the saddle equilibrium at the origin in chart K_24: {p_2=-1}.


Dynamics on the hemisphere: From our analysis of the dynamics in the charts K_21, K_22, and K_23, together with the transition maps φ_ij between them, we conclude the following. 


  * Solutions emanate from the point ℒ_-^u ∩𝕊^2 on the equator of the blown-up hemisphere, travel over the top of the hemisphere in the region enclosed by Γ_0^- and Γ_0^+, and terminate at the point ℒ_-^s ∩𝕊^2 on the equator of the blown-up hemisphere. These are the class 1 heteroclinics.

  * Solutions emanate from the point ℒ_+^u ∩𝕊^2 on the equator of the blown-up hemisphere, travel over the top of the hemisphere in the region enclosed by Γ_0^+ and the equator, and terminate at the point ℒ_-^s ∩𝕊^2 on the equator of the blown-up hemisphere. These are the class 2 heteroclinics. 

  * The unstable manifold, Γ_0^+, of the equilibrium point corresponding to the intersection of the positive p_2-axis with 𝕊^2 is the separatrix that divides between class 1 and class 2 heteroclinics. 

  * Solutions emanate from the point ℒ_-^u ∩𝕊^2 on the equator of the blown-up hemisphere, travel over the top of the hemisphere in the region enclosed by Γ_0^- and the equator, and terminate at the point ℒ_+^s. These are the class 3 heteroclinics. 

  * The stable manifold, Γ_0^-, of the equilibrium point corresponding to the intersection of the negative p_2-axis with 𝕊^2 is the separatrix that divides between class 1 and class 3 heteroclinics.









§ ANALYSIS OF THE GOVERNING EQUATIONS IN THE ENTRY/EXIT CHART LG
 


In this appendix, we present the analysis of (<ref>) in the sequence of invariant hyperplanes:
{δ_1=0}∩{ a_1=0}, 
{ r_1=0}∩{ a_1 = 0 }, 
{ r_1 = 0 }∩{δ_1 = 0 }, 
{δ_1=0 }, 
and { r_1=0},

respectively.
These are the intermediate results used in Section <ref> to obtain the invariant sets and dynamics of the full system (<ref>).

In the invariant hyperplane
{δ_1=0}∩{ a_1=0 }, system (<ref>) reduces to

    ṙ_1    = 1/2√() p_1 r_1 
    ṗ_1    = 1 - v_1 - 3/2√() p_1^2
                                  + 1/3√()r_1^2 
    v̇_1    = -2 √() p_1 v_1 
    q̇_1    = -3/2√() p_1 q_1.

The first three components of this vector field are independent of q_1, and hence the fourth equation decouples from the others.
The line I given by (<ref>)
is also invariant for this larger system,
and E_± are again fixed points on I,
recall (<ref>).
The equilibrium E_+ is a saddle with three stable eigenvalues
-√(6)^1/4,
-2√(2/3)^1/4,
and -√(3/2)^1/4,
and one unstable eigenvalue
1/√(6)^1/4.
The stable subspace is given by

    span{[ [ 0; 1; 0; 0 ]],
      [ [        0;    -√(3); √(2)^1/4;        0 ]],
      [ [ 0; 0; 0; 1 ]]
    },

and the unstable subspace is in the r_1-direction,

    span{[ [ 1; 0; 0; 0 ]]
            }.


The equilibrium E_- is a saddle with three unstable eigenvalues
√(3/2)^1/4,
2√(2/3)^1/4,
and √(6)^1/4,
and one stable eigenvalue
-1/√(6)^1/4.
The unstable subspace is given by

    span{[ [ 0; 0; 0; 1 ]],
      [ [        0;     √(3); √(2)^1/4;        0 ]],
      [ [ 0; 1; 0; 0 ]]
    },

and the stable subspace is in the r_1-direction,

    span{[ [ 1; 0; 0; 0 ]]
            }.


There is also a two-dimensional surface of equilibria
in the hyperplane {δ_1=0 }∩{a_1=0}.
It is given by

    𝒮={ r_1 ∈ℝ,  δ_1=0,   p_1=0,  
                           v_1=1 + 1/3√()r_1^2,  
                           q_1 ∈ℝ,   a_1=0}.

The eigenvalues
of the Jacobian at points on 𝒮
are

    λ_s = -√(2√() + r_1^2 ),   λ_u = √(2√() + r_1^2 ),   λ_c = 0,0.

Hence, 𝒮
is a surface of saddle points,
and it corresponds to the saddle branches
S_s^± of the critical manifold S.
The associated stable, unstable, and center subspaces are

    𝔼_s = span{[ [              -√() r_1;     2√(2√() + r_1^2 ); 4/3√()(3 + r_1^2 √());             3 √() q_1 ]] },
      𝔼_u = span{[ [              -√() r_1;    -2√(2√() + r_1^2 ); 4/3√()(3 + r_1^2 √());             3 √() q_1 ]] },
      𝔼_c 
                 = {[ [ 0; 0; 0; 1 ]],
                 [ [         3;         0; 2 √() r_1;         0 ]] }.

The center manifold W^c(𝒮) is two-dimensional
in this hyperplane, and it emanates from ℓ.



Next, we examine the invariant hyperplane
{r_1=0}∩{a_1=0}.
Here, system (<ref>) reduces to

    δ̇_1    = -√() p_1 δ_1 
    ṗ_1    = 1 - v_1 - 3/2√() p_1^2 
    v̇_1    = √() (-2 p_1 v_1 + δ_1 q_1)
    q̇_1    = -3/2√() p_1 q_1 + δ_1.

This system is fully coupled,
in contrast to the above system. The line I given by (<ref>) is also invariant for this larger system, and E_± are again fixed points on I, recall (<ref>).
The equilibrium E_+ is stable with eigenvalues
-√(6)^1/4,
-2√(2/3)^1/4,
-√(3/2)^1/4,
and -√(2/3)^1/4.
The associated eigenvectors are

    [ [ 0; 1; 0; 0 ]],
      [ [        0;    -√(3); √(2)^1/4;        0 ]],
      [ [ 0; 0; 0; 1 ]],
      [ [ ^1/4;    0;    0; √(6) ]].


The equilibrium E_- is unstable with eigenvalues
√(2/3)^1/4,
√(3/2)^1/4,
2√(2/3)^1/4,
and √(6)^1/4, and eigenvectors 

    [ [ -^1/4;     0;     0;  √(6) ]],
      [ [ 0; 0; 0; 1 ]],
      [ [        0;     √(3); √(2)^1/4;        0 ]],
      [ [ 0; 1; 0; 0 ]].

Also, the line ℓ (recall (<ref>))
is a line of saddle equilibria
of (<ref>),
since the eigenvalues are

    λ_s = -√(2)^1/4,   λ_u = √(2)^1/4,   λ_c = 0,0.

The associated stable, unstable, and center subspaces are

    𝔼_s = span{[ [          0; 2√(2)^-1/4;          4;      3 q_1 ]] },
      𝔼_u = span{[ [           0; -2√(2)^-1/4;           4;       3 q_1 ]] },
      𝔼_c = span{[ [ 0; 0; 0; 1 ]],  [ [    4/4-3√()q_1^2; 2q_1/4-3√()q_1^2;                0;                0 ]] }.

Hence, (<ref>) has a two-dimensional center manifold

    N ≡ W^c(ℓ) ∩{ r_1=0 }.

Moreover, for δ_1>0, N is unique in the half-space { p_1 < 0 }.
Also, δ_1 increases on N
in this half-space.

Proceeding, we present the main results for the dynamics in the { r_1=0}∩{δ_1=0} hyperplane. The equations here are given by system (<ref>) with the equation ȧ_1= -3/2√() p_1 a_1 appended. Since the (p_1,v_1,q_1) subsystem is independent of a_1 on this hyperplane, the dynamics in these three variables are the same as the dynamics for (<ref>), and then one may solve for a_1 using quadrature.
The line I given by (<ref>) is invariant, and E_± are fixed points on it.
The equilibrium E_+ is stable with eigenvalues
-√(6)^1/4,
-2√(2/3)^1/4,
-√(3/2)^1/4, and -√(3/2)^1/4, 
and eigenvectors

    w_1^+ = [ [ 1; 0; 0; 0 ]],
       w_2^+ = [ [    -√(3); √(2)^1/4;        0;        0 ]],
       w_3^+ = [ [ 0; 0; 1; 0 ]],
                    w_4^+ = [ [ 0; 0; 0; 1 ]].

The equilibrium E_- is unstable with eigenvalues
√(3/2)^1/4,
√(3/2)^1/4, 
2√(2/3)^1/4, √(6)^1/4,
and eigenvectors

    w_1^- = [ [ 0; 0; 1; 0 ]],
       w_2^- = [ [ 0; 0; 0; 1 ]],
       w_3^- = [ [     √(3); √(2)^1/4;        0;        0 ]],
       w_4^- = [ [ 1; 0; 0; 0 ]].


There is also a two-dimensional surface of equilibria

    𝒮̂={ r_1=0,  δ_1=0,   p_1=0,  
                           v_1=1,  
                           q_1 ∈ℝ,   a_1 ∈ℝ}.

The eigenvalues of the Jacobian at points on 𝒮̂ are

    λ_s = -√(2)^1/4,   λ_u = √(2)^1/4,   λ_c = 0,0.

Hence, 𝒮̂
is a surface of saddle points.
The associated stable, unstable, and center subspaces are

    𝔼_s = span{[ [      2√(2);      4^1/4;  3 ^1/4q_1; 3 ^1/4 a_1 ]] },
      𝔼_u = span{[ [     -2√(2);      4^1/4;  3 ^1/4q_1; 3 ^1/4 a_1 ]] },
      𝔼_c 
                 =  span{[ [ 0; 0; 1; 0 ]],
                 [ [ 0; 0; 0; 1 ]] }.

Therefore, this subsystem has a two-dimensional center manifold W^c(𝒮̂) ∩{ r_1=0 }∩{δ_1=0 }. 
It emanates from ℓ, and W^c(𝒮̂) ∩{ r_1=0 }∩{δ_1 =0 } contains
W^c(ℓ) ∩{ r_1=0 }∩{δ_1 =0 }∩{ a_1=0 }.

Up next is the invariant hyperplane {δ_1 = 0 }.
The full system (<ref>) is fifth-order here:

    ṙ_1    = 1/2√() p_1 r_1 
    ṗ_1    = 1 - v_1 - 3/2√() p_1^2
                                  + 1/3√() r_1^2 
    v̇_1    = -2 √() p_1 v_1 
    q̇_1    = -3/2√() p_1 q_1 
    ȧ_1    = 3/2√() p_1 a_1.

The line I given by (<ref>) is invariant, and E_± are fixed points on it.
The equilibrium E_+ is a saddle with stable eigenvalues
-√(6)^1/4,
-2√(2/3)^1/4,
-√(3/2)^1/4, and -√(3/2)^1/4, and unstable eigenvalue 1/√(6)^1/4.
The associated stable and unstable eigenspaces are

    𝔼_s = span{[ [ 0; 1; 0; 0; 0 ]],
      [ [        0;    -√(3); √(2)^1/4;        0;        0 ]],
      [ [ 0; 0; 0; 1; 0 ]],
      [ [ 0; 0; 0; 0; 1 ]]
    },   𝔼_u = span{[ [ 1; 0; 0; 0; 0 ]] }.

The equilibrium E_- is a saddle with stable eigenvalue
-1/√(6)^1/4
and unstable eigenvalues √(3/2)^1/4,
√(3/2)^1/4, 
2√(2/3)^1/4, √(6)^1/4,
The associated stable and unstable eigenspaces are 

    𝔼_s =span{[ [ 1; 0; 0; 0; 0 ]] },
      𝔼_u  = span{[ [ 0; 0; 0; 1; 0 ]],
      [ [ 0; 0; 0; 0; 1 ]],
      [ [        0;     √(3); √(2)^1/4;        0;        0 ]],
      [ [ 0; 1; 0; 0; 0 ]]
    }.


There is also a two-dimensional surface of equilibria

    𝒮̃={ r_1=0,  δ_1=0,   p_1=0,  
                           v_1=1 + 1/3√() r_1^2,   
                           q_1 ∈ℝ,   a_1 ∈ℝ}

The eigenvalues of the Jacobian at points on 𝒮̃ are

    λ_s = -√(2√() +  r_1^2),   λ_u = √(2√() +  r_1^2),   λ_c = 0, 0, 0.

Hence, 𝒮̃
is a surface of saddle points. 
The associated stable and unstable subspaces are

    𝔼_s = span{[ [           -√() r_1;  2√(2√() +  r_1^2); 4/3√()(3+√()r_1^2);           3 √()q_1;          3 √() a_1 ]] },
      𝔼_u = span{[ [             -√() r_1;   -2√(2√() +  r_1^2); 4/3√()(3 + √()r_1^2);            3 √() q_1;            3 √() a_1 ]] }.

The associated center subspace is 

    𝔼_c 
                 =  span{[ [ 0; 0; 0; 1; 0 ]],
                 [ [ 0; 0; 0; 0; 1 ]], 
            [ [       3;       0; 2√()r_1;       0;       0 ]] 
           }.

Therefore, this subsystem has a three-dimensional center manifold W^c(𝒮̃) ∩{ r_1=0 }∩{δ_1=0 }. 
It emanates from ℓ, coincides with the critical manifold S, and W^c(𝒮̃) ∩{δ_1 =0 } contains
W^c(ℓ) ∩{δ_1 =0 }∩{ a_1=0 }.

In the invariant hyperplane { r_1=0 }, system (<ref>) is the following fifth-order set of equations:


    d δ_1/dy_1   = -√()p_1 δ_1 
    dp_1/dy_1   = 1-v_1-3/2√()p_1^2 
    dv_1/dy_1   = √()(-2p_1 v_1 + δ_1 q_1 ) 
    dq_1/dy_1   = -3/2√() p_1 q_1 +δ_1 
    da_1/dy_1   = -3/2√() p_1 a_1.

The line I given by (<ref>) is still invariant for this system, and E_± remain as fixed points on I, recall (<ref>).
The equilibrium E_+ is stable with spectrum 

    σ^s_+ = { -√(6)^1/4,
    -2√(2/3)^1/4,
    -√(3/2)^1/4,
    -√(3/2)^1/4,
    -√(2/3)^1/4}

At E_+, the stable subspace is 

    𝔼^s =  span  {[ [ 0; 1; 0; 0; 0 ]],
    [ [        0;    -√(3); √(2)^1/4;        0;        0 ]],
    [ [ 0; 0; 0; 1; 0 ]],
    [ [ 0; 0; 0; 0; 1 ]],
    [ [ ^1/4;    0;    0; √(6);    0 ]]
    }


The other equilibrium, E_-, is unstable with spectrum given by

    σ_-^u = {√(2/3)^1/4,
        √(3/2)^1/4,
        √(3/2)^1/4,
    2√(2/3)^1/4, √(6)^1/4}.

At E_-, the unstable eigenspace is 

    𝔼^u =  span {[ [ -^1/4;     0;     0;  √(6);     0 ]],
    [ [ 0; 0; 0; 1; 0 ]],
    [ [ 0; 0; 0; 0; 1 ]],
    [ [        0;     √(3); √(2)^1/4;        0;        0 ]],
    [ [ 0; 1; 0; 0; 0 ]]
    }.


There is also a manifold of equilibria

    𝒮_0 = { r_1=0,  δ_1=0,   p_1=0,   v_1=1,   q_1 ∈ℝ,   a_1 ∈ℝ}.

It is a manifold of saddle fixed points, since the eigenvalues are

    λ_s = -√(2)^1/4,   λ_u = √(2)^1/4, 
              λ_c = 0, 0, 0.

The associated eigenspaces are

    𝔼^s =  span {[ [          0;     2 √(2);      4^1/4; 3 ^1/4 q_1;   3^1/4a_1 ]]
                 },
                   𝔼^u =  span {[ [         0;  - 2 √(2);     4^1/4; 3^1/4 q_1;  3^1/4a_1 ]] },
    
    
    
      𝔼^c =  span {[ [ 0; 0; 0; 1; 0 ]],
    [ [ 0; 0; 0; 0; 1 ]],
    [ [    4/4-3√()q_1^2; 2q_1/4-3√()q_1^2;                0;                0;                0 ]]
    }.

Therefore, in the hyperplane { r_1=0 }, there is 
a three-dimensional
center manifold

    N = W^c(ℓ).

This manifold contains the surface 𝒮_0 of equilibria.
Moreover, as shown in Section <ref>, these key manifolds persist in the full system (<ref>) for 0< r_1 ≪ 1.








§ THE PROOF OF LEMMA LG
 



   The proof is split into two steps. First, we calculate a regular perturbation expansion in powers of r_2 of the solution Γ_0 in the H_2=0 level set. This expansion is valid in K_2 on arbitrary finite intervals of the independent variable y_2. Second, we identify the fixed points on the equator of the hemisphere to which this perturbed solution limits as y_2 →±∞, thereby establishing the persistence of the heteroclinic connection, which corresponds to the intersection of the two invariant manifolds. 
   
   All calculations in these steps were performed both by hand and symbolically.

Step 1.
Let 

    Γ_δ =  Γ_0 + r_2 Γ_1 + r_2^2 Γ_2 + r_2^3 Γ_3 + 𝒪(r_2^4),    with  Γ_k=(u_2k,p_2k,v_2k,q_2k),   k=1,2,….

We substitute this expansion into (<ref>) and solve for Γ_k, order by order in powers of r_2. At 𝒪(r_2^0), we recover Γ_0, as given by (<ref>). 
At 𝒪(r_2), 
the equations are

    u_21'    = √() p_21
    p_21'    = 2 u_20u_21 - v_21
    v_21'    = √() q_21
    q_21'    = u_21 - a_2.

In the algebraic solution, there is a free parameter χ_21∈ℝ: 

    u_21(y_2)    = χ_21 y_2 + 3/2 a_2,   
    p_21(y_2)=1/√()χ_21, 
    
    v_21(y_2)   =1/6√()χ_21y_2^3 + 1/4√() a_2 y_2^2,   
    q_21(y_2) =1/2χ_21 y_2^2 +1/2 a_2 y_2.


At 𝒪(r_2^2),  the equations are

    u_22'    = √() p_22
    p_22'    = 2 u_20u_22 + u_21^2 - v_22 + 1/3√() u_20^3 
    v_22'    = √() q_22
    q_22'    = u_22.

Here, there is also a free parameter χ_22∈ℝ
in the algebraic solution. We find 

    u_22(y_2)    = 3/√()χ_21^2  + 5/96√() + χ_22 y_2 - 5^3/2/3,456 y_2^4, 
    
    p_22(y_2)   =1/√()χ_22 - 5/864 y_2^3, 
    
    v_22(y_2)   =9/4 a_2^2 + 3 a_2 χ_21 y_2 + ( 3/2χ_21^2 + 5/192) y_2^2 + 1/6√()χ_22y_2^3 - ^2/20,736 y_2^6, 
    
    q_22(y_2)    =3/√() a_2 χ_21 
    +(  
    3/√()χ_21^2 + 5√()/96) y_2 + 1/2χ_22 y_2^2 - ^3/2/3,456 y_2^5.


At 𝒪(r_2^3),  the equations are

    u_23'    = √() p_23
    p_23'    = 2 u_20u_23 + 2u_21u_22 - v_23 + √() u_20^2 u_21
    v_23'    = √() q_23
    q_23'    = u_23.

The algebraic solution at this order has the free parameter χ_23∈ℝ:

    u_23(y_2)    = 6/√()χ_21χ_22 + χ_23 y_2 - 7/96 a_2 y_2^2 - 5/144χ_21 y_2^3,
    
    p_23(y_2)   =1/√()χ_23 - 7√()/48 a_2 y_2 - 5√()/48χ_21 y_2^2, 
    
    v_23(y_2)   = a_2 ( 9/√()χ_21^2 + 29 √()/96) 
    +( 3 a_2 χ_22 + 6/√()χ_21^3 + 5√()/16χ_21) y_2 
       + 3 χ_21χ_22 y_2^2 + √()/6χ_23 y_2^3 - 7 ^3/2/1,152a_2 y_2^4
    -^3/2/576χ_21 y_2^5, 
    
    q_23(y_2)    =1/√()( 3 a_2 χ_22 + 6/√()χ_21^3 + 5√()/16χ_21)
    +6/√()χ_21χ_22 y_2 
       + 1/2χ_23 y_2^2 - 7/288 a_2 y_2^3 - 5/576χ_21 y_2^4.


Step 2.
We substitute the expansions for (u_2,p_2,v_2,q_2) from Step 1 into the Hamiltonian H_2 given by (<ref>). After lengthy calculations,  we find that all of the terms that depend on powers of y_2 vanish, and one is left with:

    H_2 |_Γ_δ =
        -1/12 a_2 r_2 - 5√()/576 r_2^2
        -1/6√()a_2^2 r_2^4
        +𝒪(r_2^5).

Next, we multiply H_2 by 12√() r_2^2, recall that
ã=√()r_2^3 a_2 by (<ref>), 
and impose H_2 |_Γ_δ = 0. This yields

    ã = - 5/48 r_2^4 + 𝒪(r_2^6).

Finally, we recall that δ=r_2^2 in chart K_2 and a=1+ã by  (<ref>).
Therefore, we have derived the formula a_c(δ) = 1 - 5/48δ^2 + 𝒪(δ^3), given in (<ref>).








99







AACS2024
A. Asch, M. Avery, A. Cortez, and A. Scheel,
Slow passage through the Busse balloon
–Predicting steps on the Eckhaus staircase,
Eur. J. Appl. Math. (2024) in press.

ADK2017 D. Avitabile, M. Desroches, and E. Knobloch, Spatiotemporal canards in neural field equations, Phys. Rev. E 95 (2017), 042205. 

ADKK2017
D. Avitabile, M. Desroches, E. Knobloch, and M. Krupa, Ducks in space: From nonlinear absolute instability to noise-sustained structures in a pattern-forming system, Proc. Roy. Soc. A 473 (2017), 20170018.

ADVW2020
D. Avitabile, M. Desroches, R. Veltz, and M. Wechselberger, Local Theory for Spatio-Temporal Canards and Delayed Bifurcations,
SIAM J. Math. An. 52 (2020), 5703–5747.

BCDD1981
E. Benoit, Callot, F. Diener, and M. Diener,
Chasse au canard,
Collectanea Mathematica 31-32 (1981), 37–119.

B2013
M. Brøns,
An iterative method for the canard explosion in general planar systems,
Discrete
and Continuous Dynamical Systems Supplement 2013 (2013), 77–83.

BBE1991
M. Brøns and Bar-Eli,
Canard explosion and excitation in a model of the Belousov-Zhabotinskii reaction,
J. Phys. Chem., 95 (1991), 8706-8713.

Brons2006
M. Brøns, M. Krupa, and M. Wechselberger,
Mixed Mode Oscillations Due to the Generalized Canard Phenomenon,
Fields Institute Communications 49 (2006), 39–63.

BDHL2023
C. Brown, G. Derks, P. van Heijster, and D.J.B. Lloyd,
Analysing transitions from a Turing
instability to large periodic patterns
in a reaction-diffusion system,
Nonlinearity 36 (2023), 6839–-6878.

Buchholtz1995
F. Buchholtz, M. Dolnik, and I.R. Epstein,
Diffusion-induced instabilities near a canard,
J. Phys. Chem. 99 (1995), 15093–15101.

Buric2006
L. Buřič, A. Klíč, and L. Purmová,
Canard solutions and travelling waves in the spruce budworm population model,
Appl. Math. Comput. 183 (2006), 1039–1051.

B1978
F.H. Busse,
Nonlinear properties of thermal convection,
Rep. Prog. Phys. 41 (1978), 1929–1967.

BC1979
F.H. Busse and  R.M. Clever
Instabilities of convection rolls in a fluid of moderate Prandtl number,
J. Fluid Mech. 1 (1979), 319–-335.

BW1971
F.H. Busse and J.A. Whitehead,
Instabilities of convection rolls in a high Prandtl number fluid,
J. Fluid Mech. 47 (1971), 305-–320.


C1981
J. Carr,
Applications of centre manifold theoy, Applied Mathematical Sciences series, 35, Springer Verlag, 
New York (1981). 






CKW2017
P. Carter, E. Knobloch, and M. Wechselberger,
Transonic canards and stellar wind,
Nonlinearity 30 (2017), 1006–1033.

CH1993
M.C. Cross and P.C. Hohenberg,
Pattern formation outside of equilibrium,
Rev. Mod. Phys. 65 (1993) 851–1112.


CN1984
M.C. Cross and A.C. Newell
Convection patterns in large aspect ratio systems,
Physica D 10 (1984), 299–328.

DPK2009
P. de Maesschalck, N. Popovic, and T.J. Kaper,
Canards and bifurcation delays of spatially homogeneous and inhomogeneous types in reaction-diffusion equations, 
Advances Differential Equations 14 (2009) 943–962.

DKO2010 
M. Desroches, B. Krauskopf, and H. M. Osinga, Numerical continuation
of canard orbits in slow-fast dynamical systems, 
Nonlinearity 23 (2010), 739–-765.

D1984
M.J. Diener, The canard unchained or how fast/slow dynamical systems
bifurcate, 
Math. Intell. 6, (1984), 38–49.

D2019
A. Doelman, 
Pattern formation in reaction-diffusion systems – an
explicit approach, 
in Complexity Science,  M. Peletier, R. van
Santen, and E. Steur (eds.), 129-182  (2019)
World Scientific. 

DRS12
A. Doelman, J.D.M. Rademacher and S. van der Stelt,
Hopf dances near the tips of Busse balloons,
Disc. Cont. Dyn. Syst. Series S
5 (2012) 61–92.

AUTO 
E.J. Doedel, A.R. Champneys, T.F. Fairgrieve, Y.A. Kuznetsov, K.E. Oldeman,
R.C. Paffenroth, B. Sandstede, X.J. Wang, and C. Zhang, 
AUTO-07P: Continuation and bifurcation software for ordinary differential equations, Technical
Report, Concordia University, Montreal, Canada (2007).

DR1996
F. Dumortier and R. Roussarie, 
Canard cycles and center manifolds, Memoirs of the AMS 557 (1996), Amer. Math. Soc, Providence, Rhode Island. 

E1965
W. Eckhaus,
Studies in Non-Linear Stability Theory,
Springer Tracts in Natural Philosophy 6,
Springer, Berlin (1965).

E1983 
W. Eckhaus, 
Relaxation oscillations including a standard chase on French ducks, in Asymptotic Analysis II, F. Verhulst (ed.) Springer, Berlin (1983), 449–497.

E1993
W. Eckhaus,
The Ginzburg-Landau manifold is an attractor,
J. Nonlin. Sci. 3 (1993), 329–348.


EK2005
L. Edelstein-Keshet,
Mathematical Models in Biology, Classics in Applied Mathematics, Series Number 46, Society Industrial Applied Mathematics, Philadelphia (2005).

EHKPPZ2022
M. Engel, F. Hummel, C. Kuehn, N. Popović, M. Ptashnyk and T. Zacharis,
Geometric analysis of fast-slow PDEs with fold singularities,
preprint (2022),  arXiv:2207.06134
	
EP1998
I.R. Epstein and J.A. Pojman,
An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos,
Oxford University Press, Oxford, UK (1998).

Fenichel1979
N. Fenichel,
Geometric singular perturbation theory for ordinary differential equations, 
J. Differ. Eq.
31 (1979), 53–98. 

GV2007
J. Galan-Vioque and A. Vanderbauwhede, Continuation of periodic orbits in symmetric Hamiltonian systems,
 in Numerical Continuation Methods for Dynamical Systems, 
 B. Krauskopf, H.M. Osinga, and J. Galan-Vioque (eds.), Springer (2007) 269–299.

GZK2018
P. Gandhi, Y.R. Zelnik, and E. Knobloch,
Spatially localized structures in the Gray–Scott model, 
Philos. Trans. A: Math. Phys. Eng. Sci.
376
(2018), 20170375.

GKSV2023
R. Goh, T.J. Kaper, A. Scheel, and T. Vo,
Fronts in the wake of a parameter ramp: slow passage through pitchfork and
fold bifurcations,
SIAM J. Appl. Dyn. Sys. 22 (2023), 
2312–2356.

GKS2024
R. Goh, T.J. Kaper, and A. Scheel,
Pitchfork bifurcation along a slow parameter
ramp: Coherent structures in the critical scaling,
Studies Appl. Math. (2024), 1–21. 

GKV2022
R. Goh, T.J. Kaper, and T. Vo,
Delayed Hopf Bifurcation and Space–Time Buffer Curves in the Complex Ginzburg–Landau Equation,
IMA J. Appl. Math. 87 (2022), 131-–186.


H1980
J. Hale,
Ordinary Differential Equations,
Krieger Pub., Malabar, Florida (1980).

HI2011
M. Haragus and G. Iooss,
Local Bifurcations, Center Manifolds, and Normal Forms in Infinite-dimensional Dynamical Systems,
Springer (2011). 

H1991
A. van Harten,
On the validity of the Ginzburg-Landau equation,
J. Nonlin. Sci. 1 (1991), 397–422.

Hasan2018
C.R. Hasan, B. Krauskopf and H.M. Osinga,
Saddle Slow Manifolds and Canard Orbits in ℝ^4 and Application to the Full Hodgkin-Huxley Model,
J. Math. Neurosci. 8 (2018), 5.
	
HJK2022
F. Hummel, S. Jelbart and C. Kuehn,
Geometric blow-up of a dynamic Turing instability in the Swift-Hohenberg equation, preprint 
(2022) arXiv:2207.03967.

IMD1989
G. Iooss, A. Mielke, and Y. Demay,
Theory of Ginzburg-Landau problems in hydrodynamic stability problems,
European J. Mech. B Fluids 8
(1989), 229–268.

IP1993
G. Iooss and M.C. Peroueme, 
Perturbed homoclinic solutions in reversible 1:1 resonance vector fields,
J. Diff. Eq. 103 
(1993), 62–88.

JK2024
S. Jelbart and C. Kuehn, 
A formal geometric blow-up method for pattern forming systems,
Contemp. Math. (2024), to appear.

Jones1995
C.K.R.T. Jones,
Geometric singular perturbation theory, 
in Dynamical Systems, R. Johnson, ed., Lecture Notes in Math. 1609, Springer, New York, 1995, pp. 44–118. 

KV2018
T.J. Kaper and T. Vo,
Delayed loss of stability due to the slow passage through Hopf bifurcations in reaction-diffusion equations.  
Chaos 28 (2018), 091103.

KV2021
T.J. Kaper and T. Vo,
A new class of chimeras in locally coupled oscillators with small-amplitude, high-frequency asynchrony and large-amplitude, low-frequency synchrony,
Chaos 31 (2021), 123111.


KLSBE2021
C. Konow, Z. Li, S. Shepherd, D. Bullara, and I.R. Epstein,
Influence of survival, promotion, and growth on pattern formation in zebrafish skin, Scientific Reports 11 (2021), 9864.

KS2001
M. Krupa and P. Szmolyan,
Extending geometric singular perturbation theory to nonhyperbolic points–fold and canard points in two dimensions,
SIAM J. Math. Anal. 33 (2001), 286–314.

KW2010
M. Krupa and M. Wechselberger,
Local analysis near a folded saddle-node singularity, 
J. Diff. Eq.
248 (2010), 2841–2888.

M1982
H. Meinhardt,
Models of Biological Pattern Formation,
Academic Press, New York (1982).

MG2000
H. Meinhardt and A. Gierer,
Pattern formation by local self-activation and lateral inhibition,
Bioessays 22 (2000), 753–760.

MKKR1984
E.F. Mishchenko, Yu.S. Kolesov, A.Yu. Kolesov, N.Kh. Rhozov; 
Asymptotic Methods in Singularly Perturbed Systems, 
Monographs in Contemporary Mathematics; Consultants Bureau, Plenum Publishing, NY (1984).

Mitry2017
J. Mitry and M. Wechselberger,
Folded saddles and faux canards,
SIAM J. Appl. Dyn. Syst. 16 (2017), 546–596.

Moehlis
J. Moehlis, 
 Canards in a surface oxidation reaction,
 J. Nonlinear Sci. 12 (2002), 319–345. 

MDK2001
D.S. Morgan, A. Doelman, and T.J. Kaper,
Stationary periodic patterns in the 1D Gray–Scott model,
Meth. Appl. An. 7 (2001), 105–150.

M1993
J.D. Murray,
Mathematical Biology,
Biomathematics Texts, Springer, Berlin 19
(1993).

P1921
B. van der Pol, 
A theory of the amplitude of free and forced triode vibrations,
Rad. Rev. 1 (1920), 701–710.

P1926
B. van der Pol, 
On relaxation-oscillations, 
Philos. Mag. Ser. VII 2 (1926), 978. 

RKZE2003
H.G. Rotstein, N. Kopell, A.M. Zhabotinsky, and I.R. Epstein,
Canard phenomenon and localization of oscillations in the Belousov-Zhabotinsky reaction with global feedback,
J. Chem. Phys. 119 (2003), 8824–8832.

Roussel1990
M.R. Roussel and S.J. Fraser,
Geometry of the steady-state approximation: Perturbation and accelerated convergence methods,
J. Chem. Phys. 93 (1990), 1072–1081.

S2003
W. van Saarloos,
Front propagation into unstable states,
Physics Reports 
386 (2003), 29–222.

SS01
B. Sandstede and A. Scheel,
Essential instabilities of fronts: bifurcation, and bifurcation failure,
Dynamical Systems 
16 (2001) 1–28. 

SU2017
G. Schneider and H. Uecker,
Nonlinear PDEs: A Dynamical Systems Approach,
Grad. Studies Math. 182,
American Mathematical Society, Providence (2017).

SW2022
G. Schneider and M. Winter, 
The amplitude system for a simultaneous short-wave Turing and long-wave Hopf instability,
Disc. Cont. Dyn. Sys. S
15 (2022) 2657–2672.

Szmolyan2001
P. Szmolyan and M. Wechselberger, 
Canards in ℝ^3,
J. Differ. Eq.
177 (2001) 419–453.

T1952
A.M. Turing,
The chemical basis of morphogenesis,
Phil. Trans. Royal Soc. London B Bio. Sciences 237 (641) (1952), 37–72.

VSC2023 
E. Vilar-Sepulveda and A. Champneys, Degenerate Turing bifurcation and the birth of localized patterns in activator-Inhibitor systems, SIAM J. Appl. Dyn. Sys.  22 (2023),  1673–1709.

VBK2020
T. Vo, R. Bertram, and T.J. Kaper,
Multi-mode attractors and spatio-temporal canards,
Physica D 411 (2020), 132544.



W1997
D. Walgraef,
Spatio-temporal Pattern Formation: With examples from physics, chemistry, and material science,
Springer, New York  (1997).

ZS1984
A.K. Zvonkin and M.A. Shubin,
Non-standard analysis and singular perturbations of ordinary differential equations,
Russ. Math. Surveys
39 (1984), 69–132.








