Bifurcation of critical sets and relaxation oscillations in singular fast-slow systems

Fast-slow dynamical systems have subsystems that evolve on vastly different timescales, and bifurcations in such systems can arise due to changes in any or all subsystems. We classify bifurcations of the critical set (the equilibria of the fast subsystem) and associated fast dynamics, parametrized by the slow variables. Using a distinguished parameter approach we are able to classify bifurcations for one fast and one slow variable. Some of these bifurcations are associated with the critical set losing manifold structure. We also conjecture a list of generic bifurcations of the critical set for one fast and two slow variables. We further consider how the bifurcations of the critical set can be associated with generic bifurcations of attracting relaxation oscillations under an appropriate singular notion of equivalence.


Introduction
Many natural systems are characterized by interactions between dynamical processes that run at very different timescales. These can often be modelled as fast-slow systems, where system dynamics can be separated into interacting fast and slowly changing variables. This has been applied to a wide range of natural phenomena, including electrical circuits [39], plasma oscillations [34], surface chemistry [27] and cell physiology [24] to ecology [41] and climate [3]. The dynamical behavior of such systems can often be understood in a common mathematical framework. See [29] for a recent monograph that summarizes both techniques and applications, and [2,6,7,12,18,21,23,45] for examples of related work.
The analysis of fast-slow systems is built around a geometric singular perturbation theory (GSPT) approach [13], perturbing from a singular limit where timescales decouple: see also [14,25]. In the singular limit, on the slow timescale there are instantaneous jumps (determined by the fast dynamics) between periods of slow evolution. The slow evolution typically takes place on stable sheets of a critical set arXiv:1902.09203v3 [math.DS] 18 Sep 2019 where the fast dynamics is in stable balance, interspersed by fast transitions. If the fast dynamics is one dimensional, then it is typically confined to a manifold (and hence the set is often called a critical manifold), though at bifurcation it may lose its manifold structure at isolated points. The fast transitions are determined by what we call the umbral map defined as the map from a fold point on the critical set to another part of the critical set. In the case of stable periodic behaviour in the singular limit, the resulting limit cycles are referred to as singular relaxation oscillations. Many examples of bifurcations of relaxation oscillations have been considered [2], including some associated with bifurcations of the critical set [3,15,16] although it seems that no exhaustive list of scenarios has been proposed.
Although singular perturbation theory has been developed to explain many aspects of behaviour near the singular limit, there is still no full understanding of generic bifurcations of limit cycles even in the singular limit. Guckenheimer [18,19,20] suggests an approach and several results along these lines, but, as far as we are aware, these conjectures are yet to be framed, let along proved, in rigorous terms. The main aim of this paper is to present an approach to doing precisely this, using singularity theory with distinguished parameters and appropriate notions of equivalence. We classify local and global transitions in the critical set by codimension and consider the consequences for the umbral map. In doing so, we find a variety of scenarios that give bifurcation of attractors in such singular systems.
We show that it is possible to split the problem of bifurcations of relaxation oscillations into two aspects: bifurcations of the critical set, and bifurcations caused by singularities of the slow flow (possibly interacting with the critical set). In the simplest case of one fast and one slow variable, bifurcations of the critical set can be directly tackled using a global version of the singularity theory with distinguished parameter in [17]. We highlight that this theory needs extension to make it suitable for systems with multiple fast and slow variables. We are able to identify a large number of scenarios that can lead to bifurcation of relaxation oscillations. Note that we only consider fast-slow systems where the fast dynamics is constrained to a subset of the variables; following [13] it is possible to extend the theory developed here to more general fast-slow systems that are not in standard form: see for example [28,30,48].
We structure the paper as follows: in Section 2 we briefly introduce the singular limit of fast-slow systems, critical sets, singular trajectories and global equivalence of critical sets. In Section 3 we explore persistence and bifurcation of critical sets by examining versal unfoldings of the fast dynamics parametrized by the slow variables, using a notion of global equivalence of the fast dynamics. In the case of one fast and one slow variable we classify persistence (Proposition 1) and bifurcations of the critical set up to codimension one (Proposition 2) using the theory of [17]. For one fast and several slow variables we highlight the need for an improved theory of bifurcations with multiple distinguished parameters. We present conjectured statements of persistence and of codimension one bifurcations of the critical set (Conjectures 1 and 2 respectively) for one fast and two slow variables. We find a rich variety of distinct mechanisms for typical codimension one bifurcations of the critical set which includes local and global bifurcations in the fast variable.
Section 5 turns to the question of bifurcation of attractors in fast-slow systems and in particular bifurcation of stable relaxation oscillations. We introduce a global singular equivalence for the singular trajectories and use this to classify persistence (Proposition 4) and codimension one bifurcations (Proposition 5) of these simple relaxation oscillations. These codimension one bifurcations naturally split into those caused by bifurcations of the critical set, and those caused by interaction of singularities of the slow flow with the critical set: in Section 5.3 we present some numerical examples of various types. Finally, Section 6 is a discussion of some of the challenges for GSPT to describe the unfolding of such bifurcations, and relation to other singularity theory approaches. We include several Appendices that give more details of the tools used for the classification and the examples.

Singular trajectories of fast-slow systems limitintro
A fast-slow system is a system of coupled ODEs for z = (x, y) ∈ R m+n of the form ẋ = g(x, y, ) y = h(x, y, ) (1) eq:mainsystem where x ∈ R m and y ∈ R n , > 0 is a small constant and t is a time measured relative to the slow dynamics. The functions g(x, y, ) and h(x, y, ) are C ∞ functions of their arguments (they are well approximated by Taylor series to arbitrary order). We refer to x as the fast and y as the slow subsystems; the dynamics in these subsystems are referred to as fast and slow dynamics respectively. These systems have a singular limit → 0, where a typical trajectory can remain close to an equilibrium of the fast system, except at isolated points where it "jumps" along a trajectory of the fast subsystem. The singularly perturbed system with > 0 will have trajectories that typically remain near a trajectory of the singular limit, although especially near bifurcations, trajectories may also explore unstable parts of the slow dynamics along so-called canard solutions (see e.g. [32,6,2,45] and [29,Chapter 8]).
In order to understand such systems it is useful to consider the reduced or slow equations in slow time t: 0 = g(x, y, 0) y = h(x, y, 0), (2) eq:reduced describing the slow dynamics in the singular limit → 0. Solutions of (2) are constrained to the critical set C[g] = (x, y) ∈ R m+n : g(x, y, 0) = 0 .
(3) eq:critman Note that by the regular value theorem [33], this critical set is a manifold at all points where the derivative of g has maximal rank. For an open dense set of g ∈ C ∞ this is true for an open and dense set (x, y) ∈ C[g]. The set is often called a critical manifold, however we do not use this notation as at bifurcation points the set may lose its manifold structure.
The flow g may have singular equilibria of the fast dynamics within C[g], where the term singular here means non-regular. The regular points of the critical set C[g] we define as The remaining non-hyperbolic (fold) points, also called singular or limit points, form the fold set of the critical set : ∂ x g(x, y, 0) is non-hyperbolic} .
At all regular points the reduced equations (2) can be used to describe the flow. At fold points, however, we need to consider the fast dynamics and there will be jumps in slow time determined by the fast subsystem only. Changing variable to a fast time τ = t/ and taking the limit → 0 gives quite a different set of equations: the layer or fast equations: x = g(x, y, 0) y = 0, (6) eq:layer where we write x to denote d dτ x and note that (a) the constant slow variable y now acts as a parameter for evolution of the fast variable x and (b) the layer equation, when restricted to C[g] consists entirely of equilibria for m = 1 (for m > 1 there may be other objects, such as limit cycles). We split the regular set into a disjoint union of attracting/repelling/saddle points C att [g], C rep [g], C sad [g] so that C reg [g] = C att [g] ∪ C rep [g] ∪ C sad [g]. Note that F[g] is the union of the set of boundaries of these sets, and also that the set C sad [g] only exists for m ≥ 2. Note that the regular set is the union of all non-singular points From now on, we only concern ourselves with the system in the singular limit, and therefore drop the dependency on in our notation, such that e.g. g(x, y, 0) := g(x, y) and h(x, y, 0) := h(x, y). However, we stress once more that the system (1) may depend crucially on near, but away from, the singular limit.
We now make the notion of jumps more precise: For any point p = (x, y) ∈ F[g] we define the (possibly set-valued) umbral map to be U [g](p) = {ω-limits of a non-trivial trajectories in (6) with α-limit p}, that is, a map from every point q, which limits to p in backward time (α-limit), to the set of forward time limits (ω-limit) of every such q, excluding p itself. The ω-limits are always non-empty since we will assume bounded global attractors. However, if all ω limits equal p then the umbral map is empty. The umbra (meaning shadow [18]) or drop set [29,Ch.5,p.109] is the image of the folds under the umbral map: We define the projection onto the slow variable as π : R m+n → R n , π(x, y) = y and for p = (x, y) we define the set of all co-folds to p as i.e. all fold points with the same slow coordinate as p. Similarly, we define the set of folds sharing a given slow (y) coordinate to be:

