A LIN’S METHOD APPROACH FOR DETECTING ALL CANARD ORBITS ARISING FROM A FOLDED NODE

. Canard orbits are relevant objects in slow-fast dynamical systems that organize the spiraling of orbits nearby. In three-dimensional vector ﬁelds with two slow and one fast variables, canard orbits arise from the intersection between an attracting and a repelling two-dimensional slow manifold. Special points called folded nodes generate such intersections: in a suitable transverse two-dimensional section Σ, the attracting and repelling slow manifolds are counter-rotating spirals that intersect in a ﬁnite number of points. We present an implementation of Lin’s method that is able to detect all of these inter- section points and, hence, all of the canard orbits arising from a folded node. With a boundary-value-problem setup we compute orbit segments on each slow manifold up to Σ, where we require that the corresponding end points in Σ lie in a one-dimensional subspace known as the Lin space Z . The Lin space Z must be transverse to the slow manifolds and it remains ﬁxed during the detection of canard orbits as zeros of the signed distance along Z . During the computation, a tangency of Z with one of the intersection curves in Σ may arise. To overcome this, we update the Lin space at an intermediate continu- ation step to detect a double tangency of Z to both curves in Σ, after which the canard detection is able to continue. Our method is demonstrated with the examples of the normal form for a folded node and of the Koper model.


1.
Introduction. Slow-fast systems arise as mathematical models in various applications and they describe physical phenomena, such as chemical reactions, nonharmonic oscillations, spiking and bursting [8,21,29,30,34,42,44,51]; their study has been an active area of research. All the natural questions from classic dynamical systems theory about the local and global configuration of phase space are valid for slow-fast systems, yet the existence of different time scales gives rise to new questions in the slow-fast context. In particular, such systems are known to exhibit mixed-mode oscillations, which are orbits of a vector field characterized by an alternation of both small and large-amplitude oscillations. Mixed-mode oscillations were first discovered in the Belousov-Zhabotinskii reaction [28,50], and since then have been found in a broad range of chemical and biological systems [9,14,19,47,53].
Solutions of slow-fast systems can be thought of as concatenations of slow motion and fast segments, where the fast segments are described by a layer problem and the slow motion is organized by slow manifolds [15,16,31,41]. These are locally invariant manifolds that, together with equilibria, periodic orbits and their corresponding invariant manifolds, organize the dynamics and the slow-fast nature of the system globally. The existence of slow manifolds as perturbations of the so-called critical manifold, formed by the intersection of the nullclines of the fast variables, is guaranteed by Fenichel theory [15,16,31], provided that a condition known as normal hyperbolicity is satisfied. We are particularly interested in three-dimensional slow-fast systems with two slow and one fast variables. In this setup, the critical manifold is a two-dimensional surface. Therefore, the corresponding slow manifolds are two-dimensional surfaces in the three-dimensional phase space that can be either attracting or repelling under the dynamics; determining their geometry in R 3 may be quite challenging. Two slow manifolds of different types may interact with each other. Generically, in our setup with two slow and one fast variables, intersections of an attracting slow manifold with a repelling slow manifold are structurally stable and give rise to canard orbits. These are orbits associated with trajectories of the slow flow that connect an attracting and a repelling sheet of the critical manifold by crossing through fold curves. Canard orbits have the unusual property that they follow a repelling slow manifold for a considerable amount of time [5,22,52]; as such, they are global objects that are responsible for the creation and/or organization of complex dynamics. The literature about canard orbits in R 3 is extensive; see for instance [2,5,12,23,39,40,52,54]. Note that canard orbits cannot be found analytically and must be computed numerically, which explains the keen interest in detecting them. While they may arise more widely, canard orbits in R 3 are closely related with folded singularities, which are singularities of the slow flow that are located on fold curves of the critical manifold. In particular, a folded node generates canard orbits, the number of which is determined by the eigenvalues of the folded node [52,54]. Such canard orbits organize small-amplitude oscillations of mixed-mode oscillations in many applied models [4,6,11,49,55]. Therefore, systems with folded nodes provide ideal and relevant test-case examples.
The use of advanced numerical techniques has been very successful for the computation and visualization of slow manifolds [17,24,25,36,37]. However, these methods do not allow for the detection of canard orbits straight away; in particular, it is generally not possible to find canard orbits by integrating a single orbit segment, because canard orbits are unstable objects in both forward and backward time. Up to now, the visualization of slow manifolds has been used for the approximate detection of canard orbits arising from folded nodes. This is done with a shooting approach by integration in forward or backward time, up to a suitable cross section, of initial conditions on curves that lie far from the fold on the critical manifold. The intersection points of the two intersection curves of the attracting and repelling slow manifolds are then found in the cross section by inspection. When this computation is performed with a continuation and boundary-value-problem setup, a correction step can be used to find canard orbits accurately [7,8]. However, this approach is ad-hoc in that each intersection point must be found by inspection and then corrected individually. This makes the detection of all canard orbits, for example, arising from a given folded singularity, cumbersome and time consuming.
In this paper we present an approach for the systematic detection of canard orbits as intersection curves of two-dimensional slow manifolds. We stress that this automatic detection step is the crucial first step. Once they are detected, canard orbits can then be continued in system parameters to study their properties and bifurcations; see, for example, [5,7,8,9,11,19,23,34,54]. We focus here on slow-fast systems where the canard orbits are generated by a folded node. More specifically, our method is able to detect in a systematic and automated fashion all of the canard orbits that (are know to) exist. It is based on the continuation of orbit segments that are solutions of two-point boundary value problems (2PBVP); these are solved with the package AUTO [10], which uses collocation in combination with pseudo-arclength continuation where the size of the continuation step is determined by considering the change along the entire orbit segment instead of merely the initial condition. This computational setup is able to cope very well with sensitive systems such as slow-fast systems [13]. Specifically, we simultaneously compute and continue two orbit segments -one on the attracting slow manifold and one on the repelling slow manifold -up to a section Σ through the folded node. These two orbit segments are coupled by requiring that their end points lie along a fixed direction. This Lin's method approach [38,43,46] is a way of defining a suitable test function that can be monitored: its zeros indicate that a canard orbit has been found; see already Fig. 1. Any or all of the canard orbits that are found in this way can then be continued in systems parameters in a straightforward manner: either as two orbit segments while keeping the test function at zero, or by first concatenating the two orbit segments to a single orbit segment. Such subsequent computations are beyond the scope of this paper; see [7,8,9,26] for examples and more details.
Lin's method was developed as an analytic technique for finding periodic or aperiodic solutions near heteroclinic or homoclinic cycles [27,32,43,48,56]. It considers one or several suitable sections transverse to the flow, and defines algebraic bifurcation equations from orbits with gaps in lower-dimensional subspaces. Simultaneous zeros of the gaps correspond to the global objects sought. Lin's method has been implemented also as a numerical technique for finding heteroclinic and homoclinic connections; see for instance [18,33,38,46]. More recently, the approach has been successfully applied in the slow-fast context; specifically, for detecting so-called connecting canard orbits arising as codimension-zero intersections between the twodimensional unstable manifold of a saddle-focus equilibrium and a two-dimensional repelling slow manifold in a model with a singular Hopf bifurcation [45]. Here, we follow a similar approach to compute canard orbits, which are codimension-zero objects since they are structurally stable in R 3 .
Lin's method requires the genericity condition that the Lin space is transverse to both the attracting and repelling slow manifolds; this is not a problem locally near the intersection point. However, finding all the canard orbits arising from a folded node is more of a global problem. As we show, tangencies of the Lin space with the intersection curves of the slow manifolds in Σ do occur and create a problem for the detection of all the intersection points; this is due to the counter-rotating nature of attracting and repelling slow manifolds near a folded node. In our approach, we detect these tangencies and, through an intermediate adjustment step, we update the Lin space in a suitable way. This allows the systematic detection of all the canard orbits arising from a folded node singularity in a sequence of continuation steps.
This paper is organized as follows. Section 2 gives the necessary background on slow-fast systems and canard orbits. Section 3 describes the general numerical setup for the computation of slow manifolds and the implementation of Lin's method. Section 4 shows our approach implemented for a normal form of a folded node; here we use the symmetry of the system to define the Lin space Z. Section 5 describes the situation for Z in general position, when during the detection of canard orbits one encounters tangencies of Z with the intersection curves of the slow manifolds with the section. In section 6 we present our overall method that deals with tangencies and allows for the systematic detection of all canard orbits. Section 7 demonstrates our approach with the Koper model for idealized chemical reactions. We end with a discussion in section 8.

