Validated numerics for period-tupling and touch-and-go bifurcations of symmetric periodic orbits in reversible systems

TL;DR

This paper develops a computer-assisted framework to verify symmetry breaking and bifurcations of symmetric periodic orbits in reversible systems, with applications to celestial mechanics.

Contribution

It introduces a general method for verifying bifurcations in reversible maps and applies it to prove the existence of complex bifurcations of halo orbits in the Three Body Problem.

Findings

Proved existence of wide branches of halo orbits bifurcating from Lyapunov families.

Demonstrated multiple period doubling and quadrupling bifurcations.

Validated bifurcations for specific physical parameters.

Abstract

We propose a general framework for computer-assisted verification of the presence of symmetry breaking, period-tupling and touch-and-go bifurcations of symmetric periodic orbits for reversible maps. The framework is then adopted to Poincar\'e maps in reversible autonomous Hamiltonian systems. In order to justify the applicability of the method, we study bifurcations of halo orbits in the Circular Restricted Three Body Problem. We give a computer-assisted proof of the existence of wide branches of halo orbits bifurcating from $L_{1,2,3}$-Lyapunov families and for wide range of mass parameter. For two physically relevant mass parameters we prove, that halo orbits undergo multiple period doubling, quadrupling and third-order touch-and-go bifurcations.

Tables

Table 1. Table 1: Bound on Jacobi constant for period k 𝑘 k -tupling and touch-and-go bifurcations of L 1 subscript 𝐿 1 L_{1} halo families for the mass parameters μ E ​ M subscript 𝜇 𝐸 𝑀 \mu_{EM} and μ S ​ J subscript 𝜇 𝑆 𝐽 \mu_{SJ} . Multiplicity k = 1 𝑘 1 k=1 corresponds to symmetry breaking bifurcations. Multiplicity k = 3 𝑘 3 k=3 stands for touch-and-go bifurcations. In the case j = 7 𝑗 7 j=7 and μ = μ S ​ J 𝜇 subscript 𝜇 𝑆 𝐽 \mu=\mu_{SJ} , denoted by a star, the convexity condition in Theorem 15 has not been checked. We could not obtain bounds on third order derivatives so sharp, which would guarantee that the second derivative of Jacobi integral along bifurcation curve is non-zero.

j	k	bound on Jacobi constant (${\mu }_{EM}$)	bound on Jacobi constant (${\mu }_{SJ}$)
0	1	3.17435[03, 36]	3.03588[11,51]
1	4	3.08384097317512242430038839[79,91]	3.01979992774088569676042962[45,59]
2	3	3.058886412529835176423178819[4,6]	3.015412945713342018918305256[4,5]
3	2	3.0216192479264201986830801047[6,8]	3.00909434962110748263791018519[3,9]
4	3	2.999986911642456326104353768[3,8]	3.00589961847845578773000478[64,75]
5	3	2.997919501600216512922520[69,71]	3.00584514988954760426886596[07,61]
6	3	2.94132864491556775199[28,32]	2.9941342902214648929134[15,26]
7	4	2.940683922766931384[68,79]	2.9940756819941148370203478568[3,4]*
8	2	-0.986509091038502895183600231[6,9]	-0.99874596801401641137454972[14,21]
9	3	-0.996795335128162658942078721[1,5]	-1.0006666037342203004067305[44,52]
10	4	-1.004727349648878143369879[07,22]	-1.00227845861825336488127[00,44]
11	1	-1.016[09,14]	-1.00460[55,77]

Table 2. Table 2: Bound on Jacobi constant for period k 𝑘 k -tupling and touch-and-go bifurcations of L 2 subscript 𝐿 2 L_{2} halo families for the mass parameters μ E ​ M subscript 𝜇 𝐸 𝑀 \mu_{EM} and μ S ​ J subscript 𝜇 𝑆 𝐽 \mu_{SJ} . Multiplicity k = 1 𝑘 1 k=1 corresponds to symmetry breaking bifurcation. Multiplicity k = 3 𝑘 3 k=3 stands for touch-and-go bifurcations. In two cases j = 8 𝑗 8 j=8 and μ = { μ S ​ J , μ E ​ M } 𝜇 subscript 𝜇 𝑆 𝐽 subscript 𝜇 𝐸 𝑀 \mu=\{\mu_{SJ},\mu_{EM}\} , denoted by a star, the convexity condition in Theorem 15 has not been checked. We could not obtain bounds on third order derivatives sharp enough, which would guarantee that the second derivative of Jacobi integral along bifurcation curve is non-zero.

j	k	bound on Jacobi constant (${\mu }_{EM}$)	bound on Jacobi constant (${\mu }_{SJ}$)
0	1	3.15211[87,91]	3.03413[67,72]
1	4	3.071869946146936057096946[46,52]	3.01887850935398942392012584[04,38]
2	3	3.05083503865220863946099978[87,93]	3.014802161770970024871372364[53,86]
3	2	3.0229911591596336379776897124[1,5]	3.0091557491590008844773187520[3,6]
4	4	3.0158978159595140970471[78,93]	n/a
5	2	3.017143662542155479781194411[85,94]	n/a
6	4	3.0190035761795315910790[19,29]	n/a
7	3	3.077836508502735302[69,73]	3.019542146861187965925873562[20,98]
8	4	3.104289978034786552471088092[8,9]*	3.024350541079342605537853984[5,6]*

Equations

$$\mathrm{Fix}(f,U,\mathcal{C})=\{x\in U:\left(f(x)=x\right)\wedge \mathcal{C}(x)\}.$$

$$(R^{-1}\circ f\circ R\circ f)(x)=x.$$

$$S(f(x))=f(S(x)).$$

$${\mathcal{C}}_S(x)=\left(x\in \mathrm{Fix}(S,X)\right).$$

$$f_{\nu }^k(x)$$

$$x_{fp}:[{\nu }_1,{\nu }_2]⟶\mathrm{int}U,x_{b_1},x_{b_2}:[{\nu }^{\ast },{\nu }_2]⟶\mathrm{int}U,$$

$$\begin{array}[]{rcll}\mathrm{Fix}(f^{p}_{\nu},U,\mathcal{C})&=&\{x_{fp}(\nu)\},&\nu\in[\nu_{1},\nu^{*}],\ 1\leq p\leq k,\\ \mathrm{Fix}(f^{p}_{\nu},U,\mathcal{C})&=&\{x_{fp}(\nu)\},&\nu\in[\nu^{*},\nu_{2}],\ 1\leq p<k,\\ \mathrm{Fix}(f^{k}_{\nu},U,\mathcal{C})&=&\{x_{fp}(\nu),x_{b_{1}}(\nu),x_{b_{2}}(\nu)\},&\nu\in[\nu^{*},\nu_{2}],\\ \#\mathrm{Fix}(f^{k}_{\nu},U,\mathcal{C})&=&3,&\nu\in(\nu^{*},\nu_{2}].\end{array}$$

$$f_{\nu }^{k/2}(x_{b_1}(\nu ))=x_{b_2}(\nu ),f_{\nu }^{k/2}(x_{b_2}(\nu ))=x_{b_1}(\nu ).$$

$$x_{fp},x_b:[{\nu }_1,{\nu }_2]⟶\mathrm{int}U,$$

$$\begin{array}[]{rcll}\mathrm{Fix}(f^{p}_{\nu},U,\mathcal{C})&=&\{x_{fp}(\nu)\},&\nu\in[\nu_{1},\nu_{2}],\ 1\leq p<k,\\ \mathrm{Fix}(f^{k}_{\nu},U,\mathcal{C})&=&\{x_{fp}(\nu),x_{b}(\nu)\},&\nu\in[\nu_{1},\nu_{2}].\end{array}$$

$$f_{\nu }^k(x)$$

$$\mathrm{Fix}(R,U)=\left\{(p,q)\in U:p,q\in {\mathbb{R}}^n,q=0\right\}.$$

$${\pi }_q(f_{\nu }^k(p,0))=0$$

$$\mathcal{J}\ni \nu \to x_{fp}(\nu )=(p_1(\nu ),p_2(\nu ),0)\in \mathbb{R}\times {\mathbb{R}}^{n-1}\times {\mathbb{R}}^n,$$

$$\mathrm{Fix}(f_{\nu }^{2i},U,{\mathcal{C}}_R)=\{x_{fp}(\nu )\}.$$

$${\pi }_{q_2}(f_{\nu }^k(p_1,p_2,0))=0.$$

$$G_k(\nu ,p_1)={\pi }_{q_1}(f_{\nu }^k(p_1,p_2(\nu ,p_1),0)).$$

$$G_k(\nu ,p_1)=(p_1-p_1(\nu ))g_k(\nu ,p_1),$$

$$g_k(\nu ,p_1)={\int }_0^1\frac{\partial G_k}{\partial p_1}(\nu ,t(p_1-p_1(\nu ))+p_1(\nu ))dt.$$

$$\left\{(\nu ,p_1)\in \mathcal{J}\times P_1:g_k(\nu ,p_1)=0\right\}=\left\{(\nu (p_1),p_1):p_1\in P_1\right\},$$

$$0\notin \frac{{\partial }^2{G}_k}{\partial \nu \partial p_1}(\mathcal{J}\times P_1)+\frac{{\partial }^2{G}_k}{\partial p_1^2}(\mathcal{J}\times P_1)p_1^{\prime }(\mathcal{J}).$$

