\section{Introduction} \label{sec: intro}

The spatial three body problem studies the motion of three point masses in $\mathbb{R}^3$ under Newton's universal gravitational law. When two of the masses are equal, there is an invariant subsystem, where there is an fixed axis in $\mathbb{R}^3$ (without loss of generality let's assume it is the $z$-axis), such that two bodies with equal masses are symmetric with respect to the axis all the time and the third body always moves on the axis. This subsystem is called the isosceles three body problem, as the three bodies always form an isosceles triangle.

\begin{figure}[h]
\centering
\includegraphics[width=0.7 \textwidth]{fig_of_isosceles_3BP2.pdf}
\caption{The isosceles three body problem.}
\label{fig:IsoThree}
\end{figure}

Let $m_i$ and $q_i$ be the mass and position of the $i$-th body, $i=1,2,3$. Set $m=m_1=m_2$, then $q_1$ and $q_2$ are symmetric with respect to the $z$-axis and $q_3$ belongs to it. Let $q_1=(x, y,z) \in \mathbb{R}^3$ and $\dot{q}_1=(\dot x, \dot y, \dot z)$, then $q_2=(-x, -y, z)$ and $\dot q_2=(-\dot x, -\dot y, \dot z)$. Since the center of mass is fixed at the origin, $q_3=(0, 0, -2 \alpha z)$ and $\dot q_3=(0, 0, - 2\alpha \dot z)$, where $\alpha = m/m_3$.

The system is entirely determined by the motion of the first body and can be described as an Lagrange system with the following Lagrangian
\begin{equation}
\label{eq_Lagxyz} \begin{aligned}
L(q, \dot q) & = K(\dot q) + U(q)  \\
             & = m \big( \dot x^2 + \dot y^2 + (1+2\alpha) \dot z^2 \big) \\
             &  + \frac{1}{2} m^2\left( \frac{1}{\sqrt{x^2 +y^2}} + \frac{4\alpha^{-1}}{\sqrt{x^2 +y^2 +(1+2\alpha)^2 z^2}} \right).
\end{aligned}
\end{equation}


Introduce the matrix $M=\text{diag}(2m, 2m, 2m(1+2\alpha))$ and the new variable
$$  \xi = \left( \begin{array}{c}
\xi_1 \\ \xi_2 \\ \xi_3
\end{array} \right)
= M^{\frac{1}{2}} \left( \begin{array}{c}
x \\ y \\ z
\end{array} \right) \; \text{ and } \; \dot{\xi} = \left( \begin{array}{c}
\dot \xi_1 \\ \dot \xi_2 \\ \dot \xi_3
\end{array} \right)
= M^{\frac{1}{2}} \left( \begin{array}{c}
\dot x \\ \dot y \\ \dot z
\end{array} \right).
$$
We transfer the previous Lagrangian to
\begin{equation}
\label{eq_LagXi} L(\xi, \dot \xi) = \frac{1}{2} |\dot \xi|^2 + \frac{m^{5/2}}{\sqrt{2}} \left( \frac{1}{\sqrt{\xi_1^2 + \xi_2^2}} + \frac{4 \alpha^{-1}}{\sqrt{\xi_1^2 +\xi_2^2+ (1+2\alpha)\xi_3^2}} \right)
\end{equation}


Due to the rotation symmetry with respect to the $z$-axis, it is more convenient to work with the cylindrical coordinates $(r, \theta,z) \in [0, \infty) \times [-\pi, \pi] \times \mathbb{R}$ defined by
$$ \xi_1= r\cos \theta, \; \xi_2 = r  \sin \theta, \; \xi_3 =z. $$
This further transfers the Lagrangian to
$$ L(r, \theta, z, \dot{r}, \dot{\theta}, \dot{z}) = \frac{1}{2} ( \dot r^2 +\dot z^2 + r^2 \dot \theta^2 ) + U(r,z), $$
where the potential function $U$ is independent of $\theta$
\begin{equation}
\label{eq_Urz} U(r,z) = \frac{1}{\sqrt{2}} m^{5/2} \left( \frac{1}{r} + \frac{4 \alpha^{-1}}{\sqrt{r^2 + (1+2\alpha)z^2 }} \right).
\end{equation}
By the Legendre transformation
\begin{equation}
\label{eq_Legendre} p_r = \dot r, \; p_z = \dot{z} , \; p_{\theta} = r^2 \dot \theta,\end{equation}
we get the corresponding Hamiltonian
\begin{equation}
\label{eq_Ham}  H(p_r, p_z, p_{\theta}, r, z, \theta) = \frac{1}{2} \left(p_r^2 + p_z^2 + \frac{p_{\theta}^2}{r^2 } \right)- U(r,z).
\end{equation}


Recall that the angular momentum of the three body problem is a three dimensional vector and does not change along a given solution. Here it is always parallel to the $z$-axis. As a result, we will treat it as a scalar function and denote it by $\varpi$. By a straight forward computation (up to multiplication of a constant)
\begin{equation}
\label{eq: AngMom} \varpi = p_{\theta} = r^2 \dot \theta.
\end{equation}

Since the Hamiltonian $H$ is independent of $\theta$, for a fixed angular momentum $\varpi$, we get a reduced Hamiltonian with only two degrees of freedom
\begin{equation}
\label{eq_HamOm} H_{\varpi}(\zeta) = \frac{1}{2} \left(p_r^2 + p_z^2  \right) + V_{\varpi}(r, z), \;\; \zeta = (p_r, p_z, r, z) \in \mathbb{R}^2 \times \mathbb{R}_+ \times \mathbb{R},
\end{equation}
where $V_{\varpi}$ can be seen as the potential energy of the reduced Hamiltonian
\begin{equation}
\label{eq_Vrphi} V_{\varpi}(r, z) = \frac{1}{2} \frac{\varpi^2}{r^2 } - U(r,z).
\end{equation}

\begin{figure}[H]
\centering
\includegraphics[width=0.7\textwidth]{fig_of_V2.pdf}
\caption{ The level sets of $V_{\varpi}$ with $\varpi=\sqrt{5/3}, \alpha=6$.}
\label{picture of V}
\end{figure}


The corresponding reduced Hamiltonian system is
\begin{equation}
\label{eq: RedHamSys} \dot{\zeta} = X_{H_{\varpi}}(\zeta)= J \nabla H_{\varpi}(\zeta),
\end{equation}
where $X_{H_{\varpi}}$ is the associated Hamiltonian vector field, $J$ is the standard symplectic matrix and $\nabla H_{\varpi}(\zeta) = (p_r, p_z, -\frac{\varpi^2}{r^3} - \partial_r U, -\partial_z U)$ with
$$ -\partial_r U(r,z) = \frac{m^{5/2}}{\sqrt{2}}(r^{-2}+\frac{4r\alpha^{-1}}{(r^2+(1+2\alpha)z^2)^{3/2}}); $$

$$ -\partial_z U(r,z) = \frac{m^{5/2}}{\sqrt{2}} \frac{4\alpha^{-1}(1+2\alpha)z}{(r^2+(1+2\alpha)z^2)^{3/2}}. $$


Given a solution $\zeta(t)$ of \eqref{eq: RedHamSys}, by the angular momentum identity \eqref{eq: AngMom}, we can find the corresponding $p_{\theta}(t)=\varpi$ and $\theta(t)=\int^t_0\frac{\varpi}{r^2(t)}dt+\theta(0)$. Together with $\zeta(t)$, they form a solution of the original Hamiltonian system as well as the isosceles three body problem.



The isosceles three body problem plays an important role in understanding the overall dynamics of the three body problem and has been studied by many authors. The restricted case with the third body having zero mass is usually known as the Sitnikov problem, where oscillatory motions were first found and proven by Sitnikov \cite{Sit60} and Alekseev \cite{Alk72} (an alternative proof can also be found in the celebrated book \cite{Ms73} by Moser). For the unrestricted case, such a result was first proven by Moeckel \cite{Mk84}. The behaviors of the bodies near triple collision were studied in \cite{Dev80}, \cite{Mk81} and \cite{SM87} using McGehee coordinates, and played a key role in Xia's construction of the first example of non-collision singularity in the five body problem \cite{Xia92}. In this paper, we will focus on (relative) periodic orbits of this problem.

In recently years, many new (relative) periodic orbits of the $n$-body problem have been found  using action minimization methods, see \cite{CM00}, \cite{Yu21} and the references within. For the isosceles three body problem, some results including numerical ones using these methods were obtained in \cite{Shi09} and \cite{Yan15b}. However up to our knowledge the action minimization methods have only been successfully applied to the fixed-period problem and in general one can not control the energy (except when it is non-negative) or the angular moment of the obtained orbit.


To study periodic orbits at a given energy surface, Poincar\'e introduced the concept of an annulus-like global surface of section in \cite{Pc12} and stated a theorem about the existence of periodic orbits without giving a proof (this is usually known as \emph{Poincar\'{e}'s last geometric theorem}). Shortly afterwards it was proven by Birkhoff \cite{Bk13}. Since then many generalizations of this result have been obtained, see \cite{Fr90,Fr92,Calvez06} and it is impossible to give a complete list of references. To apply these results one needs to find a global surface of section first, which usually is not an easy task.  For convex energy surfaces of a Hamiltonian system with two degrees of freedom, a celebrated result regarding this was obtained by Hofer, Wysocki and Zehnder \cite{HWZ98}.  Since then, there have been many interesting works in the dynamics of three dimensional contact manifolds, please refer \cite{SH18, HMS15} and the references within.


In celestial mechanics, up to our knowledge such results seem to be only available for the restricted three body problem. First by Conley \cite{Con63} and McGehee \cite{McG69} using perturbation methods and more recently in \cite{AFFHO12,MvK} using holomorphic curve techniques from \cite{HWZ98}. The main purpose of this paper is to show that for the isosceles three body problem there are some natural global surfaces of section for certain choices of energy and angular momentum. Then by investigating the Poincar\'e return maps on these global surfaces of section, we are able to prove the existence of all kinds of different periodic solutions.

To state our results in more details, let's fix an energy $h$ of the reduced Hamiltonian and denote the corresponding energy surface as
$$ \mathfrak{M}(h, \varpi)= \{(p_r, p_z, r, z) \in \mathbb{R}^2 \times \mathbb{R}_+ \times \mathbb{R}: \; H_{\varpi} = h \}. $$
Our goal is to understand the dynamics on the energy surface, which not only depends on $h$ and $\varpi$, but also the masses. Notice that it is the mass ratio $\alpha$ that really affects the dynamics of the system, so for simplicity, in the rest of the paper, we assmue $$m=2^{1/5}.$$

The following result can be found in a paper by Moeckel \cite[Proposition 2.1]{Mk84}.
\begin{prop}
\label{prop_Mcomp} For $h < 0$ and $\varpi \in \mathbb{R}$ satisfying
\beg{eqnarray} \label{energy sphere} 1 < 2 \varpi^2 |h| < \beta^{-2}, \;\text{ where } \beta = ( 1 + 4 \alpha^{-1})^{-1},
\end{eqnarray}
$\mathfrak{M}(h, \varpi)$ is a compact manifold homeomorphic to $S^3$.

Moreover, $\mathfrak{M}(h,\varpi)$ is empty if $|h|>(2\varpi^2 \beta^{2})^{-1}$, $\mathfrak{M}(h,\varpi)=\{(0,0,\varpi^2\beta,0)\}$ if $|h|=(2\varpi^2 \beta^{2})^{-1}$ and $\mathfrak{M}(h,\varpi)$ is unbounded in $z$ direction if $|h|<(2\varpi^2 \beta^{2})^{-1}$.
\end{prop}



\begin{rem} Although $\mathfrak{M}(h, \varpi)$ is convex for certain choices of $(h,\varpi, \alpha)$, it is difficult to show $\mathfrak{M}(h,\varpi)$ is a contact manifold. If this is true, one may use powerful tools from symplectic dynamics \cite{BH12}. As a result, we use more direct approaches in our proofs without assuming the contactness of $\mathfrak{M}(h,\varpi)$.
\end{rem}


\begin{rem} \label{rem:TwoParameters}
In the rest of the paper, we always assume $h<0$ with $(h, \varpi)$ satisfying \eqref{energy sphere}. As a result, $\varpi \ne 0$. Since the results for $\varpi>0$ and $\varpi <0$ are essentially the same, we will further assume $\varpi >0$.

Although the dynamics on the energy surface $\mathfrak{M}(h, \varpi)$ is determined by $\alpha, h$ and $\varpi$. Because of the homogeneity of the Newtonian potential, $h$ and $\varpi$ is essentially one parameter given as $2h \varpi^2$ (see the hand written note by Moeckel \cite{Mk94}).
\end{rem}





Consider the following two closed subsets of the energy surface $\mathfrak{M}(h, \varpi)$:
\begin{equation}
\label{eq Sigma pm}
\begin{aligned}
\Sigma_- & =\{(p_r,p_z,r,z)\in \mathfrak{M}(h, \varpi): \; z=0,\ p_z \le 0 \},\\
\Sigma_+ & =\{(p_r,p_z,r,z)\in \mathfrak{M}(h, \varpi): \; z=0,\ p_z \ge 0 \}.
\end{aligned}
\end{equation}
We will show they are two disk-like global surfaces of section (see Definition \ref{defi: GlobalSurface}).

\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{fig_of_Sigma2.pdf}
\caption{  $\Sigma_\pm$ with $\beta = \mathfrak{e}=0.6$.}
\label{picture of Sigma}
\end{figure}




\begin{thm} \label{thm: GlobalSurface}
 $\Sigma_{\pm}$ are two disk-like global surfaces of section of the Hamiltonian flow $\varphi_t$ corresponding to the reduced Hamiltonian $H_{\varpi}$.
\end{thm}

Notice that $\Sigma_{\pm}$ shares the same boundary. In Subsection \ref{sub:Euler}, we will show it is nothing but the Euler orbit, which will be denoted as $\zeta_e(t)$. For such an Euler orbit the trajectories of the masses are ellipses with the same eccentricity $\mathfrak{e}$ satisfying
$$ \mathfrak{e} = (1 + 2 h \varpi^2\beta^{2})^{1/2}. $$
Because of this, by Remark \ref{rem:TwoParameters}, the dynamics on $\mathfrak{M}(h, \varpi)$ is determined by $\beta$ and $\mathfrak{e}$, and by Proposition \ref{prop_Mcomp}, we will only consider parameters from the following region
\begin{eqnarray}\label{equ. the region of (beta,e)}
\mathcal{D}:=\{(\beta,\mathfrak{e})| 0<\mathfrak{e}<(1-\beta^2)^{\frac{1}{2}}, 0<\beta<1\},
\end{eqnarray}
for which the corresponding energy surface $\mathfrak{M}(h, \varpi)$ is homeomorphic to $S^3$.


Since $H_{\varpi}$ is a \emph{mechanical Hamiltonian}, i.e., it is a summation of the kinetic and potential energy, $\zeta(t) = (p_r, p_z, r, z)(t)$ is an $h$-energy solution of \eqref{eq: RedHamSys} if and only if $x(t) = (r, z)(t)$ is a $h$-energy solution of the following Euler-Lagrange system,
\begin{equation}
\label{eq: SecOrdSys} \ddot{x}= -\nabla V_{\varpi}(x); \quad \frac{1}{2} |\dot{x}|^2 + V_{\varpi}(x) = h.
\end{equation}

Let $\Pi$ denote the projection of $\mathfrak{M}(h, \varpi)$ to the configuration space $(r, z)$. Then
\begin{equation*}
    \mathcal{H}(h) = \Pi(\mathfrak{M}(h, \varpi)) =\{ (r, z): V_{\varpi}(r, z) \le h \}, \;\;
 \end{equation*}
$$ \mathcal{B}(h) = \partial \mathcal{H}(h) =\{ (r, z): V_{\varpi}(r, z)=h \},$$
are the corresponding \emph{Hill's region} and \emph{Hill's boundary}. Notice that for any $(r, z) \in \mathcal{B}(h)$, there is a unique point
$\Pi^{-1}(r, z) = (0,0, r,z)  \in \mathfrak{M}(h, \varpi)$. For any subset $\mathcal{A} \subset \mathcal{B}(h)$, in the following we use $\mathcal{A}$ to denote both $\mathcal{A}$ and $\Pi^{-1}(\mathcal{A})$, when there is no confusion.


	Since $\mathcal{B}(h)$ is the zero velocity curve of $h$-energy orbits, if an $h$-energy orbit $\zeta(t)$ satisfies $\zeta(t_0) \in \mathcal{B}(h)$, for some $t_0 \in \mathbb{R}$, then
\begin{equation}
\label{eq: brake}  \zeta(t_0 + \delta t) = \begin{pmatrix}
-I & 0\\
 0 & I
\end{pmatrix} \zeta(t_0 - \delta t) , \;\; \forall \delta t \in \mathbb{R}.
\end{equation}


Moreover if there is another $t_1>t_0$ (smallest after $t_0$) also satisfying $\zeta(t_1) \in \mathcal{B}(h)$, then $\zeta(t)$ must be a $2(t_1 -t_0)$-periodic orbit traveling back and forth between $\zeta(t_0)$ and $\zeta(t_1)$ along the same curve. This will be called a \emph{brake orbit} as the velocities are zero at $t=t_0$ and $t=t_1$.

\begin{rem}
The study of brake orbits goes back to the work of Seifert \cite{Sft48}. For recent advances we refer the readers to \cite{LLZ21} and the references within.
\end{rem}

To distinguish different types of brake orbits, we set
$$ \mathcal{B}_+(h) = \mathcal{B}(h) \cap \{(r, z): z \ge 0 \}, \;\; \mathcal{B}_-(h) = \mathcal{B}(h) \cap \{(r, z): z \le 0 \}. $$

\begin{defi}
\label{defi: BrakeOrbits} Let $\zeta(t)$ be a brake orbit with $\zeta(t_0), \zeta(t_1) \in \mathcal{B}(h)$ and $\zeta(t) \notin \mathcal{B}(h)$, for all $t \in (t_0, t_1)$.
\begin{enumerate}
\item We say $\zeta(t)$ is \textbf{a type-I brake orbit}, if $\zeta(t_0) \in \mathcal{B}_+(h)$ and $\zeta(t_1) \in \mathcal{B}_-(h)$, or $\zeta(t_0) \in \mathcal{B}_-(h)$ and $\zeta(t_1) \in \mathcal{B}_+(h)$.
\item We say $\zeta(t)$ is \textbf{a type-II brake orbit}, if $\zeta(t_0) \in \mathcal{B}_+(h)$ and $\zeta(t_1) \in \mathcal{B}_+(h)$, or $\zeta(t_0) \in \mathcal{B}_-(h)$ and $\zeta(t_1) \in \mathcal{B}_-(h)$.
\end{enumerate}
\end{defi}
\begin{rem}
Notice that a brake orbit can only be one of the above two types, except the Euler orbit, which is both.
\end{rem}

\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{fig_symmetric_orbits2.pdf}
 \caption{Different types of periodic orbits for $\beta = \mathfrak{e}=0.6$.}
 \label{picture of symmetric orbits}
\end{figure}

Since the potential energy $V_{\varpi}$ satisfies the following $z$-symmetry property,
\begin{equation}
\label{eq:Vz-sym} V_{\varpi}(r, z) = V_{\varpi}(r, -z),
\end{equation}
if an orbit $\zeta(t)$ satisfies
\begin{equation*} \zeta(t_0) \in \mathfrak{M}(h, \varpi) \cap \{z=0, p_r=0\}, \; \text{ for some } t_0 \in \mathbb{R},
\end{equation*}

	we say it is \emph{$z$-symmetric}, as \eqref{eq:Vz-sym} implies,
\begin{equation} \label{eq: z-sym}
\zeta(t_0 + \delta t) = \begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(t_0 - \delta t), \;\; \forall \delta t \in \mathbb{R}, \; \text{ where } A = \begin{pmatrix}
-1 & 0 \\
0 & 1
\end{pmatrix}.
\end{equation}



A brake orbit which is also $z$-symmetric will be called a \emph{$z$-symmetric brake orbit}. We point out that every $z$-symmetric brake orbit (except the Euler orbit) must be type-I (see Proposition \ref{prop:z-symm brake}). 



\begin{thm}\label{th: z symm brake}
There exists a $z$-symmetric brake orbit $\zeta_{bz}(t)= (p_r^*, p_z^*, r^*, z^*)(t)$ satisfying $\zeta_{bz}(0) \in \mathcal{B}_+(h)$ and
$$ p^*_r(t^*)=0, \; z^*(t^*) =0, \;\text{ and } \; z^*(t)>0, \; \forall t \in [0, t^*), \; \text{ for some } t^* >0.$$
Moreover $\zeta_{bz}(t)$ forms a hopf link with the Euler orbit inside $\mathfrak{M}(h, \varpi)$.
\end{thm}



By Theorem \ref{thm: GlobalSurface}, there is a Poincar\'e return map $g: \Sigma_-^o \to \Sigma_-^o$ given by the Hamiltonian flow $\varphi_t$. In Section \ref{sec:GlobalSurfSec}, we will show $g$ can be extended continuously to the boundary of $\Sigma_-^o$ to get an orientation preserving homeomorphism $g: \Sigma_- \to \Sigma_-$ (see Proposition \ref{prop bound cont}). By Theorem \ref{th: z symm brake}, $\zeta_{bz}(t^*) \in \Sigma_-^o$ is a fixed point of $g$. After blowing this point up (see Subsection \ref{sub: AnnulusMap} for the details), we can get an orientation preserving homeomorphism of the annulus
$$ f: \mathbb{A} \to \mathbb{A}, \; \text{ where } \; \mathbb{A} = \mathbb{S} \times [0, 1] = \mathbb{R} /\mathbb{Z} \times [0, 1]. $$
Since $f|_{\mathbb{S} \times \{0 \}}$ and $f|_{\mathbb{S} \times \{1 \}}$ are orientation preserving homeomorphisms, they have well-defined \emph{rotation numbers} $\rho(f|_{\mathbb{S} \times \{0 \}})$ and $\rho(f|_{\mathbb{S} \times \{ 1 \}})$ (see Proposition \ref{prop_RotationNumber}).
For simplicity, we will also call them the rotation numbers of the Euler orbit $\zeta_e$ and the $z$-symmetric brake orbit $\zeta_{bz}$, and set
$$ \rho(\zeta_e) = \rho(f|_{\mathbb{S} \times \{ 1 \}}) \; \text{ and } \; \rho(\zeta_{bz}) = \rho(f|_{\mathbb{S} \times \{0 \}}). $$


