# KAM and geodesic dynamics of blackholes

Jinxin Xue  Jingzhai 310, Department of mathematics and Yau Mathematical
Sciences Center, Tsinghua University, Beijing,100084<EMAIL_ADDRESS>
###### Abstract.

In this paper we apply KAM theory and the Aubry-Mather theory for twist maps
to the study of bound geodesic dynamics of a perturbed blackhole background.
The general theories apply mainly to two observable phenomena: the photon
shell (unstable bound spherical orbits) and the quasi-periodic oscillations.
We discover there is a gap structure in the photon shell that can be used to
reveal information of the perturbation.

###### Contents

  1. 1 Introduction
    1. 1.1 The photon shell dynamics
    2. 1.2 Quasi-periodic oscillations (QPO)
  2. 2 Introduction to nearly integrable Hamiltonian dynamics
    1. 2.1 Liouville-Arnold theorem
    2. 2.2 The Kolmogorov-Arnold-Moser theorem
    3. 2.3 Energetic reduction
    4. 2.4 Dynamics of twist maps
  3. 3 The Schwarzschild spacetime
    1. 3.1 Schwarzschild dynamics
      1. 3.1.1 The radial dynamics, critical points and homoclinic orbits
      2. 3.1.2 The normally hyperbolic invariant manifold
    2. 3.2 Dynamics on the photon shell
    3. 3.3 KAM and QPO for Schwarzschild
      1. 3.3.1 KAM and twist maps theory for the QPO
      2. 3.3.2 KAM nondegeneracy and further derivatives
      3. 3.3.3 Three fundamental frequencies
  4. 4 The Kerr spacetime
    1. 4.1 The vertical dynamics
    2. 4.2 The photon shell and the radial dynamics
      1. 4.2.1 Determining the photon shell
      2. 4.2.2 The Hamiltonian restricted to the photon shell
      3. 4.2.3 The radial dynamics
    3. 4.3 KAM and the photon shell dynamics
      1. 4.3.1 KAM and twist map nondegeneracy
      2. 4.3.2 The point of discontinuity
    4. 4.4 KAM and QPO
  5. A The theorem of normally hyperbolic invariant manifold
  6. B Aubry-Mather theory for twist maps

##  1\. Introduction

This paper is a companion paper of [X]. The main purpose is to study the
stability of the geodesic dynamics of a blackhole background under
perturbations.

The important metrics such as Schwarzschild and Kerr are both integrable, each
with four independent conserved quantities. If we consider the bounded motions
in an integrable system, it is known from Liouville-Arnold theorem that the
dynamics can be treated as toral linear flows. When an integrable system is
perturbed slightly, the integrability is generically broken and chaotic
motions occur. The remarkable discovery of the Kolmogorov-Arnold-Moser theorem
shows that most motions are still toral linear flows after a coordinate
change. In reality, the blackholes modeled on either Schwarzschild and Kerr
always undergo some perturbations. We try to understand the dynamics of a
particle, massless or mass, moving in a perturbed Schwarzschild or Kerr
background, using the theory of Hamiltonian dynamical systems developed for
nearly integrable systems.

The mathematical model that we use to study the geodesic dynamics of the
blackholes is the following map called integrable twist map:

(1.1) $\phi_{0}:\
\mathbb{T}^{1}\times[0,1]\to\mathbb{T}^{1}\times[0,1],\quad(\theta,I)\mapsto(\theta+\nu(I),I)$

where $\nu$ is assumed to be smooth and strictly monotone in $I$. The dynamics
of this example is as follows. Each $I$-circle is invariant and the dynamics
on it is a rotation by $\nu(I)$. We call this rotation rate $\nu(I)$ the
rotation number. When $\nu(I)$ is rational, then each orbit on the $I$-circle
is periodic, and when $\nu(I)$ is irrational, then each orbit is dense on the
$I$-circle. When $\phi_{0}$ is slightly perturbed within the class of smooth
symplectic maps, we have the following picture for the dynamics. First,
Moser’s version of KAM theory implies that each $I$-circle with $\nu(I)$ a
Diophantine number is perturbed to a nearby smooth loop, called invariant
circle, that remains invariant under the perturbed map and the dynamics on it
can be conjugate to the original unperturbed rotation, provided the
perturbation is sufficiently small in a smooth enough topology. For a fixed
small perturbation of size $\varepsilon$, the set of remaining invariant
circles has rotation numbers lying in a Cantor set (a nowhere dense closed
set) of measure $1-O(\sqrt{\varepsilon})$. Between two nearby invariant
circles (not necessarily KAM curve), there is a gap region called Birkhoff
instability region where the dynamics is very chaotic. There is a Aubry-Mather
theory developed for general twist maps which gives the existence of some
special orbits in the Birkhoff instability region, such as periodic orbits,
heteroclinic orbits corresponding to rational rotation numbers and Cantor set
like orbit corresponding to irrational rotation numbers. So KAM theory allows
us to find a large measure set of the phase space where the motion is regular
(quasiperiodic), while Aubry-Mather theory allows use to find some special
orbits in the gaps. We give an outline of the two theories in Section 2.4 and
more details in Appendix B.

We first locate the part of phase space with bounded motions. For massless
particles moving on null geodesics in Schwarzschild or Kerr background, these
bounded motions are called bound photon orbits, or fundamental photon orbits
in literature. Each such orbit is moving on a sphere with a fixed radius and
is unstable under radial perturbations. This is an observable feature of the
blackhole, which lies on the edges of the blackhole shadows [M87, N]. If we
also consider the motion of massive particles on timelike geodesics in
Schwarzschild or Kerr background, again we have unstable bound spherical
orbits, similar to the lightlike case. For simplicity, we use the term photon
shell to call the set of bound spherical orbits that are unstable under radial
perturbations, for both Schwarzschild and Kerr, and both lightlike or timelike
geodesics. The photon shell when considered in the phase space corresponds to
a geometric object called normally hyperbolic invariant manifold in dynamical
systems, which is known to be robust under perturbations. This means that an
orbit on the photon shell remains a bound orbit in a perturbed Schwarzschild
or Kerr metric, though it may become chaotic. When a background is a
perturbation of Schwarzschild or Kerr, we still use the terminology photon
shell to call the union of bound orbits that are unstable under radial
perturbations. Note that in this case, the radial component of each orbit may
oscillate slightly.

When studying Arnold diffusion, people often need some information on the
normally hyperbolic invariant manifold, which in our case is the photon shell
(c.f. [X]). We show in this paper that both KAM and the twist map theory apply
to the study of the dynamics on the photon shell in a perturbed Kerr
spacetime, but neither applies to perturbed Schwarzschild a priori without
further information of the perturbation. In particular, KAM implies in the
Kerr case the existence of a gap structure (many Birkhoff instability regions)
on the photon shell exhibiting coexistence of regular and chaotic motions,
which we believe is observable and can reveal information of the perturbation.

Furthermore, in the massive case, in addition to the photon shell, there are a
lot more bound timelike geodesics. It can be observed that outside a blackhole
that is a slowly rotating accretion disk emitting X-rays exhibiting some
quasi-periodic oscillations (QPO) in its frequency (c.f. [DKN, IM, JDW] etc).
In the phase space this is a region around a stable circular orbit. We show
both KAM and twist map theory apply in this region for both Schwarzchild and
Kerr cases to yield lots of quasiperiodic orbits as well as other special
orbits.

Both Schwarzschild and Kerr metrics are written in coordinates
$(\tau,r,\theta,\varphi)$ where $\tau$ is the coordinate time, $r$ is the
polar radius, $\theta$ is the latitudinal angle and $\varphi$ is the azimuthal
angle. We will use the following definition throughout the paper.

######  Definition 1.1.

A metric is called stationary if it is independent of the coordinate time
$\tau$ and axisymmetric if it is independent of the azimuth angle $\varphi$.

###  1.1. The photon shell dynamics

We consider perturbations stationary perturbation so that the particle’s
energy $E$, that is dual to $\tau,$ remains to be a constant of motion, so we
fix a value $E$ and restrict to the photon shell. The resulting system has two
degrees of freedom with coordinates $(\theta,p_{\theta},\varphi,L_{z})$ and we
only need to study the latitudinal and azimuthal motions, where $p_{\theta}$
and $L_{z}=p_{\varphi}$ are generalized momentum variables due to $\theta$ and
$\varphi$ respectively, and $L_{z}$ in particular is a constant of motion
having the physical meaning of the $z$-component of the angular momentum. This
subsystem describes a particle’s motion on the photon shell that has a
vertical oscillation within a spherical band symmetric around the equator when
orbiting around the center of the blackhole.

Our general strategy is to introduce action-angle coordinates
$(\alpha_{\theta},J_{\theta},\alpha_{\varphi},L_{z})\in\mathbb{T}\times\mathbb{R}_{+}\times\mathbb{T}\times\mathbb{R}$
for the $\theta$\- and $\varphi$-motions. In these coordinates the unperturbed
Schwarzschild or Kerr Hamiltonian has the form $H(J_{\theta},L_{z})$, so the
dynamics leaves $J_{\theta},L_{z}$ invariant and the angles
$(\alpha_{\theta},\alpha_{\varphi})$ move on the torus $\mathbb{T}^{2}$
linearly as
$(\alpha_{\theta},\alpha_{\varphi})\mapsto(\alpha_{\theta},\alpha_{\varphi})+t(\partial_{J_{\theta}}H,\partial_{L_{z}}H)\
\mathrm{mod}\ \mathbb{Z}^{2},\ t\in\mathbb{R}$. For this system, we further
fix a Hamiltonian level set $H=-\frac{1}{2}\mu^{2}$ with $\mu^{2}=0$ in the
massless case and $\mu^{2}=1$ in the massive case, and look at the return map
to the section $\\{\alpha_{\theta}=0\\}$ (or to say $\\{\theta=\pi/2\\}$).
Then the map has the form of (1.1) with
$\nu:=\partial_{L_{z}}H/\partial_{J_{\theta}}H$ parametrized by $L_{z}$. To
apply KAM theory and twist map theory, we have to find a parameter range where
the ratio has no critical point. The main difficulty is that the Hamiltonian
cannot be written explicitly as the action variables $(J_{\theta},L_{z})$.

The Schwarzschild case is relatively simple, since each photon spherical orbit
moves on a plane passing through the origin so is periodic and all the photon
spherical orbits have the same period. In this case the ratio $\nu$ is
constantly 1, so neither KAM nor twist map theory applies without further
knowledge of the perturbations.

It turns out that both theories apply to Kerr. For Kerr, we can reparametrize
the ratio $\nu$ by the radius $r$ of the photon shell orbits and denote by
$\mathcal{R}$ the domain of definition for $r$ (c.f. Definition 4.2). So the
photon shell of the unperturbed Kerr is foliated by $S_{r}$, that is a
spherical strip symmetric around the equator consisting of all spherical
orbits of radius $r\in\mathcal{R}$.

######  Theorem 1.2.

Consider null geodesics in a Kerr spacetime with mass $M$ and angular momentum
$a$. Then there exists an open set of parameter $a/M$ such that the frequency
ratio $\nu=\omega_{\varphi}/\omega_{\theta}$ as a function of
$r\in\mathcal{R}$ is smooth, strictly monotone and has a discontinuity point
at some $r_{\star}\in\mathcal{R}$ with $\nu(r_{\star-})-\nu(r_{\star+})=2$.
Orbits on $S_{r}$ with $r<r_{\star}$ are all prograde and that with
$r>r_{\star}$ are retrograde. Denote by $U$ a small neighborhood of
$r_{\star}$.

Let $\varepsilon h_{\mu\nu}dx^{\mu}dx^{\nu}$ be a stationary perturbation of
the above Kerr metric. Then if $\varepsilon$ sufficiently small, we have the
following for the dynamics of null geodesics in the perturbed system

  1. (1)

There is a closed subset $\mathcal{C}$ of $\mathcal{R}\setminus U$ such that
for each $r\in\mathcal{C}$, in the perturbed system, there is a set
$\tilde{S}_{r}$ that is a small deformation of $S_{r}$, such that
$\tilde{S}_{r}$ consists of periodic or quasiperiodic orbits.

  2. (2)

Each interval in $\mathcal{R}\setminus(U\cup\mathcal{C})$ corresponds to a
Birkhoff instability region where the spheres $S_{r}$ are broken (c.f.
Definition B.2) and twist map theory $($c.f. Theorem B.1 $(2.b),(3.b)$ and
Theorem B.3$)$ applies.

Figure 1. The frequency ratio $\omega_{\varphi}/\omega_{\theta}$ with
parameters $a/M=0.8,\mu^{2}=0$

In particular, the gaps (instability regions) often occur around rational
frequency ratios of the form $p/q$, due to the existence of Fourier modes
$\exp(\sqrt{-1}k(p\alpha_{\theta}+q\alpha_{\varphi})),$ $k,p,q\in\mathbb{Z},$
in the perturbation. In general, the gap is large when the denominator of the
rational number is small. We expect that this is an observable phenomenon, in
some sense an analogue of the Kirkwood gap in the asteroid belt in the solar
system, and can be used to conclude the information of the perturbation.

In Figure 1, we see that there is a discontinuity in the frequency ratio
$\nu$, which occurs at a point where $L_{z}$ changes sign and the motion of
the particle becomes from prograde to retrograde. This phenomenon was observed
in [Wi, T]. We give a elementary mathematical explanation in Section 4.3.2. In
the proof, we deal with a particularly chosen parameter value $a/M=0.8$, but
it is clear that our formulas are applicable to any other parameter values of
interest.

For the geodesic motion of massive particles, we have a similar theorem (c.f.
Theorem 4.4). The main difference of the massive case and the massless case is
that the dynamics on the photon shell of the massive case depends on the
particle’s $E$, while the massless case does not. So in the massive case, we
need to fix an $E$ ahead of time, then there is the similar gap structure for
the photon orbits of all massive particles with the same energy $E$.

Note that our result can be combined with the result of Johnson, Lupsasca,
Strominger et al [J] to give a delicate picture of the fine structure of the
neighborhood of the photon shell. It is called the photon ring the image on
the observer’s screen of photon’s on nearly bound null geodesics. So the
photon ring contains the photon shell. It was discovered in [J] that there is
a subring structure outside the photon shell. Indeed, if we shoot a light ray
very near, a distance from the shadow edge at $\delta r_{0}$, it will circle
$n\sim-\frac{1}{\gamma}\ln|\delta r_{0}|$ times before falling into the black
hole or escaping to infinity, due to the instability in the direction normal
to the photon shell, where $\gamma$ is the Lyapunov exponent. A light ray that
completes $n$ half orbits collects $n$ times more photons along its path. This
gives a ring on the screen subdivided into subrings labeled by the number of
circulations $n$, exponentially narrower when approaching the photon shell.

###  1.2. Quasi-periodic oscillations (QPO)

In the case of massive particles, the geodesic dynamics has lots of bounded
motions other than the bound spherical orbits. These bounded motions occur
around the stable spherical orbits, thus they undergo both radial oscillations
and the latitudinal oscillations when spinning around the blackhole. They are
supposed to be responsible for the QPO behavior in some theories.

We study the dynamics of these bound orbits using similar approaches as
before. For simplicity, we consider perturbations that are both stationary and
axisymmetric hence $E$ and $L_{z}$ remain to be constants of motion. Then we
can use coordinates $(r,p_{r},\theta,p_{\theta})$ to study the radial and
latitudinal dynamics of the Hamiltonian system. We introduce action-angle
coordinates
$(\alpha_{r},J_{r},\alpha_{\theta},J_{\theta})\in\mathbb{T}\times\mathbb{R}_{+}\times\mathbb{T}\times\mathbb{R}_{+}$
so that the Hamiltonian is a function of $J_{r}$ and $J_{\theta}$ only and the
radial frequency $\partial_{J_{r}}H$ and vertical frequency
$\partial_{J_{\theta}}H$ are called fundamental frequencies in [S]. The ratio
$\nu=\partial_{J_{r}}H/\partial_{J_{\theta}}H$ plays a similar role to the
photon shell case as above. Again we can compute $\nu$ to verify that the
assumptions for both KAM and twist map theory apply, so we have an analogue of
Theorem 1.2 in both Schwarzschild and Kerr cases. We refer readers to Theorem
3.2 and 4.5 for more details.