$$u(\tau )=g(w(\tau ),p_1(w(\tau )).$$

$$S={\int }_0^1\frac{\partial g}{\partial \nu }(u(\tau ))+\frac{\partial g}{\partial p_1}(u(\tau ))p_1^{\prime }(w(\tau ))d\tau .$$

$$\frac{\partial g}{\partial p_1}(\nu ,p_1)$$

$$\frac{\partial g}{\partial \nu }(\nu ,p_1)$$

$$\frac{\partial g}{\partial p_1}(u(\tau ))={\int }_0^1\frac{{\partial }^{2}G}{\partial p_1^2}(u(\tau ))tdt=\frac{1}{2}\frac{{\partial }^{2}G}{\partial p_1^2}(u(\tau )).$$

$$\begin{gathered} \frac{\partial g}{\partial \nu }(u(\tau ))={\int }_0^1\frac{{\partial }^{2}G}{\partial \nu \partial p_1}(u(\tau ))+\frac{{\partial }^{2}G}{\partial p_1^2}(u(\tau ))p_1^{\prime }(w(\tau ))(1-t)dt= \\ \frac{{\partial }^{2}G}{\partial \nu \partial p_1}(u(\tau ))+\frac{1}{2}\frac{{\partial }^{2}G}{\partial p_1^2}(u(\tau ))p_1^{\prime }(w(\tau )). \end{gathered}$$

$$S={\int }_0^1\frac{{\partial }^{2}G}{\partial \nu \partial p_1}(u(\tau ))+\frac{{\partial }^{2}G}{\partial p_1^2}(u(\tau ))p_1^{\prime }(w(\tau ))d\tau \in \frac{{\partial }^2{G}_k}{\partial \nu \partial p_1}(\mathcal{J}\times P_1)+\frac{{\partial }^2{G}_k}{\partial p_1^2}(\mathcal{J}\times P_1)p_1^{\prime }(\mathcal{J}).$$

$$0\notin \frac{{\partial }^2{G}_k}{\partial p_1^2}(\mathcal{J}\times P_1)$$

$$\frac{\partial g_k}{\partial p_1}({\nu }^{\ast },p_1^{\ast })={\int }_0^1\frac{{\partial }^2{G}_k}{\partial p_1^2}({\nu }^{\ast },p_1^{\ast })tdt=\frac{1}{2}\frac{{\partial }^2{G}_k}{\partial p_1^2}({\nu }^{\ast },p_1^{\ast })\ne 0.$$

Full text

In order to justify the applicability of the method, we study bifurcations of halo orbits in the Circular Restricted Three Body Problem. We give a computer-assisted proof of the existence of wide branches of halo orbits bifurcating from $L_{1,2,3}$-Lyapunov families and for wide range of mass parameter. For two physically relevant mass parameters we prove, that halo orbits undergo multiple period doubling, quadrupling and third-order touch-and-go bifurcations.

††journal: Communications in Nonlinear Science and Numerical Simulation

1 Introduction

In the past 40 years there were proposed very efficient methods for numerical continuation of periodic orbits for general ODEs, or periodic orbits satisfying certain symmetries with special focus on applications to Hamiltonian systems [19, 37, 22]. They are implemented in very efficient packages, such as AUTO [18], MATCONT [17] and CONTENT [27].

Most of the methods are similar in the spirit. The family of periodic orbits satisfies certain finite-dimensional implicit equations, which is solved by a Newton-like scheme. The equations usually involve period (return time) of the orbit, space variables and/or value of the Hamiltonian. The Newton-like iteration requires computation of monodromy matrix, which is not an issue nowadays, where many efficient tools for integration of ODEs along with variational equations are available [4, 23, 10].

Bifurcations of periodic orbits can be detected by looking at changes of determinant of monodromy matrix [54], inspecting so-called stability parameter in low dimensional systems [21, 22], analysis of normal form after Lyapunov–Schmidt reduction [16] or solving an equation specific for the type of bifurcation [38]. Literature on the topic is really wide, and the list of methods and references mentioned above is clearly incomplete.

The present paper is complementary to the above results and methods. We focus on reversible maps and reversible Hamiltonian systems. The primary result of the paper is a general framework for rigorous computer-assisted verification, that a branch of symmetric periodic orbits undergoes period-tupling and/or touch-and-go bifurcation. Finding an approximate numerical candidate for bifurcation point is just a preliminary step of the validation algorithm and for this purpose we can use any of the methods mentioned above [16, 22, 21, 38, 54]. Then, checking several inequalities on the (Poincaré) map under consideration and its derivatives on an explicit neighbourhood (input to the algorithm) of the approximate bifurcation point we can prove, that it contains a bifurcation point and we can conclude about type of bifurcation. As an output of the algorithm we obtain guaranteed bounds on both the bifurcation point and the parameter of the system, at which bifurcation occurs.

In order to justify applicability of the framework we apply it to the Circular Restricted Three Body Problem (CR3BP) [44]. We give a computer-assisted proof of the existence of wide branches of so-called halo orbits bifurcating from the families of planar Lyapunov orbits around $L_{1,2,3}$ libration points. We also prove, that for some physically relevant mass parameters of the system, these branches undergo period-tupling and touch-and-go bifurcations. These rigorous results are justification of some numerical observations from previous articles, in particular [21, 19]. We also would like to emphasize, that we have found a new phenomenon regarding bifurcations of halo orbits near $L_{3}$ libration point, when the mass parameter of the system tends to zero. We have observed, that the energy at the symmetry breaking bifurcation, which creates halo orbit, as a function of the mass parameter is not monotone, when the relative mass tends to zero — see Remark 19. Although not validated, this observation has been made using nonrigorous numerics with very high accuracy.

Rigorous numerical investigation of periodic orbits to ODEs became quite standard [2, 9, 11, 24, 25, 41, 45, 48, 49, 50]. To the best of our knowledge, there are very few results regarding computer-assisted verification of bifurcations of periodic orbits of ODEs, and the field remains widely open. Validation, that a family of periodic orbits undergoes a bifurcation usually involves computation of rigorous bounds on higher order derivatives of the trajectories with respect to initial condition (except some special cases when the system admits additional structure). The algorithm capable to do that appeared just in 2011 [53] (implemented as a part of publicly available CAPD library [10]) and to the best of our knowledge there are no other publicly available algorithms for rigorous integration of higher order variational equations. The $\mathcal{C}^{r}$-Lohner algorithm from [53] has been already applied to study period-dubling bifurcations of periodic orbits [51] in the Rössler system [43], homoclinic tangencies of periodic orbits in a time-periodic forced-damped pendulum equation [52] and non-local cocoon bifurcations [28] in the Michelson system [32]. Very recently [6], another approach to compute periodic orbits for ODEs without integration of the system has been proposed. Periodic solutions are approximated via piecewise Chebyshev polynomials and then their existence is validated by analysis of certain nonlinear operator on a Banach space of Chebyshev coefficients. The efficiency of the method is illustrated on the example of Equilateral Circular Restricted Four Body Problem. Similar approach has been successfully applied for validated computation of bifurcations of equilibria for ODEs and steady states for PDEs — see for example [8, 30, 31]. A different, geometric method for validation of bifurcations of steady states in the Kuramoto-Sivashinsky PDE is also proposed in [55].

Our algorithm for proving the existence of period-tupling and touch-and-go bifurcation is similar in the spirit to that proposed in [51, 55], but exploits the presence of a reversing symmetry of the system. After fixing an appropriate coordinate system in a neighbourhood of an apparent bifurcation points, we perform validated Lyapunov-Schmidt reduction. In this way, the analysis of bifurcation is transformed to zero-finding problem of a bivariate scalar-valued function, called bifurcation function. Then, checking some inequalities on derivatives of the bifurcation function we can prove, that its set of zeroes is the union of two smooth curves which intersect at an unique point — the bifurcation point. Finally, some non-degeneracy conditions (inequalities on higher order derivatives of the bifurcation function) let us to conclude about the type of bifurcation. In this paper we restrict to symmetry breaking, period-tupling and touch-and-go bifurcations.

The article is organized as follows. In Section 2 we introduce notation and main definitions used in the paper. Theoretical results, which constitute a basis of the computational framework for general reversible (Poincaré) maps are presented in Section 3 and then adopted to autonomous reversible Hamiltonian systems in Section 4. Finally, in Section 6 we give formal statements of theorems regarding continuation and bifurcations of halo orbits in the CR3BP.

In the Appendix we present two auxiliary yet new results, which are included mainly for self-consistency of the paper. Frequently we have to show, that a solution to an implicit equation is defined over an explicit domain. This step appears for example in continuation of periodic orbits or in validated Lyaunov-Schmidt reduction. Although there are available methods for this purpose (see for example [7, 20, 53, 9]), in A we provide an improvement, which takes advantage from higher order derivatives and flattening the implicit function by a smooth, well chosen substitution. In B we propose own and short algorithm for finding an approximate bifurcation point. It takes advantage from the Lyapunov–Schmidt reduction and computation of higher order derivatives of the (Poincaré) map under consideration. Using it, we could easily localize approximate bifurcation points in the CR3BP with the accuracy $10^{-60}$.

2 Preliminaries

2.1 Notation and basic definitions

For a map $f\colon D\subset X\to X$, a predicate $\mathcal{C}\colon X\to\{\texttt{true},\,\texttt{false}\}$ and a set $U\subset D$ we introduce the following notation

[TABLE]

Thus, $\mathrm{Fix}(f,U,\mathcal{C})$ is the set of fixed points of $f$ in $U$ satisfying constraint $\mathcal{C}$. Although there is a natural correspondence between subsets of $X$ and univariate predicates defined on $X$, we will separate $U$ and $\mathcal{C}$ to emphasize rather rare property $\mathcal{C}$ in a larger (usually open) set $U$. We will also write $\mathrm{Fix}(f,U)$ if $\mathcal{C}(x)\equiv\texttt{true}$. For a set $M\subset\mathbb{R}\times X$ and $\nu\in\mathbb{R}$, we define its slice by $M_{\nu}=\{x\in X:(\nu,x)\in M\}$. For a map $f\colon M\subset\mathbb{R}\times X\to X$ and fixed $\nu\in\mathbb{R}$ we define $f_{\nu}\colon M_{\nu}\to X$ by $f_{\nu}(x)=f(\nu,x)$.

Definition 1

A homeomorphism $R\colon X\to X$ is called a reversing symmetry for $f\colon M\subset X\to X$ if $(R\circ f)(M)\subset M$ and for $x\in M$ there holds

[TABLE]

It is easy to see that any reversing symmetry $R$ for $f$ is also a reversing symmetry for all iterations $f^{n}$ defined on their proper domains.

Definition 2

A homeomorphism $S\colon X\to X$ is called a symmetry for $f\colon M\subset X\to X$ if $S(M)\subset M$ and for $x\in M$ there holds

[TABLE]

For a map $S\colon X\to X$ we define a predicate $\mathcal{C}_{S}\colon X\to\{\texttt{true},\,\texttt{false}\}$ by

[TABLE]

We will use this notation to select points satisfying certain symmetries.

2.2 Geometric definitions of two types of bifurcations

Our primary object of interest is a $\mathcal{C}^{3}$-smooth function $f\colon M\to X$ defined on an open set $M\subset\mathbb{R}\times X$, where $X$ is a smooth manifold. In this article we make the following standing assumption:

• C1:

for $\nu\in\mathbb{R}$, if $M_{\nu}\neq\emptyset$ then the function $f_{\nu}\colon M_{\nu}\to X$ is a diffeomorphism onto image.

The above assumption fits the applications we have in mind, i.e. $f_{\nu}$ will be a one-parameter family of Poincaré maps. Thus, the domain of each map $f_{\nu}$ may vary with the parameter $\nu$.

Bifurcations are usually defined by their normal forms [13]. The following

[TABLE]

describe period-tupling and transcritical bifurcations, respectively. When an eigenvalue of the derivative of a reversible planar map at a symmetric fixed point crosses $1:k$ resonance, $k\geq 2$, generically period-tupling or touch-and-go bifurcation occurs (notation following [5]) — two types of generic bifurcations for strong resonances are illustrated on Figure 1. The branches of symmetric periodic points near bifurcation point look similarly to those for pitchfork and transcritical bifurcations, respectively.

This observation lead us to the alternative definitions, in which these two types of bifurcations are described by some geometric conditions on the mutual position of two intersecting curves, that solve equation $f^{k}_{\nu}(x)-x=0$ — see Figure 2. These definitions are motivated by our algorithm for validation of bifurcations (Section 3), which is geometric in its spirit. Then, in Section 3.4 we will show, that these geometric definitions along with some non-degeneracy conditions imply standard unfolding (2)–(3) of these bifurcations.

Definition 3

Let $\mathcal{C}\colon X\to\{\emph{{true},\,{false}}\}$ be a predicate and let $k$ be a positive integer. We say that $f_{\nu}\colon M_{\nu}\to X$ has period $k$–tupling bifurcation at $(\nu^{*},x^{*})\in M$, if there exists $V=[\nu_{1},\nu_{2}]\times U\subset M$, such that $(\nu^{*},x^{*})\in\mathrm{int}\,V$ and there are continuous and smooth in the interior of their domains functions

[TABLE]

such that $x_{b_{1}}(\nu^{*})=x_{b_{2}}(\nu^{*})=x_{fp}(\nu^{*})=x^{*}$ and the following conditions are satisfied.

Periodic points:

[TABLE] 2. 2.

If $k$ is even, then for $\nu\in(\nu^{*},\nu_{2}]$ there holds

[TABLE]

Definition 4

Let $\mathcal{C}\colon X\to\{\emph{{true},\,{false}}\}$ be a predicate and let $k$ be a positive integer. We say that $f_{\nu}\colon M_{\nu}\to X$ has $k^{\mathrm{th}}$–order touch-and-go bifurcation at $(\nu^{*},x^{*})\in M$, if there exists $V=[\nu_{1},\nu_{2}]\times U\subset M$, such that $(\nu^{*},x^{*})\in\mathrm{int}\,V$ and there are continuous and smooth in the interior of their domains functions

[TABLE]

satisfying the following conditions.

Periodic points:

[TABLE] 2. 2.

The curves $x_{fp}$ and $x_{b}$ intersect at exactly one point $x^{*}=x_{fp}(\nu^{*})=x_{b}(\nu^{*})$.

The above two definitions are slightly more general than the normal forms (2)–(3). Indeed, it is easy to see that the functions

[TABLE]

also satisfy geometric conditions from Definition 3 and Definition 4 for any $n,m\in\mathbb{N}$.

3 Validation of bifurcations of symmetric periodic orbits: reversible maps

Let $f\colon M\to X$ be a family of $R$–reversible maps satisfying our standing assumption C1. In this section we propose a general framework for computer-assisted verification that a branch of $R$-symmetric period-$2$ points for $f_{\nu}$ undergoes symmetry breaking, period-tupling or touch-and-go bifurcation. Let $V\subset M$ be an explicit neighbourhood of an approximate bifurcation point (both are input to the algorithm). The algorithm is split into two steps.

In the first step we validate, that the set of points $(\nu,x)\in V$, such that $x\in\mathrm{Fix}(R)$ is either period-$2k$ point or period-$2$ point for $f_{\nu}$ is a union of two regular curves intersecting at exactly one point. This is a common step for period-tupling and touch-and-go bifurcations and it will be described in Section 3.1. In the case of symmetry breaking bifurcation, we need an additional constraint, which allows to isolate the primary curve of period-$2$ points and the branching-off curve of period-$2$ points. In Section 3.3 we will show, that such a problem can be solved, if the system has an additional symmetry, which commutes with $R$.

The second step is specific for each type of bifurcation. We provide checkable by means of rigorous numerics conditions, which guarantee, that the two curves intersect in a manner given by Definition 3 and Definition 4.

Finally, in Section 3.4 we show that some additional non-degeneracy conditions lead to standard unfolding of the computed bifurcations to their normal forms (2) and (3).

3.1 Bifurcation as two intersecting curves

We are mostly interested in studying bifurcations of symmetric periodic orbits of $R$–reversible ODEs. In continuous-time $R$–reversible dynamical systems, a trajectory is $R$–symmetric, if it is either an equilibrium belonging to $\mathrm{Fix}(R)$ or it intersects twice the set $\mathrm{Fix}(R)$ — see [29]. In order to obtain isolated symmetric periodic orbits we assume, that $\mathrm{dim}\,X=2n$ and $\mathrm{Fix}(R)$ is an $n$-dimensional submanifold.

Let us fix $k\geq 1$ and assume that $(\hat{\nu},\hat{x})$ is a good numerical approximation of a bifurcation point, i.e. one of the eigenvalues of $Df_{\hat{\nu}}(\hat{x})$ is close to $1:k$ resonance. Since all types of bifurcations we are studying are local phenomena, we assume, that there is an open interval $\mathcal{J}$ and an open set $U\subset X$, such that $(\hat{\nu},\hat{x})\in\mathcal{J}\times U$ and for $\nu\in\mathcal{J}$ the function $f_{\nu}^{2k}$ is defined on $U$. We also assume that there are local coordinates in $U$ in which

[TABLE]

We will use the same symbols $(p,q)$ as a coordinate system near $f_{\hat{\nu}}(\hat{x})$. Solutions to

[TABLE]

are $R$–symmetric periodic points of $f_{\nu}$ with principal period not larger than $2k$. The solution set to (4) near $(\hat{\nu},\hat{x})$ is expected to be a union of two curves, one of which corresponds to the fixed points of $f_{\nu}^{2}$ and the second corresponding to period-$2k$ points of $f_{\nu}$.

Remark 5

There are rigorous numerical methods for validation, that a solution set to an implicit equation forms a regular curve over an explicit domain — see for example [9, 7, 51]. In this paper we propose slightly improved version of the Interval Newton Operator [1, 34, 39]. Since it is only an auxiliary step of the main algorithm for validation of bifurcations, we postpone this improvement to A.

The case $k=1$ (symmetry breaking bifurcation) will require additional constraint and will be discussed in Section 3.3. If $k>1$, we may apply the method described in A to solve (4) on an explicit range of parameter values $\nu\in\mathcal{J}$. In what follows we assume, that there is a smooth curve

[TABLE]

such that for $\nu\in\mathcal{J}$ and $i=1,\ldots,k-1$

[TABLE]

In most cases, the function $x_{fp}$ cannot be computed exactly, but using rigorous numerics we can prove, that it exists and we can find bounds for $x_{fp}(\nu)$ and its derivatives. Although it is not required, it is desirable for further numerical computation to choose the coordinate system in $\mathrm{Fix}(R,U)$ and the decomposition $p=(p_{1},p_{2})$ so that $p_{1}^{\prime}(\hat{\nu})\approx 0$.

The idea of computing the second curve of $R$-symmetric, period-$2k$ points, which solves (4) goes as in [51]. At first, we perform the Lyapunov-Schmidt reduction [13]. We split $q=(q_{1},q_{2})\in\mathbb{R}\times\mathbb{R}^{n-1}$ and using the method from A we solve for a function $p_{2}=p_{2}(\nu,p_{1})$ satisfying implicit equation

[TABLE]

Assuming that $p_{2}=p_{2}(\nu,p_{1})$ is locally unique solution to (7) in $\mathcal{J}\times U$, we can define so-called bifurcation function

[TABLE]

In this way we reduced the problem of finding zeros of (4) to a problem of finding zeros of a bivariate scalar-valued function (8). By the construction, the function $G_{k}$ vanishes at $(\nu,p_{1}(\nu))$, for $\nu\in\mathcal{J}$ — see (5). Therefore, we can factorize it in the following way

[TABLE]

where

[TABLE]

The function $g_{k}$ is called the reduced bifurcation function. Recall, that in most cases the function $\nu\to p_{1}(\nu)$ is unknown and we have only a rigorous bound on it and its derivatives. Similarly as $p_{1}$, the function $g_{k}$ cannot be computed exactly. However, it possible to bound values and partial derivatives of $g_{k}$ using integral representation (10).

In the second step, we apply the method from A to prove, that the set of zeroes of $g_{k}$ can be parametrized as a smooth curve $p_{1}\to\nu(p_{1})$. To sum up, in addition to C1, we make the following standing assumptions, which are all checkable be means of rigorous numerical methods:

• C2:

$P_{1}$ is a closed interval and $P_{2}$ is closed set, such that $p_{1}(\nu)\in\mathrm{int}\,P_{1}$ for $\nu\in\mathcal{J}$ and the solution set to (7) in $\mathcal{J}\times P_{1}\times P_{2}$ is a graph of smooth function $p_{2}\colon\mathcal{J}\times P_{1}\to\mathrm{int}\,P_{2}$;

• C3:

there exists a smooth function $P_{1}\ni p_{1}\to\nu(p_{1})\in\mathcal{J}$, such that

[TABLE]

where $g_{k}$ is defined by (10);

• C4:

there holds

[TABLE]

We will finish this section by showing, that assumptions C1–C3 imply, that the two curves, which solve (4) intersect at some point, while C4 implies that this intersection is unique.

Lemma 6

Let $s:A\to B$ and $t:B\to A$ be continuous, where $A$ is an open interval, $B$ is a closed interval and $s(A)\subset\mathrm{int}\,B$. Then there exists $(a^{*},b^{*})\in\mathrm{int}\left(A\times B\right)$, such that $(a^{*},s(a^{*}))=(t(b^{*}),b^{*})$.

Proof: Since $s(A)\subset\mathrm{int}B$, the set $\left(A\times B\right)\setminus\{(a,s(a)):a\in A\}$ has two connected components containing the points $(t(\min{B}),\min{B})$ and $(t(\max{B}),\max{B})$, respectively. Since $t$ is continuous, there is $b^{*}\in\mathrm{int}\,B$ such that $(t(b^{*}),b^{*})=(a^{*},s(a^{*}))$ for some $a^{*}\in A$. \qed

Lemma 7

Under assumptions C1–C4, there is a unique point $(\nu^{*},p_{1}^{*})\in\mathcal{J}\times P_{1}$, such that $\nu^{*}=\nu(p_{1}^{*})$ and $p_{1}^{*}=p_{1}(\nu^{*})$.

Proof: From C2–C3 and Lemma 6 the curves $\mathcal{J}\ni\nu\to p_{1}(\nu)\in P_{1}$ and $P_{1}\ni p_{1}\to\nu(p_{1})\in\mathcal{J}$ intersect at some point $(\nu^{*},p_{1}^{*})$. Assume that $(\hat{\nu},p_{1}(\hat{\nu}))$ is another intersection point for some $\hat{\nu}\in\mathcal{J}$. For $\tau\in[0,1]$ we define $w(\tau)=\nu^{*}+\tau(\hat{\nu}-\nu^{*})$ and

[TABLE]

From C3 we have $u(0)=u(1)=0$ and $0=u(1)-u(0)=S(\hat{\nu}-\nu^{*})$, where

[TABLE]

We will argue, that $S\in\frac{\partial^{2}G_{k}}{\partial\nu\partial p_{1}}(\mathcal{J}\times P_{1})+\frac{\partial^{2}G_{k}}{\partial p_{1}^{2}}(\mathcal{J}\times P_{1})p_{1}^{\prime}(\mathcal{J})$ and thus by C4 there must be $\hat{\nu}=\nu^{*}$. For $(\nu,p_{1})\in\mathcal{J}\times P_{1}$ we have

[TABLE]

In particular, for fixed $\tau\in[0,1]$ and $(\nu,p_{1})=u(\tau)$, both second partial derivatives of $G$ are evaluated at a point which does not depend on $t$. Therefore

[TABLE]

Similarly

[TABLE]

Using (11)–(13) we obtain

[TABLE]

\qed

3.2 Validation of period-tupling and touch-and-go bifurcations

Assumptions C1–C4 are easily checkable by means of rigorous numerics. The curve of $R$-symmetric period-$2k$ points is defined on an explicit domain $P_{1}$, which makes it possible for further continuation of this branch by the standard methods [9, 7, 20, 51] (see also  A). On the other hand, Lemma 7 guarantees, that there is an unique bifurcation point in the domain under consideration. The type of bifurcation, however, is unknown.

In what follows, we derive conditions, which along with C1–C4 guarantee, that period tupling or touch-and-go bifurcation occurs in the sense of Definition 3 and Definition 4, respectively.

Theorem 8

Under assumptions C1–C4, if $k>1$ is odd and

[TABLE]

then $f_{\nu}^{2}$ has $k^{\mathrm{th}}$-order touch-and-go bifurcation at some point $(\nu^{*},x^{*})\in\mathcal{J}\times P_{1}\times P_{2}\times\{0\}$.

Proof: From Lemma 7 there is a unique intersection point $(\nu^{*},p_{1}^{*})$ of two curves of zeroes of $G_{k}$ in $\mathcal{J}\times P_{1}$. From (14) it follows that

[TABLE]

Differentiating the identity $g_{k}(\nu(p_{1}),p_{1})\equiv 0$ we obtain

[TABLE]

From the above and (15) we conclude, that $\nu^{\prime}(p_{1}^{*})\neq 0$. Therefore, there is an interval $[\nu_{1},\nu_{2}]\subset\mathcal{J}$, such that $\nu^{*}=\nu(p_{1}^{*})\in(\nu_{1},\nu_{2})$ and the inverse function $\widetilde{x}_{b}=\nu^{-1}$ is defined on it. Since $p_{1}^{*}\in\mathrm{int}\,P_{1}$, the interval $[\nu_{1},\nu_{2}]$ can be chosen so that $\widetilde{x}_{b}([\nu_{1},\nu_{2}])\subset\mathrm{int}\,P_{1}$.

Define $U_{\varepsilon}=P_{1}\times P_{2}\times(-\varepsilon,\varepsilon)^{n}$. We assumed, that $\mathrm{Fix}(R)$ is a submanifold given locally by $q=0$. Therefore, there exists sufficiently small $\varepsilon>0$, such that for $\nu\in[\nu_{1},\nu_{2}]$ we have

[TABLE]

where

[TABLE]

From (6) it follows that

[TABLE]

for $i=1,\ldots,k-1$. Thus, all requirements from Definition 4 are satisfied on the set $[\nu_{1},\nu_{2}]\times U_{\varepsilon}$ and for the point $x^{*}=(p_{1}^{*},p_{2}(\nu^{*},p_{1}^{*}),0)$. \qed

The next theorem provides a framework for computer-assisted verification of the presence of period $k$-tupling bifurcation.

Theorem 9

Assume C1–C4. If $k\geq 2$ is even and

[TABLE]

then $f_{\nu}^{2}$ has period $k$-tupling bifurcation at some point $(\nu^{*},x^{*})\in\mathcal{J}\times P_{1}\times P_{2}\times\{0\}$.

Proof: From Lemma 7 we know, that the curves $\mathcal{J}\ni\nu\to p_{1}(\nu)\in P_{1}$ and $P_{1}\ni p_{1}\to\nu(p_{1})\in\mathcal{J}$ intersect at exactly one point $(\nu^{*},p_{1}^{*})$. We will prove, that $p_{1}^{*}$ is a proper local minimum of $p_{1}\to\nu(p_{1})$. First we will show, that $\nu^{\prime}(p^{*}_{1})=0$.

Define $x^{*}=(p_{1}^{*},p_{2}(\nu^{*},p_{1}^{*}),0)$. Let us choose a sequence $\{p_{1}^{n}\}_{n\in\mathbb{N}}$, such that $p_{1}^{n}\neq p_{1}^{*}$ and $\lim_{n\to\infty}p_{1}^{n}=p_{1}^{*}$. Define $\nu^{n}=\nu(p_{1}^{n})$ and $x^{n}=(p_{1}^{n},p_{2}(\nu^{n},p_{1}^{n}),0)$. Since $p_{1}^{n}\neq p_{1}^{*}$ and the intersection point is unique it follows, that $x^{n}$ is period-$2k$ point. Since $p_{1}^{*}\in\mathrm{int}\,P_{1}$, $k$ is even and $x^{*}$ is a fixed point for $f_{\nu^{*}}^{2}$ we conclude, that for $n$ large enough the point $y^{n}=f_{\nu^{n}}^{k}(x^{n})$ satisfies $\pi_{p_{1}}y^{n}\in\mathrm{int}\,P_{1}$ and

[TABLE]

The point $y^{n}$ is an $R$-symmetric, $2k$-periodic point for $f_{\nu^{n}}$, so it solves (4). Therefore $y^{n}\in P_{1}\times P_{2}\times\{0\}$. Since $x^{n}\neq y^{n}$, assumption C2 implies, that also $p_{1}^{n}\neq\pi_{p_{1}}y^{n}$. By the Rolle’s Theorem there is a point in the interval joining $p_{1}^{n}$ and $\pi_{p_{1}}y^{n}$ at which $\nu^{\prime}$ vanishes. In consequence, arbitrary close to $p_{1}^{*}$ there is a point at which $\nu^{\prime}$ is zero, which by the smoothness of $\nu^{\prime}(p_{1})$ implies $\nu^{\prime}(p_{1}^{*})=0$.

We will show that $\nu^{\prime\prime}(p_{1}^{*})>0$. Differentiating $g_{k}(\nu(p_{1}),p_{1})\equiv 0$ we obtain

[TABLE]

On the other hand

[TABLE]

Therefore

[TABLE]

We also have

[TABLE]

Differentiating twice $g_{k}(\nu(p_{1}),p_{1})\equiv 0$ with respect to $p_{1}$, we obtain

[TABLE]

Using $\nu^{\prime}(p_{1}^{*})=0$ and (16)–(19) we conclude, that

[TABLE]

This proves, that the function $p_{1}\to\nu(p_{1})$ has a proper local minimum at $p_{1}^{*}$.

Now, we will construct the functions $x_{b_{1}},x_{b_{2}}$ as required by Definition 3. Since $\nu^{\prime\prime}(p_{1}^{*})>0$, we can shrink $P_{1}$ in such a way, that $p_{1}^{*}\in\mathrm{int}\,P_{1}$, $\nu(\min P_{1})=\nu(\max P_{1})$ and $\nu$ is strictly convex on $P_{1}$. Given that $p_{1}(\nu)$ is smooth (so it has bounded slope near $\nu^{*}$) and $\nu^{\prime}(p_{1}^{*})=0$ we conclude, that for $\nu\in[\nu^{*},\nu(\min P_{1})]$ there holds $p_{1}(\nu)\in\mathrm{int}\,P_{1}$. Shrinking $\mathcal{J}$ from below, if necessary, we may also assume that $p_{1}(\nu)\in\mathrm{int}\,P_{1}$ for $\nu<\nu^{*}$ and $\nu\in\mathcal{J}$.

Let us fix $\nu_{2}\in(\nu^{*},\nu(\min P_{1}))$. The function $\nu(p_{1})$ is monotone on $[\min P_{1},p_{1}^{*}]$ and $[p_{1}^{*},\max P_{1}]$. Therefore, for $\widetilde{\nu}\in[\nu^{*},\nu_{2}]$

[TABLE]

and the functions $\widetilde{x}_{b_{i}}$ can be chosen to be continuous and smooth in $(\nu^{*},\nu_{2})$. Define $U_{\varepsilon}=P_{1}\times P_{2}\times(-\varepsilon,\varepsilon)^{n}$. For some sufficiently small $\varepsilon>0$ the functions

[TABLE]

given by

[TABLE]

satisfy conditions from Definition 3. From (6) it follows that

[TABLE]

for $i=1,\ldots,k-1$, which completes the proof. \qed

3.3 Symmetry breaking bifurcations

A special type of bifurcation covered by Definition 3 is when $k=1$. In the case of reversible maps, we are looking for the set of fixed points of $f_{\nu}^{2}$, which near the bifurcation point is expected to be a union of two intersecting curves. Thus, the principal branch of fixed points of $f^{2}_{\nu}$ cannot be isolated by applying Interval Newton Operator (A) to (4). Therefore, we need an extra constraint in order to isolate two curves that intersect at the bifurcation point. A possible scenario for this kind of bifurcation is breaking some symmetry $S$, which provides desired additional constraint $\mathcal{C}_{S}$, as defined by (1).

Let $f\colon M\subset\mathcal{J}\times X\to X$ be a family of $R$–reversible and $S$–symmetric maps. As in Section 3.1, we focus on computation of the set of $R$–symmetric fixed points for $f^{2}_{\nu}$.

Let $(\hat{\nu},\hat{x})\in M$ be an apparent bifurcation point and let $U\subset X$ be an open set, such that $\hat{x}\in U$ and $\mathcal{J}\times U\subset M$. We assume, that there are local coordinates in $U$, such that

[TABLE]

We impose that $\mathrm{Fix}(R,U)$ and $\mathrm{Fix}(S,U)$ intersect transversally. In the above settings, we expect that for $\nu\in\mathcal{J}$ the set of double-symmetric fixed points of $f_{\nu}^{2}$ in $U$ is a single point

[TABLE]

Thus, using the method described in A, it is possible to isolate the curve

[TABLE]

from the branching-off curves $x_{b_{1,2}}(\nu)\in\mathrm{Fix}(R,U)$, which intersect $x_{fp}(\nu)$ at exactly one point located in $\mathrm{Fix}(R,U)\cap\mathrm{Fix}(S,U)$.

The remaining steps in validation of the existence of a symmetry breaking bifurcation go as described in Section 3.1 and Section 3.2. We perform the Lyapunov-Schmidt reduction (7) and we define the bifurcation functions $G_{1}$ and $g_{1}$ by (8) and (10), respectively.

In Theorem 9 we assumed that the period of branching-off orbits is even. This lead us to a conclusion, that after half of iterations, the points sufficiently close to the bifurcation point come back to its vicinity. This allowed us to prove, that the reduced bifurcation function $g_{k}$ has a local minimum at the bifurcation point, which was the crucial step in the proof of the presence of period-tupling bifurcation.

In the case of the symmetry breaking bifurcation we cannot use this argument, because $k=1$ is odd. The next theorem says, that commutativity of $R$ and $S$ yields to the same conclusion.

Theorem 10

Let $S$ and $R$ be a symmetry and a reversing symmetry for $f_{\nu}$, respectively. Assume that C1–C4 are satisfied for $k=1$ and with $x_{fp}$ given by (20). If the symmetries $R$ and $S$ commute and (16) is satisfied, then $f_{\nu}^{2}$ has symmetry breaking bifurcation at some point $(\nu^{*},x^{*})\in\mathcal{J}\times P_{1}\times P_{2}\times\{0\}$.

Proof: From Lemma 7 there is an unique $(\nu^{*},p_{1}^{*})\in\mathrm{int}\,\mathcal{J}\times P_{1}$, such that $p_{1}(\nu^{*})=p_{1}^{*}$ and $\nu(p_{1}^{*})=\nu^{*}$. Let us set $x^{*}:=(p_{1}(\nu^{*}),p_{2}(\nu^{*}),0)\in\mathrm{Fix}(S)\cap\mathrm{Fix}(R)$.

We will show, that $\nu^{\prime}(p_{1}^{*})=0$. Take a sequence $\left\{p_{1}^{n}\right\}_{n\in\mathbb{N}}$ of points from $P_{1}$, which converges to $p_{1}^{*}$ and such that $p_{1}^{n}\neq p_{1}^{*}$ for $n\in\mathbb{N}$. Let us define $\nu^{n}=\nu(p_{1}^{n})$, $x^{n}=(p_{1}^{n},p_{2}(\nu^{n},p_{1}^{n}),0)\in\mathrm{Fix}(R)$ and $y^{n}=S(x^{n})$. By the commutativity of $R$ and $S$ we have

[TABLE]

hence $y^{n}\in\mathrm{Fix}(R)$. Similarly $f_{\nu^{n}}(y^{n})\in\mathrm{Fix}(R)$, because

[TABLE]

This shows, that $y^{n}$ solves (4) with $\nu=\nu^{n}$ and $\nu(\pi_{p_{1}}(y^{n}))=\nu^{n}$. Since

[TABLE]

we conclude, that $y^{n}\in P_{1}\times P_{2}\times\{0\}$ for $n$ large enough.

Since $p_{1}^{n}\neq p_{1}^{*}$ and the intersection point is unique, it follows that $x^{n}\notin\mathrm{Fix}(S)$, so $x^{n}\neq y^{n}$. This and C2 imply that also $p_{1}^{n}\neq\pi_{p_{1}}y^{n}$. By the Rolle’s Theorem there is a point in the interval joining $p_{1}^{n}$ and $\pi_{p_{1}}y^{n}$ at which $\nu^{\prime}$ vanishes. Given that both $p_{1}^{n}$ and $\pi_{p_{1}}(y^{n})$ (see (21)) converge to $p_{1}^{*}$ we conclude, that arbitrarily close to $p_{1}^{*}$ there is a point at which $\nu^{\prime}$ is zero and thus $\nu^{\prime}(p_{1}^{*})=0$.

The remaining part of the proof goes as in Theorem 9. Condition (16) implies, that $\nu^{\prime\prime}(p_{1}^{*})>0$, which proves that $p_{1}^{*}$ is a proper minimum of $\nu(p_{1})$. This makes it possible to define two branches $x_{b_{1}}$ and $x_{b_{2}}$, as required in Definition 3. \qed

3.4 General unfolding of period-tupling and touch-and-go bifurcations

Theorem 8 and Theorem 9 provide frameworks for computer-assisted verification, that the two types of bifurcations occurs in the sense of geometric conditions given by Definition 4 and Definition 3, respectively. In this section we will show, that the same geometric assumptions lead to standard unfolding of these bifurcations.

Theorem 11

Under assumptions of Theorem 9 (C1–C4 and (16)), there is a smooth $\nu$-dependent substitution $p_{1}=p_{1}(\nu,w)$, which brings the bifurcation function (8) to the normal form

[TABLE]

where $\alpha=\alpha(\nu)$ does not vanish at $\nu=0$ and $\beta=\beta(\nu)$ is positive at $\nu=0$.

Proof: Let $\nu^{*}$ be the parameter value at which bifurcation occurs. Without loosing generality we may assume, that $\nu^{*}=0$. From the proof of Theorem 9 we have that $\nu^{\prime}(p_{1}^{*})=0$. Differentiation of the identity $g_{k}(\nu(p_{1}),p_{1})\equiv 0$ gives

[TABLE]

First, a $\nu$-dependent substitution $p_{1}=z+p_{1}(\nu)$ brings $G_{k}$ to the form

[TABLE]

Put $z^{*}=0$. Differentiation of the product $G_{k}(\nu,z)=zg_{k}(\nu,z)$ gives

[TABLE]

because $g_{k}(\nu^{*},z^{*})=0$ and from (23) also $\frac{\partial g_{k}}{\partial z}(\nu^{*},z^{*})=\frac{\partial g_{k}}{\partial p_{1}}(\nu^{*},p_{1}^{*})=0$. The non-degeneracy condition (16) and C4 imply that

[TABLE]

Using (24)–(28) we can write $G_{k}$ in the form

[TABLE]

where $c_{i}$, $i=1,2,3$ are $\nu$-dependent coefficients with $c_{1}(0)\neq 0$ and $c_{3}(0)\neq 0$. Consider a $\nu$-dependent substitution of the form

[TABLE]

Then

[TABLE]

In order to annihilate the term $\nu w^{2}$ we should take $h=c_{2}/c_{1}$, which is well defined near $\nu^{*}=0$, because $c_{1}(0)\neq 0$. With such choice of $h$, we have

[TABLE]

where $d_{3}=c_{3}-2c_{2}^{2}\nu/c_{1}$ and $d_{3}(0)=c_{3}(0)\neq 0$. Put $\beta=-d_{3}/c_{1}$. From (16) and (28) we have $c_{1}(0)c_{3}(0)<0$. Therefore $\beta(0)=-c_{3}(0)/c_{1}(0)>0$. Setting $\alpha=c_{1}$ we conclude, that $G_{k}$ takes the required form (22). \qed

Theorem 12

Under assumptions of Theorem 8 (C1–C4 and (14)), there is a smooth $\nu$-dependent substitution $p_{1}=p_{1}(\nu,w)$ which brings the bifurcation function (8) to the normal form

[TABLE]

where $\alpha=\alpha(\nu)$, $\beta=\beta(\nu)$ and they are non-zero at $\nu=0$.

Proof: Let $\nu^{*}$ be the parameter value at which bifurcation occurs. Without loosing generality we may assume that $\nu^{*}=0$. Reasoning as in the proof of Theorem 11, substitution $z=p_{1}-p_{1}(\nu)$ brings $G_{k}$ to the form $G_{k}(\nu,z)=zg_{k}(\nu,z)$. Put $z^{*}=0$. From C1–C4 and (14) we have that

[TABLE]

Hence, in these coordinates we have

[TABLE]

with $c_{1}(0)\neq 0$ and $c_{2}(0)\neq 0$. Setting $\alpha=c_{1}$ and $\beta=c_{2}/c_{1}$ we obtain the required normal form (29). \qed

4 Validation of bifurcations of symmetric periodic orbits: reversible Hamiltonian systems

In this section we show, how to adopt the general framework described in Section 3 to autonomous Hamiltonian systems. We consider an $R$–reversible Hamiltonian system

[TABLE]

where $J$ is the standard symplectic matrix and $H\colon\mathbb{R}^{2n+2}\to\mathbb{R}$ is $\mathcal{C}^{4}$-smooth. First, we choose a $\mathcal{C}^{3}$-smooth Poincaré section $\Pi\subset\mathbb{R}^{2n+2}$. If $\mathrm{Fix}(R)\subset\Pi$, by [48, Lemma 3.3] the Poincaré map $\mathcal{P}\colon\Pi\to\Pi$ is $R|_{\Pi}$–reversible, too. In what follows we will study bifurcations of $R$-symmetric fixed points for $\mathcal{P}^{2}$.

In autonomous Hamiltonian systems, a natural choice of the bifurcation parameter is the value of $H$, because it is a constant of motion. On the other hand, from the point of view of rigorous numerics, it is much easier and efficient to work in the phase space coordinates and parametrize periodic orbits for $\mathcal{P}$ by one of them.

This dissonance between numerical efficiency and formal description of a bifurcation is solved in the following way. We expect that two families of periodic points for $\mathcal{P}$ intersect at a bifurcation point. This geometric condition can be checked in the phase-space coordinates. Additional conditions on $H$, which will be given in this section, will guarantee, that $H$ can be used (locally) as a bifurcation parameter.

Let us give a brief overview of the construction we are going to perform. Period-tupling and touch-and-go bifurcations are local phenomena. Thus, we can choose local coordinates $(p_{0},p,q_{0},q)\in\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{n}$ near an apparent bifurcation point, such that

$\Pi=\left\{(p_{0},p,q_{0},q)\in\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{n}:q_{0}=0\right\}$ and 2. 2.

$\mathrm{Fix}(R)=\{(p_{0},p,q_{0},q)\in\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{n}:q_{0}=0,q=0\}\subset\Pi$.

In order to simplify further notation we will use $(p_{0},p,q)$ coordinates in $\Pi$ and we will always skip $q_{0}=0$ as an argument of $\mathcal{P}$ and $H$. We also assume, that the vector field (31) is transverse to $\Pi$ near an apparent bifurcation point, so that the Poincaré map is well defined and smooth. Define $\Pi_{h}=\left\{(p_{0},p,q)\in\Pi:H(p_{0},p,q)=h\right\}$.

Assuming that near an apparent bifurcation point there holds

[TABLE]

we can conclude that the projection $\Pi_{h}\ni(p_{0},p,q)\to(p,q)$ is a local diffeomorphism. This allows us to parametrize $\Pi_{h}$ locally by $(p_{0}(h,p,q),p,q)$ for $h$ belonging to some interval $\mathcal{H}$. Our goal is to specify conditions on $\mathcal{P}$ and $H$, which will guarantee, that the map

[TABLE]

is well defined on some domain and it has period-tupling and/or touch-and-go bifurcations in the sense of Definition 3 and Definition 4, respectively.

Let us fix $k>1$. Following Section 3, we split $p=(p_{1},p_{2})$ and $q=(q_{1},q_{2})$. We assume, that the set of $R$–symmetric fixed points of $\mathcal{P}^{2}$ near an apparent bifurcation point forms a regular curve, which is parametrized by the coordinate $p_{0}$

[TABLE]

and defined on an explicit, open interval $P_{0}$. We also assume, that for $i=1,\ldots,k-1$ there holds

[TABLE]

First we perform the Lyapunov-Schmidt reduction. We assume that there is a set $P_{0}\times P_{1}\times P_{2}$ and a smooth function $p_{2}^{H}:P_{0}\times P_{1}\to P_{2}$, such that for $(p_{0},p_{1},p_{2})\in P_{0}\times P_{1}\times P_{2}$

[TABLE]

Using this implicit function we can define the bifurcation function by

[TABLE]

and factorize it as $G_{k}^{H}(p_{0},p_{1})=(p_{1}-p_{1}(p_{0}))g_{k}^{H}(p_{0},p_{1})$, where the reduced bifurcation function reads

[TABLE]

Now we can specify the standing assumptions in the context of Hamiltonian systems:

• HC2:

$P_{1}$ is a closed interval, and $P_{2}$ is a closed set, such that $p_{1}(p_{0})\in\mathrm{int}\,P_{1}$ for $p_{0}\in P_{0}$, where $p_{1}(p_{0})$ is given by (33) and there is a smooth function $p_{2}^{H}:P_{0}\times P_{1}\to\mathrm{int}\,P_{2}$ solving (35);

• HC3:

there exists a smooth function $P_{1}\ni p_{1}\to p_{0}(p_{1})\in P_{0}$ such that

[TABLE]

where $g_{k}^{H}$ is defined by (37);

• HC4:

there holds

[TABLE]

Assumptions HC2–HC4 and Lemma 7 imply, that the set of $R$–symmetric fixed points for $\mathcal{P}^{2k}$ in $P_{0}\times P_{1}\times P_{2}\times\{0\}$ is the union of two regular curves, which intersect at unique point $x^{*}=(p_{0}^{*},p_{1}^{*},p_{2}^{H}(p_{0}^{*},p_{1}^{*}),0)$, where $p_{2}^{H}$ is defined by (35). These parametric curves are defined on explicit range $p_{0}\in P_{0}$ and $p_{1}\in P_{1}$, respectively, which makes it possible for further continuation of these branches by the method described in A. On the other hand, we have no information about the type of bifurcation at the intersection point.

Generically we expect, that the value of $H$ can be used as the bifurcation parameter. It may happen, however, that $H$ is not monotone (or even constant) along one or both curves of periodic points. In the remaining part of this section we derive (easily checkable by means of rigorous numerics) conditions, which guarantee, that

$H$ can be used as the parameter near $x^{*}$ and 2. 2.

the mapping (32) undergoes one of the bifurcations defined by Definition 3 and/or Definition 4.

For further use we define four functions

[TABLE]

Lemma 13

If $\mathcal{P}^{2}(x^{*})=x^{*}$ and

[TABLE]

then there is an open interval $\mathcal{H}$ and an open set $U$, such that

$s(x^{*})\in\mathcal{H}\times U$, 2. 2.

$s^{-1}|_{\mathcal{H}\times U}$* is well defined diffeomorphism onto image and* 3. 3.

$\mathcal{P}^{k}\circ s^{-1}$* is defined on $\mathcal{H}\times U$.*

Proof: From (38) it follows, that $\det\left(Ds(x^{*})\right)=\frac{\partial H}{\partial p_{0}}(x^{*})\neq 0$ and thus $s$ is a local diffeomorphism near $x^{*}$. The inverse function takes the form

[TABLE]

Since $x^{*}$ is a fixed point of $\mathcal{P}^{2}$, we can also find open sets $\mathcal{H}$ and $U$, such that $s(x^{*})\in\mathcal{H}\times U$ and such that $\mathcal{P}^{k}\circ s^{-1}$ is defined on $\mathcal{H}\times U$. \qed

Theorem 14

Assume HC2–HC4 and let $x^{*}=(p_{0}^{*},p_{1}^{*},p_{2}^{H}(p_{0}^{*},p_{1}^{*}),0)$ be the unique intersection point from Lemma 7. If $k>1$ is odd, $\frac{\partial H}{\partial p_{0}}(x^{*})\neq 0$, $h_{fp}^{\prime}(p_{0}^{*})\neq 0$ and $h_{b}^{\prime}(p_{1}^{*})\neq 0$, then $f^{2}_{h}$ defined by (32) has $k^{\mathrm{th}}$-order touch-and-go bifurcation at $s(x^{*})$.

Proof: Take $\mathcal{H}\times U$ from Lemma 13. From (33)–(34) the point $x_{b}^{H}(p_{1})$, $p_{1}\in P_{1}$ is of principal period $2k$ for $\mathcal{P}$, except the intersection $p_{1}^{*}$. Since $h_{fp}^{\prime}(p_{0}^{*})\neq 0$ and $h_{b}^{\prime}(p_{1}^{*})\neq 0$, we can choose $[h_{1},h_{2}]\subset\mathcal{H}$, $h^{*}\in(h_{1},h_{2})$, such that $h_{fp}^{-1}$ and $h_{b}^{-1}$ are defined on $[h_{1},h_{2}]$ and both functions

[TABLE]

have range in $U$. By the construction they satisfy conditions from Definition 4 for the function $f$ defined by (32). \qed

Theorem 15

Assume HC2–HC4 and let $x^{*}=(p_{0}^{*},p_{1}^{*},p_{2}^{H}(p_{0}^{*},p_{1}^{*}),0)$ be the unique intersection point from Lemma 7. If $k\geq 2$ is even, $\frac{\partial H}{\partial p_{0}}(x^{*})\neq 0$, $h_{fp}^{\prime}(p_{0}^{*})>0$ and $h_{b}^{\prime\prime}(p_{1}^{*})>0$, then $f^{2}_{h}$ defined by (32) has period $k$-tupling bifurcation at $s(x^{*})$.

Proof: Take $\mathcal{H}\times U$ from Lemma 13. We will show that $h_{b}^{\prime}(p_{1}^{*})=0$. Let us fix a sequence $p_{1}^{n}\to p_{1}^{*}$, such that $p_{1}^{*}\neq p_{1}^{n}$ and $x_{b}^{H}(p_{1}^{n})$ is defined for all $n\in\mathbb{N}$. Let us set $x^{n}=x_{b}^{H}(p_{1}^{n})$ and $y^{n}=\mathcal{P}^{k}(x^{n})$. We also have $y^{n}\neq x^{n}$, because by (34) the principal period of $x^{n}$ is $2k$. This and HC3 imply, that $p_{1}^{n}\neq\pi_{p_{1}}y^{n}$. Since $k$ is even and $x^{*}$ is a fixed point for $\mathcal{P}^{2}$, we have that $\lim_{n\to\infty}y^{n}=x^{*}$, $y^{n}\neq x^{n}$ and $H(y^{n})=H(x^{n})=h_{b}(p_{1}^{n})$. Hence, in the interval joining $p_{1}^{n}$ and $\pi_{p_{1}}y^{n}$ there is a point at which $h_{b}^{\prime}$ vanishes and in consequence $h_{b}^{\prime}(p_{1}^{*})=0$. The assumption $h_{b}^{\prime\prime}(p_{1}^{*})>0$ implies, that $p_{1}^{*}$ is a proper local minimum of $h_{b}$.

The remaining part of the construction goes as in Theorem 9. The function $h_{b}$ is convex near $p_{1}^{*}\in\mathrm{int}\,P_{1}$ which allows us to define two continuous branches

[TABLE]

of $h_{b}^{-1}$, for some $h_{2}\in\mathcal{H}$. The function $h_{fp}$ is invertible near $p_{0}^{*}$, because we assumed that $h_{fp}^{\prime}(p_{0}^{*})\neq 0$. Restricting the domain, if necessary, we may assume, that the inverse is defined on $[h_{1},h_{2}]$ with $h^{*}\in(h_{1},h_{2})$. Define

[TABLE]

Shrinking again the interval $[h_{1},h_{2}]$ if necessary, we may assume, that the range of the above three functions is in $U$. Thus $f_{h}^{2}$ defined by (32) has period $k$-tupling bifurcation at $s(x^{*})$. \qed

We end this section by a version of Theorem 10 for Hamiltonian systems. Proceeding as in Section 3.3 we assume, that (31) admits a symmetry $S$ and a reversing symmetry $R$. We also assume that the function defined by (33) satisfies

[TABLE]

Then, we define $G_{1}^{H}$ and $g_{1}^{H}$ by (36) and (37), respectively. We state the following result without a proof, as it is similar to that of Theorem 10 and Theorem 15.

Theorem 16

Assume HC2–HC4 with $k=1$ and let $x^{*}=(p_{0}^{*},p_{1}^{*},p_{2}^{H}(p_{0}^{*},p_{1}^{*}),0)$ be the unique intersection point from Lemma 7. If $S$ and $R$ commute, $\frac{\partial H}{\partial p_{0}}(x^{*})\neq 0$, $h_{fp}^{\prime}(p_{0}^{*})>0$ and $h_{b}^{\prime\prime}(p_{1}^{*})>0$, then $f^{2}_{h}$ defined by (32) has symmetry breaking bifurcation at $s(x^{*})$.

5 Bifurcations of odd periodic solutions in the Falkner-Skan equation

The Falkner-Skan equation [35] is a third order ODE given by

[TABLE]

It is low-dimensional and without singularities. Therefore it is relatively easy for rigorous numerical investigation. Although for some physical reasons solutions of certain BVP for (39) are relevant, we use the system to test the methodology introduced in Section 3 and thus we will focus on periodic solutions. Using this example we will also provide the reader with some details regarding validation of the presence of period-tupling and touch-and-go bifurcations. A higher-dimensional hamiltonian system (Circular Restricted Three Body Problem), which is also more computationally demanding, will be studied in Section 6.

Our aim is to prove that some family of odd periodic solutions of (39) parametrized by $\beta$ undergoes period-doubling, third order touch-and-go and period-quadrupling bifurcations.

The equation (39) can be rewritten as a system of first order equations

[TABLE]

and in the sequel we will work with (40). The system (40) is reversible with respect to $R(x,y,z)=(-x,y,-z)$. For all $\beta>0$ the points $(0,\pm 1,0)$ are $R$-symmetric equilibria. Numerical simulation shows, that there is a family of $R$-symmetric periodic orbits $u_{\beta}(t)=(x_{\beta}(t),y_{\beta}(t),z_{\beta}(t))$ parametrized by $\beta>1$. These periodic orbits intersect the symmetry line $\mathrm{Fix}(R)=\{(0,y,0):y\in\mathbb{R}\}$ at exactly two points, which approach $(0,\pm 1,0)$, respectively, when $\beta\to 1^{+}$ — see Fig. 3. Thus, the period of these orbits goes to infinity when $\beta\to 1^{+}$. It is also observed, that $\max_{t\in\mathbb{R}}\|x_{\beta}(t)\|$ goes to infinity when $\beta\to 1^{+}$.

For small $\beta>1$, these periodic orbits are hyperbolic. For larger $\beta$ they become elliptic and crossing strong $1:k$ resonances, $k=2,3,4$ at approximate parameter values $\hat{\beta}_{k}$, respectively, where

[TABLE]

Approximate initial conditions for resonant $R$-symmetric periodic orbits are $(0,\hat{y}_{k},0)$, where

[TABLE]

and their shapes are shown in Fig. 4. It seems, the family $u_{\beta}$ continues to exists for all $\beta>1$ approaching $1:1$ resonance when $\beta\to\infty$.

Let us define a Poincaré section

[TABLE]

We will use $(x,y)$ coordinates to describe points in $\Pi$. By $\mathcal{P}_{\beta}:\Pi\to\Pi$ we denote the associated Poincaré map for the system (40) with the parameter value $\beta$. The map $\mathcal{P}_{\beta}$ is reversible with respect to involution $R(x,y)=(-x,y)$. We will also use the notation $\mathcal{P}(\beta,x,y)=\mathcal{P}_{\beta}(x,y)$.

The aim of this section is to give a computer-assisted proof of the following theorem.

Theorem 17

There is a smooth family $u_{\beta}=(x_{\beta},y_{\beta})$ of $R$-symmetric period-two points for $\mathcal{P}_{\beta}$ parametrized by $\beta\in\mathcal{J}=[\frac{9}{8},100\,000]$. This family undergoes period-doubling, third order touch-and-go and period-quadrupling bifurcations at some points $(\beta_{k}^{*},0,y_{k}^{*})$, $k=2,3,4$, respectively, with

[TABLE]

where approximate bifurcations points are listed in (41)–(42).

Proof: The existence of a smooth branch of $R$-symmetric period-two points for $\mathcal{P}_{\beta}$ has been proved by means of the method described in A. We used an adaptive cover of the parameter range $\mathcal{J}=\bigcup_{j=1}^{10478}\mathcal{J}_{j}$, where the diameters of intervals $\mathcal{J}_{j}$ are smaller (approximately $2\cdot 10^{-5}$) if $\mathcal{J}_{j}$ is close to $\frac{9}{8}$ and quite large (above $141$) in the second end of the parameter range. Then, for each subinterval $\mathcal{J}_{j}$ we check the assumptions of the parametrized Interval Newton Method (Lemma 25). If succeed, we prove that the segments $y(\mathcal{J}_{j})$ glue into a smooth curve by checking the the bounds on $y(\mathcal{J}_{j})$ resulted from the Interval Newton operator overlap, when their domains do.

We will give more details regarding validation of bifurcations. The bifurcation function is

[TABLE]

for $k=2,3,4$. Hence, the Lyapunov-Schmidt reduction is not needed. In order to apply the general framework introduced in Section 3 we need first to check conditions C2–C4. We have the following bounds on $G_{1}$.

[TABLE]

We see that in each case $\hat{y}_{k}-G_{1}(\mathcal{J}_{k}\times\{\hat{y}_{k}\})\cdot\frac{\partial G_{1}}{\partial y}(\mathcal{J}_{k}\times Y_{k})^{-1}\subset Y_{k}$, which proves that, there is a branch $(\beta,y(\beta))$ of zeroes of $G_{1}$ parametrized by $\beta\in\mathcal{J}_{k}$, for $k=2,3,4$.

In order to check C3 we need bounds on the reduced bifurcation function $g_{k}$ – see (10). We have

[TABLE]

Note, that in order to obtain tiny bounds on $g_{k}(\{\hat{\beta}_{k}\}\times Y_{k})$ we used high precision interval arithmetic [36]. Again, in each case we have $\hat{\beta}_{k}-g_{k}(\{\hat{\beta}_{k}\}\times Y_{k})\cdot\frac{\partial g_{k}}{\partial\beta}(\mathcal{J}_{k}\times Y_{k})^{-1}\subset\mathcal{J}_{k}$, which proves that for $k=2,3,4$ the function $g_{k}$ has branch of zeroes $(\beta(y),y)$ parametrized by $y\in Y_{k}$ and thus C3 is satisfied.

From the following estimates

[TABLE]

it follows that

[TABLE]

for $k=2,3,4$ and the condition C4 is also satisfied.

There remains to check conditions specific for each type of bifurcation. In (43) we have already computed bound on $\frac{\partial^{2}G_{2}}{\partial y^{2}}(\mathcal{J}_{2}\times Y_{2})$, which proves that the assumptions of Theorem 8 are satisfied. Thus, the proof of the existence of third order touch-and-go bifurcation for $\mathcal{P}^{2}$ in $\mathcal{J}_{k}\times\{0\}\times Y_{k}$ is completed.

For period-doubling and period-quadrupling bifurcations we have to check non-degeneracy condition (16). From (43) we already have, that for $k=2,4$ there holds $\frac{\partial^{2}G_{k}}{\partial\beta\partial y}(\mathcal{J}_{k}\times Y_{k})>0$. Thus, it suffices to check that $\frac{\partial^{3}G_{k}}{\partial y^{3}}(\mathcal{J}_{k}\times Y_{k})$ is non-zero. We have the following bounds

[TABLE]

Eventually, we have to check that the principal periods of bifurcating orbits. The cases $k=2,3$ do not require any computation, as both numbers are primes. For the case $k=4$ we computed

[TABLE]

This completes the proof. \qed

6 Halo orbits in the Circular Restricted Three Body Problem

In this section we apply the general framework described in Section 4 to the Circular Restricted Three Body Problem. First, we will give a short overview of the CR3BP and we list some of its relevant properties. Then, we will give a computer-assisted proof, that the well known families of halo orbits undergo various types of bifurcations.

6.1 Equations of motion

The CR3BP is a mathematical model, that describes the motion of a small body with negligible mass under the gravitational influence of two point like big bodies, called primaries, which rotate around their common mass centre on a circle.

Denote by $\mu$ the relative mass of the smaller primary. In a rotating coordinate system centred at the common mass centre of two big primaries, the dynamics of the small particle is governed by the following system of second-order differential equations [26, 44]

[TABLE]

where

[TABLE]

The system is Hamiltonian and it admits a first integral, called the Jacobi constant, which is given by

[TABLE]

The hyperplane $\{(x,y,z=0,\dot{x},\dot{y},\dot{z}=0)\}$ is invariant under the local flow induced by (44) and the corresponding four-dimensional Hamiltonian system is called the Planar Circular Restricted Three Body Problem (PCR3BP).

6.2 Symmetries of the CR3BP

The system possesses two main symmetries

[TABLE]

It is important to note that $R$ and $S$ commute. This property is required for our method of validation of the existence of symmetry breaking bifurcations — see Theorem 18.

6.3 Poincaré map in the CR3BP

Let us define the following Poincaré section

[TABLE]

and denote by $\mathcal{P}_{\mu}:\Pi\to\Pi$ the associated Poincaré map. We will skip the dependency on $\mu$ and write $\mathcal{P}$, if it will be clear from the context. Since the $y$ variable is fixed and equal to zero on the section, we will use $(x,z,\dot{x},\dot{y},\dot{z})$ coordinates to describe points in $\Pi$. With some abuse of notation on $R$, we will denote by the same letter the reversing symmetry of the system restricted to points on $\Pi$, i.e.

[TABLE]

Since

[TABLE]

by [48, Lemma 3.3] the mapping $\mathcal{P}$ is reversible, too. Thus, the frameworks for computer-assisted verification of various types of bifurcations introduced in Section 4 can be applied to $\mathcal{P}$.

6.4 Periodic orbits near libration points

The CR3BP possesses five equilibrium points, called the libration points. All of them are located in the $\{z=0,\dot{z}=0\}$ invariant hyperplane and thus they are equilibrium points for the PCR3BP, as well. Three of libration points, commonly denoted by $L_{1}$, $L_{2}$ and $L_{3}$, are collinear and are located on the $x$-axis. They are of saddle-centre type for the PCR3BP. It is well known [44], that for all $\mu\in(0,1)$ there exists a family of $R$–symmetric periodic orbits, called planar Lyapunov orbits (PLO), that surround these libration points — see also Figure 5. In [11, 12] the existence of Lyapunov orbits around $L_{1}$ and $L_{2}$ libration points for selected mass parameters has been proved in an explicit domain. Computer-assisted methods have been used to prove [2, 45, 49, 50], that for the mass parameter $\mu=0.0009537$ corresponding to the Sun-Jupiter system there are Lyapunov orbits around $L_{1}$ and $L_{2}$ for certain energy level, and there is countable infinity of connecting orbits between them in both directions.

For the full system (CR3BP) the libration points $L_{1,2,3}$ are of saddle-centre-centre type, for all $\mu$. In additional direction $z$ there exists a second family of vertical Lyapunov orbits (VLO), which are double symmetric both with respect to symmetry $S$ and reversing symmetry $R$ defined by (45). These orbits intersect twice the $x$-axis. Therefore, an object (a spacecraft) located nearby one of those orbits will be periodically collinear with the two main primaries. Thus eclipses or shadows are unavoidable for trajectories approaching planar or vertical Lyapunov orbits.

There is a numerical evidence [19, 21, 22], that a branch of out-of-plane $R$-symmetric orbits, called halo orbits, bifurcates from the Lyapunov family. In opposite to planar and vertical Lyapunov orbits, they never cross the $x$-axis, except at the bifurcation point. For large vertical amplitudes $z$ the halo orbits become more and more aligned to the $(y,z)$ plane allowing continuous observation of both primaries without eclipses. Parts of $L_{1}$-Lyapunov and $L_{1}$-halo families are shown in Figure 5 for the relative mass corresponding to the Sun-Jupiter system.

6.5 Symmetry breaking bifurcations of halo orbits in the CR3BP

The halo orbits in the CR3BP are $R$-symmetric, out-of-plane periodic orbits, which bifurcate from the planar Lyapunov family — see Figure 5. Although, there were extensive numerical study of these orbits and their bifurcations (just to mention few papers [22, 21, 19]), to the best of our knowledge, they were never proved to exist.

The best theoretical result in this direction has been done in [14, 15]. The authors consider a certain normal form, which approximates the CR3BP. They prove, that in the normal form system there is a branch of $R$-symmetric halo orbits bifurcating from $R$-symmetric planar Lyapunov orbits. Moreover, they give an explicit expression for the bifurcation point. This result is valid for all $L_{1,2,3}$ points and for all mass parameters.

The results of this section are complementary to those from [14]. We give a computer-assisted proof of the existence of halo orbits for the original CR3BP system for all libration points $L_{1,2,3}$ and for some (not full) range of $\mu$. We will also study continuation and bifurcations of $L_{1,2}$-halo families.

Theorem 18

There are continuous functions

[TABLE]

with $[\mu_{*},\mu^{*}]=[9.5\cdot 10^{-4},0.5]$ and $Z=[-1,1]\cdot\Delta_{z}$, where $\Delta_{z}=10^{-7}$, such that

$h_{i}(\mu,0)$* is a point of symmetry breaking bifurcation of out-of-plane family of $R$–symmetric periodic orbits from the planar family of Lyapunov orbits near $L_{i}$ libration point and* 2. 2.

for $\mu\in[\mu_{*},\mu^{*}]$ and $z\in Z$ the point $h_{i}(\mu,z)$ is an initial condition for an $R$–symmetric periodic out-of-plane (halo) orbit.

Proof: We applied the framework from Section 4 to the family of Poincaré maps $\mathcal{P}_{\mu}$. For fixed $i=1,2,3$ we proceed as follows (we skip dependencies on $i$ to simplify notation). The range of parameters $[\mu_{*},\mu^{*}]$ is initially subdivided into $N$ smaller overlapping subintervals $[\mu_{*},\mu^{*}]=\bigcup_{j=1}^{N}\mu_{j}$, where $\mu_{j}=[\underline{\mu_{j}},\overline{\mu_{j}}]$. For each $j=1,\ldots,N$ we execute in parallel the following (independent) tasks.

We find an approximate bifurcation point $\hat{u}_{j}=(\hat{x},0,0,\hat{\dot{y}},0)$ for the mass parameter $\hat{\mu}_{j}=\frac{1}{2}(\underline{\mu_{j}}+\overline{\mu_{j}})$. For this purpose we use the scheme described in B — see also Remark 26. 2. 2.

The planar double-symmetric Lyapunov orbits can be easily isolated by restriction to the planar system. Using Lyapunov-Schmidt reduction and the method from A we validate the existence of smooth, two-parameter families of periodic orbits

[TABLE]

which correspond to Lyapunov orbits and halo orbits, respectively. The diameters of $X$ and $\dot{Y}$ were hand-optimized to speed-up computation. 3. 3.

The set $W=X\times\{0\}\times\{0\}\times\dot{Y}\times\{0\}$ is a bound for the bifurcation point for each $\mu\in\mu_{j}$. Then we check the assumptions of Theorem 16, i.e. $\frac{d}{dx}C_{\mu_{j}}(W)\neq 0$, $\frac{d}{dx}C_{\mu_{j}}(u_{fp}(\mu_{j},X))\neq 0$ and $\frac{d^{2}}{dz^{2}}C_{\mu_{j}}h_{i}(\mu_{j},z=0)\neq 0$. Note, that in Theorem 16 we required, that the branching-off curve is convex, but switching of either sign does not change the geometry of overall picture.

If any of the above steps fails, the interval $\mu_{j}$ is being subdivided and we repeat computation for each element of subdivision. Finally, if all tasks are completed, using the methods from [20, 7] we check, that the pieces of $h_{i}$ glue into a smooth function. \qed

Graphs of $\mu\to h_{i}(\mu,0)$ are shown in Figure 6. We would like to emphasize, that they match quite well bifurcation points coming from the normal forms found in [14]. Our validation algorithm used to prove Theorem 18, by its construction, cannot continue with $\mu\to 0$. The threshold value $\mu_{*}=9.5\cdot 10^{-4}$ as well as size of out of plane amplitude $\Delta_{z}=10^{-7}$ are by our choice a compromise between CPU time needed to obtain the result and the range of $\mu$ and $z$ we can cover. In particular, the range of mass parameter $[\mu_{*},\mu^{*}]$ contains two relevant values

[TABLE]

corresponding to Sun-Jupiter and Earth-Moon systems, respectively. The values listed in (46) are taken as the nearest IEEE-754 double precision numbers to the recent mass measurements reported in [40, 3].

Remark 19

The $L_{3}$ is a special case, because for small $\mu$ the maximal order of normal form constructed in [14] is $2$, and it is divergent when $\mu\to 0$ [42]. In this case we observe an interesting phenomenon. If $\mu\to 0$, the Jacobi integral seems to be not monotone along the curve of bifurcation points — see Figure 6 right-bottom panel. The computation, which indicate the presence of a local minimum is non-rigorous but performed in high accuracy (400 bits of mantissa) floating point arithmetic [36] and using $80^{\mathrm{th}}$ order ODE solver with the tolerance per time step set to $10^{-60}$.

6.6 Continuation and bifurcations of halo orbits

Theorem 18 guarantees that for $\mu\in[9.5\cdot 10^{-4},0.5]$ a branch of halo orbits near $L_{i}$ can be parametrized by out of plane amplitude $z$ to at least $|z|\leq\Delta_{z}$. Numerical simulations [19, 21] strongly indicate, that these branches continue to exist for much larger amplitudes and that they undergo period-doubling and period-quadrupling and third-order touch-and-go bifurcations. The next theorem addresses this issue.

Theorem 20

Consider the CR3BP with $\mu\in\{\mu_{SJ},\mu_{EM}\}$ as defined in (46). There exists a smooth function $h_{\mu}:\mathcal{S}^{1}\to\mathbb{R}^{6}$ such that for $\tau\in[0,2\pi]$ the following holds true.

$h_{\mu}(\tau)=(x(\tau),0,z(\tau),0,\dot{y}(\tau),0)$* is an initial condition for an $R$–symmetric periodic (halo) orbit.* 2. 2.

$h_{\mu}(\tau)=S(h(2\pi-\tau))$* – the family is $S$ symmetric.*

This closed loop of halo orbits intersect the invariant subspace $\{z=0,\dot{z}=0\}$ at exactly two points $h(0)$ and $h(\pi)$, at which an symmetry breaking bifurcation occurs. Moreover, the branch $h_{\mu}$ undergoes period doubling, period quadrupling and third order touch-and-go bifurcations as listed in Table 1.

Proof: The branch of halo orbits is split into four pieces.

Proceeding as in the proof of Theorem 18 we validate the existence of two symmetry breaking bifurcations: one in the region $x>0$ and the second in $x<0$ — see Figure 7. From HC3 we also have, that there is a branch of out-of-plane $R$-symmetric periodic orbits parametrized by $z\in[-\Delta_{z},\Delta_{z}]$, for an explicit $\Delta_{z}>0$. 2. 2.

Both out-of-plane families are then continued (see A) and parametrized by $z$ variable until hand-chosen threshold value $z=0.625$, as shown in Figure 7. 3. 3.

The upper arc joining two halo orbits with $z=0.625$ is parametrized by $x$ variable.

By the well known techniques [7] we can check, that these pieces glue into a smooth curve. By the symmetry we obtain the lower branch of halo orbits. Summarizing, we obtained, that the branch of halo orbits is a compact, smooth, one-dimensional manifold without boundary. It is well known [33], that such a manifold is diffeomorphic to a circle. By the construction the manifold is $S$-symmetric, thus the parametrization $h$ can be chosen to preserve this symmetry, as well.

Approximate bifurcation points listed in Table1 (except for $k=1$ were found with very high accuracy (of order $10^{-60}$) using high-order ODE solvers from the CAPD library [10] based on high precision floating-point arithmetic [36]. Then we checked assumptions of Theorem 14 and Theorem 15 on very small sets (of the size about $10^{-30}$) centred at these approximate bifurcation points. Notice, that in one case $j=7$ and $\mu=\mu_{SJ}$ we could not obtain bounds on third order derivatives sharp enough to check convexity of Jacobi integral along bifurcation curve. We checked, however, the conditions HC2–HC4, which in particular means, that there are two curves of halo orbits of period $1$ and $4$ intersecting at a single point. \qed

Remark 21

Although we did not prove it, there is a strong numerical evidence that the points $h_{\mu}(0)$ belongs to the family of $L_{1}$-Lyapunov orbits. A possible method to close this gap is to apply the method for computation of invariant manifolds of Lyapunov orbits, as proposed in [12].

Projections of the two curves $h_{\mu}(\mathcal{S}^{1})$ for $\mu\in\{\mu_{SJ},\mu_{EM}\}$ resulting from Theorem 20 onto $(x,z)$ plane are shown in Figure 7. These families form a $2D$–tori in the full phase space. Such torus for the mass $\mu=\mu_{SJ}$ is shown in Figure 8.

Numerical simulation [22, 19] shows that $L_{2,3}$-halo families continue to exists until a collision with one of the primaries — see Figure 9. Hence, on the Poincaré section $\Pi$, $(x,z)$ coordinates of these orbits approach $(1-\mu,0)$ and $(\mu,0)$ respectively, while $|\dot{y}|$ tends to infinity.

We have the following partial result for the $L_{2}$-halo family.

Theorem 22

Consider the CR3BP with $\mu\in\{\mu_{EM},\mu_{SJ}\}$ as defined in (46). There is a smooth function $h_{\mu}:[-1,1]\to\mathbb{R}^{6}$ such that for $\tau\in[-1,1]$ the following statements hold true.

$h_{\mu}(\tau)=(x_{\mu}(\tau),0,z_{\mu}(\tau),0,\dot{y}_{\mu}(\tau),0)$* is an initial condition for an $R$–symmetric periodic (halo) orbit.* 2. 2.

$h_{\mu}(\tau)=S(h(-\tau))$* – the family is $S$ symmetric.* 3. 3.

$h_{\mu}(0)$* is a point of an symmetry breaking bifurcation.* 4. 4.

The function $\dot{y}_{\mu}(\tau)$ has a unique local minimum at $\tau=0$ satisfying

[TABLE] 5. 5.

The branches continue to at least

[TABLE]

Moreover, these branches undergo period doubling, period quadrupling and third order touch-and-go bifurcations as listed in Table 2.

Proof: The validation is split into three steps.

From Theorem 18 we know, that there is an symmetry breaking bifurcation of halo orbits from $L_{2}$-Lyapunov family. The estimates (47) are taken from this proof. The branching of family of halo orbits $h_{\mu}(z)=(x_{\mu}(z),0,z,0,\dot{y}_{\mu}(z),0)$ is parametrized by $z\in[-1,1]\cdot\Delta_{z}$, with $\Delta_{z}=10^{-7}$. We also check that $\dot{y}_{\mu}^{\prime\prime}(z)>0$ for $|z|\leq\Delta_{z}$. 2. 2.

The branch is then rigorously continued using A and parametrized by $h_{\mu}(z)=(x_{\mu}(z),0,z,0,\dot{y}_{\mu}(z),0)$, $z>0$ until some hand-chosen threshold value of $\hat{z}$ (dependent on $\mu$). We also checked that for $z\in[\Delta z,\hat{z}]$ there holds $\dot{y}_{\mu}^{\prime}(z)>0$. 3. 3.

Further continuation of the branch starting from $h_{\mu}(\hat{z})$ is parametrized by $\dot{y}$ until hand-chosen threshold values (48).

Summarizing, in each segment the variable $\dot{y}$ is increasing along the branch of periodic orbits which makes it possible to re-parametrize the curve as a function

[TABLE]

defined on $\tau\in[0,1]$. From the symmetry $S$ we obtain the second branch for $\tau\in[-1,0]$.

The bifurcations listed in Table 2 (except $j=0$) were validated using Theorem 15 and Theorem 14 in high-precision interval arithmetics. Notice, that in two cases $j=8$ and $\mu=\{\mu_{SJ},\mu_{EM}\}$ we could not obtain bounds on third order derivatives sharp enough to check convexity of Jacobi integral along bifurcation curve. We checked, however, the conditions HC2–HC4, which in particular means, that there are two curves of halo orbits of period $1$ and $4$ intersecting at a single point. \qed

Remark 23

The threshold values in (48) have been chosen so that the validated arc of $L_{2}$–halo orbits contains all strong resonances at which period-tupling and touch-and-go bifurcations occurs for both values of $\mu\in\{\mu_{SJ},\mu_{EM}\}$. Numerical simulation shows, that for larger values of $\dot{y}$ strong resonances are not present. However, we did not validate this conjecture. A possible approach to close this gap and obtain a full picture of what happens to $L_{2,3}$–halo families is to perform Levi-Civita regularisation [44] and continue branches of halo orbits in these coordinates.

6.7 Implementation notes

The source code of all programs is available to download from the web page of the corresponding author [47]. The programs are written in C++-11 and use rigorous ODE solvers and algorithms for computation of Poincaré maps and their derivatives [46, 53] from the CAPD library [10]. All programs the from the were compiled using g++-4.9.2 and executed on a computer equipped with Intel Xeon E7-8867 v4 2.40GHz processors (64 cores).

Appendix A Interval Newton method for implicit equations

In this section we provide a method for validated computation of implicit functions. Such an algorithm is needed to check assumptions C2–C3 or HC2–HC3, as well as to compute wide branches of periodic orbits far from bifurcation points. The method is an adaptation of the well known Interval Newton Method (INO) [34, 39] to the case of implicit equations. The main modification which significantly improves the method is the use of higher order derivatives. The method itself is quite straightforward but to the best of our knowledge, it has not appeared in the literature. First we recall the Interval Newton Method [34, Thm. 8.4], [39, Thm. 5.1.7].

Theorem 24** ([34, 39])**

Let $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be a $\mathcal{C}^{1}$ map and let $X\subset\mathbb{R}^{n}$ be a convex, compact set. For $x_{0}\in\mathrm{int}X$ we define the interval Newton operator by

[TABLE]

where by $[Df(X)]_{I}$ we mean a convex hull of the set of matrices $\{Df(x):x\in X\}$.

If $N(f,x_{0},X)\subset\mathrm{int}\,X$ then $f$ has a unique zero in $X$ that belongs to $N(f,x_{0},X)$. 2. 2.

If $N(f,x_{0},X)\cap X=\emptyset$ then $f$ has no zeros in $X$.

The above theorem can be used to validate the existence of solutions to implicit equations over an explicit domain.

Lemma 25

Let $f:\mathbb{R}^{m}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ be a $\mathcal{C}^{1}$ map, $Z\subset\mathbb{R}^{m}$ be the closure of an open set and let $X\subset\mathbb{R}^{n}$ be a convex, compact set with non-empty interior. For $x_{0}\in\mathrm{int}\,X$ we define the interval Newton operator by

[TABLE]

If $N(f,x_{0},X,Z)\subset\mathrm{int}\,X$ then there exists a $\mathcal{C}^{1}$ smooth function $g:Z\to X$ such that the set of zeroes $\{(z,x)\in Z\times X:f(z,x)=0\}$ coincides with the graph of the function $g$, i.e.

[TABLE]

Proof: Let us fix $z\in Z$ and put $f_{z}=f(z,\cdot)$. Then we have

[TABLE]

From Theorem 24 for all $z\in Z$ there exists a unique $x=g(z)\in X$ such that $f(z,x)=0$. Let us observe, that the condition (49) implies, that $D_{x}f(z,x)$ is invertible at each point $(z,x=g(z))$. By the implicit function theorem the function $g$ is smooth, as it solves an implicit equation $f(z,g(z))\equiv 0$ at every point $z\in Z$. \qed

The efficiency of the INO (Theorem 24) is hidden in the fact, that we usually have a very good approximation (from numerical experiments) for zero, i.e. $f(x_{0})\approx 0$. Then the quantity $[Df(X)]_{I}^{-1}f(x_{0})$ can be very tight, even if the computed bound on derivative $Df(X)$ is overestimated. This is not the case for parametrized maps, as equation (49) contains the term $f(Z,x_{0})$. If $Z$ is large then having $f(z,x)\approx 0$ for all $z\in Z$ is rather unlikely.

A straightforward way to overcome this problem is to make a substitution $(z,x)=s(z,w)$, such that in the new coordinates the function $z\to w(z)$, which solves the implicit equation $f(s(z,w))=0$, is flat. Let us fix $(z_{0},x_{0})\in Z\times X$, such that $f(z_{0},x_{0})\approx 0$ and assume $D_{x}f(z_{0},x_{0})$ is non-singular. We define an affine substitution by

[TABLE]

where $A=D_{x}f(z_{0},x_{0})^{-1}D_{z}f(z_{0},x_{0})$. The map $s$ is invertible because its linear part has determinant equal to one and thus zeroes of $f$ are in one-to-one correspondence with zeroes of $g:=f\circ s$. In the new coordinates the point $(z_{0},w_{0}=0)$ is an approximate zero of $g$. Moreover, $D_{z}g(z_{0},w_{0})=0$ (provided $A$ is computed exactly). The above idea of changing linearly coordinate system has been proposed in [7], but the method for validation of a branch of zeroes was different and based on so-called radii polynomial approach [20, 7].

In the remaining part of the section we will show, how to efficiently evaluate all the terms, which appear in the INO for the mapping $g$. Our test show, that using this approach we could significantly reduce overestimations in evaluation of the INO, which lead to significant advantage of this approach in comparison to direct evaluation of INO for $g$.

The INO for the mapping $g$ on the set $Z\times W$, which contains $(z_{0},w_{0})$ reads

[TABLE]

In what follows, we will show how we can bound all the terms that appear in this expression. Let $X$ be such that $s(Z,W)\subset(Z,X)$ and denote $\Delta Z=Z-z_{0}$. The term $g(Z,w_{0})$ can be bounded by means of the mean value theorem

[TABLE]

and the set of matrices $D_{z}g(Z,w_{0})$ can be bounded by

[TABLE]

By the choice of $A$, we have $D_{z}g(z_{0},w_{0})=0$. Therefore we expect that for not very large parameter radius $\Delta Z$, the term $D_{z}g(Z,w_{0})\cdot\Delta Z$ in (50) is a small box around zero, as desired. Note, that the above considerations hold true for any matrix $A$. Therefore the quantities $D_{x}f(z_{0},x_{0})^{-1}$ and $D_{z}f(z_{0},x_{0})$ can be computed using just floating point arithmetic, however their product $A$ must be bounded rigorously, if we want to take the intersection in (51).

The quantity $g(Z,w_{0})$ can be also bounded using second order Taylor expansion

[TABLE]

Recall, the point $(z_{0},x_{0})$ is chosen so that $f(z_{0},x_{0})\approx 0$ and the substitution $s$ is chosen so that $D_{z}g(z_{0},w_{0})\approx 0$ (note, usually it cannot be computed exactly). Since both quantities are evaluated at a single point, we can take advantage (if necessary) of high-precision interval arithmetic and make these terms as close to zero as desired. Therefore, the bound on $g(Z,w_{0})$ is practically quadratic in the radius of the parameter range $\Delta Z$. One can take the intersection of direct evaluation of $g(Z,w_{0})$ in interval arithmetic with the bounds obtained from (50) and (52).

In order to compute the Interval Newton Operator for $g$ we have to bound $D_{w}g(Z,W)$. Direct evaluation gives

[TABLE]

Using second order derivatives of $f$ we can obtain another enclosure

[TABLE]

which can be intersected with (53).

Numerical experiments we performed show that using the above approach we can usually validate the existence of solution to the implicit equation on much wider domain $Z$ (without subdivision) than in the original coordinates.

In principle, both $g(Z,w_{0})$ and $D_{w}g(Z,W)$ can be bounded using higher order Taylor expansions. This should come along with nonlinear (usually polynomial) substitution $s$, such that all derivatives of $g=f\circ s$ with respect to $z$ vanish at $(z_{0},w_{0})$ up to desired order $r$. Then the bound on $g(Z,w_{0})$ can be made of order $O(\|\Delta Z\|^{r+1})$.

In the context of the CR3BP we have found, that the second order expansion is very efficient. Note, that computation of higher order derivatives of Poincaré maps in a high dimensional system is costly — the complexity of the $\mathcal{C}^{r}$–Lohner algorithm [53] used to integrate variational equations is $O(n^{3}r^{n}s^{2})$, where $n$ is the dimension, $s$ is the order of Taylor method and $r\geq 1$ is the largest order of derivative of Poincaré map, we request. Clearly, increasing $r$ in a high-dimensional system is very expensive. Secondly, the bounds on higher order derivatives of Poincaré map are usually much overestimated than those of lower order.

Appendix B Newton like scheme for finding bifurcation points

A straightforward way to localize period-tupling and touch-and-go bifurcation points of a family of reversible maps $f_{\nu}(x)$ is to follow the branch of period-$2$ points $x(\nu)$ and look for resonant eigenvalues of $Df^{2}_{\nu}(x(\nu))$. In low-dimensional systems one can look at the stability parameter [21, 22]. This method is very efficient when we want to find a rough approximation to the bifurcation point. In reversible or hamiltonian case multiple eigenvalues may occur making computation of eigenvalues with high-accuracy quite non-trivial task.

In this short section we propose eigenvalue-independent yet efficient scheme for finding very accurate approximation to bifurcation points. In what follows we focus on reversible Hamiltonian systems and use the notation from Section 4, but the idea applies to any family of reversible maps.

Following Section 4 we assume that $(p_{0},p_{1},p_{2},q=0)\in\mathrm{Fix}(R)$ is an approximate period-$2$ point for the Poincaré map $\mathcal{P}$, which is close to $1:k$ resonance. We would like to refine it by a Newton-like scheme. Since we have $n+1$ unknowns $(p_{0},p_{1},p_{2})$ we need $n+1$ equations with expected isolated zero. We impose

[TABLE]

where $g_{k}^{H}$ is defined by (37). The first equation selects the points from the curve of fixed points, while the second equation guarantees that the solution is also on the bifurcation curve. In order to apply the Newton method to the system of equations (54), we need to compute derivatives of $g_{k}^{H}$. From (37) we have

[TABLE]

If the seed point for the Newton method is quite close to the bifurcation point, we may assume that $p_{1}(p_{0})$ almost constant and thus

[TABLE]

The second partial derivative reads

[TABLE]

Again, assuming $p_{1}\approx p_{1}(p_{0})$ we can approximate

[TABLE]

In order to compute $D^{2}G_{k}^{H}$ we need second order derivatives of $\mathcal{P}^{k}$ and of the function $p_{2}^{H}(p_{0},p_{1})$ obtained from the Lyapunov-Schmidt reduction — see (35). The latest can be computed by differentiation of the identity

[TABLE]

Summarizing, the Newton-like iteration for equation (54) is given by

[TABLE]

where $r$ is the solution to the linear equation $M\cdot r=b$ with

[TABLE]

Using the above scheme, finding approximate bifurcation points of halo orbits with accuracy $10^{-60}$ was not a difficult task.

Remark 26

In the computer-assisted proof of Theorem 18 we used similar strategy to localize approximate points of symmetry breaking bifurcations. We solved for zeroes of the function

[TABLE]