2.
Background: Fenichel theory and canard orbits in R 3 . We now present some background on slow-fast systems in R 3 with two slow and one fast variables, as needed for this paper. For further details, we refer the interested reader to, for example, [1,9,15,16,31,41].
We consider a slow-fast vector field of the form where f , g 1 and g 2 are smooth functions and λ ∈ R k is a vector of parameters. Here, 0 < ε 1 represents the ratio of time scales, so that the variable x ∈ R is fast and the variables y, z ∈ R are slow. (In this paper we only encounter the case that f , g 1 and g 2 do not depend on ε.) Solutions of slow-fast systems can be thought of as a concatenation of slow motion with fast segments. Considering system (1) for ε = 0 gives the slow flow or reduced system    0 = f (x, y, z, λ), y = g 1 (x, y, z, λ), z = g 2 (x, y, z, λ), for the limiting slow motion. It is a differential-algebraic equation (DAE) where the constraint on the first equation defines the critical manifold which is the nullcline of the fast variable x. The dot in system (1) represents differentiation with respect to time on the slow time scale τ . One can rewrite system (1) with respect to the fast time scale t via a time rescaling by 1/ε to obtain    x = f (x, y, z, λ), y = εg 1 (x, y, z, λ), z = εg 2 (x, y, z, λ), where the prime denotes the derivative with respect to time on the fast time scale. Fast segments of solutions of (1) are approximated by solutions of the layer equations which is a family of differential equations on the fast time scale, obtained as the singular limit of (3) for ε = 0. Here the x -equation depends on y and z, which are now treated as parameters. Note that the critical manifold S is a manifold of equilibria for the fast subsystem.
The properties of the critical manifold S come from the fast subsystem. Accordingly, we say that a subset N ⊂ S is normally hyperbolic if all its points are hyperbolic equilibria of the fast subsystem x = f (x, y, z, λ). In other words, N ⊂ S is normally hyperbolic if, for all points (p x , p y , p z ) ∈ N , the Jacobian D x f (p x , p y , p z , λ) has no eigenvalues with zero real part. Since system (1) has a single fast variable, the normal hyperbolicity is reduced to f x (p x , p y , p z , λ) = 0; hence, system (4) implies that the critical manifold S may have parts that are either attracting or repelling. More precisely, the attracting sheet is S a := S ∩ {f x (x, y, z, λ) < 0}, and the repelling sheet is The sheets S a and S r of S may meet at fold curves that are defined by Note that normal hyperbolicity of S occurs away from the set F .
The fast dynamics of (1) is well understood by analyzing system (4); its solutions are one-dimensional fast fibers that are attracted to or repelled from S. On the other hand, the slow dynamics deserves a more detailed analysis. Since the reduced system (2) is restricted to its critical manifold, one can use the normal hyperbolicity of S away from fold curves and apply the Implicit Function Theorem to describe S locally as a graph x = φ(y, z) and, thus, obtain a two-dimensional system projected onto the plane of slow variables Unfortunately, S is not a graph over the slow variables near F . Alternatively, to study the dynamics on S, one can choose, say, x and z as the defining variables, and use the constraint f = 0 and the equations forẏ andż in (2) to obtain the system Here, the functions g 1 , g 2 , f x , f y , and f z all depend on the point (x, y, z) and the parameter λ, and (x, y, z) ∈ S, that is, f (x, y, z, λ) = 0. This formulation is valid everywhere on S, although the system is singular at fold curves. System (6) can be desingularized by rescaling time by the factor −1/f x . This way, one obtains which allows the extension of (6) to fold curves. Note that in system (7) the flow on the repelling sheet S r is reversed. Generically along a fold curve, trajectories of (2) approach F in either forward or backward time on both the attracting and repelling sheets S a and S r of S. Singularities of the desingularized system (7) lie on F and are known as folded singularities. At such points trajectories of the slow flow (2) pass from S a to S r . A point q = (q x , q y , q z ) ∈ F is a folded singularity if The stability of a folded singularity q comes from the analysis of q as a singularity of (7). Let λ 1 and λ 2 denote the eigenvalues of the Jacobian matrix of the desingularized system (7) at q. We call q • a folded saddle, if λ 1 λ 2 < 0 and λ 1 , λ 2 ∈ R.
Note that folded singularities are typically not singularities of the full system (1), and they are only defined for the desingularized system (7). In this paper we focus on folded nodes, which are associated with the existence of a finite number of canard orbits in the full system (1) for 0 < ε 1; see [52,54] for a detailed analysis of folded singularities. These results allow us to determine the number of canard orbits that exist in a particular system for specific parameter values, so that we can then apply our method and test whether it computes them all.
For 0 < ε 1, Fenichel Theory [15,16] guarantees the existence of attracting and repelling smooth slow manifolds S a ε and S r ε in the full system (1) that lie at distance O(ε) away from S a and S r where S is normally hyperbolic, that is, away from fold curves. Trajectories of system (1) with ε > 0 are attracted to S a ε and repelled from S r ε in forward time at fast exponential rates; trajectories that lie on a slow manifold remain slow for an O(1) time on the slow time scale. Slow manifolds are not unique, but the distance between a pair of slow manifolds of the same type is of order O(exp(− c ε )) for some c > 0 [31]. Slow manifolds can be extended in forward and backward time by the flow; however, their behavior is not controlled by the singular limits (2) and (4). In particular, one can extend slow manifolds close to folded singularities, where Fenichel theory does not apply and slow manifolds are no longer approximations of the corresponding sheets of the critical manifold; attracting and repelling slow manifolds may exhibit complex oscillations in a neighborhood of a folded node and start interacting. In this paper, the slow manifolds are two-dimensional surfaces that intersect in canard orbits, which remain on S r ε for an O(1) time, in contrast to most trajectories of (1), which jump at folds along fast fibers.
Canard orbits in R 3 have been classified and analyzed in [2,4,52,54] by using Geometric Singular Perturbation Theory and blow-up techniques. Generically, for a folded node q, one has an inequality of the form λ s := |λ 1 | > λ w := |λ 2 | for its eigenvalues. The corresponding eigendirectionsγ s andγ w are referred to as the strong and weak singular canards, respectively. The ratio λ w /λ s < 1 between the weak and the strong eigenvalues of q determines the number of secondary (maximal) canard orbits that arise as additional transverse intersections between S a ε and S r ε for 0 < ε 1; see [52,54].
3. Numerical setup. We now describe the numerical techniques we use for the computation of slow manifolds via the continuation of solutions to a 2PBVP implemented in AUTO [10], as well as the numerical setup for the detection of canard orbits with Lin's method. For further background information we refer the interested reader to, for example, [8,9,36,38,46]. Throughout, instead of a slow-fast system of the form (1), we consider its equivalent version (3) written in the fast time scale, which contains the ratio of time scales ε as part of the right-hand side of the equation. As is standard in AUTO, we rescale time and write system (3) in the forṁ where u = (x, y, z) ∈ R 3 , the function F : R 3 × R k → R 3 corresponds to the righthand side of system (3), and λ ∈ R k is a vector of parameters. Here, any orbit segment is parameterized over the unit interval [0, 1] and T is the actual integration time, which is considered as a separate parameter. We assume that system (1) and its rescaled version (9) has a folded node at the point p 0 ∈ R 3 for λ = λ 0 , so that we can calculate the eigenvalue ratio to know beforehand the number of canard orbits that are expected to exist. The goal is to find intersections between S a ε and S r ε by looking at their corresponding intersection curves with a cross-section Σ that is transverse to the critical manifold S. Since S a ε and S r ε are expected to spiral and interact near the folded node p 0 , we stipulate that p 0 ∈ Σ and write Σ as where Y Σ is a two-dimensional subspace of R 3 that is normal to the fold curve containing p 0 . We consider the intersection sets S a ε := S a ε ∩ Σ and S r ε := S r ε ∩ Σ. To obtain S a ε and S r ε we need to compute the slow manifolds up to Σ and track their intersection sets. In this setup, the relevant part of the attracting slow manifold S a ε is a family of orbit segments that are solutions to the 2PBVP   u where the line L a lies on the attracting sheet S a of S, parallel to and sufficiently far away from the fold curve F , that is, L a is transverse to the flow. Solutions of (11) provide an accurate approximation of S a ε , since Fenichel theory [15,16,31,41] ensures, for ε > 0 small enough, the existence of slow manifolds as O(ε) perturbations of the corresponding sheets of the critical manifold away from fold curves. The end point u a (1) of a solution of (11) lies on the intersection curve S a ε , which is the diffeomorphic image of L a under the flow. We need to find a 'good' initial orbit segment that satisfies the 2PBVP (11), which is to be continued along L a to compute S a ε and S a ε . We follow [8,9] and compute an initial orbit in two homotopy steps, which we explain briefly here. Starting from the folded node p 0 , that is, from the trivial solution u a ≡ p 0 for T a = 0, we keep the condition u a (1) ∈ Σ throughout but allow u a (0) to move along the fold curve F while T a varies. When it is far enough from p 0 , we allow u a (0) to move on the attracting sheet S a away from F until the chosen line L a is reached.
Similarly, the 2PBVP   u has as its solutions a family orbit segments that approximate S r ε , where the line L r lies on the repelling sheet S r of S away from F and is transverse to the flow; the end points u r (0) describe the intersection curve S r ε , which is now the diffeomorphic image of L r under the backward-time flow; see [7,8,9,36] for further details on the computation of slow manifolds.
Having initial solutions to the 2PBVPs (11) and (12), we can set up Lin's method to find intersections between S a ε and S r ε on Σ and, therefore, intersections between S a ε and S r ε . The main idea of Lin's method is that one couples the 2PBVPs (11) and (12) to construct a new 2PBVP that has as its solution two orbit segments, one from L a to Σ and the other from Σ to L r . The difference between their two end points in Σ is chosen in a codimension-one subspace, called the Lin space Z; this gives rise to a well-defined test function, called the Lin gap η. A zero of the Lin gap η corresponds to a canard orbit. In other words, a canard orbit can be thought of as a connection from L a to L r , which is found as a concatenation of an orbit segment from L a to Σ and another orbit segment from Σ to L r .
To get started, we choose initial orbit segments u a and u r on S a ε and S r ε that are solutions of (11) and (12), respectively, and consider their corresponding end points p a := u a (1) and p r := u r (0) in Σ. We then use p a and p r to define the unit vector as well as its normal unit vector n Z ∈ Σ, with n Z ⊥ v Z . The vector v Z spans the Lin space, that is, Note that there are many ways to define Z, given by different choices for p a and p r . The genericity condition for the Lin space Z requires that Z is transverse to the intersection curves S a ε and S r ε , which are unknown at the beginning of the calculation. However, Z will generically be transverse to both S a ε and S r ε . Once it is defined, the vector v Z remains fixed during the continuation and its normal vector n Z ∈ Σ is obtained from the equation v Z · n Z = 0.
In order to cope with the counter-rotating spirals S a ε and S r ε it is convenient to represent v Z and n Z as follows. One can always find an orthonormal basis on the subspace Y Σ that defines Σ in (10), so that v Z := cos(2πα) sin(2πα) and n Z := − sin(2πα) cos(2πα) (13) in that basis. Here, the single parameter α ∈ [0, 1] parameterizes v Z and n Z , and it is initialized by the choices of p a and p r ; fixing v Z and n Z means fixing α. We stress that, for the detection of a single canard orbit with Lin's method, the vectors v Z , n Z and the parameter α remain fixed during the continuation. We now look for solutions to the coupled problems (11) and (12) with the additional boundary conditions (u a (1) − u r (0)) · n Z = 0, v Z · (u a (1) − u r (0)) = η.
The boundary condition (14) ensures that the end points u a (1) and u r (0) lie along the Lin space Z during the continuation of the overall 2PBVP, and (15) defines the signed Lin gap η. The 2PBVP (11), (12), (14) and (15) is well defined and the test function η depends on a single internal parameter, which can be thought of as identifying the end point u a (0) ∈ L a or, alternatively, u r (1) ∈ L r . Here, T a , T r and η are free parameters that vary as the end points of the corresponding orbit segments move along L a and L r . Again, once chosen, the vectors v Z and n Z remain fixed during the continuation of (11), (12), (14) and (15). Continuing this 2PBVP and monitoring η allows us to detect a canard orbit automatically as η = 0. Figure 1 shows the setup of our implementation of Lin's method for detecting canard orbits in system (16) that we present in the next section. Panel (a1) shows S a ε (red surface) and S r ε (blue surface) up to Σ from the line segments L a (red) and L r (blue), respectively. The initial orbit segments u a and u r that define p a and p r and, hence, the Lin space Z, are shown as the highlighted red and blue trajectories, respectively; the Lin space Z is the dark-gray vertical line in Σ (green plane). Note that the curves S a ε and S a ε intersect on both sides of Z in Σ. Panel (a2) shows Figure 1. Lin's method setup for finding canard orbits in system (16) with µ = 8.5. Panel (a1) shows S a ε (red surface) and S r ε (blue surface) computed from L a and L r , respectively, up to Σ. The initial orbit segments u a (red curve) and u r (blue curve) are each other's symmetric counterparts and define the Lin space Z = span{(0, 0, 1)} (dark-gray line) that defines the Lin gap η. Panel (b1) shows the situation when the Lin gap is closed and the canard orbit ξ 1 (orange) is detected. The relevant objects in Σ are shown in panels (a2) and (b2), respectively. the intersection curves S a ε (red) and S r ε (blue) in Σ together with Z, as defined by p a and p r , and the Lin gap η. Panels (b1) and (b2) show what we would like to achieve: the same objects when η = 0 and the canard orbit ξ 1 (orange) is detected; here Z is shown in orange through ξ 1 and the orange orbits on S a ε and S r ε connect as ξ 1 from L a to L r .
Overall, we want to find all of the zeros of the Lin gap, and hence, the canard orbits as the intersections of S a ε and S r ε in Σ. The intersection curves S a ε and S r ε , as well as the slow manifolds S a ε and S r ε , are computed as part of the process.
4. Systematic detection of canard orbits in a normal form of a folded node. As an illustrative example, we consider the normal form vector field for a folded node, introduced by Wechselberger [54], which is given by Here, the variable z is fast, x and y are the slow variables, and the parameter µ ∈ R is such that µ −1 corresponds to the eigenvalue ratio of the folded node singularity of the reduced flow. The normal form (16) does not have the ratio of time scales ε as part of the equations, because it has been obtained via ε-dependent blow-up and rescaling. The remainder terms are O( √ ε); see [54] for details. In spite of the absence of ε in system (16) and in order to be consistent with the standard notation, we denote the attracting and repelling slow manifolds by S a ε and S r ε and their intersection sets in Σ by S a ε and S r ε , respectively. One advantage of using system (16) is that it has an easy representation of the folded node and the relevant objects that are necessary for the slow-fast analysis. The critical manifold of (16) is which has the attracting sheet which implies that S a ε and S r ε are related by this symmetry; this means that one can obtain S r ε by reflecting S a ε with respect to the x-axis in Σ. Hence, the intersection points between S a ε and S r ε in Σ, that is, the canard orbits, occur exactly when S a ε and/or S r ε cross the line z = 0 in Σ. The geometry of the slow manifolds in (16) and how it depends on the eigenvalue ratio µ was studied in [7], which showed how canard orbits organize the number of small-amplitude oscillations near the folded node. Slow manifolds were computed with a 2PBVP and the symmetry was used to detect canard orbits, for different values of µ, as zeros of the z-coordinates of points on S a ε . We now demonstrate how our Lin's method approach is able to detect canard orbits in (16) without monitoring the z-coordinate along S a ε and/or S r ε . Throughout this section we consider system (16) with µ = 8.5, which implies the existence of  Figure 2. Illustration of the Lin's method approach to detect canard orbits for (16) with µ = 8.5 in the section Σ, which is the (x, z)-plane. Shown are the intersection sets S a ε (red curve) and S r ε (blue curve), together with the Lin space Z (vertical dark-gray line). Panels (a1) and (a2) show the detection of the canard orbit ξ 0 (cyan), and panels (b1) and (b2) show the detection of ξ 1 (orange). five canard orbits (two primary and (µ − 1)/2 = 3 secondary canard orbits) [54]. We define L a := S a ∩ {x = −3} and L r := S r ∩ {x = −3}.
Due to the symmetry of system (16), it is natural to choose the initial orbit segments u a and u r defining p a and p r , respectively, as symmetric counterparts. This means that Z is the vertical direction of the (x, z)-plane Σ. Figure 2 illustrates the continuation runs to find the first two canard orbits ξ 0 in row (a) and ξ 1 in row (b), respectively, which are located on either side of the initial choice of Z (thick dark gray line). Shown are the intersection sets S a ε and S r ε in Σ and the points p a and p r defining Z. The light-gray vertical lines are translations of Z at selected points on S a ε (and S r ε ) during the computation, so that we can appreciate how the Lin gap η is changing; the arrows indicate the direction of the end points u a (1) and u r (0) along S a ε and S r ε during the respective continuation run. Figure 2(a1) shows the first continuation run, when the corresponding end points tracing out S a ε and S r ε initially move to the right. Here, we monitor η and detect the canard orbit ξ 0 when η = 0. Panel (a2) illustrates that the end points u a (1) and u r (0) continue past a simultaneous tangency of Z with both S a ε and S r ε at points denoted p t a and p t r , respectively. While the computation continues, no further canard orbits exist past the double tangency. Figure 2(b1) shows the second continuation run, when u a (1) and u r (1) initially move from p a and p r to the left. The Lin section Z becomes again tangent to both S a ε and S r ε simultaneously, at points also denoted p t a and p t r . Panel (b2) shows that the end points on S a ε and S r ε are continued past the double tangency. Monitoring η, as in each continuation run, results in the detection of the canard orbit ξ 1 when η = 0. Note that panels (b1) and (b2) are part of the same continuation run, which continues past further double tangencies and results in the detection of all the subsequent canard orbits ξ 2 -ξ 4 . In between any two consecutive canard orbits, that is, zeros of η, the normal form (16) has a double tangency of Z with both S a ε and S r ε . Figure 3 represents a three-dimensional view of the result of the two continuation runs. Shown are the critical manifold S with the fold curve F , the slow manifolds S a ε and S r ε computed from the line segments L a ⊂ S a and L r ⊂ S r up to Σ, the intersection curves S a ε and S r ε and the canard orbits ξ 0 -ξ 4 . Note that ξ 1 -ξ 3 lie between ξ 0 and ξ 4 on S a ε and S r ε , which correspond to the so-called strong and weak maximal canards, respectively. Furthermore, from ξ 0 to ξ 3 the number of twists of the canard orbits around ξ 4 increases, as predicted by the theory [52,54]. Note that for part of S r ε and S r ε and the respective canard orbits shown in Fig. 3, the flow on L a and L r initially points away from the fold curve F . However, this does not affect their crossing through Σ and the detection of the canards with our approach. 5. Lin space in general position. The choice of Z as orthogonal to the symmetry axis is natural, but nongeneric. This makes the detection of canard orbits in Fig. 2 quite special. Figure 4 shows the setup of Lin's method for a different, more typical Lin space Z, which comes from a different choice of the points p a and p r . Panel (a1) shows the first continuation run, when u a (1) and u r (0) move initially to the right. The Lin section Z becomes tangent to S a ε at the point p t a before the canard orbit ξ 0 is detected. Panel (a2) shows that the continuation goes on and u a (1) keeps tracing S a ε further, but u r (0) turns around and goes back over the same part of S r ε . Even though the tangency point p t a is extremely close to the canard orbit ξ 0 , it is not possible to detect ξ 0 by Lin's method for this choice of Z. Similarly, row (b) of Fig. 4 illustrates the second continuation run; this time the end points on S a ε and S r ε move initially to the left. Panel (b1) shows a tangency of Z with S a ε at the point p t a , and panel (b2) shows that u a (1) continues along S a ε , while u r (0) goes back along S r ε as in row (a), so no canard orbits are found for this choice of Z.
Overall, we conclude that a single tangency of the fixed Lin section Z with one of the curves S a ε or S r ε poses a problem for the detection of canards orbits. The curve that has the tangency is traced past the fold, as was shown. However, as x y z Figure 3. Three-dimensional view of the slow manifolds S a ε and S r ε , and all canard orbits ξ 0 -ξ 4 of the normal form (16) for µ = 8.5.
a result, the continuation of the other curve also enforces a 'fold', meaning that the curve is traced back along the part that was already computed. Because of this, the next canard orbit cannot be found. This artificial turn-around is a result of the coupling between the two intersection curves by the boundary condition (14). This also explains why our approach works well when a double fold occurs; a double tangency of Z with both S a ε and S r ε does not have this issue and both  Figure 4. Illustration of the Lin's method approach for (16) with µ = 8.5 and a Lin space Z in general position. Panel (a1) shows when Z becomes tangent to S a ε as the end points tracing S a ε and S r ε move to the right, and panel (a2) shows that it is not possible to detect a canard orbit by keeping Z fixed. Panels (b1) and (b2) show a similar situation when the end points tracing S a ε and S r ε move to the left.
curves are continued correctly past the double fold. Since there is an alternation of canard orbits and (double) tangencies due to the counter-rotating spiral of S a ε and S r ε , a general algorithm will need to deal with the issue of single tangencies. The underlying idea behind our approach for finding the next and, hence, all of the canard orbits, is to move the Lin direction Z in such a way that the nearby double tangency is found. This is achieved by continuing the detected single tangency while varying Z until the second tangency is also detected. Once the double tangency has been found, Z is again fixed and we resume the computation of finding the next canard orbit. This intermediate step of adjusting Z may seem overly complicated compared to, for example, re-adjusting boundary conditions in a different way.
However, its strength lies in the fact that Z is adjusted only when necessary and fully automatically as part of a single overall computation, without prior knowledge of the geometry of the intersection curves or input from the user. 6. Detecting a double tangency and overall method. Instead of defining a fixed Lin section Z that will become tangent to both S a ε and S r ε simultaneously, which is generally impossible, one can first detect a single tangency of Z with either S a ε or S r ε , which is generic. To this end, we introduce the projections β a = n Z · u a (1), and of the end points u a (1) and u r (0) in Σ onto n Z , respectively. Specifically, the boundary conditions (17) and (18) introduce the parameters β a and β r that represent these projections.
We continue the overall 2PBVP (11), (12), (14) and (15) with the additional boundary conditions (17) and (18), where v Z and n Z are parameterized by α as defined in (13); the parameter α is fixed, while T a , T r , η, β a and β r are free continuation parameters. We monitor β a and β r and detect their folds; namely, a fold of β a or β r indicates that Z has become tangent to S a ε or S r ε , respectively, since it corresponds to an extremum of the projection on n Z . AUTO is able to detect such a fold as a limit point; once it is detected, the fold can be continued. To this end, α in (13) becomes an extra free parameter. Hence, the Lin direction Z rotates while the fold (tangency) that was found is also continued, that is, its existence is preserved. A second fold of β r or β a , respectively, can be detected at, say α = α * . This new fold then corresponds to a simultaneous double fold, where Z is tangent to both S a ε and S r ε . We then fix Z by setting α = α * and resume the detection of canard orbits as zeros of the Lin gap η for this adjusted Lin space.
Our general approach for the systematic detection of all the canard orbits via Lin's method, therefore, consists of the following steps: (I) Initialization. Define a Lin space Z, that is, define vectors v Z and n Z , by computing p a and p r . (II) Main continuation. Continue the overall 2PBVP (11), (12), (14) and (15), with the extra boundary conditions (17) and (18) for fixed α, with T a , T r , η, β a and β r as continuation parameters, while detecting canard orbits as η = 0 and monitoring for fold points in β a and β r . (III) Intermediate step. If a fold of β a , that is, a tangency of Z with S a ε , is detected, continue this fold with α as a continuation parameter until a fold point of β r is detected for α = α * ; and similarly for a first fold of β r . Then fix α = α * and return to step (II). We briefly mention a few technical implementation issues. First of all, when one activates the fold detection in AUTO, the continuation routines possibly detect folds with respect to different parameters. Therefore, it is necessary to be particularly careful to ensure that detected folds indeed concern either β a or β r . Secondly, at the moment of detecting the double tangency, one needs to take care that the end points u a (1) and u r (0) on S a ε and S r ε are continued in the correct direction when resuming step (II). This can be ensured by continuing u a (1) and u r (0) ever so slightly past the double tangency, meaning that the detection of canard orbits with the Lin space given by α = α * is not started right at the tangency points.
We now show that this overall implementation works for a practical example. 7. Systematic detection of canard orbits in the Koper model. We consider the idealized model of a chemical reaction introduced by Koper [34,35]: where x, y and z are real numbers representing chemical concentrations. System (19) was analyzed by Koper to study complex oscillations of chemical systems, where the mixed-mode oscillations are related to the presence of a Shilnikov homoclinic bifurcation [42]. This system has very rich dynamics and different sources for smallamplitude oscillations, such as folded nodes and singular Hopf bifurcations [3,9,20]. We use (19) here as a more typical test-case example of a three-dimensional slow-fast system with two slow variables. To this end, we set λ = 7, k = −10, ε 2 = 1 and ε 1 = 0.01, (20) so that x is fast and y and z are slow. The critical manifold of (19) is which is a cubic surface that has two attracting sheets System (19) with the parameters given by (20) has a folded node at the point that lies on the fold curve F − , where S a,− ε interacts with S r ε . Therefore, we look for canard orbits that arise from the intersections between these two slow manifolds. To simplify the notation, from now on we denote S a,− simply by S a and F − by F . As usual, in order to implement the coupled 2PBVP (11), (12), (14) and (15), we consider the section Σ ⊂ {z = −0.8}, which contains the folded node p 0 and is transverse to both the flow and the critical manifold, together with the lines L a := S a ∩ {x = −1.5} and L r := S r ∩ {x = −0.2}. In the same way as for the general setup described in section 3, we use homotopy steps for finding initial orbits on S a ε and S r ε that define the points p a and p r in Σ and the Lin space Z. The situation here is generic in the sense that there is no symmetry as for the normal form of a folded node singularity. Figure 5 shows a  (19) for the parameters values given by (20).
three-dimensional view of the parts of S a ε and S r ε computed up to Σ. Here L a and L r are not shown, since they are too far from F , outside the frame of the figure, which focuses on the spiraling of S a ε and S r ε near F . The inverse of the eigenvalue ratio λ w /λ s for p 0 is approximately 10.6. Hence, there exist six canard orbits, ξ 0 , . . . , ξ 5 (two primary and (10.6 − 1)/2 = 4 secondary canard orbits) [54]. Figure 6 shows the continuation runs for the detection of the canard orbits ξ 0 , ξ 1 and ξ 2 , in the same layout as in Fig. 2. Panel (a1) shows the first continuation run. Here, the points p 0 a and p 0 r in Σ define the Lin space Z 0 for α = α 0 , and the corresponding end points tracing out S a ε and S r ε move initially to the right. The canard orbit ξ 0 (cyan dot) is detected when η = 0. Then the continuation detects a fold in β a , meaning that the Lin section Z 0 becomes tangent to S a ε at the point p t a ; the intermediate step (III) from section 6 changes the Lin section so that the computation can continue (not shown), but there are no further canard orbits past this point. Figure 6(a2) illustrates the second continuation run, when u a (1) and u r (0) move initially to the left. Here, a fold in β r is detected, where the Lin space Z 0 is tangent to S r ε at the point p t r . This single tangency is very close to the double tangency, yet we know from section 5 that the continuation cannot go past it to find the next canard orbit. The intermediate step (III) detects the double tangency as a fold also on β a for α = α 1 near α 0 ; the points p 1 a and p 1 r of double tangency of the Lin space with S a ε and S r ε , respectively, define the new Lin Figure 7. Intermediate step (III) for the detection of a simultaneous tangency of the Lin space with S a ε and S r ε , for the Koper model (19) with parameters as in (20). Panel (a1) shows the detection of the points defining the Lin space Z 1 and panel (a2) shows the corresponding fold of β a . Panels (b1) and (b2) show step (III) for the detection of the points defining the Lin space Z 2 .
process continues with steps (II) and (III) until all other canard orbits ξ 3 -ξ 5 are also detected.
As illustrated, finding a double tangency of the Lin section with S a ε and S r ε through the intermediate step (III) is the key for our systematic detection of canard orbits. Figure 7 illustrates in rows (a) and (b) the continuation runs of the intermediate step (III) for finding the updated Lin spaces Z 1 and Z 2 , respectively. Panel (a1) shows an enlargement near the double tangency of Z 1 with S a ε and S r ε . The continuation run starts with the Lin space Z 0 (gray line) through the point u a (1) ∈ S a ε and the first detected tangency point u r (0) = p t r ∈ S r ε , which is the situation just after step (II). Since we are continuing a fold in β r , the intermediate end points u i a (1) and u i r (0) move along S a ε and S r ε , respectively, in such way that the line through them, that is, the moving Lin space (light-gray line), is always tangent to S r ε . The continuation then detects a fold in β a as well, that is, the two points p 1 a ∈ S a ε and p 1 r ∈ S r ε of double tangency that define the Lin space Z 1 (dark-gray line). Figure 7(a2) shows the fold in β a during step (III), which is best observed when it is plotted against the (AUTO) L 2 -norm · 2 of the corresponding orbit segment. Similarly, Fig. 7(b1) shows an enlargement near the double tangency of the Lin space Z 2 with S a ε and S r ε , and illustrates step (III) for finding and updating Z 2 when starting from Z 1 . Notice that the angle between the Lin sections Z 1 and Z 2 is now much larger, so that it is easier to see the process. The fold of β r , again best observed when plotted against the (AUTO) L 2 -norm · 2 , is shown in Fig. 7(b2). Figure 8 shows a three-dimensional view of S a ε and S r ε together with the section Σ and all canard orbits ξ 0 -ξ 5 , which were detected systematically as just described. Note the spiraling of the slow manifolds and the different numbers of small oscillations performed by ξ 0 -ξ 4 around ξ 5 ; here, ξ 0 and ξ 5 are again the so-called strong and weak maximal canards, respectively. The situation is qualitatively the same as for the normal form in Fig. 3, as one would expect, but without exact symmetry between S a ε and S r ε . Importantly, ξ 0 -ξ 5 are nevertheless found systematically via the alternation of steps (II) and (III).
8. Discussion. We presented a Lin's method approach for the detection of all the canard orbits arising from a folded node in slow-fast systems in R 3 with two slow and one fast variables, based on the continuation of solutions to a suitable 2PBVP. We showed a new context where the underlying idea of Lin's method is utilized to achieve an efficient and accurate numerical implementation. It employs two independent well-posed 2PBVPs that are coupled through additional boundary conditions.
Our overall implementation allows for the systematic detection of all of the canard orbits arising as intersection points of an attracting and a repelling slow manifold in a two-dimensional section Σ that contains a folded node. It has the novel feature that the Lin space is updated in an suitable way during an intermediate continuation step, which overcomes a geometric obstruction that arises due to the counter-rotating nature of slow manifolds near a folded node. The spiraling nature of the slow manifolds is reflected by their intersection curves in Σ. Our method uses the geometry of these intersection curves, in particular the alternation between their intersection points and tangencies with the Lin space. Specifically, we vary the Lin space until we detect double tangencies with both slow manifolds, and then fix the Lin space again. This approach has been demonstrated for the automatic detection and computation of canards in the Koper model.
Once a canard orbit is detected, it is natural to continue it in system parameters to study how the interaction between the attracting and repelling slow manifolds changes. This is particularly useful, for example, to identify the creation or destruction of regimes with different numbers of small-amplitude oscillations. To this end, one can fix the Lin gap at η = 0 in the overall 2PBVP, as was done in [45]. This has the advantage that one can still keep track of the intersection of the canard orbits with the section Σ. Alternatively, one can concatenate the two orbit segments that form the canard orbit, so that it becomes a single orbit segment that is then continued in parameters.
The overall 2PBVP implementation of Lin's method presented in this paper may also be useful for the detection of connecting orbits in a broader context, especially when the intersection curves of the manifolds of interest generate folds with any Lin space, for example, due to their spiraling.