Trajectories in the singular limit
Note that typical points in R m+n are not on C[g]: starting at (x, y) ∈ C[g] there will be fast motion governed by the layer equations (6). If this settles to a limit we will typically have arrived at a point on the critical set that is a stable equilibrium, i.e. on C att [g]. The slow dynamics then carries the trajectory around C att [g] until (possibly) it hits a fold point p = (x, y) ∈ F[g]. If U [g](p) is a single point then there is a unique non-trivial trajectory of the layer equations from this point, the fast motion will take the dynamics to U [g](p). Hence typical trajectories in the singular limit (which can be viewed as trajectories of a constrained differential equation [42]) are composed of segments of slow trajectories on C att [g] interspersed with fast jumps. Depending on the nature of the slow segments and fast jumps, characterised by the umbral maps, there may be a trajectory of the > 0 system that remains close to the singular trajectory. More precisely, we define a singular trajectory as follows (this idea is widely used and sometimes called a candidate trajectory, for example [44,32,4,5,35] and [29,Ch.3,p.64]): • The interval is partitioned as a = s 0 < s 1 < . . . < s n−1 < s m = b into a finite number of subintervals.
• The image of each subinterval γ 0 (s j−1 , s j ) is a trajectory of either the fast subsystem or the slow subsystem.
• The image γ 0 (a, b) is oriented consistently with the orientations on each subinterval induced by the fast or slow flows.
Note that s is a parametrisation of the curve rather than fast or slow time. Consequently, the image of a subinterval can be complete homoclinic or heteroclinic orbits of the fast or slow subsystem. In typical cases, the attractor will consist of subintervals that alternate between fast and slow segments, but this may not be the case at bifurcation. In the case that all slow segments are on C att [g], this will typically perturb [35] to similar trajectories → 0. If the slow segments explore other hyperbolic points on the critical manifold, canard trajectories may appear. Several possible cases of fast/slow trajectories are considered in [35].

Global equivalence of critical sets
In order to define persistence and bifurcation of critical sets, we need a notion of equivalence between critical sets. We define an equivalence of critical sets, called global equivalence, through the parametrized fast vector fields that generate them. We adopt the equivalence in [17] for the case m = n = 1, and conjecture that that equivalence or a similar one may work also for the case m = 1, n = 2. A local version of this equivalence due to [38] is presented at the end of this section. Classification of local bifurcations (classification of bifurcation of germs) has been worked out for the local equivalence [38], but it remains to be shown if the equivalence is also suitable for global classification of bifurcations.
We consider fast and slow vector fields g : R m+n → R m and h : R m+n → R n , with the generic assumptions of smoothness and being in the singular limit → 0. More precisely we consider (9) eq:restrictedsy We furthermore consider only systems (1) with bounded global attractors. We restrict ourselves to vector fields defined on some fixed compact regions M ⊂ R m and N ⊂ R n with open interior and smooth boundaries having outward normals m(x) and n(y) respectively. Implicitly fixing these M and N , we define the open sets g(x, y) · n(x) < 0 ∀(x, y) ∈ (∂M, N ) and g(x, y) = 0 ⇒ g x (x, y) = 0 ∀(x, y) ∈ (M, ∂N ) (10) eq:Xf Note that if (g, h) ∈ M × N then the forward dynamics must remain in M × N and that there are no tangencies of the flow, (or the critical set) with the boundary. These conditions ensure that the properties persist under small perturbations of the vector field.
Since global equivalence should imply local equivalence we expect the classification of local bifurcations under a global equivalence, perhaps (12), to coincide with the classification under a local equivalence such as (13). The latter was worked out in [38].
We leave the generalisation of the global equivalence (12) to the case m > 1, n > 1 open.

Persistence and bifurcation of critical sets lmanifolds
Assume we define V f as in (10) for some compact regions M and N . In order to define persistence of the critical sets we define the unfolding of the slow dynamics following [17, Section III]. We say a smooth function G(x, y, λ) for λ ∈ R r is an r-parameter unfolding of g(x, y) if G(x, y, 0) = g(x, y) for all (x, y) ∈ M × N . Reference [17] mostly assumes G and g are germs of vector fields, though in [17, Theorem III.6.1] the equivalence is global within a compact region, as we consider here.
If G and H are both unfoldings of g, we say that H factors through G if there exist smooth mappings S, X, Y, L and W ⊂ R r , a neighbourhood of 0, such that H(x, y, λ) = S(x, y, λ)G(X(x, y, λ), Y (y, λ), L(λ)), for all λ ∈ W and (x, y) ∈ (M, N ), where S(x, y, 0) = 1, X(x, y, 0) = x, Y (y, 0) = y, L(0) = 0 (see [17]). We define G to be a versal unfolding if every unfolding H of g factors through G. We say g is persistent if it is its own unfolding, i.e. for any unfolding G ∈ C ∞ (M ×N ×R r ) such that G(x, y, 0) = g(x, y), on M ×N there is a neighbourhood W of 0 in R r such that where, as before, ∼ denotes global equivalence on M × N .
If the unfolding is versal and contains a minimum number of parameters, we call it a universal unfolding [17]. The number of parameters λ in such a universal unfolding G is the codimension of g. In particular, if g is persistent then g is its own universal unfolding, and in this case we say it has codimension zero. We say that a bifurcation of g occurs if g is non-persistent, in which case the codimension of the bifurcation is that of the universal unfolding of g. We emphasise once more that the equivalence relation (12) concerns the zero sets of g, i.e. the critical set C[g]. Hence, persistence and bifurcation of g under this equivalence is identified with persistence and bifurcation of the critical set.
Note that the aforementioned meaning of persistence does not concern persistence to perturbations involving the scale separation parameter , which is sometimes the case in the fast-slow literature, e.g. [13,35]. Here, we study exclusively the system (1) in the singular limit → 0.

