Interpolating vector fields for near identity maps and averaging

For a smooth near identity map, we introduce the notion of an interpolating vector field written in terms of iterates of the map. Our construction is based on Lagrangian interpolation and provides an explicit expressions for autonomous vector fields which approximately interpolate the map. We study properties of the interpolating vector fields and explore their applications to the study of dynamics. In particular, we construct adiabatic invariants for symplectic near identity maps. We also introduce the notion of a Poincar\'e section for a near identity map and use it to visualise dynamics of four dimensional maps. We illustrate our theory with several examples, including the Chirikov standard map and a symplectic map in dimension four, an example motivated by the theory of Arnold diffusion.


Introduction
Near identity maps naturally appear in several branches of mathematics and mathematical physics. They are often used for numerical integration of ordinary differential equations and for studying the evolution under fast time-periodic forcing. In the perturbation theory of integrable Hamiltonian systems with three degrees of freedom, near identity maps are used to describe dynamics in a small neighbourhood of a double resonance, see e.g. [6]. The last problem is one of the central topics in the study of the Arnold diffusion. In these problems some important details of the discussion depend on the smoothness class of the map and, for the sake of simplicity, we mainly consider the infinitely differentiable case and, where appropriate, assume analyticity. 1 It is well known that a near identity map F , where | | 1 and F 0 = I, is formally embedded into an autonomous flow [13]. Indeed, a formal construction can be used to find a vector field in the form of a formal power series in such that the formal series of the corresponding time-map coincides with the Taylor series of the original map. Truncated series define a vector field which interpolates the map with an error of order of the first omitted term. Note that this construction is formal only as in general the map cannot be embedded into an autonomous flow, i.e. generic F cannot be represented as a time-map of an autonomous vector field [7]. Nevertheless, the difference between them can be made a flat function in at = 0 [11].
An alternative approach to the study of near identity maps is based on averaging (see e.g. [7]). A near identity map F can be written as a time-map of a time-periodic vector field called a suspension of the map. Then the averaging procedure can be used recursively to eliminate the time-dependence up to any order of the small parameter. In the analytic case, a suspension can be found in the form of an analytic autonomous vector field plus an exponentially small non-autonomous term [7] (see also [2] for a concrete example).
Both these methods are constructive, nevertheless, none of them provides an easy way to find an explicit expression for a vector field which approximately interpolates the map. In this paper we propose a new method for construction of such vector fields with the help of the Lagrange interpolation of iterates of the map, henceforth referred to as interpolating vector fields. We analyse the errors of the approximation of F by the time-map of the interpolating vector fields.
The interpolating vector fields provide a useful tool for analysis of dynamics. In particular, we introduce the notion of a "Poincaré section" for the map F . Similarly to classical Poincaré section, this construction reduces the dimension of the phase space by one and, consequently, can be used to visualise trajectories of four-dimensional maps. Compared to the method of phase-space slices traditionally used for this purpose (see e.g. [9]) our method apparently shows higher resolution of details.
In the case when the map F is symplectic, an interpolating vector field can be used to construct an adiabatic invariant h n where n stands for the order of the interpolating vector field. The function h n is approximately preserved for many iterates of the map and can be used for a further reduction of the dimension.
In order to illustrate the effectiveness of the method and before discussing technical details involved in the construction, let us consider a 4-dimensional near identity Froeschlélike symplectic diffeomorphism T defined on T 2 ×R 2 . In coordinates (ψ 1 , ψ 2 , J 1 , J 2 ) this map is defined by the equation (3.5) but the precise expression is not important for the current discussion. Let Σ be a three-dimensional cylinder defined by the equality ψ 1 = ψ 2 . Fig. 1(a) shows points from several trajectories of T projected along an interpolating vector field onto Σ. We call this picture a Poincaré section for the map T .
In Fig. 1 all initial conditions are chosen from a torus defined by the intersection of Σ with a level set of the adiabatic invariant h 10 . Since h 10 is approximately preserved by the map, all points shown on Fig. 1(a) are in a small neighbourhood of the torus. Then a clearer image of the dynamics is obtained with the help two angular coordinates, namely, ψ and Note that all points are in a small neighbourhood of a two-dimensional torus. (b) A two dimensional projection of (a) using the coordinate φ = arg(J 1 + iJ 2 ) ∈ (−π, π] for the vertical axis. The analogy with the phase space of a two-dimensional area-preserving map is apparent. (c) A magnification of (b) where further details can be observed. φ = arg(J 1 + iJ 2 ). Fig. 1(b) shows the projection of Fig. 1(a) onto the torus (ψ 1 , φ), and Fig. 1(c) is a magnification of a strip located near the top of Fig. 1(b). These picture reveal that the four-dimensional dynamics of T resemble a two-dimensional area-preserving map. This similarity can be explained with the help of the normal form theory [8]. The latter also provides an alternative visualisation method (see [4]). We stress that the method of this paper is much easier to implement and provides higher accuracy in representation of details of the dynamics as our pictures represent true trajectories of T since the interpolation is used only for selecting initial conditions and projecting points onto Σ. The details of the model and the method are described in Section 3.
To give more formal definitions for the objects discussed in this paper, consider a smooth one-parameter near identity family of maps F : D → R m where m ≥ 1, D ⊂ R m is an open domain and | | < 0 . Then the map can be written in the form The function G 0 obtained by setting = 0 is often called a limit vector field. Its time-map is 2 -close to F .
The central object of this paper is an interpolating vector field for the near identity map F defined in the following way. Fix a positive integer n. Let x ∈ D and suppose that the iterates x k = F k (x) ∈ D are defined for −n ≤ k ≤ n. Then there is a unique polynomial p n of degree 2n in t such that for every t k = k with |k| ≤ n. The coefficients of p n depend on x and . The interpolating vector field is defined as the velocity vector of the interpolating curve at t = 0, i.e., where ∂ t p n denotes the derivative of p n with respect to t. In Section 2 we will use the Lagrange interpolation to show that the interpolating vector field X n is a linear combination of the iterates of the initial point x = x 0 , where Note that the coefficients p nk are independent of the map F , e.g., for n = 1 we get We will show that the time-map of this vector field is | | 3 -close to F . Thus X 1 can provide a more accurate interpolation for F than the limit flow.
If F : D → D is a diffeomorphism, then the domain D is invariant and the equation (1.3) defines X n on D for every n ∈ N. In general, we do not assume that D is invariant, and iterates of a point x 0 ∈ D may leave the domain. Nevertheless, it is not difficult to show that for any compact subset D 0 ⊂ D, there is a constant r 0 > 0 such that x k = F k (x) is defined for every x ∈ D 0 and every integer k provided | k| ≤ r 0 . Therefore X n is defined on D 0 for every n ≤ r 0 | | −1 .
If the map F is a lift of a diffeomorphism of a torus or cylinder, then the interpolating vector field is also a lift of a vector field defined on the same manifold. It is possible to generalize the definition onto a general manifold but this discussion is beyond the aims of the present paper.
We note that if F and F −1 are both polynomial (e.g. the classical Hénon map and its generalisations), then the interpolating vector field is also polynomial.
In the symplectic case, the interpolating vector fields remove the need of computing an averaged Hamiltonian. This opens an opportunity of performing massive numerical explorations of a multi-dimensional phase space, reducing the computation time significantly. Possible dynamical applications of this visualizing techniques include, among others, studying of area-preserving maps in C 2 and of Arnold diffusion for 4-dimensional symplectic maps near double resonances.
The rest of this paper is organised as follows. In Section 2 we briefly summarise the Lagrangian interpolation theory and derive some properties of the interpolating vector fields. Then we analyse the accuracy of interpolation. We discuss connections between the interpolating vector field and suspensions of the map F . In particular, we discuss connections to averaging methods and the optimality of choosing n as a function of . We show that in the analytic case the interpolation error can become exponentially small in . Finally, in subsection 2.5 we consider the symplectic setting. Although the interpolating vector field X n is not symplectic, it still can be used to define an adiabatic invariant h n which remains approximately constant for about −2n iterates of the map F .
In Section 3 we present results of our numerical experiments to illustrate the use of the interpolating vector fields for exploring the dynamics of maps in regions where they are near the identity. First, in Section 3.1 we consider the 2-dimensional symplectic case and perform experiments with the Chirikov standard map as a model. This is a preliminary example which illustrates various useful features of the interpolating vector fields. Moreover, one can compare the interpolating vector field directly with the map as the two dimensional dynamics is easy to visualize.
The main advantage of interpolating vector fields is that they provide a useful tool to investigate higher dimensional phase spaces. In particular, they provide a powerful tool for computing projections of the discrete dynamics onto subspaces of codimension-one. We use this tool to define Poincaré sections for maps, an extension of the classical construction traditionally restricted to systems with continuous time only, see details in Section 3.2. We illustrate this construction for the 4-dimensional symplectic map T to obtain representations of different types of its dynamics with the help of plots in dimensions two and three.