If we consider non axisymmetric perturbations, then we should consider a
Hamiltonian system of three degrees of freedom with coordinates
$(r,p_{r},\theta,p_{\theta},\varphi,L_{z})$. The KAM nondegeneracy condition
can be studied similarly. We shall outline how this can be done without
working on details. Action-angle coordinates for Kerr has been constructed in
[S]. In general, we do not expect that KAM nondegeneracy condition holds
globally in the part of phase space with bounded motions for Kerr, considering
the isofrequency pairing results (c.f. [WBS] etc).

Organization of the paper The paper is organized as follows. The main body of
the paper consists mainly of introductory and conceptual arguments and
statements, and we put all the technical ingredients to the appendix. In
Section 2, we give a brief introduction to nearly integrable Hamiltonian
dynamics including Liouville-Arnold theorem, KAM theorem, and twist map
theory. In Section 3, we study the perturbed Schwarzschild dynamics. In
Section 4, we study the dynamics of perturbed Kerr. Finally, we have two
appendices. In Appendix A we introduce the theorem of normally hyperbolic
invariant manifolds, and in Appendix B, we introduce the Aubry-Mather theory
of twist maps.

##  2\. Introduction to nearly integrable Hamiltonian dynamics

In the Hamiltonian formalism of the classical mechanics, a smooth Hamiltonian
function $H$ on a symplectic manifold $(M,\omega)$ is given, and defines a
vector field $X_{H}$ through $\omega(\cdot,X_{H})=dH$ which defines a
dynamical system by solving the ODE $\dot{x}=X_{H}(x)$.

###  2.1. Liouville-Arnold theorem

The classical Liouville-Arnold theorem gives important characterization of the
integrable systems.

######  Theorem 2.1 (Liouville-Arnold).

Let $H_{1}=H:\ M^{2n}\to\mathbb{R}$ be a Hamiltonian and suppose there are
$H_{2},\ldots,H_{n}:\ M\to\mathbb{R}$ satisfying

  1. (a)

$\\{H_{i},H_{j}\\}\equiv 0$, for all $i,j=1,\ldots,n$, where
$\\{H_{i},H_{j}\\}:=\omega(X_{H_{i}},X_{H_{j}})$ is the Poisson bracket;

  2. (b)

the level set $M_{\mathbf{a}}:=\\{(q,p)\in M\ |\ H_{i}(q,p)=a_{i},\
i=1,\ldots,n\\}$ is compact;

  3. (c)

At each point of $M_{\mathbf{a}}$, the $n$ vectors $DH_{i},\ i=1,\ldots,n$ are
linearly independent.

Then

  1. (1)

$M_{\mathbf{a}}$ is diffeomorphic to
$\mathbb{T}^{n}=\mathbb{R}^{n}/\mathbb{Z}^{n}$ and is invariant under the
Hamiltonian flow of each $H_{i}$.

  2. (2)

in a neighborhood $U$ of $M_{\mathbf{a}}$, there is a symplectic transform
$\Phi(q,p)=(\theta,I)$ such that
$\Phi(U)=\mathbb{T}^{n}\times(-\delta,\delta)^{n}$ for some $\delta>0$.

  3. (3)

In the new coordinates, each $K_{i}:=H_{i}\circ\Phi^{-1}$ is a function of $I$
only so the Hamiltonian equation is

$\dot{\theta}=\omega_{i}(I):=\frac{\partial K_{i}}{\partial
I},\quad\dot{I}=0.$

We refer readers to Section 50 of [A89] for the proof and more details. The
coordinates $(\theta,I)$ are called action-angle coordinates and can be
constructed by standard recipe in [A89].

###  2.2. The Kolmogorov-Arnold-Moser theorem

The Kolmogorov-Arnold-Moser theorem is an important stability result on the
dynamics of nearly integrable systems.

######  Definition 2.2.

A vector $\alpha\in\mathbb{R}^{n},\ n\geq 2$ is called _Diophantine_ if there
exist some $C>0,\tau>n-1$ such that for all
$k\in\mathbb{Z}^{n}\setminus\\{0\\}$, we have

$|\langle\alpha,k\rangle|\geq\frac{C}{|k|^{\tau}}.$

A number $\alpha\in\mathbb{R}$ is called a Diophantine number, if
$(1,\alpha)\in\mathbb{R}^{2}$ is a Diophantine vector. In both cases, we
denote $\alpha\in\mathrm{DC}(C,\tau)$.

For each $\tau>n-1$, the set $\cup_{C>0}\mathrm{DC}(C,\tau)$ has full Lebesgue
measure in $\mathbb{R}^{n}$. For a fixed $C$, the set $\mathrm{DC}(C,\tau)$ is
a Cantor set of positive measure.

A version of the KAM theorem states as follows.

######  Theorem 2.3.

Suppose $\det D^{2}h(I)\neq 0$ for $I\in D$. Then for all $C>0,\tau>n-1$ and
all $I\in D$ with $\omega(I)\in\mathrm{DC}(C,\tau)$, there exist
$\varepsilon_{0}>0$ and $\ell$ sufficiently large such that when
$\varepsilon\|f\|_{C^{\ell}}<\varepsilon_{0}$, we have that there is a torus
$\mathcal{T}_{I}$ that is invariant under the Hamiltonian flow generated by
the Hamiltonian $H(\theta,I)=h(I)+\varepsilon f(\theta,I)$. Moreover, there is
a symplectic transform $\Phi$ whose restriction to $\mathcal{T}_{I}$ satisfies
$\Phi(\mathcal{T}_{I})=\mathbb{T}^{n}$ and the conjugated flow
$\Phi\circ\phi^{t}_{H}\circ\Phi^{-1}$ is the translation
$\theta\mapsto\theta+\omega(I)t$ when restricted to $\mathbb{T}^{n}$.

We refer readers to [P] for more details. When applying KAM theory to an
energy level, we may replace the nondegeneracy condition $\det D^{2}h(I)\neq
0$ by Arnold’s isoenergetic nondegeneracy condition
$\det\left[\begin{array}[]{cc}D^{2}h(I)&Dh(I)\\\
(Dh(I))^{t}&0\end{array}\right]\neq 0$ (c.f. Appendix 8 of [A89] and Chapter
6.3 of [AKN] for various versions of nondegeneracy conditions).

###  2.3. Energetic reduction

Let $H:\ T^{*}\mathbb{T}^{2}\to\mathbb{R}$ be a Hamiltonian system of two
degrees of freedom of the form $H(x_{1},x_{2},y_{1},y_{2})$ where
$(x_{1},x_{2})\in\mathbb{T}^{2}$ and $(y_{1},y_{2})\in\mathbb{R}^{2}$. Suppose
we have $\frac{\partial H}{\partial y_{1}}\neq 0,$ then we can fix an energy
level $E$ and solve the equation $H(x_{1},x_{2},y_{1},y_{2})=E$ for
$y_{1}=y_{1}(x_{1},x_{2},y_{2})$. We can then treat $-y_{1}$ as a
nonautonomous Hamiltonian defined on $T^{*}\mathbb{T}^{1}\times\mathbb{T}^{1}$
where $x_{1}\in\mathbb{T}^{1}$ plays the role of new time. By implicit
function theorem, we have the Hamiltonian equations

(2.1) $\begin{cases}\frac{dx_{2}}{dx_{1}}&=\frac{\partial H}{\partial
y_{2}}\Big{/}\frac{\partial H}{\partial
y_{1}}=\frac{\partial(-y_{1})}{\partial y_{2}},\\\
\frac{dy_{2}}{dx_{1}}&=-\frac{\partial H}{\partial x_{2}}\Big{/}\frac{\partial
H}{\partial y_{1}}=-\frac{\partial(-y_{1})}{\partial y_{2}}.\end{cases}$

This procedure is called the energetic reduction (Section 45 of [A89]).

###  2.4. Dynamics of twist maps

Twist map is defined as follows.

######  Definition 2.4.

A diffeomorphism $\phi:\mathbb{T}^{1}\times[0,1]\to\mathbb{T}^{1}\times[0,1]$
preserving the boundary is called an area preserving twist map if the
following holds:

  1. (1)

$\phi$ is symplectic, i.e. $\phi^{*}\omega=\omega,\ \omega=dx\wedge dy;$

  2. (2)

for each $x\in\mathbb{T}$, the map $\pi_{1}\phi:\
\\{x\\}\times[0,1]\to\mathbb{T}^{1}$ is a local homeomorphism, where $\pi_{1}$
means the projection to the first $(\mathbb{T}^{1})$ component.

One example of a twist map is the map (1.1). There is a version of KAM theorem
developed for twist maps by Moser [M2] that we shall mainly apply later in the
paper.

######  Theorem 2.5.

Let $\phi_{0}:\ \mathbb{T}^{1}\times[0,1]\to\mathbb{T}^{1}\times[0,1]$ be a
twist map as in (1.1) and $\phi_{\varepsilon}:=\phi_{0}+\varepsilon f:\
\mathbb{T}^{1}\times[0,1]\to\mathbb{T}^{1}\times[0,1]$ be a perturbation of
$\phi_{0}$. For all $C>0,\tau>1$, and all $I\in[0,1]$ with
$\nu(I)\in\mathrm{DC}(C,\tau)$, there exists $\varepsilon_{0}>0$ and $\ell$
sufficiently large such that when $\|\varepsilon
f\|_{C^{\ell}}<\varepsilon_{0}$, we have that there is a homotopically
nontrivial circle $\gamma_{I}$ that is invariant under the map
$\phi_{\varepsilon}$ and the dynamics $\phi_{\varepsilon}:\
\gamma_{I_{*},\varepsilon}\to\gamma_{I_{*},\varepsilon}$ can be conjugated to
the rotation $\theta\mapsto\theta+\nu(I_{*})$ on $\mathbb{T}^{1}$.

Another important example of twist maps is the standard map defined as
follows: $\Phi:\ \mathbb{T}^{2}\to\mathbb{T}^{2}$ via

$(x,y)\mapsto(x+y+k\sin x,y+k\sin x)$

where $k$ is a parameter.

For a general nonlinear twist map (not necessarily nearly integrable), we can
introduce a variational principle and define globally minimizing orbits which
also has the notion of rotation number playing a similar role as $\nu(I)$ as
in (1.1). In Aubry-Mather theory, the globally minimizing orbits are
classified according to their rotation numbers, denoted by
$\mathcal{M}_{\rho}$. We give an overview of this theory in Appendix B.

The Aubry-Mather theory gives us the following picture for the dynamics of a
twist map. Let $\phi_{\varepsilon}$ be a smooth twist map and when
$\varepsilon=0$, the map $\phi_{0}$ is of the form (1.1). Then for
$\varepsilon$ sufficiently small, we have

  1. (1)

KAM implies that for Diophantine irrational rotation numbers ($\rho\in
DC(C,\tau)$ for $C>0$ independent of $\varepsilon$) the set
$\mathcal{M}_{\rho}$ is a invariant curve restricted on which the map $\phi$
is conjugate to an irrational rotation, i.e. there exists a diffeomorphism
$h:\ \mathcal{M}_{\rho}\to\mathbb{T}$ such that $h\phi=h+\rho.$

  2. (2)

For each irrational rotation number $\rho$, the set $\mathcal{M}_{\rho}$ is
either a homotopically nontrivial invariant curve, or a Denjoy minimal set.

  3. (3)

For rational rotation numbers $\rho$, it may happen that there exists a
homotopically nontrivial curve consists of periodic orbits with the same
rotation number $\rho\in\mathbb{Q}$. However, typically, such curve does not
exist and $\mathcal{M}_{\rho}$ consists of one single periodic orbit with two
homoclinic orbits approaching it in both the future and the past.

  4. (4)

Whenever there is a region bounded by two neighboring invariant circles, there
exist orbits crossing the gap and visiting any two neighborhoods of the
boundary circles.

  5. (5)

In a neighborhood of an elliptic periodic point of least period $q$
($D\phi^{q}$ has no eigenvalues off the unit circle), the above
(1),(2),(3),(4) apply also when the twist condition is satisfied. Such a
neighborhood is called an elliptic island.

##  3\. The Schwarzschild spacetime

###  3.1. Schwarzschild dynamics

The Schwarzschild metric is as follows

(3.1)
$ds^{2}=-\left(1-\frac{2M}{r}\right)d\tau^{2}+\frac{1}{1-\frac{2M}{r}}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\varphi^{2})=g_{\mu\nu}dx^{\mu}dx^{\nu}$

with $x^{0}=\tau$ the time coordinate, $x^{1}=r$ is the polar radius, and
$x^{2}=\theta\in[0,\pi],x^{3}=\varphi\in[0,2\pi)$ are the spherical
coordinates where $\theta$ is the latitudinal angle and $\varphi$ is the
azimuthal angle. Here $M$ is the mass of the blackhole, and $r=2M$ is the
event horizon. We are only interested in the dynamics of a particle moving
outside of the event horizon. In the following, we denote
$\alpha=\left(1-\frac{2M}{r}\right)$.

We view the metric as twice of a Lagrangian
$L(x,\dot{x})=\frac{1}{2}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}$ (here
$\dot{x}=\frac{dx}{dt}$ means the derivative with respect to the particle’s
proper time $t$) and obtain its conjugate Hamiltonian after a formal Legendre
transform as follows

(3.2) $\displaystyle\begin{cases}\frac{d\tau}{dt}&=-\alpha^{-1}p_{\tau},\\\
\frac{dr}{dt}&=\alpha p_{r},\\\
\frac{d\theta}{dt}&=\frac{p_{\theta}}{r^{2}},\\\
\frac{d\varphi}{dt}&=\frac{p_{\varphi}}{r^{2}\sin^{2}\theta},\\\ \end{cases}$
$\displaystyle 2H(x,p)$
$\displaystyle=2(\dot{x}^{\mu}p_{x^{\mu}}-L(x,\dot{x}))=g^{\mu\nu}p_{\mu}p_{\nu}$
$\displaystyle=-\alpha^{-1}p_{\tau}^{2}+\alpha
p_{r}^{2}+r^{-2}\left(p_{\theta}^{2}+\frac{1}{\sin^{2}\theta}p_{\varphi}^{2}\right).$

The Hamiltonian system has four conserved quantities as
$H,p_{\tau},p_{\varphi},L^{2}=p_{\theta}^{2}+\frac{1}{\sin^{2}\theta}p_{\varphi}^{2}$
with the physical meanings: the total Hamiltonian, the particle energy, the
third component of the angular momentum and the square of the total angular
momentum respectively. In the following, we also use $L_{z}$ for $p_{\varphi}$
and $E$ for $p_{\tau}$.

We consider a particle moving in the Schwarzschild spacetime along geodesics
with invariant mass $\mu$, so along a geodesic, we have
$g_{\mu\nu}\frac{dx^{\mu}}{dt}\frac{dx^{\nu}}{dt}=-\mu^{2}.$ We choose
$\mu^{2}=0$ for massless particles and $\mu^{2}=1$ for massive particles, and
the geodesics in the former case is called a lightlike geodesic and the latter
timelike.

####  3.1.1. The radial dynamics, critical points and homoclinic orbits

We next analyze the radial dynamics. Setting $p_{\tau}=E$,
$H=-\frac{1}{2}\mu^{2}$ and $\frac{dr}{dt}=\alpha p_{r}$, we obtain (3.2)

(3.3) $\displaystyle\frac{1}{2}E^{2}$
$\displaystyle=\frac{1}{2}\left(\frac{dr}{dt}\right)^{2}+V(r),$ $\displaystyle
2V(r)$
$\displaystyle=\left(\mu^{2}+\frac{L^{2}}{r^{2}}\right)\alpha=\mu^{2}-\frac{2\mu^{2}M}{r}+\frac{L^{2}}{r^{2}}-\frac{2ML^{2}}{r^{3}}.$

Thus we may visualize the radial dynamics in the $(r,\frac{dr}{dt})$-plane as
a particle moving in the potential well of $V$. The critical points of $V$
corresponds to orbits with $r=\mathrm{const}$. Setting $\frac{dV}{dr}=0$ we
obtain