Persistence and codimension one bifurcation of critical sets for one fast and one slow variable
The case m = n = 1 can be directly treated using the global bifurcation theory with distinguished parameter approach of [17,Section III]. Here, a distinguished parameter is a parameter which is considered integral to the model and separate from unfolding parameters, which represent model perturbations. In the fast-slow setting we consider the slow variable y to be a distinguished parameter from the point of view of the layer equations (6). Consider some g ∈ V f and note that the critical set is and that the fold set is Table 1 lists the three degenerate fold sets D i [g], i = {1, 2, 3} for m = n = 1: fold tangency, hysteresis point, and multiple limit point. The term limit point is a historical term for fold point. The set of all degenerate folds is then  [17] refers to the fold tangency as a "simple bifurcation" and a multiple limit point as a "double limit point" but our notation offers easier generalization to higher n. The following theorem characterizes the persistent critical sets, using a result from [17].  (14) is empty then g ∈ V f is persistent on M × N Fold tangency: : g y (p) = 0}, Hysteresis point: : g xx (p) = 0}, Multiple limit point: Table 2. Subsets of D i [g] for m = n = 1 whose union contains all codimension one bifurcations. Note that det(D 2 g) = g xx g yy − g 2 xy and that the first two are local degeneracies Quadratic fold tangency: : |Π(p)| = 1 and det(D 2 g(p)) = 0}, Cubic hysteresis point: : |Π(p)| = 1 and g xxx (p) = 0}, Double limit point: To aid the classification of codimension one bifurcations, we define D 1 i [g] in Table 2. These are open dense subsets of D i [g] that avoids obvious further degeneracies. We then subdivide these cases further in Table 3. A vector field g is degenerate at codimension one if exactly one of these degeneracies D 1 i,j [g, y] occur for exactly one slow coordinate y (in exactly one fast fibre) which means we only need to compare points in P (y) for some y. We avoid higher codimension fold tangency by precluding the cases det(D 2 g) = g xx g yy − g 2 xy = 0 or higher order hysteresis g xxx = 0. Note that D 1 i,j [g, y] depends explicitly on y. This choice makes it easier notation-wise to preclude the critical set from being degenerate at multiple y, which would raise codimension. The next result shows that Table 3 and Figure 1 give a complete list of codimension one bifurcations for this case. Table 3. Complete list of degeneracies that lead to codimension 1 bifurcations listed in Proposition 2. We write P (y) = {p i } as the set of distinct singular points p i = (x i , y) of the vector field g with slow coordinate y. Note that local degeneracies have |P (y)| = 1 Hyperbolic fold tangency : |P (y)| = 1 and det(D 2 g(p)) < 0} Opposed umbrafold double limit: : |P (y)| = 2 and U [g](p 1 ) = p 2 and ν[g](p 1 ) · ν[g](p 2 ) < 0, for some p 1 , p 2 ∈ P (y)} Fig. 3 d,e,f) Aligned umbraumbra double limit: : |P (y)| = 2 and U [g](p 1 ) = U (p 2 ) and ν[g](p 1 ) · ν[g](p 2 ) > 0, for some p 1 , p 2 ∈ P (y)} for g ∈ V f are characterised in Fig. 1, such that one of the sets D j,k [g, y] in Table 3 is non-empty for precisely one y. At such a bifurcation, precisely one of the following occurs: (i) Two folds merge at a fold tangency (e.g. Fig. 2 a,b,c) or d,e,f )).
(iii) Two fold points share the same slow coordinate: there are six distinct ways this can occur (e.g. Fig. 3) Proof: To avoid persistence, at least one of the the degeneracies D i [g] listed in Table 1 Table 3). Note that ν(p) is the direction vector of a fold at a point p and det(D 2 (g(p))) is the Hessian of g at p. Similarly, f means fold, f u means fold umbra and f x means non-interacting fold. For ease of notation, we suppress explicit dependence on g, such that e.g. ν[g](p) = ν(p). For a persistent codimension one bifurcation exactly one branch must be followed for exactly one fast fibre (a single y), leading to one of the red boxes. Overlapping red boxes symbolise that the degeneracy can be of either aligned or opposed type. See the text for details acytable11 must occur for some i ∈ {1, 2, 3}: as these are independently defined we can assume that only one will occur for an open dense set of unfoldings. Without loss of generality we can assume that the open conditions in Table 2 apply.
The subcases of D 1 1 [g] follow from examining the sign of det(D 2 g): the hyperbolic fold tangency D 1,1 [g, y] is the simple bifurcation of [17] while the elliptic fold tangency D 1,2 [g, y] is also called the isola. Similarly, the cubic hysteresis D 1 2 [g] can be either stable or unstable, depending on the sign of the leading order term. These cases can be transformed into the normal forms of Table 4 (these are given in [17]). The cases D 1 1 [g] unfold on varying a typical parameter λ as shown in Fig. 2 a,b,c) and d,e,f) respectively, while D 1 2 [g] unfold as shown in Fig. 2 g,h,i) and j,k,l). The double limit point degeneracy D 1 3 [g] can be split into several subsets according to the direction of the folds given by the signs of ν[g](p) = g xx (p)g y (p) at the two limit points, and k, the number of regular sheets that separate them. The number k determines whether the umbrae and folds intersect. If k = 0, then the umbra of one fold intersects the other fold, if k = 1 then the umbrae of the folds intersect, and if k ≥ 2 then the umbrae and folds do not intersect. The six distinct subcases of D 1 3 [g] are shown in Figure 3.
We conclude this section with a few comments. First, it may appear as if the noninteracting double limit point degeneracy is not a degeneracy. However, the critical sets in e.g. Figure 3 m,o) cannot be equivalent to that in Figure 3 n) since the equivalence (12) preserves the number of zeros in fast fibres: at bifurcation there is a y for which C[g] has five zeros, while the perturbed diagrams have at most four zeros for any given y. However, non-interacting double limit point degeneracy does not cause bifurcation of singular relaxation oscillations, as is shown in Section 5.2.
It may seem from Figure 2 that the hyperbolic fold tangency is of codimension higher than one. That this is not the case is shown algebraically in [17] using a considerable technical machinery. However, in Figure 4 we aim to provide some intuition that the bifurcation does not require fine tuned perturbations, but rather arises generically as two branches of the critical set approach under one-parameter perturbation.
It might seem puzzling why some combinations of x and y are left out from the normal forms in Table 4. In general, which terms can be included in a certain normal form is a non-trivial question, treated in e.g. [17]. However, for the particular case of fold tangency the transformation X(x, y) = x − ay, Y (y) = y/ |1 − a 2 | turns the polynomial g(x, y) = x 2 + y 2 + 2axy + λ into one of the normal forms in Table 4 as long as g is non-degenerate. This transformation clearly preserves equivalence class under (12).

Persistence of critical sets for one fast and two slow variables
In analogy with the m = n = 1 case we give a conjectured list of all degeneracies that can cause nonpersistency of vector fields with one fast and two slow variables, up to codimension one. First, we introduce some notation. For any g ∈ V f we write g x = ∂g ∂x , ∇ y g = (g y 1 , g y 2 ) and ∇ ⊥ y g = (−g y 2 , g y 1 ).
By u||v we mean that vectors u and v are parallel. The non-zero vector u rescaled to unit length is denoted u = u/|u|. D 2 (g) is the Hessian of g with respect to all components Table 5. Singularities of the critical set for one fast and two slow variables Quadratic fold: The slow Hessian D 2 y (g) is defined analogously, but with i, j ∈ {2, 3}. Recall that the critical set is and the fold set is As folds are not typically isolated in this case, we also need to distinguish between quadratic folds, cubic cusps and higher order cusps (Table 5).
We believe that the list of degenerate sets D i [g] of F[g] given in Table 6 is an exhaustive list of degeneracies under a suitable equivalence. These degeneracies are natural extensions from the degeneracies for m = n = 1; generic objects (quadratic fold lines and cubic cusps) can intersect (D 1 [g], D 2 [g], D 3 [g]), their projections onto the slow variables can become tangent D 4 [g] or intersect D 5 [g] and D 6 [g]. More precisely, we define the set of degenerate points Note that the umbral map is single valued for then it can be zero, one or two-valued (see Figure D1). We now conjecture a persistence criterion for m = 1, n = 2 that is analogous to the m = n = 1 case in Proposition 1.
conj:CM120 Conjecture 1 (Codimension zero, m = 1, n = 2.) For one fast and two slow variables, the critical set Unfortunately, the method of proof in [17, Thm III.6.1] used for Proposition 1 does not easily generalize to this case of "multiple distinguished parameters" (two slow variables in our case). This is because degenerate cases appear with codimension infinity [36], at least if we consider the restricted global equivalence where we require that q(y) is the identity in (12). Hence proof of this result will require a less stringent (but still natural) form of global equivalence. Table 6. Possible degeneracies of the critical set for one fast and two slow variables, m = 1 and n = 2. As before, the co-fold set Π(p) is the subset in F[g] sharing slow coordinate with the point p Fold tangency Triple fold projection intersection A variety of degeneracies can persistently occur for one parameter families, i.e. at codimension one. Table 8 lists degeneracies that we believe contain all persistent codimension one bifurcations for a suitable notion of global equivalence. We divide these into local or global degeneracies. The local degeneracies (Table 9) , the others involve interaction of two or more points in the same fast fibre of y: these global degeneracies are denoted D 1 i,j [g, y] for i ∈ {4, 5, 6} and are listed in detail in Appendix D.
We note that our conjectured local codimension one degeneracies coincide with those in the classification of [38], which is found for the local equivalence (13); it might be possible to extend these results to a related global equivalence.
We first state necessary conditions for the degeneracies D i [g] in Table 6 to be codimension one. In this process we introduce some geometric notions, listed in Table  7.
Starting with quadratic fold degeneracies, we note that for typical fold tangency we require that the Hessian D 2 (g) has no zero eigenvalues. For typical fold projection tangency, we require that the two folds sharing slow coordinate (y 1 , y 2 ) are not aligned and with the same quadratic curvature. This is guaranteed if the quadratic fold curvature vectors κ[g](p) ( Table 7 and Figure 5 c,d)) at fold points p 1 and p 2 ∈ Π(p 1 ) are distinct: For typical triple quadratic folds we require that exactly three fold points occur for every slow coordinate (y 1 , y 2 ) and that all folds are quadratic.
Turning to cusps, there is degeneracy if the cubic cusp direction vector µ[g](p) at a cubic cusp point p (Table 7 and Figure 5 b)) is either zero or undefined. The case µ[g](p) = 0 corresponds to g xxx (p) = 0, resulting in a swallowtail bifurcation. The Table 7. Quantities used in the classification of degeneracies for one fast and two slow variables, m = 1 and n = 2. The last column refers to appendices where all definitions except ν[g](p) are motivated (we do not motivate the natural definition of ν[g](p)). See Figure 5 for illustrations of all quantities except W [g](p) Quadratic fold direction vector Scalar quadratic fold curvature  Table 7  efinitions tangency. However, for typical cusp tangency we require that the quadratic quantity Table 7) is non-zero. Furthermore, the sign of W [g](p) separates subcases of typical cusp tangency. Finally, for typical cusp-fold projection intersection, we require that exactly one cubic cusp and one quadratic fold share slow coordinate (y 1 , y 2 ).
The quantities ν[g](p), µ[g](p), W [g](p) and κ[g](p) are discussed more in Appendix A, Appendix B and Appendix C. Hypothesised normal forms for the local codimension one bifurcations are listed in Table 4.
We now go trough the persistent subcases of codimension one degeneracies listed in Table 8. Table 8. Subsets of D[g], for one fast and two slow variables, whose union we conjecture contains all codimension one bifurcation sets. Note that D 1 i [g] are local for i = 1, 2, 3 and global for i = 4, 5, 6. Π(p) is the set of all singular points sharing slow coordinate with p. The scalars a 1 , a 2 , a 3 are coefficients in equation (15) : |Π(p)| = 1 and det(D 2 g(p)) = 0} Typical cusp tangency : |Π(p)| = 1 and g xxxx = 0} Typical double quadratic-fold projection tangency Table 9. Subsets of local degeneracies, for one fast and two slow variables parametrized by the slow coordinate, conjectured to include all codimension one bifurcations. Note that P (y) is the set of all singular points of the vector field g with slow coordinate y. The number of positive eigenvalues of the Hessian sign (g xx )D 2 [g] at p is written as |Σ + |. Note that no eigenvalues are zero, since det(D 2 g(p)) = 0 by assumption. Figure  D1 is in Appendix D. See the text for details Wormhole fold tangency Fig. 6 g,h,i) Stable lips cusp tangency  Fig. 6 p,q,r)   Table 9 for a) to r), and Table D1, Table D2 and Table  D3 in the Appendix for s) to dd) ). Bifurcation occurs when the bifurcation parameter λ equals the critical value λ = λ 0 . Solid/dashed black lines show the stable/unstable sheets of the critical set while red lines show the image of the fold under the umbral map. Blue lines indicate special points of intersection 121DSample 3.3.1. Fold tangency Fold tangency occurs when a pair or continuum of folds intersect. Typical fold tangency D 1 1 [g] is classified by |Σ + |, the number of positive eigenvalues in Σ, the spectrum of sign (g xx (p))D 2 (g)(p). The cases of wormhole (|Σ + | = 1), tube (|Σ + | = 2) and isola (|Σ + | = 3) are shown in Figure 6 a) to i). Note that g xx (p) det(D 2 (g)(p) = 0 implies that all non-positive eigenvalues are negative and that at least one eigenvalue is positive.

Cusp tangency
At a cusp tangency, two cusps meet locally along a line. The subsets beaks and lips are distinguished by whether cusps are directed away from or towards each other before bifurcation (see Figure 6 and Figure D1). These degeneracies are named after the appearance of their projections onto the slow variables (Figure 10), and the type is determined by the cusp quantity W [g](p). The case W [g] > 0 gives "lips" and W [g] < 0 gives "beaks". We further subdivide these cases depending on their stability, determined by the sign of g xxx .

Swallowtail
The swallowtail ( Figure 6 p,q,r)) is well known from catastrophe theory and occurs when a fold "folds over itself" to create a degenerate fold that splits up into a pair of cusps.
3.3.4. Fold projection tangency At a fold tangency, the projection of two curves of folds onto the slow variables are tangent. We divide fold projection tangency D 1 4 [g] into subcases depending on whether the folds at points p 1 and p 2 approach each other from the same direction (aligned) or opposite directions (opposed), captured by the sign of the inner product of the quadratic fold direction vectors ν[g](p 1 ) · ν[g](p 2 ). We further subdivide the opposed cases depending on the sign of the sum where K[g](p) is the scalar quadratic fold curvature at a fold point p (see Table 7 and Figure 5 c,d)). K[g](p) > 0 corresponds to a quadratically convex fold with respect to the fold direction and K[g](p) < 0 corresponds to a quadratically concave fold. Hence means that the concave curvature dominates, and the degeneracy is called a covering fold projection tangency since the folds locally cover the slow plane (see Figure 7). Similarly, if then the degeneracy is called a non-covering fold projection tangency. Accounting for whether the fold umbrae interact with each other, or one fold umbra interacts with a fold, or neither, we get six subcases of opposed fold projection tangency (Table D1). If the two fold projections are aligned the total curvature does not matter as long as K[g](p 1 ) = K[g](p 2 ), with one exception. This exceptional case occurs if a fold umbra hits a fold, in which case it matters if the curvature of the umbral fold dominates the destination fold or not (Figure 7 a.i) and a.ii)). More details are listed in Appendix D.

