Efficient and Reliable Algorithms for the Computation of Non-Twist Invariant Circles

This paper presents a methodology to study non-twist invariant circles and their bifurcations for area preserving maps, which is supported on the theoretical framework developed in Gonzalez-Enriquez et al. (Mem. Amer. Math. Soc. 227:vi+115, 2014). We recall that non-twist invariant circles are characterized not only by being invariant, but also by having some specified normal behavior. The normal behavior may endow them with extra stability properties (e.g., against external noise), and hence, they appear as design goals in some applications, e.g., in plasma physics, astrodynamics and oceanography. The methodology leads to efficient algorithms to compute and continue, with respect to parameters, non-twist invariant circles. The algorithms are quadratically convergent and, when implemented using FFT, have low storage requirement and low operations count per step. Furthermore, the algorithms are backed up by rigorous a posteriori theorems, proved and discussed in detail in Gonzalez-Enriquez et al. (Mem. Amer. Math. Soc. 227:vi+115, 2014), which give sufficient conditions guaranteeing the existence of a true non-twist invariant circle, provided an approximate invariant circle is known. Hence, one can compute confidently even very close to breakdown. With some extra effort, the calculations could be turned into computer-assisted proofs, see Figueras et al. (Found. Comput. Math. 17:1123–1193, 2017) for examples of the latter. The algorithms are also guaranteed to converge up to the breakdown of the invariant circles, and then, they are suitable to compute regions of parameters where the non-twist invariant circles exist. The calculations involved in the computation of the boundary of these regions are very robust, and they do not require symmetries and can run without continuous manual adjustments, largely improving methods based on the computation of very long period periodic orbits to approximate invariant circles. This paper contains a detailed description of our algorithms, the corresponding implementation and some numerical results, obtained by running the computer programs. In particular, we include calculations for two-dimensional parameter regions where non-twist invariant circles (with a prescribed frequency) exist. Indeed, we present systematic results in systems that do not contain symmetry lines, which seem to be unaccessible for previous methods. These numerical explorations lead to some open questions, also included here.

ment and low operations count per step. Furthermore, the algorithms are backed up by

Introduction
The standard KAM theorem (see [18] and references there) shows that quasi-periodic solutions (geometrically invariant tori) in Hamiltonian systems which satisfy some non-degeneracy conditions are persistent under perturbations.
In such cases, where the non-degeneracy conditions of standard KAM theory are not satisfied, it is possible to establish the persistence by including extra parameters. The inclusion of extra parameters is natural in several problems, such as the design of accelerators or of confinement devices. In these problems, the designer can choose (within certain limits) the system that needs to be considered. In these applications, where stability is a main goal, one wants to modify the systems so that there are KAM tori guaranteeing the long-term stability. In design problems, by adjusting parameters, one aims to obtain not only the existence of KAM tori, but-if sufficient parameters are available-also extra properties that make the invariant tori effective barriers for transport.
With the above motivations, in [37] we introduced a rigorous methodology to study the persistence of non-twist tori, prescribing the qualitative properties of the degeneracies in the normal behavior. The main points of strength of the methodology include: -It is based on the theory of counterterms, also referred to as the method of parameters: One fixes the rotation and finds a modification (counterterm) of the system so that the modified system has an invariant torus, with the prescribed rotation. This technique was first introduced in [54] for close to integrable systems. See also [28], which includes a review of Herman's approach. -It uses symplectic geometry, so that the normal behavior of an invariant torus is described by the jet of a function (the potential) rather than a Birkhoff normal form. -The system is not assumed to be written as a close to integrable system. -It is formulated in an a posteriori format: Under the hypothesis of the existence of an approximate solution of the invariance equation which satisfies certain sufficient conditions, the existence of true solutions is proven in a constructive way. The sufficient conditions guaranteeing the existence of true solutions are made explicit in the form of conditions numbers (measuring the quality of the approximation). In this way, the classical close-to-integrable hypothesis is replaced by the assumption of the existence of an approximate solution of the invariance equation. Of course, for close-to-integrable systems, the exact solution of the integrable system can be taken as the approximate solution. Hence, we recover the customary quasiintegrable results as particular cases of our approach. Moreover, the a-posteriori format can also be used to validate the formal expansions obtained by matching powers. -As it is shown here, our theory leads to efficient and reliable algorithms for the computation of non-twist tori and detect bifurcations of invariant tori. -The method is numerically very robust.
-The theory guarantees that it can continue till the breakdown. -Given the a posteriori format of the rigorous results, we are assured that the computations are correct even close to breakdown. -The method does not rely on existence of symmetries or indicator points. -The method does not require any fine-tuning (as often happens in the methods based on approximation by periodic orbits).
One can see the results in Sect. 4, in particular in Fig. 8 The goal of this paper is to turn the rigorous results in [37] into practical algorithms, implement them and report some results and conjectures obtained by running the programs. The presentation of the algorithms is mainly for two-dimensional maps, that is, for computing invariant circles. The assumption of two-dimensional systems appear frequently in practice, and it allows some simplifications that were not used in the general results of [37], which were formulated for arbitrary dimension. Some interesting features of this paper are: -The algorithms we present are highly efficient: -They are based on Newton method and, hence, quadratically convergent. The non-degeneracy condition required to perform one iterative step is just the invertibility of a matrix, which is defined explicitly (see Algorithm 2). -To compute d-dimensional tori, one deals only with ddimensional functions (in the present paper, d = 1). This can be contrasted with methods based on normal forms, which rely on transformation theory that requires dealing with functions depending on as many variables as the dimension of the phase space (in the present paper, n = 2). It is well known that the computational effort grows exponentially fast with the dimension of the computed objects. -Low operation counts and storage requirements. The algorithms described here involve sequences of steps consisting on rather simple operations with periodic functions represented either in grid space (with N discretization points) or in Fourier space (with N Fourier coefficients). Each step has cost O(N ), in an appropriate representation, or O(N log(N )), when transforming one representation into another by using fast Fourier transform (FFT) algorithm. The storage is O(N ). This is a considerable improvement over the methods based on a direct discretization of the invariance equations that give rise to large systems for the Fourier coefficients, which require O(N 3 ) operations and O(N 2 ) storage.
-The algorithms are backed up by rigorous a posteriori theorems. This systematic assessment of the reliability of the computations is particularly important when studying invariant tori close to the breakdown. In these regimes, of course, the calculations become delicate and it is found empirically that there are many spurious solutions (i.e., solutions of the truncated equation which are not truncations of true solutions). With some extra level of effort (which we have not undertaken here), the a posteriori theorems allow to obtain computer-assisted proofs. It suffices to obtain rigorous estimates on the error of invariance, and the-rather mildnon-degeneracy assumptions. These can be obtained using interval arithmetic in function space. See [29] for a recent implementation. -Keeping certain parameters fixed, the continuation algorithms are guaranteed to converge (in machines with unlimited resources) as close as desired to values of parameters where the theoretical results (or the non-degeneracy conditions) fail to hold. This gives a practical criterion to find the limits of validity of the existence of non-twist invariant tori, which is based on the blow-up of Sobolev norms of the parameterizations at breakdown (see also [10]). In today's desktop machines (which are not infinite) in reasonable problems, such as the ones discussed in this paper, we estimate that one can get two or three decimal figures of the breakdown values in calculations of one hour or less. -The present method can be programmed once and for all: It does not require manual adjustments. This is an advantage over methods based on computing periodic orbits with large period, which often require considerable hand-tweaking of the algorithm to success to continue the right periodic orbits. An important example is the study of the breakdown of invariant circles, in which Greene's method [19,38,58,60] requires searching for periodic orbits. In practice, this requires the system to have some sym-metries, indicator points, etc. It has been observed that in systems with two harmonics, the computation of periodic orbits near breakdown is rather delicate [27,49]. In contrast, the Fourier methodology presented here is of general purpose and it is supported by rigorous results which apply, in principle, to regular regions of phase/parameter space (the hypotheses of our theorems include analyticity properties, which can be relaxed to finite differentiability properties). In these regions, the number of Fourier coefficients required by our methods to approximate an invariant circle with a good level of accuracy is relatively small (usually from 128 to 1024), while approximating the circles by periodic orbits would require periods of order 10 9 (for an accuracy 10 −9 ). In critical regimes (i.e., close to breakdown), this number increases abruptly: to obtain three or four digits of the breakdown parameter the required number of the Fourier coefficients easily increases to 10 5 or 10 6 , and several observables can be used to extrapolate the breakdown. In this paper, we provide numerical evidence that in absence of symmetries and close to breakdown, the algorithms presented here produce rather accurate results. Our methods are guaranteed to run till breakdown and, given the a posteriori results in [37], one can be confident that the computations are correct even close to breakdown. -Our method is based on computing a potential, in a parameter space, whose critical points correspond to the invariant tori. The potential is defined and can be calculated for any given symplectic map, regardless of the existence of invariant tori. Moreover, the critical points of the potential provide parameter values for which invariant tori (with a prescribed frequency) are likely to appear. This might be useful in some applications when the parameters are not subject to control, but they fluctuate (for example, in oceanography, we can predict where the barriers are likely to form).
Besides describing the algorithms, we also present implementation details: used discretizations, chosen linear solvers (not all linear solvers are practical for the problems considered here due to degeneracies in the spectrum), etc. In a final layer, we present the results of running the implementation in some problems that we believe would have been difficult by other methods: -The bifurcation of meandering non-twist circles; -Phenomena at breakdown; -The computation of the breakdown of non-twist invariant circles in families with several harmonics, where the periodic orbit theory seems complicated.
Aims of the Paper This paper requires working with different tools and standards: -Rigorous results are presented, mainly adapted from [37], but also new results in the area preserving case; -Numerical algorithms are presented; -Details of implementation of the algorithms are provided; -Numerical results obtained by running the algorithms are described; -Conjectures on the mathematical phenomena happening at breakdown obtained, from the numerical results, are presented.
We hope that the different topics covered in the different sections are sufficiently clear and that the reader appreciates having in a single paper the different tools which can bear on a rich problem coming from applied sciences. We certainly hope that this paper could inspire new further applications (e.g., the higher-dimensional case, more applied problems and more detailed studies of the breakdown phenomena). Organization of the Paper Section 2 tailors the geometric and analytic set up introduced in [37] to the area preserving case. In particular, we adopt a few simplifying assumptions: The phase space is assumed to be a two-dimensional cylinder, and the symplectic form is assumed to be the standard one. These assumptions make the theoretical results in [37] ready for the numerical implementation. We recall that a rotational circle in a cylinder is a no homotopically trivial circle. Section 2 can also serve as a guide to interpret the results of [37] in a concrete simplified case. Moreover, additional results are here shown to hold in the two-dimensional case. These are included in Appendix B. Section 3 presents some of the algorithms that can be derived from the methodology to compute invariant rotational circles. The rotation number and the type of degeneracy of a rotational invariant circle within a family of area preserving maps are assumed to be given and fixed. Section 4 reports numerical results obtained with the implemented algorithms for the cases mentioned above (fold bifurcation, meandering circles, etc.).

