Blowup Analysis of a Hysteresis Model Based Upon Singular Perturbations

In this paper, we provide a geometric analysis of a new hysteresis model that is based upon singular perturbations. Here hysteresis refers to a type of regularization of piecewise smooth differential equations where the past of a trajectory, in a small neighborhood of the discontinuity set, determines the vector-field at present. In fact, in the limit where the neighborhood of the discontinuity vanishes, hysteresis converges in an appropriate sense to Filippov’s sliding vector-field. Recently (2022), however, Bonet and Seara showed that hysteresis, in contrast to regularization through smoothing, leads to chaos in the regularization of grazing bifurcations, even in two dimensions. The hysteresis model we analyze in the present paper—which was developed by Bonet et al in a paper from 2017 as an attempt to unify different regularizations of piecewise smooth systems—involves two singular perturbation parameters and includes a combination of slow–fast and nonsmooth effects. The description of this model is therefore—from the perspective of singular perturbation theory—challenging, even in two dimensions. Using blowup as our main technical tool, we prove existence of an invariant cylinder carrying fast dynamics in the azimuthal direction and a slow drift in the axial direction. We find that the slow drift is given by Filippov’s sliding vector-field to leading order. Moreover, in the case of grazing, we identify two important parameter regimes that relate the model to smoothing (through a saddle-node bifurcation of limit cycles) and hysteresis (through chaotic dynamics, due to a folded saddle and a novel return mechanism).