Cusp-fold projection intersection
At a cusp-fold projection intersection, the projections of a cusp and a fold line coincide in their projection onto the slow variables. We classify the intersection D 5 [g] of a cusp and a fold projection into ten cases (Table  D2), depending on the stability of the cusp (determined by g xxx ), the direction from which the cusp approaches the fold (determined by the sign of ν[g](p 2 ) · µ[g](p 2 )), and k, the number of regular sheets of equilibria separating the fold and cusp. If k = 0, then one umbra intersect directly with a fold or cusp point (e.g. Figure 6). If k = 1, then two umbrae intersect, and if k ≥ 2 then none of the umbrae or folds intersect. The middle columns of Figures D4 and D5 show typical cases of these degeneracies. Note that no degeneracies involving the umbrae of a stable cusp exist, since stable cusps have no umbrae.

Aligned fold projection tangency
Opposed covering fold projection tangency Opposed non-covering fold projection tangency [g] onto the slow variables can intersect transversally in two ways: as a covering triple limit or as a non-covering triple limit (see Figure 8 b)). For brevity we write ν i := ν[g](p i ). In the covering case, all folds are opposed in the sense that their direction vectors span a convex cone covering all of R 2 . Therefore, the zero vector can be written as a linear combination of the direction vectors using only non-negative coefficients a i ≥ 0, not all zero: In the non-covering case the convex cone of the direction vectors does not cover R 2 , meaning that at least one coefficient has to be negative in order for the vector sum to be zero (see Figure 8 a)). Therefore, the two subcases are defined by the signs of the coefficients in (15) Non-covering triple limit if ± sign (a 1 , a 2 , a 3 ) = (+, +, −) for some choice of prefactor sign. Note that a higher codimension degeneracy will occur if a i = 0 for at least one i. Interactions of umbrae of the folds with other folds or umbrae give additional subclasses of triple limit points: these cases are detailed in Appendix D in Table D3. Note that it is not possible for all three fold umbrae to intersect. We summarise the discussion in this section with the following classification of codimension one bifurcations for the case of one fast and two slow variables, analogous to Proposition 2. such that precisely one of the sets D j,k [g, y] in Table 3 is non-empty, for precisely one y ∈ R 2 . At such a bifurcation, precisely one of the following occurs: (i) A loop or pair of hyperbolae appears in the fold projections at a fold tangency D 1,k .
• The function S(x, y) > 0 is smooth and positive on M × N .
• The function T (x, y) > 0 is smooth and positive on M × N .
Note that because we are only interested in equivalence of the singular systems, we allow independent re-parametrization of the fast and slow timescales. Note that T (x, y) is globally defined but only evaluated on C[g]. Clearly, if {g, h} is globally singularly equivalent to {g,h} then g is globally equivalent tog in the sense of (12). One can check that this is an equivalence relation -it is transitive and reflexive, and one can check it is symmetric by noting that if (17) holds then g(x, y) =S(x, y)g(X(x, y),Ỹ (y)) for all (x, y) because (x, y) ∈ C[g] if and only if Φ(x, y) ∈ C[g], and one can verify that: • The mapΦ(x, y) = (X(x, y),Ỹ (y)) is a diffeomorphism that is the inverse of Φ on M × N .
Note that singular trajectories are mapped onto each other by global singular equivalence as expressed in the following result.

