Singularly Perturbed Boundary-Equilibrium Bifurcations

Boundary equilibria bifurcation (BEB) arises in piecewise-smooth systems when an equilibrium collides with a discontinuity set under parameter variation. Singularly perturbed BEB refers to a bifurcation arising in singular perturbation problems which limit as some $\epsilon \to 0$ to piecewise-smooth (PWS) systems which undergo a BEB. This work completes a classification for codimension-1 singularly perturbed BEB in the plane initiated by the present authors in [19], using a combination of tools from PWS theory, geometric singular perturbation theory (GSPT) and a method of geometric desingularization known as blow-up. After deriving a local normal form capable of generating all 12 singularly perturbed BEBs, we describe the unfolding in each case. Detailed quantitative results on saddle-node, Andronov-Hopf, homoclinic and codimension-2 Bogdanov-Takens bifurcations involved in the unfoldings and classification are presented. Each bifurcation is singular in the sense that it occurs within a domain which shrinks to zero as $\epsilon \to 0$ at a rate determined by the rate at which the system loses smoothness. Detailed asymptotics for a distinguished homoclinic connection which forms the boundary between two singularly perturbed BEBs in parameter space are also given. Finally, we describe the explosive onset of oscillations arising in the unfolding of a particular singularly perturbed boundary-node (BN) bifurcation. We prove the existence of the oscillations as perturbations of PWS cycles, and derive a growth rate which is polynomial in $\epsilon$ and dependent on the rate at which the system loses smoothness. For all the results presented herein, corresponding results for regularized PWS systems are obtained via the limit $\epsilon \to 0$.


Introduction
This manuscript concerns the unfolding of singularities in planar singular perturbation problems which limit to piecewise-smooth (PWS) systems. The underlying PWS system is assumed to have a smooth codimension-1 discontinuity set, or switching manifold Σ ⊂ R 2 , which has an isolated boundary equilibrium (BE). BEs are PWS singularities which unfold generically in a codimension-1 bifurcation known as a boundary equilibrium bifurcation (BEB), whereby an isolated equilibrium collides with the switching manifold Σ under parameter variation.
A first classification of planar BE singularities appeared in Filippov's seminal work on discontinuous PWS systems [10]. Here it was shown that generically, there are 8 topologically distinct classes of BE singularities, comprised of 2 boundary-saddle (BS), 2 boundary-focus (BF) and 4 boundary-node (BN) singularities. A treatment of the unfolding of these singularities came later in [32], where the authors identify 10 topologically distinct unfoldings, and provide 'prototype systems' for each. Subsequently in [14], two more unfoldings were identified, bringing the total count to 12. Here we present a single prototype system capable of generating all 12 unfoldings, and a completeness theorem [14,Theorem 2] ruling out the possibility of additional missing cases. Explicit local normal forms (as opposed to 'prototypes') for a large number of BE singularities have been derived in [6], but their work did not treat the unfolding via BEB.
The notion of singularly perturbed BEB was developed more recently in [19], for the analysis of smooth singular perturbation problems limiting to PWS systems with a BF bifurcation. The motivation to study smooth perturbations of PWS systems arises from the observation that PWS systems often serve as approximations for smooth dynamical systems with abrupt transitions in phase space. Hence, it is natural to consider a class of smooth singular perturbation problems, which limit to PWS systems that are discontinuous along a switching manifold Σ as a perturbation parameter → 0. Abrupt dynamical transitions in such systems occur within an −dependent neighbourhood U ⊂ R 2 about Σ known as the switching layer, which satisfies U → Σ as → 0. It is important to note that singular perturbation problems in this class can arise either (i) naturally, or (ii) by a process of regularization whereby a modeller 'smooths out' discontinuity in a PWS system. In the former case, the problem is given as a smooth singular perturbation problem with a PWS singular limit; see e.g. [18,27] for applications of this kind. In the latter case, the PWS system is given, and the modeller introduces a method of regularization based on the characteristics of the problem at hand; many examples of this kind can be found in [16]. In both cases, analytical techniques from PWS systems and Geometric Singular Perturbation Theory (GSPT) [21,30,38], in combination with a method of geometric desingularization known as the blow-up method [9,29], provide a powerful analytical framework; see e.g. [4,5,13,16,24,25,26,22,27,34,36,37]. It is worthy to note that the authors in [7] consider a large number of BEBs in the context of regularized PWS systems, however the degeneracy associated with the BE singularity is not fully resolved.
The present manuscript provides a classification and detailed dynamical study of singularly perturbed BEBs in the plane. The work can be seen as a continuation of recent work in [19], see also the PhD thesis [17], where the analysis was restricted to a subset of singularly perturbed BF bifurcations, treating 3 of the total 12 BE unfoldings in detail, and successfully resolving the degeneracies associated with these cases. This manuscript aims to complete the project, by providing a 'complete' description for all 12 unfoldings. Similarly to [17,19], emphasis is placed on understanding the smooth dynamics with 0 < 1. This allows for the treatment of problems arising either naturally or via regularization simultaneously, since the corresponding results for (regularized) PWS systems are easily obtained upon taking the non-smooth singular limit → 0. First, we show that the C r≥1 local normal form derived for singularly perturbed BF bifurcations in [19] is in fact sufficient to generate all 12 unfoldings. A corresponding PWS local normal form is obtained from this expression in the non-smooth singular limit → 0.
We then study all 12 unfoldings for 0 < 1. As found in the analysis of singularly perturbed BF bifurcations in [19], each unfolding typically involves singular bifurcations, in some cases codimension-2, occurring within an −dependent domain which shrinks to zero as → 0 at a rate which can be quantified explicitly in terms of rate at which the system loses smoothness. We present 2-parameter bifurcation diagrams for a desingularized system with = 0 which determines the qualitative dynamics for 0 < 1. It is worthy to note that within the class of smooth monotonic regularizations considered, the dynamics are shown to be qualitatively determined by the underlying PWS problem, i.e. the bifurcation structure is qualitatively independent of the choice of regularization, and determined by the type of PWS unfolding in the limit → 0. It is shown how the choice of regularization does, however, effect the dynamics quantitatively, particularly due to its determination of the rate at which the system loses smoothness as → 0.
Following an analysis of the unfoldings, we present new results on the asymptotics of distinguished homoclinic solutions corresponding to boundaries between singularly perturbed BF 1 and BF 2 bifurcations. Finally, special attention is devoted to the singularly perturbed BN 3 bifurcation, which provides the necessary local mechanisms for the onset of relaxation-type oscillations. This was first observed in [27] in the context of substrate-depletion oscillations. Whereas emphasis there was on the existence of relaxation-type oscillations for the specific model, we will in the present work identify and describe the explosive onset of oscillations in the case of generic singularly perturbed BN 3 bifurcations. Similarly to the canard explosion phenomena known to occur in classical slow-fast systems [9,30,29], we will show that limit cycles in the singularly perturbed BN 3 bifurcation perturb from a continuous family of singular cycles, i.e. closed concatenated orbits having segments along a critical manifold. However, the local mechanism for the onset of explosive dynamics differs from that of classical canard explosion, and functions without the need for canard solutions. In contrast to the exponential growth rate associated with classical canard explosion, we show that the growth rate of the cycles arising in singularly perturbed BN 3 is polynomial in . We quantify this growth rate in terms of properties of the regularization.
The manuscript is structured as follows: Basic definitions and setup are introduced in Section 2. The C r≥1 local normal form capable of generating all 12 unfoldings for 0 < 1, as well as the resulting smooth and PWS classifications are also given in Section 2. Main results are presented in Section 3. Specifically, the blow-up analysis is outlined in Section 3.1, unfoldings and corresponding 2-parameter bifurcation diagrams are presented in Section 3.2, asymptotic results on boundary separatrices are presented in Section 3.3, and results on the singularly perturbed BN 3 explosion are given in Section 3.4. Main results on the unfolding and boundary separatrices are proved in Section 4, and a proof for the results pertaining to BN explosion are presented in Section 5. Finally in Section 6, we conclude and summarise our findings.