\begin{thm}\label{thm infy z brake}
If $\rho(\zeta_e) \ne \rho(\zeta_{bz})$, then
\begin{enumerate}
\item[(a).] there are infinitely many different $z$-symmetry brake orbits;
\item[(b).] there are infinitely many different type-I brake orbits which are not $z$-symmetric;
\item[(c).] there are infinitely many different type-II brake orbits;
\item[(d).] there are infinitely many z-symmetry periodic orbits which are not brake.
\end{enumerate}
\end{thm}



\begin{rem}
Previous results about the existence of brake orbits in the three body problem can be found in \cite{MMV12} and \cite{Chen13}.
\end{rem}


In general it is difficult to verify the condition required in the above theorem analytically. In Section \ref{sec: RotNum and MeanIndex}, we prove the following theorem, which shows connections between the rotation numbers of $\zeta_{e}$ and $\zeta_{bz}$ and their mean indices. The mean index of a periodic orbit is defined as the mean index of the fundamental solution of the corresponding linearized system along the periodic orbit (the details can be found in the Appendix \ref{sec: AppIndex}).

 \begin{thm}\label{thm: rot number and mean index gen} For any periodic orbit  $\zeta$ of \eqref{eq: RedHamSys}, let $\hat{i}(\zeta)$ be its mean index. Then
 \begin{equation}
 \label{rot mean index bz} \rho(\zeta_{bz})=\frac{\hat{i}(\zeta_{bz})}{2}\; (\text{mod} \; \mathbb{Z}),
 \end{equation}
 \begin{equation}
 \label{rot mean index euler } \rho(\zeta_e)=(\frac{\hat i(\zeta_e)}{2}-1)^{-1}.
 \end{equation}
 \end{thm}

Particularly in  Proposition \ref{prop z symm brake is squa}, we will show $\rho( \zeta_{bz})\neq 1/2$. As a result,
\begin{cor} \label{cor case rot nequal} $\rho(\zeta_e) \ne \rho(\zeta_{bz})$, when $\hat{i}(\zeta_e)=6$.
\end{cor}

In Section \ref{sec: Inf Period Obrits}, we will show the set of parameters $(\beta, \mathfrak{e})$ with $\hat{i}(\zeta_e)=6$ is a smooth curve in $\mathcal{D}$, which is just $\Gamma_4(1)$ in Theorem \ref{curves}. Unfortunately we are not able to verify the condition required in Theorem \ref{thm infy z brake} for other values of the parameters analytically (although some numerical results will be given in Section \ref{sec: numerical}).

On the other hand, in Section \ref{sec: Inf Period Obrits} we are able to prove that there are infinitely many $z$-symmetric periodic orbits when $\rho(\zeta_e)$ is rational, and this is true for a dense subset of $\mathcal{D}$ with positive measure as stated in the following two theorems.

\begin{thm} \label{thm: inf z-symm orbits}
When $\rho(\zeta_e) \in \mathbb{Q}$, there are infinitely many $z$-symmetric periodic orbits in $\mathfrak{M}(h, \varpi)$.
\end{thm}



\begin{thm} \label{thm:DensePosMeas}
$\mathcal{D}_\mathcal{R}=\{(\beta, \mathfrak{e}) \in \mathcal{D}: \rho(\zeta_e) \in \mathbb{Q} \}$ is a dense subset of $\mathcal{D}$ with positive measure. \end{thm}

Our paper is organized as following: in Section \ref{sec:Prepare}, we recall some basic results of the annulus maps and the Euler orbits; in Section \ref{sec:GlobalSurfSec}, the existence of global surfaces of section will be established; in Section \ref{sec:z-symm}, we give a proof of Theorem \ref{thm infy z brake}, which implies the existence of infinitely many different types of periodic orbits; in Section \ref{sec: RotNum and MeanIndex}, the connection between the rotation number and the mean index will be established which proves Theorem \ref{thm: rot number and mean index gen} and Corollary \ref{cor case rot nequal}; in Section \ref{sec: Inf Period Obrits}, we give a proof of Theorem \ref{thm: inf z-symm orbits} and \ref{thm:DensePosMeas}; in Section \ref{sec: numerical}, we give some numerical results; in Section \ref{sec: AppIndex}, some basics of the index theory for symplectic paths are given as an appendix.



\section{Preparation} \label{sec:Prepare}

\subsection{The annulus map and rotation number} \label{sub: AnnulusMap}

Consider an orientation preserving homeomorphism of the annulus, $f: \mathbb{A} \to \mathbb{A}$, we recall some basic results of $f$ that will be needed.

\begin{defi} \cite[Definition 1.1]{Fr90} \label{def: definition of rotation number}
Let $\pi: \widetilde{\ba} \to \mathbb{A}$ ($\widetilde{\ba} = \mathbb{R} \times [0, 1]$) be a universal cover and let $F: \widetilde{\ba} \to \widetilde{\ba}$ be a lift of $f$. If $x \in \widetilde{\ba}$, the \textbf{rotation number} of $x$ relative to $F$ is defined as below, if the limit exists
$$ \rho(x, F) = \lim_{n \to \infty} \frac{(F^n(x) - x)_1}{n},$$
where $(\cdot)_1$ denotes projection on the first factor of $\widetilde{\ba}$.

If $\rho(x, F)$ exists and $y = \pi(x) \in \mathbb{A}$, we define the rotation number of $y$ as $\rho(y, f)$, which is the fractional part of $\rho(x, F)$.
\end{defi}