Persistence under global singular equivalence
If system (1) is its own universal unfolding under global singular equivalence then we say the system is persistent. Clearly, in such a case the fast vector fields indexed by the slow variables must be persistent, but also we cannot have degeneracies of the slow system on the critical set. We expect (1) to be persistent under global singular equivalence if the fast subsystem is persistent and, in addition, the slow system has persistent behaviour on the critical set. For the case m = n = 1 we can make this statement more precise. We define the slow nullcline For any regular point p 0 = (x 0 , y 0 ) ∈ C reg [g] there will be a curve (X p 0 (y), y) ∈ C[g] with X p 0 (x 0 ) = y 0 such that g(X p 0 (y), y) = 0. Implicitly differentiating this gives for y close to y 0 . Then we can locally reduce (2) to an equation on the critical set of the formẏ = H p 0 (y) := h(X p 0 (y), y).
then H p 0 (y 0 ) = 0 is an equilibrium and its linear stability is determined via This highlights that the slow dynamics are essentially one-dimensional when restricted to C reg [g]. Before stating a result on persistence of fast-slow systems, we give some definitions. We define the restriction of the slow nullcline onto the critical manifold as and define the slow degenerate set E[g, h] as the union of two subsets, defined shortly: . The slow locally degenerate set is which occurs when the fast and the slow nullclines intersect tangentially. Note that this condition is equivalent to the determinant of the Jacobian of the full system being zero, and for p 0 ∈ C reg [g] this implies that H p 0 (y) = 0.
Recalling that π(p) is the projection onto the slow variables, we define the set of slow co-equilibria which we use to define the multiple slow equilibrium set We define the (mixed) fold-equilibrium multiple projection set as that is, the set of equilibria that share slow coordinate with at least one fold of the critical manifold. This finally allows us to define the fast-slow degenerate set as Theorem 1 below establishes that this set contains all degeneracies under global singular equivalence.
ingpersist Theorem 1 In the case m = n = 1, if g ∈ V f and h ∈ V s then (1)  and a y such that H p 0 (y) = H p 0 (y) = 0. But this is a non-hyperbolic equilibrium, and thus is not persistent to perturbation. Hence, for persistence we need For the "only if" part, we need to argue that there is no other way for (1) to be nonpersistent than if G[g, h] = ∅. Assume that G[g, h] = ∅. Then there is a neighbourhood of every singularity of g and equilibrium of {g, h} that has only one fold or equilibrium in the fast fibre. These are either quadratic folds of g or hyperbolic equilibria of {g, h}, both which are persistent under perturbation. Hence, {g, h} must be persistent.
Note that the global singular equivalence (17) does not depend on the nullcline N (x, y) away from the critical manifold. Bifurcation occurs if one of the assumptions in Theorem 1 is broken. Note that the assumptions that g ∈ V f implies the critical set does not intersect ∂M × N , and that h ∈ V s implies that the nullcline does not intersect M × ∂N ; more generally there will be additional persistence conditions that require persistent intersection with these boundaries.

Generic bifurcations in singular fast-slow systems
We can understand generic codimension one bifurcations of fast-slow systems (1) by examining the ways that the persistence conditions of Proposition 1 are violated. For m = 1 and n = 1, 2 this means that the codimension one bifurcations of the critical manifold are possible bifurcations under global singular equivalence. In addition, there are many ways that a change in the slow subsystem can lead to a bifurcation.
For m = n = 1 we define, as for the critical manifold, subsets of E 1 [g, h], E 2 [g, h] and M[g, h] containing all codimension one degeneracies of G[g, h], which are not only due to degeneracy of the fast subsystem, in Table 10. We further define subsets of these, which give codimension one bifurcation if all except for one of the subsets are empty, and if the nonempty subset is nonempty for only one slow coordinate y (Table 3). Equipped with the subsets in Table 11, Proposition 3 lists the codimension one degeneracies of (1) for m = n = 1.

Proposition 3
In the case m = n = 1, codimension one bifurcation of the fast-slow system (1) for {g, h} occurs due to exactly one of the following reasons, for exactly one slow coordinate y ∈ N .  for m = n = 1 whose union contains all codimension one bifurcations, not only due to bifurcation in the fast subsystem. Note that det(D{g, h}) = g x h y − g y h x is the Jacobian of the full system (1) and that the first subset is a local degeneracy Saddle-node: : |Ξ[g, h](p)| = 1 and det(D{g, h}(p)) = 0} Double slow equilibrium: Fold-equilibrium double projection set: Table 11. Degeneracies that lead to codimension one bifurcation of the singular fast-slow system (1) up to global singular equivalence. The sets F[g] and U[g] are the fold and the umbral sets and π(p) is the projection map. R(y) is the set of all equilibria sharing slow coordinate y. In saddle-node non-degeneracy condition, r and q are eigenvectors of the Jacobian of the full system and its adjoint respectively, and B = (v) There are exactly two hyperbolic equilibria that share the same slow coordinate y (i.e. there are exactly two points p 1 , p 2 ∈ E 2 [g, h] for which π(p 1 ) = π(p 2 ) = y, and Proof: Note that all degeneracies are contained in the set . Because of this, and because the defining conditions are independent, codimension one degeneracy will occur at a point that is in exactly one of those sets. Furthermore, bifurcation must occur for exactly one y since otherwise more than one equality constraint is imposed, raising the codimension. Case i) describes the only subset D 1 1 [g] of D 1 [g] containing codimension one degeneracies exclusively in D 1 [g], and therefore it produces codimension one degeneracy of {g, h}. The same is true for cases ii) and iii). Case iv) is codimension one since we impose just one equality condition and exclude higher codimension degeneracy with the non-degeneracy condition in Table 11. Case v) is codimension one since hyperbolic equilibria are persistent, and more than two hyperbolic equilibria sharing slow coordinate would impose more than one equality constraint. Case vi) is codimension one for the same reason.
Note that codimension two bifurcations may combine degeneracies in more than one of these sets.

Generic bifurcations of relaxation oscillations in singular fast-slow systems
ifurcation Not all bifurcations of the singular fast-slow system listed in Proposition 3 will lead to bifurcation of singular relaxation oscillations, as the degeneracy in the fast-slow system must interact with a limit cycle. We focus on bifurcation of singular relaxation oscillations and simple relaxation oscillations, a generic subclass due to [18]. Several cases of these bifurcations have been considered in the literature, see for example [35].

Singular relaxation oscillations
Consider a fast-slow system with m = 1 fast variables (1) in the singular limit = 0. A relaxation oscillation is a singular periodic trajectory γ : [a, b] → M × N (i.e. such that γ(b) = γ(a)) where the slow segments are in C att . If the oscillation consists of alternating stable slow segments on C att [g] up to non-degenerate folds, fast segments from these folds to their umbra, and satisfies certain other non-degeneracy conditions then we say it is a simple relaxation oscillation. These are called strongly common slow-fast cycles in [35] where it is shown that these singular trajectories will be shadowed by a stable periodic orbit for small enough . Guckenheimer stated a similar persistence theorem in [18]; in Section 5.2 we state and prove a version of it.
We say a continuous curve s k : [0, 1] → M × N is a slow segment of a singular trajectory if there is a continuous and monotonic increasing θ : [0, 1] → R such that s k (θ(s)) is a trajectory of (2). We say a slow segment s k has slow time duration T k > 0 if it can be parameterised by θ(s) = s/T k . If not, and if θ(0) = θ(1) then T k = 0, otherwise T k = ∞.
Up to equivalence of the fast segments joining the slow segment end-points, we define a relaxation oscillation in terms of its slow segments as a sequence of continuously parametrized slow segments s k : [0, 1] → M × N . We will assume that either A is a loop entirely within C att [g] or • s k (θ) ⊂ C att [g] for all θ ∈ (0, 1) • There is a trajectory φ(t) of the fast system such that α(φ(0)) = s k (1) and ω(φ(0)) = s k+1 (0) for k modulo d.
for each k, where ω(p) and α(p) are the omega and alpha limits of a point p respectively. This equivalence class of relaxation oscillations has more than one member if there is more than one fast segment joining two consecutive slow segments. We define the slow period P(A) of a relaxation oscillation A to be the total slow time duration of its slow segments. This is which may be infinite, where orbits in the equivalence class of A clearly all have the same slow period. We allow the possibility that A is a loop on C att [g] without jumps, in which case d = 1 and s 0 (1) = s 0 (0), or that the jumps are trivial and on C att [g]. Infinite slow period relaxation oscillation (P(A) = ∞) of a variety of types are covered by this definition. We define a simple relaxation oscillation (cf Guckenheimer [19,22]) as follows: Definition 2 (Simple relaxation oscillation (m = 1, n ≥ 1)) A relaxation oscillation A in (19) is simple if all of the following hold: i) The slow period P(A) is finite.
ii) The slow segments are on C att [g], except possibly the last point.

itm:sro4
iv) The slow segments are not tangent to either fold set or umbral set.
v) The singular return map local to s k (1) is well-defined with a hyperbolic equilibrium at s k (1).
itm:sro3 f:simplero Note that the assumption P (A) < ∞ implies that the slow segments do not limit to any equilibria of the slow flow. For one fast and one slow variable, our definition of a simple relaxation oscillation can be expressed in a simpler way: Definition 3 (Simple relaxation oscillation (m = n = 1)) A relaxation oscillation (19) with one fast and one slow variable (m = n = 1) is simple if i) The slow period P(A) is finite.