Acknowledgement
The first author acknowledges partial funding from the SFB/TRR 109 Discretization and Geometry in Dynamics grant, and together with the third author, partial funding from the ARC Discovery Project grant DP180103022.
The second author is grateful for his discussions with Peter Szmolyan on the explosive growth of limit cycles for the boundary node case in the context of the substratedepletion oscillator.
where the vector fields Z ± : R 2 × I → R 2 are smooth.
Assumption 2. The smooth 'regularization function' φ : R → R satisfies the monotonicity condition ∂φ(s) ∂s > 0, for all s ∈ R and, moreover, It follows from Assumption 1 and the form of the (non-uniform) limit in (3) that system (1) is (generically) PWS in the singular limit → 0. In particular, the limiting system is PWS and (generically) discontinuous along the switching manifold Remark 2.1. The more general scenario where Σ = {z ∈ R 2 : f Σ (z) = 0} for any smooth function f Σ : R 2 → R such that Df Σ | Σ = (0, 0), can easily be incorporated into the preceding formalism by replacing y with f Σ (z) in system (1) and adjusting Assumptions 1 and 2 accordingly. Since we restrict to a local analysis throughout, we may assume that f Σ (z) = y without loss of generality.
Notice that system (4) can be 'regularized' via (2) with p = φ(y −1 ). Hence, system (1) can be viewed as either of the following: • A smooth singularly perturbed system with a PWS singular limit; • A smooth regularization of the PWS system (4).
In this work we shall prioritise the former interpretation, since (i) this case is treated in less detail so far in the literature, and (ii) findings pertinent to the latter case can be immediately inferred from the dynamics of the nearby smooth system upon taking the limit → 0.
We impose one more technical assumption, which restricts the class of regularization functions φ: Assumption 3. The regularization function φ(s) has algebraic decay as s → ±∞, i.e. there exist k ± ∈ N + and smooth functions φ ± : [0, ∞] → [0, ∞) such that Assumption 3 restricts to the class of regularization functions with algebraic decay toward 0, 1, and is natural in the context of general systems (1) with analytic or sufficiently smooth right-hand-side. Specifically, it follows that both mappings u → φ(±u −1 ) for u > 0 have well-defined Taylor expansions at u = 0, which are each nondegenerate in the sense that there are leading nonzero terms (1−u k + β + and u k − β − , respectively) at order k ± , respectively. Note this assumption precludes regularizations like φ(s) = tanh(s) or φ(s) = e s /(1 + e s ), which have exponential decay toward 0, 1 and thus k = ∞. We omit the rigorous treatment of these cases, but refer to [18,22] for details on how to handle non-algebraic asymptotics using an adaptation of the blow-up method.
Remark 2.2. Regularisation functions φ which satisfy Assumptions 2 and 3 can be analytic, and should be distinguished from the well known class of non-analytic Sotomayor-Teixeira (ST) regularizations. In particular, the regularizations considered herein do not feature an artificial cutoff at the boundary to the switching layer.

PWS preliminaries.
It follows from our assumptions that the PWS system (4) is Filippov-type [10]. In particular, sliding and crossing regions of Σ can be determined in accordance with their usual definitions. Definition 2.3. Given system (4) and a point p ∈ Σ. Then p ∈ Σ is called a crossing (sliding) point if the quantity is positive (negative), where Z ± f (·) = ∇f (·), Z + (·, α) denotes a Lie derivative. We denote the set of crossing (sliding) points by Σ cr (Σ sl ).
It follows from our assumptions on φ that the sliding/Filippov vector field described as a convex combination in [10] can be written as where [Z + , Z − ] denotes a Lie bracket. If f Σ (x, y) = y as in (5), then the sliding/Filippov vector field is given in the x−coordinate chart by Figure 1. The 8 BE singularities arising in Filippov's topological classification [10]. The switching manifold Σ = {y = 0} is shown in green, with sliding/crossing submanifolds in bold/dashed respectively. We adopt the labelling convention in [14] with S, n, N, F denoting saddle, stable node, unstable node, focus respectively, and I/O denoting inward/outward flow along Σ sl . We have chosen an orientation such that the Σ sl always lies to the left.
where det(Z + (x, 0)|Z − (x, 0)) denotes the determinant of the 2×2 matrix with columns Z + (x, 0), Z − (x, 0). Sliding trajectories can leave Σ sl at a point of tangency with either vector field Z ± . Depending on the order of the tangency, such a point may also separate sliding and crossing regions of Σ sl . The following definition characterises the least degenerate case, i.e. quadratic tangency with either Z ± . Definition 2.4. Given system (4) and a point F ∈ Σ. Then F is a fold point if either It remains to review the notion of BE singularities and BEB. BE singularities arise when one or both of the vector fields Z ± (z be , α be ) = (0, 0) T for some z be ∈ Σ and parameter value α = α be . We consider the least degenerate case in which z be ∈ Σ is a hyperbolic equilibrium of Z + (·, α be ), and Z − is locally transverse to Σ. Filippov showed in [10], see also [14], that there are 8 topologically distinct cases depending on: • The type of equilibrium (focus, node or saddle); • The orientation of the sliding dynamics (towards or away from z be ); • In the case that z be is a node of Z + (·, α be ), its asymptotic stability (stable or unstable); see Figure 1. As described in [14], the 8 cases can be neatly categorised if we let S, n, N, F denote 'saddle', 'stable node', 'unstable node', 'focus' respectively, and let I/O define inward/outward sliding flow (i.e. towards or away from z be ). Then the possible cases are: SO, SI, nO, nI, NO, NI, FO and FI. BE singularities unfold generically under parameter variation in a BEB. Below we provide a formal definition for BEB in general PWS systems (4).  Figure 1) in a BN 3 bifurcation.
The topological classification in [14] shows that generically, the 8 BE bifurcations in Figure 1 unfold in 12 topologically distinct BEBs. Specifically, there are 5 BF bifurcations, 4 BN bifurcations and 3 BS bifurcations. We shall label these by BF i , BN i and BS i for i ∈ {1, . . . , 5}, i ∈ {1, . . . , 4} or i ∈ {1, 2, 3} respectively, in accordance with the notational conventions from [32]. The two unidentified BN bifurcations later described in [14] will be denoted BN 3 and BN 4 . The BN 3 unfolding is of particular interest in this work and shown in Figure 2. The fact that there may be more than one unfolding per BE is a consequence of the relative positioning of separatrices; topologically distinct BEBs can be separated by so-called 'double separatrices' [14] which connect equilibria and points of tangency on Σ. The role of separatrices is also discussed in e.g. [3,12].
Remark 2.6. The determinant condition in (12) ensures that the equilibrium of Z + (·, α) collides with Σ transversally under variation in α. To see this, notice that in the extended (z, α)−space, the vector T bf := (∇Z + 1 ∧ ∇Z + 2 )| (z bf ,α bf ) is tangent to the curve defined implicitly by Z + (z, α) = (0, 0) T . The stated determinant condition follows by the requirement that T bf has a non-zero y−component.
Finally, we introduce the notion of singularly perturbed BEB.
Definition 2.7. We say that system (1) under Assumptions 1 and 2 has a singularly perturbed BEB if the PWS system (4) obtained in the singular limit → 0 has a BEB. Notions of singularly perturbed BF bifurcation, singularly perturbed BN bifurcation and singularly perturbed BS bifurcation are similarly defined.
By definition, the existence of 12 BEBs implies the existence of 12 singularly perturbed BEBs.