Interpolating vector fields 2.1 Lagrange interpolating polynomials
An interpolating polynomial can be explicitly written with the help of Lagrange interpolating polynomials (see e.g. [12]). Suppose that nodes are located at integer points symmetrically with respect to τ = 0. Then Lagrange's fundamental interpolating polynomial is and Lagrange's basis polynomials are .
Proposition 2.1. For every n fixed, (−1) k+1 p nk is a monotone decreasing positive sequence for 1 ≤ k ≤ n, and p n,−k = −p nk . Moreover, where H n = n k=1 k −1 is the n-th harmonic number.
Proof. The first equality can be derived from the following observation. Let 0 ≤ j ≤ 2n.
Interpolating any polynomial of degree j with the help of a polynomial of degree 2n, we recover the original polynomial. Applying this rule to τ j we get that n k=−n π nk (τ )k j = τ j .
Let r > 0 and consider a continuous curve γ : [−r, r] → R m . For any = 0 and n ∈ N such that n| | < r, the equation defines the unique polynomial of degree 2n such that p n (k , ) = γ(k ) for |k| ≤ n. Let v n (γ, ) = ∂ t p n (0, ). Taking into account the expression for p n we get v n (γ, ) = −1 n k=−n p nk γ(k ). (2.7) Note that although this definition contains a division by , the following lemma implies that v n can be extended to = 0 by continuity provided the curve is sufficiently smooth.
Lemma 2.2. If γ ∈ C 2n+1 [−r, r] and | |n < r, then Proof. The standard theory of interpolating polynomials, see e.g. [12], states that where the function w n can be expressed as w n (t, ) = γ (2n+1) (τ * (t, )) with τ * (t, ) ∈ [−n , n ], an intermediate value of the time variable. Since the zeroes of π are all simple, w is continuous at t = k . Differentiating with respect to t at t = 0 and taking into account that π n (0) = 0, we getγ The estimate of the lemma immediately follows from (2.3). This inequality implies that v n (γ, ) →γ(0) when → 0.
where H n = n k=1 k −1 is the n-th harmonic number. Proof. Using (2.7) and the triangle inequality we get Then the first equality in (2.6) implies the estimate of the lemma.