Methodology
This section describes our methodology to study degenerate invariant circles for area preserving maps. As mentioned before, this is an adaptation to the two-dimensional case of the methodology introduced in [37]. It turns out that the geometric properties of area preserving maps lead to simplifications with respect to the case of arbitrary dimension and arbitrary symplectic form. With the aim of implementing the algorithms produced by the methodology and their application to concrete problems, explicit formulas are included here, although details regarding the convergence are not given in this paper. (The latter are presented in full detail in [37].) To deal with degenerate rotational circles, some parameters (acting as translations) are introduced. In this way, the existence of invariant rotational circles is reduced to finding zeros of a function, which acts as a translation. It is worth mentioning that the translation can be computed even when the invariant rotational circle does not exist. It was proved in [37] that the translation is the negative gradient of another function, which we called potential, so existence of invariant tori is reduced to finding critical points of the potential. See Appendix B for the result tailored to the setting of this paper.
In the present paper, invariant rotational circles are found by finding zeros of the translation rather than finding critical points of the potential because in the planar case both problems are equivalent (the functions are univariate) and, from the numerical point of view, the former is less expensive and allows greater accuracy than the latter.

Classification of Invariant Rotational Circles
Let T×R denote the two-dimensional cylinder, with coordinates z = (x, y), endowed with the canonical symplectic form: Given z ∈ T × R, let z x = x and z y = y denote the projections over the x and y components, respectively. Let Φ : R × T × R → T × R denote the Hamiltonian action on T × R given by the horizontal translations: Notice that Φ is the flow generated by the Hamiltonian vector field with Hamiltonian which is homotopic to the identity and has the following form: where f , g : T×R → R are 1-periodic in x. The exactness condition on F guarantees the existence of a function S : T × R → R, known as the primitive function of F, such that: In particular, F is a zero net flux area preserving map. To avoid a large amount of notation, in referring to a symplectic map F, no notational distinction is made between the map and its lift (i.e., F is often assumed to be defined on R × R). We say that a circle K ⊂ T × R is a rotational circle if it is homotopic to the zero section, i.e., K is parameterized by an embedding K : T → T × R defined as follows: where x, y are 1-periodic in θ . Let · denote the average of a periodic function, given a rotational circle K , parameterized by K , the averages x and y are called, respectively, the phase and the momentum of K .
A rotational circle K is F-invariant, with rotation number ω ∈ R, if for a parameterization K of K , the following equation holds: Notice that if K is F-invariant with rotation number ω, then, for any parameterization K and any ϕ ∈ R, K ϕ (θ ) = K (θ + ϕ) is also a parameterization of the same invariant circle K . The phase of K ϕ is x ϕ = ϕ + x . Hence, by choosing ϕ = − x , it is always possible to obtain a parameterization with phase equal to zero. Hereafter, it is assumed that the parameterization of the circles has phase equal to zero. This normalization eliminates the ambiguity in the choice of the origin of coordinates on the circles, and it makes the translated circles locally unique under mild non-degenerate condition, as it is shown in [37].
Let U ⊂ R be an open interval, a family of F-translated rotational circles with frequency ω and labeled by the momentum p ∈ U is a map and such that the following equations hold: Notice that if K and p satisfy equations (2), with p ∈ U , then the exactness of F implies σ ( p) = 0 for all p (an area preserving map with zero net flux cannot translate a circle neither upward nor downward) and hence the following holds: with Φ defined in (1), thus obtaining a family of (horizontally) translated rotational circles. Hence, if K and τ = (λ, τ ) satisfy Eq. (2), with p ∈ U , then λ is referred to as the (horizontal) translation function of the family {K (·, p)} p∈U , with p ∈ U . Moreover, any zero of λ corresponds to a F-invariant rotational circle with frequency ω: if λ( p 0 ) = 0, then K (·, p 0 ) is the parameterization of an F-invariant rotational circle K p 0 . In the latter case, we say that K p 0 is non-degenerate if p 0 is a simple zero of the equation λ( p) = 0, and we say that it is degenerate if p 0 is a multiple zero. The multiplicity of the invariant rotational circle is the minimum integer k ≥ 1 such that λ( p 0 ) = 0, λ (1)

Remark 1
In [37], we proved that the translation function λ( p) is the negative gradient of a potential function V ( p) (see Appendix B.1), which is explicitly given by We emphasize that the construction works for families of translated Lagrangian tori for symplectic maps in arbitrary dimension. An important consequence is that invariant rotational tori correspond to critical points of the potential, and hence, a bifurcation theory of KAM tori can be derived from singularity theory of critical points of functions. In the two-dimensional case ( p ∈ R), rather than the potential V ( p), we use the translation function λ( p) because it is less expensive from the numerical point of view and gives greater accuracy.