Normal form and classification.
We show that the normal form derived for singularly perturbed BF bifurcations in [19] generalises to a single normal form capable of generating all 12 singularly perturbed BEBs.
Theorem 2.8. Consider system (1) under Assumptions 1, 2, and assume that the PWS system (4) obtained in the limit → 0 has a BEB of type BF, BN or BS at z bf ∈ Σ when α = α bf . Then there exists constants such that system (1) can be smoothly transformed, up to a reversal of orientation, into the local normal form where θ i (x, y, µ), i = 1, 2 are real-valued smooth functions such that and µ is a new bifurcation parameter related to α via µ = g(α), where g : I α → R is a smooth function such that g(α bf ) = 0 and g (α bf ) = 0.
The PWS system obtained from (13) in the limit → 0 + has a BEB at the origin for µ = 0, and a Filippov/sliding vector field given by Proof. The proof is similar to derivation of the normal form for singularly perturbed BF bifurcations presented in [19, p.38], and deferred to Appendix A for brevity.
Remark 2.9. Note the qualifier "up to a reversal of orientation" in Theorem 2.8. Orientation should be reversed if the vector field component Z − 2 (z, α) in system (1) satisfies Z − 2 (0, 0) < 0. A classification of singularly perturbed BEBs with 0 < 1 can be given via the classification of the underlying PWS system for → 0. This approach is similar to the classification of singularities in slow-fast systems in terms of their 'singular imprint' for = 0.
Similarly to the prototype system given in [14], the PWS normal form (14) can be used to generate all 12 BEBs by a suitable restriction of parameters in the PWS normal form (14). Each unfolding can be identified with an open region in (τ, δ, γ)−parameter space determined by the quantities τ, δ, ∆ := τ 2 − 4δ, γ.
Double-separatrices which connect a visible fold point with an equilibrium on Σ sl also play a role in separating regions corresponding to BF 1,2 , and regions corresponding to BS 1,2 . Here the distinction lies in whether or not the separatrix emanating from the  Table 1. Classification for the singularly perturbed BEBs generated by the local normal form (13), given in terms of a PWS classification for the PWS local normal form (14) obtained in the singular limit → 0. The classification is equivalent to the PWS classification in [14, Table  1]. Here ± denotes the sign of the corresponding quantity.
fold point connects to the regionΣ sl ⊂ Σ sl which is bounded between the fold point and the equilibrium. The resulting classification, which is equivalent to that in [14, Table 1], is given in Table 1.

Main results
In this section we present our main results. We begin in Section 3.1 with an outline of the sequence of blow-up transformations necessary to resolve all degeneracy associated with singularly perturbed BEB in system (13). This allows for the identification of a desingularized system governing the unfolding of the singularity. In Section 3.2, we present the unfolding for all 12 singularly perturbed BEBs. In Section 3.3 we present results on the asymptotics of a homoclinic double-separatrix which separates singularly perturbed BF 1,2 bifurcations. The BS 1,2 boundary is also discussed. Finally in Section 3.4, we present results on an observed 'explosion' in the case of singularly perturbed BN 3 bifurcations.
3.1. Resolution via blow-up. We describe the blow-up analysis used to resolve degeneracy in system (13) due to either (i) the loss of smoothness along Σ, or (ii) the loss of hyperbolicity at fixed points. The sequence of blow-up transformations is the same as in [19], so we restrict ourselves here to an overview. System (13) loses smoothness along Σ in the singular limit → 0. To describe this, we follow [23,26] and others and consider extended system {(x , y ) = X(x, y, µ, ), = 0} , with respect to a fast time, recall (13). For this system Σ × {0} is a set of equilibria with a loss of smoothness. We gain smoothness via a homogeneous cylindrical blow-up  The loss of smoothness along Σ × {0} has been resolved, but a degenerate point Q (also in orange) stemming from the tangency point persists. An attracting critical manifold S terminating at Q is identified on the cylinder, and shown here in blue. The local projective coordinates (x, r 1 , 1 ) defined in (17) and centered at Q are also shown.
transformation of the form which replaces Σ × {0} by the cylinder {r = 0} × R × S 1 , see Figure 3. The subspace {r = 0} corresponding to the blow-up cylinder is invariant. After a suitable desingularization amounting to division by¯ , the dynamics within {r = 0} are governed by a slow-fast system with a normally hyperbolic and attracting critical manifold, denoted S in Figure 3(b) [4,5,23,26,33,34]. Moreover, there is a reduced flow on S which is topologically conjugate to the sliding/Filippov flow induced by (15).
The critical manifold S terminates tangentially to the fast flow at a degenerate point Q ∈ {r = = 0}, which is also a point of tangency with the outer dynamics induced by the vector field X + within { = 0}; see again Figure 3(b). Choosing local coordinates of the form with x unchanged, this degeneracy is identified as a fully nonhyperbolic (i.e. no eigenvalues with non-zero real part) equilibrium at (x, r 1 , 1 ) = (0, 0, 0).
The point Q is degenerate for all µ ∈ R, however degeneracy stemming from the presence of the tangency is resolved via the weighted spherical blow-up where k := k + ∈ N + is the decay exponent associated with the regularization function φ, see equation (6). After another desingularization (division by ρ k(1+k) ), nontrivial dynamics are identified within the invariant subspace {ρ = 0} corresponding to the blow-up sphere {ρ = 0} × S 2 . The critical manifold S now connects to a partially hyperbolic and (partially) attracting (i.e. there is an eigenvalue with negative real part) equilibrium p a contained within the intersection of the blow-up cylinder and blow-up sphere, see Figure 4(a). An attracting center manifold W ∈ {ρ = 0} emanates from p a , thereby 'extending' S. Whether or not equilibria are also identified along the intersection of the blow-up sphere with { = 0} depends on whether the corresponding BEB is type BS, BN or BF, as well as on the sign of µ; see Figure 4(b), (c) and (d) (additional equilibria arising in cases BN and BS are denoted q w and q o as in Figure  4(c)).
It follows from previous work [19,23] that for each fixed µ = 0, the blow-up transformations (16) and (18) are sufficient to resolve all degeneracies in system (13). For µ = 0, an additional degeneracy persists due to the BE singularity. In this case, W becomes an attracting critical manifold W 0 , and connects to another degenerate point Q bf b ∈ {ρ = = 0} at the top of the blow-up sphere [19]. This case is shown in Figure  4(a). The point Q bf b is located at the origin in local coordinates (x 2 , ρ 2 , 2 ) defined by for µ = 0 only. Appending the trivial equation µ = 0 to the system obtained in these coordinates, Q bf b is identified as a nonhyperbolic equilibrium within the extended (x 2 , ρ 2 , 2 , µ)−space. Finally, degeneracy at Q bf b is resolved via the weighted spherical blow-up which replaces Q bf b with the 3-sphere {ν = 0} × S 3 . Following this spherical blowup, and a desingularization amounting to division by ν k(1+k) , the critical manifold W 0 terminates at a partially hyperbolic and (partially) attracting equilibrium q a contained within {ν =ρ =μ = 0}, see    (18). The critical manifold S in blue connects to the blow-up sphere at an attracting, partially hyperbolic point p a , and an attracting center manifold W, also in blue, emanates from p a over the blow-up sphere shown in orange. If µ = 0, all degeneracy is resolved. For µ = 0, the case shown here, W is a critical manifold W 0 which connects to the degenerate point Q bf b (magenta), which corresponds to the BE singularity. Local coordinates (x 2 , ρ 2 , 2 ) centered at Q bf b are also shown. (b) Dynamics and geometry following spherical blow-up of Q bf b via (20) in case BF. By restricting to the invariant set defined by the scaling (22), the blow-up 3sphere (magenta) can be projected into into (x,ρ,ˇ )−space as described in the text, and plotted in 3D. Following blow-up, W 0 connects to an attracting, partially hyperbolic point q a . An attracting center manifold J contained within {ν = 0}, also in blue, extends from q a onto the new blow-up sphere. Local coordinates (x 1 , ρ 1 , ν 1 ) defined via (24) used to describe the dynamics on the sphere are also shown. (c) resp. (d) Dynamics and geometry after blow-up in cases BN resp. BS. Here one identifies additional equilibria within {ν =ˇ =μ = 0}.
The sequence of blow-up transformations (16), (18) and (20) can be written in the following form upon composition: Remark 3.1. Note that the µ−coordinate is not shown in Figure 4 Due to the conservation of µ and the original small parameter it follows that is also a conserved quantity, even for ν = 0. This conserved quantity induces a foliation of the blow-up 3-sphere by lower-dimensional 2-spheres parameterized byμ ∈ R, thereby permitting a 3-dimensional representation as in Figure 4. In the following we will, when it is convenient to do so, viewμ as our bifurcation parameter on the sphere.
Lemma 3.2. A desingularized system governing the singular limit dynamics in the scaling regime defined by µ =μ k/(1+k) can be obtained from the doubly extended system (23) by an application of the coordinate transformation , followed by the desingularization and finally, restriction to the invariant subspace {ν 1 = 0} corresponding to = 0. The resulting system is where we write β := β + = φ + (0). Moreover, system (26) is topologically equivalent to Proof. The transformation (24) is simply obtained from (21) by working in the charť = 1 with chart-specific coordinates (x 1 , ν 1 , ρ 1 , µ 1 ) defined by , recall (22), which gives the desired result upon using this expression to eliminate µ 1 . From this, we obtain (27)  Both systems (26) and (27) are useful for describing the unfolding of singularly perturbed BEB in system (13). System (26) arises from a central projection of the final blow-up transformation, and is preferred for purposes of global computations within the blown-up space. System (27) is derived by a direct parameter rescaling, and although it is preferred for local computations pertaining to e.g. bifurcations, it is less suited to global analyses. 1 Notice however, that (27) can also be obtained more directly by composing (28) with (24). This gives after eliminating ν 1 . Inserting this into (13) gives (27) for → 0 upon desingularization.
It is also possible to scale x and y by µ for µ > 0 instead of ; in fact, this is more well-suited forμ → ∞. Therefore if we definê transforms (13) into following system for µ → 0 upon desingularization. System (31) is smoothly topologically equivalent to (27) onμ > 0 through the transformation The limitμ → −∞ can be studied via an analogous scaling by −µ for µ < 0, but we will not need this in our analysis.

