Geometric Blow-Up for Folded Limit Cycle Manifolds in Three Time-Scale Systems

Geometric singular perturbation theory provides a powerful mathematical framework for the analysis of ‘stationary’ multiple time-scale systems which possess a critical manifold, i.e. a smooth manifold of steady states for the limiting fast subsystem, particularly when combined with a method of desingularisation known as blow-up. The theory for ‘oscillatory’ multiple time-scale systems which possess a limit cycle manifold instead of (or in addition to) a critical manifold is less developed, particularly in the non-normally hyperbolic regime. We use the blow-up method to analyse the global oscillatory transition near a regular folded limit cycle manifold in a class of three time-scale ‘semi-oscillatory’ systems with two small parameters. The systems considered behave like oscillatory systems as the smallest perturbation parameter tends to zero, and stationary systems as both perturbation parameters tend to zero. The additional time-scale structure is crucial for the applicability of the blow-up method, which cannot be applied directly to the two time-scale oscillatory counterpart of the problem. Our methods allow us to describe the asymptotics and strong contractivity of all solutions which traverse a neighbourhood of the global singularity. Our main results cover a range of different cases with respect to the relative time-scale of the angular dynamics and the parameter drift. We demonstrate the applicability of our results for systems with periodic forcing in the slow equation, in particular for a class of Liénard equations. Finally, we consider a toy model used to study tipping phenomena in climate systems with periodic forcing in the fast equation, which violates the conditions of our main results, in order to demonstrate the applicability of classical (two time-scale) theory for problems of this kind.


Introduction
Many physical and applied systems featuring multiple time-scale dynamics can be mathematically modelled by singularly perturbed systems of ordinary differential equations in the standard form where x ∈ R m , y ∈ R n , (•) ′ = d/dt, 0 < ε ≪ 1 is a small perturbation parameter and the functions f, g : R m ×R n ×[0, ε 0 ] → R m ×R n are at least C 1 -smooth.If the limiting system obtained from (1) when ε → 0 possesses a critical manifold, i.e. if the set of equilibria S = {(x, y) : f (x, y, 0) = 0} forms an n-dimensional submanifold of the phase space R m ×R n , then system (1) can be analysed using a mathematical framework known as geometric singular perturbation theory (GSPT) [17,30,45,64].Typical GSPT analyses consist of two parts, depending on whether S is normally hyperbolic, i.e. depending on whether the linearisation with respect to the fast variables D x f | S when ε = 0 is hyperbolic.The local dynamics near normally hyperbolic submanifolds of S can be accurately described using Fenichel-Tikhonov theory [17,60] (see also [30,45,64,65]), which ensures that solutions are well approximated by regular perturbations of singular trajectories constructed as concatenations of trajectory segments from the distinct limiting problems which arise when ε → 0 in system (1) on the fast and slow time-scales t and τ = εt, respectively.However, Fenichel-Tikhonov theory breaks down in the nonnormally hyperbolic regime.The dynamics near non-normally hyperbolic points or submanifolds of S, which generically correspond to bifurcations in the layer problem (1)| ε=0 , can be studied using various techniques.A particularly powerful approach is the so-called blow-up method, which was pioneered for fast-slow systems in [15] and later in [42,43,44].In these works and many others (see [27] for a recent survey), the authors showed that the loss of hyperbolicity at a non-normally hyperbolic point or submanifold Q can often be 'resolved' after lifting the problem to a higher dimensional space in which Q is blown-up to a higher dimensional manifold Q.Following a suitable desingularization, which amounts to a singular transformation of time, a non-trivial vector field with improved hyperbolicity properties can be identified on Q.This allows one to study the dynamics in the non-normally hyperbolic regime using classical dynamical systems methods like regular perturbation and center manifold theory.We refer to [35,36,37,39,42,43,44,47,58,59] for seminal works in both applied and theoretical contexts, based on this combination of Fenichel-Tikhonov theory and the blow-up method.
The mathematical theory and methods described above are local in the sense that they apply specifically to fast-slow systems possessing an n-dimensional critical manifold S. A corresponding global theory in which the layer problem (1)| ε=0 possesses a limit cycle manifold in place of (or in addition to) a critical manifold is less developed.Following [24], we shall refer to the former class as stationary fast-slow systems, and the latter class as oscillatory fast-slow systems.A program for the development of a global GSPT which is general enough to encompass oscillatory fast-slow systems was initiated by J. Guckenheimer in 1996 [20], however a number of key analytical and methodological obstacles to its development remain.
One such obstacle concerns the development of Fenichel-Tikhonov theory for oscillatory fast-slow systems.A number of authors have made important contributions in this direction.Anosova showed that normally hyperbolic limit cycle manifolds in oscillatory fast-slow systems of the form (1) persist as O(ε)-close locally invariant manifolds for the perturbed system [1,2], similarly to the persistence of a normally hyperbolic critical manifold in Fenichel-Tikhonov theory.For oscillatory fast-slow systems with one slow variable, the persistence and contractivity properties of the center-stable/unstable manifolds have been described in detail in [28] using properties of the Poincaré map and a discrete GSPT framework developed therein.This study also yielded an asymptotic formula for the slow drift along the manifold, which agrees with the formula predicted by classical averaging theory [55].
Global theory for oscillatory fast-slow systems in the non-normally hyperbolic regime is less developed, despite the ubiquity of fast-slow extensions of global bifurcations in applications, for example in biochemical models exhibiting bursting [7,16,22,56].Notable exceptions include the study of a dynamic saddle-node of cycles bifurcation in a variant of the FitzHugh-Nagumo equations in [34], and the rigorous topological classification of so-called torus canards in [63], which occur near folded cycle bifurcations in the layer problem in oscillatory fast-slow systems [63].See e.g.[5,6,13,38,57] for more (predominantly numerical) important work on torus canards.The results in [63] were obtained using averaging theory, Floquet theory and (stationary) GSPT.Indirect contributions to the non-normally hyperbolic theory for oscillatory fast-slow systems have also been made via the study of non-normally hyperbolic singularities in fast-slow maps, since these can be used to infer dynamical properties of corresponding limit cycle bifurcations (or fast-slow extensions thereof) in one greater dimension.The slow passage through a flip/period-doubling bifurcation (and even through an entire period-doubling cascade) was treated in [3,4], and further results on the slow passage through discrete transcritical, pitchfork and Neimark-Sacker/torus type bifurcations have been derived using non-standard analysis; we refer to the review paper [18] and the references therein.
In general, the development of mathematical methods for handling non-hyperbolic dynamics in oscillatory fast-slow systems is complicated by the fact that local techniques like the blow-up method rely upon near-equilibrium properties possessed by stationary but not oscillatory fast-slow systems.The desingularization step in blow-up analyses, for example, relies upon a formal division of the blown-up vector field by zero.This step is crucial for obtaining a desingularized vector field with improved hyperbolicity properties, but it is only well-defined if the original (blown-up) vector field is in equilibrium wherever it is formally divided by zero.For oscillatory fast-slow systems of the form (1), a 'typical' non-normally hyperbolic cycle has non-equilibrium dynamics, so it cannot be desingularized and the blow-up method does not apply.
The aim of this work is to show that stationary methods (in particular the blow-up method) which may not apply for oscillatory fast-slow systems, may be applicable in the study of oscillatory multiple time-scale systems with at least three distinct time-scales.Specifically, we consider systems of the form where (r, θ, y) ∈ R ≥0 × R/Z × R are cylindrical coordinates, 0 < ε 1 , ε 2 ≪ 1 are singular perturbation parameters, and f, g, h : R ≥0 × R/Z × R × [0, ε 1,0 ] × [0, ε 2,0 ] → R are sufficiently smooth for our purposes.Depending upon the relative magnitude of ε 1 and ε 2 , system (2) has either two or three distinct timescales.Although we present results on a range of different cases, we are primarily interested in the case 0 which defines a class of three time-scale semi-oscillatory systems that are in a certain sense 'in between' the classes of stationary and oscillatory fast-slow systems described above.Heuristically, this is because under suitable assumptions (to be outlined in detail in Section 2), system ( 2) is an oscillatory fast-slow system with respect to the limit ε 1 > 0, ε 2 → 0, and a stationary fast-slow system with respect to the double singular limit (ε 1 , ε 2 ) → (0, 0) [46].It is worthy to note that multiple time-scale systems with more than two time-scales appear frequently in applications.The long term dynamics of a forced van der Pol oscillator with three time-scales was studied as early as 1947 in [11].A theoretical basis for normally hyperbolic theory for stationary multiple time-scale systems with three or more time-scales has appeared more recently in e.g.[9,10,40,50].Progress has also been made in the non-normally hyperbolic setting, particularly via the study of three time-scale applications and 'prototypical systems' inspired by applications; we refer to [12,14,26,31,32,33,41,49,52].We present a detailed analysis of the 'jump-type' transition near a non-normally hyperbolic cycle of regular fold type using geometric blow-up.More precisely, we assume that the limiting system (2) ε1>0,ε2=0 undergoes a type of folded cycle bifurcation under variation in y which is common in applications with periodic forcing in the slow equation.This global bifurcation is closely related to (and can in many ways be seen as an oscillatory analogue to) the regular fold or jump point in stationary fastslow systems, which has been studied in R 2 and R 3 in particular using blow-up techniques in [42,59]; see also [51] for a detailed treatment of the planar case using classical asymptotic methods.After deriving a prototypical system by imposing a number of defining assumptions on system (2), we show that the blow-up method can be applied, even though a rigorous reduction to the stationary setting is not possible due to angular coupling.The formal division by zero which is necessary to obtain a desingularized vector field with improved hyperbolicity properties is possible if the time-scale associated to the rotation is sufficiently slow relative to the fast radial dynamics, i.e. if ε 1 is sufficiently small.
The blow-up analysis allows for the detailed characterisation of the transition map induced by the flow, including the asymptotic and contractivity properties of the transition undergone by solutions traversing the neighbourhood of the global singularity.If ε 2 /ε 1 ∼ 1 or ε 2 /ε 1 ≫ 1 (so that (3) is not satisfied), then system (2) is a stationary fast-slow system.In this case, the local dynamics can (for the most part) be described in detail using the results established for two time-scale systems in [42,51,59], after Taylor expansion about a given jump point.However, these results do not apply directly to the scaling regimes we consider which satisfy (3), i.e. they do not apply in the semi-oscillatory case.One important reason for this is that the singularity is 'global' in the angular coordinate θ.Consequently, it does not suffice to blow-up at a single point on the non-hyperbolic cycle.Rather, it is necessary to blowup the entire non-hyperbolic cycle to a 'torus of spheres' S 1 × S 2 .A similar approach is adopted in the geometric analysis of the periodically forced van der Pol equation in [8], however in our case, the leading order equations derived on the blown-up sphere may depend upon the angular variable θ, which remains non-local.As a consequence, the local dynamics cannot be analysed with a straightforward adaptation of the arguments used to study the dynamics near a regular fold point/curve in [42,59].Rather, new arguments are needed.We derive results for two different scaling regimes defined by (ε 1 , ε 2 ) = (ε α , ε 3 ), where α ∈ {1, 2} and 0 < ε ≪ 1.In each case, the size of the leading order term in the asymptotics for the parameter drift in y is shown to agree with the known results for the stationary regular fold point [42,59].In contrast to the classical fold, however, the leading order coefficient is shown to depend on θ, and we provide an explicit formula for this dependence in the case α = 2.We also provide asymptotics for the angular coordinate θ as a function of the initial conditions and small parameters, and an asymptotic estimate for the number of rotations about the y-axis over the course of the transition.The results obtained are shown to depend on the relative magnitude of ε 1 and ε 2 (i.e. on α), with the main qualitative difference pertaining to the asymptotic estimates for θ and the corresponding number of rotations.
Finally, we apply our results in order to derive detailed asymptotic information near folded limit cycle manifolds of periodically forced Liénard equations, and we consider a simple model proposed in [66] as a toy model for the study of tipping phenomena in climate systems.Our main results do not apply directly to the latter problem, due to periodic forcing in the fast equation.Rather, we aim to demonstrate with a partial but illustrative geometric analysis that problems of this kind can be treated using classical approaches based on established results for two time-scale systems.
The manuscript is structured as follows.In Section 2 we introduce defining assumptions and present the prototypical system for which our main results are stated.The singular dynamics and geometry, which differ in different scalings, are presented in Section 2.1.The main results are presented and described in Section 3, and the blow-up analysis and proof of the main results are presented in Section 4. The applications are treated in Section 5. Specifically, in Section 5.1 we apply our main results to periodically forced Liénard equations, and in Section 5.2 we consider the toy model for the study of tipping phenomena from [66] which cannot be treated directly with the results from Section 3. We conclude with a summary and outlook in Section 6.