In this paper, we consider piecewise smooth (PWS) systems of the following form: where z = (x, y) ∈ R n+1 , Z ± (z) = (X ± (z), Y ± (z)). The set Σ : y = 0 is called the discontinuity set or switching manifold. In a more general setting, one could define the switching manifold Σ as a smooth hypersurface h(z) = 0 for some regular function h : R n+1 → R. Locally, however, we can always introduce coordinates (x, y) so that h(x, y) = y. We will suppose that Z ± are smooth vector-fields, each defined in a neighborhood of Σ.
The basic problem of (1) is how to define solutions of (1) near Σ. The case when Y + (x, 0) < 0 and Y + (x, 0) > 0 is most interesting from a technical point of view, because in this case orbits of either systemż = Z ± (z) reach Σ in finite time, see Fig. 1. This is known as (stable) sliding. To be able to define a forward flow, a vectorfield must be assigned on Σ. The most common way to do this, is through the Filippov vector-field defined by The PWS systems, where (2) is assigned along the subset of the switching manifold with Y + (x, 0)Y − (x, 0) < 0, are called Filippov systems. These type of systems occur naturally in mechanics, e.g. in friction modelling [3,14]. However, even such mechanical models may suffer from nonuniqueness of solutions, and a meaningful forward flow may not be defined at all points [3].
From a modelling perspective, nonuniqueness may be interpreted as an insufficient model where additional information or complexity has to be added in order to select a unique forward trajectory. From this point of view, it is therefore important to study regularizations of (1). There are two basic regularizations of (1), one is by smoothing and another one is by hysteresis. In this paper, we shall -following Sotomayor and Teixeira [22] -define regularization by smoothing as replacing (1) with an -family of Figure 1. Illustration of Filippov's sliding vector-field X sl in the case of stable sliding. smooth systems: with Z(z, p) = (X(z, p), Y (z, p)) defined by: for Z ± = (X ± , Y ± ). Regarding the function φ in (3) , we assume the following assumption, so that (3) approaches (1) pointwise for → 0 for y = 0. Assumption 1. The smooth 'regularization function' φ : R → R satisfies the monotonicity condition ∂φ(s) ∂s > 0, for all s ∈ R and, moreover, On the other hand, in hysteresis, solutions ofż = Z + (z) are extended to y = −α before switching toż = Z − (z). Here α > 0 is some small parameter. Solutions oḟ z = Z − (z) are similarly extended to y = α before switching occurs, see Fig. 2.
In smoothing, we basically introduce a boundary layer of order O( ) around y = 0 where p = φ(y −1 ) changes by an O(1)-amount. From this point of view, it is also useful to think of hysteresis as introducing a "negative" boundary layer around y = 0 of size 2α.
In both types of regularizations, forward solutions can be uniquely defined for all > 0, α > 0, respectively, and in some cases, , α → 0 can be analyzed. For example, in [15] the authors studied the regularization by smoothing of a visible-invisible twofold in R 3 . The two-fold is a well-known singularity of Filippov systems that give rise to nonuniqueness of solutions. The results of [15] showed that the smooth system has a well-defined limit for → 0 which selects a distinguished forward trajectory through the two-fold. In this way, the nonuniqueness has (in a certain sense) been resolved.
The references [13] also studied regularization by smoothing but considered the planar grazing bifurcation scenario, where a limit cycle of Z + grazes Σ while Z − remains transverse and points towards Σ. The results showed, in line with [1] and analysis based upon the associated Filippov system [18], that in the case of a repelling limit Figure 2. Illustration of hysteresis in the case of stable sliding. In comparison with regularization by smoothing, hysteresis can be interpreted as introducing a negative boundary layer, the size of which is given by 2α. Figure 3. Illustration of the planar visible fold, which is important for the grazing bifurcation. The visible fold separates the switching manifold Σ into stable sliding (x < 0) and crossing points (x > 0). cycle, the smooth system has a locally unique saddle-node bifurcation of limit cycles. The analysis rested upon a careful description of the local dynamics, through a local transition map P loc : Π in → Π out , near the grazing point which is given by a visible fold, see illustration in Fig. 3. These results, together with [9,10,11,16] working on similar systems, were obtained by adapting methods from Geometric Singular Perturbation Theory [12,6]. In particular, these references use a modification of the blowup method [5,17] to gain smoothness of systems of the form (3).
Recently, in [21] the authors performed a related study but using regularization by hysteresis instead. Interestingly, the results are completely different in this case. In fact, hysteresis leads to chaotic dynamics for any 0 < α 1 under the same assumptions.
In this paper, we consider a new regularization of (1) [2]: (6)ẋ = X(z, p), y = Y (z, p), for 0 < , |α| 1. Notice that the dimension of (6) is one greater than the dimension of (1). The connection between (6) and (1) at the pointwise level is as follows: By assumption 1, (6) converges pointwise to (7)ż = Z(z, p), p = 1 y > 0 0 y < 0 , for , α → 0 for y = 0, which upon using (4) projects to (1). This model was introduced by [2], in a slightly more general framework where Z(z, p) depends nonlinearly p, with the purpose to incorporate smoothing and hysteresis in one single unified framework. The authors present asymptotic results for both α < 0 and α > 0 connecting the dynamics of (6) with Filippov's sliding vector-field in certain parameter regimes. Since trajectories in hysteresis cross each-other in the "negative boundary layer", recall Fig. 2, it makes sense that the smooth model (6) is defined in an extended space.
In [2], the authors consider functions φ(s), see assumption 1, that reach 0 and 1 at finite value. Such functions have -following [22] -been called Sotomayor-regularization functions. In this paper, also to exemplify the power of our approach, we will consider general regularization functions that are truly asymptotic, like analytic ones, e.g. φ(s) = 1 2 + 2 π arctan(s).
For this purpose, we add the following technical assumption:  There could be different k-values k ± for s = ±∞, respectively, but for simplicity we take these to be identical. In the following, β − will play little role so we will therefore for simplicity write β + as β.
Assumption 2 ensures the folllowing (by Fenichel's theory [6]): Lemma 1.1. Consider ( 6) on a compact subset U + upon which y > 0. Then by assumption 2, we haveẋ = X(z, p), y = Y (z, p), This system has an attracting slow manifold S ,α taking the graph form p = 1 + O( k |α| k ) which carries the reduced problem: This holds uniformly and smoothly on the compact subset U + and on this set ( 10) is therefore a smooth O( k |α| k )-perturbation ofż = Z + (z).
A similar result clearly holds within y < 0. The objective of our analysis is to uncover what occurs near y = 0.
It is clear that (6) involves a combination of slow-fast and nonsmooth effects. The analysis of such system seem to be rare. The reference [14] offers an exception. This manuscript studied a model of a friction oscillator, also of the form (3) but with φ nonmonotone, in the presence of a time scale separation. The combination of slow-fast and nonsmooth effects was shown to lead to chaotic dynamics through a horseshoe obtained through a folded saddle singularity [23] and a novel return mechanism. We will obtain something similar for (6) in the case of the grazing bifurcation. However, the analysis of (6) is more involved, and as opposed to the system in [23], the slow-fast and nonsmooth effects are more combined. On top of that, (6) involves two small parameters 0 < , |α| 1. In this paper, we will focus on α > 0; although the case α < 0 is somewhat different, it can be handled by the same methods and we therefore leave out its details of the present manuscript. In comparison with [2], we present a more detailed geometric picture of the dynamics and do not rely on any particular order dependency between the small parameters > 0 and α > 0 in our description of sliding. Moreover, we will also consider the grazing bifurcation in the context of (6) and obtain a (partial) unification of the results by smoothing and hysteresis.
1.1. Overview. The remainder of the paper is organized as follows. In Section 2, we present the blowup approach which will form the basis for our analysis of (6) with α > 0. This will include a review of this method in the context of regularization by smoothing. In Section 3, we then present the first main results, summarized in Theorem 3.1, on the dynamics of (6) for , α > 0 both sufficiently small in the case of stable sliding. In proving this result, we also lay out the geometry of the dynamics using our blowup approach. In Section 4, we prove an important lemma used to prove Theorem 3.1. In Section 5, we then turn our attention to the grazing bifurcation. The main results of this section are stated in Theorem 5.1. Here we focus on two parameter regimes of the ( , α)-plane with distinct dynamics. In one regime, we obtain a locally unique saddle-node bifurcation, as in regularization by smoothing, while in the other regime, we obtain chaotic dynamics through a horseshoe. The proof of these results rely on separate blowup transformations to resolve further degeneracies due to the grazing bifurcation.

Blowup
The blowup approach [5,17], which in its original framework was developed as a method to deal with lack of hyperbolicity, has recently been adapted [15] to deal with smooth systems approaching nonsmooth ones. Within this framework, we gain smoothness rather than hyperbolicity by applying blowup.

2.1.
A blowup approach for regularization by smoothing. A particular emphasis in the development of blowup for smooth systems approaching nonsmooth ones, has been on systems of the form (3). Within our context these systems correspond to regularization of the PWS system (1) by smoothing and the blowup approach proceeds as follows: Firstly, we work in the extended space (x, y, ) by adding the trivial equation˙ = 0. At the same time, to ensure that = 0 is well-defined, we consider this extended system in terms of a fast time: (11) x = X z, φ y , Then (x, y, 0) is a set of equilibria, but (x, 0, 0) is extra singular due to the lack of smoothness there. This set is therefore blown up through a cylindrical blowup transformation defined by for r ≥ 0, (ȳ,¯ ) ∈ S 1 , leaving x fixed. Here S 1 is the unit circle in R 2 . Notice that r = 0, (ȳ,¯ ) ∈ S 1 maps to (y, ) = (0, 0) and by the inverse process the points (x, 0, 0) are therefore blown up to a cylinder. See Fig. 4.
Under the assumption 2, it follows that (12) gives a smooth vector-field V on (x, r, (ȳ,¯ )) ∈ R n × [0, ∞) × S 1 by pullback of the vector-field associated with (11). In fact, due to the multiplication by , V has¯ as a common factor. However, V :=¯ −1 V is well-defined, smooth on R n × [0, ∞) × S 1 , and nontrivial for¯ = 0, see [15]. Whenever either Z ± is transverse to Σ then V has good hyperbolic properties. This holds true even if Z(z, p) depends nonlinearly on p, but we shall henceforth suppose the following: Assumption 3. Z(z, p) is affine with respect to p as in ( 4).
Suppose also that: Assumption 4. The PWS system ( 1) has stable sliding along the discontinuity set Σ: We also assume that Σ is a compact domain in R n .
Then we have the following.
Proposition 2.1. [15,19,22] Consider ( 3) and suppose that assumptions 3 and 4 both hold. Then V has a normally hyperbolic critical manifold, carrying a reduced slow flow that coincides on x with the Filippov flow.
The result is illustrated in Fig. 4, see figure caption for further details. In practice, the analysis of V is performed in directional charts. These charts correspond to settinḡ y = ±1 and¯ = 1 in such a way that (12) takes the following forms: . Illustration of blowup in the case of regularization by smoothing. Upon blowup we gain smoothness and hyperbolicity along the edges of the cylinder (indicated by the double-headed arrows) as well as a critical manifold (in pink) when the associated PWS system has stable sliding. The most fundamental result, see Proposition 2.1, is then that the slow-flow on this critical manifold is given by Filippov's sliding vector-field.
The blowup approach has also been used to study bifurcations of systems of the form (3), including boundary equilibrium bifurcations where equilibria of either Z ± collide with Σ upon parameter variation [10,11]. For this, the blowup transformation (12) is combined with additional blowup transformations in order to deal with lack of hyperbolicity. In this paper, we will adapt the blowup approach to study (6).

2.2.
A blowup approach for (6). To study (6) with α > 0, we now proceed in a similar manner. First, however, due to the time scale separation of (6), we introduce a fast time and augment trivial equations for and α ≥ 0: (14) x = αX(z, p), (6) is PWS with respect to both → and α → 0, we anticipate that we will need to perform two blowup transformation. In light of this, Section 2.1 suggests that we should consider (14) with respect to an even faster time-scale, corresponding to multiplying the right hand side by α again. But notice, despite the similarities, there is also a fundamental difference between (14) and (11) insofar that the discontinuity set of (14) for , α → 0 is y = 0, x ∈ Σ, p ∈ R, but the discontinuity only enters the p-equation. To avoid too many multiplications and subsequent divisions by the same quantities, we will therefore proceed more ad hoc in the following; in fact, the analysis will show that it is only necessary to multiply the right hand side of (14) by in order gain smoothness. Apriori it is not obvious how the two blowup transformations should be organized and whether the order is important, but leaving and α as independent small parameters, we will show that it is convinient to first blowup with respect to α. (See the end of the section for a further discussion of this.) We therefore first apply the following transformation where r ≥ 0, (ȳ,ᾱ) ∈ S 1 , leaving all other variables -including -untouched. In this way, we gain smoothness with respect to α ≥ 0 for any > 0. Indeed, the transformation (15) gives a smooth vector-field V for > 0 on (x, p, , (r, (ȳ,ᾱ)), with r ≥ 0, (ȳ,ᾱ) ∈ S 1 ,ᾱ ≥ 0, by pull-back of (14) (without the need for further transformation of time). Notice specifically, that using (15) the first term in the p-equation in (14) Figure 6. Illustration of the two consecutive blowup transformations relating to (6). The first cylinder corresponds to (15). The second one corresponds to (17). becomes: which for each > 0 is smooth on (ȳ,ᾱ) ∈ S 1 ∩ {ᾱ ≥ 0}. However, there is still a lack of smoothness along (ᾱ,ȳ) = (1, 0) as → 0. To deal with this, we perform a second blowup transformation: where ν ≥ 0, (ȳ,¯ ) ∈ S 1 . Indeed, in this way, (16) becomes regular under assumption 4. We illustrate the blowup transformations in Fig. 6. As described in Section 2.1 in the context of regularization by smoothing, we will also use different directional charts in the analysis of (6) to cover the two cylinders. In particular, to cover the first cylinder, defined by (15) and (ȳ,ᾱ) ∈ S 1 , we (re-)consider the two charts defined by: We will refer to these charts by (ȳ = 1) 1 and (ᾱ = 1) 2 , respectively, henceforth. In principle, we will also need the chart (ȳ = −1) 3 that coversȳ < 0 of cylinder, but since the analysis there is identical to the analysis in the (ȳ = 1) 1 -chart we skip this. The change of coordinates between the charts (ȳ = 1) 1 and (ᾱ = 1) 2 are given by the expressions: Subsequently, to cover the second cylinder due to (17), we notice that in (ᾱ = 1) 2 , (17) becomes for ν ≥ 0, (ȳ,¯ ) ∈ S 1 . Therefore we define the following charts (ᾱ = 1,ȳ = 1) 21 : (ᾱ = 1,¯ = 1) 22 : In both charts, we have α = r 2 . (The chart corresponding toȳ = −1 is again similar tō y = 1 and therefore left out.) The change of coordinates are given by the expressions valid for y 22 > 0.
The two blowup transformations relate to two important scalings. Firstly, in the (ᾱ = 1) 2 -chart, we have (24) y = −αp + αy 2 , upon eliminating r 2 and consequently Through the coordinate y 2 , we therefore zoom in on a O(α)-neighborhood of y = 0. From (16), we understand that the resulting vector-field V 2 in terms of (x, y 2 , p, r 2 , ) is itself PWS in the limit → 0. Consequently, following Section 2.1 and the results for gaining smoothness of (11), we see that through (17), we obtain a smooth vector-field V 2 on (x, p, α, ν, (ȳ,¯ )), ν ≥ 0,¯ ≥ 0, by pullback of V 2 . This system has¯ as a common factor and it is therefore V 2 :=¯ −1 V 2 that we will study.
Next, we emphasize that in the (ᾱ = 1,¯ = 1) 22 -chart, we have y 2 = y 22 upon eliminating and consequently Therefore we also have that and we see that coordinate y 22 provides a zoom on a O(α )-neighborhood of y = −αp.
It is obvious that the scaling defined by (26) is important; this captures the region where the first term in the p-equation (14) changes by an O(1) amount with respect to , α → 0. It also seems reasonable that the scaling (24) is useful, but it not obvious why the scaling defined by seem to play no role. To see this we have to insert this expression into (14). This gives Here we would like to divide by on the left hand side, but for this we will have to make assumptions on relative to α (i.e. whether −1 α is small, moderate or large). If we insert (24) instead, then we obtain Here α is a common factor on both sides which can therefore be divided out. This explains why (24) and (26) are both important in our analysis and why (28) will not be used.
2.3. Notation. Throughout the paper we follow the convention that a set S in the blowup space is given a subscript when viewed in a chart. I.e. the subset of a set S, which is visible in the chart (ȳ = 1) 1 , will be called S 1 . Similarly, S 2 in the chart (ᾱ = 1) 2 . In the charts, (ᾱ = 1) 2 and (ᾱ = 1,¯ = 1) 22 , r 2 = α and ν 22 = are constants, so when working in these charts, it is most convinient to eliminate r 2 and ν 22 , respectively, and return to treat and α as parameters. The only important thing to keep in mind in regards to this, is that when we change coordinates (e.g. through (20) and (23)) then this has to be viewed in the appropriate space. For example, in the (ᾱ = 1) 2 -chart, we will obtain a slow manifold S ,α,2 in the (x, y 22 , p)space. When writing this in the (ȳ = 1) 1 -chart, we first have to embed S ,α,2 in the extended (x, y 22 , p, , α)-space in the obvious way. We can then apply the change of coordinates (20) with r 2 = α and obtain S ,α,1 . We will henceforth perform similar change of coordinates without further explanation, moving back and forth between different spaces, treating and α as parameters whenever it is convinient to do so.

Main results in the case of stable sliding
In this section, we will use the blowup approach, outlined in the previous section, to describe the dynamics of (6) under the assumption 4 of stable sliding. More specifically, we will provide a detailed study of the dynamics in each of the charts (ȳ = 1) 1 , (ᾱ = 1) 2 , (ᾱ = 1,ȳ = 1) 21 , (ᾱ = 1,¯ = 1) 22 . In summary, this analysis reveals the existence of two critical manifolds C and M . Whereas C extends onto the first blowup cylinder, obtained by (15), M lies on the subsequent blowup cylinder, obtained by (17). Moreover, C is normally attracting and enables an extension of S ,α in Lemma 1.1 up to y = cα, for c > 0 and , α > 0 small enough. On the other hand, M is normally repelling. Using the geometric representation used in Fig. 6, we illustrate the findings in Fig. 7. On both C and M , we obtain a desingularized slow flow; the direction of this flow is also indicated in the figure but we emphasize that x (not shown) is a constant for this reduced flow. This leads to a singular cycle Γ x for each x ∈ Σ, which we indicate in Fig. 7 using curves of increased thickness. Due to the desingularization along C, Γ x is akin to a relaxation cycle in slow-fast systems. In the full blowup space, the curves Γ x make up a singular cylinder Γ = {Γ x } x∈Σ of dimension n + 1. The first main result basically says that this singular cylinder persists for 0 < , α 1 and that this manifold carries a reduced flow, which can be related to the Filippov sliding vector-field (2). Theorem 3.1. Suppose that assumptions 1-4 all hold true and let K > 1. Then there exists a δ > 0 such that for any 0 < , α < δ, ( 6) has an invariant cylinder C ,α of dimension n + 1, contained within y ∈ (−Kα, Kα). C ,α is uniformly Lipschitz in the blowup space and converges to Γ in the Hausdorff-distance as , α → 0. Figure 7. Illustration of the dynamics on the blowup system. Our analysis reveals two normally hyperbolic critical manifolds C and M . In case of stable sliding, the reduced slow flow on these invariant manifolds reveals closed singular cycles Γ x for each x ∈ Σ. This cycle does not have completely desirable hyperbolicity properties due to the degeneracy at the point Q.
Let Π 0 be a local section on y = −αp transverse to C ,α and define x → x + (x, , α) and x → T (x, , α) to be the corresponding return map and the transition time, respectively. Then where the order of the remainder remains unchanged upon differentiation with respect to x. Specifically, lim ,α→0 Theorem 3.1 generalizes Proposition 2.1 to the framework of (6) with 0 < , α 1 without any order dependency on and α.
We prove Theorem 3.1 in the following. In Section 3.1-Section 3.4, we first analyze the dynamics in each of the charts (ȳ = 1) 1 , (ᾱ = 1) 2 , (ᾱ = 1,¯ = 1) 22 , (ᾱ = 1,ȳ = 1) 21 , respectively. In Section 3.5, we then collect the findings in the local charts into a global result, see Fig. 7. This includes a detailed description of Γ x . Following this in Section 3.6, we first present a description of the return map defined on the section Π 22 : y 22 = 0 transverse to Γ in the (ᾱ = 1,¯ = 1) 22 -chart, see Lemma 3.8. The description of this mapping rests upon a subsequent blowup transformation of the degenerate point Q, which sits at the interface between C and M , with the purpose of gaining hyperbolicity. The details of this blowup analysis of Q and the proof of Lemma 3.8 are delayed to Section 4. Before this in Section 3.7, we show how Lemma 3.8 implies Theorem 3.1. Here we rely on a general result [24, Theorem A.1] on the existence of an invariant curve for a return mapping.
3.1. Analysis in the (ȳ = 1) 1 -chart. In this chart, we insert (18) into (14) and obtain (31) , and = 0, using assumption 2. This system is the local form of V in the (ȳ = 1) 1chart. As already advertised above, we will henceforth treat as parameter in this chart. In (31), we have defined with y = −αp + r 1 and α = r 1 α 1 on the right hand sides. The system (31) is a slowfast system in nonstandard form with respect to the small perturbation parameter . Indeed for = 0, the set C 1 defined by p = 1 is a critical manifold of the layer problem: Fig. 8. The linearization around any point in C 1 produces −1 as the only nonzero eigenvalue. C 1 is therefore normally attracting.
Lemma 3.2. Consider any compact submanifold S 0,1 of C 1 , defined as the graph p = 1 over a compact domain D 1 in the (x, r 1 , α 1 )-space. Then for all 0 < 1, there exists a locally invariant slow manifold S ,1 , which is also a smooth graph over D 1 : For any α > 0 small enough, we let S ,α,1 denote the constant α-section, defined by α = r 1 1 , of the center manifold S ,1 . The resulting invariant manifold S ,α,1 provides an extension of the slow manifold S ,α in Lemma 1.1 into the (ȳ = 1) 1 -chart.
On S ,1 , we have a reduced flow defined by , upon desingularization, corresponding division of the right hand side by α 1 .  . Dynamics in the (ȳ = 1) 1 -chart. The manifold C 1 is normally hyperbolic. C 1 actually extends to any α = r 1 α 1 but in this picture we illustrate the α = 0 limit. Figure 9. Reduced dynamics on the normally hyperbolic critical manifold C 1 (blue in Fig. 8) in the case when Y + (x, 0) < 0. Within α 1 = 0 the system is equivalent to z = Z + upon time reparametrization for r 1 = y > 0. The line r 1 = α 1 = 0 is normally hyperbolic, each point having stable and unstable manifolds (green and orange, respectively) under the assumption Y + (x, 0) < 0. In particular, the former invariant manifold lies within r 1 = 0, and along this set x is a constant.
The reduced problem is illustrated in Fig. 9. Notice it is identical to what is found by smoothing the PWS system, recall (11) and Fig. 4, near the edge of the blowup cylinder defined by (12).

3.2.
Analysis in the (ᾱ = 1) 2 -chart. In this chart, we insert (18) into (14) and obtain the following equations with = α = 0. Within = 0, we re-discover the manifold of equilibria C 1 from the (ȳ = 1) 1 -chart, in the following form: Notice that the dependency on α is regular. In particular, note that C 2 is a critical manifold for any α ≥ 0. We will often view it within α = 0 (as in Fig. 8 since α = r 1 α 1 in the (ȳ = 1) 1 -chart).
The manifold C 2 is also normally attracting for (33) and carries the following reduced problem upon passage to the slow time for = α = 0.
In further details, let S 0,α,2 ⊂ C 2 be a compact submanifold contained within y 2 > 0 for any α ≥ 0. Then S 0,α,2 perturbs to a slow manifold S ,α,2 by Fenichel's theory for 0 < 1 and an easy calculation shows that it takes the following graph form: We may choose S ,α,2 to coincide with S ,α,1 where these overlap upon change of coordinates. The reduced problem on S ,α,2 is given by

Analysis in the
In this chart, we insert (22) into (14) and obtain the following equations: and = α = 0, upon multiplication of the right hand side by . The system (35) is now a slow-fast system with respect to ≥ 0 in standard form, x and p being slow while y 22 is fast. For = 0, we obtain the following layer problem: (36)ẋ = 0, and consequently the set M 22 defined by (x, y 22 , φ(y 22 )) is a critical manifold, even for α > 0. As with C, we will often think of M 22 within α = 0. The manifold M 22 is normally repelling, since the linearization of (36) around any point (x, y 22 , φ(y 22 )) produces φ (y 22 ) > 0 as a single nonzero eigenvalue, see assumption 1.
The reduced problem N ,α,22 is given by in terms of a slow time (that corresponds to dividing the right hand side by 2 ).
Proof. For the reduced problem, we first use thaṫ on N ,α, 22 . Then upon realizing that p = φ(y 22 )+O( ), we obtain the desired result.
Notice that for = α = 0, we can also write (37) as which is more convinient. It is normally repelling (attracting) if the sliding is stable (unstable). The reduced problem on K 22 is given byẋ recall ( 2), with respect to the original (slow) time for = α = 0.
Proof. Follows from a simple calculation and the details are therefore left out.
In Fig. 10 we summarize the findings. 3.4. Analysis in the (ᾱ = 1,ȳ = 1) 21 -chart. In this chart, we obtain the following equations: , upon desingularization through division of the right hand side by 21 . Here we treat α as parameter and have introduced the following quantities  The set B 21 defined by ν 21 = 21 = 0 is a set of equilibria for any α ≥ 0. The linearization about any point in this set has two nontrivial eigenvalues: ±(1 − p). Consequently, the subset Q 21 ⊂ B 21 defined by p = 1 is fully nonhyperbolic, also for any α ≥ 0.
In the present chart, the subset of M 22 with y 22 > 0 becomes M 21 defined by (23). The corresponding graph therefore ends in Q 21 for 21 = 0. There is obviously another critical set C 21 , given by 21 = 0, p = 1, ν 21 > 0, emanating from Q 21 . It corresponds to C 1 from the (ȳ = 1)-chart, see Section 3.1. Both sets, M 21 and C 21 are normally hyperbolic, M 21 being repelling whereas C 21 is attracting; the set Q 21 therefore acts like a fold point where attracting and repelling critical branches of critical manifolds come together.
Let 21 = 0 in (39):ẋ = 0, It follows that each point on the critical set (x, 0, p, 0) ∈ B 21 with p < 1 is connected by a heteroclinic orbit through the dynamics of (39) to a point on C 21 . In particular, we have the following result, which follows from a simple calculation.

3.5.
Collecting the local results into a global picture. Fig. 7 summarizes the findings in the local charts. Notice specifically, that while we have focused on the upper part of the cylinders, the analysis of the lower part is identical and therefore skipped. In conclusion, we obtain a singular cycle Γ x for each x ∈ Σ sl , x being a constant on the two cylinders. Γ x is the union of six pieces γ xi , i = 1, . . . , 6 where: (1) γ x1 is a heteroclinic connection on the first cylinder. It is described in the coordinates of the (ᾱ = 1,ȳ) 21 -chart in Lemma 3.7 (corresponding to p = 0 in this result). In particular, x is constant along γ x1 and its α-limit set is given by (ν 21 , p, 21 ) = (0, 0, 0) on B 21 , whereas the ω-limit set is given by (ν 21 , p, 21 ) = (1, 1, 0), belonging to the normally attracting set C 21 .
Lemma 3.8. The mapping P 22 is given by The remainder terms remain unchanged upon differentiation with respect to x and p.  Upon applying the flow map to the invariant curve of P 22 in Proposition 3.9, we obtain the desired invariant cylinder C ,α in Theorem 3.1. To finish the proof of Theorem 3.1, we just have to prove (29). For this, we reduce the mapping P 22 to the invariant manifold to C ,α . From the previous analysis, we obtain This completes the proof of the expression for x + in (29). The expression for the transition time is similar; in fact, it can be obtained from the expression for x by setting X + = X − = 1 (sinceṫ = 1).

Proof of Lemma 3.8
To prove Lemma 3.8, we will chop the return map P 22 into several local pieces. However, to describe the local transition near the degenerate set Q, we have to perform an addition blowup step. In the following, we first analyze this blowup transformation and the associated dynamics in separate local charts. In this way, we obtain singular cycles Γ x with improved hyperbolicity properties.

21
= ρ 213 ν −1 213 . The advantage of working with this chart, is that in these coordinates and ρ 213 is therefore conserved. In comparison, we have in the other charts. Notice that we also haveν 21¯ 213 , which is why we only use these coordinates to cover a compact subset ofν 21 ,¯ 21 > 0. The coordinate changes between the different charts are given by the following expressions: 4.2. Entry chart (ᾱ = 1,ȳ = 1,ν 21 = 1) 211 . In this chart, we obtain the following equations: (45)ẋ = ρ k+1 211 211 αX 211 (x, ρ 211 , p 211 , α), Consequently, (x, 0, 0, 0) and (x, 0, 1, 0) are both partially hyperbolic. The former allows us to extend the critical manifold C 21 in chart (ᾱ = 1,ȳ = 1) 21 onto the blowup sphere as a normally hyperbolic invariant manifold C 211 . In fact, within ρ 211 = 0 we have that p 211 = −β k 211 is a manifold of equilibria R 211 and C 211 will therefore include these points, at least locally. We will see the resulting slow-fast structure more clearly in the chart (ᾱ = 1,ȳ = 1,ν 21¯ 21 = 1) 213 which we analyze in the following section. The hyperbolicity of C 211 allows us to extend the slow manifold S ,α as a constant -section S ,α,211 , defined by (43), of a center manifold S α,211 .
Lemma 4.2. The transition map P 4 211 from Π 4 11 to Π 5 11 takes the following form for some c > 0. The order of the remainders remain unchanged upon differentiation with respect to x andp 211 .
Proof. The proof is standard using Fenichel's theory and normal forms, see e.g. [12].
In particular, sincep 211 = 0 is invariant, we havė , upon dropping the tildes. Then upon invoking Fenichel's normal form [12], we straighten out the stable fibers of S α,211 by settingx = x + O(ρ k+1 211 211 α). Then the (x, ρ 211 , 211 )system is independent of p 211 and described by (46) upon dropping the tilde. We then simply integrate the ρ 211 and 211 -equations in (46), insert the resulting expressions into the x-equation and estimate x. On the other hand, on the time scale of (46), p 211 is given byṗ   To the left, we illustrate the dynamics in the (ᾱ = 1,ȳ = 1,ν 21¯ 21 = 1) 213 -chart using a projection onto the (p 213 , ν 213 )-coordinate plane. Here we find a critical manifold R 213 which has a regular fold jump point (in green). On the left, we summarize the findings on the blowup of Q. The local diagram on the left covers the subset of the sphere that is bounded away from the edges (purple).
Proof. The statement about the jump points, follow from an analysis of the reduced problem on R 213 : This can be obtained from [25] or more directly by writing the slow manifold approximation as where Y 213 (x, ν 213 , −βν −k , 0, 0) = Y + (x, 0), inserting the result into the (x, ν 213 )-subsystem, writing the system in terms of the slow time and then letting ρ 213 → 0.
Proof. We solve (52) for ν 212 and ρ 212 , so that ν 212 (t) = c in e −(k+1)t , ρ 212 (t) = e t ρ 2,in and define u 2 (t) by p 212 (t) = c in e −(k+1)t + (u 2 (t) − c in )e −kt . Inserting this into the p 212 equation giveṡ The transition time is T = log c out ρ −1 2,in . Notice that quantities Z 212 (· · · ), X 212 (· · · ) are uniformly bounded on this domain. By integrating the equations we therefore obtain Recall that = ρ k+1 212 ν 212 in this chart. 4.5. Completing the proof of Lemma 3.8. In Fig. 13, we summarize the findings from our analysis of the two cylindrical blowups and the blowup of Q. In particular, the blowup of Q gives rise to an improved singular cycle.
In Fig. 13, we also indicate different sections Π i , i = 1, . . . , 8, that are each transverse to Γ, that we use to decompose the return mapping P 22 in Lemma 3.8. We describe each of the local mappings Π i−1 → Π i , i = 1, . . . , 8 in the following. We try to stroke the balance between including a complete, rigorous analysis while avoiding an unnecessary amount of details, that can be found elsewhere in a similar context. Π 0 → Π 1 . The transition from Π 0 and Π 1 is regular in the (ᾱ = 1,¯ ) 22 -chart. We therefore leave out further details.
On the other hand, the transition map from Π 1 to Π 2 is described in the coordinates (x, ν 21 , p, 21 ) of the (ᾱ −1ȳ = 1,ᾱ = 1) 21 -chart. We therefore consider (39) and define the sections as follows Π 1 21 : 21 = c in , p ∈ I in to Π 2 21 : ν 21 = c out , p ∈ I out , with I in and I out open neighborhoods of p = 0. Notice, for these values of p, the set B 21 is normally hyperbolic, see Fig. 11.
Proof. The proof is standard, see e.g. [4, Proposition 2.1]. Π 2 → Π 3 . The transition map from Π 2 → Π 3 is regular in the (ᾱ = 1,ȳ = 1) 21 -chart and further details are therefore left out. Π 3 → Π 4,5 . The transition map from Π 3 to Π 4 is obtained from Fenichel's theory near the normally attracting manifold C e.g. by working in the (ᾱ = 1) 2 -chart. In fact, by working in chart (ᾱ = 1,ȳ = 1) 21 and using the blowup transformation (41) this result can be extended all the way up to the section Π 5 on the blowup of Q 21 . The details were given in Lemma 4.2.
Π 5 → Π 6 . The transition map from Π 5 to Π 6 is best described in the chart (ᾱ = 1,ȳ = 1,ν 21¯ 21 = 1) 213 where the equations are slow-fast. The transition map is then given as a regular fold, jump set with ρ 213 = as the small parameter. See e.g. [24] for further details.
The exit from the blowup sphere, that we describe by a transition map from Π 6 to Π 7 is given by the transition near a resonance saddle. The details were given in Lemma 4.4. Π 7 → Π 8 . The transition map from Π 7 → Π 8 is regular in the (ᾱ = 1,¯ = 1) 22 -chart and further details are therefore left out.

Main results in the case of grazing
In this section, we consider (6) under the following assumption (which replaces assumption 4 henceforth): Assumption 5. The PWS system Z ± is planar z = (x, y) ∈ R 2 and each Z ± depends smoothly on an unfolding parameter µ ∼ 0. In particular, for µ = 0, Z + has a hyperbolic and repelling limit cycle γ 0 that has a quadratic tangency with Σ at x = 0. Z − , on the other hand, is assumed to be transverse to Σ.
Since γ 0 is hyperbolic for Z + , we have a repelling limit cycle γ µ of Z + for every µ ∼ 0. Let Y (µ) = min t y(t) along γ µ so that Y (0) = 0 by assumption 5. We assume the following degeneracy condition.
We illustrate the setting in Fig. 14. Under these assumptions, the reference [13] proved that the system obtained from regularization by smoothing (3), has a locally unique saddle-node bifurcation of limit cycles at µ = o(1) with respect to → 0. On the other hand, the reference [21] also showed that the system obtained from regularization by hysteresis has chaotic dynamics (through a Baker-like map) for all α > 0 sufficiently small provided µ ∼ 0 is sufficiently small. In this section, we try to bridge these two results by working on (6).
To present the result, we define two wedge-shaped regions in the ( , α)-plane. Firstly, for 0 > 0, α 0 > 0, let W 1 ( 0 , α 0 ) be the region defined by 0 < α ≤ 2k α 0 for 0 < ≤ 0 . On the other hand, let W 2 ( 0 , 1 , α 0 ) be the region defined by 1 , for 0 < α ≤ α 0 and 0 < 0 < 1 . We illustrate the two regions in Fig. 15. Figure 14. The grazing bifurcation. We assume that the smooth vector-field Z + has a repelling limit Γ 0 for µ = 0 having a quadratic tangency with Σ. Under a further degeneracy condition, which ensures that the perturbation Γ µ of Γ 0 as a limit cycle of Z + for µ ∼ 0 transverses Σ with nonzero speed, see assumption 6, the reference [21] have shown that, while regularization by smoothing leads to a saddle-node bifurcation of limit cycle [13], regularization by hysteresis leads to chaotic dynamics. Theorem 5.1 is an attempt to bridge these two regimes by working on (6). In the following, we will sometimes write W 1 ( 0 , α 0 ) and W 2 ( 0 , 1 , α 0 ) as W 1 and W 2 for simplicity.
(2) For any ( , α) ∈ W 2 ( 0 , 1 , α 0 ), there exists a µ ∼ 0 such that there is a return map defined by the system ( 6) having an invariant cantor set upon which the map is homeomorphic to a full shift on n symbols.
To prove the theorem, we have to describe the local transition near the grazing point with Σ. By assumption 5, this point is a visible fold point [8,13], see T in Fig. 3. In fact, by the implicit function theorem, Z + has visible fold point for each µ ∼ 0 and this point depends smoothly on µ. Then using [1], we transform the system locally into (6) with Here f and g are smooth functions with f (0) = 0 that also depend smoothly on µ. This is the system that we will use to study the local dynamics; the global dynamics, away from Σ, is regular on an attracting slow manifold by Lemma 1.1. In our analysis, we can reuse some of the results obtained in the proof of Theorem 3.1. We still have the existence of C, M and R on the blowup of the degenerate point Q, recall Fig. 13. The essential differences lie in the reduced problems on these manifolds.
Before going into details, we first emphasize that Z ± in (54) has stable sliding for x < 0 and crossing for x > 0 along Σ. Therefore the blowup dynamics for x < 0 in a compact interval is covered by Theorem 3.1 and the blowup dynamics for = α = 0 is therefore as in Fig. 13 in this case. The blowup dynamics for x > 0 on the other hand, where assumption 4 is violated and crossing occurs, is shown in Fig. 16. This follows from the blowup analysis with Y + > 0. In each of the two diagrams, Fig. 13 and Fig. 16, x is constant on the cylinders and there is only slow flow in theȳ-direction. In order to describe the details of the dynamics associated with the visible fold, we will need to zoom in on x = 0 so that the dynamics in this direction for 0 < , α 1 becomes comparable with the dynamics in theȳ-direction. We achieve this zoom through blowup. In particular, in Section 5.1, we first reduce to the slow manifold S ,α obtained as a perturbation of the critical manifold C in the (ȳ = 1) 1 -chart and then perform two separate blowup transformations. In the parameter regime ( , α) ∈ W 1 , this is sufficiently to prove Theorem 5.1 (1). The details are similar to those in [13] covering the grazing bifurcation in the case of regularization by smoothing.
In order to prove Theorem 5.1 (2) in the regime ( , α) ∈ W 2 we have to follow dynamics that becomes unbounded in the chart (ȳ = 1) 1 . In Section 5.3, we will specifically work on the blowup of Q. Here we will study the reduced problem on the critical manifold R 213 in the (ᾱ = 1,ȳ = 1,ν 21¯ 21 = 1) 213 -chart for x ∼ 0 using a separate blowup transformation. This gives rise to a folded saddle and a canard orbit for ( , α) ∈ W 2 along which the extended versions of the slow manifolds S ,α and N ,α of C and M , obtained in chart (ȳ = 1) 1 and (ᾱ = 1,¯ = 1) 22 , respectively, intersect transversally. This provides the main mechanism for the chaotic dynamics in Theorem 5.1 (2). In fact, the geometric construction is similar to [14], which proved existence of chaos in a friction oscillator in the presence of slow-fast and nonsmooth effects. We therefore complete the proof of Theorem 5.1 (2) in Section 5.4 by exploiting this connection.

5.1.
Analysis of the slow flow on S ,α in the case of the visible fold. In this section, we work in the (ȳ = 1) 1 -chart and consider the reduced flow on S ,1 , recall Figure 16. Blowup dynamics in the case of crossing upwards where Y + (x, 0) > 0, Y − (x, 0) > 0 (corresponding to x > 0 for the visible fold in Fig. 3). In this case, the flow on C (blue) moves upwards and there is no equilibrium of the reduced flow on M (moving downwards. Lemma 3.2, in the case of (54). For this, we use (32) with Z ± as in (54): The dynamics of this system within the invariant subspaces α 1 = 0 and r 1 = 0 are illustrated in Fig. 17. Notice that r 1 = x = 0, α 1 ≥ 0 is a line of degenerate singularities for = 0. We will again need to perform consecutive blowup transformations to resolve the degeneracy stemming from the terms of the form k α k 1 . For this, we first blowup with respect to α 1 and then subsequently blowup with respect to . In further details, we first apply the transformation: Since r 1 , α 1 ≥ 0 we are only interested in the quarter sphere wherer 1 ,ᾱ 1 ≥ 0. This gives a vector-field V 1 on σ ≥ 0, (x,r 1 ,ᾱ 1 ) ∈ S 2 by pull-back of (55). It has σ k as common factor and it is therefore V := σ −k V that has improved hyperbolicity properties. The analysis of V is straightforward in the directional charts. In particular, in the chart defined by (ȳ = 1,r 1 = 1) 11 : where α = σ 2k+1 11 α 11 , we obtain the following equations (57)ẋ where [· · · ] = (2x 11 + σ k 11 g 11 (x 11 , σ k 11 ))(1 + O( k σ k 11 α k 11 )) + β k α k 1 + O( k σ 11 α k 11 ), with σ k 11 f 11 (x 11 , σ k 11 ) := f (x, r 1 ), and g 11 (x 11 , σ k 11 ) := g(x, r 1 ). In particular, we find two hyperbolic equilibria: for any ≥ 0. The eigenvalues of the linearization around these points are with x 11 = ±1 at the two points q ± 11 , respectively. Whereas the point q − 11 has a two-dimensional unstable manifold within σ 11 = 0, and a one-dimensional unstable manifold within α 11 = 0 (corresponding to the grazing orbit of the PWS system (54) within x < 0), the point q + 11 has a two-dimensional stable manifold within σ 11 = 0 and a one dimensional unstable space within α 11 = 0 (corresponding to the grazing orbit of the PWS system (54) within x > 0). Compare also Fig. 17 and Fig. 18. The dynamics on the sphere is given by (60)ẋ 13 = 1, r 13 = 2x 13 , upon using the coordinates (x 13 , r 13 ) defined by (ȳ = 1,r 1ᾱ1 = 1) 13 : where α = r 1 α 1 = σ 2k+1
To gain hyperbolicity and resolve the dynamics near the degenerate line (pink in Fig. 18), we proceed to augment˙ = 0 and then apply the following cylindrical blowup transformation: leaving σ 12 untouched. Let V 12 be the vector-field in the (ȳ = 1,ᾱ 1 = 1) 12 -chart witḣ = 0 augmented. The blowup transformation (61) then gives a vector-field V 12 by pull-back of V 12 . It has ξ k as a common factor and it is therefore V 12 := ξ −k V 12 that we shall study.
The eigenvalues of the linearization around these points are given by with x 121 = ±1 at the two points z ± 121 , respectively. The unstable manifold for z − 121 is three-dimensional and contained within ξ 121 = 0. However, for z − 121 it is the stable manifold that is three dimensional; in fact, z + 121 will be the ω-limit set of all points with ξ 121 = 0, 121 = 0. Notice, that since To describe the dynamics in further details, we focus on the cylinder ξ = 0, (x 12 ,r 12 ,¯ ) ∈ S 2 , σ 12 ≥ 0, and the two invariant subspace of V 12 | ξ=0 given by σ 12 = 0 and¯ = 0. The reason for doing so, is that these invariant spaces capture different scaling regimes of and α. In particular, within the (ȳ = 1,ᾱ 1 = 1,r 12 = 1) 121 -chart, (63) holds and on σ 12 = const. we therefore have by (66) that 2k ∼ α 2k 121 .
(Here we have used ∼ to indicate that two quantities differ by a constant that only depends upon the constant value of σ 12 .) Consequently, the regular dynamics follow close to¯ = 0 (i.e. 121 = 0) provided that Notice also that on σ 12 = const. we have upon eliminating ξ 121 . This will be important later on.
Dynamics of V 12 | ξ=0 in the invariant subspace σ 12 = 0. In the (ȳ = 1,ᾱ 1 = 1,¯ = 1) 122chart, we obtain the following local form of V 12 within ξ 122 = σ 12 = 0: Within r 122 = 0 we find two equilibria, one given by x 122 = 0 and another given by The first point is hyperbolic and repelling for (72) whereas the second one is partially hyperbolic, the linearization having a single nonzero and negative eigenvalue. A simple calculation, reveals the following: There exists a unique, attracting center manifold G 122 for ( 72) of the point (x 122 , r 122 ) = (− β 2 , 0). G 122 is its (nonhyperbolic) unstable manifold, along which r 122 is increasing.
We summarize the findings in the two charts in Fig. 21.

5.2.
Proof of Theorem 5.1 (1). For the proof Theorem 5.1 (1) we work on the slow manifold S ,α that has been extended, through the blowup approach in Section 3, to the first blowup cylinder. On this manifold, using the (x, y)-coordinates and the system (54) locally near (x, y) = 0, we then consider the return map P on a section Π in = {(x, y) : y = c in , x ∈ I in }, for some appropriate closed interval I in ⊂ (−∞, 0) so that Π in is transverse to γ 0 . We then decompose P into a local transition map P loc : Π in → Π out , with Π out = {(x, y, p) : y = c out , x ∈ I out }, see Fig. 3, and a global map P glo : Π in → Π out . The latter is regular on the attracting slow manifold, and we therefore turn our attention to P loc .  In order to describe P loc , we use the chart (ȳ = 1) 1 and the blowup transformations (56) and (61), that resolve the degeneracy of x = r 1 = 0, α 1 ≥ 0 for = 0, and chop the mapping into separate transition maps, see Fig. 22: P 0 : Π 0 → Π 1 near q − , a regular map P 1 : Π 1 → Π 2 being a regular perturbation of (60), P 2 : Π 2 → Π 3 near z − , a regular map P 3 : Π 3 → Π 4 being a perturbation of the map in Lemma 5.4, P 4 : Π 4 → Π 5 near z + , a regular map P 5 : Π 5 → Π 6 being a regular perturbation of (60), and finally P 6 : Π 6 → Π 7 near q + .
Although the eigenvalues near the points q ± , z ± are resonant, it is possible, following [13], to achieve a (suitable) linearization near each of this points. We will only present the details near q − and z − . Figure 21. The reduced flow on C 1 in the (ᾱ = 1) 1 -chart upon two consecutive blowup transformations of the degenerate set α 1 ≥ 0, r 1 = x = = 0. The dynamics on the cylinder obtained by the blowup transformation (61) (its boundary¯ = 0 being indicated in pink) breaks up into different regimes, depending on the ratio of and α. For example, whenever ( , α) ∈ W 2 then the dynamics near¯ = 0 (pink) becomes relevant, whereas within ( , α) ∈ W 1 the green region where¯ > 0, described by the Chini-equation (74), becomes relevant. In this region, which is more visible in Fig. 22, the attracting center manifold G produces a contraction -which is absent for ( , α) ∈ W 2 , see Lemma 5.2 -of the return map P loc , see Lemma 5.4. It is the balance of this contraction and the expansion along γ 0 that gives rise to the saddle-node bifurcation in Theorem 5.1 (1). Local transition map near q − . Consider (57) and divide the right hand side − 1 2 [· · · ], which is ≈ −x 11 and therefore positive near q − 11 . This gives (76)ẋ 11 = x 11 − 2(1 + σ k 11 f 11 (x 11 , σ k 11 )) 2x 11 + σ k 11 g 11 (x 11 , σ k 11 ) + k α k 11 A 11 (x 11 , σ 11 , α 11 , ), for A 11 smooth.
Proof. The proof can be found in [13], see Lemma 3.5 and Lemma 3.6 in this reference, but essentially we use that the α 11 = 0 subsystem is equivalent to z = Z + (z) which is regular. This enables a linearization within α 11 = 0 through the flow box theorem. Subsequently, we linearize the non-resonant system within σ 11 = 0. Consider (77) and notice that α =σ 2k+1 11α 11 is still conserved in the tilde variables. We therefore drop the tildes and describe the transition map P 0 11 from Π 0 11 : σ 11 = c in to Π 1 11 : α 11 = c out by integrating these equations. This produces the following result. with c > 0 fixed small enough and given by ( The order of the remainder terms does not change upon differentiation with respect to x 11 .
The analysis near q + is almost identical. In particular, although the local mapping near q − is expanding, the local mapping near q + contracts by the same order. with c > 0 and α 0 small enough and given by (x 121 , c in , σ 12 12σ (x 121 , σ 12 , 121 )). The order of the remainder terms does not change upon differentiation with respect to x 121 . Moreover, by ( 80) Proof. Simple calculation.
The analysis near z + is almost identical. In particular, although the local mapping near z − is expanding, the local mapping near z + contracts by the same order.
The local map P loc . Let x in denote the value of x on Π 0 = Π in of the grazing orbit of Z + . Similarly, let x out be the corresponding value on Π 7 = Π out . ( , α) ∈ W 1 ( 0 , α 0 ) implies that in the (ȳ = 1,ᾱ 1 = 1,r 12 = 1) 121 -chart and it is therefore consistent with (80). Consequently, by Lemma 5.7 and Lemma 5.9, we consider any ( , α) ∈ W 1 ( 0 , α 0 ) with α 0 > 0 small enough and x in a small neighborhood of x in : , for some c > 0. This leads to the following.
The following can be said about φ ± : For any n ∈ N, we can take φ ± as close to 1 as desired and φ (k) ± as close to 0 as desired for any k = 2, . . . , n, by adjusting the domains of the mappings (i.e. by adjusting c in , c out , c > 0).
Moreover, the remainder term o(1) is bounded by a constant c m (α 0 ) → 0 for α 0 → 0 in C m , m ∈ N fixed.
Proof. The proof is similar to [13,Lemma 4.3]. In particular, we write P loc as the composition of the maps P 0−6 and the result then follows from Lemma 5.4 andLemma 5.7 and Lemma 5.9, near q − and z − , along with similar results (these maps are basically the inverses (to leading order) of those in Lemma 5.7 and Lemma 5.9) near q + and z + . The fact that the remainder term can be bounded by a constant c n (α 0 ) follows from (81).
From this lemma, it follows that P 3 122x also satisfies the estimates (75) on x 2 ∈ [−c, c]. In fact, one can show (see [13,Theorem 1.3]) that for any l ∈ (0, 1), there exists constants, including c > 0, such that ( P 3 122x )(x 2 ) can be extended in such a way that (82) holds and such that ( We now write the regular map P glo in a similar way. In fact, we focus on P −1 glo . Let P −1 glox (x, µ) be the x-component of P 1 glo . Since it is regular it depends smoothly on x and on the unfolding parameter µ. By assumption 5, we have that P −1 glox (x in , 0) = x out . Consequently, we obtain the following expansion: There exists ν 0 ∈ (−1, 0) and ν 1 > 0 such that The fact that ν 0 ∈ (−1, 0) follows from the fact that γ 0 is repelling, see [13, Lemma 1.6]. Moreover, ν 1 > 0 follows by assumption 6.
To solve the fixed-point equation P x (x, µ) = x, we therefore solve P locx = P −1 glox . By (82) and (83) this gives (1), attains all values in [−1 + l, −l] with 0 < l < 1 + ν 0 < 1, we obtain a locally unique saddle-node of the fixed point by applying the implicit function theorem, see [13,Lemma 4.5]. The proof in the present case is identical. In this way, we have completed the proof of Theorem 5.1 (1).

5.3.
Dynamics on the blowup of Q. To prove Theorem 5.1 (2), we consider the regime ( , α) ∈ W 2 ( 0 , 1 , α 0 ) where 0 < α k+1 k 0 ≤ ≤ α k+1 k 1 . In this case, the dynamics within¯ = 0 becomes relevant, recall (68). We therefore decompose P loc in a different way, replacing Π 3 and Π 4 with Π 3 and Π 4 , respectively, see Lemma 5.2. In this way, since the mapping from Π 3 and Π 4 within 121 = 0 is completely "neutral" with no contraction, see (71), it follows that for all x ∈ I in and α 0 small enough, so that the dynamics is uniformly bounded in the (ȳ = 1) 1 -chart, then P loc is as close Figure 23. Illustration of the maps P loc and P −1 glo restricted to the slow manifold in the case ( , α) ∈ W 1 , see Theorem 5.1 (1). Here x in is the x-values of the orbit of Z + that grazes Σ on Π 0 . In the parameter regime ( , α) ∈ W 1 , the mapping P loc is then dominated by the attraction towards the attracting center manifold G 122 on one side (green) x x out and the dynamics of Z + (which itself is as close to x → −x as desired upon adjusting the domains) on the other side x x in . The transition inbetween (in purple), which extends over a O( 2k 2k+1 α 2k 2k+1 )-neighborhood of x in , is described by the Chini-equation, see (74), and it is concave cf. Lemma 5.4. On the other hand, since γ 0 is repelling, it follows that P glo is expanding. In particular, P glo moves with nonzero speed for µ ∼ 0 by assumption 6 and this therefore gives the saddle-node bifurcation of limit cycles as solutions of P −1 glo = P loc when the two graphs are tangent at a point.
as desired (upon adjusting the domains) to a reflection x → −x. Consequently, there can be no saddle-node bifurcations in this chart within this parameter regime.
In order to prove Theorem 5.1 (2) and describe the chaotic dynamics, we have to follow the set L − . Recall that this set is unbounded in the (ȳ = 1,ᾱ 1 = 1,r 12 = 1) 121chart, so we follow it across the first blowup cylinder and towards the blowup of the point Q.

5.4.
Completing the proof of Theorem 5.1 (2). Our strategy for completing the proof of Theorem 5.1 is as follows: Let µ ∼ 0 and consider 1 > 0 small enough, so that the system has a repelling limit cycle, that when written in the (ȳ = 1,ᾱ 1 = 1,r 12 = 1) 121 -chart intersects σ 12 = 1 transversally for each ( , α) ∈ W 2 (0, 1 , α 0 ) provided that α 0 > 0 is small enough. The fact that this is possible follows from the analysis above in the (ȳ = 1,ᾱ 1 = 1,r 12 = 1) 121 -chart, see the start of Section 5.3. Next, by decreasing α 0 > 0 if necessary there exists an 0 < 0 < 1 such that there is a canard trajectory for each ( , α) ∈ W 2 . In fact, the canard has an unstable foliation on the repelling side, which -when carried across N ,α near M -gives a twist-like return to the slow manifold S ,α . In particular, we obtain an associated foliation of points on S 0,1 ⊂ C 1 . This is the curve F . Since the limit cycle is repelling, we can track the canard backwards, and conclude that the limit cycle is the α-limit of the canard. The canard therefore transversally intersects the curve F at least n points for all ρ 213 > 0 small enough. The proof of the theorem then follows [14], which is inspired by [7] in a similar setting. In particular, we define a return map in the (ᾱ = 1,¯ = 1) 22 -chart using the scaling (86), with = ρ k+1 213 , defined on a section Π 22 transverse to M 22 and the canard. Since the expansion along M is greater than the contraction along C, recall Remark 3.6, we will study this mapping in backward time (so that M becomes attracting and C repelling). By flowing N ,α ∩ Π 22 backwards near the canard, we obtain -due to the transverse intersection of S ,α and N ,α along the canard -a stable foliation of the canard on the S ,α -side. For each transverse intersection i = 1, . . . , n of the canard with F , we then further obtain a small subset of this foliation which, upon extension by the backward flow, eventually returns to Π 22 in a "horizontal" curve H i that extends an O(1) distance in the direction tangent to N ,α ∩ Π 22 at the canard. At the same time, H i is exponentially close to N ,α ∩ Π 22 . This gives n disjoint horizontal curves H 1 , . . . , H n , whose preimages are n disjoint exponentially small intervals I 1 , . . . , I n on N ,α ∩ Π 22 . By the unstable foliation of N ,α , we obtain n "vertical strips" Figure 25. Blowup dynamics for ( , α) ∈ W 2 . On the left, we show the reduced, desingularized dynamics on C 1 in the (ȳ = 1) 1 -chart upon application of the two consecutive blowup transformations, see (56) and (61). In red we have indicated the repelling limit cycle γ 0 . It extends onto the cylinder¯ = 0, due to the blowup (61), in the singular limit , α → 0, with the limit understood within this parameter regime. In comparison with Fig. 21, we leave out the dynamics along¯ > 0 since this is not relevant for the regime ( , α) ∈ W 2 . At the same time, we also indicate the canard in cyan, see also Fig. 24, and the set F (black) which is the set of base points on C 1 , obtained by following the unstable foliation of the canard along M . On the right, we illustrate dynamics in the projection also used in Fig. 13, where the fast dynamics are also visible. Here we specifically indicate how the canard (cyan) extends across the two cylinders following C on top and M below. The section Π, transverse to M and the canard, is used in the proof of Theorem 5.1 (2).
In future work, it would be interesting to perform the same analysis for α < 0, but also, in the case of the grazing bifurcation, to explore the transition between the two regimes of Theorem 5.1. One possible explanation is that there is actually a curve in the ( , α)-plane along which saddle-node limit cycles "touch" or "grazes" the foliation of points, described by F in the singular limit and bounded by α 1 = 1 from above, due to the twist and return to S ,α away from the canard. An analysis of such a bifurcation scenario is interesting in its own right and in future work we aim to describe this in a simpler setting.