Remark 3.3.
In this work we focus on the qualitative dynamics near a nondegenerate BE bifurcation. For general systems (1) with a BE bifurcation at (z, α) = (z bf , α bf ), however, Lemma 3.2 offers a direct route to obtain quantitative information about the dynamics without the need for bringing the system into normal form, by first shifting (z,α) = (z − z bf , α − α bf ), and then applying the coordinate transformation (24) and desingularization by (25) withz = (x, y).

3.2.
Unfolding all 12 singularly perturbed BEBs. The limiting bifurcation structure can be derived for each case using either of the systems (26) or (27). We may consider system (27) for simplicity, which by [19, Lemma 3.5] has either 0, 1 or 2 equilibria. The corresponding bifurcation diagrams with 1 are obtained after lifting results for = 0.
desingularized system (27), for which saddle-node, Andronov-Hopf and homoclinic bifurcations are identified along parameter space curves given by (35) respectively, whereμ hom (γ) is defined for γ < δ/τ and τ > 0 in a neighbourhood of the Bogdanov-Takens point (μ bt , γ bt ) : Remark 3.5. The corresponding statement for regularized PWS systems is obtained by taking the singular limit → 0 in Theorem 3.4. In this context results should be stated in terms of the desingularized system (26) (or (27)) alone, for which the identified bifurcations are described by equations (35), (36) and (37).
Theorem 3.4 yields four qualitatively distinct 2-parameter bifurcation diagrams. These are shown for the desingularized system (27), i.e. in the limit → 0, in Figure  5. Theorem 3.4 asserts that the corresponding diagrams for 1 sufficiently small are qualitatively similar. We make the following observations with respect to Figure  5: (i) All 12 singularly perturbed BEBs are represented: BF 1,2,3 in (a), BN 2,3 in (b), BN 1,4 in (c), and BS 1,2,3 in (d). Cases BF 4,5 are qualitatively similar to BN 1,4 respectively in (c). (ii) Cases for which the underlying PWS BE has an incoming (outgoing) Filippov flow, see again Figure 1 and Table 1, are contained within γ < 0 (γ > 0). (iii) Cases BF 1,2 are both contained within γ > 0 in (a). The homoclinic branch represents the continuation of the separatrix which constitutes a boundary between the two cases, with BF 1 (BF 2 ) lying the the left (right) of this curve. Theorem 3.4 only provides a local parameterisation of the homoclinic curve. A global parameterisation is not given in this work; homoclinic curves in Figure  5 have been obtained by numerical continuation using MatCont [8]. However, additional properties of the homoclinic branch in Figure 5(a) are also described in Section 3.3. (iv) Cases BS 1,2 are both contained within γ < 0, and separated by a distinguished solution which connects the unstable manifold of the saddle along the (unique) trajectory which is tangent to the strong eigendirection of the stable node. This is also discussed in Section 3. (vi) All bifurcations are 'singular' in system (13) in the sense that they occur within an −dependent domain which shrinks to zero as → 0, at a rate prescribed by the scaling (22). (vii) The BN 2 bifurcation in (b) features 'hidden oscillations', i.e. oscillations which cannot be identified in the PWS system (14), within the wedge-shaped region bounded by the Andronov-Hopf and homoclinic curves. (viii) The decay coefficient k ∈ N + associated with the regularization does not effect the topology of the bifurcation diagrams. It follows that within the class of regularizations defined by Assumptions 2-3, the observed dynamics are qualitatively independent of the choice of regularization. (ix) Each of (non-equivalent) 2-parameter bifurcation diagram in Figure 5 can be obtained from any of the others by suitably varying the additional parameters (τ, δ), either across one of the boundaries δ = 0, τ = 0 or ∆ = 0, or through the origin τ = δ = 0; see again Table 1. A complete description of the dynamics involves the unfolding a (singular) codimension-4 bifurcation at (τ, δ, γ,μ) = (0, 0, 0, 0). This unfolding is expected to involve (singular) codimension-3 bifurcations, and the unfolding of these bifurcations should involve the diagrams in Figure 5.
3.3. Separatrices: The boundaries between BF 1,2 and BS 1,2 . In this section we present a result on the homoclinic double-separatrix which constitutes a boundary between singularly perturbed BF 1,2 bifurcations. A heteroclinic double-separatrix forming a boundary between singularly perturbed BS 1,2 bifurcations is also discussed.
Remark 3.7. Notice that the outer regime covers µ ∼ whereas the inner regime covers µ ∼ k/(1+k) . The two regimes do not overlap for → 0. In principle, we should be able to cover the gap using our method, but we leave that for future work.
Combining Theorem 3.4(iv) and Proposition 3.6, we have asymptotic information about the branch of homoclinic solutions in Figure 5(a) forμ ∼μ bf ,μ 1 and µ ≥ 0. Our findings are represented schematically in Figure 6(a), which shows the expected global bifurcation diagram in (γ, µ, )−space after a weighted cylindrical blow-up  rescaling chart˜ = 1, the bifurcation diagram in Figure 5(a) is identified on the cylinder, i.e. within {η = 0}, which is invariant. The bifurcation diagram for µ > 0 is bounded above the cylinder in Figure 6(a).
Remark 3.8. The BS 1,2 boundary is formed by the distinguished heteroclinic connection which connects saddle and node equilibria, tangentially to the strong eigendirection of the node. In the PWS normal form (14) obtained in the in the dual limit → 0 + , µ → 0 − , this distinguished heteroclinic connection occurs for It is straightforward to obtain an analagous result to Proposition 3.6, describing inner and outer expansions of such a heteroclinic connection, see the illustration in Figure 6(b). For simplicity, we have decided against including this result. Furthermore, numerical investigations (see Figure 5(d)) support the existence of a simple (transverse) connection to the branch of saddle-node bifurcations with base along {(γ het (0), µ, 0) :μ <μ sn (γ het (0))} as shown in Figure 6(b). 3 . The case BN 3 in Figure 2 is somewhat special, due to the existence of a continuous family of PWS homoclinic cycles for µ = = 0. As indicated in Figure 7, we parametrize this family using the negative x−coordinate:

Explosion in case BN
for any s ∈ (0, s 0 ) with s 0 > 0 sufficiently small, where Γ X + (s) is the backward orbit of (−s, 0) following X + | µ=0 while Γ sl (s) is the forward orbit of (−s, 0) following X sl | µ=0 . Note that the orbits Γ(s) are homoclinic to a BN 3 singularity, and should not therefore be confused with homoclinic orbits Γ hom from Proposition 3.6 that are homoclinic to a hyperbolic sliding equilibrium on Σ.
Since Γ(s) only exists for parameter values µ = = 0 corresponding to a BE singularity, we are motivated to consider the problem within the blown-up space described in Section 3.1. Recall that on the sphere ν = 0 there exists an attracting two-dimensional center manifold J of a partially hyperbolic point q a . Essentially, this manifold provides an extension of the critical manifold onto the sphere ν = 0. At the same time, for the present case BN 3 , there is also a hyperbolic point q w on the sphere ν = 0, alongˇ =μ = 0 with a two-dimensional stable manifold S := W s (q w ), see Figure 4(c). Let Jμ and Sμ denote the manifolds obtained from J and S after restriction to the invariant subsets {μ = const.} defined via the scaling (22).
The following result identifies the existence of a heteroclinic connection between q a and q w which will play an important role in the unfolding of the PWS cycles. The situation is sketched in Figure 8. Lemma 3.9. For each fixed γ < 0, J and S intersect in a unique heteroclinic orbit connecting q a with q w .
A similar result was proven in [27,Proposition 2] in the context of the substratedepletion oscillator, which is degenerate as a BN 3 bifurcation (see Section 6.1 for further details). Nevertheless, the proof of Lemma (3.9), which will be given in Section 5, in the course of proving Theorem 3.10 below, will follow the proof of [27,Proposition 2]. Using the parameterμ defined in (22), the heteroclinic will be obtained for a unique valueμ =μ het (γ) corresponding to an intersection of manifolds Sμ and Jμ obtained as intersections of S and J with invariant level sets defined by (22). The existence of a heteroclinic connection produces a family of heteroclinic cycles {Γ(s)} s∈(0,s 0 ) with improved hyperbolicity properties, see Figure 9. In turn, this enables a perturbation of the PWS homoclinic (42) into limit cycles for 0 < 1.
Theorem 3.10. Consider system (13). Let and fix any ν ∈ (0, 1). Then for any c > 0 sufficiently small, there exists an 0 > 0 and an s 0 > 0 such that the following holds for each ∈ (0, 0 ): There exists a parameterized family of stable limit cycles A proof is given in Section 5. The limit cycles described in this theorem are O(1) with respect to . Although it is straightforward to use our method to connect these cycles with o(1) cycles (essentially taking c = K k/(1+k) in (43) with K > 0 sufficiently large, see also Remark 3.11) that are obtained as perturbations of the heteroclinic cycle on the sphere {ν = 0}, we have decided to focus on the O(1) cycles since (i) the result is easier to state, and (ii) we have not been able to connect the cycles all the way down to the Hopf bifurcation anyways (colloquial), recall Theorem 3.4. Such a connection requires global information of the limit cycles on the sphere, which we have not been able to obtain. Figure 8. Dynamics on the the blow-up sphere in casesμ =μ − , µ =μ het andμ =μ + , whereμ − <μ het <μ + . A 3-dimensional representation is possible after restricting to invariant subspaces defined by level sets (22). Part of the path followed by the equilibrium q n (μ) under µ−variation is shown in green. Jμ and Sμ denote manifolds obtained from J and S after restriction to {μ = const.} via (22). By Lemma 3.9, Sμ and Jμ intersect for a unique parameter valueμ =μ het , providing a heteroclinic connection from q a to q w , shown here in blue. This connection breaks regularly asμ is varied overμ het . Dynamics on each side of the connection are also shown, in dark blue and red.
Remark 3.11. Figure 9 also indicates the existence of a family of nondegenerate singular cyclesΓ i , bounded between 'small' and 'large' heteroclinic cyclesΓ s andΓ l respectively. Such a construction is straightforward, and similar to the construction of singular cycles given in [19,Section D.1]. By analogy to the arguments presented in [19] in the case of singularly perturbed BF 3 bifurcation, it is possible to prove a connection between limit cycles that are O(1) with respect to and limit cycles on the (second) blow-up sphere, facilitated by a family of limit cycles obtained as perturbations of the singular cyclesΓ i .

Proof of the theorem 3.4 and proposition 3.6
In this section we prove Theorem 3.4 and Proposition 3.6. We begin with a proof of Theorem 3.4.