\begin{rem}
Notice that if $\rho(x, F)$ exists, its fractional part is independent of the choice of the lift $F$ and $\rho(x', F)=\rho(x, F)$, for any $x' \in \widetilde{\ba}$ with $\pi(x') = \pi(x)$.
\end{rem}

Unlike an annulus homeomorphism, for an orientation preserving circular homeomorphism, the corresponding rotation number always exists.

\begin{prop} \cite[Proposition 4.3.5]{HK03} \label{prop_RotationNumber}
Let $f: \mathbb{S} \to \mathbb{S} $ be an orientation preserving homeomorphism and $F:\mathbb{R}\to \mathbb{R}$ a lift of $f$. Then for any $x \in \mathbb{R}$, the rotation number of $x$ relative to $F$
\beg{eqnarray} \rho(x, F)= \lim_{n \to \infty}\frac{F^n(x) -x}{n}  \end{eqnarray}
always exists, and $\rho(x, F) = \rho(x', F)$, for any $x, x' \in \mathbb{R}$.

For any $y \in \mathbb{S}$ with $y = \pi(x)$, define the rotation number $\rho(y,f)$ of $y$ as the fractional part of $\rho(x, F)$. Then $\rho(y,f)$ is independent of the choice of the lift $F$, and $f$ has a periodic point if and only if $\rho(y, f)$ is rational.
\end{prop}
\begin{rem}
For the map $f$ and $F$ in the above proposition, since $\rho(y, f)$ and $\rho(x, F)$ always exist and are independent of the choice of $y \in \mathbb{S}$ and $x \in \mathbb{R}$, we will define them as the rotation numbers of $f$ and $F$, and denote then as $\rho(f)$ and $\rho(F)$ correspondingly.
\end{rem}



By Theorem \ref{thm: GlobalSurface}, the global surface of section we obtained is a $2$-dim disk rather than an annulus. Meanwhile a fixed point can be found in the interior of the disk from Theorem \ref{th: z symm brake}. Following \cite{Fr90} we will blow up the fixed point to obtain an annulus.


Let $D\subset\mathbb{R}^2$ be a smooth contractible closed region, and $g: D \to D$ be an area preserving homeomorphism (diffeomorphic in $\mathbb{D}^o$) with an interior fixed point $x_0 \in \mathbb{D}^o$. If there is an diffeomorphism $\psi: \mathbb{D} \to D$ satisfying $\psi(0) = x_0$, where $\mathbb{D}= \{ z \in \mathbb{R}^2: \; |z| \le 1 \}$ is the standard $2$-dim disk, then
$$ \tilde{g} = \psi^{-1} \circ g \circ \psi: \mathbb{D} \to \mathbb{D}$$
is an area preserving homeomorphism with the origin fixed.

Using the following symplectic diffeomorphism
\begin{equation} \label{defi: AnnulustoDisk}
\phi: \mathbb{S} \times (0, 1] \to \mathbb{D} \setminus \{0\}; \;\; (\theta, r) \to \sqrt{\frac{r}{\pi}} e^{2\pi \sqrt{-1} \theta},
\end{equation}
we can define $f_g = \phi^{-1} \circ \tilde{g} \circ \phi: \mathbb{S} \times (0, 1] \to \mathbb{S} \times (0, 1]$ as
\begin{equation}
\label{eq: f_g} f_g(\theta, r) = \left( \frac{\tilde{g}(\sqrt{\frac{r}{\pi}} e^{2\pi \sqrt{-1} \theta})}{|\tilde{g}(\sqrt{\frac{r}{\pi}} e^{2\pi \sqrt{-1} \theta})|}, |\tilde{g}(\sqrt{\frac{r}{\pi}} e^{2\pi \sqrt{-1} \theta})|^2 \right).
\end{equation}
As $r$ goes to zero, we can extend $f_g$ continuously to $\mathbb{S} \times \{0\}$ by defining
$$ f_g(\theta, 0) =  \left( \frac{d\tilde{g}|_{0}(e^{2\pi \sqrt{-1} \theta})}{|d\tilde{g}|_{0}(e^{2\pi \sqrt{-1} \theta})|}, 0 \right). $$
This allows us to define the \emph{infinitesimal rotation number of $x_0$ with respect to $g$} as
$$ \rho(x_0, g) = \rho(f_g|_{\mathbb{S} \times \{0 \}}). $$


\begin{lem} \label{lem:RotNumx0}
Let $\lambda$ be an eigenvalue of $d\tilde{g}|_{0}$, we have
$$ \rho(x_0, g) = \begin{cases} \theta/2\pi, & \text{ if } d\tilde{g}|_0\ \text{is\ symplectically\ similar\ to}\ R(\theta) \\
1/2\ (0), & \text{ if } \lambda \text{ is negative (positive)},
\end{cases}
$$
where $R(\theta)$ represents the rotation matrix with angle $\theta$ (see \eqref{definition of R(theta)}).
\end{lem}
\begin{proof}
Let $e^{\sqrt{-1}\theta_0}$ be the eigenvector of $\lambda$. Since $\mathrm{\rho}(x_0,g)=\rho(f_g|_{\mathbb{S} \times \{0 \}})$ is constant on $\mathbb{S}$, by choosing $\theta=\theta_0$, we have $f_g(\theta_0,0)=(e^{\sqrt{-1}(\arg(\lambda)+\theta_0)},0)$ and the lemma follows from Definition \ref{def: definition of rotation number}.
\end{proof}



\subsection{The Eulerian orbit} \label{sub:Euler}

For the three body problem, an Euler orbit is a self-similar solution, where the three bodies form a Eulerian configuration (a collinear configuration) all the time, with each body moving on a Keplerian orbit around the center of mass.

In general there are three different Euler configurations depending on which body is located between the other two. For the isosceles three body problem, there is only one with $m_3$ between $m_1$ and $m_2$. As a result, when the energy $h$ and angular momentum $\varpi$ are fixed, there is essentially one Euler orbit given as 
\begin{equation}
\label{eq: EulerOrbit} \zeta_e(t) = (p_{r, e}, p_{r, z}, r_e, z_e)(t) = (\dot{r}_e(t), 0, r_e(t), 0)
\end{equation}
Notice that for the Euler orbit, $z_e(t) \equiv 0$, for all $t$. This implies $\dot{p}_{z,e}(t) = \dot{z}_e(t)\equiv 0$. Meanwhile $r_e(t)$ satisfies
\begin{equation} \label{eq:kepler}
\ddot{r}_e=-\frac{1}{\beta r_e^{2}}+ \frac{\varpi^2}{r_e^3 }.
\end{equation}
To fix the initial condition, we further assume $m_1$ is at its aphelion when $t=0$, i.e.
\begin{equation}
\label{eq: aphelion} r_e(0) = \max \{r(t): \; \forall t \in [0, T_e] \},
\end{equation}
where $T_e$ is the minimal period of the Euler orbit.

Let $\theta_e(t)$ be the corresponding function of the $\theta$ variable of $m_1$ in $\mathbb{R}^3$. Recall that $\varpi > 0$, by \eqref{eq: AngMom}, $\theta_e(t)$ is a strictly increasing. If we assume
\begin{equation}
\label{eq: tht0} \theta_e(0) =0,
\end{equation}
then $\theta_e$ can be used as a new time variable. As a result
\begin{equation}
\label{eq_re} r_e(\theta_e) =  \frac{\varpi^2\beta}{1 + \mathfrak{e} \cos \theta_e},
\end{equation}
where $\mathfrak{e}$ is the eccentricity of the corresponding Keplerian orbit satisfying
\begin{equation}
\label{eq: eccentricity} \mathfrak{e} = \sqrt{1 + 2 h \varpi^2 \beta^2}.
\end{equation}
By Proposition \ref{prop_Mcomp}, $\mathfrak{M}(h, \varpi) \simeq \mathbb{S}^3$, when
\beg{eqnarray} \label{e energy sphere}
0<\mathfrak{e}<\sqrt{\frac{8(2+\alpha)}{(4+\alpha)^2}} = \sqrt{1 - \beta^2}.
\end{eqnarray}

Consider the following linearization of \eqref{eq: RedHamSys} along the Eulerian orbit $\zeta_e(t)$,
\begin{equation}
\label{linear euler d} \dot{\xi}(t)=J B_e(t)\xi(t),
\end{equation}
where
$$ B_e(t)=\text{diag} \left(1, 1, \frac{3 \varpi^2 - 2 \beta^{-1} r_e(t)}{r_e^4(t)},  (\beta^{-1}+ 7) \frac{1}{r^{3}_e(t)} \right).$$
We may decompose \eqref{linear euler d} into the following two linear sub-systems
 \beg{eqnarray} \label{linear euler r} \dot{\xi}_1 (t)=J \left(
\begin{array}{cccc}
 1 & 0\\
 0 &3\varpi^2 r_e^{-4}(t)-2 \beta^{-1} r_e^{-3}(t) \end{array} \right)\xi_1(t); \end{eqnarray}
 \beg{eqnarray} \label{linear euler z} \dot{\xi}_2(t)=J \left(
\begin{array}{cccc}
 1 & 0\\
 0 & (7+\beta^{-1})  r^{-3}_e(t) \end{array} \right)\xi_2(t).\end{eqnarray}

Change the time variable from $t$ to $\theta_e$ in \eqref{linear euler z}. We get 
\begin{equation} \label{linear euler z theta}
\xi'_2(\theta_e)=J \begin{pmatrix}
\frac{ r_e^2(\theta_e)}{\varpi} & 0\\
 0 & \frac{7+\beta^{-1}}{\varpi r_e(\theta_e)}
\end{pmatrix} \xi_2(\theta_e).
\end{equation}
To further simplify the above linear system, let
\begin{equation}
 \label{equ: translation matrix} \mathcal{R}(\theta_e)=\left(
\begin{array}{cccc}
  \frac{r_e}{\sqrt{\varpi}} &  \frac{-r'_e}{\sqrt{\varpi} r_e^2} \\
 0 & \frac{\sqrt\varpi}{ r_e} \end{array} \right)=\begin{pmatrix}
 \frac{r_e}{\sqrt{\varpi}}  &   - \frac{\sqrt{\varpi}}{r_e} \frac{\mathfrak{e} \sin \theta_e}{1 + \mathfrak{e} \cos \theta_e} \\
 0 & \frac{\sqrt\varpi}{ r_e} \end{pmatrix}= \begin{pmatrix}
 \frac{r_e}{\sqrt{\varpi}}  &   - \frac{\beta\mathfrak{e} \sin \theta_e}{\sqrt{\varpi}^3} \\
 0 & \frac{\sqrt\varpi}{ r_e} \end{pmatrix}.
\end{equation}
By a direct computation $\eta_2(\theta_e)=\mathcal{R}(\theta_e)\xi_2(\theta_e)$ satisfies
\beg{eqnarray} \label{linear euler z new co}
\eta'_2(\theta_e)=J \left(
\begin{array}{cccc}
 1 & 0\\
 0 & 1+7\beta(1+ \mathfrak{e}\cos\theta_e)^{-1} \end{array} \right)\eta_2(\theta_e) .
 \end{eqnarray}
This is equivalent to the following second order equation
\begin{equation}
\label{second Euler z theta} \Xi''(\theta_e) =-\Xi(\theta_e)-7\beta(1+ \mathfrak{e}\cos\theta_e)^{-1}\Xi(\theta_e).
\end{equation}
Let $\Xi^D(\theta_e)$ and $\Xi^N(\theta_e)$ be solutions of the above equation under the following initial conditions correspondingly
\begin{equation}
\label{eq: DiriBC}   \Xi(0) = 0, \; \Xi'(0) = 1;
\end{equation}
\begin{equation}
\label{eq: NeumBC}   \Xi(0) = 1, \; \Xi'(0) = 0.
\end{equation}


\begin{lem} \label{lem:tht_ei} If $0 < \theta^D_{e,1} < \theta^D_{e,2} < \cdots$ and $0 < \theta^N_{e,1} < \theta^N_{e,2} < \cdots$ are the zeros of $\Xi^D(\theta_e)$ and $\Xi^N(\theta_e)$ correspondingly. Then
$$ \theta_{e,1}^D < \pi \; \text{ and } \theta^D_{e,i+1} - \theta^D_{e,i} < \pi, \; \forall i \geq 1; $$
$$ \theta_{e,1}^N < \frac{\pi}{2} \; \text{ and } \theta^N_{e,i+1} - \theta^N_{e,i} < \pi, \; \forall i \geq 1. $$
\end{lem}

\begin{proof}
Since the solutions of $\Xi''(\theta_e) = - \Xi(\theta_e)$ under the initial conditions \eqref{eq: DiriBC} and \eqref{eq: NeumBC}  are correspondingly
$$ \Xi(\theta_e) = \sin \theta_e \; \text{ and } \Xi(\theta_e) = \cos \theta_e,$$
the desired results follow directly from the Sturm-Picone comparison theorem.
\end{proof}

If we go back to the original time variable $t$, Lemma \ref{lem:tht_ei} can be used to control the zeros of solutions of the following first order linear equation, which is equivalent to \eqref{linear euler z},
\begin{equation} \label{eq: LinearEulerZ}
 \ddot{Z}(t)=-(7+\beta^{-1})\frac{Z(t)}{r_e(t)^3}.
\end{equation}
\begin{lem} \label{lem:tHatTilde} Let $\hat{Z}(t)$ be the unique solution of \eqref{eq: LinearEulerZ} under the initial condition
\begin{equation}
\label{eq: DiriBdryZ}  Z(0) = 1, \;\; \dot{Z}(0) = 0,
\end{equation}
and let $\tilde{Z}(t)$ be the unique solution of \eqref{eq: LinearEulerZ} under the initial condition
\begin{equation}
\label{eq: DiriBdryZ2}  Z(T_e/2) = 1, \;\; \dot{Z}(T_e/2) = 0.
\end{equation}
If $\hat{t}_1 >0$ is the first zero of $\hat{Z}(t)$ after $t=0$ and $\tilde{t}_1 > T_e/2$ is the first zero of $\tilde{Z}(t)$ after $t = T_e/2$, then
$$ 0 < \hat{t}_1 < T_e/2 < \tilde{t}_1 < \frac{3T_e}{4}. $$
\end{lem}
\begin{proof}
By the assumption \eqref{eq: tht0}, $\theta_e(0)=0$. Then $\theta_e(T_e/2) = \pi$. Now if we view $\Xi^N(t) =\Xi^N(\theta_e(t))$ as a function of the original time $t$, because of the connection between \eqref{linear euler z theta} and \eqref{eq: LinearEulerZ},
$$ \Xi^N(\theta_e(t)) = 0 \; \text{ if and only if } \; \hat{Z}(t)=0. $$
Therefore by \eqref{lem:tht_ei}, $ 0 < \theta_e(\hat{t}_1) = \theta^N_{e, 1}< \pi$, which implies $0 < \hat{t}_1 < \pi$.

The fact that $T_e/2 < \tilde{t}_1 < \frac{3T_e}{4}$ follows from the assumption that the start point $r_e(0)$ is at aphelion. By Kepler's second law, $\theta(\frac{3T_e}{4})>\frac{3\pi}{2}$ and then the lemma follows.


\end{proof}









\section{The Global Surfaces of Section} \label{sec:GlobalSurfSec}

In this section, we show $\Sigma_{\pm}$ are two disk-like global surfaces of section of the reduced Hamiltonian system \eqref{eq: RedHamSys}.

\begin{defi} \label{defi: GlobalSurface} \cite[Chaper 9]{FvK18} Let $X$ be a non-vanishing vector field on  $S^3$ and $\varphi^X_t$ is the corresponding flow.
A disk-like global surface of section is an embedded closed $2$-dim disk $D \subset S^3$ satisfying
\begin{enumerate}
\item[(i.)] $X$ is tangent to $\partial D$, the boundary of $D$;
\item[(ii.)] $X$ is transverse to $D^{\circ}$ the interior of $D$;
\item[(iii.)] for any $x \in S^3 \setminus \partial D$, there exist $t^+>0$ and $t^- <0$, such that both $\varphi_{t^+}^X(x)$ and $\varphi_{t^-}^X(x)$ are contained in the interior of $D$.
\end{enumerate}

\end{defi}

It is well known the topology of the energy surface $\mathfrak{M}(h, \varpi)$ depending on the values of $\varpi$ and $h$, and is not always homeomorphic to $S^3$. We reproduce the proof of Proposition \ref{prop_Mcomp} below for the seek of completeness.

\begin{proof}[Proof of Proposition \ref{prop_Mcomp}]
The energy identity $H_{\varpi}= h$ can be written as
$$ p_r^2 + p_z^2 = 2 ( h - \frac{\varpi^2}{2r^2} + U(r, z)). $$
This shows the projection of $\mathfrak{M}(h, \varpi)$ onto the $(r, z)$-plane is
$$ \{(r, z): U(r, z) \ge \frac{\varpi^2}{2r^2}- h \}.$$
Notice that for a fixed $r>0$, $U(r, z)$ is a even function of $z$, which is strictly increasing when $z \le 0$. Hence $U(r, 0)$ is the global maximum of $U$ and
$$ \mathfrak{M}(h, \varpi) = \emptyset\ \Leftrightarrow\ U(r, 0)=\frac{1}{r\beta} < \frac{\varpi^2}{2r^2} - h,\ \forall r. $$

Rewrite the above inequality according to \eqref{eq_Urz}, we have
\begin{equation}
 Q(1/r)=\frac{\varpi^2}{2} \frac{1}{r^2} - \frac{1}{r\beta} - h >0, \label{eq_IneqRZ}
\end{equation}
where $Q$ is a quadratic function of $1/r$. As a result,
$$ \mathfrak{M}(h, \varpi) \ne \emptyset \; \text{ if and only if } \; 2\varpi^2 |h| \le \beta^{-2}. $$
Moreover the projection of $\mathfrak{M}(h, \varpi)$ to the $r$-line is the same as the solutions of $Q(1/r) \le 0$, which is a interval $[r_{min}, r_{max}]$, when $0 < 2\varpi^2 |h| < \beta^{-2}$ and a point, when $2\varpi^2 |h| =\beta^{-2}$. In the later case, the entire integral manifold $\mathfrak{M}(h, \varpi)$ is just a point corresponding to the circular Eulerian orbit.

Assume $2\varpi^2|h|\leq 1$, it means that
$$\frac{\varpi^2}{2}\frac{1}{r_0^2}-\frac{1}{\beta r_0}-h<\frac{\varpi^2}{2}\frac{1}{r_0^2}-\frac{1}{r_0}-h\leq0$$
for some $r_0\in(r_{min},r_{max})$. Then
$$\frac{\varpi^2}{2}\frac{1}{r_0^2}-h\leq\frac{1}{r_0}=U(r_0,\infty)<U(r_0,z),\ \forall z.$$
Therefore, $\{r_0\}\times \mathbb{R}\in\Pi(\mathfrak{M}(h,\varpi))$, which means $\mathfrak{M}(h,\varpi)$ is unbounded in $z$ direction.

Now we assume $1<2\varpi^2|h|<\beta^{-2}$, then we have $$U(r_0,\infty)<\frac{\varpi^2}{2}\frac{1}{r_0^2}-h<U(r_0,0),$$
for all $r_0\in[r_{min},r_{max}]$. There exists $z(r_0)>0$ such that when $|z|>z(r_0)$, $\{r_0,z\}\notin \Pi(\mathfrak{M}(h,\varpi))$. Therefore, $\mathfrak{M}(h,\varpi)$ is compact. Furthermore, since the Hamiltonian $H_{\varpi}$ has only one non-degenerate critical point $(0,0,\varpi^2\beta,0)$, by Morse theory, $\mathfrak{M}(h,\varpi)$ is homotopic to $S^3$.
\end{proof}


Recall that $\varphi_t$ is the flow of the reduced Hamiltonian system \eqref{eq: RedHamSys}.
\begin{proof}[Proof of Theorem \ref{thm: GlobalSurface}]
When $z=0$, the energy identity $H_{\varpi} = h$ becomes
$$ p_r^2 + p_z^2 = 2 (h - \frac{\varpi^2}{2r^2} + U(r, 0)). $$
Since $r \in [r_{min}, r_{max}]$ (from the proof of Proposition \ref{prop_Mcomp}), the above identity implies $\mathfrak{M}(h, \varpi) \cap \{z = 0\} \simeq \mathbb{S}^2$ (the 2-dim sphere) with $\Sigma_-$ and $\Sigma_+$ homeomorphic to the lower and upper half sphere. As $ \partial \Sigma_- = \partial \Sigma_+ =\{z=0, p_z =0 \} \cap \mathfrak{M}(h, \varpi)$ is nothing but the Eulerian orbit, it must be tangent to the vector field $X_{H_\varpi}$.

Choose an arbitrary $\zeta=(p_r, p_z, r,z) \in \Sigma^o_-$. The last component of $X_{H_{\varpi}}(\zeta) = J \nabla H_{\varpi}(\zeta)$ is just $p_z$, which is negative by the definition of $\Sigma^o_-$. Hence $X_{H_{\varpi}}(\zeta)$ is not tangent to $\Sigma^o_-$, as $z \equiv 0$ in $\Sigma^o_-$.

Let $\zeta$ be a solution start from $\Sigma^o_-$. We claim there must be a $t_1>0$, such that $p_z(t_1)=0$ and $z(t_1)<0$. If not, $p_z(t) = \dot{z}(t) < 0$, for all $t >0$. This implies $z(t)$ is monotonically decreasing and $z(t) < 0$, for all $t >0$. Since $\mathfrak{M}(h, \varpi)$ is compact, $\lim_{t \to \infty} z(t) = z_0$ for some finite $z_0$. As the $\omega$-limit set of $\zeta(t)$ is non-empty, it is entirely contained in $\mathfrak{M}(h, \varpi) \cap \{z = z_0 \}$. Meanwhile the $\omega$-limit set of an orbit is invariant under $\varphi_t$. This implies there exists a solution $\tilde{\zeta}(t) = (\tilde{p}_r, \tilde{p}_z, \tilde{r}, \tilde{z})(t)$ with $\tilde{z}(t) =z_0 <0$ and $\tilde{p}_z(t) = \dot{\tilde{z}}(t)=0$, for all $t$. However this can not be true, as
$$ \dot{\tilde{p}}_z(t) = \ddot{\tilde{z}}(t) = - \frac{m^{5/2}}{\sqrt{2}} \frac{4\alpha^{-1}(1+2\alpha)\tilde{z}(t)}{(\tilde{r}^2(t)+(1+2\alpha)\tilde{z}^2(t))^{3/2}} >0, \;\; \text{if } \tilde{z}(t) <0. $$

Once $t_1$ exists, since $\ddot z(t)<-cz(t)$ for some $c>0$, then by Strum-Picone comparison theorem, there exist $t_+>t_1$ such that $p_z(t_+)<0,\ z(t_+)=0$. Hence $\Sigma_-$ is a disk-like global surface of section of $\varphi$. The proof for $\Sigma_+$ is exactly the same.
\end{proof}




\begin{figure}[H]
\centering
\includegraphics[width=0.7\textwidth]{fig_of_t+.pdf}
\caption{The graph of $t^+$ over the $(p_r, r)$-plane with $\beta = \mathfrak{e}=0.6$.}
\label{fig:t+}
\end{figure}

By Theorem \ref{thm: GlobalSurface}, for each $\zeta \in \Sigma_-^o$ there is a smallest positive time $t^+(\zeta) >0$, such that $\varphi_{t^+(\zeta)}(\zeta) \in \Sigma_-^o$ (view $t^+$ as a function on $\Sigma_-^o$, its graph over the $(p_r, r)$-plane is indicated in Figure \ref{fig:t+}). Choose $\Sigma_-^o$ as the Poincar\'e section. It allows us to define a Poincar\'e return map $g : \Sigma_-^o \to \Sigma_-^o$. First we will extend $g$ continuously to the boundary of  $\Sigma_-^o$ (in general this may not always be possible, see \cite[Chapter 9]{FvK18}).


Recall that $ \partial \Sigma_- = \{ \zeta_e(t): \; t \in [0, T_e)\}$, where $\zeta_e(t)$ is the Euler orbit and $T_e$ its minimal period. For each $t_0 \in [0, T_e)$, we define $g(\zeta_e(t_0))$ as
\begin{equation}
\label{eq: ExtBdry}  g(\zeta_e(t_0)) = \varphi_{t_2(t_0)-t_0}(\zeta_e(t_0)) = \zeta_e(t_2(t_0)),
\end{equation}
where $t_2(t_0) >t_1(t_0)>t_0$ are the first two zeros of $Z(t)$ after $t_0$ with $Z(t)$ being the unique solution of the linear equation \eqref{eq: LinearEulerZ} under the initial condition
\begin{equation}
\label{eq:NeumanBdryZt0} Z(t_0)=0, \; \dot{Z}(t_0)=-1.
\end{equation}


\begin{prop} \label{prop bound cont}
The above definition extends $g$ continuously to $\Sigma_-$ and $g|_{\partial \Sigma_-}: \partial \Sigma_- \to \partial \Sigma_-$ is an orientation preserving homeomorphism.
\end{prop}

\begin{proof}

Choose an arbitrary sequence of points $\zeta_n \in \Sigma_-^o$ with $\zeta_n \to \zeta_e(t_0)$, as $n \to \infty$. For each $\zeta_n$, let $\zeta_n(t) = (p_{r,n}, p_{z, n}, r_n, z_n)(t)= \varphi_{t-t_0}(\zeta_n)$. From the proof of Theorem \ref{thm: GlobalSurface}, there exist $t_{n,2}>t_{n,1}> t_0$, which are the first two zeros of $z_n(t)$ after $t_0$. In particular $t_{n, 2} - t_0 = t^+(\zeta_n)$ and $\varphi_{t_{n, 2} -t_0}(\zeta_n) \in \Sigma_-$. As a result, it is enough to show $t_{n,2} \to t_2(t_0)$, as $n \to \infty$.

Since $\zeta_n(t)$ is a solution of \eqref{eq: RedHamSys}, its last component $z_n(t)$ must satisfy
\begin{equation*}
\label{second z} \ddot{z}_n = \dot{p}_{z,n}=-(7+\beta^{-1}) \frac{z_n}{(r^2_n+(1+2\alpha)z^2_n)^{3/2}}.
\end{equation*}
Then $z_n(t)/|\dot{z}_n(t_0)|$ is a solution of the second order linear equation
\begin{equation} \label{eq: SecLinZ}
\ddot{Z}(t)=-(7+\beta^{-1}) f_n(t) Z(t);\;  f_n(t)=\left( r_n^2(t)+(1+2\alpha)z_n^2(t) \right)^{-\frac{3}{2}}
\end{equation}
under the initial condition \eqref{eq:NeumanBdryZt0}.


Since $(r_n, z_n)(t)$ converges to $(r_e, z_e)(t)$ uniformly on any compact time interval, for any $\varepsilon>0$ small enough, there exists a $n_{\varepsilon} >0$, such that for all $n>n_{\varepsilon}$,
\begin{equation}
\label{control r}   \frac{1-\varepsilon}{r_e(t)^{3}}< f_n(t)<  \frac{1+\varepsilon}{r_e(t)^{3}}, \; \forall t\in[t_0, T].
\end{equation}
where $T$ is a finite time large than every $t_{n, 2}$.

Let $Z_{\varepsilon^{\pm}}(t)$ be the solution of the following equation under the initial condition \eqref{eq:NeumanBdryZt0}
$$  \ddot{Z}(t)=-(7+\beta^{-1})\frac{1\pm\varepsilon}{r_e(t)^{3}}Z(t).$$
and let $t_{\varepsilon^{\pm}, 2} >t_{\varepsilon^{\pm}, 1}>t_0$ be the first two zeros of $Z_{\varepsilon^{\pm}}(t)$ after $t_0$. Since $t_{n, 2} > t_{n, 1} > t_0$ are the first two zeros of $z_n(t)/|\dot{z}_n(t_0)|$ after $t_0$, by \eqref{control r} and the Sturm-Picone comparison theorem, we have
$$ t_{\varepsilon, i}\leq t_{n,i}\leq t_{-\varepsilon, i},\; \text{ when } i =1,2. $$
Meanwhile $t_{\pm\varepsilon, i}\to t_i$, as $n \to \infty$, for $i=1,2$, because $f_n(t)$ converges uniformly to $ r_e(t)^{-3}$ on any compact time interval. This gives us the desired result.
\end{proof}

For simplicity, we omit the $h$ in $\mathcal{B}(h), \mathcal{B}_{\pm}(h)$ and the interior $\mathcal{B}^o_{\pm}(h)$ in the rest of paper.
\begin{prop} \label{exis z-symm brake}  There exist at least one $z$-symmetric brake orbit which is different from the Euler orbit.
\end{prop}
\begin{proof}
Let $\zeta(t)= (p_r, p_z, r, z)(t) = \varphi_t(\zeta)$ be an orbit with the initial condition $\zeta \in \mathcal{B}_+^o$. Then $z(t)$ satisfies
$$ \ddot{z}=-(7+\beta^{-1}) \frac{z}{(r^2+(1+2\alpha)z^2)^{3/2}}.  $$
Since $z(0) >0$ and $\dot{z}(0)=0$, there must be a finite $t_1(\zeta) >0$ with $ \zeta(t_1(\zeta)) \in \Sigma_-^o$.  Moreover if $p_r(t_1(\zeta)) =0$, then $\varphi_t(\zeta)$ must be a $z$-symmetric brake orbit different from the Euler orbit.

To obtain the desired $z$-symmetric brake orbit, by the continuity of solutions with respect to the initial condition, it is enough to show there are two different $\zeta$'s from $\mathcal{B}^o_+$, such that the corresponding $p_r(t_1(\zeta))$'s having opposite signs (see Figure \ref{fig:angle} for a numerical picture, where the red curves are different $z$-symmetric brake orbits started from the upper half Hill's boundary).


\begin{figure}[H]
\centering
\includegraphics[width=0.6\textwidth]{fig_of_the_angle2.pdf}
\caption{$z$-symmetric brake orbits with $\beta = \mathfrak{e}=0.6$.}
\label{fig:angle}
\end{figure}


Choose a sequence of points $\zeta_n=(0, 0, r_n, z_n) \in \mathcal{B}^o_+$ satisfying $(r_n, z_n) \to (r_{max}, 0)$ as $n \to \infty$. Let $\zeta_n(t) = (p_{r, n}, p_{z, n}, r_n, z_n)(t) = \varphi_t(\zeta_n)$, for all $t \in \mathbb{R}$. Denote the first zero of $z_n(t)$ after $t=0$ as $t_{n,1}$.

Notice that $z_n(t)/|z_n(0)|$ is a solution of the linear equation \eqref{eq: SecLinZ} under the following initial condition
\begin{equation}
\label{eq: DiriBdryZ11}  Z(0) = 1, \;\; \dot{Z}(0) = 0.
\end{equation}
If $\hat{Z}(t)$ is the unique solution of the linear equation \eqref{eq: LinearEulerZ} under the initial condition \eqref{eq: DiriBdryZ11} with $\hat{t}_1$ being the first zero of $\hat{Z}(t)$ after $t=0$. Then Lemma \ref{lem:tHatTilde} implies $0< \hat{t}_1 < T_e/2$, and hence $p_{r, e}(\hat{t}_1) < 0$.

Meanwhile by a similar argument as in the proof of Proposition \ref{prop bound cont}, one can show $t_{n, 1} \to \hat{t}_1$, as $n \to \infty$. This implies $p_r(t_{n,1}) < 0$, when $n$ is large enough.

On the other hand if we assume the sequence of points $\zeta_n= (0, 0, r_n, z_n) \in  \mathcal{B}_+^o$ satisfying $(r_n, z_n) \to (r_{min}, 0)$ as $n \to \infty$. By a similar argument as above, one can show $p_r(t_{n, 1})>0$.

\end{proof}

For the rest of the paper, denote the $z$-symmetric brake orbit obtained in the above proposition as $\zeta_{bz}(t)= (p_r^*, p_z^*, r^*, z^*)(t)$ satisfying $\zeta_{bz}(0) \in \mathcal{B}_+(h)$ and
$$ p^*_r(t^*)=0, \; z^*(t^*) =0, \;\text{ and } \; z^*(t)>0, \; \forall t \in [0, t^*), \; \text{ for some } t^* >0.$$


\begin{prop} \label{hopf link} $\zeta_{bz}(t)$ forms a Hopf link with the Euler orbit $\zeta_e(t)$ in $\mathfrak{M}(h, \varpi)$.
\end{prop}
\begin{proof} Notice that $\zeta_{bz}(t)$ intersects $\Sigma^o$ only once during a minimal period $T_{bz}=4t^*$. As a result, the linking number of $\zeta_{bz}$ and $\zeta_{e}$ is $\pm1$. Hence it is enough to show  $\zeta_{bz}$ is not self-knotted.

For this, we define an isotopic $F_s(t), (s,t)\in[0,1]\times[0,T_{bz}]$ as
$$ F_s(t) = \begin{cases} (\sqrt{1-s}p_r^*,-(\sqrt{s}(p_r^*)^2+ (p_z^*)^2)^{1/2}, r^*, z^*)(t), & \text{ if } t\in[0,T_{bz}/2],  \\
(\sqrt{1-s}p_r^*, (\sqrt{s}(p_r^*)^2+ (p_z^*)^2)^{1/2}, r^*, z^*)(t), & \text{ if } t\in[T_{bz}/2, T_{bz}].
\end{cases}
$$
Since $-p_z^*(T_{bz}-t)=p_z^*(t)<0$ for $t\in(0,T_{bz}/2)$, $F_s$ is an isotopy with $F_0(t)=\zeta_{bz}(t)$. Then $\zeta_{bz}$ is unknotted if and only if $F_1$ unknotted.

Recall that  $\mathcal{H}(h) $ is the Hill's region and let
 $$    \hat{\Sigma}_\pm =\{ (0, \pm(h-V_{\varpi}(r, z))^{1/2}, r, z):  (r, z)\in \mathcal{H}(h)  \}.$$
Then $\hat{\Sigma}_\pm \simeq \mathbb{D}$ with $\hat{\Sigma}_+\cap \hat{\Sigma}_-=\mathcal{B}(h)$, and $\hat{\Sigma}_+\cup \hat{\Sigma}_-$ is a $2$-dim sphere embedded in   $\mathfrak{M}(h, \varpi)$. Since $F_1$ is an embedded loop in $\hat{\Sigma}_+\cup \hat{\Sigma}_-$, then it is unknotted by the Jordan curve theorem.


\end{proof}

Theorem \ref{th: z symm brake} now follows directly form Proposition \ref{exis z-symm brake} and \ref{hopf link}.
 
Since $\Sigma^o_\pm$ is a global surface of section, $X_H$ must be transversal to $\Sigma^o_\pm$. Then $(\Sigma^o_\pm,\Omega|_{\Sigma^o_\pm})$ is a symplectic surface with $\Omega=dp_r\wedge dr+dp_z\wedge dz$. A direct computation shows the return map $g:(\Sigma^o_\pm,\Omega|_{\Sigma^o_\pm})\rightarrow (\Sigma^o_\pm,\Omega|_{\Sigma^o_\pm})$ is symplectic. 

By arguments similarly to the above proofs, we can obtain two more homeomorphisms $g_1: \Sigma_-\rightarrow \Sigma_+$ and $\ g_2: \Sigma_+\rightarrow \Sigma_-$ with
$$ g = g_2 \circ g_1. $$

Let $\Upsilon$ be a $2$-dim connected region defined as below with area form $dp_r\wedge dr$
\begin{equation}
\label{eq: Upsiion} \Upsilon=\{(p_r,r)\in\mathbb{R}^2|\ p_r^2\leq2h-\frac{\varpi^2}{r^2}+\frac{2}{r\beta}\}.
\end{equation}
The natural projections $\mathcal{P_{\pm}}: \Sigma_\pm\to \Upsilon $ given below  are symplectic homeomorphisms
$$\mathcal{P_{\pm}}: \Sigma_{\pm} \rightarrow \Upsilon; \; (p_{r},p_{z},r,0) \mapsto  (p_{r},r). $$

Then the following maps are area preserving homeomorphisms from $\Upsilon$ to itself:
\begin{equation}
\label{eq:gbar}  \bar{g} = \mathcal{P}_- \circ g \circ \mathcal{P}_-^{-1}, \;\; \bar{g}_1= \mathcal{P}_+ \circ g_1 \circ \mathcal{P}_-^{-1}, \;\; \bar{g}_2= \mathcal{P}_- \circ g_2 \circ \mathcal{P}_+^{-1}.
\end{equation}
Obviously $\bar{g} = \bar{g}_2 \circ \bar{g}_1$. Meanwhile the $z$-symmetric property of the system \eqref{eq: RedHamSys}, see \eqref{eq: z-sym}, implies $\bar{g}_1 = \bar{g}_2$. Therefore
\begin{equation}
\label{eq: gbar} \bar{g}=\bar g_i^2, \;\; \forall i =1, 2.
\end{equation}


\begin{prop}\label{prop z symm brake is squa} $ \rho(\zeta_{bz}) \neq 1/2. $
\end{prop}
\begin{proof}
Recall that
$$ \rho(\zeta_{bz}) = \rho(f|_{\mathbb{S} \times \{0\}}) = \rho(\zeta_{bz}(t^*), g),$$
where $\zeta_{bz}(t^*)$ is a fixed point of $g$ as given in Theorem \ref{th: z symm brake}, and $\rho(\zeta_{bz}(t^*), g)$ is the infinitesimal rotation number of $\zeta_{bz}(t^*)$ with respect to $g$.

Since $\zeta_{bz}$ a $z$-symmetric brake orbit, $\mathcal{P}_-(\zeta_{bz}(t^*))$ is a fixed point of $\bar{g}_1$. Then \eqref{eq: gbar} implies
\[ d\bar{g}|_{\mathcal{P}_-(\zeta_{bz}(t^*))}=\left(d \bar{g}_1|_{\mathcal{P}_-(\zeta_{bz}(t^*))}\right)^2 , \]
which implies that no eigenvalue of $d\bar{g}|_{\mathcal{P}_-(\zeta_{bz}(t^*))}$ on the negative real line. Then so is $dg|_{\zeta_{bz}(t^*)}$. The desired result now follows from Lemma \ref{lem:RotNumx0}.
\end{proof}



\section{Existence of different types of periodic orbits} \label{sec:z-symm}

The purpose of this section is to give a proof of Theorem \ref{thm infy z brake}. Recall that $\zeta_{bz}(t)$ is the $z$-symmetric brake orbit obtained in Proposition \ref{exis z-symm brake} with $\zeta_{bz}(0) \in \mathcal{B}_+$ and $\zeta_{bz}(t^*)= \Sigma_-$. As a result, $\zeta_{bz}(2t^*) \in \mathcal{B}_-$, and $\zeta_{bz}(0)$, $\zeta_{bz}(2 t^*)$ divide $\mathcal{B}_{+}$, $\mathcal{B}_-$ into the following closed curves correspondingly
$$ \mathcal{B}_+ = \mathcal{B}_{+, 1} \cup \mathcal{B}_{+, 2}, \;\; \mathcal{B}_-= \mathcal{B}_{-, 1} \cup \mathcal{B}_{-,2},$$
with $\mathcal{B}_{+, i}$ and $\mathcal{B}_{-, i}$ symmetric with respect to the $z$-axis (see Figure \ref{fig: CurvesB}).
\begin{figure}
 \centering
\includegraphics[width=0.6\textwidth]{fig_of_B1.pdf}
\caption{ $\mathcal{B}_{\pm, i}$ with $\beta = \mathfrak{e}=0.6.$ }
\label{fig: CurvesB}
\end{figure}

Denote the image of $\mathcal{B}_{\pm, i}$ under the flow $\varphi_t$ at the first positive time it reaches $\Sigma_-$ as $\mathcal{C}_{\pm, i}$, for $i=1,2$, and $\mathcal{C}_{\pm} = \mathcal{C}_{\pm, 1} \cup \mathcal{C}_{\pm, 2}$. Notice that $\mathcal{C}_{\pm, i}$ starts from $\zeta_{bz}(t^*)$ and ends in the Hill's boundary $\mathcal{B}$ and $\zeta_{bz}(t^*)$ divides $\mathcal{D} =\Sigma_- \cap \{p_r=0\}$ into two closed sub-curves, which will be denoted as $\mathcal{D}_{1}, \mathcal{D}_2$ (see Figure \ref{fig:CurvesCD}).


\begin{figure}
\centering
\includegraphics[width=0.6\textwidth]{fig_of_ACD2.pdf}
\caption{ $\mathcal{C}_{\pm, i}$ and $\mathcal{D}_i$ with $\beta = \mathfrak{e}=0.6$. }
\label{fig:CurvesCD}
\end{figure}



\begin{prop}
\label{prop:z-symm brake} If $\zeta(t) = (p_r, p_z, r, z)(t)$ is a $z$-symmetric brake orbit different from $\zeta_{bz}(t)$ or the Euler orbit $\zeta_e(t)$, then it must be a type-I brake orbit and it can not intersection both $\mathcal{B}_{+, i_1}$ and $\mathcal{B}_{-, i_2}$, when $i_1 \ne i_2$.
\end{prop}

\begin{proof}
Since $\zeta$ is a brake orbit, there must exist $t_1 < t_2$ with $\zeta(t_1) \ne \zeta(t_2)$,
$$ \zeta(t_1), \zeta(t_2) \in \mathcal{B} \; \text{ and } \zeta(t) \notin \mathcal{B}, \; \forall t \in (t_1, t_2). $$

First let's assume $\zeta(t)$ is type-II, then either both $\zeta(t_1)$ and $\zeta(t_2)$ belong to $\mathcal{B}^o_+$, or to $\mathcal{B}^o_-$. Without loss of generality, let's say the former holds.

As $\zeta(t)$ is also $z$-symmetric, there is a $t_0 \in (t_1, t_2)$ with $p_r(t_0)=z(t_0)=0$. Then \eqref{eq: z-sym} implies
$$ \zeta(2t_1 -t_0), \zeta(2t_2 -t_0) \in \mathcal{B}_-^o \; \text{ and } \; \zeta(2t_1 -t_0) \ne \zeta(2t_2 -t_0). $$
This means $\zeta(t)$ have four different intersections with $\mathcal{B}$, which is absurd.

For the second part, without loss of generality let's assume $\zeta(t_1) \in \mathcal{B}_{+, i_1}$ and $\zeta(t_2) \in \mathcal{B}_{-, i_2}$. Then \eqref{eq: z-sym} implies
$$ \zeta(2t_1 -t_0) \in \mathcal{B}_{-, i_1} \; \text{ and }  \zeta(2t_2 -t_0) \in \mathcal{B}_{+, i_2}. $$
As $i_1 \ne i_2$, and $\zeta(t)$ is neither $\zeta_{bz}(t)$ nor $\zeta_e(t)$, all these four intersections with $\mathcal{B}$ are different from the other, which is absurd.
\end{proof}
Let $g: \Sigma_- \to \Sigma_-$ be the homeomorphism given in the previous section.
\begin{lem} \label{lem:PerOrbitSigma}
Let $\zeta(t)=\varphi_t(\zeta_0)$ be an orbit with the initial condition $\zeta_0 \in \Sigma_-^o$. If there is an positive integer $n \in \mathbb{Z}^+$ satisfying
\begin{enumerate}
\item[(a).] $\zeta_0 \in g^n(\mathcal{C}_+) \cap \mathcal{D}$, then $\zeta(t)$ a $z$-symmetric brake orbit;
\item[(b).] $\zeta_0 \in g^{n}(\mathcal{C}_{+, i_0}) \cap \mathcal{C}_{-, i_1}$ with $i_0 \ne i_1$, then $\zeta(t)$ is a type-I brake orbit which is not $z$-symmetric;
\item[(c).] $\zeta_0 \in g^n(\mathcal{C}_{+, i_0}) \cap \mathcal{C}_{+, i_1}$ with $i_0 \ne i_1$, then $\zeta(t)$ is a type-II brake orbit;
\item[(d).] $\zeta_0 \in g^n(\mathcal{D}_{i_0}) \cap \mathcal{D}_{i_1}$ with $i_0 \ne i_1$, then $\zeta(t)$ is a $z$-symmetry periodic orbit which is not brake.
\end{enumerate}
\end{lem}

\begin{proof}
(a). By the definitions of $\mathcal{C}_+$ and $\mathcal{D}$, there exist $t_0 < 0 < t_1$, such that
$$ \zeta(t_0) \in \mathcal{B}_+, \; \zeta(t_1) \in \mathcal{D}, \; \text{ and } \zeta(t) \notin \mathcal{B}, \; \forall t \in (t_0, t_1). $$
the desired result then follows from \eqref{eq: brake} and \eqref{eq: z-sym}.

(b). By the definition of $\mathcal{C}_{\pm, i}$, there must exist $t_0 < t_1<0$, such that $\zeta(t_0) \in \mathcal{B}_{+, i_0}$ and $\zeta(t_1) \in \mathcal{B}_{-, i_1}. $ This shows $\zeta(t)$ must be a type-I brake orbit.

To see it is not $z$-symmetric, we just need to show $\zeta(t) \cap \mathcal{D} = \emptyset$, for all $t \in (t_0, t_1)$. By a contradiction argument, let's assume there is a $t_2 \in (t_0, t_1)$ with $\zeta(t_2) \in \mathcal{D}$. Then the $z$-symmetric property implies $\zeta(2 t_2 -t_0) \in \mathcal{B}_{-, i_0}$. Since $i_0 \ne i_1$, $\zeta(2 t_2 -t_0) \ne \zeta(t_1)$. As a result the brake orbit $\zeta(t)$ have three different intersections with $\mathcal{B}$, which is absurd.

(c). Like above we can find $t_0 < t_1 < 0$ satisfying $ \zeta(t_0) \in \mathcal{B}_{+, i_0}$ and $\zeta(t_1) \in \mathcal{B}_{+, i_1}.$ Since $i_0 \ne i_1$, we must have $\zeta(t_0) \ne \zeta(t_1)$. As a result, one can not find another moment $t_2$ with $\zeta(t_2) \in \mathcal{B}_{-}$. As otherwise, $\zeta(t)$ will have three different intersection points with $\mathcal{B}$, which is absurd. This means $\zeta(t)$ can not be a type-I orbit, so it must be a type-II orbit.

(d). As $\zeta_0 \in g^n(\mathcal{D}_{i_0}) \cap \mathcal{D}_{i_1}$, we can find a $t_0 < 0$ with
$$ \zeta(t_0) \in \mathcal{D}_{i_0} \; \text{ and } \; \zeta(t) \notin \mathcal{D}, \; \forall t \in (t_0, 0). $$
By the definition of $\mathcal{D}_{i_0}$ and $\mathcal{D}_{i_1}$, $\zeta(t)$ must be $z$-symmetric and periodic. To show it is not brake, it is enough to show $\zeta(t) \notin \mathcal{B}$, for any $t \in (t_0, 0)$.  By a contradictory argument, assume there is a $t_1 \in (t_0, 0)$ with $\zeta(t_1) \in \mathcal{B}$. Then \eqref{eq: brake} and \eqref{eq: z-sym} imply
$$ \zeta(4t_1 -3t_0) = \zeta(t_0) \in \mathcal{D}_{i_0} \; \text{ and } \zeta(t) \notin \mathcal{D}, \forall t \in (t_0, 4t_1-3t_0). $$
As a result, we must have
$$ 4t_1- 3t_0 =0 \; \text{ and } \; \zeta(4t_1- 3t_0) = \zeta(0), $$
which is impossible, as $i_0 \ne i_1$ and $\zeta(t) \ne \zeta_{bz}(t)$.
\end{proof}




Recall that $\zeta_{bz}(t^*)$ is a fixed point of $g: \Sigma_- \to \Sigma_-$. Following the blow-up and extension process given in Subsection \ref{sub: AnnulusMap}, we can get an orientation preserving diffeomorphism $ f: \mathbb{A} \to \mathbb{A}.$ By abuse of notation, the corresponding images of $\tilde{\mathcal{C}}_{\pm, i}$ and $\tilde{\mathcal{D}}_i$ in $\mathbb{A}$ will be denoted by the same symbol. The following result is a simple corollary of Lemma \ref{lem:PerOrbitSigma}.




\begin{cor}
\label{cor:fCurveIntsect} Given an $x \in \mathbb{A}^o$ and a positive integer $n \in \mathbb{Z}^+$.
\begin{enumerate}
\item[(a).] Every $x \in f^n(\tilde{\mathcal{C}}_+) \cap \tilde{\mathcal{D}}$ corresponds to the initial condition of a $z$-symmetry brake orbit of \eqref{eq: RedHamSys}.
\item[(b).] Every $x \in f^n(\tilde{\mathcal{C}}_{+, i_0}) \cap \tilde{\mathcal{C}}_{-, i_1}$ $(i_0 \ne i_1)$ corresponds to the initial condition of a type-I brake orbit of \eqref{eq: RedHamSys} which is not $z$-symmetric.
\item[(c).] Every $x \in f^n(\tilde{\mathcal{C}}_{+, i_0}) \cap \tilde{\mathcal{C}}_{+, i_1}$ $(i_0 \ne i_1)$ corresponds to the initial condition of a type-II brake orbit of \eqref{eq: RedHamSys}.
\item[(c).] Every $x \in f^n(\tilde{\mathcal{D}}_{i_0}) \cap \tilde{\mathcal{D}}_{i_1}$ $(i_0 \ne i_1)$ corresponds to the initial condition of a $z$-symmetry orbit of \eqref{eq: RedHamSys} which is not brake.
\end{enumerate}
\end{cor}
Now we are ready to prove Theorem \ref{thm infy z brake}.
\begin{proof} [proof of Theorem \ref{thm infy z brake}]
Recall that $\pi: \widetilde{\ba} \to \mathbb{A}$ is the universal cover. Denote the preimages of $\tilde{\mathcal{C}}_{\pm, i}, \tilde{\mathcal{D}}_i$, $i =1,2$, in $\tilde{\mathbb{A}}$ as
$$\pi^{-1}(\tilde{\mathcal{C}}_{\pm, i}) = \cup_{j \in \mathbb{Z}} (\tilde{\mathcal{C}}_{\pm, i}+j), \;\;\pi^{-1}(\tilde{\mathcal{D}}_i) = \cup_{j \in \mathbb{Z}} (\tilde{\mathcal{D}}_i+j).  $$


Choose an arbitrary pair of $(i_1, i_2) \in \{1, 2\}^2$ with $i_0 \ne i_1$, by Corollary \ref{cor:fCurveIntsect}, it is enough to show
$$ \#(f^{n}(\tilde{\mathcal{C}}_{+,i_0}) \cap \tilde{\mathcal{D}}_{i_1} \cap \mathbb{A}^o) \to \infty, \; \text{ as } n \to \infty;$$
$$ \#(f^{n}(\tilde{\mathcal{C}}_{+,i_0}) \cap \tilde{\mathcal{C}}_{-,i_1} \cap \mathbb{A}^o) \to \infty, \; \text{ as } n \to \infty;$$
$$ \#(f^{n}(\tilde{\mathcal{C}}_{+,i_0}) \cap \tilde{\mathcal{C}}_{+,i_1} \cap \mathbb{A}^o) \to \infty, \; \text{ as } n \to \infty;$$
$$ \#(f^{n}(\tilde{\mathcal{D}}_{+,i_0}) \cap \tilde{\mathcal{D}}_{i_1} \cap \mathbb{A}^o) \to \infty, \; \text{ as } n \to \infty. $$
We show how the first limit can be obtained in details, while the others can be proven similarly. Let $F: \widetilde{\ba} \to \widetilde{\ba}$ be a lift of $f$. The assumption $ \rho(f|_{\mathbb{S} \times \{0\}}) \ne \rho(f|_{\mathbb{S} \times \{1\}})$ implies
\begin{equation}
\label{eq: DiffRotNumber} \rho_0 = \rho(F|_{\mathbb{R} \times \{0\}}) \ne \rho_1 = \rho(F|_{\mathbb{R} \times \{1\}}).
\end{equation}

If $x_k \in \mathbb{R} \times \{k\}$, $k=0,1$, are the two end-points of $\tilde{\mathcal{C}}_{+,i_0}+0$, then
$$ \rho(x_k, F) = \lim_{n \to \infty} \frac{(F^n(x_k) -x_k)_1}{n} = \rho_k, \;\; k =1,2.$$
This implies
$$ \lim_{n \to \infty} \frac{(F^n(x_1)- F^n(x_0))_1 - (x_1 -x_0)_1}{n} = \rho_1 -\rho_0. $$
Since $\rho_1 \ne \rho_0$ and $|(x_1 -x_0)_1|$ is finite,
\begin{equation}
 \label{eq: FirstFactorUnbdd} |(F^n(x_1)- F^n(x_0))_1| \to \infty, \; \text{ as } n \to \infty.
 \end{equation}
As a result,
$$ \# (F^n(\tilde{\mathcal{C}}_{+,i_0}+0) \cap \pi^{-1}(\tilde{\mathcal{D}}_{i_1}) \cap \widetilde{\ba}^o) \to \infty, \; \text{ as } n \to \infty. $$
Notice that if two different points $x, x^* \in F^n(\tilde{\mathcal{C}}_{+,i_0}+0) \cap \pi^{-1}(\tilde{\mathcal{D}}_{i_1}) \cap \widetilde{\ba}^o$, then $\pi(x)$ and $\pi(x^*)$ can not be the same point in $\mathbb{A}$. As otherwise the curve $f^n(\tilde{\mathcal{C}}_{+,i_0})$ will have a self-intersection point, which is a contradiction to the existence and uniqueness theorem of ODE. As a result,
$$ \#(f^{n}(\tilde{\mathcal{C}}_{+,i_0}) \cap \tilde{\mathcal{D}}_{i_1} \cap \mathbb{A}^o) \to \infty, \; \text{ as } n \to \infty.$$
\end{proof}

\section{Connection between the rotation number and mean index} \label{sec: RotNum and MeanIndex}

We give a proof Theorem \ref{thm: rot number and mean index gen} in this section, which establish connections between the rotation numbers of the $z$-symmetric brake orbit $\zeta_{bz}(t)$ and the Euler orbit $\zeta_e(t)$ and their mean indices.  It will follow directly from Theorem \ref{thm: RotNumMeanIndex1} and \ref{thm: RotNumMeanIndex2}.
 

\subsection{The mean index of the z-symmetric brake orbit }
By Lemma \ref{lem:RotNumx0}, the rotation number of $f|_{\mathbb{S} \times \{0 \}}$ only depends on the normal form of  $dg(\zeta_{bz}(t^*))$, which is a symplectic matrix. Recall that $\zeta_{bz}(t^*) \in \Sigma_-^o$ is a fixed point of $g$. Denote the fundamental solution of the linearized Hamiltonian equation \eqref{eq: RedHamSys} along $\zeta_{bz}(t)$ as $\gamma_{bz}(t)$ with $\gamma_{bz}(t^*) =I$, the identity matrix. 


Let $\{e_1, e_2, e_3, e_4\}$ be a symplectic basis of $T_{\zeta_{bz}(t^*)} \mathfrak{M}$ defined as: $e_1$ is the normal vector pointing outward of the energy sphere at $\zeta_{bz}(t^*)$; $e_2=\frac{X_{H_{\varpi}(\zeta_{bz}(t^*))}}{\|X_{H_{\varpi}(\zeta_{bz}(t^*))}\|}$; $T_{\zeta_{bz}(t^*)}\Sigma^-=\{e_3, e_4\}$. As a result, $\Omega(e_1, e_2)=\Omega(e_3, e_4)=1$ and other $\Omega(e_i, e_j)=0$.


 \begin{lem} Under the above symplectic basis,
 \begin{eqnarray}\label{equ: form of M}
N : =\gamma_{bz}(t^*)=\begin{pmatrix}1&0&0&0\\a_1&1&b_1&b_2\\ a_2&0&M_1&M_2\\ a_3&0&M_3&M_4\end{pmatrix},
\end{eqnarray}
 where $M:= dg(\zeta_{bz}(t^*))=\left(
\begin{array}{cccc}
M_1 & M_2\\
M_3 & M_4\\
\end{array} \right)\in\mathrm{Sp}(2)$ and
\beg{eqnarray} b_1-a_3M_1+a_2M_3=b_2-a_3M_2+a_2M_4=0. \label{abM} \end{eqnarray}
\end{lem}
\begin{proof}   Let $\xi_i=N e_i$ for $i=1,\cdots,4$. $N e_2=e_2$ get the form of the second column.

For any $x \in \Sigma^o_-$, $t_+(x)$ denote the first time the flow $\varphi_t$ comes back to $\Sigma^o_-$, then
\beg{eqnarray} dg(\zeta_{bz}(t^*))\xi=N\xi+(dt_+(\zeta_{bz}(t^*)\xi)X_{H_\varpi},  \end{eqnarray}
this gives the form of  $\xi_3, \xi_4$ and $dg(\zeta_{bz}(t^*))=M$.
$\Omega(\xi_3, \xi_4)=1$ implies that $M\in\mathrm{Sp}(2)$. Moreover, $\Omega(\xi_1, e_2)=1$ gives the form of $\xi_1$, and \eqref{abM} is from $\Omega(\xi_1,\xi_i)=0$ for $i=3,4$.
\end{proof}

Let $\sigma(N)$, $\sigma(M)$ be the spectra of $N, M$ correspondingly. The above lemma immediately implies
\begin{cor}  $\sigma(N)=\sigma(M)\cup\{1\}. $
\end{cor}

Recall that the mean index of $\zeta_{bz}$ is just $\hat{i}(\gamma_{bz})$.

\begin{thm}\label{thm: RotNumMeanIndex1} $\rho(\zeta_{bz}) = \rho(f|_{\mathbb{S} \times \{0\}})=\frac{\hat{i}(\gamma_{bz})}{2} \; (\text{mod} \; \mathbb{Z}).$
\end{thm}


\begin{proof}
By Lemma \ref{lem: mean index form}\begin{equation}
\hat{i}(\gamma_{bz})=i(\gamma_{bz})+S_{N}^{+}(1)-C(N)+\sum_{\theta \in(0,2 \pi)} \frac{\theta}{\pi} S_{N}^{-}\left(e^{\sqrt{-1} \theta}\right),  \label{mean index expre}
\end{equation}
where $C(N)=\sum_{0<\theta<2 \pi} S_{N}^{-}\left(e^{\sqrt{-1} \theta}\right)$. For any $\omega \in \mathbb{U}$, the unit circle in $\mathbb{C}$, $S^{\pm}_N(\omega)$ is the splitting number of $N$ at $\omega$ (see \eqref{eq: SplitNumber} for the details).

Notice that $i(\gamma_{bz})+S_{N}^{+}(1)=i_{e^{\sqrt{-1}\varepsilon}}(\gamma_{bz})$ for $\varepsilon>0$ small enough. If $\sigma(N)\in\mathbb{R}^+$, $N\in\mathrm{Sp}(2n)^+_{e^{\sqrt{-1}\varepsilon}}$. By Lemma \ref{lem: index even}, $i_{e^{\sqrt{-1}\varepsilon}}(\gamma_{bz})$ is even. Hence $\hat{i}(\gamma_{bz})$ is even and
$$\hat{i}(\gamma_{bz})/2=0 \;(\text{mod} \; \mathbb{Z}). $$

In this case, $\sigma(M)\in \mathbb{C}\setminus \mathbb{R}^+$, then
$N\in\mathrm{Sp}(2n)^-_{e^{\sqrt{-1}\varepsilon}}$ which shows that
$i_{e^{\sqrt{-1}\varepsilon}}(\gamma_{bz})$ is odd. If $\sigma(M)\in \mathbb{R}^{-}$, then  $\hat{i}(\gamma_{bz})=i_{e^{\sqrt{-1}\varepsilon}}(\gamma_{bz})$ is odd and we get
$$ \hat{i}(\gamma_{bz})/2=1/2 \; (\text{mod} \; \mathbb{Z}). $$

In the case, $\sigma(M)\in\mathbb{U}\setminus \{\pm 1\}$,  then $M\approx R(\theta)$. Since $i_{e^{\sqrt{-1}\varepsilon}}(\gamma_{bz})$ is odd, $C(M)=1$ and $\displaystyle \sum_{\theta \in(0,2 \pi)} \frac{\theta}{\pi} S_{N}^{-}\left(e^{\sqrt{-1} \theta}\right)=\frac{\theta}{\pi} $.  We have
$$\hat{i}(\gamma_{bz})/2=\theta/2\pi \; (\text{mod} \; \mathbb{Z}). $$
The desired result now follows from Lemma \ref{lem:RotNumx0}.
\end{proof}



\subsection{The mean index of the Euler orbit} \label{sec: mean index euler}


Recall that \eqref{linear euler d} is the linearization of the Hamiltonian system \eqref{eq: RedHamSys} along the Euler orbit $\zeta_e(t)$ and it can be decomposed into two linear sub-systems \eqref{linear euler r} and \eqref{linear euler z}. Denote the fundamental solution of \eqref{eq: EulerOrbFundSol} as $\gamma_e(t)$:
\begin{equation}
\label{eq: EulerOrbFundSol} \dot{\gamma}_e(t)=JB_e(t)\gamma_e(t), \; \; \gamma_e(0)=I.\end{equation}
Similarly the fundamental solutions of \eqref{linear euler r} and \eqref{linear euler z} will be denoted as $\gamma_1(t)$ and $\gamma_2(t)$ correspondingly.
Then the mean index of the Euler orbit is just $\hat{i}(\gamma_e)$ with $\gamma_e=\gamma_1\diamond \gamma_2$,  then by the  property  of mean index \begin{equation}
\label{eq: SumMeanIndex} \hat{i}(\gamma_e)=\hat i(\gamma_1)+\hat i(\gamma_2).
\end{equation}



\begin{lem}\label{mean index ga1}  $ \hat{i}(\gamma_{1})=2.$
\end{lem}
\begin{proof}
Recall that $\zeta_e(t) = (\dot{r}_e(t), 0, r_e(t), 0)$ and $r_e(t)$ satisfies (\ref{eq:kepler}). By differentiating both sides of (\ref{eq:kepler}), we get
\beg{eqnarray} \label{linear euler1} \frac{d^2}{dt^2}(\dot r_e)=(-3\varpi^2 r_e^{-4}+2\beta^{-1} r_e^{-3})\dot{r}_e.\end{eqnarray}
Then $\xi_1(t)=c\cdot(\ddot{r}_e(t),\dot{r}_e(t))$ is a periodic solution of \eqref{linear euler1}.  After fixing $t=0$ at the Aphelion and choosing $c$ properly, we have $\xi_1(0)=(1,0)^T$ and solution $\xi_1$ becomes the first column of $\gamma_1$.  Let $\xi_2$ be the second solution of  \eqref{linear euler1} which satisfied $\xi_2(0)=(0,1)^T$, then $\gamma_1=(\xi_1, \xi_2)$. Since $\gamma_1$ is symplectic, we have
$$\gamma_1(T)=\left(\begin{array}{cccc}
 1 & b\\
 0 &1 \end{array} \right)  $$ with some constant $b$.
Let $V_d=\mathbb{R}\times\{0\}$, From the above discussion and the Morse index theorem, we have
\beg{eqnarray}
\mu(V_{d},\gamma_{1}(t)V_{d}: [0,kT])&=&\#\{t:\dot{r}(t)=0,t\in[0,kT]\}\nonumber\\
&=&\#\{t:\sin\theta(t)=0,t\in[0,kT]\}\nonumber \\
&=&2k+1.
\end{eqnarray}
\beg{eqnarray}
\hat{i}(\gamma_{1})=\lim_{k\rightarrow+\infty}\frac{\mu(V_{d},\gamma_{1}(t)V_{d}: [0,kT])}{k}=2.\nonumber
\end{eqnarray}
\end{proof}

For $k \in \mathbb{Z}^+$, let $\theta^D_{e,k}$ and $\theta^N_{e,k}$ be defined as in Subsection \ref{sub:Euler}.
\begin{lem} \label{mean ind euler2}
$ \hat{i}(\gamma_2)=\lim_{k\to\infty}\frac{2\pi k}{\theta^D_{e,k}}=\lim_{k\to\infty}\frac{2\pi k}{\theta^N_{e,k}}>2.$
\end{lem}
\begin{proof} Let $\Lambda_D=\mathbb{R}\oplus\{0\}$ be the Lagrangian subspace of $\mathbb{R}^2$. By \eqref{mean index D},
\beg{eqnarray} \hat{i}(\gamma_2)=\lim_{m\to \infty}\frac{1}{m}\sum_{0<t_{0}<mT} \operatorname{dim} \gamma_2\left(t_{0}\right) \Lambda_{D} \cap \Lambda_{D}. \nonumber \end{eqnarray}
Take $m$ as $[\theta_k^D/2\pi]$, then $k=\sum_{0<t_{0}<mT} \operatorname{dim} \gamma_2 \left(t_{0}\right) \Lambda_{D} \cap \Lambda_{D}$, which implies the first equality. The second equality is from \eqref{mean index N} with same discussion.

The inequality follows from Lemma \ref{lem:tht_ei}, which implies $\theta^D_k<k\pi$.
\end{proof}

\begin{lem}\label{rotation number at Euler}
$ \rho(\zeta_e)= \rho(f|_{\mathbb{S} \times \{1\}})=2/\hat{i}(\gamma_2)$.
\end{lem}
\begin{proof} Notice that $\rho(f|_{\mathbb{S} \times \{1\}}) =\rho(g|_{\partial \Sigma_-})$. If we project $\Sigma^-$ to the $(P_r, r)$ plane, then $\partial_{P_r}, \partial_r$ is a symplectic basis of it and the Euler orbits $\zeta_e$, which is its boundary, rotates along the positive direction.

Recall that $g|_{\partial \Sigma_-}$ is given by \eqref{eq: ExtBdry}. Then by the definition of rotation number,
\beg{eqnarray} \rho(g|_{\partial \Sigma_-})=\lim_{k\to+\infty}\frac{\theta^D_{2k}}{2k\pi}=\frac{2}{\hat{i}(\gamma_2)}. \end{eqnarray}
where the second equality from Lemma \ref{mean ind euler2}.
\end{proof}
Combining the above lemma with \eqref{eq: SumMeanIndex}, immediately we get
\begin{thm} \label{thm: RotNumMeanIndex2} $ \rho(\zeta_e)=(\frac{\hat i(\gamma_e)}{2}-1)^{-1}.$
\end{thm}


\section{generic existence of infinitely many periodic orbits} \label{sec: Inf Period Obrits}
In Subsection \ref{z-symm orb and rev map}, we give a proof of Theorem \ref{thm: inf z-symm orbits} based on a result by Kang \cite{kang18} on reversible maps.
In Subsection \ref{dense euler meanindex}, by studying the mean index of the Euler orbit using the Maslov-type index theory for symplectic paths from \cite{Lon02}, a proof of Theorem \ref{thm:DensePosMeas} will be given.

\subsection{z-symmetric periodic orbits and the symmetric periodic point  }\label{z-symm orb and rev map}
Recall that $\Upsilon$ defined in \eqref{eq: Upsiion} is a $2$-dim disk. Let $\mathcal{N}$ be the reflection on $ \Upsilon$ given by
 $$(p_r,r) \mapsto (-p_r, r). $$
Following \cite{kang18}, we call a homomorphism $f: \Upsilon \to \Upsilon$ \emph{a reversible map}, if
$$ f\mathcal{N}=\mathcal{N}f^{-1},$$
and call a point $z\in \Upsilon$  \emph{a symmetric periodic point}, if
$$ f^k(z)=z, \; f^l(z)=\mathcal{N}(z) \; \text{ for some }  k,l\in \mathbb{N}. $$

\begin{lem}{\cite[Corollary 1.2]{kang18}}
\label{Lem:Kang} Every area-preserving reversible map on $\Upsilon$ has either precisely one interior symmetric fixed point with no other periodic points or infinitely many interior symmetric periodic points.
\end{lem}

\begin{rem} A famous theorem of Franks \cite{Fr92} states that an area preserving homeomorphism of the open or closed annulus which has at least one periodic point must in fact have infinitely many interior periodic points. Kang \cite{kang18} extend this result to symmetric periodic points of reversible maps. Recently, Liu. et. \cite{LWY} give a refined result.
\end{rem}



Recall that $\bar{g}$  is a homeomorphism of $ \Upsilon $ (see \eqref{eq:gbar}).
\begin{prop} \label{prop reversible g bar} $\bar{g}$ is a reversible map of $\Upsilon$, i.e., $\bar{g}\mathcal{N}\bar{g}=\mathcal{N}.$
\end{prop}
\begin{proof}
Let's first consider the interior of $\Upsilon$, take $q\in\Upsilon\setminus\partial\Upsilon$, consider equation
\beg{eqnarray}
\left\{\begin{array}{ll}
\dot{\zeta}(t)= J \nabla H_{\varpi}(\zeta(t)), &   \\
\zeta(0)=\mathcal{P}^{-1}_{-}q. &
\end{array}\right.\nonumber
\end{eqnarray}
Let $t_{1}>0$ be the first time such that $\zeta(t_{1})\in\Sigma_-$, i.e $\zeta(t_{1})=g(\zeta(0))$, define
\begin{equation}
\zeta_{1}(t) = \begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(t_1 -t), \; \text{ where } A = \begin{pmatrix}
-1 & 0 \\
0 & 1
\end{pmatrix}.\nonumber
\end{equation}
It's east to check that $\zeta_{1}(t)$ satisfy
\beg{eqnarray}
\left\{\begin{array}{ll}
\dot{\zeta_{1}}(t)= J \nabla H_{\varpi}(\zeta_{1}(t)), &   \\
\zeta_{1}(0)=\mathcal{P}^{-1}_{-}\mathcal{N}\mathcal{P}_{-}g(\zeta(0)). &
\end{array}\right.\nonumber
\end{eqnarray}
and
\begin{equation}
\zeta_{1}(t_{1}) = \begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(0)=\mathcal{P}^{-1}_{-}\mathcal{N}\mathcal{P}_{-}\zeta(0)
,\nonumber
\end{equation}
$t_{1}$ is also the first time of $\zeta_{1}(t)$ such that $\zeta_{1}(t_{1})\in\Sigma_-$, we get
\beg{eqnarray}
\zeta_{1}(t_{1})=g(\zeta_{1}(0))=g\mathcal{P}^{-1}_{-}\mathcal{N}\mathcal{P}_{-}g(\zeta(0)).\nonumber
\end{eqnarray}
hence
\beg{eqnarray}
g\mathcal{P}^{-1}_{-}\mathcal{N}\mathcal{P}_{-}g(\zeta(0))=\mathcal{P}^{-1}_{-}\mathcal{N}\mathcal{P}_{-}\zeta(0),\nonumber
\end{eqnarray}
this implies
$$
\bar{g}\mathcal{N}\bar{g}(q)=\mathcal{N}q,\ \ \forall q\in\Upsilon\setminus\partial\Upsilon,
$$
since $\bar{g}$ is continuous at $\Upsilon$, we still have
$$
\bar{g}\mathcal{N}\bar{g}(q)=\mathcal{N}q,\ \ \forall q\in\partial\Upsilon,
$$
hence $\bar{g}$ is a reversible map on $\Upsilon$.
\end{proof}
From the proposition above, $\bar{g}$ becomes a reversible map on $\Upsilon$, it's natural to consider the
relation between the symmetric point of $\bar{g}$ and the $z$-symmetric orbits of the hamiltonian flow, this relation will be given
in the following propositions.
\begin{prop} \label{z-symm oribit imples symm point}
Let $\zeta(t)$ be a z-symmetric orbit,
\begin{equation*} \zeta(t_0) \in \mathfrak{M}(h, \varpi) \cap \{z=0, p_r=0\}, \; \text{ for some } t_0 \in \mathbb{R},
\end{equation*}
and let $t_{-1}<t_{0}$ be the first time in the opposite direction of time such that $\zeta(t_{-1})=g^{-1}(\zeta(t_{0}))$,
then $\mathcal{P}_{-}\zeta(t_{-1})\in\Upsilon$ is a symmetric point of $\bar{g}$. Moreover, if $\zeta(t)$ is
a $z$-symmetric periodic orbit, but not Euler orbit, then $\mathcal{P}_{-}\zeta(t_{-1})$ will be a symmetric periodic point of $\bar{g}$.
\end{prop}
\begin{proof}
We first consider the case, $\zeta(t)$ is not the Euler orbit, i.e $ \zeta(t_0)\not\in\partial\Upsilon$.
Let $t_{1}$ be the first time of $\zeta(t)$ after $t_{0}$ such that $\zeta(t_1)\in\Sigma_-$, we have $\zeta(t_1)=g(\zeta(t_{-1}))$ if
$\zeta(t_0)\not\in\Sigma_-$ or $\zeta(t_1)=g^{2}(\zeta(t_{-1}))$ if $\zeta(t_0)\in\Sigma_-$ , define
\begin{equation}
\zeta_{1}(t)=\zeta(t_{0}+t),\ \
\zeta_{2}(t) = \begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(t_0-t),
\nonumber
\end{equation}
then $\zeta_{1}(t)$ satisfies
\beg{eqnarray}
\left\{\begin{array}{ll}
\dot{\zeta_{1}}(t)= J \nabla H_{\varpi}(\zeta_{1}(t)), &   \\
\zeta_{1}(0)=\zeta(t_0). &
\end{array}\right.\nonumber
\end{eqnarray}
since
$$
\begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(t_0)=\zeta(t_0),
$$
it's easy to check that $\zeta_{2}(t)$ also satisfies
\beg{eqnarray}
\left\{\begin{array}{ll}
\dot{\zeta_{2}}(t)= J \nabla H_{\varpi}(\zeta_{2}(t)), &   \\
\zeta_{2}(0)=\zeta(t_0). &
\end{array}\right.\nonumber
\end{eqnarray}
hence
\begin{equation}
\zeta(t_{0}+t) = \begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(t_0-t),
\nonumber
\end{equation}
since $t_{-1}<t_{0}$ be the first time in the opposite direction of time such that $\zeta(t_{-1})\in\Sigma_-$, combine with the above equation, it's easy to know that $t_{0}+(t_{0}-t_{-1})$ will be the first time of $\zeta(t)$ after $t_{0}$ such that $\zeta(t_{0}+(t_{0}-t_{-1}))\in\Sigma_-$, hence $t_{1}=t_{0}+(t_{0}-t_{-1})$, that is $t_{0}=\frac{t_{-1}+t_{1}}{2}$, we obtain
\begin{equation}
\zeta(t_{1}) = \begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(t_{-1}),
\nonumber
\end{equation}
this implies
\begin{equation}
\mathcal{P}_{-}\zeta(t_{1}) =
\mathcal{P}_{-}\begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix}\mathcal{P}^{-1}_{-}\mathcal{P}_{-}\zeta(t_{-1})=\mathcal{N}\mathcal{P}_{-}\zeta(t_{-1}),
\nonumber
\end{equation}
since
$\zeta(t_1)=g(\zeta(t_{-1}))$ or $\zeta(t_1)=g^{2}(\zeta(t_{-1}))$, we get
$$
\bar{g}(\mathcal{P}_{-}\zeta(t_{-1}))=\mathcal{N}\mathcal{P}_{-}\zeta(t_{-1}),\ \ \text{or} \ \ \bar{g}^{2}(\mathcal{P}_{-}\zeta(t_{-1}))=\mathcal{N}\mathcal{P}_{-}\zeta(t_{-1}).
$$
then $\mathcal{P}_{-}\zeta(t_{-1})\in\Upsilon$ is a symmetric point of $\bar{g}$, since $\zeta(t_{-1})\in\Sigma_-\setminus\partial\Sigma_-$,
we have $\mathcal{P}_{-}\zeta(t_{-1})\in\Upsilon\setminus\partial\Upsilon$.
Moreover, if $\zeta(t)$ is a $z$-symmetric periodic orbit, then $g^{l}(\zeta(t_{-1}))=\zeta(t_{-1})$ for some $l\in\mathbb{N}$, hence
$\bar{g}^{l}(\mathcal{P}_{-}\zeta(t_{-1}))=\mathcal{P}_{-}\zeta(t_{-1})$, we get $\mathcal{P}_{-}\zeta(t_{-1})\in\Upsilon$ is a periodic symmetric point of $\bar{g}$.

For the Euler orbit $\zeta(t)=\zeta_{e}(t)$, take a family $z$-symmetric solution $\widetilde{\zeta}_{n}(t)$ such that
\begin{equation*} \widetilde{\zeta}_{n}(t_0) \in \mathfrak{M}(h, \varpi) \cap \{z=0, p_r=0\}\cap\Sigma_-\setminus\partial\Sigma_-,\ \ \lim_{n\rightarrow\infty}\widetilde{\zeta}_{n}(t_0)=\zeta_{e}(t_{0}),
\end{equation*}
take $\widetilde{t}_{-1,n}<t_{0}$ be the first time in the opposite direction of time such that
$\widetilde{\zeta}_{n}(\widetilde{t}_{-1,n})=g^{-1}(\widetilde{\zeta}_{n}(t_{0}))$, and
$\widetilde{t}_{1,n}>t_{0}$ be the first time such that $\widetilde{\zeta}_{n}(\widetilde{t}_{1,n})=g(\widetilde{\zeta}_{n}(t_{0}))$,
then similar above analysis, we have
\begin{equation}
\widetilde{\zeta}_{n}(\widetilde{t}_{1,n}) = \begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \widetilde{\zeta}_{n}(\widetilde{t}_{-1,n}),
\nonumber
\end{equation}
this implies
\begin{equation}
\bar{g}^{2}(\mathcal{P}_{-}\widetilde{\zeta}_{n}(\widetilde{t}_{-1,n}))=\mathcal{N}\mathcal{P}_{-}\widetilde{\zeta}_{n}(\widetilde{t}_{-1,n}).\nonumber
\end{equation}
Since
$$
\lim_{n\rightarrow\infty}g^{-1}(\widetilde{\zeta}_{n}(t_0))=g^{-1}(\zeta_{e}(t_{0})),
$$
we have
$$
\lim_{n\rightarrow\infty}\widetilde{\zeta}_{n}(\widetilde{t}_{-1,n})=\zeta_{e}(t_{-1}),
$$
hence
$$
\bar{g}^{2}(\mathcal{P}_{-}\zeta_{e}(t_{-1}))=\mathcal{N}\mathcal{P}_{-}\zeta_{e}(t_{-1}),\nonumber
$$
$\mathcal{P}_{-}\zeta_{e}(t_{-1})$ is a symmetric point of $\bar{g}$, we must mention that the Euler orbit $\zeta_{e}(t)$ is always a periodic solution,
this does not imply $\mathcal{P}_{-}\zeta_{e}(t_{-1})$ is a periodic point, this is different from the interior points of $\Sigma_-$, in fact
the periodicity of $\mathcal{P}_{-}\zeta_{e}(t_{-1})$ depends on the rationality of the rotation number $\rho(\zeta_{e})$.
\end{proof}
\begin{prop} \label{symm point imples z-symm oribit}
Let $q\in\Upsilon$ be a symmetric point, i.e
$$ \bar{g}^l(q)=\mathcal{N}(q) \text{\ for\ some\ } l\in \mathbb{N},$$
and let $\zeta(t)$ solve equation
\beg{eqnarray}
\left\{\begin{array}{ll}
\dot{\zeta}(t)= J \nabla H_{\varpi}(\zeta(t)), &   \\
\zeta(0)=\mathcal{P}^{-1}_{-}q. &
\end{array}\right.\nonumber
\end{eqnarray}
then $\zeta(t)$ is a $z$-symmetric orbit. Moreover, if $q$ is a symmetric periodic point, then $\zeta(t)$ must be a $z$-symmetric periodic orbit.
\end{prop}
\begin{proof}
If $q\in\partial\Upsilon$, then $\zeta(t)$ is the Euler orbits, it is naturally a $z$-symmetric orbit, so we consider $q\in\Upsilon\setminus\partial\Upsilon$ in the following cases.

Case 1: $l=2n$ for some $n\in\mathbb{N}$, let $p=\bar{g}^{n}(q)$, then
$
\bar{g}^{2n}(q)=\mathcal{N}(q)
$
together with the fact $\bar{g}^{n}\mathcal{N}\bar{g}^{n}=\mathcal{N}$ implies that
$
\bar{g}^{2n}(q)=\bar{g}^{n}\mathcal{N}\bar{g}^{n}(q),
$
hence
$$
p=\bar{g}^{n}(q)=\mathcal{N}\bar{g}^{n}(q)=\mathcal{N}p,
$$
then the coordinate $p_{r}$ of $p$ must be $0$ and we have
$$
g^{n}(\zeta(0))=g^{n}(\mathcal{P}^{-1}_{-}q)=\mathcal{P}^{-1}_{-}\bar{g}^{n}(q)=\mathcal{P}^{-1}_{-}p\in \mathfrak{M}(h, \varpi) \cap \{z=0, p_r=0\},
$$
so $\zeta(t)$ is a $z$-symmetric orbit.

Case 2: $l=2n+1$ for some $n\in\mathbb{Z}$, then let $p=\bar{g}^{n}(q)$, similar to Case 1, we get
$$
\bar{g}(p)=\mathcal{N}p,
$$
hence
$$
g(\mathcal{P}^{-1}_{-}(p))=\mathcal{P}^{-1}_{-}\mathcal{N}\mathcal{P}_{-}(\mathcal{P}^{-1}_{-}p).
$$
Note that $\zeta(t)$ satisfies
\beg{eqnarray}
\left\{\begin{array}{ll}
\dot{\zeta}(t)= J \nabla H_{\varpi}(\zeta(t)), &   \\
\zeta(0)=\mathcal{P}^{-1}_{-}q, &
\end{array}\right.\nonumber
\end{eqnarray}
let $t_{n}>0, n\in\mathbb{Z}_{+},$ be the $n$th time such that $\zeta(t_{n})\in\Sigma_-$, then
$$
\zeta(t_{n})=g^{n}(\zeta(0))=g^{n}(\mathcal{P}^{-1}_{-}q)=\mathcal{P}^{-1}_{-}\bar{g}^{n}(q)=\mathcal{P}^{-1}_{-}p,
$$
$$
\zeta(t_{n+1})=g(\zeta(t_{n}))=g(\mathcal{P}^{-1}_{-}p)=\mathcal{P}^{-1}_{-}\bar{g}(p)=\mathcal{P}^{-1}_{-}\mathcal{N}p.
$$
Define
\begin{equation}
\zeta_{3}(t) = \begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(t_{n}+t_{n+1} -t), \nonumber
\end{equation}
then direct computation shows that $\zeta_{3}(t)$ satisfies
\beg{eqnarray}
\left\{\begin{array}{ll}
\dot{\zeta_{3}}(t)= J \nabla H_{\varpi}(\zeta_{3}(t)), &   \\
\zeta_{3}(t_n)=\mathcal{P}^{-1}_{-}p. &
\end{array}\right.\nonumber
\end{eqnarray}
Note that $\zeta(t)$ satisfies the same equation,
\beg{eqnarray}
\left\{\begin{array}{ll}
\dot{\zeta}(t)= J \nabla H_{\varpi}(\zeta(t)), &   \\
\zeta(t_{n})=\mathcal{P}^{-1}_{-}p, &
\end{array}\right.\nonumber
\end{eqnarray}
and we have
\begin{equation}
\zeta(t)=\zeta_{3}(t) = \begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(t_{n}+t_{n+1} -t). \nonumber
\end{equation}
By taking $t=\frac{t_{n}+t_{n+1}}{2}$, we have
\begin{equation}
\zeta(\frac{t_{n}+t_{n+1}}{2})=\begin{pmatrix}
A & 0 \\
0 & -A
\end{pmatrix} \zeta(\frac{t_{n}+t_{n+1}}{2}), \nonumber
\end{equation}
and this implies that the coordinate $p_{r}$ and $z$ of $\zeta(\frac{t_{n}+t_{n+1}}{2})$ must be $0$. Then
$$
\zeta(\frac{t_{n}+t_{n+1}}{2})\in\mathfrak{M}(h, \varpi) \cap \{z=0, p_r=0\},
$$
and $\zeta(t)$ is a $z$-symmetric orbit. Moreover, if $q$ is a symmetric periodic point, then $\bar{g}^{k}(q)=q$ for some $k\in\mathbb{N}$, hence
$$
g^{k}(\zeta(0))=g^{k}(\mathcal{P}^{-1}_{-}q)=\mathcal{P}^{-1}_{-}\bar{g}^{k}q=\mathcal{P}^{-1}_{-}q=\zeta(0),
$$
then $\zeta(t)$ must be a $z$-symmetric periodic orbit.
\end{proof}


\begin{proof}[Proof of Theorem \ref{thm: inf z-symm orbits}]
From Proposition \ref{prop reversible g bar}, \ref{z-symm oribit imples symm point}, we know that $\bar{g}$ is an area-preserving reversible map on $\Upsilon$, and we can get one interior symmetric fixed point of $\bar{g}$ from $\zeta_{bz}(t)$. If $\rho(\zeta_{e})\in\mathbb{Q}$, then $\bar{g}$ has an periodic point on $\partial\Upsilon$. Combining this with lemma \ref{Lem:Kang}, we get infinitely many interior symmetric periodic points in $\Upsilon$, hence from Proposition \ref{symm point imples z-symm oribit}, we obtain infinitely many $z$-symmetric periodic orbit.
\end{proof}


In Figure \ref{fig:z-sym}, we give some additional examples of $z$-symmetric orbits. All these orbits are found numerically with initial conditions from $\mathfrak{M}(h,\varpi)\cap\{p_r=0,\ z=0\}$, so they could be brake orbits as well. One can distinguish them from the number of intersection with the global surface of section $\Sigma_-$, denoted by $S$. For any integer $S \ge 1$, the numerical computation indicates the existence of at least one $z$-symmetric orbit, which intersects $\Sigma_-$ $S$ times in one period.
\begin{figure}[H]
\centering
\includegraphics[width=0.7\textwidth]{fig_of_symmetric_orbits2.pdf}
\caption{More examples on $z$-symmetric brake orbits.}
\label{fig:z-sym}
\end{figure}

\subsection{ The distribution of $\hat i(\gamma_2)$ in the parameter set}\label{dense euler meanindex}
This subsection will be devoted to the proof of Theorem \ref{thm:DensePosMeas}.

For this we need to study the distribution of $\hat i(\gamma_2)$ for $(\beta,\mathfrak{e})\in\mathcal{D}$, where $\gamma_2(t)$ is the fundamental solution of \eqref{linear euler z}. To indicate the dependence of the system on the parameters $\beta$ and $\mathfrak{e}$, we will write $\gamma_2(t)$ as $\gamma_{2, \beta, \mathfrak{e}}(t)$ from now on.

Compare to \eqref{linear euler z}, it is more convenient to study the linear system (\ref{linear euler z new co}). As we explained in Subsection \ref{sub:Euler}, $\mathcal{R}(\theta_{e})\gamma_{2, \beta, \mathfrak{e}}(t(\theta_{e}))$ is the fundamental solution of \eqref{linear euler z new co} when $\gamma_{2, \beta, \mathfrak{e}}(t)$ is the fundamental solution of \eqref{linear euler z}. By abuse of notation, we will denote $\mathcal{R}(\theta_{e})\gamma_{2, \beta, \mathfrak{e}}(t(\theta_{e}))$ by $\gamma_{2, \beta, \mathfrak{e}}(\theta_{e})$. 

Although it is enough to consider the problem for  $ (\beta, \mathfrak{e}) \in \mathcal{D}$, we shall study the distribution of $\hat{i}(\gamma_{2, \beta, \mathfrak{e}})$ for all $(\beta, \mathfrak{e})\in [0,+\infty)\times [0,1]$ for its mathematical interest.
\begin{prop}\label{index for e=0}
\noindent i) When $\beta=0$,
\beg{eqnarray}
&&i_{\omega}(\gamma_{2,0,\mathfrak{e}})=\left\{\begin{array}{ll}
1, & \text { for } \omega=1  \\
2, & \text { for } \omega\in\mathbb{U}\setminus\{1\}
\end{array}\right.\nonumber\\
&&v_{\omega}(\gamma_{2,0,\mathfrak{e}})=\left\{\begin{array}{ll}
2, & \text { for } \omega=1  \\
0, & \text { for } \omega\in\mathbb{U}\setminus\{1\}
\end{array}\right.\nonumber
\end{eqnarray}
\noindent ii) When $\mathfrak{e}=0$ and $ \omega\neq\pm1$,
\beg{eqnarray}
v_{\omega}(\gamma_{2,\beta,0})=\left\{\begin{array}{ll}
1, & \text { for } \omega=e^{\pm2\pi \sqrt{-1}\sqrt{1+7\beta}}. \\
0, & \text { otherwise.}
\end{array}\right.\nonumber
\end{eqnarray}
When $\omega=\pm 1$,
\beg{eqnarray}
&&i_{1}(\gamma_{2,\beta,0})=\left\{\begin{array}{ll}
1, & \text { if } \beta=\beta_{1}  \\
2n+1, & \text { if }\beta\in(\beta_{n},\beta_{n+1}]\text { for } n\in\mathbb{Z}^+
\end{array}\right.\nonumber\\
&&v_{1}(\gamma_{2,\beta,0})=\left\{\begin{array}{ll}
2, & \text { if } \beta=\beta_{n}\text { for } n\in\mathbb{Z}^+  \\
0, & \text { otherwise.}
\end{array}\right.\nonumber
\end{eqnarray}
\beg{eqnarray}
&&i_{-1}(\gamma_{2,\beta,0})=\left\{\begin{array}{ll}
2, & \text { if } \beta\in[0,\beta_{\frac{3}{2}}]  \\
2n+2, & \text { if } \beta\in(\beta_{n+\frac{1}{2}},\beta_{n+\frac{3}{2}}]\text { for } n\in\mathbb{Z}^+
\end{array}\right.\nonumber\\
&&v_{-1}(\gamma_{2,\beta,0})=\left\{\begin{array}{ll}
2, & \text { if } \beta=\beta_{n+\frac{1}{2}}\text { for } n\in\mathbb{Z}^+  \\
0, & \text { otherwise.}
\end{array}\right.\nonumber
\end{eqnarray}
where $$\beta_{n}=\frac{n^2-1}{7},\ \beta_{n+\frac{1}{2}}=\frac{(n+1/2)^2-1}{7},\ \ n\in\mathbb{Z}^+.$$


\end{prop}
\begin{proof}
When $\beta=0$,  the fundamental solution $\gamma_{2,\beta,\mathfrak{e}}(\theta_{e})$ has expression $e^{J_{2}\theta_{e}}=R(\theta_{e})$ ($R(\theta_e)$ is the rotation matrix, see \eqref{definition of R(theta)}).

 By \cite{Lon02} , we can
compute its Maslov-type index directly and we get conclusion $i)$.

When $\mathfrak{e}=0$, A direct computation shows
\beg{eqnarray} \label{gamm for e=0}
\gamma_{2,\beta,0}(\theta_{e})=\left(
                     \begin{array}{cc}
                       \cos\sqrt{1+7\beta}\theta_{e} & -\sqrt{1+7\beta}\sin\sqrt{1+7\beta}\theta_{e} \\
                       \frac{1}{\sqrt{1+7\beta}}\sin\sqrt{1+7\beta}\theta_{e} & \cos\sqrt{1+7\beta}\theta_{e} \\
                     \end{array}
                   \right),
\end{eqnarray}
and
$$\sigma(\gamma_{2,\beta,0}(2\pi))=\{e^{2\pi \sqrt{-1}\sqrt{1+7\beta}}, e^{-2\pi \sqrt{-1}\sqrt{1+7\beta}}\}.$$
Let $\beta_{n}=\frac{n^{2}-1}{7}$ and $\beta_{n+\frac{1}{2}}=\frac{(n+1/2)^{2}-1}{7}$, $\forall n\in\mathbb{Z}^+$. Then
\beg{eqnarray}
&&\gamma_{2,\beta,0}(2\pi)=I_{2}, \text{ if } \beta=\beta_{n}, n\in\mathbb{Z}^+,\nonumber \\
&&\gamma_{2,\beta,0}(2\pi)=-I_{2}, \text{ if } \beta=\beta_{n+\frac{1}{2}}, n\in\mathbb{Z}^+.\nonumber
\end{eqnarray}
Combining these with the expression of fundamental solution of $\gamma_{2,\beta,0}(\theta_{e})$. Following \cite{Lon02}, we can compute its Maslov-type index
directly and get conclusion $ii)$.
\end{proof}
Define
$$\beta_n(\mathfrak{e},\omega):=\inf\{\beta>0|\ i_\omega(\gamma_{2,\beta,\mathfrak{e}})\geq n\}, \forall (\mathfrak{e},\omega)\in [0,1)\times\mathbb{U}.$$
For the whole region $(\beta, \mathfrak{e})\in[0,+\infty)\times[0,1)$, we further have the following result regarding the distribution of the index and the symplectic normal form of $\gamma_{2,\beta,\mathfrak{e}}(2\pi)$.
\begin{thm}\label{curves}
For any fixed $\mathfrak{e}\in[0,1), \omega\in\mathbb{U}$, the $\omega$-index $\beta\mapsto i_{\omega}(\gamma_{2,\beta,\mathfrak{e}})$ is non-decreasing, and  $i_{\omega}(\gamma_{2,\beta,\mathfrak{e}})$ is only strictly increasing
when it passes $\beta_{n}(\mathfrak{e},\omega),n\in\mathbb{Z}^+$. Furthermore, the curves $$\Gamma_{n}(\omega)=\{(\beta_{n}(\mathfrak{e},\omega),\mathfrak{e})| \mathfrak{e}\in[0,1)\},$$ possess the following properties:

i) For any $k\in\mathbb{Z}^+, \mathfrak{e}\in(0,1), \omega\in\mathbb{U}\setminus\{1, -1\}$,
\begin{eqnarray*}
&&\beta_{2k+1}(\mathfrak{e},1)<\beta_{2k+1}(\mathfrak{e},\omega)<\beta_{2k+1}(\mathfrak{e},-1)\\
&&<\beta_{2k+2}(\mathfrak{e},-1)<\beta_{2k+2}(\mathfrak{e},\omega)<\beta_{2k+2}(\mathfrak{e},1)=\beta_{2k+3}(\mathfrak{e},1).
\end{eqnarray*}

ii) For any $\mathfrak{e}\in[0,1)$, we have
$$\beta_{1}(\mathfrak{e},1)=\beta_{2}(\mathfrak{e},1)=\beta_{3}(\mathfrak{e},1)=0,$$
especially, for $\mathfrak{e}=0$, any $k\in\mathbb{Z}^+$, we have
$$\beta_{2k}(0,1)=\beta_{2k+1}(0,1)=\beta_{k},\ \beta_{2k+1}(0,-1)=
\beta_{2k+2}(0,-1)=\beta_{k+\frac{1}{2}}.$$

iii) For any $\mathfrak{e}\in(0,1), n, k\in\mathbb{Z}^+$, $\gamma_{2,\beta_{n}(\mathfrak{e},1),\mathfrak{e}}(2\pi)= I_{2}$ and
$$ \gamma_{2,\beta_{2k+1}(\mathfrak{e},-1),\mathfrak{e}}(2\pi)\approx N(-1,-1),\ \gamma_{2,\beta_{2k+2}(\mathfrak{e},-1),\mathfrak{e}}(2\pi)\approx N(-1,1),$$
where $\approx$ means symplectic similarity.

iv) When $\gamma_{2,\beta,\mathfrak{e}}(2\pi)\approx R(\theta)$ for some $\theta\in(0,\pi)$ and thus it is strong linear stable on
the segment $\beta_{2k+1}(\mathfrak{e},1)<\beta<\beta_{2k+1}(\mathfrak{e},-1)$ and $$i_{1}(\gamma_{2,\beta,\mathfrak{e}})=2k+1, i_{-1}(\gamma_{2,\beta,\mathfrak{e}})=2k.$$

v) When $\gamma_{2,\beta,\mathfrak{e}}(2\pi)\approx D(\lambda)$ for some $0>\lambda\neq 1$ and thus it is hyperbolic on
the segment $\beta_{2k+1}(\mathfrak{e},-1)<\beta<\beta_{2k+2}(\mathfrak{e},-1)$ and
$$i_{1}(\gamma_{2,\beta,\mathfrak{e}})=2k+1, i_{-1}(\gamma_{2,\beta,\mathfrak{e}})=2k+1.$$

vi) When $\gamma_{2,\beta,\mathfrak{e}}(2\pi)\approx R(\theta)$ for some $\theta\in(\pi,2\pi)$ and thus it is strong linear stable on
the segment $\beta_{2k+2}(\mathfrak{e},-1)<\beta<\beta_{2k+2}(\mathfrak{e},1)$ and $$i_{1}(\gamma_{2,\beta,\mathfrak{e}})=2k+1, i_{-1}(\gamma_{2,\beta,\mathfrak{e}})=2k+2.$$

vii) $\Gamma_{n}(\omega)$ is real analytic curves,  i.e.  $\beta_{n}(\mathfrak{e},\omega)$ is analytic function of $\mathfrak{e}$. Moreover,
$\beta'_{n}(0,\omega)=0$.

\end{thm}

\begin{rem}\label{mean index equa 6} When $(\beta, \mathfrak{e})\in \Gamma_4(1)=\Gamma_5(1)$,  $i_1(\gamma_{2,\beta, \mathfrak{e}})=3$. By \eqref{mean index 2}, $\hat{i}(\gamma_2)=4$. This then implies $$\hat{i}(\gamma_{\zeta_e})=6, $$
which is precisely the condition in Corollary \ref{cor case rot nequal}.
\end{rem}


In order to proof Theorem $\ref{curves}$, let's consider a family of operator,
\beg{eqnarray}
A(\beta,\mathfrak{e})=-\frac{d^{2}}{d\theta_{e}^{2}}-1-\frac{7\beta}{1+\mathfrak{e}\cos\theta_{e}}, \ \ (\beta, \mathfrak{e})\in\mathbb{R}\times[0,1)
\end{eqnarray}
it is a self-adjoint operator defining on domain
\beg{eqnarray}
\bar{D}(\omega,2\pi)=\{W^{2,2}([0,2\pi], \mathbb{C})|y(2\pi)=\omega y(0), \dot{y}(2\pi)=\omega \dot{y}(0)\},
\end{eqnarray}
we set $v_{\omega}(A(\beta,\mathfrak{e})) =\dim\ker(A(\beta,\mathfrak{e}))$ and
denote by $\mathcal{M}_\omega(A(\beta,\mathfrak{e}))$ its Morse index, then from \eqref{more equa mas omega}, we have the fact
\beg{eqnarray}\label{index equ}
\mathcal{M}_{\omega}(A(\beta,\mathfrak{e}))=i_{\omega}(\gamma_{2,\beta,\mathfrak{e}}), \ \ v_{\omega}(A(\beta,\mathfrak{e}))=\nu_{\omega}(\gamma_{2,\beta,\mathfrak{e}}).
\end{eqnarray}
We call $A(\beta,\mathfrak{e})$ is $\omega$ degenerate, if $v_{\omega}(A(\beta,\mathfrak{e}))\neq 0$, and define
\beg{eqnarray}
&&B(\mathfrak{e},\omega)=A(-\frac{2}{7},\mathfrak{e})^{-\frac{1}{2}}\frac{1}{1+\mathfrak{e}\cos\theta_{e}}A(-\frac{2}{7},\mathfrak{e})^{-\frac{1}{2}},\nonumber\\
&&A(-\frac{2}{7},\mathfrak{e})=-\frac{d^{2}}{d\theta_{e}^{2}}-1+\frac{2}{1+\mathfrak{e}\cos\theta_{e}},\nonumber
\end{eqnarray}
on domain $\bar{D}(\omega,2\pi)$,
then we have following lemmas, the idea of proof is from \cite{HLS14}. In \cite{Zhou17a}, similar ideas have also been used to study the linear stability and index distributions of elliptic Euler solutions. We mention that for $\beta\in[-\frac{1}{7},0]$, this operator has been studied in
\cite{Zhou17b}.
\begin{lem}\label{oper decreasing}
For any $\mathfrak{e}\in[0,1)$, $A(\beta,\mathfrak{e})$ is non-increasing with respect to $\beta\in[0,+\infty)$, for any fixed $\omega\in\mathbb{U}$,
there exist $\varepsilon_{0}=\varepsilon_{0}(\beta,\mathfrak{e})>0$ small enough such that for all $\varepsilon\in(0,\varepsilon_{0})$, there
holds
\beg{eqnarray}
i_{\omega}(\gamma_{2,\beta+\varepsilon,\mathfrak{e}})-i_{\omega}(\gamma_{2,\beta,\mathfrak{e}})=\nu_{\omega}(\gamma_{2,\beta,\mathfrak{e}}),\nonumber
\end{eqnarray}
\end{lem}
\begin{proof}
The non-increasing property is easy to get from the definition of $A(\beta,\mathfrak{e})$, moreover, if $\gamma_{2,\beta,\mathfrak{e}}(2\pi)$ is non-degenerate,
then there exists an $\varepsilon_{0}=\varepsilon_{0}(\beta,\mathfrak{e})>0$ small enough such that for all $\varepsilon\in(0,\varepsilon_{0})$, there
holds $\gamma_{2,\beta+\varepsilon_{0},\mathfrak{e}}(2\pi)$ is non-degenerate, this implies
$i_{\omega}(\gamma_{2,\beta+\varepsilon,\mathfrak{e}})-i_{\omega}(\gamma_{2,\beta,\mathfrak{e}})=\nu_{\omega}(\gamma_{2,\beta,\mathfrak{e}})=0$. If
$\gamma_{2,\beta,\mathfrak{e}}(2\pi)$ is degenerate, then $\nu_{\omega}(\gamma_{2,\beta,\mathfrak{e}})\neq0$, the index relation (\ref{index equ}) implies that
there exist unit $x_{\beta,\mathfrak{e}}\in \bar{D}(\omega,2\pi)$, and $\lambda_{\beta,\mathfrak{e}}=0$ such that
$A(\beta,\mathfrak{e})x_{\beta,\mathfrak{e}}=\lambda_{\beta,\mathfrak{e}}x_{\beta,\mathfrak{e}}$, direct computation show that $\frac{\partial\lambda_{\beta,\mathfrak{e}}}{\partial\beta}<0$,
then there exists $\varepsilon_{0}=\varepsilon_{0}(\beta,\mathfrak{e})>0$ small enough such that for all $\varepsilon\in(0,\varepsilon_{0})$, there
holds $\mathcal{M}_{\omega}(A(\beta+\varepsilon,\mathfrak{e}))-\mathcal{M}_{\omega}(A(\beta,\mathfrak{e}))=v_{\omega}(A(\beta,\mathfrak{e}))$, use the index relation (\ref{index equ}) again, we have $i_{\omega}(\gamma_{2,\beta+\varepsilon,\mathfrak{e}})-i_{\omega}(\gamma_{2,\beta,\mathfrak{e}})=\nu_{\omega}(\gamma_{2,\beta,\mathfrak{e}})$.
\end{proof}

\begin{lem}\label{omega degenerate}
For any $\omega$ boundary condition and $\mathfrak{e}\in[0,1)$, $A(\beta,\mathfrak{e})$ is $\omega$ degenerate if and only if $\frac{-1}{2+7\beta}$ is an
eigenvalue of $B(\mathfrak{e},\omega)$.
\end{lem}
\begin{proof}
Suppose $\mathfrak{e}\in[0,1), A(\beta,\mathfrak{e})x=0$ for some $x\in\bar{D}(\omega,2\pi)$. Let $y=A(\beta,\mathfrak{e})^{\frac{1}{2}}x$, then we have
\beg{eqnarray}
&&A(-\frac{2}{7},\mathfrak{e})^{\frac{1}{2}}(\frac{1}{-2-7\beta}+B(\mathfrak{e},\omega))y\nonumber\\
&&=(\frac{1}{-2-7\beta}A(-\frac{2}{7},\mathfrak{e})+\frac{1}{1+\mathfrak{e}\cos\theta_{e}})x\nonumber\\
&&=\frac{1}{-2-7\beta}A(\beta,\mathfrak{e})x\nonumber\\
&&=0\nonumber
\end{eqnarray}
\end{proof}
\begin{lem}\label{1,-1-dege coexistence}
For any fixed $\mathfrak{e}\in(0,1)$, if $A(\beta,\mathfrak{e})$ is $1$-degenerate, then we have $v_{1}(A(\beta,\mathfrak{e}))=2$, and if
$A(\beta,\mathfrak{e})$ is $-1$-degenerate, then we have $v_{-1}(A(\beta,\mathfrak{e}))=1$.
\end{lem}
\begin{proof}
In fact, the above statement is just called the ``coexistence" problem (see chapter VII of \cite{MW66}). By taking
$\theta_{e}=2u, \widetilde{x}(u)=x(2u)$, the following equation
$$
A(\beta,\mathfrak{e})x(\theta_{e})=-\ddot{x}(\theta_{e})-x(\theta_{e})-\frac{7\beta x(\theta_{e})}{1+\mathfrak{e}\cos\theta_{e}}=0,
$$
is equivalent to
\beg{eqnarray}\label{ince equation}
(1+\mathfrak{e}\cos2u)\ddot{\widetilde{x}}(u)+(4(1+7\beta)+4\mathfrak{e}\cos2u)\widetilde{x}(u)=0,
\end{eqnarray}
it is a special case of ``Ince's equation", from Theorem 7.1 in \cite{MW66}, we have following discriminant polynomial
\beg{eqnarray}
Q(\mu)=2\mathfrak{e}\mu^{2}-2\mathfrak{e},\ \ Q^{*}(\mu)=\mathfrak{e}(2\mu-1)^{2}-4\mathfrak{e},
\end{eqnarray}
which show that if equation (\ref{ince equation}) has two linearly independent solution of period $\pi$, then $Q(\mu)$ has
zero at one of the points $\mu=0, \pm1, \pm2,\cdots$, and equation (\ref{ince equation}) has two linearly independent solution of period $2\pi$,
then $Q^{*}(\mu)$ vanishes for one of the values of $\mu=0, \pm1, \pm2,\cdots$. Direct computation shows that for $\mathfrak{e}\in(0,1)$,
$Q(\mu)$ has zero points $\pm1$ and $Q^{*}(\mu)$ has zero points $-\frac{1}{2}, \frac{3}{2}$, hence for $\mathfrak{e}\in(0,1)$, equation
(\ref{ince equation}) does not has two linearly independent solution of period $2\pi$, this implies if
$A(\beta,\mathfrak{e})$ is $-1$-degenerate, then $v_{-1}(A(\beta,\mathfrak{e}))=1$. Also, if equation
(\ref{ince equation}) has one solution of period $\pi$, from the definition of `order' in page98 of \cite{MW66}, it's easy to check that the order of this solution will bigger than $1$, then Theorem 7.3 of \cite{MW66} tells us that the existence of one solution of $\pi$ implies there exists another linearly independent solution of period $\pi$, it means that if
$A(\beta,\mathfrak{e})$ is $1$-degenerate, then $v_{1}(A(\beta,\mathfrak{e}))=2$. This completes the proof.
\end{proof}

From above lemmas, we can prove the properties i)-vi) in Theorem \ref{curves},
\begin{proof}
We know $\beta_{1}=0$, then for $\beta=\beta_{1}$
\beg{eqnarray}
A(\beta,\mathfrak{e})=-\frac{d^{2}}{d\theta_{e}^{2}}-1,\nonumber
\end{eqnarray}
it's easy to show that
\beg{eqnarray}
\text{ ker }(A(\beta,\mathfrak{e}))=\Big\{\cos\theta_{e},\sin\theta_{e}\Big\} \text{ on } \bar{D}(1,2\pi),
\end{eqnarray}
Let $\beta_{1}(\mathfrak{e},1)=\beta_{1}=0$, then $\Gamma_{1}(1)$ is real analytic degenerate curves.

Further, it's not hard to find out that for any fixed $\mathfrak{e}\in[0,1)$, $A(\beta,\mathfrak{e})$ is non-decreasing with respect to $\beta$, and
$$\lim_{\beta\rightarrow+\infty}\mathcal{M}_{\omega}(A(\beta,\mathfrak{e}))=\lim_{\beta\rightarrow+\infty}i_{\omega}(\gamma_{2,\beta,\mathfrak{e}})=+\infty.$$
hence there exits a series of values $\beta_{j}(\mathfrak{e},\omega)\in[0,+\infty), j\in\mathbb{Z}^+$ such that
$A(\beta_{j}(\mathfrak{e},\omega),\mathfrak{e})$ is  $\omega$-degenerate, from lemma \ref{omega degenerate}, this equivalent to $\frac{1}{-2-7\beta_{j}(\mathfrak{e},\omega)}$ is an
eigenvalue of $B(\mathfrak{e},\omega)$. Note that $B(\mathfrak{e},\omega)$ is a compact operator and self-adjoint when $\mathfrak{e}\in[0,1)$, by using similar idea in \cite{HLS14}, we get
$\frac{1}{-2-7\beta_{j}(\mathfrak{e},\omega)}$ is real analytic in $\mathfrak{e}$, this implies $\beta_{j}(\mathfrak{e},\omega)$ are real analytic functions of $\mathfrak{e}$.

Furthermore, when
$\mathfrak{e}=0$, from Proposition \ref{index for e=0} and the index relation (\ref{index equ}), we have $A(\beta,0)$ is $1$-degenerate at $\beta_{n}$, and
$-1$-degenerate at $\beta_{n+\frac{1}{2}}$, and
\beg{eqnarray}
&&i_{1}(\gamma_{2,\beta,0})=\left\{\begin{array}{ll}
1, & \text { if } \beta=\beta_{1}  \\
2n+1, & \text { if }\beta\in(\beta_{n},\beta_{n+1}]\text { for } n\in\mathbb{Z}^+
\end{array}\right.\nonumber\\
&&v_{1}(\gamma_{2,\beta,0})=\left\{\begin{array}{ll}
2, & \text { if } \beta=\beta_{n}\text { for } n\in\mathbb{Z}^+  \\
0, & \text { other cases}
\end{array}\right.,\nonumber
\end{eqnarray}
\beg{eqnarray}
&&i_{-1}(\gamma_{2,\beta,0})=\left\{\begin{array}{ll}
2, & \text { if } \beta\in[0,\beta_{\frac{3}{2}}]  \\
2n, & \text { if } \beta\in(\beta_{n+\frac{1}{2}},\beta_{n+\frac{3}{2}}]\text { for } n\in\mathbb{Z}^+
\end{array}\right.\nonumber\\
&&v_{-1}(\gamma_{2,\beta,0})=\left\{\begin{array}{ll}
2, & \text { if } \beta=\beta_{n+\frac{1}{2}}\text { for } n\in\mathbb{Z}^+  \\
0, & \text { other cases}
\end{array}\right.\nonumber
\end{eqnarray}
this combines with lemma \ref{1,-1-dege coexistence}, we can order the analytic function $\beta_{j}(\mathfrak{e},-1), j\in \mathbb{Z}^+$ in following way,
\beg{eqnarray}
&&0=\beta_{1}(\mathfrak{e},-1)=\beta_{2}(\mathfrak{e},-1)<\beta_{3}(\mathfrak{e},-1)\nonumber\\
&&\beta_{2k+1}(\mathfrak{e},-1)<\beta_{2k+2}(\mathfrak{e},-1),\ \ k\in\mathbb{Z}^+, \mathfrak{e}\in(0,1)\nonumber\\
&&\beta_{2k+1}(0,-1)=\beta_{2k+2}(0,-1)=\beta_{k+\frac{1}{2}},\ \ k\in\mathbb{Z}^+.\nonumber
\end{eqnarray}
Based on the face $1$-degenerate curves, $-1$-degenerate curves and $\omega$-degenerate curves $(\omega\in\mathbb{U}\setminus{1,-1})$ are disjoint with each other, this combines with lemma \ref{1,-1-dege coexistence}, we can order $\beta_{j}(\mathfrak{e},1), j\in \mathbb{Z}^+$ in following way,
\beg{eqnarray}
&\beta_{2k}(\mathfrak{e},1)=\beta_{2k+1}(\mathfrak{e},1)<\beta_{2k+1}(\mathfrak{e},-1)<\beta_{2k+2}(\mathfrak{e},-1)<\beta_{2k+2}(\mathfrak{e},1)
=\beta_{2k+3}(\mathfrak{e},1),\ \nonumber\\
&\beta_{1}(0,1)=\beta_{2}(0,1)=0,\ \ \beta_{2k}(0,1)=\beta_{2k+1}(0,1)=\beta_{k},\ \ k\in\mathbb{Z}^+\nonumber
\end{eqnarray}
Moreover, from Lemma \ref{oper decreasing}, we get the estimation of $1$-index, $-1$-index for $\mathfrak{e}\in(0,1), k\in\mathbb{Z}^+$,
\beg{eqnarray}
&&i_{1}(\gamma_{2,\beta,\mathfrak{e}})=2k+1,\ \ i_{-1}(\gamma_{2,\beta,\mathfrak{e}})=2k,
\ \beta\in(\beta_{2k+1}(\mathfrak{e},1), \beta_{2k+1}(\mathfrak{e},-1)], \nonumber\\
&&i_{1}(\gamma_{2,\beta,\mathfrak{e}})=2k+1,\ \ i_{-1}(\gamma_{2,\beta,\mathfrak{e}})=2k+1, \ \beta\in(\beta_{2k+1}(\mathfrak{e},-1), \beta_{2k+2}(\mathfrak{e},-1)], \nonumber\\
&&i_{1}(\gamma_{2,\beta,\mathfrak{e}})=2k+1,\ \ i_{-1}(\gamma_{2,\beta,\mathfrak{e}})=2k+2, \ \beta\in(\beta_{2k+2}(\mathfrak{e},-1), \beta_{2k+3}(\mathfrak{e},1)], \nonumber\\
\end{eqnarray}
Let $M=\gamma_{2,\beta,\mathfrak{e}}(2\pi)$, then from the index formula (See \cite{Lon02})
\beg{eqnarray}
i_{-1}(\gamma_{2,\beta,\mathfrak{e}})&=&i_{1}(\gamma_{2,\beta,\mathfrak{e}})+S^{+}_{M}(1)\nonumber\\
&+&\sum_{0<\theta<\pi}(S^{+}_{M}(e^{\sqrt{-1}\theta})
-S^{-}_{M}(e^{\sqrt{-1}\theta}))
-S^{-}_{M}(-1)\nonumber
\end{eqnarray}
for $\theta\in(0,\pi)$ and $e^{i\theta}\in\sigma(M)$, we get the following results,
\beg{eqnarray}
&&(S^{+}_{M}(e^{\sqrt{-1}\theta}), S^{-}_{M}(e^{\sqrt{-1}\theta}))=(0,1),\ \ \beta\in(\beta_{2k+1}(\mathfrak{e},1), \beta_{2k+1}(\mathfrak{e},-1))\nonumber\\
&&(S^{+}_{M}(e^{\sqrt{-1}\theta}), S^{-}_{M}(e^{\sqrt{-1}\theta}))=(0,0), \ \ S^{-}_{M}(-1)=1,\ \ \beta=\beta_{2k+1}(\mathfrak{e},-1),\nonumber \\
&&(S^{+}_{M}(e^{\sqrt{-1}\theta}), S^{-}_{M}(e^{\sqrt{-1}\theta}))=(0,0),\ \ \beta\in(\beta_{2k+1}(\mathfrak{e},-1), \beta_{2k+2}(\mathfrak{e},-1))\nonumber \\
&&(S^{+}_{M}(e^{\sqrt{-1}\theta}), S^{-}_{M}(e^{\sqrt{-1}\theta}))=(0,0),\ \ S^{-}_{M}(-1)=0,\ \ \beta=\beta_{2k+2}(\mathfrak{e},-1),\nonumber \\
&&(S^{+}_{M}(e^{\sqrt{-1}\theta}), S^{-}_{M}(e^{\sqrt{-1}\theta}))=(1,0),\ \ \beta\in(\beta_{2k+2}(\mathfrak{e},-1), \beta_{2k+3}(\mathfrak{e},1))\nonumber
\end{eqnarray}
from the splitting number, we have the symplectic normal form of $\gamma_{2,\hat\alpha,e}(2\pi)$ for $e\in(0,1)$,
\beg{eqnarray} \label{spitt region}
&&\gamma_{2,\beta,\mathfrak{e}}(2\pi)\approx R(\theta),\ \ \theta\in(0,\pi),\ \ \beta\in(\beta_{2k+1}(\mathfrak{e},1), \beta_{2k+1}(\mathfrak{e},-1))\nonumber\\
&&\gamma_{2,\beta,\mathfrak{e}}(2\pi)\approx N(-1,-1),\ \ \beta=\beta_{2k+1}(\mathfrak{e},-1)\nonumber\\
&&\gamma_{2,\beta,\mathfrak{e}}(2\pi)\approx D(\lambda),\ \ -1\neq\lambda<0,\ \ \beta\in(\beta_{2k+1}(\mathfrak{e},-1), \beta_{2k+2}(\mathfrak{e},-1))\nonumber\\
&&\gamma_{2,\beta,\mathfrak{e}}(2\pi)\approx N(-1,1),\ \ \beta=\beta_{2k+2}(\mathfrak{e},-1)\nonumber\\
&&\gamma_{2,\beta,\mathfrak{e}}(2\pi)\approx R(\theta),\ \ \theta\in(\pi,2\pi),\ \ \beta\in(\beta_{2k+2}(\mathfrak{e},-1), \beta_{2k+3}(\mathfrak{e},1))
\end{eqnarray}
from the above normal form, it's easy to know that there exists only one $\omega$-degenerate curve $(\omega\in\mathbb{U}\setminus\{1,-1\})$ on each interval
 $$\big(\beta_{2k+1}(\mathfrak{e},1), \beta_{2k+1}(\mathfrak{e},-1)\big),
\big(\beta_{2k+2}(\mathfrak{e},-1), \beta_{2k+3}(\mathfrak{e},1)\big),$$
let's order the analytic functions $\beta_{j}(\mathfrak{e},\omega), \omega\in\mathbb{U}\setminus\{1,-1\}, j\in \mathbb{Z}^+$ in following way,
\beg{eqnarray}
&\beta_{2k}(\mathfrak{e},1)=\beta_{2k+1}(\mathfrak{e},1)<\beta_{2k+1}(\mathfrak{e},\omega)<\beta_{2k+1}(\mathfrak{e},-1),\nonumber \\
&\beta_{2k+2}(\mathfrak{e},-1)<\beta_{2k+2}(\mathfrak{e},\omega)<\beta_{2k+2}(\mathfrak{e},1)=\beta_{2k+3}(\mathfrak{e},1).\nonumber
\end{eqnarray}
This completes the proof of properties i)-vi) in Theorem \ref{curves}.
\end{proof}

For property vii) of Theorem \ref{curves},
we first prove the case $\omega\neq -1$.  Although $A(\beta,\mathfrak{e})=-\frac{d^{2}}{d\theta_{e}^{2}}-1-\frac{7\beta}{1+\mathfrak{e}\cos\theta_{e}}$ is defined for  $(\beta, \mathfrak{e})\in[0,+\infty)\times[0,1)$, we can extend $\mathfrak{e}\in (-1,1)$. Recall that the $\omega$ boundary condition is $x(t)=\omega x(t+2\pi)$,
$\dot{x}(t)=\omega \dot{x}(t+2\pi)$. From \cite{HLS14}, let
$\psi(x)(t)=x(t+\pi)$, direct computation show that
\[\psi^*A(\beta,\mathfrak{e})\psi=A(\beta,-\mathfrak{e}). \]
This property shows that the $\omega$ degeneracy curve must be symmetric with respect to $\mathfrak{e}=0$. Since $v_{\omega}(\gamma_{2,\hat\alpha,0})=1$ for $\omega\neq\pm1$, then $\beta_k(\mathfrak{e},\omega)$ is even function of $\mathfrak{e}$, which implies $\beta'_k(0,\omega)=0$. This is also true for $\omega=1$, since $\beta_{2k-1}(\mathfrak{e},1)=\beta_{2k}(\mathfrak{e},1)$.


The proof of $-1$ degenerate curve is different.
Let consider the $-1$-degenerate curves at the bifurcation point $\beta_{n+\frac{1}{2}}$, we have
\begin{lem}
The tangent directions of curves $\Gamma_{n}(-1)$ at bifurcation
point $(\beta_{n+\frac{1}{2}},0)$ are given by
$$\frac{\partial\beta_{n}(\mathfrak{e},-1)}{\partial \mathfrak{e}}\Big|_{\mathfrak{e}=0}=0, \ \ n\in\mathbb{Z}^+.$$
\end{lem}
\begin{proof}
For the $-1$-degenerate curves $\Gamma_{n}(-1)$, we have
\beg{eqnarray}
\text{Ker}A(\beta_{n+\frac{1}{2}},0)=\Big\{\cos(n+\frac{1}{2})\theta_{e}, \sin(n+\frac{1}{2})\theta_{e}\Big\},\nonumber
\end{eqnarray}
take $x_{\mathfrak{e}}\in \bar{D}(-1,2\pi)$ such that
\beg{eqnarray}
A(\beta_{n}(\mathfrak{e},-1),\mathfrak{e})x_{\mathfrak{e}}=0, \nonumber
\end{eqnarray}
then
\beg{eqnarray}
<A(\beta_{n}(\mathfrak{e},-1),\mathfrak{e})x_{\mathfrak{e}},x_{\mathfrak{e}}>=0,
\end{eqnarray}
Differentiate $\mathfrak{e}$ on both sides, we get
\beg{eqnarray}
\beta_{n}^{\prime}(\mathfrak{e},-1)<\frac{\partial}{\partial\beta}A(\beta,\mathfrak{e})x_{\mathfrak{e}},x_{\mathfrak{e}}>\big|_{(\beta_{n}(\mathfrak{e},-1),\mathfrak{e})}
+<\frac{\partial}{\partial \mathfrak{e}}A(\beta,\mathfrak{e})x_{\mathfrak{e}}, x_{\mathfrak{e}}>\big|_{(\beta_{n}(\mathfrak{e},-1),\mathfrak{e})}\nonumber\\
+2<A(\beta_{n}(-1,\mathfrak{e}),\mathfrak{e})x_{\mathfrak{e}},\frac{\partial}{\partial \mathfrak{e}}x_{\mathfrak{e}})>=0\nonumber.
\end{eqnarray}
For $\mathfrak{e}=0$, $x_{\mathfrak{e}}\in\text{Ker}A(\beta_{n+\frac{1}{2}},0)$, hence $x_{0}=\eta_{1}\cos(n+\frac{1}{2})\theta_{e}+\eta_{2}\sin(n+\frac{1}{2})\theta_{e}$ with $\eta_{1}, \eta_{2}\in\mathbb{R}, \eta_{1}^{2}+\eta_{2}^{2}\neq0$, and we have
\beg{eqnarray}
\beta_{n}^{\prime}(\mathfrak{e},-1)<\frac{\partial}{\partial\hat\alpha}A(\beta_{n+\frac{1}{2}},0)x_{0},x_{0}>
+<\frac{\partial}{\partial \mathfrak{e}}A(\beta_{n+\frac{1}{2}},0)x_{0},x_{0}>=0\nonumber.
\end{eqnarray}
Direct computation shows that
\beg{eqnarray}
&&\frac{\partial}{\partial\beta}A(\beta_{n+\frac{1}{2}},0)=-\frac{7}{1+\mathfrak{e}\cos\theta_{e}}\big|_{(\beta_{n+\frac{1}{2}},0)}=-7,\nonumber\\
&&\frac{\partial}{\partial \mathfrak{e}}A(\beta_{n+\frac{1}{2}},0)=\frac{7\beta\cos\theta_{e}}{(1+\mathfrak{e}\cos\theta_{e})^{2}}\big|_{(\beta_{n+\frac{1}{2}},0)}
=7\beta_{n+\frac{1}{2}}\cos\theta_{e},\nonumber
\end{eqnarray}
this implies that
\beg{eqnarray}
<\frac{\partial}{\partial\beta}A(\beta_{n+\frac{1}{2}},0)x_{0},x_{0}>&=&-7\int_{0}^{2\pi}
(\eta_{1}\cos(n+\frac{1}{2})\theta_{e}+\eta_{2}\sin(n+\frac{1}{2})\theta_{e})^{2}d\theta_{e}\nonumber\\
&=&-(\eta_{1}^{2}+\eta_{2}^{2})\pi,\nonumber
\end{eqnarray}
\beg{eqnarray}
<\frac{\partial}{\partial \mathfrak{e}}A(\beta_{n+\frac{1}{2}},0)x_{0},x_{0}>&=&\int_{0}^{2\pi}\beta_{n+\frac{1}{2}}\cos\theta_{e}
(\eta_{1}\cos(n+\frac{1}{2})\theta_{e}+\eta_{2}\sin(n+\frac{1}{2})\theta_{e})^{2}d\theta_{e}\nonumber\\
&=&0\nonumber
\end{eqnarray}
hence,
$$\frac{\partial\beta_{n}(\mathfrak{e},-1)}{\partial \mathfrak{e}}\Big|_{\mathfrak{e}=0}=0, \ \ n\in\mathbb{Z}^+,$$
this completes the proof.
\end{proof}



\begin{proof}[Proof of Theorem \ref{thm:DensePosMeas}]
From Lemma \ref{rotation number at Euler}, $\rho(x_e)\in\mathbb{Q}$ equivalent to
$\hat{i}(\gamma_2)\in \mathbb{Q}$.   By \eqref{mean index 2}, we know that $\hat{i}(\gamma_2)\in \mathbb{Q}$ if and only if
$$ \sigma( \gamma_{2,\beta,\mathfrak{e}}(2\pi))\subset \mathbb{R}^+\cup \mathbb{U}_Q, $$
where $\mathbb{U}_Q:=\{\omega=e^{i\theta}, \theta/\pi\in\mathbb{Q}  \}$.
Let $D_k(-1):=\{ (\beta,\mathfrak{e})\in\mathcal{D}, \beta_{2k+1}(\mathfrak{e},-1)<\beta< \beta_{2k+2}(\mathfrak{e},-1)\}$,
we have
\beg{eqnarray}  \mathcal{D}_R=\left(\cup_{\omega\in \mathbb{U}_Q, n\in\mathbb{Z}^+}  \Gamma_n(\omega) \bigcup \cup_{k\in\mathbb{Z}^+} D_k(-1)  \right)\bigcap \mathcal{D}. \end{eqnarray}
Please note that $D_k(-1)$ with positive measure, we only need to prove $\mathcal{D}_R$ is dense.

From \eqref{spitt region}, we have
\beg{eqnarray} \gamma_{2,\beta,\mathfrak{e}}(2\pi)\approx R(\theta),\ \ \theta\in(0,\pi),\ \ \beta\in(\beta_{2k+1}(\mathfrak{e},1), \beta_{2k+1}(\mathfrak{e},-1)).  \end{eqnarray}
Fixing $\mathfrak{e}_0$, then $\theta=\theta(\beta)$ is continuous and monotone on $(\beta_{2k+1}(\mathfrak{e},1), \beta_{2k+1}(\mathfrak{e},-1))$, which implies
$ \theta^{-1}(\mathbb{Q}\cap(0,1)) $ is dense on $e=e_0$.  The proof of other cases is similar.
\end{proof}
Since the map $\beta\mapsto i_\pm(\gamma_{2,\beta,\mathfrak{e}})$ is nondecreasing with discrete jumping by Theorem \ref{curves}. We compute the index $i_\pm(\gamma_{2,\beta,\mathfrak{e}})$  numerically for all $(\beta,\mathfrak{e})\in\mathcal{D}$ and show these jumping as the following picture


\begin{figure}[h]
\centering
\includegraphics[width=0.6\textwidth]{fig_of_Gamma.jpg}
\caption{The figure of $\Gamma_k(\pm 1)=(\beta_k(\mathfrak{e},\pm 1),\mathfrak{e})$ in $(\beta,\mathfrak{e})\in [0,1]\times [0,1]$.}
\label{picture of Gamma}
\end{figure}

\noindent where every time $\gamma_{2,\beta,\mathfrak{e}}(2\pi)$ go through $\mathrm{Sp}^0_{\pm1}(2n)$, the color changes. The border lines from left to right are $$\Gamma_2(1)=\Gamma_3(1),\Gamma_3(-1),\Gamma_4(-1),\Gamma_4(1)=\Gamma_5(1),\Gamma_5(-1)\thickapprox\Gamma_6(-1),\Gamma_6(1)\cdots.$$
It can be shown that the 'green' region is a hyperbolic region which has positive measure. Moreover, whenever $\beta\mapsto \gamma_{2,\beta,\mathfrak{e}}(2\pi)$ passes $\mathrm{Sp}^0_{\pm1}(2n)$, the index $i_1(\gamma_{2,\beta,\mathfrak{e}})$ will increase by 2 and $i_{-1}(\gamma_{2,\beta,\mathfrak{e}})$ will increase by 1. Otherwise, $i_{\pm 1}(\gamma_{2,\beta,\mathfrak{e}})$ will be constance.



\section{Numerical results on rotation numbers} \label{sec: numerical}

As we mentioned the condition required in Theorem \ref{thm infy z brake} is difficult to verify and we are only able to show it is true for a smooth curve in the parameter space $\mathcal{D}$ (see Corollary \ref{cor case rot nequal}). Despite of this, we suspect the condition holds for a much larger subset of the parameter space $\mathcal{D}$. To support such a speculation, in this section we compute the rotation numbers of the Euler orbit $\zeta_e$ and the $z$-symmetric brake orbit $\zeta_{bz}$ numerically and show they are different, when the parameters satisfy
$$ \beta=\frac{\alpha}{4+\alpha}=0.6, \; \mathfrak{e}=\sqrt{1+2h\varpi^2\beta^2} =0.6,$$
which corresponds to
$$ m= 2^{\frac{1}{5}}, \; \alpha=6, \; \varpi=\sqrt{5/3}, \; h = -8/15. $$





We first compute $\hat i(\zeta_e)$ and $\rho(\zeta_e)$ numerically. On the one hand, according to the simplified linear system (\ref{linear euler z theta}) and the period of $\zeta_e$, i.e. $$T_e=2\pi\sqrt{\beta}(\frac{\varpi^2\beta}{1-\mathfrak{e}^2})^{\frac{3}{2}}=2\pi\sqrt{0.6}\cdot 0.8^{-3}=9.5057,$$
we obtain the following numerical picture of curve $a(t)=(r(t),\theta(t),z(t)),\ t\in[0,T_e]$, which satisfies
$$\gamma_{2,0.6,0.6}(t)=\begin{pmatrix}r(t)&z(t)\\ z(t)&\frac{1-z(t)^2}{r(t)}\end{pmatrix}\begin{pmatrix}\cos\theta(t)&-\sin\theta(t)\\ \sin\theta(t)&\cos\theta(t)\end{pmatrix}\in \mathrm{Sp}(2).$$
Since $a(t)$ intersects five times with the singular set $\mathrm{Sp}^0_1(2)$ (see section \ref{sec: AppIndex}) including $a(0)$, we have $i_1(\gamma_{2,0.6,0.6})=5$.
\begin{figure}[H]
\centering
\includegraphics[width=0.6\textwidth]{fig_of_fund_solu_Euler2.pdf}
\caption{The cylindrical coordinate $(r(t),\theta(t),z(t)),\ t\in [0,T_e]$.}
\end{figure}
Furthermore, we obtain $-0.9437\pm0.3307\sqrt{-1}\in\sigma(\gamma_{2,0.6,0.6}(T_e))$. Then by (\ref{eq: SumMeanIndex}) and Lemma \ref{mean index ga1},
$$\hat i(\gamma_e)=\hat i(\gamma_{2,0.6,0.6})+2=i_1(\gamma_{2,0.6,0.6})+1+\frac{\theta(T_e)}{\pi}=6.8928.$$

On the other hand, by computing the numerical solutions of (\ref{eq: RedHamSys}) with initial point in $\Sigma_-$, we can obtain the picture of return time function as Figure \ref{fig:t+}. Finally, by computing the position of discrete orbit $\{g^k(\zeta_e(0))\}$ in $\Sigma_-$ and using the definition of rotation number in Section \ref{sub: AnnulusMap}, we have
$$\rho(\zeta_e)=0.4090\thickapprox 0.4088=\frac{2}{\hat i(\gamma_{2,0.6,0.6})}=(\frac{\hat i(\gamma_e)}{2}-1)^{-1}$$ with 0.0002 error, which fits (\ref{rot mean index euler }) of Theorem \ref{thm: rot number and mean index gen} very well.

Now, we are in the position to compute $\hat i(\gamma_{bz})$ and $\rho(\zeta_{bz})$. As the proof of Theorem \ref{th: z symm brake}, we can find a $z$-symmetric brake orbit numerically as  Figure \ref{fig: CurvesB} with initial condition $x_{bz}=(0,-0.5234,1.4150,0)$ and period $T_{bz}=9.7067$.


On the one hand, by solving the linearization of (\ref{eq: RedHamSys}) at $\zeta_{bz}$, we can obtain $\gamma_{bz}\in\mathrm{Sp}(4)$ numerically. Since $\gamma_{bz}(T_{bz})$ has the normal form $N_1(1,1)\diamond R(\pm\theta)$, by a numerical computation, we get
 $$0.6354\pm0.7722\sqrt{-1},\ 0.9984\pm0.0564\sqrt{-1}\in\sigma(\gamma_{bz}(T_{bz})).$$
This implies $\theta=0.8823$. However, the eigenvector $\xi$ of $e^{-\sqrt{-1}\theta}$ satisfies $\langle \sqrt{-1}J\xi,\xi \rangle=-0.5958<0$, i.e. $e^{-\sqrt{-1}\theta}$ is Krein negative, we have
$$\gamma_{bz}(T_{bz})=N_1(1,1)\diamond R(-\theta).$$
Finally, by Lemma \ref{lem: mean index form}, the mean index $$\hat i(\gamma_{bz})=\frac{2\pi-\theta}{\pi}\ mod\ 2\mathbb{Z}=1.7189\ mod\ 2\mathbb{Z}.$$

On the other hand, we can also compute the rotation number of $\zeta_{bz}$ numerically as the Euler case and we have
$$\rho(x_{bz})=0.8599\thickapprox 0.8594 \ mod\ \mathbb{Z}=\frac{\hat i(\gamma_{bz})}{2}\ mod\ \mathbb{Z}$$ with error 0.0005, which also meets (\ref{rot mean index bz}) of Theorem \ref{thm: rot number and mean index gen} very well.

In conclusion, even we can not confirm that $(\beta,\mathfrak{e})\in\mathcal{D}_{\mathcal{R}}$ from the numerical results above, however, $\rho(x_{bz})\neq\rho(x_e)$ is obviously true, which supports our speculation at the beginning.


\section{Appendix: the index theory for symplectic paths} \label{sec: AppIndex}

In this section, we briefly review the index theory for symplectic path.  The detail  could be  found  in \cite{Lon02}. Let $(\mathbb{R}^{2n}, \Omega)$ be the standard symplectic vector space with coordinates $(x_1,\cdots,x_n, y_1,\cdots,y_n)$, then $\Omega=\sum_{i=1}^n dx_i\wedge dy_i$. Let $J=\left(\begin{array}{cc}0 & -I_{n} \\ I_{n} & 0\end{array}\right)$, where $I_n$ is the identity matrix on $\mathbb{R}^{n}$.
The symplectic group $\mathrm{Sp}(2n)$ is defined by
\[    \mathrm{Sp}(2n)=\{M\in GL(2n, \mathbb{R}) | M^TJM=J \},   \]
with topology induced from $\mathbb{R}^{4n^2}$.  For $\tau>0$ we are interested in paths in $\mathrm{Sp}(2n)$:
$$
\mathcal{P}_{\tau}(2 n)=\left\{\gamma \in C([0, \tau], \operatorname{Sp}(2 n)) \mid \gamma(0)=I_{2 n}\right\}.
$$
For any $\omega\in\mathbb{U}$(the unite circle in $\mathbb{C}$),
the following real function was introduced in  \cite{Lon02},
$$
D_{\omega}(M)=(-1)^{n-1} \bar{\omega}^{n} \operatorname{det}\left(M-\omega I_{2 n}\right).
$$
Thus for any $\omega\in\mathbb{U}$
 the following  $\omega$-degenerate hypersurface of codimension one in $\mathrm{Sp}(2n)$ is defined  \cite{Lon02}:
\[    \mathrm{Sp}(2n)_\omega^0=\{M\in\mathrm{Sp}(2n) | D_{\omega}(M)=0 \}.   \]
Moreover, the $\omega$-regular set of $\mathrm{Sp}(2n)$ is defined by
$$
\begin{gathered}
\operatorname{Sp}(2 n)_{\omega}^{\pm}=\left\{M \in \operatorname{Sp}(2 n) \mid \pm D_{\omega}(M)<0\right\}, \\
\operatorname{Sp}(2 n)_{\omega}^{*}=\operatorname{Sp}(2 n)_{\omega}^{+} \cup \operatorname{Sp}(2 n)_{\omega}^{-}
\end{gathered}.
$$
For any two continuous path $\xi$ and $\eta:[0, \tau] \rightarrow \operatorname{Sp}(2 n)$ with $\xi(\tau)=\eta(0)$, we define their concatenation as:
$$
\eta * \xi(t)= \begin{cases}\xi(2 t), & \text { if } 0 \leqslant t \leqslant \tau / 2 \\ \eta(2 t-\tau), & \text { if } \tau / 2 \leqslant t \leqslant \tau.\end{cases}
$$


Given any  two $2m_k\times2m_k$ matrices of square block form
$M_{k}=\left(\begin{array}{ll}A_{k} & B_{k} \\ C_{k} & D_{k}\end{array}\right)$ with $k=1, 2$, the $\diamond$-product of $M_1$ and $M_2$
 is defined by the following $2(m_1+m_2)\times 2(m_1+m_2)$ matrix  $M_{1} \diamond M_{2}$:
$$
M_{1} \diamond M_{2}=\left(\begin{array}{cccc}
A_{1} & 0 & B_{1} & 0 \\
0 & A_{2} & 0 & B_{2} \\
C_{1} & 0 & D_{1} & 0 \\
0 & C_{2} & 0 & D_{2}
\end{array}\right).
$$


We define a special continuous symplectic path $\xi_n\subset \mathrm{Sp}(2n)$ by
$$
\xi_{n}(t)=\left(\begin{array}{cc}
2-\frac{t}{\tau} & 0 \\
0 & \left(2-\frac{t}{\tau}\right)^{-1}
\end{array}\right)^{\diamond n} \quad \text { for } 0 \leqslant t \leqslant \tau.
$$

For $M\in \mathrm{Sp}(2n)_\omega^0$, we define a co-orientation of  $\mathrm{Sp}(2n)_\omega^0$
at $M$ by the positive direction $\frac{d}{dt}Me^{tJ}|_{t=0}$ of the path $Me^{tJ}$ with $|t|$ sufficiently small.
Now will give the definition of $\omega$-index \cite{Lon02}.

\begin{defi}\label{def:Maslov-type index} For $\omega\in \mathbb{U}$, $\gamma\in \mathcal{P}_{\tau}$, \textbf{the $\omega$-index} of $\gamma$  can be defined as
\[i_\omega(\gamma)=[e^{-\varepsilon J}\gamma*\xi_n: \mathrm{Sp}(2n)_\omega^0], \] For $\varepsilon>0$ small enough, where $[\cdot : \cdot]$ is the intersection number. We also denote $$\nu_\omega(\gamma)=\dim_{\mathbb{C}}\ker_{\mathbb{C}}(\gamma(\tau)-\omega I).$$
\end{defi}
Obviously, If $\gamma(\tau)\in\mathrm{Sp}(2n)^*_\omega$, we can take $\varepsilon=0$. Since $\xi_n(0)\in\mathrm{Sp}(2n)^+_\omega$, from \cite{Lon02}, we have
\begin{lem} \label{lem: index even}
$i_\omega(\gamma)$ is even if and only if  $e^{\varepsilon J}\gamma(\tau)\in\mathrm{Sp}(2n)^+_\omega$ for $\varepsilon>0$ small enough.
\end{lem}


For any two paths $\gamma_{j} \in \mathcal{P}_{\tau}\left(2 n_{j}\right)$ with $j=0, 1$, let $\gamma_{0} \diamond \gamma_{1}(t)=\gamma_{0}(t) \diamond \gamma_{1}(t)$ for all $t\in[0,\tau]$. It is obvious that
\begin{equation}
\label{eq:IndSum} i_{\omega}\left(\gamma_{1} \diamond \gamma_{2}\right)=i_{\omega}\left(\gamma_{1}\right)+i_{\omega}\left(\gamma_{2}\right), \quad \forall \omega \in \mathbb{U}.
\end{equation}


For any $M\in\mathrm{Sp}(2n)$ and $\omega\in\mathbb{U}$, the splitting
number $S^\pm_M(\omega)$ of $M$ at $\omega$ are defined by
\begin{equation}
\label{eq: SplitNumber} S_{M}^{\pm}(\omega)=\lim _{\varepsilon \rightarrow 0^{+}} i_{\omega e^{\pm \sqrt{-1} \varepsilon}}(\gamma)-i_{\omega}(\gamma),
\end{equation}
for any path  $\gamma \in \mathcal{P}_{\tau}(2 n)$ satisfying $\gamma(\tau)=M$.  The splitting number is well defined, i.e. they are independent of the choice of path $\gamma$.  From the property of Maslov-type index, if $M_0, M_1\in\mathrm{Sp}(2n)$ which are symplectic similar, i.e. $M_0\approx M_1$, then
$$S_{M_{0}}^\pm(\omega)= S_{M_{1}}^\pm(\omega). $$
For any $M_i\in\mathrm{Sp}(2n_i)$ with $i=0$ and $1$, there holds
$$
S_{M_{0} \diamond M_{1}}^{\pm}(\omega)=S_{M_{0}}^{\pm}(\omega)+S_{M_{1}}^{\pm}(\omega), \quad \forall \omega \in \mathbb{U}.
$$


We only list the splitting number of the normal form of $\mathrm{Sp}(2)$, which is enough for the current paper.
For $M\in\mathrm{Sp}(2)$, it is symplectic similar to one of the following  normal forms:
\begin{equation}\label{definition of R(theta)}
\begin{array}{c}
D(\lambda)=\left(\begin{array}{cc}
\lambda & 0 \\
0 & \lambda^{-1}
\end{array}\right), \quad \lambda\in \mathbb{R}\setminus\{0\} \\
N_{1}(\lambda, a)=\left(\begin{array}{cc}
\lambda & a \\
0 & \lambda
\end{array}\right), \quad \lambda=\pm 1, a=\pm 1,0 \\
R(\theta)=\left(\begin{array}{cc}
\cos \theta & -\sin \theta \\
\sin \theta & \cos \theta
\end{array}\right), \quad \theta \in(0, \pi) \cup(\pi, 2 \pi)
\end{array}
\end{equation}

\begin{lem}\label{split norm} (\cite{Lon02}) For $M\in\mathrm{Sp}(2n)$ and $\omega\in\mathbb{U}$, $0<\theta<\pi$, there hold
\begin{equation}
\begin{aligned}
S_{M}^{\pm}(\omega) &=0, \quad \text { if } \omega \notin \sigma(M), \\
S_{M}^{\pm}(\omega) &=S_{M}^{\mp}(\bar{\omega}), \\
0 \leqslant S_{M}^{\pm}(\omega) & \leqslant \nu_\omega(M), \\
S_{M}^{+}(\omega)+S_{M}^{-}(\omega) & \leqslant \operatorname{dim} \operatorname{ker}(M-\omega I)^{2 n}, \quad \text { if } \omega \in \sigma(M), \\
\left(S_{N_{1}(1, a)}^{+}(1), S_{N_{1}(1, a)}^{-}(1)\right) &=\left\{\begin{array}{ll}
(1,1), & \text { if } a=0,1, \\
(0,0), & \text { if } a=-1,
\end{array}\right.\\
\left(S_{N_{1}(-1, a)}^{+}(-1), S_{N_{1}(-1, a)}^{-}(-1)\right) &=\left\{\begin{array}{ll}
(1,1), & \text { if } a=0,-1, \\
(0,0), & \text { if } a=1,
\end{array}\right.\\
\left(S_{R(\theta)}^{+}\left(e^{\sqrt{-1} \theta}\right), S_{R(\theta)}^{-}\left(e^{\sqrt{-1} \theta}\right)\right) &=(0,1) \\
\left(S_{R(2 \pi-\theta)}^{+}\left(e^{\sqrt{-1} \theta}\right), S_{R(2 \pi-\theta)}^{-}\left(e^{\sqrt{-1} \theta}\right)\right) &=(1,0),
\end{aligned}
\end{equation}


\end{lem}

For any symplectic path $\gamma\in\mathcal{P}_{\tau}(2n)$ and $m\in\mathbb{Z}^+$, we define the $m$-th iteration $\gamma^m: [0,m\tau]\to\mathrm{Sp}(2n)$ by
$$
\gamma^{m}(t)=\gamma(t-j \tau) \gamma(\tau)^{j}, \quad \text { for } j \tau \leqslant t \leqslant(j+1) \tau, j=0,1, \ldots, m-1.
$$

\begin{defi} \label{def:mean index} The \textbf{mean index} of $\gamma$ is defined by
\beg{eqnarray} \hat{i}(\gamma)=\lim_{m\to\infty}\frac{i_1(\gamma^m)}{m}. \end{eqnarray}

\end{defi}

Meanwhile the Bott-type iteration formula shows
$$
i_{z}\left(\gamma^{m}\right)=\sum_{\omega^{m}=z} i_{\omega}(\gamma), \;\; \forall z\in\mathbb{U},
$$



which implies
$$
\hat{i}(\gamma) =\frac{1}{2 \pi} \int_{0}^{2 \pi} i_{e^{\sqrt{-1} \theta}}(\gamma) d \theta.
$$

Y. Long has showed a more explicit form of the mean index by the Bott-type formula and the basic norm form of $\mathrm{Sp}(2n)$. We give this form as follows.
\begin{lem}\label{lem: mean index form}(\cite{Lon02}) For any symplectic path $\gamma\in\mathcal{P}_{\tau}(2n)$ and $M=\gamma(\tau)$, we have
\begin{equation}
\hat{i}(\gamma)=i(\gamma)+S_{M}^{+}(1)-C(M)+\sum_{\theta \in(0,2 \pi)} \frac{\theta}{\pi} S_{M}^{-}\left(e^{\sqrt{-1} \theta}\right), 
\end{equation}
where  $C(M)=\sum_{0<\theta<2 \pi} S_{M}^{-}\left(e^{\sqrt{-1} \theta}\right)$.

\end{lem}

In the case $\gamma\in\mathcal{P}_{\tau}(2)$, from \cite{Lon02}, we have

\begin{equation} \label{mean index 2}
\hat{i}(\gamma)=\left\{\begin{array}{ll}
i_1(\gamma)-1+\theta/\pi, & \text { if } \gamma(\tau)\approx R(\theta), \theta \in(0, \pi) \cup(\pi, 2 \pi)  \\
i_1(\gamma)+1, & \text { if }\gamma(\tau)\approx N_1(1,a), a=1,0 \\ i_1(\gamma), & other\quad cases
\end{array}\right.
\end{equation}
Since $i_1(\gamma)$ is even only in the cases $\gamma(\tau)\approx N_1(1,-1)$ or $\gamma(\tau)\approx D(\lambda)$ with $\lambda>0$, we have
\begin{equation} \label{mean index mod 2}
\hat{i}(\gamma)/2\quad (mod \quad \mathbb{Z})=\left\{\begin{array}{ll}
\theta/2\pi, & \text { if } \gamma(\tau)\approx R(\theta), \theta \in(0, \pi) \cup(\pi, 2 \pi),  \\
1/2, & \text { if } \sigma(\gamma(\tau))\in\mathbb{R}^-, \\ 0, & \text { if } \lambda\in\sigma(\gamma(\tau))\quad with \quad \lambda>0. \end{array} \right.
\end{equation}

Now we brief review the Morse index theorem.
Consider the second order system \beg{eqnarray} \ddot{x}=R(t)x, \label{second eq}\end{eqnarray}
where  $R(t+T)=R(t)$ is continuous path of symmetric  $n\times n$ matrices.   Let
$$\mathcal{A}=-\frac{d^2}{dt^2}+R(t)$$ be the corresponding  operator, it is self-adjoint operator on $L^2([0,T], \mathbb{C}^n)$ with self-adjoint boundary conditions. We always denote $\mathcal{M}^0(\mathcal{A})=\dim\ker(\mathcal{A})$ and the Morse index $\mathcal{M}^-(\mathcal{A})$ is the total number of negative eigenvalues of $\mathcal{A}$.

Obviously, the Morse index depends on the boundary conditions. We list the three kinds of condition:

i)  $\omega$-boundary conditions for $\omega\in\mathbb{U}$,
$\Lambda_\omega:=\{ x(0)=\omega x(T),\ \dot{x}(0)=\omega\dot{x}(T)  \}$;

ii) Dirichlet boundary conditions, $\Lambda_D:=\{ x(0)=0,\ x(T)=0 \}$;

iii) Neumann boundary conditions,  $\Lambda_N:=\{ \dot{x}(0)=0,\ \dot x(T)=0 \}$.

By the standard Legendre transformation, the system \eqref{second eq} becomes in to
\beg{eqnarray}  \dot{z}(t)=J\mathcal{B}(t)z(t), \label{ham equ} \end{eqnarray} where  $\mathcal{B}(t)=\left[\begin{array}{cc}I_n & 0_n \\ 0_n & -R(t)\end{array}\right]$.  Then $\Lambda_1=\mathbb{R}^n\oplus\{0\}$,  $\Lambda_2=\{0\}\oplus\mathbb{R}^n$ are Lagrangian spaces of $(\mathbb{R}^{2n}, \Omega)$.

Let $\gamma(t)$ be the fundamental solution of  \eqref{ham equ}, the well known Morse index theorem shows that
\begin{equation}
 \mathcal{M}^{-}\left(\Lambda_{D}\right)=\sum_{0<t_{0}<T} \operatorname{dim} \gamma\left(t_{0}\right) \Lambda_1 \cap \Lambda_1. \label{morse D}
\end{equation}
For the $\omega$ boundary conditions, by \cite{Lon02}
\begin{equation}
 \mathcal{M}^{-}\left(\Lambda_{\omega}\right)=i_\omega(\gamma). \label{more equa mas omega} \end{equation}
For the Neumann boundary conditions, $\mathcal{M}^{-}(\Lambda_{N})$ expressed by Maslov index, we will only introduce a simple case here which is enough for the current paper. The path $\gamma$ is called positive path if $\mathcal{B}(t)>0$ for $t\in[0,T]$.  In this case, we have
\begin{equation}
 \mathcal{M}^{-}\left(\Lambda_{N}\right)=\sum_{0\leq t_{0}<T} \operatorname{dim} \gamma\left(t_{0}\right) \Lambda_2 \cap \Lambda_2. \label{morse N p}
\end{equation}

From \cite{HWY20}, we have
\begin{equation}
\mathcal{M}^{-}\left(\Lambda_{D}\right)\leq \mathcal{M}^{-}\left(\Lambda_{\omega}\right)\leq \mathcal{M}^{-}\left(\Lambda_{N}\right)\leq \mathcal{M}^{-}\left(\Lambda_{D}\right)+2n. \label{index estimate}\end{equation}

From \eqref{morse D} and \eqref{index estimate}, we have
\beg{eqnarray} \hat{i}(\gamma)=\lim_{k\to \infty}\frac{1}{k}\sum_{0<t_{0}<kT} \operatorname{dim} \gamma\left(t_{0}\right) \Lambda_1 \cap \Lambda_1. \label{mean index D}  \end{eqnarray}
Similarly, for positive path $\gamma$, we have
\beg{eqnarray} \hat{i}(\gamma)=\lim_{k\to \infty}\frac{1}{k}\sum_{0\leq t_{0}<kT} \operatorname{dim} \gamma\left(t_{0}\right) \Lambda_2 \cap \Lambda_2. \label{mean index N} \end{eqnarray}



\hfill\newline
\noindent{\bf Acknowledgement.}
The authors thank Prof. Jian Wang for useful discussion about reversible maps.





\bibliographystyle{abbrv}