A Quasi-Newton Step for Finding Translated Rotational Circles
In this section, we introduce a quasi-Newton method to find families of translated rotational circles with a given frequency ω. Within the setting of Sect. 2.1, this means finding a solution (K , τ ) of equations (2).
As it is standard in KAM theory, we assume ω to be Diophantine: There exists γ > 0 and ν ≥ 1 such that |qω − n| ≥ γ |q| −ν for all n, q ∈ Z with q = 0.
In what follows, we describe one step of a quasi-Newton method to solve Eq. (2). Assuming that (K , τ ) is an approximate solution of (2), a new approximate solution can be defined by (K ,τ ) = (K , τ ) + (ΔK , Δτ ), if (ΔK , Δτ ) solves the linearized equation where are assumed to be "small". Our goal is to find (ΔK , Δτ ) solving the linear Eq. (6), up to an error which is quadratic with respect to (E, e) (in suitable norms). If one defines then the linearized Eq. (6) can be written as follows: Moreover, the following hold: To avoid technical details that are not be used in the definition of the approximate solution of (6), we concentrate our efforts on finding a solution (U , u) of (8) under the following assumptions: M : T → R 2×2 is symplectic, ω ∈ R is Diophantine, V : T → R 2 and v ∈ R 2 are known, and there is L : T → R 2 such that L(θ ) = 0 for all θ and moreover, the following equations hold: Remark 2 From a dynamical point of view, (9) implies that L parameterizes an invariant bundle for the bundle map (M, R ω ) : From a functional analysis point of view, (9) implies that L parameterizes an eigensection of eigenvalue 1 for the transfer operator M acting on sections V : We note that the Newton method involves solving equations involving the transfer operators as above.

Remark 3
Invertibility of transfer operators is key on defining and implementing a Newton method. It is important to notice that the spectral properties of transfer operators are rather subtle. The spectrum presents certain peculiarities affecting the feasibility of the corresponding numerical algorithms. The spectral theory of transfer operators in spaces of C 0 functions was studied in [52] for rather general dynamics, where it was shown that the spectrum is invariant under rotations. The arguments in [52] do not generalize to the spectrum in spaces of more differentiable functions. Nevertheless, for bundle maps under irrational rotations, [42] shows that the spectrum acting on C r (r ∈ N ∪ {ω}) is invariant under rotations.
The invariance of the spectrum of transfer operators over rotations has several practical consequences (see, e.g., [1] for a discussion). In particular, the iterative methods based on selecting leading eigenvalues (such as Arnoldi iteration methods [4], based on constructing orthonormal basis of Krylov subspaces), do not work well in this case. The only alternative, we aware of, to the geometric procedure presented here is the use of the full matrix method [41], which consists in discretizing the linear equation into a full matrix equation [45].

Remark 4
Notice that in the case that E = 0 in 7, the bundle map (M, R ω ) is the linearized dynamics around a translated rotational circle (K , τ ), and L gives the tangent directions to the circle. This property is satisfied approximately when E = 0 but is small.
In what follows, we give a detailed description of our method to solve equations (8), under assumptions (9); explicit formulas are provided with the aim of making the resulting algorithm ready for implementation. First, define the following matrix-valued function: where Notice that for each θ ∈ T, P(θ ) is symplectic: P(θ ) Ω 0 P(θ ) = Ω 0 . In particular, P is invertible with inverse given by: Symplectness and invertibility of matrix P, defined by (10), are guaranteed by the fact that it is 2 × 2 dimensional. In the general dimensional case, P is only approximately symplectic and approximately invertible with approximate inverse given by (11) In (x, y)-coordinates, From (9) one obtains where In particular, the linearized dynamics (M, R ω ) around a translated rotational circle is reducible by the symplectic bundle map (P, Id) to a triangular bundle map (Λ, R ω ) with the identity in the diagonal. The torsion, defined byT = T , measures how much the normal bundle is twisted by the action of (M, R ω ), in average. We say that the circle is twist ifT = 0; otherwise, we say that the circle is non-twist. It was proved in [37] that twist tori are non-degenerate and non-twist are degenerate in the sense introduced in (4) (see Appendix B.2 for details in the two-dimensional case).
Define U (θ ) = P(θ )ξ(θ ), then performing some computations, equations in (8) are transformed into the following system of equations: where Let L ω denote the left operator acting on periodic functions ξ * : T → R as follows: Then, using the following notation enhancing tangent and normal components. Making explicit the tangent and normal components of the corresponding vectors, (12) is written as follows: where the unknowns are ξ L , ξ N , u x and u y . A fundamental result in KAM theory [56] is the following. If ω satisfies Diophantine conditions, then given a periodic function η * : T → R, the equation has an unique solution. We denote this unique solution by R ω η * (θ ) = ξ * (θ ), we refer to R ω as the right operator. The paper [56] contains a detailed study of the spaces in which Eq. (18) is solvable and the estimates satisfied by their solutions. Improved estimates appear in [29,30].

Remark 5
The operators L ω and R ω introduced above are diagonal in Fourier space: If then formally one has: This is important from both theoretical and practical points of view. Convergence of Fourier series R ω η * is possible thanks to Diophantine properties on ω and smoothness properties of η * . Numerical implementation of the action of L ω and R ω on (truncated) Fourier series is straightforward. Indeed, if η * is a (truncated) Fourier series with N harmonics, then computing ξ * requires O(N ) storage size and O(N ) number of operations.
Below we find explicit formulas for the solution ξ L , ξ N , u x and u y of Eqs. 14-17 in terms of R ω (see Algorithm 2).
From (15) and our hypothesis L x (θ ) = 1, L y (θ ) = 0, we obtain the following equalities: and where u x and the average ξ N 0 have to be determined. Letξ N denote the known part of then, substituting (20) in (14), one obtains For any ξ L 0 ∈ R, and ξ N 0 and u x such that the following defines a solution of Eq. (21): Replacing (23) into (17) and using the following notation: Equations (22) and (17) can be written as the following 2 × 2 linear system: It is clear that if the determinant ofT is different from zero, then ξ N 0 and u x can be computed as follows: Hence, from any ξ L 0 ∈ R, solutions of Eqs. 14, 15, 17 are defined by (19), (20)(23) and (25).
Remark 6 We refer to matrixT , defined in (24), as the supertorsion of the circle K with respect to the map F and the Hamiltonian action Φ given by the horizontal translations (see (1)). Note that the torsion T involves geometric properties only of the circle and the map, whereas the supertorsionT also involves the whole family of translated rotational circles (see (3)). We claim thatT is a symmetric matrix and moreover the following equality holds: (this result does not appear in [37]). This follows from the fact that for any periodic functions u and v, the following holds:

Remark 7
It is worth mentioning that if in (8), M is in a reduced form (as it happens in integrable systems): then, the supertorsionT , defined in (24), is non-degenerate no matter whetherT is non-degenerate or degenerate. Indeed, it is given by:

Existence of Translated Rotational Circles
Under suitable conditions, the iterative procedure based on the quasi-Newton method sketched in Sect. 2.2 is quadratically convergent, leading to an algorithm to compute translated rotational circles (i.e., solutions of (2)). This algorithm is discussed in Sect. 3.2. The convergence result is informally stated here without proof as Theorem 1. See Chapter 6 in [37] for a rigorous and more general statement, and a complete proof, for existence of families of translated tori for exact symplectic maps (with respect to a given exact symplectic form). It is also proved in [37] that translated tori (here rotational circles) depend smoothly on the momentum.
Theorem 1 Let F : T×R → T×R be a real-analytic area preserving map, homotopic to the identity and with zero net flux. Let ω ∈ R be a frequency satisfying Diophantine conditions and let p 0 ∈ R be a momentum. Let K 0 : T → T × R be a real-analytic parameterization of a rotational circle, and τ 0 = (λ 0 , σ 0 ) be a translation. Assume that: H2. The supertorsionT of K 0 (with respect to F and ω), defined in (24), is a nondegenerate 2 × 2 matrix.
Let (E, e) be the error of the approximately translated rotational circle (K 0 , τ 0 ) of momentum p 0 : If (E, e) is sufficiently small (in analytic norms), then:

T1. there exists a locally unique real-analytic family of F-translated rotational circles
with frequency ω, labeled with the momentum in a neighborhood U of p 0 , that is, a real-analytic mapK = (K , τ ) : T2. the vertical translation is zero: σ ( p) := τ y ( p) = 0 for all p ∈ U ; T3. the supertorsion of all translated tori in the family is non-degenerate.

Remark 8
Hypothesis H2 of Theorem 1, the fact that the supertorsionT in nondegenerate, is much weaker than the usual twist condition in KAM theory, that in this setting corresponds to the fact the torsion T is non-degenerate (see Remark 7). In a sense, the condition of existence of families of invariant tori (with a fixed frequency) in families of exact symplectomorphisms is weaker than the condition of existence of an invariant torus (with that frequency) for a particular exact symplectomorphism. In other words, given an exact symplectomorphism, by adding suitable parameters on can obtain invariant tori by tunning those parameters. The extra parameters (the translations) and the conditions (the averaged momentum of the tori) are geometrically related through the momentum map, as given in thesis T.1 of Theorem 1.
A quantitative version of Theorem 1 is suitable to perform rigorous proofs of the existence of (translated) rotational circles because the smallness condition on the error (E, e) was made explicit in [37]. The paper [29] contains a methodology to perform computer-assisted proofs in the context of KAM theory. We hope that it can be possible to extend the techniques in [29] to the problems discussed in this paper.

A Methodology to Study Bifurcations of Invariant Rotational Circles
This section describes the main ingredients of our methodology to study smooth bifurcations of invariant rotational circles.
Theorem 1 and the procedure described in Sect. 2.2 can be adapted to families of area-preserving maps. To study bifurcations of rotational circles for area-preserving maps, we assume that the family of area-preserving maps depends on parameters μ and ε. Parameters μ = (μ 0 , . . . , μ k−1 ) are used to unfold a zero, of multiplicity k, of the translation function, while ε is assumed to be a perturbation parameter which preserve the multiplicity k of the zero. More concretely, if the symplectomorphism F depends smoothly on parameters μ, ε, then the families of F-translated rotational circles, with fixed rotation ω, and the translation λ = λ( p; μ, ε) depend also smoothly on μ, ε. Hence, the study of bifurcations of invariant circles (with fixed frequency ω), which in principle is an infinite-dimensional problem, can be reduced to analyze the bifurcations of the zeros of λ( p; μ, ε), which is a finite-dimensional problem.
The procedure of computing and studying the bifurcations of rotational circles consists of three following main stages: Stage 1 For given parameters μ, ε * , fix a value p for the average momentum. Apply a quasi-Newton method (which is an infinite dimensional problem) to obtain: 1. a parameterization of the translate rotational circle with momentum equal to p (as well as other normalizations in the phase); 2. the value of the horizontal translation λ, which make the momentum of the translated rotational circle equal to p. Since the value of the momentum depends on the parameters p and μ and ε * , we will write λ( p; μ, ε * ).
Stage 2 For ε * fixed, run a finite dimensional Newton method to search for a parameter value μ * for which the function λ( p; μ * , ε * ) has a zero at p * , with the prescribed multiplicity k. That is p * is such that the following holds: To do so, one applies perturbation theory to compute the k-jet of λ( p; μ, ε * ) and its first-order derivatives (with respect to μ). Stage 3 Starting at ε * , continue with respect to ε, the solution ( p(ε); μ(ε)) of From the computational point of view, the most delicate step is Stage 1, because it involves solving an infinite-dimensional problem. Moreover, it involves dealing with small divisors and it has unbounded derivatives. Hence, issues such as the truncation of Fourier series are rather delicate, see Sect. 3.

Algorithms
This section contains the explicit description of the algorithms that implement the three-stage procedure introduced in Sect. 2.4.
The algorithm which computes the solution of the linearized Eq. (8) is formulated in Section 3.1. This algorithm is then applied to perform one step of our quasi-Newton method to compute translated rotational circles (Stage 1 in Sect. 2.4), see Sect. 3.2. Section 3.3 contains the algorithms for computing and continuing invariant rotational circles, which is done by finding simple zeros of the translation (μ is fixed). Although, it would be sufficient to apply a root-finding algorithm, such as the secant method, Brent's method [8], or Steffensen's method [7], we present an algorithm based on a Newton method. The case of zeros of multiplicity k > 1 is considered in Sect. 3.4. In this case, one needs to obtain the k-jet of the translation function with respect to the momentum and first-order derivatives with respect to parameters μ. In this paper, we only consider zeros of finite multiplicity; hence, only a finite number of computations is required, see Sect. 3.5. (For examples with zeros of infinite multiplicity, see [37]).
In Sect. 3.6, we discuss some implementation details.

Approximate Solution of the Linearized Equations
Our first algorithm solves the linearized Eq. (8), providing the basis for our quasi-Newton method to solve the invariance equations.

Algorithm 2 (Approximate solution of the linearized equations)
Let ω ∈ R satisfy Diophantine conditions. Assume that M : T → R 2 is a matrixvalued function and that L : T → R 2 satisfies the following conditions: where U : T → R 2 and u ∈ R 2 , through the following sequence of operations: The output (U , u) =R M,L,ω (V , v) of Algorithm 2 satisfies the following: where is the only solution of the linearized system (8).

Remark 9 If (K , ν)
is an approximate solution of (2), is an approximate solution of the linearized Eq. (6) (which is quadratic with respect to (E, e), in suitable norms).

Remark 10
Note that the number of steps in the algorithm is not very large (under 20). Moreover, each of the steps can be implemented in a few lines in a high level language or in a sophisticated numerical library. Therefore, a bare-bones implementation is not too difficult. Of course, a professional quality implementation that monitors the quality of the results requires more considerations and more sophistication.
The same applies to other algorithms discussed later that depend on this one to perform other related calculations.

Computation of Families of Translated Rotational Circles
Implementation of Stage 1 in Sect. 2.4 is performed by applying the algorithms introduced here: Algorithm 4 computes a family of translated rotational circles with fixed parameters, and then, Algorithm 5 performs the numeric continuation of the translated rotational circle.

Newton Method for Translated Rotational Circles
To compute a translated rotational circle and its corresponding translation, with prescribed momentum p and fixed Diophantine frequency ω ∈ R, we apply the quasi-Newton method described in Sect. 2.2. In particular, we apply Algorithm 2, see Remark 9.
LetR M,L,ω denote the output of Algorithm 2 (see equality (26)). Let p ∈ R be fixed, assume that (K , τ ) is an approximate solution of: Let (E, e) be the corresponding error: which is assumed to be small (in appropriate norms). (24), is non-degenerate, Algorithm (2) provides the correction ΔK and Δτ given by: The new approximate solution is then defined as follows: In [37], it was proved (Theorem 2.7) that under suitable conditions the error, corresponding to the new approximate solution, is quadratic with respect to (E, e) in analytic norms (see Remark 9). This is an important ingredient in proving the convergence of the quasi-Newton method described in Sect. 2.2.
Algorithm 3 implements one step of the Newton method describing the mathematical theory described in this section.

Algorithm 3 (Newton step for computing a translated rotational circle)
Let the momentum p be given. Assume that (K , τ ) is an approximate solution of (28) with invertible supertorsionT . Then, a new approximate solution (K ,τ ) is produced as follows: In practice, the previous algorithm needs to be driven by an algorithm that involves choices of thresholds for success and for unrecoverable failure.