Assumptions and Setting
We consider C k -smooth multiple time-scale systems in the general form (2), restated here for convenience: where k ∈ N will be assumed to be 'sufficiently large' throughout, (•) ′ = d/dt, the variables are given in cylindrical coordinates (r, θ, y) ∈ R ≥0 × R/Z × R, and ε 1 , ε 2 are singular perturbation parameters satisfying 0 < ε 1 , ε 2 ≪ 1.Note also that smoothness implies that f, g, h are 1-periodic in θ.System (4) evolves on either two or three time-scales, depending on whether the ratio ε 1 /ε 2 is asymptotically large, constant or small.The setup and defining assumptions presented below are primarily motivated by the case 0 < ε 2 ≪ ε 1 ≪ 1, for which system (4) in a certain sense 'intermediate' between stationary and oscillatory fast-slow systems.We shall refer to this as the semi-oscillatory case.Since we also consider other possibilities, however, we leave the exact relation between ε 1 and ε 2 unspecified for now.It suffices to observe that the forward evolution of a generic initial condition is characterised by radial motion on the fast time-scale t, angular motion on a time-scale τ ε1 = ε 1 t, and vertical 'parameter drift' on a time-scale τ ε2 = ε 2 t.
In the following we impose a number of defining conditions in terms of the limiting oscillatory fast-slow system obtained in the singular limit ε 1 > 0, ε 2 → 0, i.e. on r ′ = f (r, θ, y, ε 1 , 0), We remark that the singular limit ε 1 > 0, ε 2 → 0 is only 'natural' in the semi-oscillatory case, i.e. if ε 1 /ε 2 ≫ 1 so that the rotation is fast relative to the parameter drift.
Remark 2.1.The assumption that the limit cycle S c 0 is circular and in particular the zero condition on f in (6), is natural in applications with small-amplitude external periodic forcing (amplitudes of O(ε 2 ) or smaller).However, it rules out applications with 'large' periodic forcing in the fast equation for r.One reason for imposing such a restriction is that problems of the latter kind can often be treated using classical theory for two time-scale systems.This is demonstrated for a particular application in Section 5.2.
We shall be interested in the dynamics in a neighbourhood of S c 0 .This motivates the introduction of the signed radius variable r := r − v, in which case S c 0 = {(r, θ, y) : r = 0, θ ∈ R/Z, y = 0}.We assume without loss of generality that g(v, θ, 0, ε 1 , 0) > 0, divide the right-hand side of system (4) by g(r + v, θ, y, ε 1 , ε 2 ), i.e. we apply a time-dependent transformation satisfying d t = g(r, θ, y, ε 1 , ε 2 )dt (which is positive in a sufficiently small toroidal or tubular neighbourhood V of S c 0 ), and rewrite the system in (r, θ, y)-coordinates in order to obtain r′ = F (r, θ, y, ε 1 , ε 2 ), where by a slight abuse of notation the dash now denotes differentiation with respect to t, and In the derivation of system (7) and in many of the proofs below we make use of transformations of time which are formulated in terms of differentials, e.g.d t = g(r, θ, y, ε 1 , ε 2 )dt.Strictly speaking, such an 'transformation' only defines t as a unique function of t up to an additive constant.
Since we are interested in the behaviour of solutions to autonomous ODEs, which are invariant under time translation, this additive constant can be set to zero without loss of generality.
We are interested in three time-scale systems (7) which feature a regular fast-slow fold of limit cycles with respect to the partial singular limit ε 1 > 0, ε 2 → 0. Necessary and sufficient conditions for this to occur can be given in terms of the Poincaré map induced on the transversal section ∆ obtained by intersecting the toroidal/tubular neighbourhood V ⊃ S c 0 with the half-plane defined by θ = 0 (decreasing the size of V if necessary).Proposition 2.3.System (7) with ε 1 ∈ (0, ε 1,0 ) fixed and 0 < ε 2 ≪ 1 induces a Poincaré map P : ∆ → ∆ given by we have The result follows after Taylor expanding about ε 2 = 0.
The defining conditions for S c 0 to be a regular fast-slow fold of cycles are as follows: together with for all ε 1 ∈ (0, ε 1,0 ).The conditions in ( 9)-( 10) are in 1-1 correspondence with the defining conditions for a fold bifurcation in the 1D 'layer map' r → r + P r (r, y, ε 1 , 0) (with y as a bifurcation parameter), see e.g.[48,Ch. 4.3], except for the integral condition, which can be viewed as the analogue of the slow regularity condition on the fast-slow regular fold point in planar continuous time systems in [42,45].
Remark 2.4.The defining conditions for a regular fast-slow fold of cycles in (9)- (10) do not depend on the specific form of the Poincaré map in (8).Nevertheless, we have chosen to state Proposition 2.3 prior to the conditions in (9)- (10) in order to clarify the interpretation of the slow regularity condition, i.e. the integral expression in (10), which is not common in the literature.This condition can be viewed as a condition on the 'reduced map'; we refer to [28,29] for details.
In the following, we shall actually assume stronger conditions that are sufficient but not necessary for a regular fast-slow fold of cycles, instead of those in ( 9)- (10).Similarly to Assumption 1, these conditions are expected to be satisfied in applications with small-amplitude external periodic forcing; see again Remark 2.1.