Persistence and bifurcation rsistbif11
We say that a simple relaxation oscillation undergoes bifurcation if the relaxation oscillation ceases to be simple under perturbation of the singular fast-slow system. If not, we say that the relaxation oscillation is persistent. Note that as we only consider fastslow systems on absorbing regions in R 2 , singular relaxation oscillations can bifurcate to either equilibrium points or other singular relaxation oscillations. The following proposition links bifurcation of simple relaxation oscillations to degeneracy in the fastslow system (cf [18,19,20]): Proposition 4 A simple relaxation oscillation A is persistent for n = m = 1 if the fast-slow system is persistent under global singular equivalence.
iessropers Proof: Assume {g, h} is persistent. Because the slow period is finite, no slow equilibrium can intersect a slow segment s k (θ) in a degenerate way; intersection for θ ∈ {0, 1} implies that M[g, h] = ∅ and intersection for θ ∈ (0, 1) implies E 1 [g, h] = ∅, since the flow must have the same direction on both sides of the equilibrium, which is not possible if the equilibrium is hyperbolic. Hence Condition i) is persistent. Condition ii) is persistent since s k (0) ∈ C att [g] implies that both s k−1 (1) and s k (0) are in F[g], which in turn would imply that D 3 [g] = ∅. Condition iii) is persistent since the fold set is nondegenerate everywhere, and since trivial relaxation oscillations s 0 (1) = s 0 (0) coincide with hyperbolic equilibria in the interior of C att [g].
The converse is not true, however: only some of the degeneracies in the fast-slow system will give rise to a bifurcation of a simple relaxation oscillation. The reasons are listed in the following Proposition and in Table 12. Examples are portrayed in Figures 11 and Figure 12.
Proposition 5 Bifurcation of a singular relaxation oscillation A for m = n = 1 due to codimension one bifurcation of the singular fast-slow system {f, g} occurs for exactly one of the following reasons, at exactly one point (x, y) ∈ A.
(i) There is saddle-node bifurcation of the slow subsystem in the interior of C att [g]: (x, y) ∈ E 1,1 [g, h, y].
(iii) There is a stable hysteresis of the critical manifold: (x, y) ∈ D 2,1 [g, y]. Table 12. Bifurcations of singular relaxation oscillations for m = n = 1 due to codimension one bifurcation of the fast-slow system. Definitions of the calligraphic sets are found in Table 3 and Table 11. It is understood that A(λ) is perturbed due to perturbed g(λ) = g(x, y, λ) and h(λ) = h(x, y, λ). The final column indicates whether the bifurcation is due to bifurcation of the critical set

Type: Example
Conditions for a simple singular relaxation oscillation A(λ) to bifurcate at λ = λ 0 Saddle node on invariant circle [12]: Figure 11 a,b,c) There exists a unique y ∈ N such that A lim (λ 0 ) = ∅ and E 1,

No
Singular Hopf: Figure 11 d,e,f) There exists a unique y ∈ N such that A lim (λ 0 ) = ∅ and M 1,

No
Singular homoclinic: Figure 11 g,h,i) There exists a unique y ∈ N such that A lim (λ 0 ) = ∅ and M 1,

No
Hyperbolic fold tangency: Figure
(v) There is an opposed double fold-umbra or umbra-umbra limit point of the critical manifold: (x, y) ∈ D 3,2 [g, y] or (x, y) ∈ D 3,4 [g, y]. Proof: We determine which of the codimension one degeneracies in Table 1 and Table 11 can cause bifurcation of singular relaxation oscillations by ruling out those that cannot, and by providing examples in the following section. First, two hyperbolic equilibria sharing y-coordinate E 2,1 [g, h, y], a non-interacting double limit point degeneracy D 3,5 [g, y] or D 3,6 [g, y] or an equilibrium sharing slow coordinate with a fold point but not intersecting it or its umbra M 1,5 [g, h, y] are ruled out because these degeneracies cannot break the simple property of limit cycles at codimension one.
Furthermore, codimension one bifurcation of limit cycles requires that a regular stable part of the critical manifold exists in a neighbourhood of the bifurcation point, or else the limit cycle does not generically pass through that point. This excludes elliptic fold tangency D 1,2 [g, y] and unstable hysteresis D 2,2 [g, y].
Moreover, a source equilibrium intersecting a fold M 1,2 [g, h, y] or a sink interacting with a fold umbra M 1,3 [g, h, y] are excluded since no relaxation oscillation can exist either at or in a neighbourhood of the bifurcation parameter at codimension one. Additionally, a "fold umbra -fold umbra" double limit point degeneracy cannot cause bifurcation at codimension one, since the umbra is generically on C reg [g], and not intersecting any equilibria (making the period finite). Therefore, there is no way for a simple relaxation oscillation to be lost at such a point.
The remaining codimension one bifurcations can break the simple property by violating one of the defining conditions: these are listed in Table 12. If a vector field is perturbed by a distinguished parameter λ ∈ R, a simple singular relaxation oscillation may cease to exist for some critical λ 0 where we assume the limit is from below. In such cases we define a limit relaxation oscillation A lim (λ 0 ) as the limit, in the Hausdorff distance, of a sequence of relaxation oscillations A(λ) parametrized by λ (21) eq:limitro The limit is well defined for simple relaxation oscillations but may be empty. The bifurcation of relaxation oscillations will unfold for the non-singular systems > 0 to give a variety of canards that will appear on a case to case basis: see for example [29,Chapter 8] and [35]. Outside a small (in ) range of parameters λ( ) near the critical λ c , many of the solutions will closely resemble those of the singular system. Hence, the bifurcations in Proposition 5 will give rise to a detectable qualitative change even for non-singular systems.
In this paper we have used bifurcation theory with distinguished parameters and singular equivalence to take some steps towards such a classification. Indeed, in [19] Guckenheimer gives the following list of codimension one degeneracies that we can relate to our classification: G1: A fast segment ends at a regular fold point. There are two cases depending on whether the slow flow approaches or leaves the fold near this point.
G2: A slow segment ends at a folded saddle.
G3: A fast segment encounters a saddle point.

G4:
There is a point of Hopf bifurcation at a fold.
G5: A slow segment ends at a cusp. G6: The reduced system has a quadratic umbral tangency between projections of fold and umbra.
The degeneracy G1 is a subcase of fold projection intersection for m = 1, n = 1, 2 and fold projection tangency for m = 1, n = 2. The degeneracy G2 appears when the slow flow is tangent to a fold line: this can occur for n ≥ 2. Degeneracy G3 can appear at a saddle for m ≥ 2 or at an unstable node for m = 1. Note that the formulation of G3 is slightly modified from [19]. Degeneracy G4 corresponds to a singular Hopf bifurcation, which we discussed in the context of m = n = 1. Degeneracy G5 corresponds to a hysteresis bifurcation for m = n = 1, and to a limit cycle hitting a cusp on the slow manifold for m = 1, n = 2. Finally, degeneracy G6 requires n ≥ 2. Guckenheimer states in [19] that the list is incomplete, and mentions the case that a slow segment ends at a folded node as an example. Degeneracy due to fold tangency of the critical set is missing from the list, since the slow variables were not regarded as distinguished parameters in [19]. We believe that Proposition 5 completes the list for m = n = 1.
For n = 2, the degeneracies involve tangencies of generic one-dimensional objects such as relaxation oscillations, fold lines and fold umbrae, as well as intersections of one-dimensional objects and generic zero-dimensional objects such as cusp points and equilibria in the slow subsystem. Some of these cases are in Guckenheimer's list. Note that degeneracies of the critical manifold do not cause codimension one bifurcation of relaxation oscillations since they occur at points, which do not generically intersect relaxation oscillations. They will be involved in bifurcation of invariant tori or more complex singular attractors or of relaxation oscillations at higher codimension however.

Relation to other singularity theory problems
There appears to be connections between global equivalence of critical sets with two distinguished parameters (for m = 1 and n = 2) and the equivalence of vector fields under projection to the slow plane. In particular, several hypothesised degeneracies under strong equivalence also appear as degeneracies of orthogonal projections of vector fields, see for instance [1,9,37,43,47] and references. Crucially, however, such equivalences do not produce degeneracies such as the fold tangency where the manifold structure is lost.
Our approach to bifurcation of the critical manifold in fast-slow systems uses a singularity theory approach with distinguished parameters from [17]. In the following we briefly indicate how this approach is related to singularity theory, catastrophe theory, the theory of constrained equations, and projections from manifolds to manifolds.
Much of singularity theory concerns the stability of zero sets of smooth functions g(x) = 0 under perturbation [46]. In the singularity theory approach to bifurcation theory of [17], there is a distinguished (bifurcation) parameter y that is not "mixed up" with the unfolding parameters. Hence, e.g. the quadratic fold g = x 2 + y is a codimension one degeneracy of g(x) = x 2 in singularity theory, but is codimension zero in bifurcation theory with one distinguished parameter y. For our interpretation, y is identified with the slow variables.
Catastrophe theory [2,40] classifies the changes to stationary points of potentials V (x), with ∇V (x) = g(x), by codimension of deformation. Although different equivalences and objects are studied in catastrophe theory and singularity theory, for one (fast) variable, the classification of local singularities of the critical set is the same. Constrained systems [42,26] correspond to singular fast-slow systems where equivalence of critical manifolds is defined by potential functions, as in catastrophe theory, together with a slow flow local to a point. Unlike our approach, there are no distinguished parameters and the unfolding parameters are identified with the slow variables. This means that some local bifurcations (notably fold tangencies) that are present for the distinguished parameter approach are missed, because "slow variables" never appear in powers higher than one in the local normal forms.
Singularity theory, catastrophe theory and constrained equations have been framed in terms of germs, which are local notions of functions. This means that global intersections of projections of singularities (such as double limit points which are important for bifurcations of relaxation oscillations) have not been widely studied in these contexts, some exceptions being [2,8,19,35].