$\mu^{2}Mr^{2}-L^{2}r+3ML^{2}=0$

with solutions

(3.4) $r_{\pm}=\frac{L^{2}\pm\sqrt{L^{4}-12M^{2}L^{2}}}{2M}$

for $\mu^{2}=1,$ and $r=3M$ for $\mu=0$.

The case of $\mu^{2}=0$ differs drastically from the case of $\mu^{2}=1$.
Indeed, when $\mu^{2}=0$, we have
$V(r)=\frac{L^{2}}{r^{2}}-\frac{2ML^{2}}{r^{3}}\to-\infty$ as $r\to 0$,
$V(\infty)=0$, and the point $r=3M$ is the global maximum point and the unique
critical point of $V$. This means that the bound spherical orbits are the only
bound orbits in the lightlike case.

(a) $\mu^{2}=1$, $L^{2}>L^{2}_{isco}$ but not large

(b) $\mu^{2}=1$, large $L^{2}$

Figure 2. The Schwarzschild effective potential

We consider next the case $\mu^{2}=1$. For $L$ satisfying $L^{2}>12M^{2}$,
there are two solutions $r_{+}\neq r_{-}$ by (3.4). The root $r_{-}$ is a
local maximum of $V$ and $r_{+}$ is a local minimum. All orbits with initial
value $r<r_{-}$ plunge into the blackhole. When $L^{2}=12M^{2}$, the two
critical points of $V$ coincide. The value $12M^{2}$ is called $L_{isco}^{2}$
where _isco_ stands for the _innermost stable circular orbit_. For all values
$L^{2}<L^{2}_{isco}$, all orbits plunge into the central blackhole. The set
$\\{r=r_{-}\\}$ is called the _photon sphere_ , orbits on which remain on a
circle with constant radius $r_{-}$ and are unstable under radial
perturbations.

####  3.1.2. The normally hyperbolic invariant manifold

Denote by $r_{c}=3M$ in the $\mu^{2}=0$ case and $r_{c}=r_{-}$ in the
$\mu^{2}=1$ case. If we fix the constant $p_{\tau}^{2}=E^{2}=2V(r_{c})$, which
is a function of $L^{2}$ and $M$, then we can treat the Hamiltonian $H$ in
(3.2) as a system of three degrees of freedom depending on the variables
$r,p_{r},\theta,p_{\theta},\varphi,p_{\varphi}$. The four dimensional
submanifold

$N:=\left\\{r=r_{c},\ \frac{dr}{dt}=0\right\\}$

is a normally hyperbolic invariant manifold (NHIM) in the sense of Definition
A.1. One remarkable property of the NHIM is its persistence under small
perturbations, which is summarized in Theorem A.2 in Appendix A. The normal
Lyapunov exponents are $\pm\sqrt{-V^{\prime\prime}(r_{c})}$ obtained as the
eigenvalues of the linearized radial Hamiltonian equation

$\left(\begin{array}[]{c}\delta r\\\
\delta\dot{r}\end{array}\right)^{\prime}=\left(\begin{array}[]{cc}0&1\\\
-V^{\prime\prime}(r_{c})&0\end{array}\right)\left(\begin{array}[]{c}\delta
r\\\ \delta\dot{r}\end{array}\right)$

at $r=r_{c}$. The coefficient matrix has two eigenvalues of opposite signs due
to $V^{\prime\prime}(r_{c})<0$. The Hamiltonian system restricted to $N$ is a
Hamiltonian system of two degrees of freedom of the form

(3.5)
$H_{N}(\theta,p_{\theta},\varphi,p_{\varphi})=r_{c}^{-2}\left(p_{\theta}^{2}+\frac{1}{\sin^{2}\theta}L_{z}^{2}\right)-\alpha^{-1}(r_{c})E^{2}.$

###  3.2. Dynamics on the photon shell

The system (3.5) has two degrees of freedom and is integrable in the sense of
Theorem 2.1. Note that this is the case for both $\mu^{2}=0$ and $\mu^{2}=1$.

The nontrivial part of the Hamiltonian (3.5) is the square of the total
angular momentum, denoted by

$L^{2}=p_{\theta}^{2}+\frac{1}{\sin^{2}\theta}L_{z}^{2}.$

The dynamics of this Hamiltonian is as follows. First the Hamiltonian has no
dependence on $\varphi$ so $L_{z}$ is a constant of motion. For fixed $L_{z}$,
the resulting system describes only the latitudinal motion, which can be
considered as a particle moving in the potential well of
$\frac{1}{\sin^{2}\theta}L_{z}^{2}$, which is a infinitely deep single well.
Let $\theta^{-}<\theta^{+}$ be two zeros of the
$L^{2}-\frac{1}{\sin^{2}\theta}L_{z}^{2}$, then along an orbit, the
latitudinal angle moves in the spherical strip
$\theta^{-}\leq\theta\leq\theta^{+}$.

We next introduce action-angle coordinates following Section 50 of [A89]. We
introduce the action variable

(3.6) $J_{\theta}=\frac{1}{2\pi}\oint
p_{\theta}\,d\theta=\frac{1}{\pi}\int_{\theta^{-}}^{\theta^{+}}\sqrt{L^{2}-\frac{1}{\sin^{2}\theta}L_{z}^{2}}\,d\theta.$

The other action variable is taken to be $L_{z}$ and their dual angular
variables are denoted by $\alpha_{\theta},\alpha_{\varphi}\in\mathbb{T}^{1}$
respectively, introduced to preserve the symplectic form

$dp_{\theta}\wedge d\theta+dL_{z}\wedge d\varphi=dJ_{\theta}\wedge
d\alpha_{\theta}+dL_{z}\wedge d\alpha_{\varphi}.$

In the new coordinates the Hamiltonian $H:=L^{2}(J_{\theta},L_{z}):\
\mathbb{R}_{+}\times\mathbb{T}^{1}\times\mathbb{R}\times\mathbb{T}^{1}\to\mathbb{R}$
is independent of the angular variables $\alpha_{\theta},\alpha_{\varphi}$.
Then the Hamiltonian equation becomes
$\begin{cases}\dot{\alpha}_{\theta}=\frac{\partial H}{\partial J_{\theta}}\\\
\dot{\alpha}_{\varphi}=\frac{\partial H}{\partial L_{z}}\end{cases}$. Hence
the Hamiltonian flow is a translation on $\mathbb{T}^{2}$ via
$(\alpha_{r},\alpha_{\theta})\mapsto(\alpha_{\theta},\alpha_{\varphi})+t\omega,\mathrm{mod}\
\mathbb{Z}^{2},\ \ t\in\mathbb{R}$, where
$\omega=(\omega_{\theta},\omega_{\varphi})=(\frac{\partial H}{\partial
J_{\theta}},\frac{\partial H}{\partial L_{z}})$ is called the frequency.

When we restrict to the level set $H_{N}=-\frac{1}{2}\mu^{2}$ and perform the
energetic reduction procedure as we did in Section 2.3 by taking $-J_{\theta}$
as the new Hamiltonian, we see that we get a the resulting time-1 map
$\phi_{0}:\ \mathbb{T}\times I\to\mathbb{T}\times I$ via
$(\alpha_{\varphi},L_{z})\mapsto(\alpha_{\varphi}+\nu,L_{z})$, where
$\nu:=\omega_{\varphi}/\omega_{\theta}$ is exactly the rotation number in the
integrable case, and $I=[\varepsilon,|L|]$ for any small $\varepsilon>0$. In
order to apply Moser’s theorem and the twist map theory we have to verify the
twist condition (Definition 2.4(2)), which is exactly that the ratio
$\omega_{\varphi}/\omega_{\theta}$ depends on its parameter $L_{z}$ strictly
monotonically.

The following result show that the map $\phi_{0}$ does not satisfy the twist
condition, so neither Moser’s theorem nor the twist map theory applies to a
small perturbation of $\phi_{0}$ a priori.

######  Proposition 3.1.

The Hamiltonian $H_{N}$ written in action-angle coordinates satisfies

  1. (1)

$\frac{\partial H_{N}}{\partial J_{\theta}}=\frac{\partial H_{N}}{\partial
L_{z}}$;

  2. (2)

for fixed $L_{z}>0$, as a function of $J_{\theta}$, the Hamiltonian is
strictly convex in the region $J_{\theta}>0$.

###### Proof.

We first prove item (1). Let $f$ be the map relating the two sets of constants
of motion $(J_{\theta},L_{z})$ and $(L^{2},L_{z})$, i.e.
$f(J_{\theta},L_{z})=(L^{2},L_{z})$. From $Df\cdot Df^{-1}=I$, we obtain

$\left[\begin{array}[]{ccc}\frac{\partial L^{2}}{\partial
J_{\theta}}&\frac{\partial L^{2}}{\partial L_{z}}\\\ \frac{\partial
L_{z}}{\partial J_{\theta}}&\frac{\partial L_{z}}{\partial
L_{z}}\end{array}\right]\left[\begin{array}[]{ccc}\frac{\partial
J_{\theta}}{\partial L^{2}}&\frac{\partial J_{\theta}}{\partial L_{z}}\\\
0&1\end{array}\right]=I.$

This implies $Df=\left[\begin{array}[]{ccc}\left(\frac{\partial
J_{\theta}}{\partial L^{2}}\right)^{-1}&-\frac{\partial J_{\theta}}{\partial
L_{z}}\cdot\left(\frac{\partial J_{\theta}}{\partial L^{2}}\right)^{-1}\\\
0&1\end{array}\right].$ Denote by

$\omega:=DH(J_{\theta},L_{z})=\left(\left(\frac{\partial J_{\theta}}{\partial
L^{2}}\right)^{-1},-\frac{\partial J_{\theta}}{\partial
L_{z}}\cdot\left(\frac{\partial J_{\theta}}{\partial
L^{2}}\right)^{-1}\right)$

the first row of $Df$.

We next show that $-\frac{\partial J_{\theta}}{\partial L_{z}}=1$ so that
$\omega=\frac{\partial L^{2}}{\partial J_{\theta}}(1,1)$. Denoting by
$\theta^{-}<\theta^{+}$ the two roots of $p_{\theta}$, then we have from (3.6)

$\displaystyle\frac{\partial J_{\theta}}{\partial L^{2}}$
$\displaystyle=\frac{1}{2\pi}\int_{\theta^{-}}^{\theta^{+}}\left(L^{2}-\frac{L_{z}^{2}}{\sin^{2}\theta}\right)^{-1/2}\,d\theta,$
$\displaystyle\frac{\partial J_{\theta}}{\partial L_{z}}$
$\displaystyle=-\frac{1}{\pi}\int_{\theta^{-}}^{\theta^{+}}\frac{L_{z}}{\sin^{2}\theta}\left(L^{2}-\frac{L_{z}^{2}}{\sin^{2}\theta}\right)^{-1/2}\,d\theta.$

Note that the integrand of the last expression is exactly
$\frac{d\varphi}{d\theta}$ (see (3.2)). Each Schwarzschild geodesic on the
photon sphere lies on a fixed plane passing through the origin hence is a
closed orbit, thus we get

$\int_{0}^{2\pi}d\varphi=2\int_{\theta^{-}}^{\theta^{+}}\frac{L_{z}d\theta}{\sin^{2}\theta\sqrt{L^{2}-\frac{L_{z}^{2}}{\sin^{2}\theta}}}.$

This gives $\frac{\partial J_{\theta}}{\partial L_{z}}=-1$, hence proves item
(1) by applying the implicit function theorem.

We next work on item (2). We develop Mather’s variational method [Ma1] to our
new setting. We consider the $(\theta,p_{\theta})$-plane with
$\theta\in[0,\pi],\ p_{\theta}\in\mathbb{R}$. The point $(\pi/2,0)$ is an
elliptic fixed point. If we remove this point from the plane we get a
topological cylinder on which we introduce the polar coordinates
$F:(\rho,\psi)\in\mathbb{R}_{+}\times\mathbb{T}\mapsto(p_{\theta},\theta,)$,
i.e. $p_{\theta}=\sqrt{\rho}\sin\psi,\ \theta-\pi/2=\sqrt{\rho}\cos\psi$. On
$\mathbb{R}_{+}\times\mathbb{T}$, we introduce the differential form
$\eta_{c}=f(\psi)d\psi$, where $f:\ \mathbb{T}\to\mathbb{R}_{+}$, with
cohomology class $c=\int_{\mathbb{T}}f(\psi)\,d\psi$. It becomes a
differential form on $[0,\pi]\times\mathbb{R}\setminus\\{(\pi/2,0)\\}$ after
being pulled-back by $F^{*}$, still denoted by $\eta_{c}$.

We next introduce the space $\mathfrak{M}$ of probability measures defined on
the phase space that are invariant under the Hamiltonian flow. For each
$\mu\in\mathfrak{M}$, we introduce

$A_{c}(\mu)=\int\mathbf{L}-\eta_{c}\,d\mu,$

where
$\mathbf{L}=\frac{1}{2}\dot{\theta}^{2}+\frac{1}{2}\dot{\varphi}^{2}\sin^{2}\theta$
is the formal Legendre transform of the Hamiltonian $H:=\frac{1}{2}L^{2}$. For
given $c=(c_{\theta},c_{\varphi})$ with $c_{\theta}>0$, we choose the
differential 1-form $\eta_{c}$ to be

$p_{\theta}\,d\theta+L_{z}d\varphi:=\pm\sqrt{L^{2}-\frac{L_{z}^{2}}{\sin^{2}\theta}}\,d\theta+L_{z}d\varphi$

such that
$\frac{1}{\pi}\int\sqrt{L^{2}-\frac{L_{z}^{2}}{\sin^{2}\theta}}\,d\theta=c_{\theta}$
and $c_{\varphi}=L_{z}$. Note that this fixes the constant $L^{2}$. Then we
get

$\displaystyle A_{c}(\mu)$ $\displaystyle=\int\mathbf{L}-\eta_{c}\,d\mu.$
$\displaystyle=\int\mathbf{L}-p_{\theta}\cdot\dot{\theta}\,-L_{z}\dot{\varphi}d\mu$
$\displaystyle=\int\frac{1}{2}(\dot{\theta}-p_{\theta})^{2}+\frac{1}{2}(\dot{\varphi}\sin\theta-\frac{L_{z}}{\sin\theta})^{2}-\frac{L_{z}^{2}}{2\sin^{2}\theta}-\frac{1}{2}p_{\theta}^{2}\,d\mu$
$\displaystyle=-\frac{1}{2}L^{2}+\int\frac{1}{2}(\dot{\theta}-p_{\theta})^{2}+\frac{1}{2}(\dot{\varphi}\sin\theta-\frac{L_{z}}{\sin\theta})^{2}\,d\mu.$

If we minimize $A_{c}(\mu)$ among all $\mu\in\mathfrak{M}$, we get that the
only choice is
$-\alpha(c):=\inf_{\mu\in\mathfrak{M}}A_{c}(\mu)=-\frac{1}{2}L^{2}$, and
$\dot{\theta}=p_{\theta},\ \dot{\varphi}=\frac{L_{z}}{\sin^{2}\theta}$. Note
that the latter determines uniquely an invariant torus of the Hamiltonian flow
for given $L_{z}$ and $L^{2}$, on which the minimizing measure $\mu$ supports.

Since we have the cohomology class $(c_{\theta},L_{z})$ is exactly the action
variables $(J_{\theta},L_{z})$, we get that the function $\alpha$ as a
function of the cohomology classes is exactly the total energy
$\frac{1}{2}L^{2}$ as a function of the action variables $(J_{\theta},L_{z})$.
We next prove that $\alpha(c)$ is convex. Let $\mu_{1}$ and $\mu_{2}$ be such
that $-\alpha(c_{1})=A_{c_{1}}(\mu_{1}),\ -\alpha(c_{2})=A_{c_{2}}(\mu_{2})$
and $\mu_{\lambda}$ be such that $-\alpha(\lambda
c_{1}+(1-\lambda)c_{2})=A_{\lambda c_{1}+(1-\lambda)c_{2}}(\mu_{\lambda})$.
Then we have