Basic properties of interpolating vector fields
The interpolating polynomial of the equation (1.1) can be explicitly written in the form of a Lagrange interpolating polynomial, namely, Differentiating this equality with respect to t at t = 0 and using that π nk (0) = p nk , we recover the explicit expression (1.3) for the interpolating vector field (1.2): For example, for n = 2 we get The definition of the interpolating vector field involves division by , but we will see that X n can be continuously extended to = 0.
Then there is r 0 > 0 such that the interpolating vector field is defined for all x ∈ D 0 and all n ∈ N such that n| | ≤ r 0 . For every n, the interpolating vector field X n is as smooth as the map F itself and X n (x, 0) = G 0 (x).
Then the finite induction implies that Repeating the previous arguments we conclude that Using the finite induction, the equations (2.10) and (2.11) and compactness of D 0 , we prove that there is Substituting the expressions (2.10) and (2.11) into the definition (2.9) of X n , and using (2.5) with j = 0, we obtain and the smoothness of X n follows directly. Finally, we note that if = 0, then x j = x 0 = x for every j. Substituting = 0 into the equation (2.12) we get Then the equation (2.5) with j = 1 implies that X n (x, 0) = G 0 (X).
An interpolating vector field can be used to approximately restore a vector field from its time-map. Proposition 2.5. Let Φ τ be a time-τ flow of a smooth autonomous vector field Y and n ∈ N. Let X n be an interpolating vector field for the map Φ . Then The estimate is uniform on every compact subset of D.
Proof. Let x ∈ D and γ be a solution of the initial value problem for the ordinary differential equationγ = Y (γ), γ(0) = x. Then Proposition 2.5 follows from Lemma 2.2 as the solution is smooth.
We note that the following arguments can also be used to prove the proposition. Expanding the solution γ of the initial value problem into Taylor series in time we get . Substituting this expression into the definition of the interpolating vector field we get Here we used the equation (2.5) to simplify the sum.
Remark 2.7. For a fixed n and a small the interpolating vector field X n (x, ) only involves values of F from a small neighbourhood of the point x.
Remark 2.8. If F is a lift of a map defined on a cylinder (or a torus), then the interpolating vector field X n is also a lift of a vector field defined on the same manifold. This argument also applies to every manifold obtained by factorising R m with respect to action of a discrete group of (affine) linear transformations.
We recall that a map F is called reversible, if there is an involution R (i.e., R • R = I) which conjugates the map and its inverse, i.e., In this paper we consider only the case when R is a linear map. This case is often used in applications. Proposition 2.9. If F is a reversible map with a linear reversor R, then the interpolating vector field is also reversible, i.e. the application of the map R changes the direction of time: This proposition is proved by a straightforward computation using the symmetric expression for the interpolating vector field (1.3) and the fact that the reversor R maps a trajectory of F into a trajectory of F −1 .
Let F be a map defined on an one-dimensional interval. Suppose that x 0 ∈ R is a fixed point of F and let λ = F (x 0 ). Then x 0 is an equilibrium for the interpolating vector field, i.e. X n (x 0 ) = 0, and If is sufficiently small, the linear stability type of the fixed point is the same for the map and for the interpolating vector field. Moreover, since F is a monotone function, it is easy to see that the time-map of the interpolating vector field is topologically conjugate to F . On the other hand, in general, the conjugating function cannot be differentiable as the multipliers λ and exp( X n (x 0 )) of the fixed points can be different. The conclusions about topological equivalence do not require hyperbolicity of F and do not generalize to higher dimensions.