4.1.
Proof of the Theorem 3.4. We proceed by studying the dynamics of the relevant desingularized system from Lemma 3.2. Theorems 3.4 will follow immediately after lifting system (26) out of the invariant plane {ν 1 = 0} into {ν 1 ∈ [0, σ)} for sufficiently small σ > 0 and applying the blow-down transformation given by the inverse to (24) defined on { > 0} = {ν 1 > 0, ρ 1 > 0}. Figure 9. Nondegenerate singular cycles obtained whenμ =μ het by concatenating orbit segments, after the resolution of all degeneracy via the sequence blow-up transformations described in Section 3.1. The cy-clesΓ(s 1 ) andΓ(s 2 ) shown in red correspond to the PWS cycles Γ(s 1 ) and Γ(s 2 ) in Figure 7, respectively. In terms of the dynamics after blowup, Theorem 3.10 describes the existence and growth of limit cycles obtained as perturbations of singular cyclesΓ(s) with s > 0, i.e. with orbit segments bounded away from the blow-ups spheres. Perturbations of singular cyclesΓ s ,Γ l , and the family of cycles bounded between (represented here byΓ i ), are not described by Theorem 3.10, see Remark 3.11. It is possible to show as in [19, Theorems 3.11 and E.1], however, that these cycles mediate a connection to limit cycles on the (magenta) blow-up sphere. Transversal sections Σ 1 and Σ 2 used in the proof of Theorem 3.10 are also shown.
Lifting the expressions derived above for ν 1 ∈ [0, σ) with σ > 0 sufficiently small and applying the blow-down transformation, in particular the relation we obtain the desired result.
Remark 4.1. In the preceding proof σ must be sufficiently small so that (35), (36), (37) and (46) can be extended in (x 1 , ρ 1 , ν 1 , µ 1 )−space via suitable applications of the implicit function theorem. We omit this argument -which is standard -for the sake of brevity, but refer the reader to [19, eqn. (D7)] where the extended system is considered in detail.
We therefore focus on the inner expansion in the dual limit case, setting = µ (1+k)/kˆ and letting µ → 0. For this we consider system (31). The case of k = 1 is easier, so we will also focus on this case, repeated here for convenience where we have dropped the hat notation on X and Y . We leave the discussion of the general case k ∈ N to the end of the section.
The system (48) is for 0 <ˆ 1 a slow-fast system in nonstandard form [20,38]. Indeed, forˆ = 0 we obtain the layer problem for which {Y = 0} is a manifold of equilibria. Linearization of any point (X, 0) gives X as the only nonzero eigenvalue. Consequently, S a := {(X, 0) : X < 0} is normally hyperbolic and attracting, (0, 0) is fully nonhyperbolic, and S r := {(X, 0) : X > 0} is normally hyperbolic and repelling. Notice also that for Y > 0 we obtain the equivalent system upon dividing the right hand side of (49) by Y . Let φ t denote the flow of (50). We then define Γ as {φ t (0, 0)} t∈(0,t d ] where t d > 0 is the first return time to Y = 0. Notice that Γ is well-defined since (50) is just the linearization of the vector field X + having, in the BF case considered, an unstable focus at (0, δ −1 ). It is a simple calculation to show that t d > 0 is the first positive root of R(t), recall (39), and that Γ ∩ {Y = 0} = (X d , 0) with Next, setting Y =ˆ Y 2 brings (48) into a slow-fast system in standard form. Upon passing to a slow time and then settingˆ = 0, we obtain the following reduced problem on S a :Ẋ  Figure 10. Singular limit dynamics for the nonstandard form slowfast system (48) arising in case k = 1. Attracting and repelling critical manifolds S a and S r are shown in blue and red respectively. The point (0, 0), shown in orange, is a regular fold point. There is an unstable focus at (0, δ −1 ), and an equilibrium (−γ −1 , 0) which is repelling as an equilibrium for the reduced flow on S a ; both are indicated as black disks. We show the situation where γ = γ hom,0 = −X −1 d with X d given by (51), for which there is a singular homoclinic orbit Γ (shown here in black).
which has a repelling equilibrium at X = −γ −1 , seeing that γ > 0. We note that reduced problem can also be obtained from more general procedures described in [20,38].
Combining our analysis of the layer problem and the reduced problem, we obtain Figure 10. Specifically, for γ = γ hom,0 := −X −1 d we have a singular saddle-homoclinic connection. At the singular levelˆ = 0, the connection is clearly transverse with respect to γ; in fact, Γ is independent of γ so this is obvious from X d = X d (γ). For 0 <ˆ 1 we then use Fenichel theory to perturb the saddle and the result of [28] to track its unstable manifold near Γ. Defining a section Σ transverse to Γ within Y > 0, we then obtain a bifurcation equation for the saddle-homoclinic connection of the form H(γ,ˆ ) = 0, with H, which measures the separation of the stable and unstable manifolds on Σ, being at least C 1 in γ, continuously dependent onˆ ∈ [0,ˆ 0 ) and such that The existence of γ inner hom in Proposition 3.6 follows after applying the implicit function theorem to H(γ,ˆ ) = 0 at (γ,ˆ ) = (γ hom,0 , 0). For the final part of Proposition 3.6, we fixˆ small enough (i.e.μ large enough) and perturb in µ > 0 (or equivalently > 0, having fixedˆ ) small enough.
Remark 4.2. Notice for this last part that the µ-perturbation of (48) will include terms of the form using (30), which are ill-defined for µ = Y = 0. This is in the sense of which the charts (30) are more ill-suited for global computations. Here, however, fixingˆ > 0, where the saddle connection occurs within Y > 0, we just require that the perturbation is continuous with respect to µ on this domain. To cover the saddle-homoclinic case in a full neighbourhood of ( , µ) = (0, 0), we have to work with our full blowup system, tracking the saddle across the first blowup sphere. We leave the details of this to future work.
For k ≥ 2, {Y = 0} is fully nonhyperbolic forˆ = 0. We then gain hyperbolicity by blowing up the points (X, 0, 0) in the extended (X, Y,ˆ )−space via followed by a desingularization corresponding to division of the right hand side by r k−1 . Working in the directional chart corresponding toȲ = 1 using chart-specified coordinates (r 1 , Y 1 , E 1 ) defined by Y = r 1 ,ˆ = r 1 E 1 , we find a normally hyperbolic and attracting critical manifold on r 1 = 0, carrying a reduced problem given by (52).
We therefore obtain the result as in the k = 1 case, performing a separate blowup of (X, r 1 , E 1 ) = (0, 0, 0), replacing the result of [29], to track the slow manifold for Y > 0 in this case. We leave out the details for simplicity.
The desingularized equations in the chartρ = 1 take the following form: (58) The quantityμ is conserved for the flow of (58).
Proof. This follows by lengthy, but standard calculations. We defer the proof to Appendix B for expository reasons.
On the other hand, within ρ 1 = µ 1 = 0 we have (62) Here we find the fully hyperbolic equilibrium q f with x 1 = ν 1 = 0. In particular, a simple calculations shows that within ν 1 = 0, q f is a source. On the other hand, for (62) we also find x 1 = −β, ν 1 ≥ 0 as the critical manifold W 0 , see Figure 4, of partially hyperbolic points. Indeed, the linearization of any point on W 0 has a single nonzero eigenvalue −kβ, also at the point q a with coordinates (x 1 , ρ 1 , ν 1 , µ 1 ) = (−β, 0, 0, 0) ∈ W 0 . At q a , we therefore have a three-dimensional attracting center manifold. We shall denote the ν 1 = 0 subset of this manifold by J , as indicated in Figure 4. Using the parameterμ, we may foliate J into invariant subsets Jμ. For simplicity, we denote the projection of Jμ onto the (x 1 , ρ 1 )-subspace by the same symbol. Then Jμ becomes an attracting center manifold of the point (x 1 , ρ 1 ) = (−β, 0), which we for simplicity also denote by q a , for the system (26). A simple calculation shows that it takes the following smooth graph form: over ρ 1 ≥ 0 locally near q a . This gives and ρ 1 > 0 is therefore locally increasing on Jμ. In conclusion, we have the following.
Lemma 5.2. Consider (26). Then q a is a nonhyperbolic saddle on {ρ 1 ≥ 0} and the center manifold Jμ is unique on this set as the nonhyperbolic unstable manifold of q a for allμ ∈ R.
We will need the following result in our proof of Lemma 3.9. Lemma 5.3. There exists aμ − such that the ω-limit set of Jμ is q n (μ) for allμ ≤μ − .
Proof. The result follows from the center manifold theory, the fact that q n (μ) is a stable node forμ −1 and finally that q n (μ) → q a forμ → −∞.
As a corollary, the α-limit set of S ∞ is q n,∞ . But then by regular perturbation theory, and the hyperbolicity of q n,∞ , we obtain the following result, which we also need in our proof of Lemma 3.9.
Forμ ∈ R, we project Sμ onto the (x 2 , 2 )−space and denote the projection by the same symbol. A simple calculation shows that it takes the following local form: for 2 > 0 small enough. 5.3. Proof of Lemma 3.9. To prove Lemma 3.9, we combine our analyses in chartš = 1 andρ = 1 in order to show the existence of a uniqueμ het such that Jμ het intersects Sμ het , transversally with respect toμ.
• The α-limit set of Sμ is q f . Consider anyμ ≥μ + . Then: • The ω-limit set of Jμ is q o .
Following this result, we therefore fix an interval I = [μ − ,μ + ] ofμ-values, and then insert a section Σ at ρ 1 = c for c > 0 small enough. The ρ 1 -nullcline intersects Σ in a tangency point (x 1,t , c) for x 1,t := −β so thatρ 1 ≷ 0 for all points on Σ with x 1 ≷ −β. By the previous analysis any heteroclinic connection intersects Σ with x 1 > x 1,t . The center manifold calculation, see (63), shows that the manifold Jμ intersects the section Σ transversally in a point (x 1,c (μ), c) for eachμ ∈ I with x 1,c (μ) > x 1,t , for allμ ∈ I so thatẋ 1 > 0,ρ 1 > 0 in a neighbourhood of J ∩ Σ. By Lemma 5.6, we have that forμ =μ − the manifold Sμ intersects Σ in a point (x 1,s (μ − ), c) with x 1,s (μ − ) > x 1,c (μ − ). The intersection is therefore transverse and we can continue x 1,s (μ) smoothly for larger values ofμ >μ − . However, by Lemma 5.6 we know that Sμ does not intersect Σ for allμ ∈ I. The process of continuing x 1,s for larger values ofμ >μ − will therefore have to stop when either: x 1,s grows unboundedly or x 1,s → x + 1,t . We can exclude the former by the analysis in theρ = 1 chart. Therefore there is aμ t >μ − such that x 1,s (μ) → x + 1,t forμ →μ − t . With x 1,c (μ t ) > x 1,t we conclude that the smooth function:μ → x 1,c (μ) − x 1,s (μ) forμ <μ t changes sign at least once. The corresponding root corresponds to a heteroclinic connection. This connection is unique by the monotonicity of ρ 1,het (t) and the fact that the associated Melnikov integral, being the derivative of the Melnikov distance function, has one sign. To see the latter, consider (26) and notice that the derivative of the right hand side with respect toμ is (ρ , 0). Therefore the sign of the Melnikov integrand [31] is determined by This also shows that the intersection of J and S is transverse, completing the proof of Lemma 3.9.