$\displaystyle-\alpha(\lambda c_{1}+(1-\lambda)c_{2})$
$\displaystyle=A(\mu_{\lambda})=\int\mathbf{L}-\eta_{\lambda
c_{1}+(1-\lambda)c_{2}}\,d\mu_{\lambda}$
$\displaystyle=\lambda\int\mathbf{L}-\eta_{c_{1}}\,d\mu_{\lambda}+(1-\lambda)\int\mathbf{L}-\eta_{c_{2}}\,d\mu_{\lambda}$
$\displaystyle\geq\lambda\int\mathbf{L}-\eta_{c_{1}}\,d\mu_{1}+(1-\lambda)\int\mathbf{L}-\eta_{c_{2}}\,d\mu_{2}$
$\displaystyle=-\lambda\alpha(c_{1})-(1-\lambda)\alpha(c_{2}).$

This shows that $\alpha$ is convex in $c$, hence the Hamiltonian
$\frac{1}{2}L^{2}$ is convex in $J_{\theta},L_{z}$. It remains to prove the
strict convexity in the direction transverse the energy levels. We introduce
the $\beta$-function that is the Legendre dual of $\alpha$

$\beta(h)=\sup_{c}\langle h,c\rangle-\alpha(c).$

Following the argument of [Ca], we get that $\beta$ is differentiable in the
radial direction. The Legendre transform gives $h=D\alpha(c)=(\frac{\partial
H}{\partial J_{\theta}},\frac{\partial H}{\partial L_{z}})$ attaining the sup
in the definition of $\beta$. By the first item of the Proposition, we get
that $\alpha(c)=f(J_{\theta}+L_{z})$ for some convex function $f$. Then for
$L_{z}$ fixed, the differentiability of $\beta$ in the radial direction
implies that $f$ is strictly convex.

∎

The first item reflects the fact that all orbit on the photon sphere lies on a
plane passing through the origin hence is periodic and they all have the same
period on a fixed energy level. So neither KAM nor twist map theory can be
applied at this moment. However, when a perturbation is added, it may
introduce certain twistness, which require further more careful analysis
depending on the form of the perturbation.

###  3.3. KAM and QPO for Schwarzschild

In the case of $\mu^{2}=1$, we have other bound orbits than those on the
photon shell. In particular, orbits with $r=r_{+}$ are stable circular. In a
neighborhood of this stable circular orbit, there are a lot of bound orbits.

####  3.3.1. KAM and twist maps theory for the QPO

We introduce the following bounded part of the phase space where Liouville-
Arnold theorem applies. Let $C$ be a large constant and define

(3.7)
$\mathcal{B}(C):=\left\\{V(r_{+})<\frac{E^{2}}{2}<\min\\{V(r_{-}),V(\infty)\\},\quad|L_{isco}|<|L|<C\right\\}\cap
H^{-1}(-\frac{1}{2}).$

In this set $\mathcal{B}(C)$, we note that both the radial $r$-component and
the latitudinal $\theta$-component are lying in potential wells around the
local minimums, so nearby orbit will have small oscillations in the $r$\- and
$\theta$-components (c.f. Figure 2).

Suppose for simplicity that we consider a stationary and axisymmetric
perturbation of the Schwarzchild metric so that the perturbation has no
dependence on $\varphi$ and $\tau$ then $E$ and $L_{z}$ are constants of
motion. Fixing a value for $E$ and $L_{z}$ satisfying the definition of
$\mathcal{B}(C)$, then the resulting Hamiltonian system has two degrees of
freedom with coordinates $(r,p_{r},\theta,p_{\theta})$.

We next introduce the action-angle coordinates and work out the frequencies of
the oscillations. The action variables $J_{r}$ and $J_{\theta}$ are defined as
follows (c.f. (3.2))

(3.8) $\displaystyle J_{r}$ $\displaystyle=\frac{1}{2\pi}\oint
p_{r}dr=\frac{1}{\pi}\int_{r^{-}}^{r^{+}}\frac{1}{\alpha}\sqrt{E^{2}-\alpha
r^{-2}L^{2}-\alpha}\,dr,$ $\displaystyle J_{\theta}$
$\displaystyle=\frac{1}{2\pi}\oint
p_{\theta}d\theta=\frac{1}{\pi}\int_{\theta^{-}}^{\theta^{+}}\sqrt{L^{2}-\frac{L_{z}^{2}}{\sin^{2}\theta}}\,d\theta,$

where $r^{-}<r^{+}$ are zeros of the integrand of $J_{r}$, and the definition
of $J_{\theta}$ is the same as (3.6). The dual angular coordinates
$\alpha_{r},\alpha_{\theta}\in\mathbb{T}$ are introduced to preserve the
symplectic structure, i.e.

$dp_{r}\wedge dr+dp_{\theta}\wedge d\theta=dJ_{r}\wedge
d\alpha_{r}+dJ_{\theta}\wedge d\alpha_{\theta}.$

The motion of the angular coordinates $(\alpha_{r},\alpha_{\theta})$ on
$\mathbb{T}^{2}$ is given by a linear motion with frequency
$\omega=(\omega_{r},\omega_{\theta})=(\frac{\partial H}{\partial
J_{r}},\frac{\partial H}{\partial J_{\theta}})$ where $H$ is the Schwarzschild
Hamiltonian written in action-angle coordinates and with fixed parameters $E$
and $L_{z}$. We call $\omega_{r}:=\frac{\partial H}{\partial J_{r}}$ the
radial frequency and $\omega_{\theta}:=\frac{\partial H}{\partial J_{\theta}}$
the latitudinal frequency.

The ratio $\omega_{r}/\omega_{\theta}$ can be parametrized by one parameter
that we choose to be $L^{2}$. In order to apply KAM theorem and the twist map
theory, we need to find an interval where the ratio has no critical point. We
see from Figure 3 (plotted with $E=0.97$ and $L_{z}=3$) that this is indeed
the case.

######  Theorem 3.2.

Let $\varepsilon h_{\mu\nu}dx^{\mu}dx^{\nu}$ be a stationary and axisymmetric
perturbation to the Schwarzschild metric. Then there exist an open set of
parameters $(E,L_{z})$, a neighborhood $\mathcal{D}$ of the point
$(r,p_{r},\theta,p_{\theta})=(r_{+},0,\pi/2,0)$ and $\varepsilon_{0}$ such
that KAM and the twist map theory are applicable to the perturbed system in
$\mathcal{D}$ and for all $0<|\varepsilon|<\varepsilon_{0}$.

Figure 3. The ratio $\omega_{r}/\omega_{\theta}$ for Schwarzschild with
parameters $E=0.97$ and $L_{z}=3$

###### Proof.

Let $f$ be the map relating the two sets of constants of motion
$(J_{r},J_{\theta})$ and $(H,L^{2})$: $f(J_{r},J_{\theta})=(H,L^{2})$, where
$H$ is the Schwarzchild Hamiltonian, then we have $Df\cdot Df^{-1}=I$, which
is more explicitly

$\left[\begin{array}[]{ccc}\frac{\partial H}{\partial J_{r}}&\frac{\partial
H}{\partial J_{\theta}}\\\ \frac{\partial L^{2}}{\partial
J_{r}}&\frac{\partial L^{2}}{\partial
J_{\theta}}\end{array}\right]\left[\begin{array}[]{ccc}\frac{\partial
J_{r}}{\partial H}&\frac{\partial J_{r}}{\partial L^{2}}\\\ 0&\frac{\partial
J_{\theta}}{\partial L^{2}}\end{array}\right]=I.$

Here we use the fact that $J_{\theta}$ is independent of $H$ as we see from
the expression of $J_{\theta}$. This implies
$Df=\left[\begin{array}[]{ccc}\left(\frac{\partial J_{r}}{\partial
H}\right)^{-1}&-\frac{\partial J_{r}}{\partial L^{2}}\cdot\left(\frac{\partial
J_{r}}{\partial H}\frac{\partial J_{\theta}}{\partial L^{2}}\right)^{-1}\\\
0&\left(\frac{\partial J_{\theta}}{\partial
L^{2}}\right)^{-1}\end{array}\right].$ Denote by

$\omega=(\omega_{r},\omega_{\theta})=\left(\left(\frac{\partial
J_{r}}{\partial H}\right)^{-1},-\frac{\partial J_{r}}{\partial
L^{2}}\cdot\left(\frac{\partial J_{r}}{\partial H}\frac{\partial
J_{\theta}}{\partial L^{2}}\right)^{-1}\right)$

the first row of $Df$. Denoting $P(r)=E^{2}-\alpha r^{-2}L^{2}-\alpha$, we get
from (3.8)

$\displaystyle\frac{\partial J_{r}}{\partial L^{2}}$
$\displaystyle=-\frac{1}{2\pi}\int_{r^{-}}^{r^{+}}\frac{1}{r^{2}\sqrt{P(r)}}\,dr,$
$\displaystyle\frac{\partial J_{\theta}}{\partial L^{2}}$
$\displaystyle=\frac{1}{2\pi}\int_{\theta^{-}}^{\theta^{+}}\frac{1}{\sqrt{L^{2}-\frac{L_{z}^{2}}{\sin^{2}\theta}}}\,d\theta.$

So we get

$\frac{\omega_{r}}{\omega_{\theta}}=-\frac{\partial J_{\theta}}{\partial
L^{2}}\cdot\left(\frac{\partial J_{r}}{\partial
L^{2}}\right)^{-1}=\int_{\theta^{-}}^{\theta^{+}}\frac{1}{\sqrt{L^{2}-\frac{L_{z}^{2}}{\sin^{2}\theta}}}\,d\theta\Big{/}\int_{r^{-}}^{r^{+}}\frac{1}{r^{2}\sqrt{P(r)}}\,dr,$

which is a function of $L^{2}$ from the above formulas. Figure 3 is plotted
using this formula. Since we have $\frac{\partial L^{2}}{\partial
J_{r}}=0,\frac{\partial L^{2}}{\partial J_{\theta}}=(\frac{\partial
J_{\theta}}{\partial L^{2}})^{-1}$, we can then eventual compute
$\frac{\partial}{\partial J_{\theta}}\frac{\omega_{r}}{\omega_{\theta}}.$ We
see from Figure 3 that the statement holds for a particularly chosen paramter
$(E,L_{z})$. Since the action-angle coordinate change is smooth, so we get the
statement for the parameter $(E,L_{z})$ in an open set. The proof of the
statement is now complete.

∎

####  3.3.2. KAM nondegeneracy and further derivatives

The preceding Theorem 3.2 uses Moser’s version of KAM theorem. If we would
like to verify the KAM nondegeneracy conditions in Section 2.2, we have to
calculate $\nabla\omega$ where $\nabla$ is the derivative with respect to
$J_{r},J_{\theta}$. In this section, for the sake of completeness, we show how
to take further derivative of $\omega$ and clarify some subtleties. Since we
do not use it in the proof of the main theorem, readers are suggested to skip
this section when first reading.

We have

(3.9) $\displaystyle\nabla\omega$
$\displaystyle=\frac{\partial(E^{2},L^{2})}{\partial(J_{r},J_{\theta})}(\partial_{E^{2}},\partial_{L^{2}})\omega$
$\displaystyle=Df\cdot\left[\begin{array}[]{ccc}\partial_{E^{2}}\left(\frac{\partial
J_{r}}{\partial E^{2}}\right)^{-1}&\partial_{L^{2}}\left(\frac{\partial
J_{r}}{\partial E^{2}}\right)^{-1}\\\ -\partial_{E^{2}}\left(\frac{\partial
J_{r}}{\partial L^{2}}\cdot\left(\frac{\partial J_{r}}{\partial
E^{2}}\frac{\partial J_{\theta}}{\partial
L^{2}}\right)^{-1}\right)&-\partial_{L^{2}}\left(\frac{\partial
J_{r}}{\partial L^{2}}\cdot\left(\frac{\partial J_{r}}{\partial
E^{2}}\frac{\partial J_{\theta}}{\partial
L^{2}}\right)^{-1}\right)\end{array}\right]$
$\displaystyle=Df\cdot\left(\frac{\partial J_{r}}{\partial
E^{2}}\right)^{-3}\left(\frac{\partial J_{\theta}}{\partial
L^{2}}\right)^{-1}\left[\begin{array}[]{ll}\frac{\partial^{2}J_{r}}{(\partial
E^{2})^{2}}&\frac{\partial^{2}J_{r}}{\partial L^{2}\partial E^{2}}\\\
\frac{\partial^{2}J_{r}}{\partial E^{2}\partial
L^{2}}&\frac{\partial^{2}J_{r}}{(\partial
L^{2})^{2}}-\frac{\partial^{2}J_{\theta}}{(\partial L^{2})^{2}}\frac{\partial
J_{r}}{\partial L^{2}}\cdot\left(\frac{\partial J_{\theta}}{\partial
L^{2}}\right)^{-1}\end{array}\right].$

The expressions in the last matrix are all computable from the definitions of
$J_{r},J_{\theta}$.

To obtain the second order derivatives, we cannot differentiate the first
order derivatives in (3.3.1) for the technical reason that further derivative
will introduce divergent improper integrals. This problem can be resolved as
follows. We consider only the case $J_{r}$ and the case $J_{\theta}$ is
similar. Note that the singularity originates from the fact that
$P(r^{\pm})=0$ will appear in the denominator, which we should avoid. The
action $J_{r}$ is the area enclosed by the curve Graph$p_{r}$. Let $r_{*}$ be
the minimum of $\alpha r^{-2}L^{2}+2\alpha$. Then we split the integral
$J_{r}=\int_{r^{-}}^{r^{+}}\frac{1}{\pi\alpha}\sqrt{P(r)}\,dr$ into two pieces
$J_{r}^{1}=\int_{r^{-}}^{r_{*}}\frac{1}{\pi\alpha}\sqrt{P(r)}\,dr$ and
$J_{r}^{2}=\int_{r_{*}}^{r^{+}}\frac{1}{\pi\alpha}\sqrt{P(r)}\,dr$ and
consider only $J_{r}^{1}$ since the other one is completely similar. We next
pick $\tilde{r}\in(r^{-},r_{*})$ sufficiently close to $r^{-}$ and denote by
$\tilde{p}$ the positive value of $p_{r}(\tilde{r})$. We denote
$J_{r}^{1,R}=\int_{\tilde{r}}^{r_{*}}\frac{1}{\pi\alpha}\sqrt{P(r)}\,dr$ the
area enclosed by $r=\tilde{r},r=r_{*}$ and Graph$p_{r}$. For the area to the
left of the line $r=\tilde{r}$, we realize the Graph$p_{r}$ to the left of
$r=\tilde{r}$ as a graph over the $p_{r}$-axis. So we need to solve the
function $\alpha^{2}p_{r}^{2}=P(r(p_{r}))$ for $r$ as a function of $p_{r}$.
This function is solvable if $\tilde{r}$ is sufficiently close to $r_{-}$
since $P^{\prime}(r^{-})\neq 0$ hence the implicit function theorem applies.
With this function $r(p_{r})$, we introduce the integral
$J^{1,L}_{r}=\frac{1}{\pi}\int_{-\tilde{p}}^{\tilde{p}}r(p_{r})dp_{r}$. The
overlap of $J^{1,L}_{r}$ and $J^{1,R}_{r}$ is a square of area
$2\tilde{p}|r_{*}-\tilde{r}|.$ Therefore we obtain

$J^{1}_{r}=J^{1,R}_{r}+J^{1,L}_{r}-2\tilde{p}|r_{*}-\tilde{r}|=\int_{\tilde{r}}^{r_{*}}\frac{1}{\pi\alpha}\sqrt{P(r)}\,dr+\frac{1}{\pi}\int_{-\tilde{p}}^{\tilde{p}}r(p_{r})dp_{r}-2\tilde{p}|r_{*}-\tilde{r}|.$

We remark that $\tilde{p}$ and $r(p_{r})$ are both dependent on the variables
$E^{2}$ and $L^{2}$ while $\tilde{r}$ is not. The derivatives of this
expression does not involve improper integrals hence we can take its second
derivative to yield an expression for $D\omega$.