Suspensions of a map and averaging
In general it is not possible to construct an autonomous vector field such that its timemap coincides with F . On the other hand, it is possible to construct a time-periodic vector field with this property. Such vector field is called a suspension of F . Suppose that Y is a suspension of F . This vector field has the following characteristic property: if ξ = ξ(τ, x, ) is a solution of the initial value problem It is not too difficult to construct a suspension. Indeed, let χ : Then consider a curve parametrized by the function with τ ∈ [0, 1]. This curve connects the points x and F ( . This map is -close to the identity in the C 1 topology and, hence, a local diffeomorphism by the inverse function theorem. Let Ψ τ = (Φ τ ) −1 and consider a non-autonomous vector field originally defined for τ ∈ [0, 1] and extended periodically as a function of τ . It is easy to see that Y is smooth. Obviously, the function (2.14) is a solution of the initial value problem (2.13) with the vector-field defined by (2.15). Consequently, the time-τ map of the vector field Y coincides with Φ τ for τ ∈ [0, 1]. Since Φ 1 = F , the vector field Y is a suspension of the map F .
Note that in this construction solutions of the differential equations (2.13) are given explicitly for τ ∈ [0, 1]. These solutions extend recursively to other values of τ using the identity ξ(τ, x, ) = ξ(τ − 1, F (x), ). Since the function χ is flat at τ = 0 and τ = 1, the function ξ is smooth in time as long as the solution remains in D. It is also possible to use the flow Φ τ (x) = F χ(τ ) (x) instead of (2.14) in the construction of the suspension.
A suspension of F is not unique and we are interested in finding a suspension which is as close to an autonomous vector field as possible. This task can be performed with the help of averaging [1]. An averaging step consists of a substitution of the form where the function S is periodic in τ . Substituting (2.16) into the differential equation (2.13) we get where S is evaluated at (τ, ξ 1 , ). Let Y = 1 0 Y dτ be the mean value of Y over its period in τ , andỸ = Y − Y be its oscillatory part with zero mean. It is convenient to let This function is periodic in τ and S(k, ξ, ) = 0 for every k ∈ Z. Then the map ξ 1 → ξ is -close to the identity for every τ and is exactly the identity for τ ∈ Z. Therefore, the change of variables (2.16) transforms Y into another suspension of the same map F . Multiplying the equation (2.17) by the matrix (1 + ∂ ξ S) −1 , we write the new suspension in the form Taking into account that S = O(| |), the definition ofỸ and differentiability of Y and S, we can check that the oscillatory part of Y 1 vanishes at the leading order, i.e., Ỹ 1 = O( Ỹ ).
The averaging procedure can be repeated to further decrease the time-dependent part of the suspension. Note that the expression for the vector field Y 1 contains a first derivative with respect to the space variable, thus if Y ∈ C k then we can claim that Y 1 ∈ C k−1 only and the maximum number of averaging steps is bounded by the smoothness class of Y . If Y is infinitely differentiable, the time-dependent part of the suspension vector field can be made smaller than any power of | | by repeating the averaging procedure multiple times: after n steps we obtain a suspension of the form Y n (τ, x, ) = A n (x, ) + n B n (τ, x, ).
(2.18) Neishtadt [7] proved that if Y is analytic in a complex neighbourhood of D then after n ∼ | | −1 averaging steps the time-dependent part of the suspension vector field becomes exponentially small compared to | |. I.e., if F is a family of maps analytic in a complex neighbourhood of D, then there is a suspension vector field defined on D such that with B = O(exp(−c/| |)) for some c > 0.

Suspensions of a map and interpolating vector fields
Although the averaging is an effective tool for studying near identity maps, it is not usually possible to find an explicit expression for the vector fields Y n . where the C 2n norms of A n and B n are bounded uniformly with respect , then for every compact D 0 ⊂ D there is a constant C n such that where X n is the interpolating vector field for the map F .
Proof. Since D 0 is compact, there is 0 > 0 such that every solution of with initial condition x ∈ D 0 remains inside D for |τ | ≤ n and | | ≤ 0 . Since Y is a suspension of F , we have F k (x) = ξ(k, x, ) for |k| ≤ n. Repeatedly differentiating the equation (2.21) with respect to τ and taking into account the form of Y , we see that there is a constant C such that ∂ k τ ξ(τ, x, ) ≤ C| | k for k ≤ 2n + 1. Then the theorem follows from the Lemma 2.2 with γ(t) = ξ( −1 t, x, ).
This theorem has several important corollaries. The first one states that the time-shift along trajectories of the interpolating vector field approximates F . Corollary 2.11. If F ∈ C 2n+1 and D 0 is a compact subset of D then on this subset the interpolating vector field is uniformly bounded for | | < 0 and where Φ Xn is the time-map associated with the vector field X n .
If F is analytic in a complex neighbourhood of D 0 , then we can choose n ∼ | | −1 in order to obtain a vector field which interpolates F with an error exponentially small in . To prove this, let Y be the suspension of F given by (2.19). So its non-autonomous part B is exponentially small in . Suppose that a solution of the autonomous equation ∂ t η(t, x, ) = A(η, ) with η(0, x, ) = x ∈ D 0 is analytic in time for |t| < 3r because it remains inside the domain of the vector field A for these values of t. Then the Cauchy estimate implies that |∂ k+1 t η| ≤ k!r −k A for |t| ≤ r. LetX n (x, ) be the interpolating vector field for the time-map of the flow of A. Applying Lemma 2.2 to the curve γ(t) = η(t, x, ), we obtain that for n| | < r
Since B is small, the time-t map of A is close to the time-t map of Y . In particular, Lemma 2.3 implies that where C A > 0 is a suitable constant so that C A B bounds the distance between the solutions of the initial value problems related to the vector fields A and Y for |t| ≤ r. Using the triangle inequality we get Since B = O(e −c/| | ) for a suitable c > 0 it follows that, taking n = [r/| |], one obtains Thus for this n = n( ) the vector field X n interpolates the map with an error exponentially small in .

Adiabatic invariants of a symplectic map
Let m = 2d and ω be a symplectic form on R 2d . We will mainly consider the standard symplectic form A map is called symplectic if it preserves ω. Even when F is symplectic, the corresponding interpolating vector field X n does not necessary define a symplectic flow and, consequently, is not necessarily Hamiltonian. Therefore there is no reason to expect that X n has an integral. Nevertheless, X n is very close to a symplectic one and can be used to define an adiabatic invariant -a function which is approximately preserved on time scales much longer than | | −1 .
In order to construct the adiabatic invariant, we consider a differential one-form ν n = ω(X n , ·). In the case of the standard symplectic form where X i n denotes a component of the vector X n . Then we fix a base point x b ∈ D and, for every x ∈ D, consider an integral where γ(x b , x) is a curve which connects the base point x b and x. If the form ν n is closed and the domain D is simply connected, this integral does not depend on the choice of the path (yet it depends on the choice of x b ) and uniquely defines the function h n : D → R such that dh n = ν n . In general, the form ν n is not necessarily closed and the integral may depend on the choice of the path. In order to obtain a well-defined function we choose a rule for selecting γ(x b , x). For example, if the domain D is star-shaped (e.g. convex), then γ can be a straight segment connecting its end points. Another convenient choice is a path which consists of straight segments parallel to coordinate axes.
In contrast with the interpolating vector fields, a suspension vector field (2.20) for symplectic map can be chosen (locally) Hamiltonian [7,5].
Proposition 2.12. Let F be a smooth family of near identity symplectic maps defined on a domain D ⊂ R 2d with d ≥ 1 and let n ∈ N. Let C > 0 be a constant and γ(x b , x) be piecewise smooth paths such that |γ(x b , x)| ≤ C|x − x b | for every x ∈ D. If a suspension of F can be written in the form of a Hamiltonian vector field with a Hamiltonian function where the C 2n+1 norms of h a n and h b n are bounded uniformly with respect and h a n (x b , ) = 0, then for every compact D 0 ⊂ D there is a constant C n such that where h n is defined by (2.23).
Let A n be the Hamiltonian vector field defined by the Hamiltonian function h a n . Then h a n (x, ) = Using the bi-linearity of the symplectic form ω, we obtain Since |ω(v 1 , v 2 )| ≤ C ω |v 1 | |v 2 | for any two vectors v 1 , v 2 (e.g. C ω = 1 for the standard symplectic form), we get The proposition follows directly from Theorem 2.10. Proof. Since F − Φ An = O(| | 2n+1 ), we get h n (F (x), ) = h a n (F (x), ) + O( 2n ) = h a n (Φ An (x), ) + O( 2n ). The flow of A n preserves h a n , in particular h a n (Φ An (x), ) = h a n (x, ), and the desired estimate follows.
This corollary implies that h n is an adiabatic invariant of the map F , as it is approximately preserved for −2n iterations of the map, provided the corresponding segment of the trajectory remains inside D. We expect that in the case of an analytic map, the adiabatic invariant h n( ) is preserved over exponentially long times where n( ) ∼ | | −1 .
We note that the function h n can be as smooth as the interpolating vector field (provided the paths γ(x b , x) are chosen appropriately). If the domain D is not simply-connected, the function h n can become multivalued if it does not return to the original value when x makes a round along a closed non-contractible loop.
3 Numerical study of dynamics using interpolating vector fields

Two-dimensional area-preserving maps: the Chirikov standard map
Our first example is the Chirikov standard map [3] written in the form M : (x, y) → (x,ȳ) = (x + ȳ, y − sin(x)), (3.1) where is a real parameter. We consider this map on the cylinder S 1 × R. On every bounded subset of the cylinder, the map (3.1) is close to identity provided | | is small enough. In our first numerical experiment we take = 0.1, n = 5 and study the interpolating vector field X n defined by equation (1.2). We choose some initial conditions on the cylinder and compare their iterates under the original map M with the iterates under the time-map Φ Xn corresponding to the interpolating vector field X n . Fig. 2 shows the first 10 3 iterates for 200 initial conditions. Both pictures use the same set of initial conditions. The interpolating vector field is integrated up to τ = 10 3 with the help of a Runge-Kutta-Felberg 7-8 integrator with variable step size. There is no visually perceptible difference between the two images.
In order to quantitatively describe the difference between the map and its interpolating flow, we consider interpolating vector fields X n for n = 5, 10, 15 and 20. Fig. 3 shows level plots of |Φ Xn (x 0 ) − M (x 0 )| computed for 5 × 10 5 initial conditions selected on a uniform mesh in [−π, π] × [−2π, 2π]. We clearly see that the error vanishes at the fixed points of the map and that the error decreases as n increases. Of course, the interpolation error will eventually grow with n due to Runge oscillations in interpolation. Remark 3.1. For the standard map (3.1), the interpolating flow X n , for any n, defines a one degree of freedom Hamiltonian system (with a non-standard symplectic form). This follows from the fact that the map is reversible, hence so is the interpolating vector field, see Proposition 2.11. The reversibility of X n forces the phase space to be foliated by periodic orbits, see Fig. 2 right and also the right plots of Fig. 4 as illustrations. A reversible 2dimensional system having a foliation of periodic orbits is Hamiltonian (possibly with a nonstandard symplectic form).
To visually inspect the differences between the orbits for the map and the interpolating flow, we increase the parameter up to = 0.5 and show the results in Fig. 4. The left plots represent the iterates of the standard map M , while the right ones correspond to the timemap Φ Xn for n = 10. We recall that X n defines an integrable vector field, see Remark 3.1. The bottom raw of Fig. 4 shows magnifications of a part of the pictures from the top row: we can clearly see the differences between the phase portraits when magnifying a strip located near a chain of resonant islands of the map.
To study the adiabatic invariants defined by the equation (2.23), we fix a base point x b ∈ R 2d (d = 1 in this section but we will also use d = 2 later) and consider paths represented by the function γ(s, x b , x) = x b + sv where v = x − x b and 0 ≤ s ≤ 1. Then the integral (2.23) takes the form where v i are components of the vector v. In numerical experiments, we evaluate this integral using a trapezoidal rule combined with the Romberg extrapolation scheme. We accept a numerical estimate of the integral value if the difference between two consecutive approximations of the Romberg scheme is less than 10 −8 . We use this rule to evaluate the adiabatic invariants in all examples presented in the paper unless otherwise stated. In principle, this method can be used to achieve higher precision, however for the visualization purposes higher precision is not needed.
To investigate the preservation of the adiabatic invariants as a function of n and under iterates of M , we select 10 4 initial conditions which form a uniform mesh in [−π, π] 2 and for every point compute ∆h n (x 0 ) = |h n (M (x 0 ), ) − h n (x 0 , )| and use the corresponding maximum value to estimate the supremum norm ∆h n . Fig. 5(a) shows plots of ∆h n as a function of log 10 for every 1 ≤ n ≤ 20. Each line of the plot corresponds to a fixed value of n and is obtained by joining 50 points with different values of ∈ [0, 1/2]. In this plot n = 1 corresponds to the largest values of ∆h n . This plot confirms that the upper bound of Corollary 2.13 is not too far from being optimal as the lines on the plot are approximately linear with the slope about 2n + 1 (till they reach the On the other hand, we see that ∆h n is not necessarily monotone in n for larger values of as the lines have intersection. Moreover, the plot indicates that for a fixed the value of ∆h n cannot be moved below a certain threshold by increasing the value of n (see Fig 5(b)). Therefore for a fixed we can find an optimal value of n which corresponds to a point after which the adiabatic invariant is not improved when n is increased. The existence of this threshold can be attributed to the non-integrability of the map M . A similar phenomena are observed in the study of optimal truncation of asymptotic series.
Finally, we note that the methods of this section can be used to study some maps which are not a priori near identity but have an iterate which is near identity on some subset of the phase space. In particular, this situation can arise in a study of a near integrable system near a multiple resonance. For example, if is not small, the map M is no longer close to identity. Nevertheless, in a neighbourhood of a q-periodic point, the map M q becomes close to identity. Let us illustrate how the interpolating vector field can be adapted to study the dynamics near a q-resonant chain of islands. Let = 0.5. We established (see Fig. 4) that the interpolating vector field provides an accurate approximation of the dynamics of M for y ∈ [−3, 3]. Now we consider larger values of y and investigate what happens with the approximation. We take initial points with x = π and iterate them. Comparing Fig 6 (a) and (b) which represent the dynamics of M and the interpolating vector field X n for n = 5, respectively, we see that the interpolating vector field does not correctly capture the dynamics around the 2-periodic chain of islands. On the other hand, in this part of the phase space the dynamics of M 2 is sufficiently close to identity so that the interpolating vector field X 2,n , computed from iterates of M 2 , provides a good approximation of the dynamics as can be seen in Fig 6(c). Note, that only one of the two islands can be seen due to the choice of initial conditions. Fig. 7 shows similar results for q = 3. Here we take initial points with x = 0, hence only one of the 3-periodic islands appears in the picture. The interpolating vector field X 3 n computed for M 3 accurately describes the dynamics in a narrow zone around the resonant 3-periodic island. Note that the 5-periodic island observed in Fig. 7(a) is not present in Fig. 7(b).

Exploring higher dimensional phase spaces: Poincaré sections for maps
If the dimension of the phase space m ≥ 3, the visualization of the dynamics becomes more difficult. In the case of a system with continuous time, a Poincaré section provides a convenient tool to reduce the dimension: a trajectory is represented by its intersections with a codimension one surface. In the case of discrete time, the reduction of the dimension  cannot be performed in a similar way. A typical solution to this problem is either to plot a projection of the trajectory on a subset of lower dimension, or to use the method of slices [9], i.e. to plot only a part of the trajectory which consists of points from a narrow strip near a codimension one surface (called a slice). In the last case, the points are also projected on the surface to reduce the dimension. The interpolating vector fields provide a new tool for visualization of dynamics which is especially effective in dimensions three and four. Suppose that g : R m → R is a smooth function such that its zero set Σ = {x ∈ R m : g(x) = 0} is a smooth hyper-surface of codimension one. Taking an initial condition x 0 ∈ D we compute the points x k+1 = F (x k ) recursively. The surface Σ locally divides the space. So we can look for consecutive points of the trajectory which are separated by Σ, i.e., Suppose that the limit vector field G 0 is transversal to Σ (at least in a neighbourhood of the intersection of Σ with the straight segment with endpoints at x k and x k+1 ). Then, for small enough, there is a unique t k ∈ [0, ] such that g(Φ t k Xn (x k )) = 0, and we plot the point instead of x k . Obviously, y k ∈ Σ and the trajectory is represented on a manifold of lower dimension. The point y k is O( ) close to x k and the error does not accumulate when k grows as the construction of y k does not affect the computation of the trajectory (x k ) k≥0 .

Four-dimensional symplectic maps: the Froechlé-like map
We apply the construction of the previous section to a 4-dimensional symplectic map and show that the method reveals interesting details of the dynamics. We consider the Froechlélike map T : where a 1 , a 2 , a 3 , η, are real parameters. The map T is a symplectic diffeomorphism of the cylinder M = T 2 × R 2 . It was introduced in [6] to model the dynamics near a double resonance in a near-integrable Hamiltonian system with three degrees of freedom. In our numerical experiments we use The quadratic form a 1 J 2 1 +2a 2 J 1 J 2 +a 3 J 2 2 is positive definite, since a 3 −a 2 2 > 0. The map (3.5) has four fixed points. If is positive and not too large, the origin p 1 = (0, 0, 0, 0) is ellipticelliptic, p 2 = (0, π, 0, 0) and p 3 = (π, 0, 0, 0) are hyperbolic-elliptic and p 4 = (π, π, 0, 0) is hyperbolic-hyperbolic.

A Poincaré section for T
Since the map (3.5) is symplectic, its limit flow is Hamiltonian. It is easy to find the corresponding Hamiltonian function explicitly: This Hamiltonian defines a non-integrable Hamiltonian system with two degrees of freedom. The Hamiltonian h 0 has four critical points which coincide with the fixed points of T . Levels of constant energy, M 0 E = {x : h 0 (x) = E}, are smooth hyper-surfaces for every regular value of h 0 . It is natural to study the dynamics of the limit system restricted on each energy level separately as the Hamiltonian function h 0 remains constant along trajectories of the limit flow. Let Σ be the 3-dimensional hyper-surface defined by the equality ψ 1 = ψ 2 . Note that we do not call Σ a "hyper-plane" because we treat ψ 1 and ψ 2 as angular variables. The limit vector field is transversal to Σ except for points which satisfy the equation (a 1 − a 2 )J 1 = (a 3 − a 2 )J 2 where the vector field is tangent to Σ.
The intersection Σ 0 E = Σ ∩ M 0 E defines a Poincaré section for the limit flow. Outside a neighbourhood of the tangencies, the first return map of the limit flow defines an areapreserving map on Σ 0 E . The dynamics of the limit Hamiltonian system are described by a collection of the Poincaré sections for different values of E. Since the quadratic form a 1 E is diffeomorphic to a three dimensional torus T 3 for every E > h 0 (p 4 ) = 1 + η. Then Σ 0 E is diffeomorphic to a two dimensional torus T 2 . It is convenient to use ψ = ψ 1 = ψ 2 and φ = arg(J 1 + iJ 2 ) as coordinates on Σ 0 E . In contrast to the limit flow, the map T does not have a first integral. Moreover, even if x 0 ∈ Σ, it is unlikely that x k = T k (x) will ever come back to Σ (periodic points of T are obvious exceptions) so the direct implementation of the Poincaré section is not possible. In order to visualise a trajectory x k = T k (x 0 ) we implement the procedure explained in Section 3.2. The procedure consists in finding a subsequence k j such that the trajectory jumps over Σ between x k j and x k j +1 . We note that Σ does not divide the cylinder M into two subsets globally but, since the map T is near identity, we can check this condition locally. Since Σ is defined by the equality ψ 1 = ψ 2 , we look for k j such that (ψ Then a point y k j ∈ Σ is defined by projecting x k j to Σ along the interpolating vector field X n as defined by the equation (3.4).
Since the section Σ is three dimensional, the sequence y k j can be plotted and used to visually inspect the behaviour of the trajectory x k . On a moderate time scale, a further reduction of the dimension can be achieved by noting that h n (defined by the equation (3.2)) is an adiabatic invariant of the map, so the trajectory x k stays in a small neighbourhood of the set M n E = {x : h n (x, ) = E} where E = h n (x 0 , ). Since M n E is close to M 0 E and the latter is nicely described by the coordinates (ψ, φ), we can project the points y k j on the torus of the coordinates (ψ, φ). In this way a trajectory of a 4-dimensional map is represented by a sequence of points on a two dimensional torus.
We remark that this procedure relies on the closeness of the map to the identity. Similar to the standard map, acceptable values of depend on the values of the variables J 1 and J 2 . In our numerical experiments we use | | ≤ 0.5, so in practical terms the parameter does not need to be very small.

Visualization of the dynamics of T
Examples of visualization of dynamics for T are shown on Fig. 1 and 8. Some comments concerning the implementation can help the reader. First, the computation of the points y k j on Σ requires integration of X n which is performed using a RK7-8 method that only requires evaluating the vector field. The time t k in (3.4) is then computed using the Newton method in a way similar to [10].
Second, to show different trajectories on a single 2-dimensional torus we select initial conditions on Σ n E = Σ ∩ M n E for a fixed E. To find initial conditions with the same value of E, we use the following procedure: we select values of ψ = 0, 1, 2, 3 and, for each value, we compute a point p = (ψ, ψ, 0, J 0 2 ), with J 0 2 > 0, such that h n (p, ) = E (using a bisection method with respect to J 2 to get a zero of h n (p, ) − E). Since ∇h n is orthogonal to the vector (0, 0, −∂h n /∂J 2 , ∂h n /∂J 1 ), we numerically integrate the auxiliary vector fielḋ with initial condition (J 1 (0), J 2 (0)) = (0, J 2 ) (using a RK7-8 method). One obtains points x 0,i = (ψ, ψ, J 1 (t i ), J 2 (t i )) ∈ Σ n E for a sequence of t i provided by the integration method. Finally, we use x 0,i as initial conditions for the map T . Fig. 8 shows 500 projected iterates y k j , as defined by (3.4), obtained from the iterates x k j under the map T for each of around 400 different initial conditions. The parameters of the map are defined by (3.6) and = 0.2. The hyperbolic-hyperbolic fixed point p 4 is used as a base point x b in the definition of the adiabatic invariant and the initial conditions are chosen on M n E with n = 10 and E = 1. We see that all points are located near a 2-dimensional torus embedded in the 3-dimensional surface Σ. The projection of the trajectories onto the coordinates (ψ, φ) resembles the dynamics of an area-preserving map. Fig. 9(a) shows iterates of the point (ψ 1 , ψ 2 , J 1 , J 2 ) = (3, 3, −1.043523, 1.385456) that belongs to Σ n E with E = 1. This is one of the chaotic orbits which can be seen in Fig. 8. Fig. 9(b) shows a plot of the adiabatic invariant as a function of time along this orbit. For better visualization we show the value of the adiabatic invariant for one out of every 250 consecutive iterates of the point under T . Some important observations follow from this plot: • The orbit we study apparently belongs to a chaotic zone and is not located on a KAM torus.
• The adiabatic invariant is preserved (meaning that the range of oscillations remains unchanged) up to, at least, 10 6 iterates. Fig. 9(c) shows that this behaviour continues at least for 10 10 iterates of T . The dynamical interpretation of such preservation is clear: our numerical method has detected the slowest variable of the system which evolves on a very long time scale. We note that in this example one expects to see an (Arnold like) diffusion process and, in particular, the slow variable requires exponentially long (with respect to −1 ) time to be changed by order one. Our method is able to detect such a slow variable in a simple and efficient way.
• There is no systematic drift of the adiabatic invariant. Its values are distributed in a Gaussian-like way around the initial value (see Fig. 9(d)). We remark that, for T and  the chosen parameters, there is numerical evidence supporting that the detection of the expected drift in the adiabatic invariant due to the Arnold diffusion process requires a much longer time scaling.
A final illustration of the amount of details one can visualize with this methodology is presented in Fig. 1. There we show a plot similar to the one in Fig. 8 but for = 0.35 and E = 4. As before, we take the hyperbolic-hyperbolic fixed point p 4 as a base point to define the adiabatic invariant. The interpolating vector field is constructed with n = 10. One can clearly recognize the typical structures that show up in a phase space of an area-preserving map: islands of stability (including secondary islands and satellites in the chaotic zone), invariant curves and chaotic zones. Yet this is not a 2-dimensional map: we are just plotting a suitable projection (along the solutions of the interpolating vector field) of the iterates of the 4-dimensional map! 4 Acknowledgments AV has been supported by grants MTM2016-80117-P (Spain) and 2014-SGR-1145 (Catalonia). VG's research was supported by EPRC (grant EP/J003948/1). The authors are grateful to Ernest Fontich and Carles Simó for useful discussions and remarks.