5.4.
Finishing the proof of Theorem 3.10. In Figure 9 we combine our findings into a new figure illustrating the improved singular cycles Γ(s), see the figure caption for further details. We then obtain the family of attracting limit cycles in Theorem 3.10 with the prescribed growth rate by perturbing Γ(s). For this, we work nearμ ≈μ het and define two sections Σ 1 and Σ 2 as illustrated in Figure 9. We then flow points on Σ 1 forward and backward and measure their separation on Σ 2 . The sections are defined in the chartρ = 1 with coordinates (x 2 , 2 , ν 2 , µ 2 ), recall (54), as follows: for χ, ν > 0 small enough and I a small enough neighbourhood of x 2,w such that the following local arguments apply near q w = (x 2,w , 0, 0, 0), recall (65). The bifurcation equation is then given by where F and B are defined as the x 2 -coordinates of the points on Σ 2 obtained by following initial conditions (x 2 , 2 , ξ, µ 2 ) on Σ 1 forward and backward, respectively, where µ 2 = k/(1+k) 2μ . Notice, by conservation of andμ, solutions of (70) with 2 > 0 define closed orbits. Let J be a sufficiently small neighbourhood ofμ =μ het . Then we have the following.

5.5.
Proof of Proposition 5.7. The properties of F are standard; see [19]. Here we just summarize the approach. Firstly, as we flow points Σ 1 forward, following Γ X + (s), they are eventually exponentially contracted towards the slow manifold. We then guide the flow along this manifold, first along W 0 , and eventually along Jμ using the center manifold at q a . Seeing thatμ ≈μ het , we conclude that solutions reach Σ 2 so that F is well-defined. In particular, F (x 2 ,μ, 0) is the x 2 -value of the intersection Jμ ∩ Σ 2 .
For B we proceed more carefully. Firstly, we recall that q w is fully hyperbolic. Indeed, the linearization has the following non-zero eigenvalues: where λ = 2 √ ∆/(τ − √ ∆) as defined in Theorem 3.10. To describe the map Σ 1 → Σ 2 defined by the backward flow of (58), we use partial linearizations (since resonances in general will preclude a full linearization) within the invariant spaces { 2 = 0} and {ν 2 = µ 2 = 0} in order to obtain the following.
Lemma 5.8. There exists a C 1 diffeomorphism bringing (58) into the following system (72) upon a regular reparametrization of time.
Proof. See Appendix C.
We now integrate (72) backwards from ν 2 = ξ to 2 = χ. A simple calculation, based upon a Gronwall-type estimate, gives which defines B in the new local coordinates. From this, we then similarly obtain the desired estimate of the x 2 -derivative of B in Proposition 5.7 using the relevant variational equations (these are taken with respect to system (83) in Appendix C). Moreover, it is clear that B(x 2 ,μ, 0) coincides with the x 2 -coordinate of the intersection Sμ ∩ Σ 2 . Consequently, (71) holds by Lemma 3.9 and the transverse intersection of J and S.