Further perspectives
Persistence and codimension one bifurcation of the critical set for one fast (m = 1) and two slow variables (n = 2) remains to be proved. This requires a suitable equivalence, which should give rise to the degeneracies that we have listed, but possibly more.
A full investigation of bifurcations of singular relaxation oscillations for m = 1, n = 2 is outside the scope of this paper. Some specific examples have been studied by Guckenheimer [21,23,20] who outlined a scheme for investigation of bifurcation of solutions to singular fast-slow systems in [18].
The general case of two fast variables m = 2 is considerably more complicated as the vector fields cannot be written as gradients of potentials, and hence there can be other asymptotic behaviour than fixed points. If n = 2 then for generic asymptotic fast dynamics, the system will approach a critical set that is a union of all equilibria, periodic orbits and homoclinic/heteroclinic cycles of the fast system. The persistence of bifurcations on the critical set will depend on the number of slow variables. For n = 1 then we expect persistence precisely when (a) All singularities of equilibria within the critical set are quadratic folds or Hopf points. (b) All singularities of limit cycles within the critical set are one of saddle-nodes of limit cycles, saddle node on a periodic orbit, or homoclinic bifurcation. (c) The slow flow has generic intersection with umbrae of the singularities. For n = 2 we will get in addition generic local and global codimension two singularities at isolated points in the slow variables; this will include, for example, cusp points, Bogdanov Takens points and Bautin points at singular equilibria, and a wide variety of possible generic codimension two bifurcations of homoclinic orbits [10].
We have ignored the phenomena that arise when the scale separation is imperfect, that is for > 0. In that case, the fast and slow subsystems evolve at similar speeds close to singular points; this gives rise to canards and mixed mode oscillations [11]. Canard behaviour has been extensively studied, especially near regular values of the critical set see e.g. [6,29,18,45]. Canards for degenerate critical sets are discussed in [2,8,35], but we are unaware of any systematic treatment.
Appendix A. The quadratic curvature :curvature We define the scalar quadratic fold curvature (see Table 6)) at a fold point p in the direction of a fold of the critical manifold as is the Hessian of the slow subsystem, superscript T denotes transpose, ⊥ denotes perpendicular, and ∇ ⊥ y g denotes ∇ ⊥ y g scaled to unit length. For the remainder of this section we drop the dependency on p, such that e.g. g(p) is written just g. If ∇ ⊥ y g = 0, then K[g] is undefined. Note that this notion of scalar quadratic fold curvature only captures the quadratic curvature of the fold curve as a projection onto the slow (0, y 1 , y 2 ) plane.
In this section, we first motivate our definition of (A.1) and then offer an interpretation.
To motivate (A.1), we start with a quadratic fold at the origin and completely in the slow plane whose curvature in the slow plane reasonably should be considered to be along the y 1 -axis g(x, y 1 , y 2 ) = ξx 2 + c 1 y 1 + c 2 y 2 2 , where ξ, c 1 and c 2 are real constants.
Recall that the quadratic direction vector of a fold is given by ν[g] = g xx ∇ y g = (2ξc 1 , 0), so if ξc 1 > 0 then the fold is directed rightward, while if ξc 1 < 0 the direction is directed leftward. It makes sense to define the slow quadratic fold curvature of this fold as since then the curvature is independent of the magnitude of ξ, proportional to c 2 = 0, and inversely proportional to c 1 . The fold is convex in the direction of the fold if K[g](p) > 0 and concave if K[g](p) < 0; see Figure 5 for a graphical representation. We now consider a general quadratic polynomial function of a quadratic fold at the origin with terms of relevant order g(x, y 1 , y 2 ) = ξx 2 + ay 1 + by 2 + αy 2 1 + βy 2 2 + 2γy 1 y 2 + δxy 1 + ηxy 2 , (A.3) eq:generalfold where a, b, α, β, γ, δ and η are real constants unrelated to any quantities with the same names elsewhere in this text. The x term and constant term are missing because of the quadratic fold condition g = g x = 0, and higher order terms are not present since they should not matter for the quadratic curvature.
We can now define the quadratic fold curvature vector as The quadratic fold curvature vector points in the direction of the fold if the fold is convex, and against the direction if it is concave (see Figure 5). Some special cases of Equation (A.1) are insightful. First, if the fold lies locally in the slow plane (that is, if the coefficients of the xy 1 and xy 2 terms in (A.3) are zero), then a positive definite matrix sign (g xx )D 2 y (g) implies that the fold is convex and a negative definite sign (g xx )D 2 y (g) implies that it is concave. However, if sign (g xx )D 2 y (g) is indefinite or has a zero eigenvalue, then the sign of K[g] can be either positive, negative or zero depending on the direction of the fold.
On the other hand, if the xy 1 and xy 2 terms are present in (A.3), then the second term in (A.1) only serves to reduce convexity (or equivalently increase concavity). The extent to which convexity is reduced depends quadratically on the component of the gradient of ∇ y g x perpendicular to the gradient, and inversely on the magnitude of the gradient and the magnitude of the curvature in the x direction.

Appendix A.1. Persistent subcases of the fold projection tangency
We classify the fold projection tangency degeneracy (degeneracy subset D 4 [g] for m = 1 fast and n = 2 slow variables) into a number of qualitatively different subcases. The subcases are separated by the scalar quadratic fold curvatures at the points of degeneracy, whether the umbrae interact with each other or a fold, and in the case of fold-umbra degeneracy, whether the dominant curvature belongs to the fold curve with the largest x-coordinate.
Aligned folds generate four distinct subcases. In analogy with the situation for m = n = 1, we have fold-fold, fold-umbra, and non-interacting fold cases. But the foldumbra has two subcases, depending on whether the fold with interacting umbra also has the greatest scalar quadratic fold curvature K[g]. Hence, there are four subcases.
Opposed folds have six distinct subcases, three for each case that either the sum of curvatures is positive (net convex) or negative (net concave). Opposed fold projection tangency does not have the two fold-umbra subcases of aligned fold projection tangency, since the dominant x-values are reversed by one half rotation of the slow variables. Hence there are six such subcases, and ten subcases in total.
The nonpersistence condition (at codimension one) for fold projection tangency degeneracy at points p 1 ∈ D 4 [g] and p 2 ∈ Π(p 1 ) ∩ D 4 [g] is that is, the folds must have distinct quadratic curvature vectors (see Appendix A). However, to distinguish subcases of degeneracy we will later use the scalar quadratic fold curvature K[g] and information about whether folds are aligned or opposed.
We now separate subcases depending on whether folds are aligned or opposed. If ν[g](p 1 ) · ν[g](p 2 ) > 0 then the folds are aligned, and qualitatively indistinguishable unless the umbra of the fold with the larger x component hits the other fold. If the umbra U [g](p 1 ) of an aligned fold at a point p 1 hits another fold at a point p 2 and the scalar quadratic fold curvatures satisfy K[g](p 1 ) > K[g](p 2 ), then the degeneracy is called umbra dominant. Otherwise if K[g](p 1 ) < K[g](p 2 ), then the degeneracy is called fold dominant. Therefore, we get the four subcases illustrated in Figure D2 and tabulated in Table D1. (Figure 7 a,b,c) shows a slow projection sketch of aligned fold projection tangency).
If, on the other hand ν[g](p 1 )·ν[g](p 2 ) < 0, then the folds are opposed and therefore distinguishable at bifurcation. (Figure 7 d,e,f) shows a slow projection sketch of opposed fold projection tangency).
Assuming that folds are opposed, the defining condition for covering opposed fold projection tangency is that is, if at least one fold is concave, and the magnitude of the curvature of the convex fold is smaller than that of the convex curve (see Figure 7 d,e,f)). If, on the other hand then we have covering opposed fold projection tangency (see Figure 7 g,h,i) ).
Opposed fold curves do not have the umbra-dominant and fold-dominant subcases that aligned fold curves do, since a rotation of the slow variables by 180 degrees turns one such case into the other. Therefore, there are four aligned cases and six opposed cases, giving in total the ten subcases in Table D1.
Appendix B. Definition of the cusp direction vector pdirection We define the direction of a cubic cusp (see Table 6)) to be For the remainder of this section we do not explicitly write out the dependence on p, so that for instance g(p) is written g. For the example in this section we assume that p = (0, 0, 0). As for the definition of scalar quadratic fold curvature, we motivate the definition starting from a normal form of the cubic cusp g(x, y 1 , y 1 ) = ξx 3 + c 1 xy 2 + c 2 y 1 , where ξ, c 1 and c 2 are real constants. In this case we naturally let the cusp direction vector be where (0, −c 2 ) is a vector perpendicular to the gradient (c 2 , 0). It makes sense for the direction of the cusp to be along the gradient perpendicular, since it is parallel to the pair of fold curves which emanate from the cusp. Next, we consider a more general expression g(x, y 1 , y 2 ) = ξx 3 + αxy 1 + dxy 2 + ay 1 + by 2 , (B.4) eq:generalcuspn where ξ, a, b, α and β are coefficients, and where we do not keep terms of intermediate orders xy 2 1 or x 2 y 2 , since visually they do not seem to alter the direction or sharpness of the cusp. We then rotate the slow subsystem as to make the gradient in the old coordinates (a, b) directed along the positive y 1 axis. We accomplish this with a rotation R mapping the new coordinatesŷ 1 ,ŷ 2 to the old ones such that in the new coordinates Equation (B.4) becomes By reading off the coefficients of the x 3 , y 2 x and y 1 terms, we find that the cubic cusp direction vector in the new coordinates is Rotating this vector back to the original coordinates we get that which we recognise can be written (up to constant scaling) where ∇ ⊥ y g = ∇ ⊥ y g/|∇ ⊥ y g|. Hence, the cusp is always perpendicular to the slow gradient, with magnitude inversely proportional to the projection of ∇ y g x onto the gradient perpendicular. As a consequence, the magnitude of the cubic cusp direction vector blows up (becomes undefined) if ∇ y g x is parallel to the gradient (its perpendicular component vanishes).