Algorithm 4 (Newton method for computing a translated rotational circle)
Let the momentum p be given. Assume that (K , τ ) is an approximate solution of (28) with invertible supertorsionT . Then, repeat iteratively the following procedure (up to a given number of times): is too big (given a failure threshold), break and report failure; 4. If (E, e) is small enough (given a tolerance error), break and report success; 5. Otherwise, call Algorithm 3 (but, of course, steps (1) and (2) are already done).
Notice that the time taken by Algorithm 3 is completely predictable, but Algorithm 4 may require more or less iterations depending on the thresholds chosen. Then, the resources of storage are very easy to predict and they are linear in the number of discretization modes.

Continuation of the Family of Translated Rotational Circles
To continue the family of translated rotational circles (K (θ, p), τ ( p)), with respect to the parameter p (the momentum), we apply the quasi-Newton method introduced in Sect. 2.2, starting with a suitable approximate solution of 28, forp close to p, defined as follows.
Assume that the supertorsionT ( p) of K (θ, p) is non-degenerate (this is, in fact, the condition to apply the implicit function theorem to continue the family of translated rotational circles). Let and, since L(θ, p) = DK (θ, p) satisfy (9), we define Then, forp = p + δ p close to p, the approximate solution of 28 is defined by:
Notice that K p , τ p and λ p in Algorithm 5 are just the derivatives with respect to p derived above, in (29).

Computation and Continuation of Invariant Rotational Circles
This section contains the algorithms to numerically compute and continue (with respect to a perturbative parameter ε) invariant rotational circles. Later, we will discuss computation and continuation of invariant circles imposing a prescribed degeneracy.
Let F(·; ε) be a parameter family of exact symplectomorphisms. Assume that, for ε in a suitable open interval, {(K (θ, p; ε), τ ( p; ε))} p∈U is a family of F(·; ε)-translated rotational circles with Diophantine frequency ω: Let λ( p; ε) = τ x ( p; ε) denote the horizontal translation. It is clear that F-invariant rotational circles correspond to parameter values for which the following equality holds: The computation and continuation F-invariant rotational circles is done by continuing solutions of Eq. (30).

Computation of Non-Degenerate Invariant Rotational Circles
Let ε ∈ I be fixed, assume that p is an approximate solution of (30): A classical Newton method can be applied to solve (30), provided that e 0 is sufficiently small and λ p = λ p ( p; ε) = 0, i.e., the F-invariant rotational circle is non-degenerate (see Sect. 3.2.2). At each step, the correction is taken to be the solution of the corresponding linearized equation: where is computed in (29). Then, provided that λ p = λ p ( p; ε) = 0, Δp is given by: Algorithm 6 (Newton method for non-degenerate invariant rotational circles) Given a rotational circle K with momentum p and a translation τ such that (K , τ ) is approximately F(·; ε)-invariant, for ε fixed. A new approximately F(·; ε)-invariant rotational circleK with momentump and translationτ is defined by the following computations:

Remark 11
Another algorithm suitable for computing twist invariant rotational circles for exact symplectic maps follows from the KAM theory without angle-action angles introduced in [15,18]. See, e.g., [40] for actual implementations, and [29] for rigorous estimates. It is worth mentioning that Algorithm 6 can also be applied to non-exact symplectic maps to compute vertically translated rotational circles.

Remark 12
We emphasize that Algorithm 6 can be used to refine numerical results obtained by other methods, such as neighboring periodic orbits (using Fourier interpolation), or from efficient quantitative methods to compute rotation numbers [13,14,50,51].
First, we compute the partial derivatives with respect to ε of the involved functions (K , τ , p) as follows. From we then obtain that where dp dε is chosen so that dλ dε = 0. Then, a seed for the Newton method given by Algorithm 6, for finding a F(·; ε)-invariant rotational circle with frequency ω and momentum p(ε) is the triplet (K ,τ ,p) given by: Algorithm 7 (Predictor of the natural continuation method) Given a solution (K , 0, p, ε) of (32) for a value ε and an given a parameter step δε, an approximate solution (K ,τ ,p,ε) of (32) forε = ε + δε is computed as follows: We have again used the convention that subscripts p, ε indicate the variables with respect to which the derivatives are made.

Pseudo-Arclength Continuation Method
The continuation method specified in the previous section fails to continue turning points of the circle λ( p; ε) = 0. At the turning points, the invariant rotational circles are degenerate (non-twist). If the circle is regular, one of the derivatives ∂λ ∂ε , ∂λ ∂ p is different from zero. In this case, a pseudo-arclength continuation method can be used to follow turning points (and easily detect degenerate invariant rotational circles). Within the pseudo-arclength continuation method [46], the solution curve λ( p; ε) = 0 is considered to be parameterized by the arclength (notice that we consider (K , τ ) as objects depending on ( p, ε)).
The pseudo-arclength continuation method involves the following steps: 1. The corrector step: given an approximate solution, represented by a point in the parameter space that is close to the solution curve λ( p, ε) = 0, the corrector obtains the point in the solution curve which is the closest to the given approximation. This is done by applying Newton method to compute the minimum-norm solution of the corresponding linearized equations (notice that there are more unknowns than equations). 2. The predictor step: given a point ( p, ε) in the solution curve λ( p; ε) = 0 and a continuation step δs, the predictor obtains a new point in the parameter space which approximates a new point in the solution curve with arclength distance from ( p, ε) equal to δs.
Given an approximate solution ( p; ε) of (30), with λ( p; ε) = e 0 , the corrector step aims to obtain a correction (Δp; Δε) to ( p; ε) and a correction to the corresponding translated rotational circle (K , τ ), (ΔK , Δτ ). The linearized equation for (Δp; Δε) is the following: where the following notation is used: The solution that minimizes (Δp) 2 + (Δε) 2 (other weights in the Euclidean norm can also be used) is the one that solves the following linear equation: It is clear that to have a solution of 34, it is sufficient to assume λ 2 p + λ 2 ε = 0. As already showed before, the partial derivatives of K and τ are obtained as follows: and, in particular, λ p = τ x p and λ ε = τ x ε . The corrections for (K , τ ) are then computed as follows: The above procedure describes the corrector of the pseudo-arclength continuation method, which is summarized in the following algorithm.
If the initial seed is sufficiently close to regular points of the solution curve λ( p, ε) = 0, the repeated application of the previous step would converge to a point ( p, ε) in the curve, which corresponds to an invariant rotational circle for F(x, y, ε). In the predictor step, λ( p, ε) = 0 is considered to be a curve parameterized by the arclength, s, so that one has a function s → ( p(s), ε(s)), such that λ( p(s), ε(s)) = 0, dp ds (s)) 2 + ( dε ds (s)) 2 = 1.

Computation of Jets of the Translation Function
As explained before, in this paper the classification of an invariant rotational circle is based on the multiplicity of the corresponding momentum as a zero of the translation function. Hence, the continuation and unfolding of an invariant rotating circle is performed through the continuation and unfolding of the translation function with respect to the momentum. This, in particular, requires the study of the derivatives of the translation function with respect to the momentum and parameters, as it is shown in this section.
Let F(·; ε) be a parametric family of exact symplectomorphisms, with ε ∈ R (the case ε ∈ R follows the steps). Assume that for each ε in a neighborhood of parameter space, there exists a family of translated rotational circles (K (θ, p; ε), τ ( p; ε) = (λ( p; ε), 0)): We aim to compute derivatives of λ( p; ε) up to order k with respect to p and up to order 1 with respect to ε. This is equivalent to compute the k-jets of λ( p; ε) w.r.t. p (Taylor expansions up to order k) and the corresponding derivatives with respect to ε: Formally, the corresponding jets and the first-order derivatives with respect to ε of the corresponding translated rotational circles are given by: By applying the methods introduced in the previous sections, the zero-and firstorder terms of λ ≤k , λ ≤k,ε , K ≤k and K ≤k,ε can be computed (via Newton method). The computation of higher-order terms and the corresponding derivatives with respect to ε is performed by recurrence as explained below.
Given i ≥ 1, assume that the following have been already computed: -The (i − 1)-jet (K <i , τ <i ) := (K ≤i−1 , τ ≤i−1 ); -The derivatives of the (i − 1)-jet with respect to ε; -The (i −1)-jets and derivatives with respect to ε of the corresponding compositions with F: An iteration step, i.e., to compute (K ≤i , τ ≤i ), its derivatives with respect to ε, the i-jets of their compositions with F and the corresponding derivatives with respect to ε, is obtained by defining Indeed, the i-th-order terms in (36) satisfy: where δ i0 = 1 if i = 1, and δ i0 = 0 if i > 1. Then, θ, p; ε), ε) ≤i are computed. By differentiating (36) with respect to ε, and then taking terms of order i, we obtain In summary, we get

