{'fraction_non_alphanumeric': 0.05864033219761499, 'fraction_numerical': 0.027230621805792162, 'mean_word_length': 3.925144205558469, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discretized Fast-Slow Systems with Canards in Two Dimensions

We study the problem of preservation of maximal canards for time discretized fast-slow systems with canard fold points. In order to ensure such preservation, certain favorable structure preserving properties of the discretization scheme are required. Conventional schemes do not possess such properties. We perform a detailed analysis for an unconventional discretization scheme due to Kahan. The analysis uses the blow-up method to deal with the loss of normal hyperbolicity at the canard point. We show that the structure preserving properties of the Kahan discretization for quadratic vector fields imply a similar result as in continuous time, guaranteeing the occurrence of maximal canards between attracting and repelling slow manifolds upon variation of a bifurcation parameter. The proof is based on a Melnikov computation along an invariant separating curve, which organizes the dynamics of the map similarly to the ODE problem.


Introduction
In this paper, we study the effect of the time discretization upon systems of ordinary differential equations (ODEs) which exhibit the phenomenon called "canards". It takes place, under certain conditions, in singularly perturbed (slow-fast) systems exhibiting fold points. The simplest form of such a system is x = f (x, y, λ, ε), y = εg(x, y, λ, ε), (1.1) where we interpret ε > 0 as a small time scale parameter, separating between the fast variable x and the slow variable y. For λ = 0, the origin is assumed to be a non-hyperbolic fold point, Based on observations of this paper, the employment of Kahan's method for a treatment of canards can also be found in [17]. There, the simplest canonical form for folded, pitchfork and transcritical canards is studied and the focus lies on the linearization along trajectories. While it is demonstrated that explicit Runge-Kutta methods cannot provide symmetry of entry-exit relations, the linearization along the Kahan scheme and similar symmetric, A-stable methods are shown to preserve the typical continuous-time behaviour. Hence, the discussion of symmetry and linear stability in [17] supplements the paper at hand; here, we establish the existence and extension of maximal canards along parameter combinations for the nonlinear problem of folded canards with additional quadratic perturbation terms, in particular using the blow-up technique.
The paper is organized as follows. Section 2 recalls the setting of fast-slow systems in continuous time and summarizes the main result on maximal canards, Theorem 2.2, with a short sketch of the proof, as given in [30]. In Section 3, we study the problem of a maximal canard for systems with folds in discrete time. We establish the Kahan discretization of the canard problem in Section 3.1 and discuss the reduced subsystem of the slow time scale in Section 3.2. In Section 3.3, we introduce the blow-up transformation for the discretized problem. We discuss the dynamics for the entering and exiting chart in Section 3.4, and for the rescaling chart in Section 3.5. In Section 3.6, we explore the dynamical properties of the Kahan map in the rescaling chart, including a formal conserved quantity, an invariant measure and an invariant separating curve. Following this, we conduct the Melnikov computation along the invariant curve in Section 3.7, leading to the proof of the main Theorem 3.11, which is the discrete-time analogue to Theorem 2.2. Finally, we provide various numerical illustrations in Section 3.8 and conclude with an outlook in Section 4.
Thus, we succeeded in adding the problem of maximal canards to the recent list of results, where a geometric analysis shows that certain features of fast-slow systems with non-hyperbolic singularities can be preserved via a suitable discretization, including the cases of the fold singularity [39], the transcritical singularity [18] and the pitchfork singularity [1]. More broadly viewed, our results also provide a continuation of a line of research on discrete-time fast-slow dynamical systems, which includes the study of canard/delay behavior in iterated maps via normal form transformations [38], non-standard analysis [20,21], renormalization [2], Gevrey series [3], complex-analytic methods [22], and phase plane partitioning [37].

Fast-slow systems
We start with a brief review and notation for continuous-time fast-slow systems. Consider a system of singularly perturbed ordinary differential equations (ODEs) of the form ε dx dτ = εẋ = f (x, y, ε), dy dτ =ẏ = g(x, y, ε), x ∈ R m , y ∈ R n , 0 < ε 1 , where f, g, are C k -functions with k ≥ 3. Since ε is a small parameter, the variables The reduced problem (2.3) can be understood as a dynamical system on the critical manifold S 0 = {(x, y) ∈ R m+n : f (x, y, 0) = 0} .
Observe that the manifold S 0 consists of equilibria of the layer problem (2.4). S 0 is called normally hyperbolic if for all p ∈ S 0 the matrix D x f (p) ∈ R m×m has no eigenvalues on the imaginary axis. For a normally hyperbolic S 0 , Fenichel theory [19,28,35,44] implies that, for sufficiently small ε, there is a locally invariant slow manifold S ε such that the restriction of (2.1) to S ε is a regular perturbation of the reduced problem (2.3). Furthermore, it follows from Fenichel's perturbation results that S ε possesses an invariant stable and unstable foliation, where the dynamics behave as a small perturbation of the layer problem (2.4).

Main result on maximal canards in slow-fast systems with a fold
A challenging phenomenon is the breakdown of normal hyperbolicity of S 0 such that Fenichel theory cannot be applied. Typical examples of such a breakdown are found at bifurcation points p ∈ S 0 , where the Jacobi matrix D x f (p) has at least one eigenvalue with zero real part. The simplest examples are folds in planar systems (m = n = 1), i.e., points p = (x 0 , y 0 ) ∈ R 2 (without loss of generality p = (x 0 , y 0 ) = (0, 0)) where ∂f /∂x vanishes and in whose neighbourhood S 0 looks like a parabola. The left part of S 0 (with x < 0) is denoted by S a (a for "attractive"), while its right part (with x > 0) is denoted by S r (r for "repelling"). These notations refer to the properties of dynamics of the layer problem in the region y > 0 (see e.g. [35, Figure 8.1]). By standard Fenichel theory, for sufficiently small ε > 0, outside of an arbitrarily small neighborhood of p, the manifolds S a and S r perturb smoothly to invariant manifolds S a,ε and S r,ε .
In the following we focus on the particularly challenging problem of fold points admitting maximal canards. In this case, the critical curve S 0 = {f (x, y, 0) = 0} can be locally parametrized as y = ϕ(x) such that the reduced dynamics on S 0 are given bẏ In our setting, the function at the right-hand side is smooth at the origin, so that the reduced flow goes through the origin via a maximal solution x 0 (t) of (2.5) with x 0 (0) = 0. The solution (x 0 (t), y 0 (t)) with y 0 (t) = ϕ(x 0 (t)) connects both parts S a and S r of S 0 . However, there is no reason to expect that for ε > 0, the (extension of the) solution parametrizing S a,ε will coincide with the (extension of the) solution parametrizing S r,ε , unless there are some special reasons, like symmetry, forcing such a coincidence.
Definition 2.1. We say that a planar slow-fast system admits a maximal canard, if the extension of the attracting slow manifold S a,ε coincides with the extension of a repelling slow manifold S r,ε .
Example. Consider the system corresponding to f (x, y, ε) = x 2 − y and g(x, y, ε) = x. For the reduced system (ε = 0) we obtain y = ϕ(x) = x 2 and 2xẋ = x, henceẋ = 1/2 (regular at x = 0). The solution x 0 (t) is given by Observe that the system is symmetric with respect to the reversion of time τ → −τ simultaneously with x → −x. This ensures the existence of the maximal canard also for any ε > 0. In this particular example, one can easily find the maximal canard explicitly. Indeed, one can easily check that, for any ε > 0, is a solution of (2.6) which parametrizes the invariant set which consists precisely of the attracting branch S a,ε = {(x, y) ∈ S ε : x < 0} and the repelling branch S r,ε = {(x, y) ∈ S ε : x > 0}, such that trajectories on S ε go through x = 0 with the speeḋ x = ε/2. However, any generic perturbation of this example, e.g. with g(x, y, ε) = x + x 2 , will destroy its peculiarity and will not display a maximal canard.
Thus, maximal canards are not a generic phenomenon in the above setting. In order to find a context where they become generic, we have to consider families depending on an additional parameter λ: (2.8) We assume that at λ = ε = 0, the vector fields f and g satisfy the above conditions. By a local change of coordinates, the problem can be brought into the canonical form The main result on existence of maximal canards, as given in [30,Theorem 3.1], can be summarized as follows. Set  (2.9) such that the solution (x 0 (t), y 0 (t)) of the reduced problem for ε = 0, λ = 0 connects S a and S r . Assume that C = 0. Then there exist ε 0 > 0 and a smooth function The main result of this paper will be a discretized version of Theorem 2.2 restricted to quadratic vector fields, proving for Kahan maps the extension of canards along a parameter curve, as opposed to [17] where only Example (2.6) and its linearization are studied.
The proof of Theorem 2.2 is based on the blow-up technique, transforming the singular problem to a manifold where the dynamics can be desingularized and studied in two different charts. The crucial step in the second chart K 2 is the continuation of center manifold connections via a Melnikov method based on an integral of motion H. In the Appendix A, we summarize this procedure from [30], adding several observations on the dynamics; its separatrix, its invariant measure and an alternative non-Hamiltonian expression that relates to the discrete-time proof we will provide in the following.
3 Maximal canard for a system with a fold in discrete time

Kahan discretization of canard problem
We discretize system (2.9) with the Kahan method. It was introduced in [29] as an unconventional discretization scheme applicable to arbitrary ODEs with quadratic vector fields. It was demonstrated in [40,41,42] and in [5] that this scheme tends to preserve integrals of motion and invariant volume forms. There are few general results available to support this claim, in particular, two general cases of preservation of invariant volume forms in [41, Section 2] and a similar result for Hamiltonian systems with a cubic Hamilton function in [5]. However, the number of particular results not covered by any general theory and reviewed in the above references, is quite impressive. Our study here will contribute an additional evidence, as the result of Section 3.6.2 also belongs to this category, i.e., is not covered by known general statements.
Consider an ODE with a quadratic vector field: where each component of Q : R n → R n is a quadratic form, B ∈ R n×n and c ∈ R n . The Kahan discretization of this system reads as is the symmetric bilinear form such thatQ(z, z) = Q(z). Note that equation (3.2) is linear with respect toz and therefore defines a rational mapz = F f (z, h), which approximates the time h shift along the solutions of the ODE (3.1). Further note that F −1 f (z, h) = F f (z, −h) and, hence, the map is birational. An explicit form of the map F f defined by equation (3.2) is given bỹ In order to be able to apply the Kahan discretization scheme, we restrict ourselves to systems (2.1), (2.2) which are quadratic, that is, to resp.
This opens the way of extending our present results for more general (not necessarily quadratic) systems (2.9). In the present paper we restrict ourselves to the case (3.5), since the algebraic structure keeps the calculations clear and explicit and demonstrates the central methodological aspects of our proofs. However, we additionally apply the scheme (3.6) to the folded canard problem with cubic nonlinearity in Section 3.8, illustrating its numerical capacity beyond the quadratic case. A proof of maximal canards for the non-quadratic case remains an open problem for future work.

Reduced subsystem of the slow flow
Kahan discretization of (3.4) reads: (3.7) Proposition 3.2. The reduced system (3.7) with ε = 0 defines an evolution on a curve which supports a one-parameter family of solutions x h (n; x 0 ) with x h (0; x 0 ) = x 0 . For small ε > 0, this curve is perturbed to normally hyperbolic invariant curves S a,h,ε resp. S r,h,ε of the slow flow (3.8) for x < 0, resp. for x > 0.
For the simplest case a 1 = a 2 = a 4 = a 5 = 0 and λ = 0, everything can be done explicitly. Straightforward computations lead to the following results.
The reduced system has an invariant critical curve The evolution on this curve is given byx = x + h 2 , so that x h (n; x 0 ) = x 0 + nh 2 . For the full system (3.8), the symmetry x → −x, h → −h ensures the existence of an invariant curve whose parts with x < 0, resp x > 0 are the invariant curves S a,h,ε resp. S r,h,ε . This curve supports solutions with x(n) = x 0 + nh 2 . Thus, system (3.8) exhibits a maximal canard. Our goal is to establish the existence of a maximal canard for system (3.7).

Blow-up of the fast flow
Kahan discretization of the fast flow (3.5) is the system (3.7) with h → hε: (3.12) We introduce a quasi-homogeneous blow-up transformation for the discrete time system, interpreting the step size h as a variable in the full system. Similarly to the continuous time situation, the transformation reads The change of variables in h is chosen such that the map is desingularized in the relevant charts.
This transformation is a map Φ : B → R 5 . If F denotes the map obtained from the timediscretization, the map Φ induces a map F on B by Φ • F • Φ −1 = F . Analogously to the continuous time case, we are using the charts K i , i = 1, 2, to describe the dynamics. The chart K 1 (settingȳ = 1) focuses on the entry and exit of trajectories, and is given by In the scaling chart K 2 (settingε = 1) the dynamics arbitrarily close to the origin are analyzed. It is given via the mapping 14) The change of coordinates from K 1 to K 2 is denoted by κ 12 and, for ε 1 > 0, is given by Similarly, for y > 0, the map κ 21 = κ −1 12 is given by