####  3.3.3. Three fundamental frequencies

In Theorem 3.2 we consider stationary and axisymmetric perturbations so we
have only two fundamental frequencies $\omega_{r}$ and $\omega_{\theta}$ to
consider. In general, when a perturbation depends on $\varphi$, i.e.
nonaxisymmetric, but not on $\tau$, we have to treat the Hamiltonian system
with as one three degrees of freedom, in which case, we cannot apply Moser’s
theorem or twist map theory, so we have to verify the KAM nondegeneracy
assumptions in Section 2.2 if we wish to obtain KAM type dynamics. Again
readers are suggested to skip this section when first reading, since here we
only show how to get the formulas for $\nabla\omega$, where $\nabla$ means the
derivative with respect to the action variables
$J_{r},J_{\theta},J_{\varphi}$. We introduce the action variables
$J_{r},J_{\theta}$ as before and the extra $J_{\varphi}$ is chosen to be
$L_{z}$.

Then there is a map $F$ relating the two set of conserved quantities:
$F(J_{r},J_{\theta},J_{\varphi})=(H,L^{2},L_{z})$. We have $DF\cdot
DF^{-1}=I$. More explicitly we have

$\left[\begin{array}[]{ccc}\frac{\partial H}{\partial J_{r}}&\frac{\partial
H}{\partial J_{\theta}}&\frac{\partial H}{\partial J_{\varphi}}\\\
\frac{\partial L^{2}}{\partial J_{r}}&\frac{\partial L^{2}}{\partial
J_{\theta}}&\frac{\partial L^{2}}{\partial J_{\varphi}}\\\
0&0&1\end{array}\right]\left[\begin{array}[]{ccc}\frac{\partial
J_{r}}{\partial H}&\frac{\partial J_{r}}{\partial L^{2}}&0\\\ 0&\frac{\partial
J_{\theta}}{\partial L^{2}}&\frac{\partial J_{\theta}}{\partial L_{z}}\\\
0&0&1\end{array}\right]=I.$

This implies $DF=\left[\begin{array}[]{ccc}a^{-1}&-b/(ac)&bd/(ac)\\\
0&c^{-1}&-d/c\\\ 0&0&1\end{array}\right]$, where $a=\frac{\partial
J_{r}}{\partial H},\ b=\frac{\partial J_{r}}{\partial L^{2}},\
c=\frac{\partial J_{\theta}}{\partial L^{2}},\ d=\frac{\partial
J_{\theta}}{\partial L_{z}}$ can be evaluated directly from equation 3.8.

Denoting by $\omega=\omega(J_{r},J_{\theta},J_{\varphi})$ the first row of
$DF$. We thus have

$\nabla\omega=\frac{\partial(H,L^{2},L_{z})}{\partial(J_{r},J_{\theta},J_{\varphi})}(\partial_{H},\partial_{L^{2}},\partial_{L_{z}})\omega=DF\cdot\left[\begin{array}[]{ccc}\partial_{H}(a^{-1})&\partial_{L^{2}}(a^{-1})&0\\\
\partial_{H}(-b/(ac))&\partial_{L^{2}}(-b/(ac))&\partial_{L_{z}}(-b/(ac))\\\
\partial_{H}(bd/(ac))&\partial_{L^{2}}(bd/(ac))&\partial_{L_{z}}(bd/(ac))\\\
\end{array}\right]$

Since $DF$ is nondegenerate, the nondegeneracy condition is then reduced to
that of $(\partial_{E},\partial_{L^{2}},\partial_{L_{z}})\omega$.

##  4\. The Kerr spacetime

In this section, we study the dynamics of perturbed Kerr spacetime. The Kerr
spacetime in the standard Boyer-Lindquist coordinates has the form

$ds^{2}=-\left(1-\frac{2Mr}{\Sigma}\right)d\tau^{2}-2\left(\frac{2Mr}{\Sigma}\right)a\sin^{2}\theta
d\tau
d\varphi+\Sigma\left(\frac{dr^{2}}{\Delta}+d\theta^{2}\right)+\frac{\mathcal{A}}{\Sigma}\sin^{2}\theta
d\varphi^{2},$

where

(4.1) $\displaystyle\Sigma$ $\displaystyle=r^{2}+a^{2}\cos^{2}\theta,$
$\displaystyle\Delta$ $\displaystyle=r^{2}+a^{2}-2Mr,$
$\displaystyle\mathcal{A}$ $\displaystyle=(r^{2}+a^{2})^{2}-\Delta
a^{2}\sin^{2}\theta.$

Here $M$ is the mass and $a$ is the angular momentum of the blackhole with
$0\leq|a|\leq M$. When $a=0$, the Kerr metric reduces to the Schwarzschild
metric. The event horizon is where the metric coefficient
$g_{rr}=\frac{\Sigma}{\Delta}$ becomes singular. Let
$r_{h}=M+\sqrt{M^{2}-a^{2}}$ be the larger root of $\Delta$, hence the outer
event horizon is the sphere $r=r_{h}$.

The geodesic equations are of the following form (assuming $ds^{2}=-\mu^{2}$).

(4.2) $\begin{cases}\Sigma^{2}(\frac{dr}{dt})^{2}&=R,\\\
\Sigma^{2}(\frac{d\theta}{dt})^{2}&=\Theta,\\\
\Sigma\frac{d\tau}{dt}&=-a(aE\sin^{2}\theta-
L_{z})+\frac{r^{2}+a^{2}}{\Delta}P(r),\\\
\Sigma\frac{d\varphi}{dt}&=\frac{a}{\Delta}P(r)-\left(aE-\frac{L_{z}}{\sin^{2}\theta}\right),\end{cases}$

where we have

(4.3) $\displaystyle R(r)=P(r)^{2}-\Delta(\mu^{2}r^{2}+(L_{z}-aE)^{2}+Q),$
$\displaystyle\Theta(\theta)=Q-\left[(\mu^{2}-E^{2})a^{2}+\dfrac{L_{z}^{2}}{\sin^{2}\theta}\right]\cos^{2}\theta,$
$\displaystyle P(r)=E(r^{2}+a^{2})-aL_{z},$

and the constant
$Q=p_{\theta}^{2}+\cos^{2}\theta\left(a^{2}(\mu^{2}-E^{2})+\left(\frac{L_{z}}{\sin\theta}\right)^{2}\right)$
is a constant of motion called _Carter constant_.

We next introduce the Hamiltonian formalism. We treat the Kerr metric as twice
a Lagrangian and get the corresponding Hamiltonian via the formal Legendre
transform
$H(p_{x^{\mu}},x^{\mu})=p_{x^{\mu}}\frac{dx^{\mu}}{dt}-L(x^{\mu},\frac{dx^{\mu}}{dt})$
where $x^{\mu}=(\tau,r,\theta,\varphi)$ and
$p_{x^{\mu}}=(p_{\tau},p_{r},p_{\theta},p_{\varphi})$ with

(4.4)
$\begin{cases}p_{\tau}&=E=-\left(1-\frac{2Mr}{\Sigma}\right)\dot{\tau}-\left(\frac{2Mr}{\Sigma}\right)a\sin^{2}\theta\dot{\varphi},\\\
p_{r}&=\Sigma\frac{\dot{r}}{\Delta},\\\ p_{\theta}&=\Sigma\dot{\theta},\\\
p_{\varphi}&=L_{z}=-\left(\frac{2Mr}{\Sigma}\right)a\sin^{2}\theta\dot{\tau}+\frac{\mathcal{A}}{\Sigma}\sin^{2}\theta\dot{\varphi}.\end{cases}$

The Hamiltonian $H$ can be written as

(4.5) $\displaystyle H=\frac{1}{2}\Sigma^{-1}(H_{r}+H_{\theta}),$
$\displaystyle H_{r}=\Delta
p_{r}^{2}-\dfrac{1}{\Delta}((r^{2}+a^{2})E-aL_{z})^{2},$ $\displaystyle
H_{\theta}=p_{\theta}^{2}+\dfrac{1}{\sin^{2}\theta}(a\sin^{2}\theta
E-L_{z})^{2}.$

The system has four independent constants of motion $H,E,L_{z},Q$.

###  4.1. The vertical dynamics

Introducing a time reparametrization $d\lambda=\Sigma(r,\theta)^{-1}dt$ called
Mino time, we determine the vertical dynamics by the equation

(4.6) $(\frac{d\theta}{d\lambda})^{2}=\Theta=Q-U(\theta),\ \mathrm{where}\
U(\theta)=\left[(\mu^{2}-E^{2})a^{2}+\dfrac{L_{z}^{2}}{\sin^{2}\theta}\right]\cos^{2}\theta.$

So the dynamics can be visualized as a particle moving in the potential well
of $U(\theta)$.

The case of $L_{z}=0$ differs drastically from $L_{z}\neq 0$ case. Indeed, if
$L_{z}=0$ and $Q>(\mu^{2}-E)a^{2}$, the orbit goes through the south and north
poles, while $L_{z}\neq 0$, as will will see in the following, the orbit is
bounded away from the the south and north poles.

Figure 4. Graph of $U$

We first consider the case $\mu^{2}=1$ and suppose $(\mu^{2}-E^{2})>0$, then
we have $U(\theta)\geq 0$ is an infinitely deep single well potential and
$Q\geq 0$ (see Figure 4). The minimum $\min U=0$ is attained at
$\theta=\pi/2$. We also have $U(\pi/2+\theta)=U(\pi/2-\theta)$ for
$\theta\in[0,\pi/2)$ and $U(\theta)\to\infty$ for $\theta\to 0,\pi.$ Let
$\theta^{-}<\theta^{+}$ be two root of the equation $\Theta(\theta)=0$ with
$Q>0$. Then the orbit on the photon sphere with conserved quantities
$Q,E,L_{z}$ oscillates within the spherical strip
$\theta\in[\theta^{-},\theta^{+}]$. The explicit solution of $\theta(\lambda)$
is obtained as follows. Denoting by

$\ell=L_{z}/E,\ q=Q/E^{2},\ v=\cos\theta,\ \alpha^{2}=a^{2}(E^{-2}-1)>0$

and

$\Theta_{v}=q-(\ell^{2}+q+\alpha^{2})v^{2}+\alpha^{2}v^{4},$

and by $v^{2}_{\pm}$ the two roots of the equation $\Theta_{v}=0$ with
$v_{-}^{2}<v_{+}^{2}$. Then we get

$\lambda-\nu=\int_{0}^{\theta}\frac{d\theta}{\sqrt{\Theta}}=E^{-1}\int_{0}^{v}\frac{dv}{\sqrt{\Theta}_{v}}=\frac{1}{E\alpha}\int_{0}^{v}\frac{dv}{\sqrt{(v^{2}-v_{+}^{2})(v^{2}-v_{-}^{2})}}=\frac{1}{E\alpha
v_{+}}F(\psi,k),$

where $\nu$ is the integral constant, $k=\frac{v_{-}}{v_{+}}$ and
$\sin\psi=\frac{v}{v_{-}}$, and $F$ is the incomplete elliptic integral of the
first kind. This gives us the function of $\lambda$ in terms of $v$ hence of
$\theta$, then inverting it we obtain the function $\theta(\lambda)$.

We next consider the case of $\mu^{2}=0$.

######  Lemma 4.1.

  1. (1)

When $a^{2}E^{2}\leq L_{z}^{2}$, the potential $U$ is a single well with a
global minimum value at $\theta=\pi/2$;

  2. (2)

when $a^{2}E^{2}>L_{z}^{2}$, the potential $U$ is a double well with a local
max value 0 at $\theta=\pi/2$.

###### Proof.

When $a^{2}E^{2}<L_{z}^{2}$ we have $U\geq 0$ with a global mininum value 0 at
$\theta=\pi/2$. So the situation is similar to the above $\mu^{2}=1$ and
$E^{2}<1$ case. It is necessary that $Q\geq 0$ and when $Q>0$ all orbits
crosses the equator plane. On the other hand, we have that

$U^{\prime\prime}(\pi/2)=4(L_{z}^{2}-E^{2}a^{2}).$

This shows that when $L_{z}^{2}<E^{2}a^{2}$, the point $\theta=\pi/2$ is a
local max of $U$. Since we have $U(\pi/2)=0$. This implies that the potential
$U$ is negative in some region between $0$ and $\pi/2$, and between $\pi/2$
and $\pi$. Hence $U$ is a double well potential. In this case, it is possible
to have $Q<0$, in which case each orbit is confined within one well bounded
away from the equator plane. When $Q>0$, the particle can move from one well
to the other crossing the equator plane. The $Q=0$ case is critical, in which
case the point $\theta=\pi/2$ corresponds to a hyperbolic fixed point and the
orbit moving in one well approaches the fixed point in both future and pass
hence is on a homoclinic orbit. ∎

In this critical case $Q=0$ and the double well case, using the general
principle of [X], in a generically periodically perturbed Kerr metric,
separatix splitting will be created and chaotic motion will occur, which
implies that there exists orbit tunneling from one well to the other. In each
of the above cases, solutions $\theta(\lambda)$ of the equation of motion can
be found explicitly (c.f. Section 63 of [C]).

###  4.2. The photon shell and the radial dynamics

Similar to the photon sphere in Schwarzschild case, the Kerr spacetime also
admits unstable spherical orbits. However, the dynamics of these orbits in the
Kerr spacetime differs drastrically from the Schwarzschild case. Such an orbit
in the Schwarzschild case lies on a fixed plane through the origin thus is
periodic and all these orbits form a sphere. However, in the Kerr case, such
an orbit is no longer confined to a plane but may oscillate in a spherical
strip and may be quasiperiodic. Moreover, depending on the parameters
$E,L_{z},Q$, the radius varies, so the union of all such orbits forms a ring
called _photon shell_. The photon shell for Kerr was first studied by [Wi] in
the extreme case $a=M$ for timelike geodesics. The lightlike case was studied
by [T]. The equatorial case of $Q=0$ was studied by [GLP].

####  4.2.1. Determining the photon shell

An orbit in the photon shell have constant radial component. So from the
radial equation of (4.2), the following equations should be satisfied.

(4.7) $R(r)=0,\quad\frac{d}{dr}R(r)=0.$

Let $r_{c}$ be a solution. Moreover, orbits in the photon shell are unstable,
which implies the second order derivative $\frac{d^{2}}{dr^{2}}R(r_{c})>0$.

Introducing $\ell=\frac{L_{z}}{E}$ and $q=\frac{Q}{E^{2}}$, we write
explicitly equations (4.7) as

$\displaystyle r^{4}+(a^{2}-\ell^{2}-q)r^{2}+2M(q+(\ell-a)^{2})r-a^{2}q=0,$
$\displaystyle 4r^{3}+2(a^{2}-\ell^{2}-q)r+2M(q+(\ell-a)^{2})=0$

in the $\mu^{2}=0$ case and

$\displaystyle
3r^{4}+r^{2}a^{2}-q(r^{2}-a^{2})-\frac{r^{2}}{E^{2}}(3r^{2}-4Mr+a^{2})=r^{2}\ell^{2},$
$\displaystyle
r^{4}-a^{2}Mr+q(a^{2}-Mr)-\frac{r^{3}}{E^{2}}(r-M)=rM(\ell^{2}-2a\ell)$

in the $\mu^{2}=1$ case.

Each pair of equations involves four variables $r,E,\ell,q$. We solve for
$\ell$ and $q$

(4.8) $\displaystyle\ell$
$\displaystyle=\frac{1}{a(r-M)}(M(r^{2}-a^{2})-r\Delta),$ $\displaystyle q$
$\displaystyle=-\frac{r^{3}}{a^{2}(r-M)^{2}}(4M\Delta-r(r-M)^{2})$

in the $\mu^{2}=0$ case and