Outlook
The unfolding of BE singularities in BEB under parameter variation is generic in PWS systems. It follows that singularly perturbed BEB is also generic under parameter variation in singular perturbation problems losing smoothness along a codimension-1 switching manifold Σ as a perturbation parameter → 0. In this manuscript the notion of singularly perturbed BEB is formally defined (see Definition 2.5), and a classification based on known classifications for PWS BEB from [10,14,32] is given; see Table 1. We showed in Theorem 2.8 that the local normal form (13) first derived in [19] for the analysis of singularly perturbed BF bifurcations in particular, is capable of generating all 12 singularly perturbed BEBs. It is worthy to note that a corresponding PWS normal form (14) is also obtained in the limit → 0.
Following the introduction of the normal form (13), we studied its dynamics in parameter regions corresponding to each singularly perturbed BEB. Using a sequence of blow-up transformations to resolve a loss of smoothness along Σ, and subsequently, degeneracy arising from the BE singularity itself, we derived two desingularized systems (26) and (27) in Lemma 3.2. Studying the dynamics of these systems allowed for a detailed description of the unfolding for all 12 singularly perturbed BEBs. This was presented succinctly in Theorem 3.4 and Figure 5. In many cases, we were able to provide explicit paramterizations for the location of codimension 1 and 2 bifurcations involved in the unfoldings. It is worthy to note that in general, the bifurcation structure depends quantitatively, but not qualitatively, on the decay rate k determining the rate at which the system loses smoothness along Σ as → 0; see equation (6). In particular, all identified bifurcations are singular, in the sense that they occur within a parameter regime µ = O( k/(k+1) ) which shrinks to zero in the PWS limit → 0.
We then demonstrated the suitability of our framework for studying the geometry of so-called double-separatrices, which constitute non-trivial boundaries between cases BF 1,2 and BS 1,2 in parameter space. A result on the boundary between cases BF 1,2 was presented in Proposition 3.6, but a complete analysis of this and the BS 1,2 boundary is left for future work.
Finally, special attention was devoted to the so-called BN 3 -explosion, which may be considered a 'generic analogue' of the (degenerate) explosion identified already in [27], which can be considered as the 'γ = 0 case' of system (13). We showed that the continuous family of singular cycles shown in Figure 9 perturbs to a continuous family of stable limit cycles for 0 < 1. This is described in Theorem 3.10, where the growth rate of the cycles is also quantified as a function of and the parameter k determining the rate at which the system loses smoothness. We emphasize that the focus of [27] was on the existence of relaxation oscillations, whereas here we focus on the details of the explosion itself.
We conclude with some discussion on the relation to explosive onset of oscillations in classical slow-fast systems and other singularly perturbed BEBs is discussed below. Applications, more degenerate cases of interest, and singular bifurcations of higher codimension are also considered.
6.1. Relation to classical canard explosion and singularly perturbed BF 3 explosion. In Theorem 3.10 we described a novel explosion mechanism of limit cycles due to the BN 3 bifurcation. This 'explosion' is reminiscent of the canard explosion phenomenon in classical slow-fast systems. Here too, an entire family of singular cycles exists for a unique parameter value. Geometrically both families are upon blow-up identified in a similar way through heteroclinic cycles; for Theorem 3.10 the heteroclinic cycles occur due to the transverse intersection of the manifolds Jμ and Sμ. Moreover, in both cases, the singular cycles perturb to limit cycles for 0 < 1, see e.g. [9,29,30]. However, the 'explosion' described by Theorem 3.10 differs from the classical canard explosion phenomenon in a number of important respects. First, for the BN 3 bifurcation, there is only an attracting slow manifold. Since there is no repelling slow manifold, there are no canards. Instead, repulsion in the BN 3 explosion comes from the unstable node. Consequently, the limit cycles in Theorem 3.10 are also always stable since the contraction of the slow manifold dominates the hyperbolic repulsion from the node; in the canard case the limit cycles can be either attracting, repelling or neither, the details depending on a slow divergence integral [29]. These differences also manifest themselves through different growth rates. In the BN 3 case, the growth rate (44) is algebraic, whereas the classical canard explosions are characterised by exponential growth (since (44) is exponentially small in this case).
The onset of oscillations that are O(1) with respect to in the singularly perturbed BF 3 bifurcation, described in detail in [19], is again different. In the BF 3 case one also identifies a family of singular cycles, but these cycles all lie on the blow-up spheres. As a consequence, the family of limit cycles obtained after perturbation and blow-down is not explosive in any way. We will discuss this further in the following section in the context of a regularized stick-slip oscillator capable of producing both bifurcations BN 3 and BF 3 (albeit degenerate ones) upon parameter variations.
6.2. Degenerate singularly perturbed BE bifurcation in applications. Singularly perturbed BN 3 explosion of the kind described in Theorem 3.10 occurs in a regularized Gause problem; see [19] for the regularized model, and [11] for the original PWS system. As previously stated, a singularly perturbed BN bifurcation also occurs in the model for substrate-depletion in [27], where it shown to provide a mechanism for the onset of the relaxation-type oscillations. In the context of the normal form (13), this case is 'degenerate' due to γ = 0. As shown in [27], this produces a critical manifold S with no reduced flow, however with an additional (infra-)slow timescale.
Similar (degenerate) explosions may also be observed in regularized stick-slip oscillators under variation of the belt speed, see e.g. [19,Section 5.1]. This model takes the form (1) satisfying Assumption 1 with where µ describes the friction law. In [23] it is given as which was proposed by [2] and also studied in [35,39]. Here µ(0) = µ s and µ m > µ m > 0, ρ > 0, c ∈ (0, ρ(µ s − µ m )) to ensure that µ (0) < 0. The model with (73) has a BEB for α = 0 at (x, y) = (−µ s , 0). An easy calculation, see also [19], shows that this bifurcation can be either a degenerate BN 3 for µ (0) < −2, or a degenerate BF 3 for µ (0) ∈ (−2, 0). Although these bifurcations are degenerate with γ = 0, we nevertheless use this example to illustrate in Figures 11 and 12 the differences between these cases. Specifically, in Figure 11 (76), corresponding to (degenerate) BF 3 and BN 3 bifurcations, respectively. The thinner lines track the equilibrium whereas the thicker curves correspond to the limit cycles (denoted LC in the figure), using max x as a measure of the amplitude, emerging from the two Andronov-Hopf AH bifurcation points indicated by two black circles. The limit cycles corresponding to the purple, blue and red points (disks for the BF case and squares for the BN case) are illustrated in Figure 12(a) and (b), respectively. Homoclinic and heteroclinic branches appear to extend and intersect for γ = 0. The degenerate singularly perturbed BN bifurcation for γ = 0, i.e. between cases BN 2 and BN 3 , has been described in detail already in an application in [27]. Corresponding dynamics after blow-up are shown for γ < 0, γ = 0 and γ > 0.
These two cases give µ (0) = −1.10 and µ (0) = −2.90, respectively, and therefore correspond to (degenerate) BF 3 and BN 3 cases. We use a regularization (77) φ(y) = 1 2 1 + y y 2 + 1 which satisfies Assumptions 1, 2 and 3 with k = 2, and set = 0.001 in system (1). As discussed, we only see an explosive growth of the limit cycle amplitude in the BN 3 case. For further comparison, Figure 12 illustrates examples of limit cycles; (a) in case BF 3 , and (b) in case BN 3 . The colours correspond to the colours of the points in Figure 11. See figure captions for further details.
6.3. Connecting BN 2,3 across γ = 0. Since Theorem 2.8 applies for all γ ∈ R, it follows that both smooth and PWS normal forms for degenerate BE bifurcations with γ = 0 are obtained by setting γ = 0 in (13) and (14) respectively. Hence, it is expected that much of the analysis presented herein can be applied in order to study these degenerate cases, which can be thought of as 'boundary cases' separating (singularly perturbed) BF 1,3 , BF 4,5 , BN 2,3 , BN 1,4 and BS 1,3 bifurcations; see again Figure 5 and the caption. The BN 2,3 boundary is particularly interesting as (i) it arises naturally in the context of substrate-depletion oscillations as described above [27], and (ii) there is evidence that upon extension through γ = 0, the homoclinic branch identified in case BN 2 with γ > 0 connects to the heteroclinic branch for γ < 0 which is responsible for the explosion in case BN 3 . Computing the heteroclinic connection of Lemma 3.9 numerically and plotting it over the bifurcation diagram in Figure 5(b), we obtain the diagram in Figure 13. Here we see the expected transition from case BN 3 to BN 2 as γ crosses zero. Our observations also provide evidence thatμ het (γ) <μ ah (γ) on γ < 0. Note that ifμ het (γ) is an analytic continuation of µ hom (γ) through γ = 0, then despite appearances in Figure 13, it cannot be linear. This would follow from the nonlinearity in the local parameterization of µ hom (γ) near the BT point.
6.4. Higher codimensions. As highlighted in observation (ix) following the statement of Theorem 3.4, the diagrams in Figure 5 can be obtained from one another under suitable variation in two additional parameters τ and δ. In particular, it follows from the classification in Table 1 that δ = 0, τ = 0 and ∆ = 0 form boundaries between the cases represented in Figure 5. For example, panel (b) is obtained from panel (a) by crossing from ∆ < 0 to ∆ > 0 with τ, δ > 0. To see this transition in the blown-up space, compare Figures 4(b) and 4(c). For ∆ = 0, a saddle-node bifurcation occurs on the intersection of the upper blow-up sphere (shown in magenta) with the plane { = 0}, giving rise to the equilibria q w and q o . Our analysis provides a framework within which transitions such as these can be analysed, thereby providing a solid foundation and program for future work.