Dynamics in the entering and exiting chart K 1
Here we extend the dynamical equations (3.12) bỹ (3.17) and then introduce the coordinate chart K 1 by (3.13): defined on the domain where ρ, δ, ν > 0 are sufficiently small. To transform the map (3.12) into the coordinates of K 1 , we start with the particular case a 1 = a 2 = a 4 = a 5 = 0, generated by difference equations supplied, as usual, by (3.17). Written explicitly, this is the map Upon substitution K 1 , we have: where Setting we come to the following expression for the map (3.21) in the chart K 1 : Now it is straightforward to extend these results to the general case of the map (3.12) with arbitrary constants a i . For this, we observe: -in the first equation, the terms y and x 2 on the right-hand side scale as r 2 1 and r 2 1 x 2 1 , while the terms εx and xy scale as r 3 1 ε 1 x 1 and r 3 1 x 1 , respectively; -in the second equation, the terms εx and ελ on the right-hand side scale as r 3 1 ε 1 x 1 and r 3 1 ε 1 λ 1 , while the terms εy and εx 2 scale as r 4 1 ε 1 and r 4 1 ε 1 x 2 1 , respectively.
Therefore, we can treat all terms involving a 1 , a 2 , a 4 , a 5 as O(r 1 ). The resulting map is given by formulas analogous to (3.33), with We now analyze the dynamics of this map.
We have: In all other directions, the multipliers of these fixed points are equal to 1.
• Similarly, we have on {ε 1 = 0, λ 1 = 0} ∩ D 1 for small r 1 > 0: By the implicit function theorem, we can conclude that on {ε 1 = 0, λ 1 = 0} ∩ D 1 , there exist two families of normally hyperbolic (for h 1 > 0) curves of fixed points denoted as S a,1 (h 1 ) and S r,1 (h 1 ), parametrized by r 1 ∈ [0, ρ] and ending for r 1 = 0 at p a,1 (h 1 ) and p r,1 (h 1 ), respectively. For the map (3.21), corresponding to difference equation (3.20) (that is, to (3.12) with all a i = 0), the O(r 1 )-term vanishes, and the above families are simply given by • On the invariant set {r 1 = 0, λ 1 = 0} ∩ D 1 , the dynamics of x 1 , ε 1 and h 1 are given bỹ We compute the Jacobi matrices of the map (3.34) at p a,1 (h 1 ) and p r,1 (h 1 ), restricting to the invariant set The matrix A a has a two-dimensional invariant space corresponding to the eigenvalue 1, spanned by the vectors v (1) a = (0, 0, 1) and v (2) a = (−1, 4, 0) , such that Similarly, the matrix A r has a two-dimensional invariant space corresponding to the eigenvalue 1, spanned by the vectors v (1) r = (0, 0, 1) and v (2) r = (1, 4, 0) , such that It is instructive to compare this with the continuous-time case h 1 → 0 (see, e.g., [ and v (2) a are eigenvectors of the corresponding linearized system, with v (1) a being tangent to S a, 1 and v (2) a corresponding to the center direction in the invariant plane r 1 = 0 (and similarly for v (1) r and v (2) r ). We summarize these observations into the following statement. Proposition 3.3. For system (3.33), there exist a center-stable manifold M a,1 and a centerunstable manifold M r,1 , with the following properties: 1. For i = a, r, the manifold M i,1 contains the curve of fixed points S i,1 (h 1 ) on {ε 1 = 0, λ 1 = 0} ⊂ D 1 , parametrized by r 1 , and the center manifold N i,1 whose branch for ε 1 , h 1 > 0 is unique (see Figure 3 (b)). In D 1 , the manifold M i,1 is given as a graph x 1 =ĝ i (r 1 , ε 1 , λ 1 , h 1 ).

2.
For i = a, r, there exist two-dimensional invariant manifolds M i,1 which are given as graphs Proof. The first part follows by standard center manifold theory (see, e.g., [26]). There exist two-dimensional center manifolds N a,1 and N r,1 , parametrized by h 1 , ε 1 , which at ε 1 = 0 coincide with the sets of fixed points respectively (see Figure 3 (b)). Note that, by (3.34), on Hence, for δ small enough, the branch of the manifold N a,1 on On the other hand, we observe that for Thus, for x 1 from a neighborhood of 1, we see that ν < 2K 1+K 2 < 1 guarantees that, for δ small enough depending on K, the branch of the manifold N r,1 on The second part follows from the invariancesr 1λ1 = r 1 λ 1 andh 1 /r 1 = h 1 /r 1 , compare [18, Proposition 3.3 and Figure 2] for details.

Dynamics in the scaling chart K 2
Next, we investigate the dynamics in the scaling chart K 2 , in order to find a trajectory connecting M a,1 with M r,1 , or M a,1 with M r,1 respectively. Recall from (3.14) that in chart K 2 we have In this chart and upon the time rescaling t = t 2 /r 2 , equation (3.5) takes the form where the prime now denotes the derivative with respect to t 2 , compare (A.4). Since in this chart r 2 = √ ε is not a dynamical variable (remains fixed in time), we will not write down explicitly differential, resp. difference evolution equations for λ 2 = λ/ √ ε and for h 2 = h √ ε. We will restore these variables as we come to the matching with the chart K 1 . The Kahan discretization of equation (3.37) with the time step h 2 can be written as On the blow-up manifold r 2 = 0, we are dealing with the simple model system This yields the birational map This gives the following expressions for the map F = (F 1 , F 2 ) andĴ = (Ĵ 1 ,Ĵ 2 ) in (3.38): Explicit expressions for the functionsĜ 1 andĜ 2 can be easily obtained, as well, but are omitted here due to their length.

Dynamical properties of the model map in the scaling chart
For a better readability, we omit index "2" referring to the chart K 2 starting from here. In particular, we write x, y, r, λ, h for x 2 , y 2 , r 2 , λ 2 , h 2 rather than for the original variables (before rescaling). Similarly to the continuous-time case, we start the analysis in K 2 with the case λ = 0, r = 0 for h > 0 fixed. This means that we study the dynamics of the map given by F (3.41), F : which comes as the solution of the difference equation We discuss in detail the most important properties of the model map (3.43).

Formal integral of motion
Recall that, for r = λ = 0, the ODE system (A.6) in the chart K 2 has a conserved quantity (A.7). Its level set H(x, y) = 0 supports the special canard solution (A.11), In general, Kahan discretization has a distinguished property of possessing a conserved quantity for unusually numerous instances of quadratic vector fields. For (A.6), it turns out to possess a formal conserved quantity in the form of an asymptotic power series in h. However, there are indications that this power series is divergent, so that map F (3.43) does not possess a true integral of motion. Nevertheless, it possesses all nice properties of symplectic or Poisson integrators, in particular, a truncated formal integral is very well preserved on very long intervals of time. Moreover, as we will now demonstrate, the zero level set of the formal conserved quantity supports the special family of solutions of the discrete time system crucial for our main results.
We recall a method for constructing a formal conserved quantitȳ To compute the Taylor expansion of the left hand side, we observe: Here, the h and the h 2 terms vanish, as follows from (3.46) and its Lie derivative: Plugging this, as well as a Taylor expansion of H 2 (z) similar to H(z), into (3.47), we see that vanishing of the h 3 terms is equivalent to This is a linear PDE defining H 2 up to an additive term which is an arbitrary function of H.
Following terms H 4 , H 6 , . . . can be determined in a similar manner, from linear PDEs like (3.50) with recursively determined functions on the right hand side.
We now apply this scheme to obtain (the first terms of) the formal conserved quantitȳ  H(x, y, h) for (3.41). It turns out to be possible to find it in the form where  .49) gives The differential equation (3.50) reads in the present case: A solution forH 2 which is a polynomial of degree 4 reads: Hence, we obtain the approximation A straightforward computation shows that on the curve y −x 2 + 1 2 = 0 (the level set H(x, y) = 0), the functionH 2 (x, y) takes a constant value 1 8 . Therefore, the level setH(x, y, h) = 0 is given, Remarkably, we have the following statement.

Proposition 3.4. The curve (3.56) represents a zero level set of the (divergent) formal integral H(x, y, h). More precisely, on this curve
We will not prove this statement, but rather derive a different dynamical characterization of the curve (3.56).

Invariant measure Proposition 3.5. The map F given by (3.43) admits an invariant measure
with ϕ h (x, y) given in (3.56). This measure µ h is singular on the curve ϕ h (x, y) = 0.
Proof. Difference equations (3.44) can be written as a linear system for (x,ỹ): Differentiating with respect to x, y, we obtain: Computing determinants, we find: (3.58) Next, we derive from the first equation in (3.43): Since the system (3.44) is symmetric with respect to interchanging (x, y) ↔ (x,ỹ) with the simultaneous change h → −h, we can perform this operation in the latter equation, resulting in Comparing the last two formulas, we obtain: Together with (3.58), this results in which is equivalent to the statement of proposition.

Invariant separating curve
It turns out that the singular curve of the invariant measure µ h is an invariant curve under the map (3.43).

Proposition 3.6. The parabola
is invariant under the map F given by (3.43). Solutions on S h are given by For (x, y) ∈ S h , we have: , we obtain upon a straightforward computation: This proves the first two claims. As for the last claim, we compute by differentiating the first equation in (3.43): For (x, y) ∈ S h , this gives: , which implies inequalities (3.63). (We remark that the right hand side tends to infinity as  We can show the following connection to the chart K 1 : Lemma 3.7. The trajectory γ h (n), transformed into the chart K 1 via for large |n|, lies in M a,1 as well as in M r,1 .
Proof. From (3.16) there follows that for sufficiently large |n|, the component ε 1 (n) of γ 1 h (n) is sufficiently small such that γ 1 h , which lies on the invariant manifold κ 21 (S h , h), has to be in N a,1 for n < 0, and in N r,1 for n > 0 respectively, due to the uniqueness of the invariant center manifolds (see Proposition 3.3). In particular, observe that, if h is small enough, γ 1 h reaches an arbitrarily close vicinity of some p a,1 (h * 1 ) for sufficiently large n < 0 and of some p r,1 (h * 1 ) for sufficiently large n > 0, within N a,1 ⊂ M a,1 and N r,1 ⊂ M r,1 respectively (see also Figure 3 (b)). This finishes the proof.
The trajectory γ h is shown in global blow-up coordinates as γh in Figure 3 (a), in comparison to the ODE trajectoryγ 0 corresponding to γ 0,2 in K 2 .  0 (b). The figures also show the special ODE solutionγ 0 connectingp r (h) andp a (h) (a), and γ 0,1 connecting p r,1 (h * 1 ) and p a,1 (h * 1 ) for fixed h * 1 > 0 (b) respectively. In Figure (a), the fixed pointsq in (h) andq out (h), forε = 0, are added, whose existence can be seen in an extra chart (similarly to [30]). In Figure (b), the trajectory γ 1 h is shown on the attracting center manifold N a,1 ⊂ M a,1 and on the repelling center manifold N r,1 ⊂ M r,1 (see Section 3.4 and Lemma 3.7).

Melnikov computation along the invariant curve
We consider a Melnikov-type computation for the distance between invariant manifolds, which is a discrete time analogue of continuous time results in [31] and, for a more general framework, in [43].
Consider an invertible map depending on a parameter µ: where (x, y) ∈ R 2 , and F = (F 1 , F 2 ) , G = (G 1 , G 2 ) are C k , vector-valued maps, k ≥ 1. The following theory can be easily extended to µ ∈ R m , like in [43], but for reasons of clarity we formulate it for µ ∈ R.
We formulate the following assumptions: (A1) There exist invariant center manifolds M ± of the dynamical system (3.65), given as graphs of C k -functions y = g ± (x, µ) and intersecting at µ = 0 along the smooth curve where g ± (x, 0) = g(x, 0).
The following proposition is a discrete time analogue of [31, Proposition 3.1].
We now apply Proposition 3.8 (or, better to say, its generalization for the case of two parameters µ = (r 2 , λ 2 )) to the Kahan map (3.38) in the rescaling chart K 2 . First of all, we have to justify Assumptions (A1)-(A4) for this case. Assumption (A1) follows from the fact that for µ = (r, λ) = 0, the center manifolds M a,2 and M r,2 intersect along the curve S h given in (3.61). Assumption (A2) follows from the explicit formula (3.62) for the solution γ h,x 0 , as well as from formulas (3.42) for the functionsĴ and similar formulas for the functionsĜ. Assumption (A3) follows from the existence of the center manifolds away from µ = (r, λ) = 0, established in Proposition 3.3. Turning to the assumption (A4), we have the following results. Proposition 3.9. For problem (3.38), the adjoint linear system (3.67), (3.69) has the decaying solution (a(−k)) −1 for n > 0, (3.71) and (3.72) We have: Here the symbol ≈ relates quantities whose quotient has a limit as n → ±∞.
be a fundamental matrix solution of the linear difference equation with det Φ(0) = 1. The first column of the fundamental matrix solution Φ(n) can be found as ∂ x 0 γ h,x 0 . Using formula (3.62) for γ h,x 0 , we have: A fundamental solution of the adjoint difference equation is given by Its second column is a solution of the adjoint system as given in (3.70), with X(n) = det Φ(n).
To compute X(n), we observe that from and from det Φ(0) = 1, there follows a discrete analogue of Liouville's formula: for n > 0, which coincides with (3.71) with a(k) = det A(k). Expression (3.72) for these quantities follows from (3.58).
To prove the estimate (3.73), we observe: Therefore, for n > 0, Using the formula Γ(n + c) ∼ n c Γ(n) by n → +∞ (in the sense that the quotient of the both expressions tends to 1), we obtain for n → +∞: This completes the proof.
With the help of estimates of Proposition 3.9, we derive from Proposition 3.8 the following statement: where O(2) denotes terms of order ≥ 2 with respect to λ, r, and In particular, convergence of the series in equation (3.76) is obtained for any h > 0 and convergence of the series in equation (3.77) is obtained for 0 < h < 4/3.
Proof. The form of the first order separation follows from Proposition 3.8. Furthermore, recall from equation (3.42) that Using Proposition 3.9, this yields (3.76) for any h > 0. Note from equation (3.5) that the highest order n κ we can obtain in the termsĜ(γ h (n), h) is κ = 3 (coming from the term with factor a 2 ) such that for large |n| we have This means that the convergence in (3.77) is given for −4/h 2 + 2 < −1 such that the claim follows.
We are now prepared to show our main result.

The function λ h c has the expansion
where C is given as in (2.12) (for a 3 = 0).
Proof. First, we will work in chart K 2 and show that the quantities d h,x 0 ,λ , d h,x 0 ,r in (3.76), (3.77) with x 0 = 0 approximate the quantities d λ , d r in (A.13), (A.14)(up to change of sign). We prove: where, recall, the functionĴ is defined as in (3.42), and similar formulas hold true also for the functionĜ. Further recall that the Melnikov integrals can be solved explicitly, yielding where a i and C are as introduced in Section 2.2 (for a 3 = 0, see (3.4) and (3.5)). We show (3.78) -the simpler case (3.79) then follows similarly. We observe: We denote by n * the closest integer to α = 2/h 2 + 1/2, and recall that β = 2/h 2 + 3/2. Since we can write, for all n ≥ 2/h 2 + 3/2, Since, with Proposition 3.9 the summands ofŜ(N ) converge to zero even faster for smaller h, we obtain by choosing N ≥ 2/h 2 + 3/2 , and hence T ≥ 2/h + 5h/2, that 3. For T = 3/h, we get by the standard methods the estimate Hence, we can conclude that equations (3.78) and (3.79) hold, and, in particular, that d h,0,λ and d h,0,r are bounded away from zero for sufficiently small h. Recall from (3.75) that where D h,0 (0, 0) = 0. Hence, the fact that d h,0,λ and d h,r are not zero implies, by the implicit function theorem, that there is a smooth function λ h (r) such that D h,0 (r, λ h (r)) = 0 in a small neighborhood of (0, 0). Transforming back from K 2 into original coordinates then proves the first claim. Furthermore, we obtain

Transformation into original coordinates gives
Hence, the second claim follows.
Numerical computations show that h 0 in Theorem 3.11 does not have to be extremely small but that our results are quite robust for different step sizes. In Figure 4, we display such computations for the case a 1 = 1, a 2 = a 4 = a 5 = 0. In this case, the rescaled Kahan discretization in chart K 2 is given bỹ (3.82) Hence, we obtain (3.83) For different values of h and N we calculate and, for the situation of (3.82) withĜ as in (3.83), We compare these quantities with the values of the respective continuous-time integrals d λ = − √ 2π and d r = − √ 2π/2 (we have C = 1/2 in this case). We observe in Figure 4 that the sums converge very fast for relatively small hN in both cases. Additionally, we see that

Numerical illustrations for ε > 0
We illustrate the results by some additional numerics for ε > 0, supplemenenting the illustrations of the dynamics in the rescaling chart K 2 , as given by Figures 1 and 2. Firstly, we consider the simplest case where a i = 0 for all i, i.e., situation (3.8) with invariant curve S ε,h (3.11). Figure 5 shows different trajectories of the map (3.8) for ε = 0.1 and h = 0.02, illustrating the organization of dynamics around S ε,h analogously to the dynamics of (3.44) around S h (3.61) (see Figure 1). Secondly, we consider the map (3.12) with a 1 = 1, i.e., a small additional perturbation of the canonical form, similarly to the end of the previous section. We take ε = 0.1, h = 0.02 and λ = −(a 1 /2)ε, as a leading order approximation of λ h c ( √ ε) (see Theorem 3.11). In Figure 6, we observe that the numerics given by the Kahan discretization approximate very well the maximal canard, which slightly deviates from S ε,h , again illustrating the organization of dynamics into bounded and unbounded trajectories sperated by the maximal canard. Note that we have chosen ε = 0.1 to demonstrate the extension up to a relatively large ε.
In addition, we consider a model with cubic nonlinearity in order to demonstrate the application of the Kahan method beyond the purely quadratic case. Consider the equation  (3.85) such that the cubic nonlinearity does not vanish and we do not directly obtain an explicit form. However, we can use (3.85) as a numerical scheme by always taking the unique real solutionx, closest to x in absolut value, of the cubic polynomial. 2). Observe that the numerics indicate the existence of a maximal canard, also in this situation, separating bounded, now spiralling, orbits and unbounded orbits. Note that the implementation is based on the fact that the cubic polynomial inx always has exactly one real solution, which we take as the next value, plus a complex conjugate pair with non-trivial imaginary part. A more general, algebraic analysis extends beyond the scope of this work and is left for additional research.  Note that the results on maximal canards for perturbations of the canonical form are local and do not make statements on the global stability. The preservation of canards for the van der Pol equation as depicted in Figure 7 is apparently also of predominantly local nature. Hence, we take a closer look in Figure 8, zooming into a neighbourhood of the inward-spiralling orbits from Figure 7. Here, we observe that the Kahan discretization even seems to capture the occurrence of a Hopf bifurcation in a neighbourhood of the maximal canard, as we slightly vary the parameter λ. Furthermore, the scheme seems to avoid crossing trajectories near the fold, which do occur as spurious solutions for some forward numerical methods near maximal canards. Indeed, there are also robust methods from boundary value problems (BVPs) [10,23] and control theory [16,27] to track canards for the van der Pol equation. However, these approaches do not take direct advantage of the polynomial structure, nor of the particular locally approximately integrable or symmetry structures of the van der Pol equation. Hence, building on the presented insights for the Kahan method, we consider an analytical treatment of the discretized cubic canard problem an intriguing direction for future work.

Conclusion
Our results show the importance of combining geometric invariants or integrable structures hidden in blow-up coordinates with suitable discretization schemes. Although we have just treated a very low-dimensional fast-slow fold case, one anticipates similar results also to be relevant for various other higher-dimensional singularities and bifurcation points, where blow-up is a standard tool. For example, it is well-known that in the Bogdanov-Takens unfolding one obtains small homoclinic orbits via a hidden integrable structure visible only after re-scaling. A thorough discretization analysis of higher-dimensional canards, similarly to the one at hand, would also deserve further investigation.
From a numerical perspective, forward integration schemes often provide an exploratory perspective to actually detect interesting dynamics or find a suitable invariant solution for fixed parameter values. In several cases these particular forward solutions are then used as starting conditions in numerical continuation techniques [11,12] to study parametric dependence in a setting of BVPs. BVPs have also been successfully adapted to parametrically continue canardtype solutions [9,10,24,33]. In particular, BVPs for canards turn out to be well-posed with a small numerical error, yet to set up the problem purely by continuation one already needs a very good understanding of phase space for the initial canard orbits. Therefore, a direct numerical integration scheme can be very helpful to automatically yield suitable starting solutions close to a maximal canard.
The Kahan method has mainly turned out to be favorable, since explicit, for quadratic vector fields; hence, in our analysis we have focused on this situation. However, we have seen in the numerical investigations in Section 3.8 that, by using its implicit form, also non-quadratic problems can be tackled, at least numerically. A further investigation into the dynamical and algebraic properties of the scheme, in particular for cubic nonlinearities, is a highly intriguing research question for the future, in general, and also in particular with respect to geometric multiscale problems as the one presented in this work.

Appendix A Existence of maximal canards for ODEs
In order to use specific geometric methods in singular perturbation theory, we consider ε and λ as variables, writing equation (2.8) as Note that in equation (2.9) the Jacobi matrix of the vector field in (x, y, λ, ε) has a quadruple zero eigenvalue at the origin. A well established way to gain (partial) hyperbolicity at such a singularity is the blow-up technique which replaces the singularity by a manifold on which the dynamics can be desingularized. An important technical assumption for this technique is quasi-homogeneity of the vector field f : R n → R n of the ODE (cf. [35,Definition 7.3.2]), which means that there are (a 1 , . . . , a n ) ∈ N n and k ∈ N such that for every r ∈ R and each component f j : R n → R of f we have f j (r a 1 z 1 , . . . , r an z n ) = r k+a j f j (z 1 , . . . , z n ).