Assumption 2. (Sufficient conditions for S c
0 to be a regular folded cycle) The following sufficient (but not necessary) conditions for a fast-slow fold of cycles are satisfied by system (7): together with for all θ ∈ R/Z and ε 1 ∈ (0, ε 1,0 ).
It is straightforward to verify that the conditions in (11)-( 12) are sufficient to ensure that the Poincaré map (8) satisfies the fold conditions ( 9)- (10).In particular, (11)- (12) are directly analogous to the defining conditions for the (stationary) regular fold point in [42], except that we require them to hold for θ-dependent functions.
Assumptions 1-2 and the implicit function theorem imply that system (5) has a two-dimensional manifold of regular limit cycles where V ⊆ R ≥0 × R/Z × R is the toroidal/tubular neighbourhood about S c 0 introduced above.In other words, system (7) (and by a simple computation also system (4)) is an oscillatory fast-slow system with respect to the (partial) singular limit ε 1 > 0, ε 2 → 0.
The combination of signs taken by the various non-zero terms in (12) determines the orientation of the bifurcation.In the following we assume without loss of generality that which are consistent with a 'jump-type' orientation in forward time; see Figures 1 and 2. Based on Assumptions 1-2 and these sign conventions, it suffices to work with the simplified system provided in the following result.
System ( 14) is related to the local normal form for the regular fold point in [42] by a simple variable rescaling if ε 1 = 0, i.e. if θ is fixed.For ε 1 > 0, however, the dependence on the angular variable appears via the functions a(θ), b(θ), c(θ) and the higher order terms R r (r, θ, y, ε 1 , ε 2 ) and R y (r, θ, y, ε 1 , ε 2 ).In what follows, we drop the tilde notation on r and work with system (14) for the remainder of this manuscript.
Remark 2.6.In Proposition 2.5 we assert that the functions a(θ), b(θ) and c(θ) are 'smooth'.A more precise statement would be to say that they are C k -smooth, since the system obtained after Taylor expansion is precisely as smooth as the original system (7).Since we will not be interested in smoothness per se, we shall adopt a similar terminology throughout for simplicity, i.e. by 'smooth' we shall mean sufficiently smooth for the validity of our methods (e.g.Taylor expansions).
Remark 2.7.The equation for r in system (14) can be further simplified after setting r = b(θ)/a(θ) r.This leads to where a ′ := ∂a/∂θ, b ′ := ∂b/∂θ and ι(θ) := a(θ)b(θ).The 'price' of this simplification, however, is that the O(ε 1 r) term also appears in the leading order terms in the blow-up analysis in later sections.
For this reason, we continue to work with the formulation in (14).
Remark 2.8.If ε 1 ∼ ε 2 or ε 1 ≪ ε 2 , then the θ-variable is 'slow enough' to validate the Taylor expansion of system (14) about a fixed point (0, θ * , 0) ∈ S c 0 , i.e. in this case, one can also Taylor expand in the angular coordinate θ.This allows for a subsequent transformation into the simpler local normal form near a fold curve in [59,Lem. 3], thereby showing that the dynamics in these cases are governed by the well-known result for two time-scale systems in [59,Thm. 1].In the semi-oscillatory case of interest with ε 2 ≪ ε 1 , however, θ is fast relative to y, and varies over the entire domain R/Z as solutions approach S c 0 .As a consequence, one cannot Taylor expand the θ coordinate, and transformation to the local normal form in [59] is not possible.Remark 2.9.In order to sketch geometric objects like S 0 in the upcoming figures, we choose the positive, 1-periodic and smooth functions a(θ) = 2 + sin(4πθ) and b(θ) = 5 + cos(2πθ − 1).For constant functions a and b, the figures including the θ-coordinate would be rotationally symmetric.
The stability of S 0 with respect to the fast radial dynamics is determined by the unique non-trivial (i.e.not identically zero) eigenvalue of the linearisation, namely It follows that the critical manifold has the structure S 0 = S a 0 ∪ S c 0 ∪ S r 0 , where are normally hyperbolic and attracting/repelling, respectively (assuming r 0 > 0 is sufficiently small).The circular set S c 0 = {(0, θ, 0) ∈ S 0 }, which corresponds to the regular folded cycle in Assumption 2, is non-normally hyperbolic.The situation is sketched in Figures 1 and 2.
The dynamics and geometry for the layer problem (15) do not depend upon the relative magnitude of ε 1 and ε 2 .The reduced dynamics on S 0 , however, are expected to differ significantly depending on the size of ε 1 /ε 2 .As noted already in Section 1, we consider three distinct possibilities: (C1) Angular dynamics are fast relative to the parameter drift, i.e.
(C2) Angular dynamics occur on the same time-scale as the parameter drift, i.e. there is a constant σ > 0 such that Figure 1: Singular geometry and dynamics in case (C1).Fast and slow dynamics are depicted (here and throughout) by double and single arrows, respectively.The attracting and repelling normally hyperbolic submanifolds of the critical manifold S 0 , denoted by S a 0 and S r 0 , are shown in shaded turquoise and gray, respectively, and sketched for the particular choice of a(θ) and b(θ) defined in Remark 2.9.The non-normally hyperbolic folded cycle S c 0 is shown in blue.The reduced flow in case (C1) is periodic, i.e. y is a parameter and S 0 is foliated by limit cycles of period τ ε1 = 1.
. Figure 2: Singular geometry and dynamics in case (C3) sketched for the particular choice of a(θ) and b(θ) in Remark 2.9.The dynamics is distinguished from cases (C1) and (C2) by the reduced flow on S 0 .In this case, θ is a parameter and singular orbits (concatenations of solution segments of layer and reduced problem) are contained within constant angle planes {θ = const.}.An example of such an orbit is sketched in red.
(C3) Angular dynamics are slow relative to the parameter drift, i.e.
A different reduced problem is obtained on S 0 in each case.We briefly consider each case in turn.
The Reduced Problem in Case (C1) In this case we rewrite system (14) on the slow time-scale τ ε1 = ε 1 t.This leads to where the dot denotes differentiation with respect to the slow time τ ε1 .Since ε 1 /ε 2 ≫ 1, we first take ε 2 → 0, and then ε 1 → 0 (in that order).This leads to the reduced problem In this case, S 0 is foliated by limit cycles of period τ ε1 = 1, i.e. t = 1/ε 1 .An expression for the vector field on S 0 , expressed in the (r, θ)-coordinate chart, can be obtained by differentiating the constraint y = φ 0 (r, θ) with respect to τ ε1 and rearranging terms.We obtain where a ′ := ∂a/∂θ and b ′ := ∂b/∂θ.This case is sketched in Figure 1.