Appendix C. The quantity W [g]
app:W To construct the cusp quantity W [g](p) ( Table 6) at a point p = (x, y 1 , y 2 ), we start from a modified normal form for the cusp tangency bifurcation (see Table 4) g(x, y 1 , y 2 ) = δ 1 x 3 + δ 2 xy 2 1 + δ 3 x 2 y 1 + xλ + y 2 , where λ is an unfolding parameter. The sign of the cusp quantity W [g] = 6δ 1 δ 2 −2δ 3 = 0 gives the type of bifurcation (W [g] < 0 gives beaks and W [g] gives lips), see [38]. Our aim is to express W [g] for a general function equivalent to g. In this section we do not explicitly write out the dependence on p, which for the normal form is p = (0, 0, 0), so that e.g. W [g](p) is written W [g]. The parameter δ 3 does not appear in the normal form in Table 4 because there are only four subcases of cusp tangency, which can be represented with just δ 1 and δ 2 . However, all of δ 1 , δ 2 and δ 3 are needed in the nondegeneracy condition W [g] = 0 and for separating different subcases.
Rotating the coordinate system to have gradient in the positive y 1 direction by the transformation gives that the coefficient of the xŷ 2 2 term is and the coefficient of the x 2ŷ 2 term is Note that the gradients and the slow Hessian D 2 y (g x ) are expressed in the original slow coordinates (y 1 , y 2 ). Combining the expressions for δ 1 , δ 2 and δ 3 , we have that  Table 9). Each row shows unfolding with a bifurcation parameter λ. Bifurcation occurs as λ = λ 0 . Solid/dashed black lines show stable/unstable sheets of the critical set and red lines show the image of the fold under the umbral map ureCM121D2 where ∇ ⊥ y g = ∇ ⊥ y g/|∇ ⊥ y g| is the unit length slow gradient perpendicular. We repeat that if W [g] < 0 then the cusp tangency is of beaks type, whereas if W [g] > 0 it is of lips type.
(Color online) Examples of codimension one aligned fold projection tangency bifurcation (see Table D1). Each row shows unfolding with a bifucation parameter λ. Bifurcation occurs as λ = λ 0 . Solid/dashed black lines show the stable/unstable sheets of the critical set, red lines show the image of the fold under the umbral map and blue lines indicate tangency of projections of fold sheets. As a visual aid, the number of sheets of the critical set in a neighbourhood of the bifurcation is shown to the right, to be viewed as a projection onto the slow variables  Table D1. Each row shows unfolding with a bifucation parameter λ. Bifurcation occurs as λ = λ 0 . Solid/dashed black lines show the stable/unstable sheets of the critical set, red lines show the image of the fold under the umbral map and blue lines indicate tangency of projections of fold sheets. As a visual aid, the number of sheets of the critical set in a neighbourhood of the bifurcation is shown to the right, to be viewed as a projection onto the slow variables  Table D2). Each row shows unfolding with a bifucation parameter λ. Bifurcation occurs as λ = λ 0 . Solid/dashed black lines show the stable/unstable critical set, red lines show umbrae, and blue lines indicate intersection of projections of cusps and fold onto the slow plane. As a visual aid, the number of sheets of the critical set in a neighbourhood of the bifurcation is shown to the right, to be viewed as a projection onto the slow variables  Table D2). Each row shows unfolding with a bifucation parameter λ. Bifurcation occurs as λ = λ 0 . Solid/dashed black lines show the stable/unstable critical set, red lines show umbrae, and blue lines indicate intersection of projections of cusps and fold onto the slow plane. As a visual aid, the number of sheets of the critical set in a neighbourhood of the bifurcation is shown to the right, to be viewed as a projection onto the slow variables  Table D3). Each row shows unfolding with a bifucation parameter λ. Bifurcation occurs as λ = λ 0 . Solid/dashed black curves show the stable/unstable critical set, red curves show umbrae, and blue curves show intersections of at least the lower two folds, but possibly all three folds. As a visual aid, the number of sheets of the critical set in a neighbourhood of the bifurcation is shown to the right, to be viewed as a projection onto the slow variables ureCM121D6 Appendix E. Examples of bifurcations of relaxation oscillations for one fast and one slow variable plesection In this section we present the equations for the example fast-slow systems for m = n = 1, showing bifurcations of relaxation oscillations due to the critical manifold in Figure 12, how they were constructed, and how the figures were produced.
Appendix E.1. Bifurcation of relaxation oscillation due to hyperbolic fold tangency We seek fast and slow subsystems g(x, y) and h(x, y) such that (1) displays bifurcation of singular relaxation oscillations due to hyperbolic fold tangency bifurcation (Figure 12  a,b,c)).
The fast subsystem g(x, y) is written as a perturbed product of a hysteresis curve and a circle: g hyst (x, y) = x 3 − 2x + y g circ (x, y) = (x − λ) 2 + (y − y c ) 2 − R 2 g(x, y) = −(g hyst (x, y)g circ (x, y) + λx + q).  (Table D3), to be viewed almost as a projection onto the (x, y 1 ) plane. Each row shows unfolding with a bifucation parameter λ. Bifurcation occurs as λ = λ 0 . Solid/dashed black horizontal curves show stable/unstable sheets of the critical set, and dashed vertical lines indicate coordinates where slow projections of folds intersect transversally. Red, blue and magenta lines are representations of three fold lines (the row a,b,c) illustrates this representaion for the d,e,f) row). Bifurcation occurs as the slow projections of all three folds intersect. Each case can be either covering or or non-covering, as shown in Fig. D6

D6Sideview
We fix x 1 = −1,x 3 = 1/2, x 4 = 5/4 and leave a and x 2 as free parameters. We take g to be the primitive function of g x with zero constant term (default of Matlab's int command) such that and g(0, 0) = 0. We solve the linear pair of equations g(x 1 , 1) = g(x 3 , 1) = 0 for a and x 2 , giving x 2 = −13/40 and a = 640/49. Then we add a bifurcation parameter λ breaking the degeneracy, giving In Figure 12 g,h,i) λ ∈ [−0.1, 0.1]. Finally, we reverse the sign of x, such that The slow subsystem is set to be positive above the constant nullcline x = x nc = 1.5 and negative below such that h(x, y) = x − x nc .
In Figure 12 g,h,i) we solve (1) with Matlab's stiff solver ode23s for 1000 time units, starting from initial conditions x = y = 0 and with scale separation parameter = 0.01.

Appendix E.4. Opposed double limit point bifurcation of relaxation oscillations
We seek fast and slow subsystems g(x, y) and h(x, y) such that (1) displays bifurcation of singular relaxation oscillations due to opposed double limit point bifurcation ( Figure 12 j,k,l)).