Computation and Continuation of Degenerate Invariant Rotational Circles
In this section, we consider the problem of computing degenerate invariant rotational circles, of multiplicity greater than 1. Of course, for multiplicity k > 1, the system at hand has to depend on a sufficient number of parameters, at least k − 1 unfolding parameters, which are denoted by the multidimensional parameter μ = (μ 0 , μ 1 , . . . , μ k−2 ) ∈ R k−1 . We also consider the problem of continuing such a degenerate invariant rotational circle with respect to a perturbative parameter ε.
Let F(·; μ, ε) be a multi-parametric family of exact symplectomorphisms, with parameters (μ, ε). Assume that a family of F-translated rotational circles has been computed for (μ, ε) in a given neighborhood. Let (K (θ, p; μ, ε), τ ( p; μ, ε)) denote such family, then the following holds: where λ( p; μ, ε) = τ x ( p; μ, ε) and τ y ( p; μ, ε) = 0. Such a family of F-translated rotational circles can be computed by applying Algorithm 4 in Sect.3.2.1. A degenerate invariant rotational circle of multiplicity k corresponds to values of the parameters for which the following equality holds in the (k − 1)-jet space: for the unknowns ( p; μ, ε). Using the coefficients of the jet, to solve Eq. (40) it is sufficient to solve the following system of equations: where See Sect. 3.4 for a methodology to compute jets.

Computation of Degenerate Invariant Rotational Circles
Let ε be fixed, and assume that an approximate solution ( p; μ) of (41) is known, so that are small errors. In what follows, we describe one step of a quasi-Newton method to find a solution of (41). To improve the approximate solution, we define a correction (Δp; Δμ), with Δμ = (Δμ 0 , . . . , Δμ k−2 ), which satisfies the following linearized equation: where we use the notation and we neglect quadratically small terms, involving products λ j Δp = e j Δp with j < k. Sect. 3.4 contains methodology for computing such derivatives. The solubility of (42) depends on the fact that μ is the unfolding parameter of the singularity and that multiplicity is k (λ k = 0). To solve (42) we proceed as follows. First, we solve the first k − 1 equations finding Δμ = (Δμ 0 , . . . , Δμ k−2 ). Then, we find Δp from the last equation,

Continuation of Degenerate Invariant Rotational Circles
Continuation with respect to the natural parameter ε of the degenerate invariant circle, requires the computation of derivatives with respect to ε of the all involved objects. To obtain a first-order approximation of these objects for a value close to ε,ε = ε + δε, we first consider the problem from a formal point of view and assume that there is a function ε → ( p(ε); μ(ε)) satisfying the following equalities: The (formal) derivatives of p(ε) and μ(ε) with respect to ε are found by solving the following equation: The coefficients in (43) are found by computing jets and their derivatives with respect to μ and ε. The procedure to solve Eq. (43) is similar to the procedure we used to solve (42). Uniqueness of the solutions of 43 follows from the fact that μ is assumed to be the unfolding parameters of the degeneracy of multiplicity k. After computing the derivatives, with respect to ε, of the momentum p and the unfolding parameter μ, the derivatives with respect to ε of the corresponding k-degenerate invariant rotational circle, K (θ ; ε) := K (θ, p(ε); μ(ε), ε), and the translation τ (θ ; ε) := τ (θ, p(ε); μ(ε), ε) (which should be zero) are computed as follows: where (K 1 , τ 1 ), (K 0,μ , τ 0,μ ) are assumed to be known (computed by Algorithm 10) and (K 0,ε , τ 0,ε ) is computed in a similar fashion. We summarize the findings of this section in the following algorithm.

Implementation of the Algorithms
On implementing the algorithms presented in this paper, it is convenient to represent periodic functions in two different ways: The operations with cost O(N ) done in grid space are: algebraic operations between (matrices of) periodic functions, composition of known functions (e.g., functions corresponding to the family of symplectic maps and its derivatives). Here we assume the cost for evaluating known functions at a grid point is O (1). See, e.g., [40] for implementation details of related algorithms.
The algorithms were coded using the C programming language, and the C library FFTW [32]. This library was used to perform the backward and forward FFT transformations between grid space and Fourier space.
In the implementations, two different accuracy controls were fixed: one (e.g., 10 −9 ) for the Newton methods (using 1 -norms in Fourier spaces and sup-norms in grid spaces) and other for the size of the tails of the truncated Fourier series, by making the harmonics of order higher than N /2 smaller than (e.g., 10 −11 ). These and other controls can be easily tuned. See, e.g., [40].

Results of Numerical Computations
In this section, we report some results obtained by running the implementations of the algorithms detailed in Sect. 3 in some concrete examples. Some of them have been considered by using methods different to those presented here. Moreover, we include some examples without symmetry lines, which seem to be inaccessible with methods based on the calculation of periodic orbits. We consider the following phenomena: detection of fold bifurcations (Sect. 4.1), birth of meandering rotational circles (Sect. 4.2), and continuation of non-twist rotational circles up to breakdown (Sect. 4.3). Examples in Sects. 4.1 and 4.2 focus on the detection of bifurcations in regular regions of parameters, and aim to illustrate the geometrical ideas underlying our methodology. In Sect. 4.3, we report the numerical results regarding the continuation, with respect to parameters, of a non-twist invariant rotational circle in a regular region of the param-eters. Moreover, we provide a critical boundary of this regular region. This boundary corresponds to values of the parameters for which the non-twist rotational curve breaks down. This preliminary study of the breakdown reveals several interesting features.
The results presented in this section are numerical computations, for which we will give error estimates that provide a certain degree of confidence. In order to validate rigorously these numerical results, one could follow the proofs of the a posteriori results presented in [37] to obtain explicit bounds for the errors ensuring the convergence of the algorithms (and then the existence of invariant tori), and then adapt the methodology introduced in [29] to perform computer assisted proofs of existence of KAM tori.