The Reduced Problem in Case (C2)
To obtain the reduced problem in case (C2) we may write system ( 14) on either time-scale τ ε1 or τ ε2 , which are related via τ ε1 = στ ε2 .Writing the system on the τ ε2 time-scale leads to where this time, the dot denotes differentiation with respect to the slow time τ ε2 .Taking the double singular limit and using the fact that ε 1 /ε 2 = σ leads to the reduced problem In contrast to the limiting equations (17) in case (C1), none of the variables become parameters in system (19).This is natural because in case (C2), system ( 14) only has two (as opposed to three) time-scales.In the (r, θ)-coordinate chart, the dynamics are determined by In particular, solutions reach r = 0 (and therefore S c 0 ) in finite time.
The Reduced Problem in Case (C3) In this case we rewrite system ( 14) on the slow time-scale τ ε2 = ε 2 t, thereby obtaining system (18).Since ε 1 /ε 2 ≪ 1, we first take ε 1 → 0, and then ε 2 → 0 (in that order).This leads to the reduced problem This time, the angular variable θ is the slow variable to be considered as a parameter in system (19).In particular, the dynamics on S 0 can be represented by the 1-parameter family of ODEs where a = a(θ), b = b(θ) and c = c(θ) are constants paramaterised by θ ∈ R/Z.Similarly to case (C2), solutions reach S c 0 in finite time.Moreover, singular orbits obtained as concatenations of layer and reduced orbit segments are contained within constant θ planes.As a result, the singular geometry and dynamics in case (C3) is equivalent to the singular geometry and dynamics of the normal form for the planar regular fold point in [42].This case is sketched in Figure 2.
For 0 < ε 1 , ε 2 ≪ 1, Fenichel-Tikhonov theory implies that compact submanifolds of the normally hyperbolic critical manifolds S a 0 and S r 0 persist as O(l(ε 1 , ε 2 ))-close locally invariant slow manifolds S a l(ε1,ε2) and S r l(ε1,ε2) , respectively [17,30,45,64,65], where we write l(ε 1 , ε 2 ) := max{ε 1 , ε 2 } in order to keep the discussion general, i.e. so that we need not distinguish between cases (C1), (C2) and (C3).Our goal is to describe the extension of the attracting slow manifold S a l(ε1,ε2) through a neighbourhood of the non-hyperbolic cycle S c 0 corresponding to the regular folded cycle in system ( 14).Remark 2.10.In [10], the authors extend GSPT for a class of multiple time-scale systems with n ≥ 3 time-scales which feature 'nested critical manifolds'.A requirement for the application of this theory to system (14) is that the reduced problem on S 0 has a one-dimensional critical manifold.This condition is not satisfied in any case (C1), (C2) or (C3) because none of the reduced problems (17), (19) and (21) have equilibria in the neighbourhood of interest (i.e.close to r = y = 0).
Note that we have introduced an additional case (C1) * .This case is dynamically distinct from the others, but it is not distinguished in Section 2 because the geometry and dynamics in the double singular limit are the same as for case (C1).
Remark 3.1.The choice to write ε 3 instead of ε in the equation for y is made a-posteriori in order to avoid fractional exponents in the proofs.Comparisons with known results for the stationary regular fold point in [42,59] are possible via the simple relation ε KSW = ε 3 , where we denote by ε KSW the small parameter in [42] and/or [59].Similar observations motivated the use of a cubic exponent for the small parameter in other works involving folded singularities, e.g. in [53,54].
Our aim is to describe the forward evolution of initial conditions in an annular entry section where R is a small positive constant and β − < β + < 0 are two negative constants chosen such that We track solutions of system (22) up to their intersection with the cylindrical exit section for a small positive constant y 0 > 0. The critical manifold S 0 , the Fenichel slow manifolds (denoted now by S a ε and S r ε ) and the entry/exit sections are visualized in the (r, y)-plane in Figure 3 and in the three-dimensional space in Figure 4.
More precisely, ∂ ∂r h (α) y (h Theorem 3.2 characterises the asymptotic behaviour of solutions and the extension of the attracting Fenichel slow manifold S a ε through a neighbourhood of the regular folded cycle in system (22).The geometry and dynamics for all four cases (C1) * , (C1), (C2) and (C3) are sketched in Figures 5, 6, 7 and 8, respectively.The proof is broken into two parts, depending on whether or not α ∈ {1, 2}.A detailed proof based on the blow-up method will be given for the cases α ∈ {1, 2} in Section 4. Recall that these are the primary cases of interest, since they correspond to semi-oscillatory dynamics.If α ≥ 3, minor adaptations of the proof for case α = 2 can be applied.Moreover, in this case, one can show directly that S c 0 is a closed regular fold curve.From here, Assertions (a)-(c) can be proven directly using established results (or a straightforward adaptation thereof) on the passage past a fold curve in 2-fast 1-slow systems, specifically [59,Thm. 1].We shall therefore omit the details of these cases.Note that for α ≥ 4 in particular, h (α) θ (r, θ, ε) ∼ θ as ε → 0, which implies that the leading order coefficient in the expansion for h (α) y (θ, ε) is also constant in θ.This is expected, since the dynamics for α ≥ 4 should resemble the dynamics of the stationary fold point considered in [42] up to a slight perturbation.
For each α ∈ N + , the size of the leading order term in the asymptotics for the parameter drift in y is also the same as for the stationary regular fold point, at least up to O(ε 3 ln ε), since which agrees with the asymptotic estimates in [42,51,59] (recall from Remark 3.1 that ε KSW = ε 3 , where ε KSW denotes the small parameter in [42,59]).The strong contraction property in Assertion (c), which The critical manifold and its submanifolds are sketched in colours consistent with earlier figures for the particular choice of a(θ) and b(θ) defined in Remark 2.9.The entry, exit sections ∆ in , ∆ out (magenta, orange) and the (extended) Fenichel slow manifolds S a ε , S r ε are also shown as shaded regions in green and light gray, respectively.For S a ε and S a 0 additionally a sample trajectory for a fixed θ initial condition is shown.
. Figure 4: The extension of the attracting slow manifold S a ε (again in shaded green) as described by Theorem 3.2, in the full (r, θ, y)-space.The entry, exit sections ∆ in , ∆ out as well as the critical manifold and its submanifolds are also shown, with the same colouring as in Figure 3 for the particular choice of a(θ) and b(θ) in Remark 2.9.The intersection π (α) (S a ε ∩∆ in ) ⊂ ∆ out is topologically equivalent to a circle (shown in dark green), and O(ε 2 )-close to the plane {y = 0} in the Hausdorff distance.The specific behavior of solutions, which depends on α, is not shown (see however Figures 5,6,7 and 8).
does not depend on α, is also the same.These two facts explain the similarity between the dynamics observed in the (r, y)-plane and the dynamics near a stationary regular fold point; c.f. Figure 3 and [42, Figure 2.1].In contrast to the stationary fold, however, the coefficients of h (α) y (θ, ε) depend on θ, and (26) gives the precise form for the leading order coefficient if α > 1, i.e. if α ∈ N + \ {1}.Moreover, Theorem (3.2) provides detailed asymptotic information about the angular coordinate θ via (25).If α > 2 then estimates are sharp as ε → 0.
Finally, the dynamics in different cases, i.e. for differing values of α, are distinguished via the angular dynamics and in particular, the number of complete rotations about the y-axis during the transition from ∆ in to ∆ out .We can estimate this number by introducing the function We obtain the following corollary as an immediate consequence of Theorem 3.2.
complete rotations about the y-axis during its passage to ∆ out .
rot (r, θ, ε) = ⌊O(ε −2 )⌋ rotations about the y-axis before leaving the neighbourhood close to S a ε ∩ ∆ out at an angle approximated by the expression for h rotations about the y-axis before leaving the neighbourhood close to S a ε ∩ ∆ out at an angle approximated by the expression for h    (22).The solution sketched in red makes N (α) rot (r, θ, ε) = 0 rotations about the y-axis before leaving the neighbourhood close to S a ε ∩ ∆ out at an angle approximated by the expression for h Proof.This follows immediately from the asymptotic estimates for h(α) (r, θ, ε) in Theorem 3.2 and the definition of N (α) rot .

Proof of Theorem 3.2
Our aim is to investigate the three time-scale system (22) with α ∈ {1, 2} using the blow-up method developed for fast-slow systems in [15,42,43,44]; we refer again to [27] for a recent survey.Many aspects of the proof rely in particular on arguments used in the blow-up analysis of the (stationary) regular fold point in [42].However, many of the calculations are complicated by the fact that the angular variable θ cannot be treated locally.In particular, for α ∈ {1, 2}, we cannot simply Taylor expand the equations about a fixed value of θ.Consequently, a local transformation into the local normal form in [59] is not possible.A similar feature arises when blowing up the fold cycle in the periodically forced van der Pol equation in the 'intermediate frequency regime' [8], except that in our case, there is no decoupling of the angular dynamics in the leading order.We adopt the (now well established) notational conventions introduced in [42,43,44].
The blow-up transformation is defined in Section 4.1, as are the three local coordinate charts that we use for calculations.The geometry and dynamics in all three coordinate charts are considered in turn in Sections 4.2, 4.3 and 4.4.Theorem 3.2 is proved in Section 4.5 using the information obtained in local coordinate charts.

Blow-up and Local Coordinate Charts
As is standard in blow-up approaches, we consider the extended system obtained from (22) after appending the trivial equation ε ′ = 0, i.e.
In order to describe the map π (α) : ∆ in → ∆ out , we introduce extended entry and exit sections and respectively.The entry, exit sections ∆ in , ∆ out defined in Section 3 (recall equations ( 23) and ( 24)) can be viewed as constant ε sections of ∆ in ε , ∆ out ε , respectively.Based on this simple relationship, we study the map π (α) : ∆ in → ∆ out described in Theorem 3.2 via the extended map π (27).
We now define the relevant blow-up transformation.Let I = [0, ρ 0 ], where ρ 0 > 0 is fixed small enough for the validity of local computations, let and define the (weighted) blow-up transformation via where (r, ȳ, ε) ∈ S 2 .The blow-up map φ is a diffeomorphism for ρ ∈ (0, ρ 0 ], but not for ρ = 0.In particular, the preimage of the non-hyperbolic cycle S c 0 under φ is a 'torus of spheres' S 2 × R/Z × {0} ∼ = S 2 × S 1 .For calculational purposes, we introduce local coordinate charts in order to describe the dynamics on Following [42], we introduce affine projective coordinates via Here we permit a small abuse of notation by introducing a new variable ε 1 , which should not be confused with the small parameter with the same notation in Sections 1-2.This leads to the following coordinates: In the analysis, it will be necessary to change coordinates between different charts.The following lemma provides the relevant change of coordinates formulae.
Lemma 4.1.The change of coordinate maps κ ij from K i to K j are diffeomorphisms given by for y 2 > 0, Proof.This follows from the local coordinate expressions in (31).
Remark 4.2.The blow-up transformation (30) has the same form as the blow-up map used to study dynamics near a regular fold curve in R 3 in [59], except that the domain of the uncoupled variable θ is R/Z instead of R. Note also that since θ is unaffected by (30), it follows that φ can be defined more succinctly in terms of its action on the remaining variables, i.e. via the map which is precisely the blow-up map used to study the regular fold point in [42] (recall that ε KSW = ε 3 by the discussion following the statement of Theorem 3.2).
Remark 4.3.By construction, the blown-up vector field induced by the pushforward of the vector field induced by system (27) under φ is invariant in the hyperplanes {ρ = 0} and {ε = 0}.The former corresponds to blown-up preimage of the non-hyperbolic cycle S c 0 , i.e. the torus of spheres S 2 × S 1 .The latter contains the blown-up preimage of the critical manifold S 0 .Since r2 + ȳ2 = 1 in {ε = 0}, the preimage of Thus, the part of the blown-up singular cycle S c 0 within {ε = 0} is a torus.
Remark 4.4.Since ε = const.in system (27) we have constants of the motion defined by We turn now to the dynamics in charts K i , i = 1, 2, 3.

Dynamics in the Entry Chart K 1
In chart K 1 we analyse solutions which track the blown-up preimage of the attracting Fenichel slow manifold S a ε as they enter a neighbourhood of the non-hyperbolic circle S c 0 .
Lemma 4.5.Following the positive transformation of time ρ 1 dt = dt 1 , the desingularized equations in chart K 1 are given by where by a slight abuse of notation we now write (•) ′ = d/dt 1 .
Proof.This follows after direct differentiation of the local coordinate expressions in (31) and subsequent application of the desingularization ρ 1 dt = dt 1 .
The analysis in chart K 1 is restricted to the set where R is the constant which defines the entry set ∆ in in ( 23) and E = ε 0 /R > 0 due to the relationship ε = ρ 1 ε 1 (recall Remark 4.4).The set D 1 is sketched with other geometric and dynamical objects in chart K 1 in Figure 9. System (33) is well-defined within {ρ 1 = 0} (recall that α ∈ {1, 2}), which corresponds to the part of the blow-up surface in K 1 .The subspace {ε 1 = 0} is also invariant, and contains two two-dimensional critical manifolds The manifolds S a 1 and S r 1 correspond to the blown-up preimages of the critical manifolds S a 0 and S r 0 in chart K 1 , respectively.Both S a 1 and S r 1 are topologically equivalent to cylindrical segments, and they are normally hyperbolic and attracting resp.repelling up to and including the sets , θ 1 , 0, 0 : , θ 1 , 0, 0 : which intersect with the blow-up surface; see Figures 9 and 10.The linearisation of system (33) along P a is where we write a ′ = ∂a/∂θ 1 and b ′ = ∂b/∂θ 1 , and we write In both cases, the set P a is partially hyperbolic with a single non-trivial eigenvalue −2 a(θ 1 )b(θ 1 ) < 0. The remaining three eigenvalues along P a are identically zero.Thus in the blown-up space (i.e. for system (33)), we have regained partial hyperbolicity.This allows us to extend the attracting center manifold with base along S a 1 up onto the blow-up surface using center manifold theory.Lemma 4.6.Consider system (33) on D 1 with E, R > 0 sufficiently small.There exists a threedimensional center-stable manifold M a 1 such that , where S a 1 ⊂ {ε 1 = 0} is the two-dimensional attracting critical manifold in (34) and N a 1 ⊂ {ρ 1 = 0} is a unique two-dimensional center-stable manifold for the restricted system (33)| ρ1=0 emanating from P a .The manifold M a 1 is given locally as a graph , where and In both cases, there exists a constant ϱ ∈ (0, ϑ), where ϑ := min θ1∈[0,1) 2 a(θ 1 )b(θ 1 ) > 0, such that initial conditions in Σ in 1 are attracted to M a 1 along one-dimensional stable fibers faster than e −ϱt1 .
Figure 9: Geometry and dynamics within D 1 , shown in (r 1 , ρ 1 , ε 1 )-space.Projections of the twodimensional attracting and repelling critical manifolds S a 1 , S r 1 ⊂ {ε 1 = 0} (shown in blue and shaded blue) extend up to their intersection with the blow-up surface at P a , P r ⊂ {ρ 1 = ε 1 = 0}, respectively.Projections of the two-dimensional center manifolds N a 1 , N r 1 ⊂ {ρ 1 = 0} emanating from P a , P r , as described in Lemma 4.6, are sketched in red and shaded red for the case α = 2, but they look similar in the case α = 1.Entry and exit sections Σ in 1 and Σ out 1 are shown in shaded magenta and yellow, respectively.The projected three-dimensional center manifold M a 1 and its image under Π (α) 1 is described by Proposition 4.9.This is shown along with the image Π Proof.The existence of a three-dimensional center manifold M a 1 follows from the linearisation (35) and center manifold theory [48].The fact that M a 1 contains two-dimensional manifolds N a 1 and S a 1 with the properties described follows after direct restriction to the invariant hyperplanes {ρ 1 = 0} and {ε 1 = 0}, respectively.
The two-dimensional center manifold N a 1 is sketched within {ρ 1 = 0} in Figure 11.Remark 4.7.Similar to Lemma 4.6, there exists a three-dimensional center-unstable manifold M r 1 at P r that contains the repelling critical manifold S r 1 ⊂ {ε 1 = 0} and a repelling center manifold N r 1 ⊂ {ρ 1 = 0}.These objects are shown in Figures 10 and 11.In D 1 , the manifold M r 1 is given as a graph r 1 = h(α) r1 (θ 1 , ρ 1 , ε 1 ), where In contrast to N a 1 , the manifold N r 1 is not unique.We omit the explicit treatment of the dynamics near M r 1 since it is not relevant to the proof of Theorem 3.2.Lemma 4.6 can be used to describe the r 1 -component of the map Π induced by the flow of initial conditions in e. the representation of the entry section ∆ in ε defined in (28) in K 1 coordinates, up to the exit section Σ out 1 The constants β − < β + < 0 are chosen in such a way that the r 1 -coordinate of the center manifold M a 1 in D 1 is always contained in the interval (β − /R, β + /R), cf.Definition (23) of the entry section ∆ in .
The remaining θ 1 , ρ 1 and ε 1 components of the map Π (α) 1 and the transition time taken for solutions with initial conditions in Σ in 1 to reach Σ out 1 can be estimated directly.To obtain better estimates, we rewrite the positive, 1-periodic smooth function c(θ 1 ) as where c 0 := 1 0 c(θ 1 (t 1 )) dt 1 is the mean value of c(θ 1 ) over one period and c rem (θ 1 ) := c(θ 1 ) − c 0 is the smooth, 1-periodic remainder with mean zero.The estimates are provided by the following result.Lemma 4.8.Consider an initial condition (r (33).Then , , where ϕ(t 1 ) := is the transition time taken for the solution to reach Σ out 1 .
The expression for ρ 1 (t 1 ) can be obtained directly from the expression for ε 1 (t 1 ) using the fact that ρ 1 (t)ε 1 (t) = ε = Rε * 1 is a constant of the motion; recall Remark 4.4.One can show that |ψ(t 1 )| = O(R) by appealing to the fact that χ(r 1 (s 1 ), θ 1 (s 1 ), ρ 1 (s 1 ), ε 1 (s 1 )) is bounded uniformly for all It remains to estimate θ 1 (t 1 ) and the transition time T 1 .We have that The expression for θ 1 (t 1 ) is obtained by integrating the expression for ε 1 (t 1 ) and using (38) to estimate ) dt 1 = 0 and θ 1 is bounded and 1-periodic.To estimate the transition time T 1 , the boundary constraint By (38), the second term on the right-hand side is O(1) and the third term can be estimated by T 1 O(R).
Rearranging yields the desired result.
Combining Lemmas 4.6 and 4.8 we obtain the following characterisation of the transition map Π (α) 1 : Σ in 1 → Σ out 1 , which summarises the dynamics in chart K 1 .

Dynamics in the Rescaling Chart K 2
In chart K 2 we study solutions close to the extension of the center manifold M a 2 = κ 12 (M a 1 ).Lemma 4.11.Following the singular time rescaling ρ 2 dt = dt 2 or equivalently, ρ 2 t = t 2 , the desingularized equations in chart K 2 are given by where by a slight abuse of notation we now write (•) ′ = d/dt 2 .Since 0 < ρ 2 = ε ≪ 1, system (39) can also be viewed as a perturbation problem in (r 2 , θ 2 , y 2 )-space as ρ 2 → 0.
Proof.This follows immediately after differentiating the defining expressions for local K 2 coordinates in (31) and applying ρ 2 t = t 2 .
In order to understand system (39), we first consider the limiting system as ρ 2 → 0, i.e.
There are two possibilities for the angular dynamics, depending on α, since We start with the case α = 2.In this case, we obtain a θ 2 -family of planar systems where a = a(θ 2 ), b = b(θ 2 ) and c = c(θ 2 ) define (parameter dependent) positive constants.The transformation leads to which is precisely the Riccati-type equation which arises within the K 2 chart in the analysis of the regular fold point in [42].For each fixed θ ∈ R/Z, solutions to system (43) (and therefore also (41)) can be written in terms of Airy functions, whose asymptotic properties are known [19,51].The properties that are relevant for our purposes are collected in [51] and reformulated in a notation similar to ours in [42].
The following result is a direct extension of the latter formulation; we simply append the existing result with the decoupled angular dynamics induced by the equation for θ 2 .
Proposition 4.12.Fix α = 2 and consider the restricted system (39)| ρ2=0 or, equivalently, the limiting system (40).The following assertions are true: (a) Every orbit approaches a two-dimensional horizontal asymptote/plane y 2 = y r from y 2 > y r as r 2 → ∞.The value of y r depends on the initial conditions.
(b) There exists a unique, invariant two-dimensional surface where h (2) y2 (r 2 ) is smooth with asymptotics and Ω 0 is the constant defined in Theorem 3.2. .
(c) All orbits with initial conditions to the right of γ 2 in the (r 2 , y 2 )-plane are backwards asymptotic to the paraboloid {(r 2 , θ 2 , (b/a)r 2 2 , 0) : r 2 ≥ 0, θ 2 ∈ R/Z}.(d) All orbits with initial condition to the left of γ 2 in the (r 2 , y 2 )-plane are backwards asymptotic to a horizontal asymptote/plane y 2 = y l .Specifically, y 2 (t 2 ) → y l from below and r 2 (t 2 ) → −∞ as t 2 → −∞.The value of y l depends on the initial conditions, but satisfies y l > y r for each fixed orbit.
The Riccati dynamics described in Proposition 4.12 are sketched for the decoupled planar system (41) in Figure 12, and for the three-dimensional limiting system (40) in Figure 13.
Proposition 4.13.Fix α = 1 and consider the restricted system (39)| ρ2=0 or, equivalently, the limiting system (40).For sufficiently small but fixed E, R > 0, solutions with initial conditions in finite time, where ν > 0 is a constant.In particular, we have that 1 is the unique center manifold described in Lemma 4.6.Proof.The idea is to bound solutions to the Riccati equation ( 47) between upper and lower solutions with known asymptotics as r 2 → ∞.Specifically, upper and lower solutions will be obtained as concatenations of solutions to Riccati equations of the form where λ, µ ̸ = 0 are suitably chosen constants.Note that it is equivalent to compare solutions of the two planar autonomous systems (45) and which has asymptotic properties similar to those of system ( 41) with a = λ, b = µ and c = 1 (the factor of c(t 2 ) does not effect the phase portrait since dr 2 /dy 2 has no explicit time dependence).We start by considering solutions of system (45) with initial conditions (r 2 (0), Our first task is to bound the r 2 -component of solutions to (45) on the interval t 2 ∈ [0, T 0 ], where T 0 > 0 is the unique solution to φ(T 0 ) = E −2 .If a solution exists until t 2 = T 0 (and has not blown up in a finite time smaller than T 0 ), the definition of T 0 corresponds to the intersection of solutions with {y 2 = 0} (recall that y 2 (t 2 ) = E −2 − φ(t 2 ), see (46)).We get a left-bound to the Riccati equation ( 47) by identifying a lower solution r 2 (t 2 ) ≤ r 2 (t 2 ).This is obtained by letting r 2 (t 2 ) be a solution of ( 49) with (r 2 (0), y 2 (0)) = (β − /ER, E −2 ) and λ, µ replaced by respectively.This ensures that r 2 (0) ≤ r 2 (0) and for all t 2 ∈ [0, T 0 ], as φ(t 2 )−E −2 ≤ 0 for t 2 ∈ [0, T 0 ].Solutions of equation ( 49) with λ = λ and µ = µ are described by Proposition 4.12 after a simple transformation similar to that in (42).In particular, there exists an invariant two-dimensional surface γ 2 with properties similar to the surface γ 2 , which intersects {y 2 = 0}.Choosing E and R sufficiently small guarantees that the initial condition r 2 (0) = β − /ER is smaller than the r 2 -coordinate of the intersection γ 2 ∩ {y 2 = E −2 }.This implies that the solution (r 2 (t 2 ), y 2 (t 2 )) also intersects {y 2 = 0}.Thus, using (47), we conclude that r 2 (t 2 ) is a lower solution for r 2 (t 2 ), with r 2 (t 2 ) ≤ r 2 (t 2 ) for all t 2 ∈ [0, T 0 ].An upper solution r 2 (t 2 ) with the initial condition (r 2 (0), y 2 (0)) = (β + /ER, E −2 ) can be constructed in a similar way, by identifying r 2 (t 2 ) as the solution of the Riccati equation ( 49) with λ, µ replaced by respectively.We have that r 2 (0) ≥ r 2 (0) and for all t 2 ∈ [0, T 0 ].There are now two different possibilities, depending on whether r 2 (t 2 ) blows up in finite time for r 2 → ∞ (case A), or exists up to the intersection with {y 2 = 0}, i.e. for all t 2 ∈ [0, T 0 ] (case B).Blow-up for r 2 → −∞ is not possible as solutions of equation ( 49) with λ = λ and µ = µ are described by Proposition 4.12 after a simple transformation similar to that in (42).In particular, there exists an invariant two-dimensional surface γ 2 with properties similar to the surface γ 2 , which intersects {y 2 = 0}.Choosing E and R sufficiently small guarantees that the initial condition r 2 (0) = β + /ER is larger than the r 2 -coordinate of the intersection γ 2 ∩ {y 2 = E −2 }, implying that r 2 cannot blow up as r 2 → −∞.Hence it is sufficient to study case A and case B in the following.
In case B, r 2 (t 2 ) exists for all t 2 ∈ [0, T 0 ] and hence intersects {y 2 = 0}.Using (47), we can conclude that r 2 (t 2 ) ≥ r 2 (t 2 ) for all t 2 ∈ [0, T 0 ].Combining this with the results for the lower solution r 2 (t 2 ), we have that In particular, along y 2 = 0 we have The situation is sketched in Figure 14 for case A and in Figure 15 for case B. The functions r 2 and r 2 (the latter only in case B) do not necessarily define lower and upper solutions when y 2 < 0 (i.e. when t 2 > T 0 ), since the inequalities ( 51) and ( 52) may no longer be satisfied.We can, however, connect to different lower and upper solutions which we denote by r2 (t 2 ) and r2 (t 2 ) respectively, with initial conditions r2 (T 0 ) = r 2 (T 0 ) and r2 (T 0 ) = r 2 (T 0 ), obtained as segments of solutions to (49) with λ, µ replaced by λ = A − , μ = B − , respectively.This time, we may apply Proposition 4.12 for the Riccati equation with a = λ, b = μ, c = 1 and a = λ, b = μ, c = 1, respectively.We find that r2 (t 2 ), r2 (t 2 ) → ∞, and that orbits of the corresponding planar systems approach horizontal asymptotes y 2 = y r± in the (r 2 , y 2 )-plane from above.Thus, r2 (t 2 ) and r2 (t 2 ) blow up in finite time.The relation y 2 (t 2 ) = E −2 − φ(t 2 ) leads to the following implicit equations for the blow-up times t ± 2 : transversally.The lower solution r 2 intersects {y 2 = 0} at r 2 (T 0 ), where it connects to the lower solution r2 , which also converges to a horizontal asymptote intersecting Σ out 2 transversally.Solutions to the Riccati equation ( 47) with initial conditions in κ 12 (N a 1 ∩ Σ out 1 ) ⊂ Σ in 2 are sketched in shaded red and have r 2 -coordinates which are bounded between r 2 and r 2 (r 2 and r 2 ) when y 2 ≥ 0 (y 2 ≤ 0).For fixed initial condition θ * 2 , a sample trajectory is shown in red.
. transversally.The lower solution r 2 intersects {y 2 = 0} at r 2 (T 0 ), where it connects to the lower solution r2 , which also converges to a horizontal asymptote intersecting Σ out 2 transversally.Solutions to the Riccati equation ( 47) with initial conditions in κ 12 (N a 1 ∩ Σ out 1 ) ⊂ Σ in 2 are sketched in shaded red and have r 2 -coordinates which are bounded between r 2 and r 2 (r 2 and r2 ) when y 2 ≥ 0 (y 2 ≤ 0).For a fixed initial condition θ * 2 , a sample trajectory is shown in red.
Remark 4.14.As pointed out by an anonymous referee, the proof of Proposition 4.13 can be simplified by fixing ν ≥ 1, in which case the existence of lower solutions r 2 and r2 and the fact that y ′ 2 < 0 are sufficient to show that solutions intersect Σ in 2 .The existence of upper solutions r 2 and r2 allows for greater control of the solutions when ν < 1, but it is not strictly necessary for the proof of Theorem 3.2.Propositions 4.12 and 4.13 can be used to describe the limiting behaviour of the y 2 -component of the map Π induced by the flow of initial conditions in which is related to the exit section in chart K 1 via Σ in 2 = κ 12 (Σ out 1 ), up to the exit section where ν > 0 is the constant in Proposition 4.13.
Similarly to the analysis in K 1 , the remaining components of the image Π (39).Then where c(θ * 2 ) = const.> 0 and φ(t 2 ) is the function defined in (46).Proof.The solutions are obtained by direct integration.
In the following we define the function h (Proposition 4.13 guarantees the existence of this intersection), and (ii) h (2) y2 defining the special Riccati solution in Proposition 4.12 (we simply suppressed the parameter dependence on θ 2 in the notation there).Combining Propositions 4.12-4.13and Lemma 4.15, we obtain the following characterization of the transition map Π Proposition 4.16.For sufficiently small but fixed E, R > 0, Π where r 2,c is the (generally Proof.We start with the case α = 2 and let (r 2,c , θ 2 , E −2 , 0) denote the K 2 coordinates of the intersection γ 2 ∩ Σ in 2 .It follows from Proposition 4.12 and Lemma 4.15 that where ) and c = c(θ 2 ) are constant functions of the initial value θ 2 here).The transition time for an initial condition with ρ 2 = 0 can be estimated using Lemma 4.15.In particular, Since the transition time T 2 is finite and system (40) is a regular perturbation problem, a neighbourhood of (r 2,c , θ The form of the map in Proposition 4.16 follows. Now fix α = 1.In this case the transition time T 2 satisfies 0 < T − 2 ≤ T 2 ≤ T + 2 < ∞, where T ± 2 are the transition times associated to the lower/upper solutions used in the proof of Proposition 4.13.The result from here follows similarly to the case α = 2, using the fact that T 2 is finite and that system (40) is a regular perturbation problem.
Figure 16: Geometry and dynamics projected into (r 2 , y 2 , ρ 2 )-space for α = 2; c.f. Figures 12 and 13.The extension of the (projected) three-dimensional manifold M a 2 is sketched in shaded green.The image of the wedge-shaped region κ 12 (Π is also shown in shaded green. The extension of M a 2 up to Σ out 2 , as described by Proposition 4.16, is shown in projections onto (r 2 , θ 2 , y 2 )-and (r 2 , y 2 , ρ 2 )-space for α = 2 in Figures 13 and 16, respectively.
Finally we note that if α = 2, then the expression for h 2 ) can be simplified using the following result, which can be found in [42].

Dynamics in the Exit Chart K 3
In chart K 3 we study solutions close to the extension of the manifold M a 3 = κ 23 (M a 2 ) as it leaves a neighbourhood of the singular cycle S c 0 .
Lemma 4.18.Following the positive transformation of time ρ 3 dt = dt 3 , the desingularized equations in chart K 3 are given by ρ where Proof.This follows after direct differentiation of the local coordinate expressions in (31) and subsequent application of the desingularization ρ 3 dt = dt 3 .
Except for the angular coordinate θ 3 , the analysis in K 3 is entirely local.Specifically, we focus on dynamics within the set where E, R and ν are the same positive constants used to define the entry and exit sections in charts K 1 and K 2 .System (54) has a circular critical manifold with the following properties.Lemma 4.19.Consider system (54).Q is normally hyperbolic and saddle-type.More precisely, the linearization along Q has eigenvalues and corresponding eigenvectors (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0) and (0, −αρ α−1 Proof.Direct calculation.Lemma 4.19 implies the existence of local center, stable and unstable manifolds at Q. Specifically, there is a one-dimensional center manifold W c 3 (Q), a two-dimensional stable manifold W s 3 (Q), and a one-dimensional unstable manifold W u 3 (Q).The geometry is sketched in different three-dimensional subspaces and projections in Figures 17,18  If α = 2, we can show that the extension of the surface γ 2 connects to the critical manifold Q tangentially to the cylinder spanned by the center and weakly stable eigenvectors.Lemma 4.20.Fix α = 2. Then the extension of γ 2 under system (54), i.e. γ 3 = κ 23 (γ 2 ), connects to Q tangentially to the cylinder segment {(0, θ 3 , 0, ε 3 ) ∈ D 3 }.
Having described the geometry for α = 2, we turn our attention to the characteristics of the map Π induced by the flow of initial conditions in , up to the exit section which is precisely the representation of the exit section ∆ out ε defined in (29) after blow-up in chart K 3 if we set y 0 = R 2 ν.Similarly to the K 3 analysis for the regular fold point in [42], Π (α) 3 can be analysed directly after a second positive transformation of time F (ρ 3 , θ 3 , y 3 , ε 3 )dt 3 = d t3 , which amounts to division of 1 (shown in shaded green in Figure 16) is shown here in Σ in 3 (again in shaded green), along with its image in Σ out 3 under Π (2) 3 .The projection of the three-dimensional manifold M a 3 with base along γ 3 ∪ {y 3 = ε 3 = 0, ρ 3 ≥ 0} is also shown.the right-hand side in system (54) by F (ρ 3 , θ 3 , y 3 , ε 3 ) (which is strictly positive in D 3 ).This leads to the orbitally equivalent system where by a slight abuse of notation the prime notation now refers to differentiation with respect to t3 .
System ( 55) is used to obtain the following result.
One can show with direct estimates that the term in parentheses is uniformly bounded by a positive constant, which implies that 3 e t3 , it follows that θ3 can be written in terms of ρ 3 , specifically, we obtain where we permit a slight abuse of notation in switching the argument from t3 to ρ 3 .Substituting ρ 3 ( t3 ) = ρ * 3 e t3 and evaluating the resulting expression at t3 = T 3 yields the desired estimate for θ 3 (T 3 ) = h(2) θ3 (ρ 3 , θ 3 , y 3 ) mod 1.It remains to estimate y 3 (T 3 ).It will be helpful to use (58) in order to write where a 3 = a(θ * 3 ), b 3 = b(θ * 3 ) and c 3 = c(θ * 3 ).From here, we use and adaptation of the 'partial linearisation' method in the proof [42,Prop. 2.11].Define a new variable ε3 = ε 3  3 , and consider the resulting equations within the invariant hyperplane {ρ 3 = 0}: which may be considered as a θ * 3 -family of planar systems in the (y 3 , ε3 ) variables, with constants a 3 , b 3 and c 3 depending on the parameter θ * 3 .Since the Jacobian at the equilibrium (0, 0) is hyperbolic and non-resonant with the eigenvalues −2 and −3, there exists a parameter-dependent, near-identity transformation of the form ỹ3 = ψ(y 3 , ε3 , θ * 3 ) = y 3 + O(y 3 ε3 ), such that the transformed system has been linearized, i.e.

Proof of Theorem 3.2
We now combine the results obtained in the charts K 1 , K 2 and K 3 in order to prove Theorem 3.2.The idea is to describe the extended map π defined in Section 4.1, the first three components of which coincide with components of π (α) : ∆ in → ∆ out , via its preimage in the blown-up space: This can be done explicitly using the change of coordinate maps κ ij in Lemma 4.1 and the characterisation of Π in Propositions 4.9, 4.16 and 4.21, respectively.The geometry is sketched in Figure 20.
In order to prove Assertions (a)-(c), we need to derive the form of the map Π (α) .Initial conditions (r 1 , θ 1 , R, ε 1 ) ∈ Σ in 1 are mapped to Σ out 1 under Π (α) 1 as described by Proposition 4.9.Since Σ in 2 = κ 12 (Σ out 1 ), we obtain the following input for the map Π where we recall that h and c 0 := 1 0 c(θ 1 (t 1 )) dt 1 .Using (64), Lemma 4.1 and Proposition 4.16, we obtain an expression for the input to the map Π where we used the fact that r 2 − r 2,c = O(e − ρ/ε 3 1 ) in the derivation of the third component in the right-hand side, and h where c(θ 2 ) = c h (2) θ1 (r 1 , θ 1 , ε 1 ) .With regards to the third component in (66), if α = 1 then we have where ν > 0 is the constant appearing in the definition of Σ out 2 (recall Proposition 4.13).If α = 2, then by Lemma 4.17 it holds h (2)  y2 θ1 (r 1 , θ 1 , ε 1 ) .Substituting the right-most expression in (66) into the expression for Π (α) 3 (ρ 3 , θ 3 , y 3 , E) in Proposition 4.21 yields an expression for Π (α) (r 1 , θ 1 , R, ε 1 ).We obtain and specify the second and third components in the following.We have h where we used ( 65) and (67) in order to obtain the right-hand side.It remains to specify the third component in the right-hand side of (68).If α = 1, we can use the bounds in (56) to estimate for constants ν > 0 and σ − < σ + , such that h (1) where α ∈ N + and ϑ > 0 is a parameter such that ϑ = O(1) as ε → 0. Note that this is equivalent to considering the scaling regimes defined by After one more change of time-scale t = ϑ t, our system becomes where by a slight abuse of notation, we now allow the dash to denote differentiation with respect to the time t.Finally, we assume the following properties of the functions K, g and A.
Assumption 3.There exists a point x F such that Assumption 3 implies that the critical manifold, given by S 0 = {(x, θ, K(x)) : x ∈ R, θ ∈ R/Z}, contains a regular folded cycle S c 0 = {(x F , θ, K(x F )) : θ ∈ R/Z}.The set S c 0 can be moved to the 'origin' via the coordinate translation (x, ỹ) = (x − x F , y − K(x F )), which yields the system which is -up to relabelling (x, ỹ) ↔ (r, y) -in the general form (22) with all of which are strictly positive due to Assumption 3. Using the same relabelling, we define entry and exit sections ∆ in and ∆ out analogously to those in ( 23) and (24), respectively, and let π (α) : ∆ in → ∆ out denote the transition map induced by the forward flow of system (73) for a given α ∈ N + .
Note that via (71), the results in Theorem 5.2 can be reformulated and restated for the scaling regimes defined by (72).As for our main results in Section 3, the results for scalings corresponding to semi-oscillatory dynamics (α ∈ {1, 2}) cannot be obtained using established theory for two time-scale systems.These results are new, to the best of our knowledge, even for the forced van der Pol example (the intermediate frequency regime considered in [8] covers the case α = 3/2).
The reduced problem on S 0 will differ depending on the scaling, i.e. depending on the value of α.We are primarily interested in the cases α ∈ {1, 2}, for which 0 < δ ≪ ω ≪ 1; see however Remark 5.8 below for the cases α ≥ 3. Fixing α ∈ {1, 2}, rewriting system (76) on the slow time-scale τ = ε α t, and taking the limit ε → 0 yields the reduced problem This leads to the following reduced vector field on S 0 , expressed in the (z, θ)-coordinate chart: It is typical for problems of this kind to consider the so-called desingularised reduced problem [45,58,64].This may be obtained after a (singular) time transformation which amounts to multiplication of the vector field by z (see again Remark 2.2), leading to System ( 80) is orbitally equivalent to system (79) on S 0 \ S c 0 , but the orientation of orbits on S r 0 (where z < 0) is reversed.Direct calculations reveal the presence of two equilibria along S c 0 , namely p s : (0, 1/4), p c : (0, 3/4), which are of neutral folded saddle and folded center type, respectively (see [58] for definitions), see Figure 21.Moreover, system (80) is Hamiltonian, with Hamiltonian function Except for the points p s and p c , there are three different types of orbits shown in Figure 21.These are identified in the following result.79) on S 0 , projected into the (θ, z)-plane.We sketch the case with A = 2 and R = 1, where we recall that R > 0 is the constant used to define ∆ in .The normally hyperbolic and attracting (repelling) branch S a 0 (S r 0 ) is indicated in shaded turquoise (gray).The fold cycle S c 0 which separates these two branches is sketched in blue.There are two folded singularities on S c 0 : a neutral folded saddle p s , and a folded center p c .The singular true and faux canards through p s coincide in this case (denoted γ c and shown in green), forming a limit cycle of period 2, and the region bounded outside of γ c is foliated by periodic orbits of period 1.Finally, we sketch the forward evolution of an initial condition in ∆ in ∩ S a 0 with initial angle θ * ∈ I (the interval I is shown in blue), i.e. inside the region bounded by γ c .Such points reach a regular jump point on S c 0 with angle θ e ∈ (3/4, 1) ∪ [0, 1/4) in finite time.For A = 2 and R = 1 we have θ l = 5/12, θ r = 1/12.The initial condition in red has angle θ * = 5/6, for which our results imply a jump angle of θ e (θ * ) = 11/12.
We can use the information obtained above to describe the local dynamics near S c 0 .Specifically, we want to describe the transition map π (α) : ∆ in → ∆ out induced by the forward flow of system (76), where for small but fixed R, β, a 0 > 0. A singular transition map π (α) 0 : ∆ in → ∆ out for ε = 0 can be constructed explicitly, by concatenating solutions to the layer and reduced problems described above.If we choose R > 0 so that then ∆ in ∩ S a 0 intersects the singular canard solution γ c at the two points (R, θ l ) and (R, θ r ), where see Figure 21.Let I ⊂ R/Z denote the interior of the interval between θ l and θ r which is contained within the region bounded by γ c .Different types of singular orbits can be constructed depending on whether or not one starts with an initial angle θ ∈ I, θ ∈ {θ l , θ r } or θ / ∈ I. Consider an initial condition (z * , θ * , R 2 ) ∈ ∆ in with θ * ∈ I.This is mapped by the layer flow of system (77) to the point (R, θ * , R 2 ) ∈ ∆ in ∩ S a 0 .Since θ * ∈ I, the point (R, θ * , R 2 ) lies on an orbit of the desingularized reduced problem (80) which is contained within a level curve H(z, θ) = K ∈ (−A/2, A/2).In fact, one can show directly that Proposition 5.5 describes regular jump type dynamics for solutions with initial conditions in ∆in , which can be chosen arbitrarily close to ∆ in | θ∈I by decreasing ε 0 if necessary.It is worthy to note that the leading order approximation for the 'jump angle' π(α) θ (z, θ, ε) ∼ θ e (θ) can be calculated explicitly as a function of the initial angle using (84).It is also worthy to note that for α ∈ {1, 2} we have π(α) a (z, θ, ε) = O(ε 2α/3 ), in contrast to the O(ε 2 ) estimates for the parameter drift in Theorem 3.2.This is a consequence of the fact that the 'parameter drift' in system (76) in a occurs on the intermediate, rotational time-scale.This is not the case for systems satisfying Assumptions 1-2.
We may also consider the case in which initial conditions are chosen 'higher up', in an annular section where the constant ρ > 0 satisfies c.f. equation (81).A subset of solutions with initial conditions in Σ in are described by the following result.We formulate it in terms of the transition map Π (α) : Σ in → ∆ out .
Remark 5.7.Due to the transcendental nature of equation (86), a closed form expression for the function T = T (θ, ε) can only be given in terms of special functions.Nevertheless, for any given θ, T (θ, ε) be calculated to arbitrary precision with standard numerical methods (e.g.Newton or bisection methods).
Propositions 5.5 and 5.6 describe solutions of 'regular jump type', but they do not explain what happens to solutions with initial conditions in ∆ in \ J 1 .This corresponds to a set of initial conditions that is 'small' in these sense that the diameter of the interval [0, 1) \ I tends to zero as R → 0, but for fixed R > 0 it is still O(1) with respect to ε → 0.Moreover, this set contains solutions close to the singular canard γ c , which are expected to play an important role in determining the qualitative dynamics [58].A detailed description of the dynamics in this case is left for future work.
Remark 5.8.For α ≥ 3, the dynamics can be analysed and understood entirely with established theory for two time-scale systems.The case α = 3 is non-trivial, and features a degenerate bifurcation in the desingularized reduced problem as the folded saddle p s and neutral folded center p c collide and annihilate at a particular value of the forcing A = A c > 0. For A < A c , the folded cycle S c 0 is regular, i.e. there are no folded singularities.For α ≥ 4, S c 0 is always regular.A direct application of the results in [59] in this case yields an explicit expression for the transition map: π (α) z, θ, R 2 = −R, θ + O(ε α ln ε), O(ε 2α/3 ) , α ≥ 4.
Remark 5.9.The local geometry and dynamics near the fold lines/circles in the forced van der Pol system considered in e.g.[8,21,23] are characterised by alternation of folded saddles and folded focus singularities, in contrast to the alternation of folded saddle and folded center singularities in Figure 21.
Nevertheless, there are many similarities, including the fact that the image of S a ε on ∆ out is non-circular as ε → 0, due to the folded saddles.This feature, which is not present for the systems described by Theorem 3.2, is utilised in the construction of horseshoes in the forced van der Pol system [23].

Summary and Outlook
GSPT is an established and powerful tool for the analysis of stationary fast-slow systems, particularly when combined with the blow-up method.The corresponding theory for oscillatory fast-slow systems which possess a limit cycle manifold instead of (or in addition to) a critical manifold is less developed, despite the fact that oscillatory systems are common in applications.One reason for this appears to be that the theory for stationary fast-slow systems, notably Fenichel-Tikhonov theory and the blow-up method, depend on quasi-equilibrium properties possessed by stationary but not oscillatory fast-slow systems.The main purpose of this article has been to show that both Fenichel-Tikhonov theory and the blow-up method can be applied to study the dynamics of a class of multiple time-scale systems with three or more time-scales, as long as the angular/rotational dynamics are sufficiently slow relative to the radial dynamics.The class of systems considered, namely those in the general form (4), contains a 'semioscillatory' class of systems that are oscillatory with respect to the partial singular limit ε 1 > 0, ε 2 → 0, but stationary with respect to the double singular limit (ε 1 , ε 2 ) → (0, 0).Although our approach is not applicable in the purely oscillatory case, our analysis showed that Fenichel-Tikhonov theory and the blow-up method can be applied directly in the semi-oscillatory case.
More concretely, we focused on the dynamics near a (three time-scale) global singularity corresponding to a kind of folded cycle bifurcation in the layer problem obtained in the partial singular limit.In order to do so we derived a prototypical system for the three time-scale global singularity near the non-hyperbolic cycle S c 0 , recall Proposition 2.5, and studied the transition map induced by the flow.Our main result is Theorem 3.2, which is stated for system (22), which provides a rigorous characterisation of the asymptotic and strong contraction properties of the flow near S c 0 .The asymptotics h (α) y (θ, ε) = O(ε 2 ) and the strong contraction property are the same as for the stationary regular fold point studied in [42,51,59], not only in the semi-oscillatory cases α ∈ {1, 2} considered herein in detail, but for all α ∈ N + .The semioscillatory cases α ∈ {1, 2} are distinguished from the cases α ≥ 3 by the fact that the leading order estimates for h (α) y (θ, ε) are, in general, non-constant functions of θ.We derived an explicit expression for the leading order estimate in terms of the functions a(θ), b(θ) and c(θ) appearing in system (22) when α = 2 in particular.In addition, we provided detailed estimates for the exit angle, recall (25), as well as estimates for the total number of rotations about the y-axis as solutions traverse the region between ∆ in and ∆ out , recall Corollary 3.3.Theorem 3.2 is proved using the blow-up method in Section 4, thereby demonstrating the applicability of geometric blow-up techniques to study global singularities of limit cycle type in suitable classes of three time-scale systems.We focused on a detailed presentation of the proof in the semi-oscillatory cases α ∈ {1, 2}.The proof for α ≥ 3 can be obtained either via straightforward adaptations of the proof for α = 2, or via straightforward adaptations of the proof of the established results for two time-scale systems in [42,51,59].The primary complication in cases α ∈ {1, 2} stems from the fact that the dynamics remain 'global' in θ after blow-up near S c 0 .Consequently, one has to consider a blow-up of the entire cycle (as opposed to a local blow-up about a point on S c 0 ).We also demonstrated the applicability of our results by using them to derive detailed geometric and asymptotic information about the passage through a folded cycle singularity in a class of periodically forced Liénard equations; recall Section 5.1 and in particular Theorem 5.2.Finally, in Section 5.2, we considered a toy model for the study of tipping phenomena.Due to the periodic forcing in the fast equation, this model violates both of our primary Assumptions 1-2.The partial geometric analysis presented herein provides reason to argue that cases of this kind can for the most part be analysed using classical two time-scale theory.
In summary, the current manuscript provides a precedent for the effectiveness of stationary methods including Fenichel-Tikhonov theory and the blow-up method in semi-oscillatory multiple time-scale systems with at least three time-scales.There are many interesting directions which one might take from here.Our analysis of the toy model for tipping considered in Section 5.2 is incomplete, and motivates the analysis of a more general class of systems with regular folded cycle manifolds which violate our main Assumptions 1-2.One might also consider the applicability of the methods developed herein for other global cycle singularities for which there is no direct and lower dimensional analogue in stationary fast-slow systems, e.g.near a non-hyperbolic flip/period-doubling type cycle in three time-scale systems of the form (4). We also expect the detailed estimates on the local dynamics near regular folded cycles in Theorem 3.2 to play an important role in the description of multi-scale oscillations which involve regular folded cycles as a key ingredient.In [59], the authors derive a detailed characterisation of the local dynamics near a regular fold curve using geometric blow-up, and combine it with Fenichel theory in order to derive the form of the Poincaré map associated to prototypical 1-fast 2-slow systems with an S-shaped critical manifold.These results were used to prove the existence of an invariant torus for a certain parameter regime in the forced van der Pol equation, which contains an attracting relaxation oscillation if a particular circle map induced on the angular variable is contracting.We conjecture that similar constructions could be applied to the three time-scale semi-oscillatory systems considered herein, under suitable assumptions on the global geometry.The details of these and other related problems are left for future work.

Figure 3 :
Figure 3: Projected geometry and dynamics in the (r, y)-plane, as described by Theorem 3.2.The critical manifold and its submanifolds are sketched in colours consistent with earlier figures for the particular choice of a(θ) and b(θ) defined in Remark 2.9.The entry, exit sections ∆ in , ∆ out (magenta, orange) and the (extended) Fenichel slow manifolds S a ε , S r ε are also shown as shaded regions in green and light gray, respectively.For S a ε and S a 0 additionally a sample trajectory for a fixed θ initial condition is shown.

Figure 6 :
Figure 6: Case (C1).Sketch of the flow of system (22) for α = 2.The solution sketched in red makes N(2) rot (r, θ, ε) = ⌊O(ε −1 )⌋ rotations about the y-axis before leaving the neighbourhood close to S a ε ∩ ∆ out at an angle approximated by the expression for h

Figure 7 :
Figure 7: Case (C2).Sketch of the flow of system (22) for α = 3.The solution sketched in red makes N

Figure 8 :
Figure 8: Case (C3).Sketch of the flow of system (22).The solution sketched in red makes N

2 y 2 r 2 γ 2 1Figure 12 : 2 2
Figure 12: Dynamics for the Riccati equation (41) in the (r 2 , y 2 )-plane for α = 2 and a fixed choice of θ 2 (θ 2 is a parameter when α = 2 and ρ 2 = 0).A qualitatively similar figure is obtained for each θ 2 ∈ R/Z.The distinguished solution γ 2 with asymptotics described by Proposition 4.12 is shown in red.The parabola y 2 = (b(θ 2 )/a(θ 2 ))r 2 2 which separates solutions with different asymptotic properties is also shown.Projections of the entry and exit sections Σ in 2 and Σ out 2 are shown in yellow and cyan, respectively.

Figure 14 :
Figure 14: Case A of the Riccati dynamics for α = 1.The upper solution r 2 does not intersect with {y 2 = 0} and converges to a horizontal asymptote intersecting Σ out 2

Figure 15 :
Figure 15: Case B of the Riccati dynamics for α = 1.The upper solution r 2 intersects {y 2 = 0} at r 2 (T 0 ), where it connects to the upper solution r2 , which converges to a horizontal asymptote intersecting Σ out 2

Figure 17 :
Figure 17: Geometry and dynamics within D 3 , shown in (r 3 , ρ 3 , ε 3 )-space for α = 2.The (projection of the) two-dimensional surface γ 3 = κ 23 (γ 2 ) shown in red connects to the circular saddle-type critical manifold Q, which is indicated by the purple dot at the origin.Entry and exit sections Σ in 3 and Σ out 3 are shown in shaded cyan and orange, respectively.The wedge-shaped region within Σ out 2

Lemma 5 . 4 .
Consider the reduced problem (79).Orbits are contained within constant level sets H(z, θ) = K ≥ −A/2, where 1.K = −A/2 corresponds to the folded center p c .2.K ∈ (−A/2, A/2) correspond to orbits that connect to the folded cycle S c 0 in both forward and backward time.

3 .
K = A/2 corresponds to two saddle homoclinic orbits in the desingularised reduced problem (80), and a periodic orbit γ c of period τ = 2 along the true and faux canards through the neutral folded saddle p s in the reduced problem (79).

4 .Figure 21 :
Figure21: Dynamics for the reduced problem (79) on S 0 , projected into the (θ, z)-plane.We sketch the case with A = 2 and R = 1, where we recall that R > 0 is the constant used to define ∆ in .The normally hyperbolic and attracting (repelling) branch S a 0 (S r 0 ) is indicated in shaded turquoise (gray).The fold cycle S c 0 which separates these two branches is sketched in blue.There are two folded singularities on S c 0 : a neutral folded saddle p s , and a folded center p c .The singular true and faux canards through p s coincide in this case (denoted γ c and shown in green), forming a limit cycle of period 2, and the region bounded outside of γ c is foliated by periodic orbits of period 1.Finally, we sketch the forward evolution of an initial condition in ∆ in ∩ S a 0 with initial angle θ * ∈ I (the interval I is shown in blue), i.e. inside the region bounded by γ c .Such points reach a regular jump point on S c 0 with angle θ e ∈ (3/4, 1) ∪ [0, 1/4) in finite time.For A = 2 and R = 1 we have θ l = 5/12, θ r = 1/12.The initial condition in red has angle θ * = 5/6, for which our results imply a jump angle of θ e (θ * ) = 11/12.