(4.9) $\displaystyle\ell$
$\displaystyle=\frac{1}{a(r-M)}\left(M(r^{2}-a^{2})-r\Delta\left(1-\frac{1}{E^{2}}(1-\frac{M}{r})\right)^{1/2}\right),$
$\displaystyle q$
$\displaystyle=\frac{1}{a^{2}(r-M)}\left[\frac{r^{3}}{r-M}(4a^{2}M-r(r-3M)^{2})+\frac{r^{2}}{E^{2}}\left[r(r-2M)^{2}-a^{2}M\right]\right.$
$\displaystyle\left.-\frac{2r^{3}M}{r-M}\Delta\left(1-\left(1-\frac{1}{E^{2}}(1-\frac{M}{r})\right)^{1/2}\right)\right]$

in the $\mu^{2}=1$ case. The formulas can be found in Section 63 and 64 of [C]
and we have discarded an unphysical solution in each case.

The functions $\ell$ and $q$ in (4.9) are plotted in Figure 5 with data
$M=1,a=0.8,E=1$ and Figure 5 with $M=1,a=0.8,E=0.92$ in the case of
$\mu^{2}=1$. For nearby choices of parameters the picture is qualitatively
similar since the functions are continuous in $E^{2}$ and $a$. In the massless
case, the RHS of (4.8) has no dependence on $E$. The picture of $\ell$ and $q$
in (4.8) is qualitatively similar to Figure 5 for $a=0.8$.

(a) $M=1,a=0.8,E=1$

(b) $M=1,a=0.8,E=0.92$

Figure 5. Graphs of $\ell$ and $q$ in the $\mu^{2}=1$ case

Given $r=r_{c}$, it determines the values of $\ell=\ell_{c},q=q_{c}$ such that
$r_{c}$ solves equation (4.7) with parameters $\ell_{c}$ and $q_{c}$. So for
this choice of parameters $\ell_{c},q_{c}$, the radius $r_{c}$ is a double
root of $R(r)=0$, thus $R$ can be put in the form

(4.10) $R(r)/E^{2}=(r-r_{c})^{2}\left(Ar^{2}+2rB+C\right),$

where we define
$A=1-1/E^{2}<0,B=r_{c}\left(1-\frac{1}{E^{2}}\left(1-\frac{M}{r_{c}}\right)\right),C=-\frac{a^{2}}{r_{c}^{2}}q_{c}$
in the $\mu^{2}=1$ case, and $A=1,B=r_{c},C=-a^{2}q_{c}/r_{c}^{2}$ in the
$\mu^{2}=0$ case.

We next introduce the following admissible set of radii of bound spherical
orbits.

######  Definition 4.2.

  1. (1)

In the case of $\mu^{2}=1$, let $M>0,\ |a|\in(0,M]$ and
$E\in(1-Mr_{h}^{-1},1)$ be given parameters. We define the admissible set
$\mathcal{R}:=\mathcal{R}(M,a,E)$ as the radius $r_{c}$ satisfying the
following:

    1. (a)

$r_{h}<r_{c}\leq\frac{M}{1-E^{2}}$ so that $r_{c}$ lies outside the event
horizon and the square root in (4.9) makes sense;

    2. (b)

$q(r_{c})\geq 0;$

    3. (c)

inequality $Ar_{c}^{2}+2Br_{c}+C>0$ holds, so that $r_{c}$ is a local max for
$-R$.

  2. (2)

In the case of $\mu^{2}=0$, the admissible set $\mathcal{R}:=\mathcal{R}(M,a)$
radii of bound spherical orbits is defined as the set of $r_{c}$ satisfying
$q(r_{c})\geq 0$ in (4.8).

####  4.2.2. The Hamiltonian restricted to the photon shell

When restricted to the photon shell

(4.11) $N:=\\{r=r_{c},p_{r}=0,\ r_{c}\in\mathcal{R}\\}$

and fixing an admissible energy $E$, the Hamiltonian system becomes the
Hamiltonian subsystem

(4.12) $\displaystyle
H_{N}=\Sigma(r_{c},\theta)^{-1}\left(-\dfrac{((r_{c}^{2}+a^{2})E-aL_{z})^{2}}{\Delta(r_{c})}+p_{\theta}^{2}+\dfrac{1}{\sin^{2}\theta}(a\sin^{2}\theta
E-L_{z})^{2}\right),$

obtained from the Hamiltonian (4.5) by setting $r=r_{c},p_{r}=0$. We should
further eliminate the $r_{c}$-dependence by substituting $r_{c}$ as a function
of $L_{z}$ by inverting the first equation in (4.8) or (4.9) on each of its
monotone intervals. Thus we get a Hamiltonian system of two degrees of freedom
with coordinates $(\theta,p_{\theta},\varphi,L_{z})$ integrable in the
Liouville-Arnold sense.

####  4.2.3. The radial dynamics

We look at the radial motion in (4.2). After a time rescaling
$d\lambda:=\frac{1}{\Sigma}d\tau$ called the Mino time, the radial equation of
motion becomes $\left(\frac{dr}{d\lambda}\right)^{2}-R(r)=0$. Thus the
dynamics can be visualized as a particle moving in the potential well of
$-R(r)$ on the zero energy level.

(a) $\mu^{2}=0$

(b) $\mu^{2}=1$

Figure 6. Graph of $-R$

######  Lemma 4.3.

In the case of $\mu^{2}=0$, the only bound orbits outside the event horizon
are on the photon shell, i.e. we have $r=r_{c}$ along each such orbit and
$-R^{\prime\prime}(r_{c})<0$.

###### Proof.

Denote by $\bar{r}<\bar{\bar{r}}$ the two roots of the factor $Ar^{2}+2Br+C$
in (4.10). Note that $AC<0$, so we have $\bar{r}<0$. The lemma follows from
the fact that $\bar{\bar{r}}<r_{c}$ (see Figure 6). To prove it, it is enough
to show that $Ar_{c}^{2}+2Br_{c}+C>0$, which is equivalent to
$3r_{c}^{2}>\frac{a^{2}q}{r_{c}^{2}}$. Using (4.8), the inequality is further
reduced to $r_{c}(r_{c}-M)^{2}>M(r_{c}^{2}-2Mr_{c}+a^{2})$ in the case of
$r_{c}>M$, which is clearly true. If $r_{c}=M$, the inequality is still true
if $|a|<M$. In the extreme case $|a|=M$, the event horizon is $r=M$ so the
inequality is only violated when $r_{c}\leq M$ that means that the photon
shell is not outside the event horizon. It can be checked that when $|a|=M$
and $r_{c}\leq M$, there is no bound orbit lying outside the event horizon. ∎

We next consider the $\mu^{2}=1$ case. Consider an admissible
$r_{c}\in\mathcal{R}$, then we have

$R(r)=(E^{2}-1)(r-r_{1})(r-r_{c})^{2}(r-r_{3}),\ r_{1}<r_{c}<r_{3}.$

See Figure 6. Besides the bound spherical orbit, there is also a homoclinic
orbit with $r(\lambda)$ approaching $r_{c}$ as $\lambda\to\pm\infty$. We do
not treat the dynamics of the homoclinic orbit in this paper but refer readers
to [X] for more details.

###  4.3. KAM and the photon shell dynamics

We next consider the dynamics on the photon shell in Kerr spacetime.
Restricted to the NHIM $N$ (c.f. (4.11)), the Hamiltonian $H_{N}$ is (4.12)
with two degrees of freedom. Each orbit has constant radius hence lies on a
sphere. The azimuthal angle $\varphi$ rotates periodically in $[0,2\pi)$ and
the latitudinal angle $\theta$ oscillates within a spherical band symmetric
around the equator according to (4.6).

####  4.3.1. KAM and twist map nondegeneracy

We next introduce action-angle coordinates for the Hamiltonian $H_{N}$. In the
equation $\Sigma^{2}\dot{\theta}^{2}=\Theta$ (c.f. (4.2)), we substitute
$p_{\theta}=\Sigma\dot{\theta}$ (c.f. (4.4)) to yield $p_{\theta}^{2}=\Theta$.
We consider both cases $\mu^{2}=0$ and $\mu^{2}=1$ together, and always assume
$Q>0$. When $\mu^{2}=0$, the Carter constant $Q$ can be negative, which case
can be considered similar to the procedure that we are doing now.

We next introduce the action variables $J_{\varphi}:=p_{\varphi}=L_{z}$ and

(4.13) $J_{\theta}=\frac{1}{2\pi}\oint
p_{\theta}\,d\theta=\frac{1}{\pi}\int_{\theta^{-}}^{\theta^{+}}\sqrt{\Theta}\,d\theta=\frac{1}{\pi}\int_{\theta^{-}}^{\theta^{+}}\sqrt{Q-U(\theta)}\,d\theta,$

where $\theta^{-}<\theta^{+}$ are two roots of the integrand (c.f. (4.6)). We
further introduce their dual angular variables
$\alpha_{\theta},\alpha_{\varphi}\in\mathbb{T}$ respectively using the
symplectic form (c.f. Section 50 of [A89]). In the new coordinates, we get
that the Hamiltonian is a function of $J_{\theta}$ and $L_{z}$ only, and the
the dynamics of $(\alpha_{\theta},\alpha_{\varphi})$ on $\mathbb{T}^{2}$ is a
linear flow with frequency
$\omega=(\omega_{\theta},\omega_{\varphi})=(\frac{\partial H_{N}}{\partial
J_{\theta}},\frac{\partial H_{N}}{\partial L_{z}})$. Again the KAM
nondegeneracy and the twist condition amount to show that the ratio
$\nu:=\omega_{\varphi}/\omega_{\theta}$ has no critical point.

Here we point out some essential differences between the massive case
$(\mu^{2}=1)$ and massless case $(\mu^{2}=0)$. In the case of $\mu^{2}=0$, we
note that the $E$ variable can be eliminated by the rescaling $L_{z}\mapsto
L_{z}/E,\ Q\mapsto Q/E^{2},\ p_{\theta}\mapsto p_{\theta}/E$. The RHS of
equation (4.8) has no dependence on the particle’s $E$, nor does the
Hamiltonian (4.12). Thus the dynamics, in particular the frequency ratio
$\nu$, on the photon shell has no dependence on $E$ so in the following it is
enough to choose $E=1$. However, for the massive case, the energy $E$ does
play a role in the frequency ratio $\nu$.

In the proof of the next theorem, we provide explicit formula for the ratio
$\nu$. Let $a,M,E$ be fixed constants. Restricted to the energy level
$H_{N}(J_{\theta},L_{z})=-\frac{1}{2}\mu^{2}$, the ratio $\nu$ depends only on
one parameter that can be chosen to be either $J_{\theta}$ or $L_{z}$. To be
more physically relevant, we introduce a new parametrization as follows. In
the Kerr case, the photon shell is foliated by bound unstable orbits with
constant radius varying in an interval $\mathcal{R}$. Recall that in equation
(4.9) we have expressed $L_{z}/E$ and $Q/E$ as functions of $r_{c}$ which is
the radius of the photo sphere. Since both $(L_{z},Q)$ and
$(L_{z},J_{\theta})$ can be used as independent integrals of the Hamiltonian
$H_{N}$, there is a map $f$ relating the two sets of integrals:
$g(Q,L_{z})=(J_{\theta},L_{z})$. Therefore $\nu\circ g$ can be considered as a
function of $r_{c}$ using (4.9). We thus obtain a reparametrization of the
ratio $\nu$ in terms of the radius $r_{c}$ of the slice of the photon shell.

(a) $a=0.8,\ M=1,\ E=0.96$

(b) $a=0.8,\ M=1,\ E=0.92$

Figure 7. The frequency ratio $\omega_{\varphi}/\omega_{\theta}$ in the
$\mu^{2}=1$ case

######  Theorem 4.4.

Let $\varepsilon h_{\mu\nu}dx^{\mu}dx^{\nu}$ be a stationary perturbation of
the Kerr metric. In either of the following two settings

  1. (1)

In the massless case $\mu^{2}=0$, there exists an open set of parameter $a/M$;

  2. (2)

In the massive case $\mu^{2}=1$, there exists an open set of parameters
$a,M,E$.

such that the following holds. Denote by $U\subset\mathcal{R}$ a small
neighborhood (may be empty) of the point of discontinuity of $\nu\circ
g(Q(r_{c}),L_{z}(r_{c})).$ For $\varepsilon$ sufficiently small, we have

  1. (1)

There is a closed subset $\mathcal{C}$ of $\mathcal{R}\setminus U$ with
measure $|\mathcal{R}\setminus(U\cup\mathcal{C})|=O(\varepsilon^{1/2})$ such
that for each $r_{c}\in\mathcal{C}$ in the perturbed system there is a photon
sphere on which all orbit is periodic or quasiperiodic with frequency
$\nu\circ g(Q(r_{c}),L_{z}(r_{c}))$ and all orbit lies in $\mathcal{M}_{\nu}$.

  2. (2)

Each open interval in $\mathcal{R}\setminus(U\cup\mathcal{C})$ corresponds to
a Birkhoff instability region where the photon spheres are broken in the sense
that $\mathcal{M}_{\nu}$ with $\nu\circ g(Q(r_{c}),L_{z}(r_{c})),\
r_{c}\in\mathcal{R}\setminus\mathcal{C},$ are all broken (see Definition B.2),
and twist map theory $($c.f. Theorem B.1 $(2.b),(3.b)$ and Theorem B.3$)$
applies.

###### Proof.

Let $f$ be the map relating the two set of constants of motions
$(J_{\theta},L_{z})$ and $(H,Q)$, i.e. $f(J_{\theta},L_{z})=(H,Q)$. From
$Df\cdot Df^{-1}=I$, we obtain

$\left[\begin{array}[]{ccc}\frac{\partial H}{\partial
J_{\theta}}&\frac{\partial H}{\partial L_{z}}\\\ \frac{\partial Q}{\partial
J_{\theta}}&\frac{\partial Q}{\partial
L_{z}}\end{array}\right]\left[\begin{array}[]{ccc}\frac{\partial
J_{\theta}}{\partial H}&\frac{\partial J_{\theta}}{\partial Q}\\\
\frac{\partial L_{z}}{\partial H}&\frac{\partial L_{z}}{\partial
Q}\end{array}\right]=I.$

This implies $Df=\frac{1}{\det
Df^{-1}}\left[\begin{array}[]{ccc}\frac{\partial L_{z}}{\partial
Q}&-\frac{\partial J_{\theta}}{\partial Q}\\\ -\frac{\partial L_{z}}{\partial
H}&\frac{\partial J_{\theta}}{\partial H}\end{array}\right].$ Denote by

$\omega=(\omega_{\theta},\omega_{\varphi})=\frac{1}{\det
Df^{-1}}\left(\frac{\partial L_{z}}{\partial Q},-\frac{\partial
J_{\theta}}{\partial Q}\right)$

the first row of $Df$. We next provide explicit formula for the derivatives in
the above expression (noting that $U$ depends on $L_{z}$ and $L_{z}$ is a
function of $Q$)

$\displaystyle\frac{\partial J_{\theta}}{\partial Q}$
$\displaystyle=\frac{1}{2\pi}\int_{\theta^{-}}^{\theta^{+}}\frac{1}{\sqrt{Q-U(\theta)}}\,d\theta-\frac{L_{z}}{\pi}\frac{\partial
L_{z}}{\partial
Q}\int_{\theta^{-}}^{\theta^{+}}\frac{\cos^{2}\theta}{\sin^{2}\theta}\frac{1}{\sqrt{Q-U(\theta)}}\,d\theta.$

We are only interested in the ratio
$\frac{\omega_{\varphi}}{\omega_{\theta}}$, then we get

(4.14) $\displaystyle\frac{\omega_{\varphi}}{\omega_{\theta}}$
$\displaystyle=-\frac{\partial J_{\theta}}{\partial Q}/\frac{\partial
L_{z}}{\partial Q}$ $\displaystyle=-\left(\frac{\partial L_{z}}{\partial
Q}\right)^{-1}\frac{1}{2\pi}\int_{\theta^{-}}^{\theta^{+}}\frac{1}{\sqrt{Q-U(\theta)}}\,d\theta+\frac{L_{z}}{\pi}\int_{\theta^{-}}^{\theta^{+}}\frac{\cos^{2}\theta}{\sin^{2}\theta}\frac{d\theta}{\sqrt{Q-U(\theta)}}.$