Detection of Fold Bifurcations
Fold bifurcations of invariant rotational circles occur when two invariant circles collide and annihilate each other when tuning control parameters. This is similar to the fold bifurcation of zeros of a function. To illustrate the fold bifurcation mechanism for invariant rotational circles, we consider a one-parameter family of quadratic standard maps F ε : where α is a fixed parameter and ε is the perturbation parameter. In this case, the primitive function of F ε is S(x, y; ε) = 2 3ȳ 3 + ε 4 π 2 cos(2π x) .
We are interested in finding invariant rotational circles with a fixed frequency ω.
In the following, we fix α = 0.375 and look for invariant rotational circles with frequency ω = 3− √ 5 2 . The choice is selected so that ω is a Diophantine number (it is a quadratic irrational, closely related to the golden mean ratio) and, for ε = 0, the map F 0 is integrable and possesses two invariant rotational circles with frequency ω. The rotational circles of the p-parameterized family K 0, p (θ ) = (θ, p) are horizontally translated by the corresponding potential function is given by: Hence, the two invariant rotational circles of F 0 are parameterized by K 0, p ± , where p ± = ± √ ω − α ±0.08346263 are the two zeros of λ 0 ( p) or, equivalently, the two critical points of the potential V 0 . The corresponding torsion isT 0 ( p ± ) = 2 p ± ±0.1669253: The rotational circles are twist (positive and negative, respectively). We emphasize that the supertorsion iŝ and, hence, its determinant is different from zero. We continued both twist invariant rotational circles, with respect to ε, using the algorithms 6 and 7. These circles are drawn in blue and red colors in the phase portraits of the left column of Fig. 1, for different values of ε. The two rotational circles smoothly collide at ε = ε c 1.361408, and the resulting invariant rotational circle is non-twist (i.e., its torsion is zero), shown in magenta in Fig. 1.
As already mentioned, the supertorsion is non-degenerate. This allows us to apply algorithms 4 and 5 to continue also a family of translated rotational circles, (K (θ, p; ε), (λ( p; ε), 0) ) (with frequency ω), for ε > 0. Computing each translation function λ( p; ε) for p ∈ [−0.1, 0.1] (with a step size 0.01), takes a few minutes. It is worth mentioning that families of translated rotational circles can be computed for values of ε greater than the collision value. The second and third columns of Fig. 1 display the translation function and the potential, respectively, for different values of ε. The zeroes of the function λ (or critical points of the potential V , obtained from (5)) correspond to momentum p for with the map F ε has rotational invariant circles with rotation number ω and such a momentum. The double zero of λ (degenerate critical point of V ) corresponds to the non-twist rotational invariant circle, as it is illustrated in the third row of Fig. 1. Hence, the observed bifurcation of rotational invariant circles corresponds to: (a) a saddle-node bifurcation of zeroes in a one-parameter family of one-dimensional equations (λ( p; ε) = 0); (b) a fold bifurcation of critical points of a one-parameter family of one-dimensional functions (V ( p; ε)).

Birth of Meandering Rotational Circles
In this section, we illustrate the methodology in an example exhibiting meandering invariant rotational circles, which are circles that cannot be described as graphs of the y component over the x component. These circles are not present in the non-perturbed system, but appear for positive value of the perturbation parameter by means of a mechanism known as reconnection scenario [44,60]. By the so-called second Birkhoff theorem [6,43,53], meandering circles do not appear in twist maps.
The model we consider in this section is a one-parameter family of quadratic standard maps F ε : T × R → T × R, defined by:    Fig. 2, where we mimic the information provided in Fig. 1. Then, we detect the fold bifurcation at which meandering rotational invariant circles are born. These circles are referred to as meanders because of their shape [59]. At the bifurcation value ε c 0.430396, the meander is non-twist.
Once we have detected the non-twist invariant rotational circle, parameters a or b can be tuned to perform a continuation of the non-twist invariant rotational circle. In the following section, we provide some implementation details to do so.

Remark 14
The rotational invariant circle at the bifurcation value in Fig. 2 is a non-twist meander circle. In [59], the author proves the existence of meanders for certain nontwist maps. It should be noticed that in our nomenclature, the meanders in the paper [59] are twist. Indeed, after performing appropriate changes of variables in certain regions, the system considered in [59] is written as a near-integrable twist map. Then, the existence of the meanders in [59] is guaranteed by the classical Moser's twist map theorem. Of course, the twist theorem can not be applied to prove existence of non-twist invariant circles, but a parameter version can.

Continuation of Non-Twist Rotational Circles up to Breakdown
In this section, we consider the continuation of a non-twist invariant rotational circle from the integrable case up to breakdown. This problem has been previously studied [19,20,23] by using the Greene's method in which, assuming reversibility and extra symmetries properties, non-twist invariant rotational circles are continued by computing stable periodic orbits with very large period. The breakdown of the non-twist circle is investigated by studying the stability properties of the approximating periodic orbits. In our experience with reversible systems, the presence of symmetries is crucial for the computation of periodic orbits with large period (see also [58] for the use of the so-called indicators points). Even with the use of the symmetry lines, it is hard to systematize the continuation of periodic orbits with large period. Moreover, if symmetries are not present, finding periodic orbits can be very cumbersome, even for twist area preserving maps [49].
The present methodology can be used to compute breakdown of twist and nontwist invariant circles, regardless of the existence of symmetries, and does not rely on the computation of periodic orbits. It is rigorously shown [17] that the present methods of continuation, given a computer with unlimited memory and running long enough, could reach arbitrarily close to the breakdown. (In real computers, of course it can get only very close as it is standard in all convergence proofs in numerical analysis.) Furthermore, in [10] it is shown that this breakdown happens if and only if one (Sobolev) H r norm of the parameterization of the invariant circle blows up. (In such a case, all H s norms with s ≥ r also blow up.) Note that, since we are computing Fourier series, H r norms are easy to compute. The criterion of blow up of Sobolev norms is a very practical tool to ascertain the breakup, and it is often much easier to implement than the classic Greene's method. The paper [10] contains a detailed comparison of several methods to compute the breakdown.
In the present paper, the continuation of a non-twist rotational circle is done by applying the algorithms introduced in Sect. 3.5 for degenerate rotational circles. In particular, we seek degenerate zeros of the translation function. The continuation is performed with respect to a perturbation parameter. To deal with the type of degeneracy, we use an unfolding parameter, which has to be adjusted during the continuation process to fix the order of degeneracy. In the case considered here, the multiplicity of the target zeros is k = 2.
To illustrate our methodology, throughout this section, we consider, for each ϕ ∈ T, the two parameter family of quadratic standard maps F ϕ,μ,ε : T×R → T×R, defined by: where ω = √ 5−1 2 is a fixed frequency, μ is the unfolding parameter and ε the perturbation parameter. The parameter ϕ breaks reversibility (notice that for ϕ = 0, the map is reversible). For ε = 0, the rotational circle y = 0 is a non-twist invariant rotational circle with frequency ω for μ = 0. Once we fix ϕ, the goal is to continue the nontwist rotational circle with respect to ε, and the corresponding unfolding parameter μ = μ(ε), up to the breakdown.