Up to now the formula remains the same formally for the $\mu^{2}=0$ and
$\mu^{2}=1$ cases. We next express $\frac{\omega_{\varphi}}{\omega_{\theta}}$
as a function of $r_{c}$. It is enough to express $L_{z}$ and $Q$ as a
function of $r_{c}$ in (4.9) for the $\mu^{2}=1$ case and (4.8) for the
$\mu^{2}=0$ case. We get $\frac{\partial L_{z}}{\partial
Q}=\frac{d\ell}{dr}/\frac{dq}{dr}$ from (4.9) for the $\mu^{2}=1$ case (4.8)
for the $\mu^{2}=0$ case. Thus the ratio $\omega_{\varphi}/\omega_{\theta}$ is
a function of $r_{c}$ in both cases.

In the massless case, the ratio $\omega_{\varphi}/\omega_{\theta}$ is
independent of $E$. Figure 1 is plotted with choice of parameter $a/M=0.8$. In
the massive case, we plot Figure 7 with parameter $a=0.8,M=1,E=0.96$ and
Figure 7 with parameter $a=0.8,M=1,E=0.92$, corresponding to Figure 5 and 5
respectively. We note that Figure 7 is discontinuous which is related to the
fact that $L_{z}$ changes sign, and we have that 7 is continuous since $\ell$
does not change sign in Figure 5. We shall explain the discontinuity in the
next section. Also there is a critical point in 7, so we can only apply KAM
and twist map theory outside a small neighborhood of the critical point. ∎

####  4.3.2. The point of discontinuity

There is an interesting feature in the graph in Figure 1 and Figure 7 that is
the appearance of a discontinuity point. The point of discontinuity occurs at
the place where $\ell(r_{c})=0$ in (4.9). This means that the orbit becomes
from prograde to retrograde as $r$ crosses the zero of $L_{z}$ when
increasing. This phenomenon was noticed and explained in [Wi, T]. Let us give
a more mathematical explanation.

We see from Section 4.1 that the vertical dynamics in the $L_{z}=0$ case
differs drastically from the $L_{z}\neq 0$ case. When $L_{z}=0$ and
$\bar{Q}:=Q-(\mu^{2}-E^{2})a^{2}>0$, the latitudinal motion $\theta$ ranges
over $[0,\pi]$, in particular, the orbit passes through the south and north
pole. For $|L_{z}|\neq 0$ but sufficient small, the latitudinal angle
$\theta\in[\delta,\pi-\delta]$ where $\delta$ solves the equation
$Q-U(\delta)=0$. In particular, as $L_{z}\to 0$ we have
$\frac{L_{z}^{2}}{\sin^{2}\delta}\to\bar{Q}$.

The discontinuity originates from the second term on the RHS of (4.14).
Another way to see it is to evaluate the quantity $\Delta\varphi$ that is the
change of the azimuthal angle during each period of the $\theta$-motion. We
have

$\Delta\varphi=\oint\frac{d\varphi}{d\theta}d\theta=2\int^{\pi-\delta}_{\delta}\frac{\frac{a}{\Delta}P(r)-\left(aE-\frac{L_{z}}{\sin^{2}\theta}\right)}{\sqrt{Q-U(\theta)}}d\theta.$

We shall now work on $\Delta\varphi$ and it will turn out that the the
discontinuity in $\Delta\varphi$ agrees with the second term of the RHS of
4.14.

We note that
$\frac{d\varphi}{d\lambda}=\frac{L_{z}}{\sin^{2}\theta}+O(1)=\frac{\bar{Q}}{L_{z}}+O(1)$
as $\theta\to\delta$ or $\pi-\delta$. Then we have
$\frac{d\varphi}{d\theta}\simeq\frac{L_{z}}{\sin^{2}\theta}\frac{1}{\sqrt{Q-U}}\simeq\frac{L_{z}}{\sin^{2}\theta\sqrt{\bar{Q}-\frac{L_{z}^{2}}{\sin^{2}\theta}}}$
when $\theta$ is $c\sqrt{\delta}$-close to $\delta$ or $\pi-\delta$ for large
number $c$. For $\theta\in[c\sqrt{\delta},\pi-c\sqrt{\delta}]$, the $L_{z}$
dependence in the integrand of $\Delta\varphi$ is negligible so the integral
$\int_{c\sqrt{\delta}}^{\pi-c\sqrt{\delta}}\frac{d\varphi}{d\theta}d\theta$ is
almost the same for $L_{z}>0$ and $L_{z}<0$ with small modulus. The
contribution to $\Delta\varphi$ when $\theta$ varies in the interval
$[\delta,c\sqrt{\delta}]$ is given by

$\displaystyle
L_{z}\int_{\delta}^{c\sqrt{\delta}}\frac{d\theta}{\sin^{2}\theta\sqrt{\bar{Q}-\frac{L_{z}^{2}}{\sin^{2}\theta}}}\simeq
L_{z}\int_{\delta}^{c\sqrt{\delta}}\frac{d\theta}{\theta\sqrt{\theta^{2}\bar{Q}-\delta^{2}\bar{Q}}}$
$\displaystyle=$
$\displaystyle\mathrm{sgn}(L_{z})\int_{1}^{c/\sqrt{\delta}}\frac{dx}{x\sqrt{x^{2}-1}}\simeq\mathrm{sgn}(L_{z})\int_{0}^{1}\frac{dy}{\sqrt{1-y^{2}}}$
$\displaystyle=$ $\displaystyle\mathrm{sgn}(L_{z})\frac{\pi}{2}$

where we substitute $L_{z}\simeq\sqrt{\bar{Q}}\delta,\ \sin\theta\simeq\theta$
and the $\simeq$ becomes $=$ when we take limit $L_{z}\to 0$ hence $\delta\to
0$. During each period of the $\theta$ motion, this discontinuity should be
counted four times, thus the total discontinuity when $L_{z}$ changes sign is
$4\pi$. In particular, when $L_{z}$ changes sign, the orbit changes from
prograde to retrograde. Similar calculation also occurs in (4.14). Since we
use a different normalization than $\Delta\varphi$, the discontinuity in
(4.14) is 2.

We remark that for some choices of parameters $a,E$, the function $\ell(r)$ in
(4.9) may be strictly positive (see Figure 5), in which case the ratio $\nu$
does not change sign.

###  4.4. KAM and QPO

In this section, we consider only the timelike case with $\mu^{2}=1$.

Suppose $E^{2}<1$ and $-R$ attains its local max at $r^{*}(>r_{h})$ and a
local min at $r_{*}>r^{*}$. Suppose also $-R(r^{*})>0$ and $-R(r_{*})<0$
(imagine that we shift the graph of $-R$ in Figure 6 upward slightly). Let
$C>0$ be a large constant, and we define the following bounded part of the
phase space where all Kerr orbits are bounded hence Liouville-Arnold theorem
applies

(4.15) $\mathcal{B}(C):=\left\\{E^{2}<1,\ |L_{z}|<C,\ 0<Q<C\ |\ -R(r^{*})>0,\
-R(r_{*})<0\right\\}\cap H^{-1}(-\frac{1}{2}).$

This is a neighborhood of the stable circular orbit. Denote by $r^{-}<r^{+}$
the two largest roots of $-R(r)=0$, then we have $r^{*}<r_{-}<r_{*}<r_{+}$. In
the $(\theta,p_{\theta})$-component, the dynamics is again periodic around the
point $(\theta,p_{\theta})=(\pi/2,0)$ as we have analyzed in Section 4.1.

Again, similar to the Schwarzchild case, we consider for simplicity
perturbations that are stationary and axisymmetric, thus $E$ and $L_{z}$
remain constants of motion in the perturbed system. Therefore we fix values of
$L_{z}$ and $E$, and introduce the action-angle coordinates
$(\alpha_{r},J_{r},\alpha_{\theta},J_{\theta})\in\mathbb{T}\times\mathbb{R}_{+}\times\mathbb{T}\times\mathbb{R}_{+}$
as follows

$\displaystyle J_{r}$ $\displaystyle=\frac{1}{2\pi}\oint
p_{r}\,dr=\frac{1}{\pi}\int_{r^{-}}^{r^{+}}\frac{\sqrt{R(r)}}{\Delta(r)}\,dr,$
$\displaystyle J_{\theta}$ $\displaystyle=\frac{1}{2\pi}\oint
p_{\theta}\,d\theta=\frac{1}{\pi}\int_{\theta^{-}}^{\theta^{+}}\sqrt{Q-U(\theta)}\,d\theta,$

where $r^{-}<r^{+}$ and $\theta^{-}<\theta^{+}$ are zeros of the integrands
above respectively. Again we get the Hamiltonian written in these coordinates
$H=H(J_{r},J_{\theta})$ and we have the $\omega_{r}:=\frac{\partial
H}{\partial J_{r}}$ the radial frequency and $\omega_{\theta}:=\frac{\partial
H}{\partial J_{\theta}}$ the latitudinal frequency. Both the expressions of
$J_{r}$ and $J_{\theta}$ depend on $Q$, so we can write the frequency ratio
$\omega_{r}/\omega_{\theta}$ as a function of $Q$.

######  Theorem 4.5.

Let $\varepsilon h_{\mu\nu}dx^{\mu}dx^{\nu}$ be a stationary and axisymmetric
perturbation to the Kerr metric. Then there exists an open set of parameters
$a,E,L_{z}$ such that there are a neighborhood $\mathcal{D}$ of the point
$(r,p_{r},\theta,p_{\theta})=(r_{*},0,\pi/2,0)$ and $\varepsilon_{0}$ such
that KAM and the twist map theory are applicable to the perturbed system in
$\mathcal{D}$ and for all $0<|\varepsilon|<\varepsilon_{0}$.

Figure 8. The ratio $\omega_{r}/\omega_{\theta}$ for Kerr with parameters
$a=0.8,E=0.96,L_{z}=1.5$

Figure 8 is plotted with parameters $a=0.8,E=0.96,L_{z}=1.5$.

An example of such perturbations is given by the Monko-Novikov metric [GAC],
which is so complicated that we do not cite it here. In [GAC] the authors
analyzed the $1:2$ and $2:3$ resonances and discovered that the dynamics of
the perturbed Kerr is chaotic. The frequency ratio
$\omega_{r}/\omega_{\theta}$ has been studied by many authors, c.f. [BGH, S]
etc.

###### Proof of Theorem 4.5.

Let $f$ be the map relating two sets of constants of motion
$(J_{r},J_{\theta})$ and $H,Q$: $f(J_{r},J_{\theta})=(H,Q)$. From $Df\cdot
Df^{-1}=I$, we obtain

$\left[\begin{array}[]{ccc}\frac{\partial H}{\partial J_{r}}&\frac{\partial
H}{\partial J_{\theta}}\\\ \frac{\partial Q}{\partial J_{r}}&\frac{\partial
Q}{\partial
J_{\theta}}\end{array}\right]\left[\begin{array}[]{ccc}\frac{\partial
J_{r}}{\partial H}&\frac{\partial J_{r}}{\partial Q}\\\ \frac{\partial
J_{\theta}}{\partial H}&\frac{\partial J_{\theta}}{\partial
Q}\end{array}\right]=I.$

This implies $Df=\frac{1}{\det
Df^{-1}}\left[\begin{array}[]{ccc}\frac{\partial J_{\theta}}{\partial
Q}&-\frac{\partial J_{r}}{\partial Q}\\\ -\frac{\partial J_{\theta}}{\partial
H}&\frac{\partial J_{r}}{\partial H}\end{array}\right].$ Denote by

$\omega=(\omega_{r},\omega_{\theta})=\frac{1}{\det
Df^{-1}}\left(\frac{\partial J_{\theta}}{\partial Q},-\frac{\partial
J_{r}}{\partial Q}\right)$

the first row of $Df$. In our setting, we fix $H=-1/2,E<1$ and $L_{z}$ as
constants, then the only remaining parameter is $Q$. We next provide explicit
formula for the derivatives in the above expression

(4.16) $\displaystyle\frac{\partial J_{r}}{\partial H}$
$\displaystyle=\frac{1}{\pi}\int_{r^{-}}^{r^{+}}\frac{r^{2}}{\sqrt{R}}\,dr,\quad\frac{\partial
J_{r}}{\partial Q}=-\frac{1}{2\pi}\int_{r^{-}}^{r^{+}}\frac{1}{\sqrt{R}}\,dr,$
$\displaystyle\frac{\partial J_{\theta}}{\partial H}$
$\displaystyle=\frac{2}{\pi}\int_{\theta^{-}}^{\theta^{+}}\frac{a^{2}\cos^{2}\theta}{\sqrt{Q-U(\theta)}}\,d\theta,\quad\frac{\partial
J_{\theta}}{\partial
Q}=\frac{1}{2\pi}\int_{\theta^{-}}^{\theta^{+}}\frac{1}{\sqrt{Q-U(\theta)}}\,d\theta.$

Then we get

$\frac{\omega_{r}}{\omega_{\theta}}=-\frac{\partial J_{\theta}}{\partial
Q}/\frac{\partial J_{r}}{\partial
Q}=\int_{\theta^{-}}^{\theta^{+}}\frac{1}{\sqrt{Q-U(\theta)}}\,d\theta\Big{/}\int_{r^{-}}^{r^{+}}\frac{1}{\sqrt{R}}\,dr.$

Figure 8 is then plotted with this formula. This completes the proof of the
theorem.

∎

Finally, for the sake of completeness, we provide formula for the verification
of the KAM nondgeneracy conditions in Section 2.2. This part is similar to
Section 3.3.2 so readers are suggested to skip it when first reading. We need
to evaluate $\nabla\omega$ where $\nabla$ is the derivative with respect to
$J_{r},J_{\theta}$. We thus arrive at

$\displaystyle\nabla\omega$
$\displaystyle=\frac{\partial(H,Q)}{\partial(J_{r},J_{\theta})}(\partial_{H},\partial_{Q})\omega=Df\cdot(\partial_{H},\partial_{Q})\left(\frac{1}{\det
Df^{-1}}\left(\frac{\partial J_{\theta}}{\partial Q},-\frac{\partial
J_{r}}{\partial Q}\right)\right)$ $\displaystyle=\left(\frac{1}{\det
Df^{-1}}Df\right)(\partial_{H},\partial_{Q})\left(\frac{\partial
J_{\theta}}{\partial Q},-\frac{\partial J_{r}}{\partial Q}\right)$
$\displaystyle+Df\cdot\left((\partial_{H},\partial_{Q})\left(\frac{1}{\det
Df^{-1}}\right)\right)\otimes\left(\frac{\partial J_{\theta}}{\partial
Q},-\frac{\partial J_{r}}{\partial Q}\right).$

We can then compute $\nabla\omega$ similar to Section 3.3.2.

If we are interested in perturbations also depending on the angle $\varphi$,
then we have a Hamiltonian system of three degrees of freedom and we have to
consider three fundamental frequencies as in Section 3.3.3. The procedure is
similar to the Schwarzschild case in Section 3.3.3, though the computations
are more involved. We skip the details here.

##  Appendix A The theorem of normally hyperbolic invariant manifold

In this section we give the version of normally hyperbolic invariant manifold
theorem that we used in the proof of our main theorem. The standard references
are [HPS, F].

######  Definition A.1 (NHIM).

Let $N\subset M$ be a submanifold $($maybe noncompact$)$ invariant under $f$,
$f(N)=N$. We say that $N$ is a normally hyperbolic invariant manifold if there
exist a constant $C>0$, rates $0<\lambda<\mu^{-1}<1$ and an invariant $($under
$Df)$ splitting for every $x\in N$

$T_{x}M=E^{s}(x)\oplus E^{u}(x)\oplus T_{x}N$

in such a way that

$\displaystyle v\in E^{s}(x)$
$\displaystyle\Leftrightarrow\quad|Df^{n}(x)v|\leq C\lambda^{n}|v|,\quad n\geq
0,$ $\displaystyle v\in E^{u}(x)$
$\displaystyle\Leftrightarrow\quad|Df^{n}(x)v|\leq C\lambda^{|n|}|v|,\ \ n\leq
0,$ $\displaystyle v\in T_{x}N$
$\displaystyle\Leftrightarrow\quad|Df^{n}(x)v|\leq C\mu^{n}|v|,\quad
n\in\mathbb{Z}.$