Remark 15
Similar models have been considered in the literature. The one proposed in the seminal works [19,20,23] is If ϕ = 0, then (46) is conjugated to the model in (47), with conjugacy given by We first illustrate the case ϕ = 1 8 . The numerical results are summarized in Table 1, where for each value of ε, the corresponding computed values for the following quantities, involved in the continuation, are reported: -The unfolding parameter μ; -The number of Fourier nodes N F ; -The determinant of the super-torsion, detT (which should be different from zero to be able to solve approximately the linearized equations); -The torsion T (which should be zero, since we are computing non-twist invariant circles); -The translation λ (which should be zero since we are seeking invariant circles); -The derivative of the translation, λ (which should be zero, since we are computing degenerate invariant circles); -Error estimates: the invariance error E and the reducibility error E red .
As it can be noticed in third column of Table 1, the number of Fourier coefficients for accurate approximations of the circles increases abruptly when approaching to the breakdown. Consequently, the computational time also increases. In slightly out of date desktop computers, 1 Figure 3 shows the Fourier spectrum of the two components of the parameterization of the non-twist invariant circle for ε = 2.870, giving a good numerical evidence of the convergence of the Fourier series.
The dynamics of the map (46), with ϕ = 1 8 , is illustrated in Fig. 4, for values of the parameters in the regular region (i.e., far from breakdown) and in Fig. 5 for values of the parameters in the critical region (i.e., close to breakdown). These figures illustrate: (left) the continuation of the non-twist rotational circle K (θ ) and the dynamics in the two connected components into which K (θ ) separates the phase space; (right) the tangent bundle, generated by DK (θ ) (in red), and the corresponding normal invariant bundle, generated byN (θ ) = N (θ ) − DK (θ )R ω T (θ ) (in blue), both represented by the corresponding polar angle α as a function of the θ coordinate. As mentioned in Sect. 2.2, the torsionT = T measures how much, in average, the normal bundle,N (θ ), is twisted by the action of (M, R ω ).
In the case under study, the torsion is zero because we are continuing non-twist rotational curves, from which one has that the normal bundleN (θ ) is invariant.
For parameter values close to breakdown (Fig. 5), it is interesting to observe how the chaotic sea grows as the twist invariant rotational circles are destroyed, but the non-twist invariant rotational circle seems to be more robust. The strange behavior of the invariant bundles announces the breakdown of the invariant circle. Sobolev norms H r have been used to detect the breakdown of non-twist invariant circles, monitoring their blow-up [10,11,16,31,40]. It has been shown rigorously that the breakdown of analytic tori happens iff and only iff all Sobolev norms of index r > r c blow up.
Furthermore, one of the predictions of the renormalization group picture is that all Sobolev norms H r , r > r c blow-up following a power law asymptotic of the form:  (46), with ϕ = 1 8 . For each parameter value ε, μ is the unfolding parameter, the left picture shows the invariant rotational circle and the dynamics, and the right picture shows the angle representation of the tangent bundle (red) and the normal invariant bundle (blue) as a function of the circle coordinate θ (Color figure online) Furthermore, the renormalization group theory of [10] predicts that B r is affine in r : The numbers B r and α, β are universal in the sense of renormalization group theory, that is, they should be the same for families in an open set and related by explicit formulas to properties of a fixed point (or periodic orbit) of a renormalization group operator.
Note that, in contrast with the blow-up, which is an unconditional rigorous result, the power law (48) depends on the assumption that the family considered is in the domain of universality of a renormalization group fixed point (or periodic orbit). Indeed, our numerical results, in particular Fig. 8, suggest that there are families which do not lie in the standard domains of universality.
We analyzed some Sobolev norms of the invariant rotational circles for (46), with ϕ = 1 8 . Figure 6 (top) shows the corresponding Sobolev norm functions h 2 and h 2.5 of the parameterization of the non-twist circle, and the corresponding fits. Figure 6 (down-left) shows the critical exponents B r as a function of the Sobolev regularity r . The glimpsed linear behavior suggests the renormalization group scenario described in [10], and that the critical Sobolev regularity is bigger than 0.5. Figure 6 (down-right) shows the dependence of the unfolding parameter μ with respect to ε, from which the critical value of the unfolding parameter, μ c , can be extrapolated. Table 2 displays several extrapolations of the critical parameters from the asymptotic laws, from which we conjecture the following critical values ε c 2.879 and μ c 0.01987 (fits have been performed with last 10 results displayed in Table 2).

Remark 16
If ϕ = 0 in (46) (the reversible case), our estimates for the critical parameters are ε c 3.8641 and μ c 0.068011, which correspond to the critical values a c 0.686044 and b c 0.74249 for the model (47). Our estimates follow from the analysis of the observables, see Fig. 7 and Table 3. In [19], the critical parameters reported are a c 0.686049 and b c 0.742493 that correspond to ε c 3.86411 and μ c 0.068015. These values were obtained by using refined versions of the Greene's method, which is a specialized algorithm for detecting breakdown for uniharmonic reversible non-twist and symmetric area preserving maps.
Continuation of non-twist invariant curves for (46), for several values of ϕ, can easily be handled by a computer cluster 2 . We performed the computations for the values ϕ = 5 · 10 −3 i, with i ∈ 0, 1, . . . , 199, and obtained a domain in parameter plane (ε 1 , ε 2 ) = (cos(2πϕ)ε, sin(2πϕ)ε) in which there exists a non-twist rotational circle with frequency ω for (46). Figure 8 Table 2 Extrapolation of critical parameters at breakdown for the system (46) in the case ϕ = 1 8 , and fits with power laws asymptotics (48). The exponents are fit by B r αr − β, where α = 0.801536, β = 0.458888 (see Fig. 6 Table 3 Extrapolation of critical parameters at breakdown for the system (46) in the reversible case ϕ = 0, and fits with power laws asymptotics (48). The exponents are fit by B r αr −β, where α = 1.045900, β = 0.920093 (see Fig. 7) r ε c μ c A B  The picture on the right shows a zoom of one of the peaks point was computed up to 3 digits (and extrapolated to 4 digits). The computational time, parameterizing the circles with Fourier series up to hundreds of thousands of coefficients, is only a few hours. Requiring one more precision digit results into the need of using millions of Fourier coefficients, which clearly lead to much more computation time. It is important to recall that in the whole process, our methodology does not assume any reversibility property.

Conclusions and Some Open Questions
In this paper, we presented several algorithms which can be applied to study bifurcations of invariant rotational circles, with fixed rotational number. The algorithms can be used to study any kind of finitely determined degeneracy class of invariant circle. Bifurcations of invariant circles are encoded in a parameter-dependent function, which we called translation function, so that bifurcations of invariant rotational circles correspond to zeros of the translation function. Our algorithms are also suitable to continue, with respect to parameters, non-twist invariant circles (with a specified degeneracy). We reported some numerical results produced by our algorithms in the simplest class of degeneracy: the fold singularity. To the best of our knowledge, the accurateness reported here has not been obtained before and we believe that it would be extremely hard to obtain it with techniques based on periodic orbits, due to the absence of symmetries and because periodic orbits may appear in complicated ways [49]. Moreover, the numerical results could be made rigorous by entangling the a posteriori theorems and methodologies presented in [29,37].
Even though our methodology has not been tailored for studying the mechanisms happening at the breakdown, in Sect. 4 we gave numerical evidence that our methodology can be used to provide observables that announce the presence of the critical phenomena. More precisely, we monitored Sobolev norms H r and their blow-up to compute values of the parameters corresponding to the critical curves.
The numerical results in Sect. 4 lead to several interesting open questions, all of them somehow related to each other. For example: 1. The boundary of the domain of existence in families of several parameters contains smooth pieces, but it also includes sharp corners. The smooth pieces can be explained by the standard renormalization group picture.
(The breakdown corresponds to the intersection of the family with a stable manifold of renormalization.) The corners can be explained as indication that the family considered leaves the domain of universality.
The question of what possible behaviors can the renormalization group experienceit is a dynamical system in infinite dimensions-is very open. 2. What can be said about the regularity of the invariant circles as they approach to the critical curve? Some preliminary results appear in [2]. We also note that the precise regularity of the conjugacies is a prediction of the renormalization group picture. 3. Is it possible to provide estimates of the critical exponent B r in (48) in terms of the degree of the Sobolev regularity r ?
From our numerical computations, we conjecture that such behavior of B r is affine with respect to r , for r greater than certain Sobolev regularity r c . This conjecture favors a renormalization group picture [10,11,16,31,40]. It would be also very interesting to obtain formulas for the coefficients of such linear behavior, and study the dependence with respect to the families, and so the existence of universality classes. Notice that the degenerate nature of the invariant circles approaching the critical circle makes the computation of the breakdown parameters very challenging. Summarizing, a related open question is: 3 How the universality classes and the renormalization groups depend on the type of degeneracy of the invariant circles?
Notice that classes of degeneracy higher to the fold singularity correspond to higher "flatness" of the frequency map around the circle, so a natural question, which is also of interest in plasma physics [3,9,[19][20][21][22][23]33,61] and oceanography [5,39,55,57]: 4 What is the influence of the "flatness" on the robustness of the non-invariant circle?
Last but not least, what about bifurcations of non-twist invariant tori with dimension greater than 1? The results in [37] hold in arbitrary dimensions, and the algorithms in this paper can be generalized to the setting in [37]. In that case, to study bifurcations of non-twist tori the potential function (and its critical points) is more suitable than the translation function (and its zeros). Of course, from the numerical point of view, questions 1-4 become more challenging in higher-dimensional cases. from which the result follows.

Remark 17
In terms of the potential of the family, Lemma 2 gives:
If the expression (50) does not vanish, the translated rotational circles are transversal with respect to their momentum p, guaranteeing the existence of a (local) diffeomorphism (θ, p) → K (θ ; p) giving rise to a (local) Lagrangian foliation. In particular, on a non-twist torus of momentum p,