Here the Riemannian metric $|\cdot|$ can be any prescribed one, which may
change the constant $C$ but not $\lambda,\mu$.

######  Theorem A.2.

Suppose $N$ is a NHIM under the $C^{r}$, $r>1$, diffeomorphism $f:\ M\to M$.
Denote $\ell=\min\\{r,\frac{|\ln\lambda|}{|\ln\mu|}\\}$. Then for any $C^{r}$
$f_{\epsilon}$ that is sufficiently close to $f$ in the $C^{1}$ norm,

  1. (1)

there exists a NHIM $N_{\epsilon}$ that is a $C^{\ell}$ graph over $N$,

  2. (2)

$($Invariant splitting$)$ There exists a splitting for $x\in N_{\epsilon}$

$T_{x}M=E^{u}_{\epsilon}(x)\oplus E^{s}_{\epsilon}(x)\oplus T_{x}N_{\epsilon}$

invariant under the map $f_{\epsilon}$. The bundle $E^{u,s}_{\epsilon}(x)$ is
$C^{\ell-1}$ in $x$.

  3. (3)

There exist $C^{\ell}$ stable and unstable manifolds $W^{s}(N_{\epsilon})$ and
$W^{u}(N_{\epsilon})$ that are invariant under $f$ and are tangent to
$E_{\epsilon}^{s}\oplus TN_{\epsilon}$ and $E_{\epsilon}^{u}\oplus
TN_{\epsilon}$ respectively.

  4. (4)

The stable and unstable manifolds $W^{u,s}(N_{\epsilon})$ are fibered by the
corresponding stable and unstable leaves $W^{u,s}_{x,\epsilon}:=\\{y\
|d(f_{\varepsilon}^{n}(y),f_{\varepsilon}^{n}(x))\to 0,n\to+\infty\
\mathrm{for}\ s,n\to-\infty\ \mathrm{for}\ s\\}$:

$W^{u}(N_{\epsilon})=\cup_{x\in N_{\epsilon}}W_{x,\epsilon}^{u},\quad
W^{s}(N_{\epsilon})=\cup_{x\in N_{\epsilon}}W_{x,\epsilon}^{s}.$

  5. (5)

The maps $x\mapsto W_{x,\epsilon}^{u,s}$ are $C^{\ell-j}$ when
$W_{x,\epsilon}^{u,s}$ is given $C^{j}$ topology.

  6. (6)

If $f$ and $f_{\epsilon}$ are Hamiltonian, then $N_{\epsilon}$ is symplectic
and the map $f_{\epsilon}$ restricted to $N_{\epsilon}$ is also Hamiltonian.

##  Appendix B Aubry-Mather theory for twist maps

The dynamics of twist maps are studied in details in Aubry-Mather theory. Let
us summarize the main result and refer readers to [B] for more detailed
proofs. The theory can be put in a very general setting. For the sake of
concreteness, we consider the following setting. Let $H:\
T^{*}\mathbb{T}^{1}\times\mathbb{T}^{1}\to\mathbb{R}$ be a nonautonomous
1-periodic Hamiltonian satisfying

  1. (1)

$H(x,y,t)=H(x+1,y,t)=H(x,y,t+1)$;

  2. (2)

$\partial_{y}^{2}H>0$ everywhere;

  3. (3)

$\frac{H(x,y,t)}{|y|}\to\infty$ as $|y|\to\infty$ for all
$(x,t)\in\mathbb{T}^{2}$;

  4. (4)

The flow $\phi^{t}$ is complete, i.e. for each initial condition
$z_{0}=(x_{0},y_{0})$, the orbit $\phi^{t}(z_{0})$ is defined for all
$t\in\mathbb{R}$.

The time-1 map, $\phi^{1}:\ T^{*}\mathbb{T}^{1}\to T^{*}\mathbb{T}^{1}$ of the
Hamiltonian flow is symplectic and is the composition of twist maps [M1].

Let $H:\ T^{*}\mathbb{T}^{1}\times\mathbb{T}^{1}\to\mathbb{R}$ be as above and
we introduce its dual Lagrangian by taking Legendre transform
$L(x,\dot{x},t):=\sup_{y}(\dot{x}\cdot y-H(x,y,t))$, and for any curve
$\gamma:\ [0,n]\to\mathbb{T}^{1},\ n\in\mathbb{Z}$, we define its action

$A(\gamma)=\int_{0}^{n}L(\gamma(t),\dot{\gamma}(t),t)\,dt.$

A curve is a critical point of the action with fixed endpoints if and only if
it solves the Euler-Lagrange equation $\frac{d}{dt}\frac{\partial
L}{\partial\dot{x}}(\gamma(t),\dot{\gamma}(t),t)=\frac{\partial L}{\partial
x}(\gamma(t),\dot{\gamma}(t),t)$. If the curve $\gamma:\
\mathbb{R}\to\mathbb{T}^{1}$ solves the Euler-Lagrange equation, then the
Legendre transform ($\\{(\gamma(t),y(\gamma(t)))\\}\subset
T^{*}\mathbb{T}^{1}$) of ($\\{\gamma(t),\dot{\gamma}(t))\\}\subset
T\mathbb{T}^{1}$) solves the Hamiltonian equation, then the time-1 map
$\phi^{1}$ maps $(\gamma(n),y(\gamma(n)))$ to $(\gamma(n+1),y(\gamma(n+1)))$
for all $n\in\mathbb{Z}$. We say that a $C^{1}$ curve $\gamma:\
\mathbb{R}\to\mathbb{T}^{1}$ is globally minimizing, if for any
$r,s\in\mathbb{R}$ and any $\xi:\ [r^{\prime},s^{\prime}]\to\mathbb{T}^{1}$
with $r^{\prime}-r\in\mathbb{Z},s^{\prime}-s\in\mathbb{Z}$, we have

$\int_{r}^{s}L(\gamma(t),\dot{\gamma}(t),t)\,dt\leq\int_{r^{\prime}}^{s^{\prime}}L(\xi(t),\dot{\xi}(t),t)\,dt.$

Since the Lagrangian $L$ is $1$-periodic in its first entry, we may lift each
curve $\gamma:\ \mathbb{R}\to\mathbb{T}^{1}$ to a curve $\tilde{\gamma}:\
\mathbb{R}\to\mathbb{R}$ by unfolding $\mathbb{T}^{1}$ to its universal
covering space $\mathbb{R}$, so we can visualize the graph of $\tilde{\gamma}$
in $\mathbb{R}^{2}$.

The following theorem summarizes the main results of the twist map theory
(c.f. [B]).

######  Theorem B.1.

  1. (1)

Each globally minimizing curve $\gamma$ has a well-defined rotation number
$\rho=\lim_{n\to\pm\infty}\frac{\tilde{\gamma}(n)}{n}$. For each $\rho$, the
set of globally minimizing curves with rotation number $\rho$ is not empty, we
introduce the set

$\displaystyle\mathcal{M}_{\rho}$
$\displaystyle=\\{(\gamma(n),\dot{\gamma}(n))\in T\mathbb{T}^{1},\
n\in\mathbb{Z}\ |$ $\displaystyle\ \gamma:\ \mathbb{R}\to\mathbb{T}^{1}\
\mathrm{is\ globally\ minimizing\ with\ rotation\ number}\ \rho\\}.$

We also use the same notation to denote its Legendre transform to
$T^{*}\mathbb{T}^{1}$.

  2. (2)

If the rotation number $\rho$ is irrational, $\mathcal{M}_{\rho}$ is

    1. (a)

either an invariant circle (a homotopically nontrivial circle on
$T^{*}\mathbb{T}^{1}$),

    2. (b)

or a Denjoy minimal set (a Cantor subset of a homotopically nontrivial circle
on $T^{*}\mathbb{T}^{1}$).

  3. (3)

If the rotation number $\rho$ is rational, $\mathcal{M}_{\rho}$ always
contains periodic orbits.

    1. (a)

It may be an invariant circle consists of periodic orbits of the same period.

    2. (b)

If $\mathcal{M}_{\rho}$ is not an invariant circle, then in the gap between
two neighboring periodic orbits $\tilde{\gamma}^{-},\,\tilde{\gamma}^{+}$ in
$\mathbb{R}^{2}$ where
$\\{(\gamma^{\pm}(n),\dot{\gamma}^{\pm}(n))\\}\subset\mathcal{M}_{\rho}$,
there is a globally minimizing curve $\gamma$ with rotation number $\rho$
whose lift is asymptotic to $\tilde{\gamma}^{-}$ in the future and to
$\tilde{\gamma}^{+}$ in the past, and another asymptotic to
$\tilde{\gamma}^{+}$ in the future and to $\tilde{\gamma}^{-}$ in the past.

######  Definition B.2.

When $\mathcal{M}_{\rho}$ is not a homotopically nontrivial invariant circle
in $T^{*}\mathbb{T}^{1}$, we say that $\mathcal{M}_{\rho}$ is _broken_. This
includes case (2.b) and (3.b) in the above Theorem B.1.

Let $R$ be the region (called Birkhoff region of instability) bounded by
$\phi^{1}$-invariant homotopically non-trivial Jordan curves $\Gamma_{-}$ and
$\Gamma_{+}$ with $\rho(\Gamma_{-})<\rho(\Gamma_{+})$, where
$\rho(\Gamma_{\pm})$ is the rotation number of $\phi^{1}|_{\Gamma_{\pm}}$.
Note that all $\mathcal{M}_{\rho}\subset R$ with
$\rho(\Gamma_{-})<\rho<\rho(\Gamma_{+})$ is broken. We cite the following
theorem of Mather.

######  Theorem B.3 ([Ma2]).

  1. (1)

Suppose $\rho(\Gamma_{-})<\alpha,\omega<\rho(\Gamma_{+})$. Then there is an
orbit of $\phi^{1}$ whose $\alpha$-limit set lies in $\mathcal{M}_{\alpha}$
and whose $\omega$-limit set lies in $\mathcal{M}_{\omega}$. Furthermore, if
$\rho(\Gamma_{-})$ $($resp. $\rho(\Gamma_{+}))$ is irrational, then this
conclusion still holds with the weaker hypothesis
$\rho(\Gamma_{-})\leq\alpha$, (resp. $\alpha,\ \omega\leq\rho(\Gamma_{+}))$.

  2. (2)

Consider for each $i\in\mathbb{Z}$ a real number
$\rho(\Gamma_{-})\leq\omega_{i}\leq\rho(\Gamma_{+})$ and a positive number
$\varepsilon_{i}$. Then there exists an orbit $(\ldots,P_{j},\ldots)$ and an
increasing bi-infinite sequence of integers $j(i)$ such that dist.
$(P_{j(i)},),\mathcal{M}_{\omega(i)})<\varepsilon_{i}$.

## Acknowledgment

The author would like to express his deep gratitudes to Professor S.-T. Yau
for suggesting the topic. He also would like to thank Mr. Yifan Guo for
drawing all the figures. The author is supported by NSFC (Significant project
No.11790273) in China and Beijing Natural Science Foundation (Z180003).

## References

  * [1]
  * [A89] Arnol’d, V.I. Mathematical methods of classical mechanics. Vol. 60. Springer, 1989.
  * [AKN] Arnold, Vladimir I., Valery V. Kozlov, and Anatoly I. Neishtadt. Mathematical aspects of classical and celestial mechanics. Vol. 3. Springer Science & Business Media, 2007.
  * [B] Bangert, Victor. Mather sets for twist maps and geodesics on tori. Dynamics reported. Vieweg+ Teubner Verlag, 1988. 1-56.
  * [BGH] Brink, Jeandrew, Marisa Geyer, and Tanja Hinderer. The Astrophysics of Resonant Orbits in the Kerr Metric. arXiv preprint arXiv:1501.07728 (2015).
  * [C] Chandrasekhar, Subrahmanyan. The mathematical theory of black holes. Research supported by NSF. Oxford/New York, Clarendon Press/Oxford University Press (International Series of Monographs on Physics. Volume 69), 1983, 663 p. 1 (1983).
  * [Ca] Dias Carneiro, On minimizing measures of the action of autonomous Lagrangians, Nonlinearily, 8 (1995)1077-1085.
  * [DKN] Efthimia Deligianni, Jutta Kunz, Petya Nedkova, Quasi-periodic oscillations from the accretion disk around distorted black holes, Arxiv: 2003.01252.
  * [F] Fenichel N., Persistence and smoothness of invariant manifolds for flows, Indiana Univ. Math. J 21 (1971) 193-226.
  * [GAC] Lukes-Gerakopoulos, Georgios, Theocharis A. Apostolatos, and George Contopoulos. Observable signature of a background deviating from the Kerr metric. Physical Review D 81.12 (2010): 124005.
  * [GHW] Gralla, Samuel E., Daniel E. Holz, and Robert M. Wald. Black hole shadows, photon shells, and lensing rings. Physical Review D 100.2 (2019): 024018.
  * [GL] Gralla, Samuel E., and Alexandru Lupsasca. Null geodesics of the Kerr exterior. Physical Review D 101.4 (2020): 044032.
  * [GLP] Grossman, Rebecca, Janna Levin, and Gabe Perez-Giz. Harmonic structure of generic Kerr orbits. Physical Review D 85.2 (2012): 023012
  * [HF] Hinderer, Tanja, and Eanna E. Flanagan. Two-timescale analysis of extreme mass ratio inspirals in Kerr spacetime: Orbital motion. Physical Review D 78.6 (2008): 064028.
  * [HPS] Hirsch M. Pugh C. & Shub M., Invariant Manifolds, Lecture Notes Math. 583 (1977) Springer-Verlag.
  * [IM] Adam R. Ingram, Sara E. Motta, A review of quasi-periodic oscillations from black hole X-ray binaries: observation and theory, New Astronomy Reviews 85 (2019) 101524
  * [J] Michael D. Johnson, Alexandru Lupsasca, Andrew Strominger, et al, Universal interferometric signatures of a black hole’s photon shell. arXiv:1907.04329v2 [astro-ph.IM] 27 Mar 2020
  * [JDW] Chichuan Jin, Chris Done, and Martin Ward, Reobserving the NLS1 galaxy RE J1034+396 - I. The long-term, recurrent X-ray QPO with a high significance, MNRAS 495, 3538-3550 (2020).
  * [M1] Moser, J, Monotone twist mappings and calculus of variations, Erg. Theory Dyn. Sys, 1986, 6, 401-413.
  * [M2] Moser, J. On invariant curves of area-preserving mappings of an annulus. Nachr. Akad. Wiss. Göttingen, II (1962): 1-20.
  * [M87] Wielgus, M, et al, Monitoring the Morphology of M87* in 2009-2017 with the Event Horizon Telescope, Astrophysics Journal, vol 901, no 67, 2020.
  * [Ma1] Mather J., Action minimizing invariant measures for positive definite Lagrangian systems, Math. Z., 207(2)(1991) 169-207.
  * [Ma2] Mather, John N. Variational construction of orbits of twist diffeomorphisms. Journal of the American Mathematical Society (1991): 207-263.
  * [N] New York Times, Infinite Visions Were Hiding in the First Black Hole Image’s Rings, https://www.nytimes.com/2020/03/28/science/black-hole-rings.html.
  * [P] Pöschel, Jürgen. A lecture on the classical KAM theorem. arXiv preprint arXiv:0908.2234 (2009).
  * [S] Schmidt, Wolfram. Celestial mechanics in Kerr spacetime. Classical and Quantum Gravity 19.10 (2002): 2743.
  * [T] Teo, Edward. Spherical photon orbits around a Kerr black hole. General Relativity and Gravitation 35.11 (2003): 1909-1926.
  * [WBS] Warburton, Niels, Leor Barack, and Norichika Sago. Isofrequency pairing of geodesic orbits in Kerr geometry. Physical Review D 87.8 (2013): 084012.
  * [Wi] D. Wilkins, Bounded geodesics in the Kerr metric, Physical Review D, Vol 5, No. 4, 814-822, 1972,
  * [X] J. Xue, Arnold diffusion and geodesic dynamics of blackholes, preprint.