\section{Introduction}\label{sec-I}

In a few years a new era of astronomy is expected to begin when the Advanced LIGO
and other gravitational wave detectors are anticipated to receive the first signals.
The main source of gravitational wave signals are binary systems consisting of
compact objects whose mass proportion ranges from equal masses to extreme mass
ratio. The gravitational wave signals from binary systems are closely related to
the motion of the inspiraling compact objects. Thus, orbital studies of binary
systems have attracted much attention (e.g., \cite{Levinetal,Hackmann}). In general,
orbital motion can be studied by using a Hamiltonian function which splits into an
orbital and a spin part.

One possible approximation is the post-Newtonian (e.g., \cite{PostNewton} and
references therein). In the post-Newtonian approximation the Hamiltonian's
orbital part is expanded in powers of $\frac 1 {c^2}$. This approach has also
been extended to binaries that are perturbed by a third, much lighter object
(e.g., \cite{3body}). These post-Newtonian approaches have in common that
the mass ratio of the binary should not be very  large.

When the binary system consists of compact objects of extreme mass ratio,
the motion in such a system is called an Extreme Mass Ratio Inspiral (EMRI). We
expect to find such binaries at the center of galaxies, where a comparably light
compact object inspirals into a central Supermassive Black Hole (SMBH). The
lighter inspiraling compact object moves adiabatically from one geodesic orbit
to another. Thus, an EMRI can be approximated by the geodesic
motion of a test particle in the spacetime background of the central
supermassive object (e.g., \cite{Amaro12} and references therein).

One interesting aspect of EMRI is that the shape of the spacetime background is
encoded in the gravitational waves emitted by the inspiraling object. Ryan
showed that we can, in principle, extract this information from the
gravitational waves \cite{Ryan}. Therefore, when future low frequency
space-borne gravitational wave detectors are launched, we will be able to test
whether the spacetime backgrounds around SMBHs are indeed described by the Kerr
metric or by something else \cite{Amaro12}. In order to provide such a test,
Collins and Hughes \cite{Collins04} constructed a perturbed Schwarzschild
spacetime background and called it a \textit{bumpy black hole}. Since then,
several authors have suggested different tests \cite{Glampedakis06,Barausse07,
Gair08,Apostolatos09,Vigeland10a,Vigeland10b,Lukes10,JohanPsaltis,BambiBar,
Vigeland11,Gair11,Bambi,Johannsen12,Pappas12} using various forms of bumpy black
holes (e.g., \cite{ManNov92,MMS,MSM}), which are also known as non-Kerr metrics.

Usually, these non-Kerr metrics are stationary and axisymmetric, but they miss
the fourth integral which would make them integrable systems, as the Carter
constant \cite{Carter68} does in the case of the Kerr metric. Thus, in
general the non-Kerr metrics correspond to non-integrable systems, which in
practical calculations often are investigated for chaos (e.g., \cite{Lukes10,Wu,
Dubeibe07,Han08,Semerak,Lukes}). This, in turn, implies numerical simulations
over long time intervals. Therefore, high-efficient integrators with a good
long term behavior are required.

Over the last decades, part of the Numerical Analysis community has focused on
long-term integrations of differential equations. Geometric or structure
preserving algorithms, such as symplectic methods for Hamiltonian equations of
motions or symmetric integrators for reversible systems, have been developed
aplenty. Contrary to standard explicit integrators such as
Runge-Kutta schemes, these integrators nearly conserve first integrals (e.g.,
the energy) and overall integration errors increase only slowly with time.
Whereas for standard integration schemes the overall error is normally
proportional to the square of the length of the integration interval $t_i$, it
only increases linearly with $t_i$ for structure preserving integrators. A
detailed discussion of such methods can be found in \cite{hairerlubichwanner}.

The efficiency of numerical integrators can be increased by the use of adaptive
step size controllers. But, by using standard step size controllers for
geometric algorithms, we lose the integrators' good long term behavior. However,
over the last 15 years, algorithms allowing structure preserving step size
control have been developed (\cite{stoffer,soederlind}). This development has
resulted in a lot of available tools to confront problems in classical astronomy.

But, in general relativity, the structure of the equations is more complex.
The inclusion of spins to a binary system leads to a Poisson-system, which can
be transformed to symplectic form as stated by the \textit{Darboux-Lie-Theorem}.
Such an approach was proposed in \cite{wuxie}. After such a transformation,
symplectic  schemes can finally be applied, e.g. \cite{zhongwu}. But also for the
the post-Newtonian formulation as a Poisson-system, a non-canonically symplectic
integrator has been found, see \cite{lubich}. In the geodesic approach, the
Hamiltonian of the orbit has a position dependent mass-matrix. Thus, from this
point of view, the system of differential equations is non-separable. Furthermore,
the vector field of the resulting equations of motion is strongly varying in
non-integrable spacetime backgrounds. These two features pose a new set of problems
for the implementation of efficient and accurate algorithms. As an answer to these
problems we introduce a new integration scheme for such differential equations
based on Gauss-Runge-Kutta collocation methods. In particular, we use a new
structure preserving adaptive step size controller. As a result, we get a symmetric,
reversible integrator which efficiently conserves the constants of motion. We
tested this integrator's performance in the case of the Manko, Sanabria-G\'{o}mez,
Manko (MSM) metric \cite{MSM} and the results show that the new integration scheme
is much faster than an adaptive step fifth-order Cash-Karp Runge-Kutta scheme.

The paper is organized as follows: Section \ref{sec-II} summarizes some
basic elements about the MSM metric, and the geodesic motion in the MSM spacetime
background. In the subsequent Section \ref{sec-III}, we consider the topic from
a Numerical Analysis point of view and discuss the characteristics of the
geodesic equations. Our step size controller, along with the integrator, is
presented in section \ref{sec-IV}. The performance of the integrator in
numerical tests is shown in section \ref{sec-V}. Section \ref{sec-VI} summarizes
the main results of the study.

\section{Manko, Sanabri-G\'{o}mez, Manko spacetime background} \label{sec-II}

 \subsection{The metric} \label{subsec-II-a}

 We tested the various integration schemes on the five-parameter vacuum solution
 introduced by Manko et al. \cite{MSM}. The MSM solution is asymptotically flat,
 axisymmetric and stationary, and describes the ``exterior field of a charged,
 magnetized, spinning deformed mass'' \cite{MSM}. The MSM spacetime depends on
 five real parameters, namely the mass $m$, the spin (per unit mass)
 $a$, the total charge $q$, the magnetic dipole moment ${\cal M}$, and the
 mass-quadrupole moment ${\cal Q}$. ${\cal M}$ and ${\cal Q}$ are represented in
 the metric by the real parameters  $\mu$ and $b$, i.e.
 \begin{eqnarray} \label{eq:MagQuad}
  {\cal M} &=& \mu+q (a-b)~~, \nonumber \\
  {\cal Q} &=& -m(d-\delta-a~b+a^2)~~,
 \end{eqnarray}
 where
 \begin{eqnarray} \label{eq:deltad}
   \delta &:=& \frac{\mu^2-m^2 b^2}{m^2-(a-b)^2-q^2}~~, \nonumber \\
   d &:=& \frac{1}{4}[m^2-(a-b)^2-q^2]~~.
 \end{eqnarray}

 The metric for the MSM spacetime can be given by the Weyl-Papapetrou line
 element
 \begin{eqnarray}\label{eq:MNSpro}
  ds^2 & = & -f(dt-\omega d\phi)^2 \nonumber \\
       & + & f^{-1} \left[ e^{2\gamma} (d\rho^2+dz^2)+\rho^2 d\phi^2 \right]~~,
 \end{eqnarray}
 where all metric functions $f$,~$\omega$,~and~$\gamma$ are considered as
 functions of the prolate spheroidal coordinates $u,v$, while the coordinates
 $\rho,z$ are the corresponding cylindrical coordinates. The transformation
 between the two coordinate systems is
 \begin{equation} \label{eq:TrCylPS}
  \rho=\kappa \sqrt{(u^2-1)(1-v^2)},~~ z=\kappa u v~~,
 \end{equation}
 where
 \begin{equation} \label{eq:kappa}
  \kappa :=\sqrt{d+\delta}~~.
 \end{equation}

 The metric functions are
 \begin{eqnarray} \label{eq:MetrFunc}
  f &=& E/D~~, \nonumber \\
  e^{2\gamma} &=& E/16 \kappa^8 (u^2-v^2)^4~~,\\
  \omega &=& (v^2-1)F/E~~,\nonumber
 \end{eqnarray}
 where
 \begin{eqnarray} \label{eq:EDF}
  E &=& R^2+\lambda_1\lambda_2 S^2~~, \nonumber \\
  D &=& E+R P+\lambda_2 S T~~, \\
  F &=& R T-\lambda_1 S P~~,\nonumber
 \end{eqnarray}
 \begin{equation} \label{eq:lambda}
  \lambda_1=\kappa^2 (u^2-1),~~\lambda_2=v^2-1
 \end{equation}
 and
 \begin{eqnarray}\label{eq:PRST}
  P & := & 2 \{\kappa m u [(2 \kappa u+m)^2-2 v^2 (2 \delta +a b -b^2)\nonumber \\
    & - & a^2+b^2 -q^2]-2 \kappa^2 q^2 u^2-2 v^2 (4 \delta d-m^2 b^2)\}~~,
   \nonumber \\
  R & := & 4 [\kappa^2 (u^2-1)+\delta (1-v^2)]^2 \nonumber \\
    &  + & (a-b) [(a-b)(d-\delta)-m^2 b +q~\mu](1-v^2)^2~~,\nonumber \\
    & ~ &  \\
  S & := & -4 {(a-b)[\kappa^2(u^2-v^2)+2 \delta v^2]+v^2 (m^2 b-q~\mu)}~~,
  \nonumber \\
  T & := & 4(2 \kappa m b u+2 m^2 b-q~\mu)[\kappa^2 (u^2-1)+\delta (1-v^2)]
  \nonumber \\
    & + & (1-v^2)\{(a-b)(m^2 b^2 -4 \delta d)  \nonumber \\
    & - & (4 \kappa m u+2 m^2-q^2) [(a-b)(d-\delta)-m^2 b +q~\mu]\}. \nonumber
 \end{eqnarray}

 The MSM metric is of astrophysical interest as it has been suggested to be
 appropriate to model neutron stars \cite{MSM,Berti}. By investigating the
 dynamics of the geodesic orbits in the MSM spacetime, Dubeibe et al.
 \cite{Dubeibe07} and Han \cite{Han08} debated on the appearance and on the
 significance of chaos in the MSM model. Even if the main purpose of our article
 is to present a new integration scheme, in section \ref{sec-V} we contribute to
 the aforementioned debate.

 \subsection{Geodesic motion} \label{subsec-II-b}

 The equations of geodesic motion of a ``test'' particle of rest mass $m_0$ in a
 spacetime given by the metric $g_{\mu\nu}$ are produced by the Lagrangian
 function
 \begin{equation} \label{eq:LagDef}
  L=\frac{1}{2}~m_0~g_{\mu\nu}~ \dot{x}^{\mu} \dot{x}^{\nu}~~,
 \end{equation}
 where the dot denotes derivation with respect to proper time $\tau$. The
 Lagrangian (\ref{eq:LagDef}) has a constant value $L=-m_0/2$ along a
 geodesic orbit, due to the four-velocity
 $g_{\mu\nu}~\dot{x}^{\mu}\dot{x}^{\nu}=-1$ constraint.

 Since the spacetime is axisymmetric and stationary, the corresponding momenta
 \begin{equation}\label{eq:MomDef}
  p_\nu=\pd{L}{\dot{x}^\nu}
 \end{equation}
 are conserved. These are the specific energy
 \begin{equation} \label{eq:EnCon}
  E=-\frac1{m_0}\pd{L}{\dot{t}}~~,
 \end{equation}
 and the specific azimuthal component of the angular momentum
 \begin{equation} \label{eq:AnMomCon}
  L_z =\frac1{m_0}\frac{\partial L}{\partial \dot{\varphi}}~~.
 \end{equation}
 For brevity, we refer hereafter to these two integrals simply as the energy $E$
 and the angular momentum $L_z$. Due to these two integrals of motion, we can
 restrict our study to the meridian plane. We just have to re-express
 \eqref{eq:EnCon}, \eqref{eq:AnMomCon} in order to get $\dot{t}$ and $\dot{\phi}$
 as functions of $E$ and $L_z$, and then replace $\dot{t}$, $\dot{\phi}$ in
 the two remaining equations of motion. Then, from the original set of $4$
 coupled second order ordinary differential equations (ODEs), we arrive to a
 set of $2$ coupled ODEs.

 By simply applying the Legendre transform
 \begin{equation}\label{eq:LegTr}
  H=p_\mu \dot{x}^\mu-L
 \end{equation}
 to the Lagrangian function (\ref{eq:LagDef}), we get the Hamiltonian function
 \begin{equation}\label{eq:HamDef}
  H=\frac{1}{2m_0}g^{\mu\nu} p_\mu p_\nu~~,
 \end{equation}
 where the momenta $p_\mu$ are given by (\ref{eq:MomDef}) and
 $p^\nu=m_0 \dot{x}^\nu$. The Hamiltonian equations are
  \begin{align}\label{eqn-pdot_qdot}
   \dot{x}^\mu &= g^{\mu\nu}\pd{H}{p^\nu}~~, \nonumber \\
               &~                           \\
   \dot{p}^\mu &=-g^{\mu\nu}\pd{H}{x^\nu}~~. \nonumber
  \end{align}
 As the system is  autonomous
 ($\frac{d H}{d \tau}=\frac{\partial H}{\partial \tau}=0$), the Hamiltonian
 function is an integral of motion and equal to $H=-m_0/2$.

\section{Standard symplectic schemes}\label{sec-III}

 We have seen in section \ref{subsec-II-b} that the geodesic equations of motion
 can be described by a Hamiltonian formalism. Thus, symplectic schemes should be
 appropriate for integrating these equations.

 \subsection{Notation}\label{subsec not}
 In the rest of the paper we drop the Tensor Analysis covariant and
 contravariant symbolism by indices and exponents, and we adopt a symbolism more
 convenient for Numerical Analysis. By bold characters, we denote vectors, and,
 by arrows over characters, we denote vectors of vectors.

  We define
  \begin{align*}
    \mathbf y &:=\begin{pmatrix}\mathbf{p}\\ \mathbf{x}\end{pmatrix}~~, \\
    J &:=\begin{pmatrix} 0&I\\-I&0\end{pmatrix}~~, \\
    f(\mathbf{y}) &:=J^{-1}\grad{H(\mathbf{y})}~~,
  \end{align*}
   and write the given Hamiltonian system as
  \begin{align}\label{eqn-eom}
    \frac{d \mathbf{y}}{d \tau}=f(\mathbf{y})~~.
  \end{align}

 Further, $\Phi_h$ denotes an integration step of step size $h$. If the
 current position in phase space is $\mathbf{y}_n$, then $\Phi_h$
 propagates the system to the next position $\mathbf{y}_{n+1}$, i.e.
 \begin{align*}
   \mathbf{y}_{n+1}=\Phi_h(\mathbf{y}_n)~~.
 \end{align*}
 To simplify the Runge-Kutta schemes notation we define for an $s$-stage scheme
\begin{itemize}
 \item{}the vector
\begin{align}
  \vec y_n&:=\underbrace{\begin{pmatrix}\mathbf y_n&...&\mathbf y_n\end{pmatrix}^T}_{\text{s times}}~~,\nonumber
\end{align}
 \item{}the inner-stage values
\begin{align}
  \vec Y &:=\begin{pmatrix} \mathbf{Y}_1&...&\mathbf{Y}_s\end{pmatrix}^T~~,\nonumber
\end{align}
 \item{}the auxiliary variables
\begin{align}\label{not-Z_F_A}
  \vec Z &:=\begin{pmatrix} \mathbf{Z}_1&...&\mathbf{Z}_s\end{pmatrix}^T=\vec Y-\vec y_n~~,
\end{align}
  \item{}the function
\begin{align}
  F(\vec Y) &:=  \begin{pmatrix} f(\mathbf{Y}_1)&... &
  f(\mathbf{Y}_s)\end{pmatrix}^T~~,\nonumber
\end{align}
  \item{}and the coefficient matrix
\begin{align}
  A &:=\begin{pmatrix}a_{11}&\hdots&a_{1s}\\ \vdots&\ddots&\vdots
   \\ a_{s1}&\hdots&\ a_{ss}\end{pmatrix}~~.\nonumber
\end{align}
\end{itemize}

 The exclamation mark ``!'' over relation symbols denotes requirement.
 Furthermore, $\cal X$ denotes the phase space and $\vec{\cal X}$ denotes the
 space on which the $\vec Z$ are defined, i.e.
 $\vec{\cal X}:=\cal{X} \times \hdots \times \cal{X}$.

 As a norm for matrices we use the \textit{Frobenius norm}, which for a matrix
 $A$ is
  \begin{align*}
   \|A\|_{\text{Frob}}=\sqrt{\sum_{i,j}A_{ij}^2}~~.
  \end{align*}

 \subsection{The equations of motion from a numerical point of view}

 From a numerical point of view, it is easier to handle the Hamiltonian equations
 $\eqref{eqn-pdot_qdot}$ in the form $\eqref{eqn-eom}$. The Hamiltonian system
 $\eqref{eq:HamDef}$ is non-separable in the sense that the mass-matrix in the
 Hamiltonian depends on the positions $\mathbf{x}$ and, therefore, the momenta
 cannot be separated from the position coordinates. In the case of a
 non-integrable system, the set of coupled first order differential equations
 $\eqref{eqn-eom}$ can be very sensitive to changes in the arguments
 $\mathbf{y}$. This happens when we evolve  chaotic orbits by $\Phi_h$. In
 other words, if $Df(\mathbf{y})$ denotes the Jacobian of $f(\mathbf{y})$ at
 a point $\mathbf{y}$ of the phase-space $\cal X$, then there
 exists a subset ${\cal U} \subset {\cal X}$, so that
 \begin{align}\label{eqn-inequal_DF}
  \|Df(\mathbf{y})\|\gg 1 ~~~\text{$\forall~\mathbf{y}\in {\cal U}$}~~.
 \end{align}

 As the Hamiltonian system $\eqref{eq:HamDef}$ is not separable, a natural
 choice for a symmetric, symplectic integrator is a \textit{Gauss-collocation
 scheme} (see e.g., \cite{hairernorsettwanner}).
 \subsubsection{Collocation methods}

 Given an interval $[0,h]$, stages $0\le c_1<...<c_s\le 1$ and an initial-value
 problem
 \begin{align}
    &y(0) = y_0\nonumber\\
    &\dot{y} = f(\tau,y)~~,\nonumber
 \end{align}
 the polynomial $U(\tau)$ of degree $s$, satisfying
  \begin{align}
     &U(0) =y_0\nonumber\\
     &\dot{U} (c_ih) =f(c_ih,U(c_ih)),~~i=1,...,s~~,\nonumber
  \end{align}
 is called a \textit{collocation polynomial}. A \textit{collocation method}
 consists of finding such a $U(\tau)$ and then setting
  \begin{align}
    y(h)\approx U(h).
  \end{align}
 It can easily be shown (e.g., theorem 7.7. in \cite{hairernorsettwanner}) that
 a collocation method as described above is equivalent to an $s$-stage
 Runge-Kutta scheme
 \begin{align}\label{eqn-Zi}
  \mathbf{y}_{n+1} &= \mathbf{y}_n+h\sum_{i=1}^sb_if(\mathbf{Y}_i)~~,
  \nonumber\\
  \mathbf{Y}_i &=\mathbf y_n+ h\sum_{j=1}^sa_{ij}f(\mathbf{Y}_j)~~,
 \end{align}
 with coefficients
 \begin{align}
    a_{ij}=\int_0^{c_i}l_j(\tau)d\tau\nonumber\\
    b_i=\int_0^1l_i(\tau)d\tau\nonumber,
  \end{align}
 where $l_i(\tau)$ are the Lagrange-polynomials
 \begin{align*}
  l_i(\tau)=\frac{\prod_{i\neq j}(\tau-c_j)}{c_i-c_j}~~.
 \end{align*}
 If the inner stages $c_i$ are chosen as
 \begin{align*}
	c_i=\frac 1 2(1+\tilde c_i)~~,
 \end{align*}
 where $\tilde c_i$ are the roots of the Legendre-polynomial of degree $s$, one
 obtains a Gauss-collocation scheme, which can be shown to be of the highest
 possible order $\mathcal O(h^{2s})$. Furthermore, Gauss-collocation or
 \textit{Gauss-Runge-Kutta methods} are symmetric and symplectic, see
 \cite{hairerlubichwanner}, chapter V and VI.

 If we rewrite the system of implicit eqs. $\eqref{eqn-Zi}$ in matrix-vector
 notation, we get
 \begin{align}
  \begin{pmatrix}\mathbf Y_1\\ \vdots\\ \mathbf Y_s\end{pmatrix}=\begin{pmatrix}\mathbf y_n\\
  \vdots\\ \mathbf y_n\end{pmatrix}+h\begin{pmatrix}a_{11}I&...&a_{1s}I\\
  \vdots&\ddots&\vdots\\ a_{s1}I&...&a_{ss}I\end{pmatrix}
 \begin{pmatrix}f(\mathbf Y_1)\\ \vdots \\f(\mathbf Y_s)\end{pmatrix}\label{eqn-matrix-vector}~~,
 \end{align}
 with the $s\times s$-identity matrix $I$. With the help of the direct product
 of matrices $\otimes$ and the notations in $\ref{subsec not}$,
 eq. $\eqref{eqn-matrix-vector}$ becomes
 \begin{align}\label{eqn-Y}
  \vec Y=\vec y_n+h(A\otimes I)F(\vec Y)~~.
 \end{align}
 Following \cite{hairerwanner}, chapter IV.8, we use the auxiliary variables
 defined in eqs. $\eqref{not-Z_F_A}$, to get the system in a shorter form
 \begin{align}\label{eqn-Z}
  \vec Z=h(A\otimes I)F(\vec Z)~~,
 \end{align}
 where $F(\vec Z)$ is to be read as
 \begin{align}\label{eqn-FZ}
  F(\vec Z)=\begin{pmatrix}f(\mathbf Y_1)\\ \vdots \\f(\mathbf Y_s)\end{pmatrix}
  =\begin{pmatrix}f(\mathbf y_n+\mathbf Z_1)\\
   \vdots \\f(\mathbf y_n+\mathbf Z_s)\end{pmatrix}~~.
 \end{align}
 There are two possibilities to solve this implicit equation for the inner stage
 values $\vec Z$.
 \begin{itemize}
 \item The simplest and the most popular method is the
 \textit{Fixed-Point iteration}. One has to solve iteratively
  \begin{align}
   \vec Z^{k+1}=h(A\otimes I)F(\vec Z^k)~~.
  \end{align}
  The \textit{Banach fixed-point theorem} guarantees that the iteration converges
  towards the correct $\vec Z$, if it is a contraction. This requires
  \begin{align*}
    \|\vec Z^{k+2}-\vec Z^{k+1}\| \stackrel{!}{<}\|\vec Z^{k+1}-\vec Z^k\|~~.
  \end{align*}
  Because of
  \begin{align*}
   & \|\vec Z^{k+2}-\vec Z^{k+1}\|=\|h(A\otimes I)(F(\vec Z^{k+1})-F(\vec Z^k))\|\\
   & \le \max_{\vec{Z}\in \vec{\cal X}}~
   \|h(A\otimes I)DF(\vec Z)\|\|\vec Z^{k+1}-\vec Z^k\|~~,
   \end{align*}
  this implies
  \begin{align} \label{eqn-req_fpi}
   h\|(A\otimes I)DF(\vec{Z})\|\stackrel{!}{<}1,
     ~~~\forall~\vec{Z}\in \vec{\cal X}~~.
  \end{align}
 \item If the convergence of the fixed-point iteration is not ensured, we can
 look forward to get better accuracy by searching for the roots of
 \begin{align}
  \hat F(\vec Z):=\vec Z-h(A\otimes I)F(\vec Z)
 \end{align}
 via a \textit{modified Newton iteration}. In this case, we have to iterate
 \begin{align}
  \vec Z^{k+1} &= \vec Z^{k}+\Delta \vec Z^{k}\\
  \Delta\vec Z^k     &= -M^{-1}\hat F(\vec Z^k)~~\label{eqn-DeltaZ},
 \end{align}
 with
 \begin{align}
  M:=I-(A\otimes I)( I\otimes Df(\mathbf y_n))~~.
 \end{align}
 Again, convergence can only be ensured for a contraction. Therefore, we need
 \begin{align*}
  &\|\vec Z^{k+2}-\vec Z^{k+1}\|= \\
  &\|\vec Z^{k+1}-M^{-1}\hat F(\vec Z^{k+1})-\left(\vec Z^k-M^{-1}\hat
  F(\vec Z^k)\right)\|    \\
  & \le \max_{\vec{Z}\in \vec{\cal X}}~h\|I-M^{-1}D\hat F(\vec{Z})
  \|\|\vec Z^{k+1}-\vec Z^k\| \\
  & \stackrel{!}{<} \|\vec Z^{k+1}-\vec Z^k\|~~,
 \end{align*}
 or
 \begin{align} \label{eqn-req_ni}
  h\|I-M^{-1}D\hat F(\vec{Z})\|\stackrel{!}{<}1,
  ~~\forall ~\vec{Z}~\in~ \vec{\cal X}~~.
 \end{align}
 \end{itemize}
 As a single iteration step for the fixed-point iteration is computationally
 much cheaper than in the modified newton iteration case, one should prefer the
 first whenever it is applicable. In both cases, inequality
 $\eqref{eqn-inequal_DF}$ together with the respective requirements
 $\eqref{eqn-req_fpi}$ or $\eqref{eqn-req_ni}$ lead to severe restrictions on
 the possible step-size $h$. The \textit{ultima ratio} for the solution of
 eq. $\eqref{eqn-Z}$ is the use of a \textit{Newton-Raphson} method. There, $M$
 in eq. $\eqref{eqn-DeltaZ}$ is replaced by
 \begin{align*}
	M(\vec Z^k)=I-(A\otimes I)(DF(\vec Z^k))~~.
 \end{align*}
 This gives an iteration with quadratic convergence but is much more expensive
 because for every iteration, a Jacobian has to be calculated. Therefore,
 standard literature recommends to proceed without it whenever it is possible, e.g.
 \cite{hairerlubichwanner}, chapter VIII. We will see in section $\ref{sec-V}$,
 that the use of Newton-Raphson doesn't improve the performance of the integration
 scheme when it is applied in the case of geodesic orbits either.

 \subsection{Existing step-size controls} \label{subsec-existing}

 For efficiency's sake one would expect to use small time steps only for the
 subspace for which condition $\eqref{eqn-inequal_DF}$ holds. Therefore, a
 reliable step size control is necessary. A feasible step size  control has to
 preserve the geometric properties of the underlying integration scheme. As it
 was shown in \cite{stoffer88}, there cannot be an efficient symplectic
 integrator with adaptive step size. Thus, other structure preserving properties
 of the integrator, such us \textit{symmetry}, i.e.
 \begin{align}\label{eqn-symmetry1}
  \Phi_{-h}(y)=\Phi^{-1}(y)
 \end{align}
 and \textit{reversibility}, i.e.
 \begin{align}\label{eqn-reversibility1}
  \Phi^{-1}_h\circ \zeta =\zeta \circ \Phi_h
 \end{align}
 for the involution $\zeta(\mathbf{p},\mathbf{x})=(-\mathbf{p},\mathbf{x})$,
 have to be conserved. If we apply such an integrator to a symmetric, reversible
 differential equation, we would expect the same advantageous long term behavior
 as for a symplectic integrator, e.g., \cite{hairerlubichwanner}, chapter VIII.

 One of the most efficient step size controls has been proposed by Hairer and
 S\"oderlind \cite{soederlind}. They express the discrete variable time step
 sequence $\tau={\tau_1,...,\tau_n,\tau_{n+1},...}$ by means of a constant step
 size sequence $\epsilon={\epsilon_1,...,\epsilon_n,\epsilon_{n+1},...}$ and a
 scalar function $\sigma(\mathbf y)$ via $d\tau=\sigma(\mathbf y)d\epsilon$. For
 the dependence of $\mathbf{y}$ on $\epsilon$, they obtain

 \begin{align}\label{eqn-dy_dtau}
  \frac{d \mathbf{y}}{d \epsilon}=\frac{d \mathbf{y}}{d \tau} \d \tau \epsilon
  =f(\mathbf{y})\d \tau \epsilon~~.
 \end{align}

 Defining
 \begin{align}\label{eqn-w_hairer}
   w:=\frac 1{\sigma(\mathbf{y})}~~,
 \end{align}
 one finds
 \begin{align}\label{eqn-w_dot_hairer}
   \frac{d w}{d\epsilon}&=-\frac1{\sigma^2}\grad\sigma \frac{d \mathbf{y}}{d \epsilon} \\
     &=-\frac1{\sigma}\grad\sigma f(\mathbf{y})~~.
 \end{align}
 We now combine eqs. $\eqref{eqn-dy_dtau}$ and $\eqref{eqn-w_dot_hairer}$ and
 get the system
 \begin{align}\label{eqn-augmented_hairer}
  \begin{pmatrix}
   \frac{d\mathbf{y}}{d \epsilon}\\
   \frac{d w}{d \epsilon}
  \end{pmatrix}=
  \begin{pmatrix}
   \frac 1 w f(\mathbf{y})\\
   -\frac1{\sigma}\grad\sigma f(\mathbf{y})=:G(\mathbf{y})
  \end{pmatrix}~~.
 \end{align}

 Solving this system with the scheme
 \begin{align}
  w_{\text{$n+\frac 1 2$}} &=w_n+\epsilon \cdot G(\mathbf{y}_n) \label{eqn-w_hairer_num}\\
  \mathbf{y}_{\text{$n+1$}} &= \Phi_{\frac{\epsilon}{w_{n+\frac 1 2}}}(\mathbf{y}_n)\\
  w_{\text{$n+1$}} &=w_{\text{$n+\frac 1 2$}}+\epsilon \cdot G(\mathbf{y}_{n+1})~~,\\
  w_0 &=\frac1{\sigma(\mathbf{y}_0)}
 \end{align}
 we get a symmetric, reversible integrator, as it is shown in \cite{soederlind}.
 Nevertheless, this step-size algorithm has some drawbacks which can show up when
 integrating geodesic equations of motion of a non-integrable system.

 The above integration scheme is explicit in the variable $w$. Therefore, eq.
 $\eqref{eqn-w_hairer_num}$ represents a constraint on possible step size
 adapters $\sigma(\mathbf{y})$. If we don't want the underlying step size
 $\epsilon$ to be too small, $G(\sigma(\mathbf{y}))$ has to be bounded. This is
 similar to the case when we use explicit Runge-Kutta schemes for stiff
 differential equations, where, due to a very small range of stability, the
 step size has to be unfeasibly small.

 However, for the iterations to converge, we have to choose a step size small
 enough to get a contraction.
 \begin{itemize}

  \item If we want to solve equation $\eqref{eqn-Z}$ with a
  fixed-point-Iteration, the step size must be somehow proportional to
  the inverse of $\|Df(\mathbf{y})\|$, i.e. we have
  \begin{align}
   \sigma(\mathbf{y})\propto\frac1{\|Df(\mathbf{y})\|}~~,
  \end{align}
  and such a controller $\sigma$ yields
 \begin{align}
  G(\mathbf{y})=\frac{\sum_{ijk}\left(Df(\mathbf{y})\right)_{ij}
  \left(\pd{\left(Df(\mathbf{y})_{ij}\right)}
  {\mathbf{y}_k}\right)\left(f(\mathbf{y})\right)_k}
  {\|Df(\mathbf{y})\|^2}~~.
 \end{align}

 For regions of the phase space, where
 \begin{align}
  \|D^2f\cdot f\|\gg\|Df\|~~,
 \end{align}
 this would conflict with the requirement which demands $G(\mathbf{y})$ not to
 be too large. Numerical tests even show that in such cases the step size can
 become negative.

 \item In order to guarantee the convergence of the newton-iteration, i.e. to
 satisfy condition $\eqref{eqn-req_ni}$, it is required that
 \begin{align}
  & \sigma(\mathbf{y}) \propto  \nonumber \\
  & \frac1{\|I-M^{-1}D\hat F(\vec{Z})\|} = \nonumber \\
  & \frac1{\|I-(I-(A\otimes Df(\mathbf{y}_n)))^{-1}(I-(A\otimes I)DF(\vec{Z}))\|}~~.
  \nonumber \\
   & ~
 \end{align}
 It is by no means clear how one can consider the last expression as a function
 of $\mathbf y$. The only possibility would be to set
 $\vec Z=(\mathbf 0\hdots\mathbf 0)^T$ in the last expression. The calculation
 of the corresponding $G(\mathbf{y})$ would then lead to an expression with a
 factor
 \begin{align}
  \frac1{\|I-(I-(A\otimes Df(\mathbf{y}_n)))^{-1}(I-(A\otimes I)DF(\vec Z))\|^2}~~.\nonumber
 \end{align}
 Because of $DF(\vec Z)=I\otimes Df(y)$ for the $\vec Z$ chosen above, this
 factor would be
 \begin{align}\label{eqn-factor_w_hairer}
  \frac1{\|I-(I-(A\otimes Df(\mathbf{y}_n)))^{-1}(I-(A\otimes Df(\mathbf y)))\|^2}~~.
 \end{align}
 Thus, the evaluation of the function $G$ at the point $\mathbf{y}_n$ (which is
 necessary in eq. $\eqref{eqn-w_hairer_num}$) is impossible as we would have to
 divide by zero due to the factor $\eqref{eqn-factor_w_hairer}$. We see that
 it is not possible to use the modified newton iteration along with the
 presented step size control algorithm.

\end{itemize}

 \section{A new adaptive step size control of motion}
 \label{sec-IV}

 We now present a new integration scheme $\Phi_h$ with an adaptive step size
 $h(\epsilon,\mathbf{y})$. The integrator uses an $s$-stage Gauss-collocation
 method as the underlying integrator. The step size depends on an underlying
 constant step size $\epsilon$ and the actual state of the system $\mathbf{y}$.
 The integration scheme has the following properties:

 \begin{itemize}

 \item For the step size $h$ an equation similar to
  \begin{align}
   h(\epsilon,\mathbf{y})\propto \frac1{\|Df(\mathbf{y})\|}
  \end{align}
 holds, and, thus, condition $\eqref{eqn-req_fpi}$ is satisfied.

 \item The integrator is both symmetric and reversible. Hence, the theoretical
 results on the long term behavior of such integrators are applicable
 (e.g., \cite{hairerlubichwanner}, chapter XI).

 \end{itemize}

 The main idea of the step size control algorithm is to replace requirement
 $\eqref{eqn-req_fpi}$ by
 \begin{align}
  \|h\cdot Df(\mathbf{y})\|\stackrel !=\epsilon~~,
 \end{align}
 for some $\epsilon<1$. We slightly modify this formulation to
 \begin{align}\label{eqn-h}
  \|\frac h 2\left(Df(\mathbf{Y}_1)
  +Df(\mathbf{Y}_s)\right)\|\stackrel !=\epsilon~~,
 \end{align}
 with $\mathbf{Y}_{1(s)}$ being the inner-stage value at stage $1(s)$ of the
 Gauss-Runge-Kutta scheme. In the subsequent paragraphs we demonstrate that this
 modification gives the required structure preserving properties.

 For the integrator to be symmetric, we must make sure that when we propagate the
 system backwards in time (i.e. from $\mathbf{y}_{n+1}$ to $\mathbf{y}_n$) the step
 size has the same behavior as when propagating forward (i.e. from  $\mathbf{y}_n$
 to $\mathbf{y}_{n+1}$), i.e. we need
 \begin{align}
  h(\epsilon,\mathbf{y}_n)=-h(-\epsilon,\mathbf{y}_{n+1})~~.
 \end{align}

 This, together with the symmetry of the underlying integrator (i.e. the
 validity of property $\eqref{eqn-symmetry1}$ for constant steps $h$), ensures
 the symmetry of the whole integration scheme. Let $ \mathbf{\hat Y}_i$ denote
 the inner-stage values for the integration backwards in time. The symmetry of
 the Gauss-collocation method then results in
 \begin{align}\label{eqn-symmetry_Z}
  \mathbf{\hat Y}_i=\mathbf{Y}_{s+1-i}
  ~~\forall i=1,\hdots,s~~,
 \end{align}
 because the collocation polynomial for $\mathbf{y}_{n+1}=\Phi_h(\mathbf{y}_n)$,
 which is the interpolation polynomial through the points $(0,y_n)$ and
 $(c_i,\mathbf{Y}_{i})$, is the same as the
 collocation polynomial for $\mathbf{y}_{n}=\Phi_{-h}(\mathbf{y}_{n+1})$, i.e.
 the interpolation polynomial through the points $(1,y_{n+1})$ and
 $(c_i,\mathbf{\hat Y}_{i})$. Due to property
 $\eqref{eqn-symmetry_Z}$, we get
 \begin{align}
  h(-\epsilon,\mathbf{y}_{n+1})
  &= \frac{-\epsilon}{\|\frac 1 2[Df(\mathbf{\hat Y}_1)+Df(\mathbf{\hat Y}_s)]\|} \nonumber \\
  &= \frac{-\epsilon}{\|\frac 1 2[Df(\mathbf{Y}_s)+Df(\mathbf{Y}_1)]\|} \nonumber \\
  &= -h(\epsilon,\mathbf{y}_n)~~.
 \end{align}

 If we consider the new variable step size integrator
 $\Phi_h(\epsilon,\mathbf{y})$ as a constant step size integration scheme
 $\Psi_{\epsilon}$, then the reversibility condition $\eqref{eqn-reversibility1}$
 reads
 \begin{align}\label{eqn-reversibility2}
  \Psi^{-1}_\epsilon\circ \zeta =\zeta \circ \Psi_\epsilon~~.
 \end{align}

 If we have
 \begin{align}\label{eqn-h_reversible}
  h(-\epsilon,\zeta~\mathbf{y}_{n+1})=-h(\epsilon,\mathbf{y}_n)~~,
 \end{align}
 then the reversibility of the method is a direct consequence of the
 reversibility of the underlying integrator $\eqref{eqn-reversibility1}$. In
 order to prove condition $\eqref{eqn-h_reversible}$, we denote by
 $\mathbf{\hat Y}$ the inner stage values for the integration that starts at
 $\zeta~\mathbf{y}_{n+1}$. We then notice, that the reversibility of the
 Gauss-Runge-Kutta scheme implies
 \begin{align}
  \mathbf{\hat Y}_i=\zeta\mathbf{Y}_{s+1-i}~~.
 \end{align}
 This yields
 \begin{align}
  h(-\epsilon,\zeta \mathbf{y}_{n+1}) &=
  \frac{-\epsilon}{\|\frac 1 2[Df(\mathbf{\hat Y}_1)+
   Df(\mathbf{\hat Y}_s)]\|} \nonumber \\
  & =\frac{-\epsilon}{\|\frac 1 2\zeta^{-1} [Df(\mathbf{Y}_s)
    +Df(\mathbf{Y}_1)]\|} \nonumber \\
  & =-h(\epsilon,\mathbf{y}_n)~~,
 \end{align}
 where the last equality is a consequence of the orthogonality of the involution
 $\zeta$.

 Equation $\eqref{eqn-h}$, together with equation $\eqref{eqn-Z}$ for the
 inner-stage values $\mathbf{Z}_i$, leads to the system of equations
 \begin{align}\label{eqn-impl-IGEM}
  \begin{pmatrix} \vec Z\\h\end{pmatrix}=\begin{pmatrix}h(A\otimes I)F(\vec Z)
  \\\frac{\epsilon}{\|\frac 1 2[Df(\mathbf y_n+\mathbf{Z}_1)+Df(\mathbf y_n+\mathbf{Z}_s)]\|}\end{pmatrix}~~.
 \end{align}

 Inserting the second part into the first, we can easily apply a
 fixed-point-Iteration, which yields
 \begin{align} \label{eqn-FP-IGEM}
  \vec Z^{k+1}=\frac{\epsilon(A\otimes I)F(\vec Z^k)}
     {\frac 1 2\|Df(\mathbf y_n+\mathbf{Z}_1^k)+Df(\mathbf y_n+\mathbf{Z}_s^k)\|}~~.
 \end{align}

  \begin{figure*} [htp]
   \centering
   \includegraphics[width=0.8\textwidth]{FigOr}
    \caption{For $m=2.904$, $a=1.549$, $q=0$, $\mu=0$, $b=0.8$, $E=0.971$, and
   $L_z=9.3$, the left panel shows the projections on the $(\rho,z)$ plane of
   two regular orbits, and the right one projections of a chaotic orbit. In the
   left panel the initial conditions of the left orbit (magnified at the upper left
   embedded plot) are $\rho=1.7$, $z=\dot{\rho}=0$, while the initial conditions
   for the right regular orbit are $\rho=30.7$, $z=\dot{\rho}=0$. In the right
   panel, the initial conditions for the chaotic orbit are $\rho=0.7$,
   $z=\dot{\rho}=0$. A magnification of the left part of the chaotic orbit is shown
   in the embedded upper left plot.}
   \label{FigOr}
  \end{figure*}

 \section{Numerical experiments}\label{sec-V}

 In this section we test our Integrator for Geodesic Equations of Motion (IGEM)
 and compare it with several other integration schemes. These are

 \begin{itemize}

  \item A standard fifth order \textit{Cash-Karp-Runge-Kutta} scheme with
   constant step size $\epsilon$ (RK5con) as proposed by Numerical Recipes
   \cite{nr}.

  \item The Cash-Karp-Runge-Kutta scheme (RK5var) with initial step size
  $\epsilon$ and a step size algorithm as follows:

   \begin{itemize}

    \item Starting with $\mathbf{y}_n$ calculate
            $\mathbf{y}_{n+1}=\Phi_h(\mathbf{y}_n)$

    \item Calculate the Hamiltonian's relative error
  $\delta H=\abs{\frac{H(\mathbf{y}_{n+1})-H(\mathbf{y}_{n})}{H(\mathbf{y}_{n})}}$,
   between two consecutive integration points $\mathbf{y}_{n},~\mathbf{y}_{n+1}$.

    \item
     \begin{itemize}
      \item If $\delta H$ is greater than a minimum tolerance value
      $tol_1=10^{-12}$, then repeat the integration step with $h_{new}=\frac h 2$.

      \item If $\delta H$ is smaller than a second tolerance value
      $tol_2=10^{-14}$, then double the step size $h$ before calculating the
      next step.

      \item If $tol_2\le\delta H\le tol_1$, then calculate the next step.
     \end{itemize}
   \end{itemize}

  \item A standard symplectic integration scheme used in classical celestial
  mechanics given by the Gauss-collocation as underlying integrator together
  with the step size control explained in section $\ref{subsec-existing}$ with
  $\sigma(\mathbf{y})=\frac1{\|f(\mathbf{y})\|}$ (CCM).

  \item Again the standard scheme from classical celestial mechanics with
  $\sigma(\mathbf{y})$ replaced by
  $\sigma(\mathbf{y})=\frac1{\|Df(\mathbf{y})\|}$ (CCM2).

  \end{itemize}

 In all implicit schemes, we use a fixed-point iteration which we implement as
 proposed by \cite{hairer2}:
 In order to increase the accuracy, one iterates until either the difference
 between two subsequent $\vec Z$ is numerically $0$ or fluctuations due to
 round-off errors start to occur. In other words, we use the stopping criterion
 \begin{align}
  \|\vec Z^{k+1}-\vec Z^k\|=0~\text{or}~\|\vec Z^{k+1}-\vec Z^k\|> \|\vec Z^k-\vec Z^{k-1}\|
  \nonumber
 \end{align}

 As test case for the algorithms, we take the system described in section
 $\ref{sec-II}$. In particular, we use the set of MSM parameters and initial
 condition discussed in \cite{Han08}. Accordingly, we choose a central object
 with mass $m=2.904$, spin $a=1.549$, and charge $q=0$, while the
 real parameters $\mu$ and $b$ correlated with the magnetic dipole and the
 mass-quadrapole deformation are set to be $0$ and $0.8$ respectively.
 Regarding the initial conditions, the energy and the angular momentum are set
 to be $E=0.971$ and $L_z=9.3$, respectively, while we fix $\dot{\rho}=z=0$.
 The integrators are tested for 3 different values of $\rho$, which correspond
 to 3 different cases of orbits: two regular ones (left panel of Fig.
 \ref{FigOr}), each covering quasiperiodically a corresponding torus in the
 phase space, and one chaotic (right panel of Fig. \ref{FigOr}), which covers in
 an irregular manner all the available phase space. We propagate the
 system in time until the proper time reaches $\tau_f=500000$ and compare the
 cpu times $T_{\text{calc}}$  the calculation lasted. As a measure of our code's
 accuracy, we track the overall relative error
 $\Delta H=\abs{\frac{H(\mathbf{y}_{n})-H}{H}}$, where $H$ is the theoretical
 value, while $H(\mathbf{y}_{n})$ is the computed Hamiltonian value at the
 integration point $\mathbf{y}_{n}$. We abort the simulations if
 $\Delta H>10^{-6}$.

  \begin{figure*} [htp]
   \centering
   \includegraphics[width=0.8\textwidth]{energies307}
    \caption{The relative error $\Delta H$ as function of the proper time $\tau$
    in semi-logarithmic scale for the regular orbit with $\rho=30.7$.}
   \label{fig-307-energy}
  \end{figure*}

 \subsection{Case with $\rho=30.7$}

 The initial condition $\rho=30.7$ corresponds to a regular orbit (right orbit in
 the left panel of Fig. \ref{FigOr}). Fig. $\ref{fig-307-energy}$ shows that the
 error is much larger for the explicit schemes, although we have chosen a much
 smaller step size $\epsilon$. In table $\ref{tab-307}$ we see that not only are
 the explicit integrators less accurate, but they also require much more
 computation time. Moreover, if we increase the step sizes $\epsilon$ of the
 explicit schemes, the energy drift becomes much worse. On the other hand, CCM
 is the fastest algorithm in this case because the corresponding regular orbit does
 not posses the special feature discussed in section \ref{subsec-existing} and
 expressed by condition $\eqref{eqn-inequal_DF}$. Overall, the explicit
 integrators show a linear drift whereas the symmetric integrators preserve the
 Hamiltonian's constant value.

  \begin{table}[htp]
   \centering
   \begin{tabular}{ccc}
    \toprule
     Integrator & $\epsilon$ & $T_{\text{calc}}[s]$\\
    \midrule
     RK5con & $0.01$ & $222.9$\\
     RK5var & $0.01$ & $160.6$\\
     CCM    & $1.0$ & $17.3$\\
     CCM2   & $1.0$ & $46.1$\\
     IGEM   & $1.0$ & $41.3$\\
    \bottomrule
  \end{tabular}

 \normalsize

  \caption{The cpu calculation times for the proper time interval
 $\tau \in [0,500000]$ for different integration schemes in the case of the
 regular orbit with $\rho=30.7$.
  }
  \label{tab-307}
 \end{table}

 \subsection{Case with $\rho=1.7$}

 Even though the initial condition $\rho=1.7$ corresponds to a regular orbit
 (embedded plot of the left panel of Fig. \ref{FigOr}), the system of
 differential equations behaves in a more ill-mannered way than in the
 $\rho=30.7$ regular case. This is shown by tracking the step sizes of the
 IGEM-integrator along the propagation (Fig. $\ref{fig-17-h(IGEM)}$), as we have
 $h\propto\frac 1{\|Df(y(\tau))\|}$. $\|Df(y(\tau))\|$ is varying fast with time.

 \begin{figure} [htp]
  \centering
  \includegraphics[width=0.4\textwidth]{hIGEM}
  \caption{ The step size $h$ of the IGEM as function of proper time
  $\tau\in[0,500]$ for the regular orbit with $\rho=1.7$.}
  \label{fig-17-h(IGEM)}
 \end{figure}
 \begin{figure} [htp]
  \centering
  \includegraphics[width=0.4\textwidth]{energies17}
  \caption{The relative error $\Delta H$ as function of the proper time $\tau$
   in normal scale for the regular orbit with $\rho=1.7$.}
  \label{fig-17-energy}
 \end{figure}
 \begin{table}[htp]
  \centering

 \begin{tabular}{ccc}
  \toprule
   Integrator & $\epsilon$ & $T_{\text{calc}}[s]$\\
  \midrule
   RK5con & $10^{-4}$ & $23813.2$\\
   RK5var & $10^{-4}$ & $19207.7$\\
   CCM & $0.1$ & $266.9$ \\
   CCM2 & $0.1$ & $2100.9$ \\
   IGEM & $0.1$ & $2202.1$\\
  \bottomrule
 \end{tabular}

 \normalsize
 \caption{The calculation times $T_{\text{calc}}$ for the proper time interval
  $\tau \in[0,500000]$ for the different integrators in the case of the
 regular orbit with $\rho=1.7$.}
  \label{tab-17}
 \end{table}

 When we tested the RK5var scheme, it needed $1$ week of calculation time only
 to arrive at $\tau=758$. Therefore, we relaxed the restrictions on the
 acceptable relative errors per step and set $tol_1=10^{-9}$ and
 $tol_2=10^{-11}$. Now, if we compare the five integrations schemes in table
 $\ref{tab-17}$ and Fig. $\ref{fig-17-energy}$ we notice that the symplectic
 standard schemes from classical celestial mechanics -although fast- show very
 bad conservation properties for the constant of motion. The relative error is
 even worse than for the constant step size RK5con which shows a linear drift.
 On the other hand, IGEM shows the best long term behavior for affordable
 computational cost.

 \subsection{Case with $\rho=0.7$}

 \begin{figure} [htp]
  \centering
  \includegraphics[width=0.4\textwidth]{spectrum07II}
  \caption{The spectrum of frequencies
  $\abs{f(\omega)}=\abs{\int{e^{i\omega\tau}}\rho(\tau)d\tau}$,
  shows the characteristic noise of chaotic orbits. A magnification of the Fourier
 spectrum is shown in the embedded plot.}
  \label{fig-07-spec}
 \end{figure}

 \begin{table*}[htp]
  \centering

  \begin{tabular}{cccc}
  \toprule
   Integrator & $\epsilon$ & result & $T_{\text{calc}}[s]$ until abortion\\
  \midrule
   RK5con & $10^{-6}$ & aborted after propagation time of $\tau=11254.0$ because
   $\Delta H>10^{-6}$ & $40739.7$ \\
   RK5var & $10^{-4}$ & aborted after propagation time of $\tau=1087.3$ because
   $\Delta H>10^{-6}$ & $149469.9$\\
   CCM & $0.01$ & aborted after propagation time of $\tau=9559.9$ because
   $h(l,\mathbf{y})<0$ & $34.2$\\
   CCM2 & $0.01$ & aborted after propagation time of $\tau=4180.1$ because
   $\Delta H>10^{-6}$ & $1353.9$\\
   IGEM & $0.1$ & no abortion in $\tau=[0,50000]$, $T_{\text{calc}}=995.3s$.
  & $-$\\
  \bottomrule
  \end{tabular}
  \normalsize
  \caption{Results of the runs of the chaotic orbit with $\rho=0.7$ for
 different integration schemes.}
 \label{tab-07}
 \end{table*}

 The initial condition $\rho=0.7$ corresponds to an orbit evolving in a strongly
 chaotic region (right panel of Fig. \ref{FigOr}). The appearing randomness of the
 motion in a strongly chaotic region produces a Fourier spectrum which corresponds
 to noise\footnote{In a Fourier spectrum of a weakly chaotic orbit we can ``detect''
 frequencies from neighboring regular orbits embedded in noise. However, such
 frequencies disappear as the chaos gets stronger, and finally we get only noise,
 see e.g. \cite{Contop02} and references therein.} (Fig. $\ref{fig-07-spec}$). In
 this case the calculations become much more computationally expensive than in the
 two other cases. Therefore, we restrict the propagation interval to
 $\tau\in[0,50000]$. The runs using the well established integrators were aborted
 very soon due to low accuracy, or in the case of CCM because the variable step
 size indeed became negative (table $\ref{tab-07}$)!

 The only integrator to pass this test was IGEM. For IGEM, we plot the
 relative error $\Delta H$ in Fig. $\ref{fig-07-energy-IGEM}$ and observe that
 some peaks appear. This happens because, during each step, we have to solve the
 iteration $\eqref{eqn-FP-IGEM}$. For this to be a contraction, we need
 \begin{align}
  \left\|D\left(\frac{\epsilon(A\otimes I)F(\vec Z)}
  {\frac 1 2\|Df(\mathbf{Z}_1)+Df(\mathbf{Z}_s)\|}\right)\right\|\stackrel ! <1~~,
 \end{align}
  with
 \begin{align}
  D\left(\hdots\right)_{ij}
  &=\frac{\left(\epsilon(A\otimes I)DF(\vec Z)\right)_{ij}}
  {\frac 1 2\|Df(\mathbf{Z}_1)+Df(\mathbf{Z}_s)\|} \nonumber \\
  &-\frac{\left(\epsilon(A\otimes I)F(\vec Z)\right)_{i}}{\frac 1 2\
  |Df(\mathbf{Z}_1)+Df(
  \mathbf{Z}_s)\|^3}  \nonumber \\
  &\times \sum_{kl}\left(Df(\mathbf{Z}_1)_{kl}+Df(\mathbf{Z}_s)_{kl}\right)
  \nonumber \\
  &\times \pd{\left(Df(\mathbf{Z}_1)_{kl}
  +Df(\mathbf{Z}_s)_{kl}\right)}{\vec Z_j} ~~.
 \end{align}

 \begin{figure} [htp]
  \centering
  \includegraphics[width=0.4\textwidth]{deltaIGEM}
  \caption{The relative error $\Delta H$ for the IGEM method as a function of
   proper time $\tau$, semi-logarithmic scale, for the chaotic orbit with
   $\rho=0.7$.}
  \label{fig-07-energy-IGEM}
 \end{figure}

 \begin{figure} [htp]
  \centering
  \includegraphics[width=0.4\textwidth]{velocities07}
  \caption{The 'velocity' $V(\tau)^2$ as function of proper time $\tau$, for the
  chaotic orbit with $\rho=0.7$.}
  \label{fig-07-velocity}
 \end{figure}

 The norm of the first term is smaller than $1$ by construction of the algorithm
 but the second term can be large. This happens at points, where $f(\mathbf y)$
 becomes almost singular. To see when such a situation occurs, we plot the
 `velocity` $V^2(\tau)=\dot z^2+\dot\rho^2$ in Fig. $\ref{fig-07-velocity}$.
 We see that at every instant of time where the Hamiltonian deviates from its
 otherwise constant value, there is a peak in the velocity. To analyze this further,
 we notice that the first peak in $V^2$ occurs at $\tau=2296.39$. We thus analyze
 the trajectory for $\tau\in[2286,2306]$. To do so, we compare the trajectory
 obtained via IGEM with $\epsilon=10^{-1}$ and once for a more 'accurate'
 calculation. In order to do the latter, we propagate the system with IGEM until
 $\tau=2296$. Then, we reduce the underlying step size $\epsilon$ to
 $\epsilon=10^{-4}$ and propagate the system until $\tau=2297$. Afterwards, we
 relax again to $\epsilon=0.1$ and calculate the configuration for
 $\tau\in[2297,2306]$. We observe that the jump in the Hamiltonian coincides with
 an inflexion point of the trajectory (i.e. the test-particle behaves like a ball
 that is thrown upon a wall), where $f(\mathbf y)$ becomes almost singular. At this
 point, IGEM with $\epsilon=10^{-1}$ cannot dissolve the trajectory with the same
 accuracy as for the other points. But, we see clearly that after the inflexion,
 the trajectory obtained via IGEM tends to follow the more accurate calculation.
 Namely, by calculating the quantity
 $   \mathrm{Dif}^2 =\frac{(\rho_{\text{IGEM}}-\rho_{\text{ac}})^2}{\rho_{\text{ac}}^2}
  +\frac{(\dot{\rho}_{\text{IGEM}}-\dot{\rho}_{\text{ac}})^2}{\dot{\rho}_{\text{ac}}^2}
 +\frac{(z_{\text{IGEM}}-z_{\text{ac}})^2}{z_{\text{ac}}^2}+\frac{(\dot{z}_{\text{IGEM}}-
   \dot{z}_{\text{ac}})^2}{\dot{z}_{\text{ac}}^2}~$,
 we observe, that $\mathrm{Dif}$ is larger than $10^{-2}$ for times near the
 inflexion but decreases to below $10^{-6}$ as $\tau$ increases until $\tau=2306$.
 This observation is confirmed by Fig. $\ref{fig-07-energy-IGEM}$,  where we see
 that after the inflexion, the Hamiltonian recovers its previous value. Hence, the
 relative error $\Delta H$ jump is negligible compared to the time scale of the
 studied phenomenon.

 To improve the performance of IGEM even more, one can hope to erase the
 convergence problems at the singular points by using the Newton-Raphson-method to
 solve the implicit eq. (\ref{eqn-impl-IGEM}). Doing so, we notice, however,
 that the calculation time increases considerably to $T_{\text{calc}}=7521.0~s$.
 This happens, because for each iteration, the costly derivative of the rhs of eq.
 (\ref{eqn-FP-IGEM}) has to be calculated, while the average number of iterations
 per step only decreases from $3.94$ for the Fixed-Point iteration to $3.57$ for
 Newton-Raphson iteration. Therefore, even though in theory the Newton-Raphson has
 a quadratic convergence rate, it cannot catch up with the Fixed-Point iteration.

 Our last numerical analysis argument is that IGEM is efficient and accurate in
 calculating Poincar\'e-sections because of the collocation-property. It is known
 from  Numerical Analysis that the solution at point $\mathbf y_{n+1}$ calculated
 with an $s$-stage Gauss-collocation method coincides with $U(h)$, where $ U(\tau)$
 is the interpolation polynomial through the points $(0,\mathbf y_n)$ and
 $(c_1,\mathbf y_n+\mathbf Z_1)...(c_s,\mathbf y_n+\mathbf Z_s)$.
 Numerical Analysis further shows that the interpolation-polynomial $ U(\tau)$
 stays $\mathcal O(h^s)$ close to the real solution within the whole interval
 between $\mathbf y_n$ and $\mathbf y_{n+1}$. Thus, if we search for the root of
 $ U(\tau)$ (or rather the root of the component of the vector-valued $U(\tau)$
 which corresponds to $z$) with a fast bisection method, we get the location of the
 Poincar\'{e} section $\mathcal O(h^s)$ close to the real section. We use this
 advantage of the IGEM to enter the discussion about the appearance of chaos in the
 MSM model.

 \begin{figure} [htp]
  \centering
  \includegraphics[width=0.4\textwidth]{FigSS01}
  \caption{The Poincar\'{e} section $z=0$ for $m=2.904$, $a=1.549$, $q=0$, $\mu=0$
  and $b=0.8$, where $E=0.971$ and $L_z=9.3$. The embedded plot is a detail of the
  Poincar\'{e} section.}
 \label{figSS}
 \end{figure}

 In \cite{Han08}, Han showed that chaos appears in both oblate and prolate
 deformations (with respect to the Kerr metric) of the MSM metric. However, he
 states that ``it is impossible to present chaos regions for the oblate
 deformation massive bodies, because those regions found in \cite{Dubeibe07} are
 actually inside a neutron star''. Even though Han discusses the oblate case
 $m=2.904$, $a=1.549$, $q=0$, $\mu=0$, $b=0.8$, $E=0.971$, $L_z=9.3$, he doesn't
 show the respective Poincar\'{e} section. This Poincar\'{e} section is shown in
 Fig. \ref{figSS} and it is obvious that the main island of stability, which
 lies far from the central object, is embedded in a pronounced chaotic layer. Thus,
 Han's statement should rather be that ``for realistic astrophysical objects that
 are oblate'', there are chaotic orbits of ``test particles around single oblate
 deformation neutron stars described by'' the MSM model. However, orbits starting
 from the external chaotic layer follow trajectories like those shown in the right
 panel of Fig. \ref{FigOr}, thus these chaotic orbits in the relativistic case
 would probably plunge to the central compact object.

 \section{Conclusion}\label{sec-VI}

 Applying several well established standard integration schemes to the system of
 differential equations describing geodesic motion in an example of a non-Kerr
 spacetime background (i.e. MSM \cite{MSM}), we showed that these integrators
 cannot guarantee satisfactory long term behavior for orbits in these
 non-integrable systems. Therefore, we introduce a new integration scheme
 appropriate for evolving orbits of such systems.

 The new integration scheme effectively conserves the integrals of motion that
 such Hamiltonian systems possess, and it is well-behaved in the case of long
 term evolution of strongly chaotic orbits. Thus, the new integration scheme is
 well-suited for studying geodesic orbits in the case of spacetime backgrounds
 which are non-integrable perturbations of the Kerr spacetime.

 Moreover, we show that chaos can appear in oblate deformations of the MSM metric
 modeling the exterior of a neutron star.

\begin{acknowledgments}
 We would like to thank B.~Bruegmann and Ch.~Lubich for useful discussions and
 suggestions. This work was supported by the DFG grant SFB/Transregio 7.
\end{acknowledgments}


\section{Introduction}

An important source of gravitational waves for space borne detectors are E(I)MRIs, extreme (intermediate) mass ratio inspirals. Such sources will allow us to test directly the Kerr hypothesis, and allow us to map the spacetime surrounding black holes. Moreover, the detection of E(I)MRIs will allow us to determine the mechanisms that shape stellar dynamics in galactic nuclei with unprecedented precision (\cite{Amaro:2012} and references cited therein).

The orbits of E(I)MRIs are typically highly relativistic, and exhibit exciting phenomena, e.g.~extreme periastron and orbital plane precessions. Because the orbital evolution time scale (``radiation reaction time scale'') is much longer than the orbital period(s), over short time scales the orbit is approximately geodesic, yet on long time scales it deviates strongly from geodesic motion of the background. Instead, the smaller objects moves along a geodesic of a perturbed spacetime. Equivalently, one may construe the orbit as an accelerated, non-geodesic motion in the spacetime of the unperturbed central object, where the acceleration is caused by the self force  of the smaller object (\cite{Poisson-Pound-Vega}). We undertake here the latter approach.

We consider here a point particle of mass $\mu$ that moves in a quasi circular orbit (i.e., an orbit that would be circular but for the radiation reaction effects) in the space-time of a Schwarzschild black hole of mass $M$, under the assumption that $\mu\ll M$. Specifically, we take here $\mu=10^{-2}\,M$, and the orbit starts at $r_0=10\,M$ down to close to the ISCO at $r=6\,M$. The orbital evolution is driven by the particle's self force. In practice we use the self force found for circular orbits in the Lorenz gauge and interpolate (using the least number of necessary terms) the tabulated values to obtain smooth functions for the self force to the accuracy given in \cite{barack-sago:2007}. Specifically, we expand the components of the self force as follows:
\begin{eqnarray*}
f^t_{r\le 8\,M}&=&-
\frac{1}{\sqrt{1-\frac{3M}{r}}\left(1-\frac{2M}{r}\right)} \; \bigg(\frac{M}{r}\bigg)^5\; \bigg(\frac{\mu}{M}\bigg)^2\;
\Bigg[a^-_0+a^-_1\,\frac{M}{r}+a^-_2\,\bigg(\frac{M}{r}\bigg)^2+a^-_3\,\bigg(\frac{M}{r}\bigg)^3 + \cdots \Bigg]
\nonumber\\
f^t_{r\ge 8\,M}&=&-
\frac{32}{5}\,\frac{1}{\sqrt{1-\frac{3M}{r}}\left(1-\frac{2M}{r}\right)} \;\left(\frac{M}{r}\right)^5\;\left(\frac{\mu}{M}\right)^2\;
\left[ {\rm PN}_{5.5}+\left(a^+_6+a^+_{6L}\,\ln\frac{M}{r}\right)\left(\frac{M}{r}\right)^6 + \cdots \right]
\nonumber\\
f^r_{r\le 8\,M}&=&
\bigg(1-\frac{2M}{r}\bigg)\; \bigg(\frac{M}{r}\bigg)^2
\bigg(\frac{\mu}{M}\bigg)^2\;
\Bigg[b^-_0+b^-_1\,\bigg(1-\frac{6M}{r}\bigg)+b^-_2\,\bigg(1-\frac{6M}{r}\bigg)^2+b^-_3\,\bigg(1-\frac{6M}{r}\bigg)^3 + \cdots\Bigg]\nonumber\\
f^r_{r\ge 8\,M}&=& \bigg(\frac{M}{r}\bigg)^2
\bigg(\frac{\mu}{M}\bigg)^2\;
\Bigg[b^+_0+b^+_1\,\frac{M}{r}+b^+_2\,\bigg(\frac{M}{r}\bigg)^2
+b^+_3\,\bigg(\frac{M}{r}\bigg)^3 + \cdots\Bigg]\nonumber
\end{eqnarray*}
where ${\rm PN}_{5.5}$ stands for the standard $\frac{11}{2}$--post--Newtonian expression (converting Eq.~(3.1) in \cite{tanaka:1996} from luminosity to $f^t$). Fitting the free parameters, we find the values appearing in Table 1.
\begin{table}[htdp]
\caption{The fit parameters for the self force. These parameters reproduce the accuracy of \cite{barack-sago:2007} to all significant figures for all data points. Our results for $a^+_{6,6L}$ are very crude predictions for the corresponding PN parameters, as our fit ignores all higher--order terms.}
\begin{center}
\begin{tabular}{|c|c||c|c||c|c||c|c||}\hline\hline
$a^-_0$ & 4.57583 & $a^+_6$ & 331.525 & $b^-_0$ & 1.32120 & $b^+_0$ & 1.999991\\
\hline
$a^-_1$ & 31.8117 & $a^+_{6L}$ & -2081.57 & $b^-_1$ & 1.2391 & $b^+_1$ & -6.9969\\
\hline
$a^-_2$ & -267.250 & & & $b^-_2$ & -1.297 & $b^+_2$ & 6.29\\
\hline
$a^-_3$ & 1049.27 & & & $b^-_3$ & 1.07 & $b^+_3$ & -24.6\\
\hline\hline
\end{tabular}
\end{center}
\label{default}
\end{table}%
Finally, the last remaining non--zero component of the self force ($u^{\phi}$) is found from the condition that the self force must satisfy $f^{\rm SF}_{\mu}\,u^{\mu}=0$, where the 4--velocity is taken for the circular geodesic orbit for which the self force was calculated.

\begin{figure}[!h]
\includegraphics[width=10.0cm,angle=0]{conv_c.eps}
\caption{Convergence tests for the two codes. Top: the 4-position (A) and 4-velocity (B) for the osculating code as functions of the azimuthal angle $\phi$. Bottom: the 4-position (C) and 4-velocity (D) for the direct code as a function of the time $t$.}
\label{fig1}
\end{figure}

\begin{figure}[!h]
\includegraphics[width=10.0cm,angle=0]{orbit4_c.eps}
\caption{The osculating--code variables $\alpha$ and $\beta$ as functions of the time $t$ (upper two panels), and the effective eccentricity $\epsilon$ of the orbit (lower panel) as a function of $t$ in the osculating case.}
\label{fig2}
\end{figure}

We evolve the orbit using three different codes. First, we use the energy balance approach to evolve the orbit. This approach ignores conservative effects on the orbital evolution and correspondingly on the waveforms, and can be done without the self force: fluxes of otherwise conserved qualities to infinity and down the event horizon are used to update the particle's constants of motion. As expected, we show in \cite{burko:2012} that the orbit and the resulting waveforms are identical in the energy balance approach and in our self force calculations, when the conservative piece of the self force is turned off (i.e., when we take $b_i^A=0$.) We used two independent codes for the computation of the orbital evolution with the self force. First, we used the method of osculating geodesics (specifically eqs. (43)--(47) in \cite{pound:2008}). We also integrate the equations of motion directly, i.e., solve directly for $u^{\mu}$ from the equation of motion $u^{\beta}\,\nabla_{\beta}u^{\alpha}=\mu^{-1}\,f^{\alpha}_{\rm SF}$ (``Newton's second law,'' with covariant differentiation compatible with the background metric) and integrate its solution to find $x^{\mu}$. Both codes are numerically stable and convergent. Specifically, the osculating code converges with 5$^{\rm th}$ order, and the direct code converges with $4^{\rm th}$ order (Fig.~\ref{fig1}).

The integration using the method of osculating geodesics cannot keep the value of the eccentricity as precisely zero. As both variables $\alpha$ and $\beta$ (see \cite{pound:2008} for definitions) are dynamical, the eccentricity $\epsilon$ must evolve along the orbit too. This behavior is shown in Fig.~\ref{fig2}. Interestingly, the inclusion of the conservative piece of the self force amplifies the resulting effective eccentricity.

\begin{figure}[!h]
\includegraphics[width=10.0cm,angle=0]{four_pos_new_c.eps}
\caption{The Orbit. The 4--position for the three orbital evolution codes: energy balance (dotted), direct evolution (dashed), and osculating code (solid).}
\label{fig3}
\end{figure}

\begin{figure}[!h]
\includegraphics[width=10.0cm,angle=0]{four_v_new_c.eps}
\caption{The Orbit. The 4--velocity for the three orbital evolution codes: energy balance (dotted), direct evolution (dashed), and osculating code (solid).}
\label{fig3b}
\end{figure}

\begin{figure}[!h]
\includegraphics[width=10.0cm,angle=0]{new_orbit_c.eps}
\caption{The Orbit. The shape of the orbit for the three orbital evolution codes: energy balance (dotted), direct evolution (dashed), and osculating code (solid).}
\label{fig3c}
\end{figure}

\begin{figure}[!h]
\includegraphics[width=10.0cm,angle=0]{new_orbit_section_c.eps}
\caption{The Orbit. A small portion of the orbit is magnified to show detail.}
\label{fig3d}
\end{figure}

The orbit is displayed in Figs.~\ref{fig3} and \ref{fig3b} and Figs.~\ref{fig3c} \ref{fig3d} for the three codes. Notably, the two independent self force codes reproduce the orbit to high level of agreement, with a difference much smaller than the difference between either and the orbit generated in the energy balance approach. This difference is attributed to the effect of the conservative piece of the self force.

The orbit, of course, is a gauge dependent quantity. Indeed, the position vector changes trivially under gauge transformations, $x^{\mu}\to x^{\mu}+\xi^{\mu}$. We can, however, create gauge invariant quantities with a specific gauge choice (in our case, the Lorenz gauge), and then those quantities -- by virtue of their gauge invariance -- are guaranteed to remain unchanged in any other gauge. Two independent gauge invariant quantities are $u^t$ (``gravitational redshift") and the angular frequency $\Omega$ (\cite{detweiler:2008}). In Fig.~\ref{fig4} we plot $u^t$ as a function of $\Omega$ with and without the conservative piece of the self force. Notably, to the accuracy of our numerical computation the two curves overlap. That is, we find that $u^t$ as a function of $\Omega$ is insensitive to the conservative piece of the self force. This conclusion implies that when an actual data stream is used and this gauge invariant figure is plotted, one may use a simplified radiation--reaction scheme, that does not include the conservative effects in its analysis.

There is, however, an aspect of the gauge invariant figure that is sensitive to the conservative effects, specifically the speed with which the data point moves along the curve. The way the conservative effects are manifested in the gauge invariant plot is not is the shape of the curve, but in the time it takes the signal to move along it.

\begin{figure}[!h]
\includegraphics[width=10.0cm,angle=0]{gauge_inv_c.eps}
\caption{The gauge invariant figure of $u^t$ as a function of the angular frequency $\Omega$. The curve without (dashed, $\circ$) and with the conservative effects  (solid, $\square$) are shown, together with equal-$t$ spacing marks starting at $t=1,500\,M$ in increments of $1,000\,M$. The two curves are indistinguishable to the numerical accuracy of our computation
(using the osculating code in both cases).}
\label{fig4}
\end{figure}

\begin{figure}[!h]
\includegraphics[width=17.0cm,angle=0]{full_wave_c.eps}
\caption{The real part of the waveform $\psi_4$ when the conservative piece of the self force is included (solid for the osculating orbit method, and dashed for the direct integration method), and using the energy balance approach (dotted). The three waveforms are in phase at the beginning of the waveforms, but the former two dephase with time from the latter (and to a much smaller extent from each other).}
\label{fig5}
\end{figure}

After we obtain the orbit we used it in order to evolve the linearized Einstein equations, using a code for the sourced Teukolsky equation with hyperboloidal slicing (\cite{zenginoglu:2011}). We present the waveforms in Fig.~\ref{fig5}. We show the three waveforms, for the energy balance approximation --- which is equivalent to the self force waveform with the conservative piece turned off --- and the two independent waveforms obtained from the direct approach and from the osculating geodesics approach that use the full self force, including its conservative piece. All three waveforms start in phase at early times. The two full self-force waveforms are dephased from the energy balance waveforms by much more than from each other. Specifically, the total commutative dephasing of either waveform from the energy balance waveform over the entire orbit from $r_0=10\,M$ down to near the ISCO, is $\,\Delta\phi=8.4\pm 0.4\;{\rm rad}$. This dephasing corresponds to about $4/3$ of a cycle, compared with the $54.3$ cycles the particle makes in the energy balance case.

We next study the importance of the conservative piece of the self force for the waveforms. Specifically, we take a window of duration $L$ from the end of the chirp part of the energy balance waveform. and find its overlap integral with the full self force waveform, which in this case models the actual data stream. As expected, the longer $L$, the smaller the overlap integral. When we take $L=L_{\rm threshold}=816.6\,M$ we find that the overlap integral is reduced to $0.96$, which corresponds to a loss of $10\%$ in the event rate.  Therefore, for $L\lesssim L_{\rm threshold}$ ignoring the conservative piece of the self force does not lead to a significant loss of accuracy in detection of E(I)MRI events. However, longer stretches of data, $L\gtrsim L_{\rm threshold}$ leads to a significant decrease in event rates, so that the inclusion of the conservative piece of the self force is important. Further detail of this work appears in \cite{burko:2012}.

The authors are indebted to Gaurav Khanna for discussions and for use of his numerical code.
This work was supported by a NASA EPSCoR RID grant and by NSF grants PHY--0757344 and DUE--0941327. LMB is grateful to Alessandro Spallicci for hospitality.


\section{Introduction}

In an $R$-parity conserving supersymmetric theory, the
lightest supersymmetric particle (LSP) is stable and thus
provides us with a candidate for dark matter (DM).
Searching for supersymmetry at colliders always
means searching for a certain LSP scenario. The neutralino
LSP scenario is clearly the most widely studied scenario for
the LHC and it is often referred to as \emph{the} SUSY search.
However, if one extends the MSSM, other cosmologically
viable DM candidates show up.
If, for instance,  we impose the idea of SUSY on
the fundamental theory of nature we have to supersymmetrize
gravity predicting the gravitino, the superpartner of the graviton.
The gravitino is indeed a very well motivated LSP\,\cite{Fayet:1981sq}
and may even be considered as favored since it alleviates the gravitino
problem\,\cite{Ellis:1984er} allowing for larger reheating temperatures.
Another perfectly viable DM candidate is the axino,\cite{Nilles:1981py}
a supersymmetric partner of the QCD axion
which is introduced in the Peccei-Quinn mechanism.
Both candidates are very weakly interacting particles---the couplings
are suppressed by the Planck- or the Peccei-Quinn scale, respectively.
This naturally renders the next-to-LSP (NLSP) long-lived.
If the decay length of the NLSP is large
compared to the size of a detector of a collider experiment, the LSP is
typically not involved in the interactions inside the detector at all and the
NLSP determines the signature. For an electrically charged NLSP, like the stau,
this leads to a very distinct SUSY signature---no missing energy
but highly ionizing tracks leaving the detector.\cite{Fairbairn:2006gg}
Several studies for long-lived staus have been performed, assuming
a certain underlying high scale model---such as minimal gauge mediated
supersymmetry breaking or the constrained MSSM (CMSSM).\footnote{%
Currently, searches for heavy stable charged particles are performed by
ATLAS\,\cite{Aad:2011hz} and CMS\,\cite{Chatrchyan1205.0272}.}
However, long-lived staus can basically occur in all common
breaking models.

In this work we address the question whether it is possible to derive robust bounds
on the mass parameters from the LHC experiment that cover all possible
spectra accommodating a long-lived stau NLSP. We will therefore adopt the
approach of simplified models.\cite{Alwall:2008ag,Alves:2011wf}
As an important outcome, we will see that in the case of long-lived staus
the simplified models provide a very powerful description covering a large set of
possible realistic spectra. We will restrict ourselves to the production via squarks
and gluinos as well as to direct pair production of staus.

\section{Simplified models for long-lived staus}\label{sec:models}

Simplified models have successfully been used in the neutralino LSP
scenario and have lead to bounds on the gluino and squark mass beyond
those derived within the CMSSM.
However, one immediately notices some important differences between that
scenario and the one considered here. In the neutralino LSP scenario, the only
signal from the LSPs is missing energy. The SUSY search is almost
independent of the neutralino kinematics, rather it depends strongly on the
type and hardness of the SM particle radiation. In the considered long-lived
stau scenario it is---roughly speaking---the other way around.
The kinematics of the staus, in particular their velocities are crucial quantities.
For velocities well below one the signature is basically background free,
whereas, if the velocity approaches one a comparatively huge amount of
muon background is present. Consequently, the significance of the signal
will drop drastically.

When the stau is produced in a cascade decay following the
production of colored sparticles a very large number
of parameters enter the computation through the
mixing angles and masses of the intermediate sparticles that
can participate in the cascade decay. But not all parameters
are equally important. Obviously, the masses of the squarks
and gluinos play an important role in determining the production
cross section. Additionally, the mass of the stau determines
the total phase space available in the cascade decay
of the lightest colored sparticle (LCP) down to the stau and
thus dominantly determines the stau velocity. The impact
of the intermediate sparticles can successfully be described
by introducing a few distinct simplified models which represent
limiting cases of realistic spectra and allow for the free parameters
$m_{\widetilde{g}}$, ${m_{\widetilde{q}}}$ and $m_{{\widetilde{\tau}}_1}$.

Considering a fixed mass gap $m_\textnormal{\tiny LCP}-m_{{\widetilde{\tau}}_1}$,
in the limit of massless SM radiation a non-trivial extremum of
the stau velocity is given by the mass pattern
\begin{equation}
\label{eq:masspatt}
m_i=(m_\textnormal{\tiny LCP})^{\frac{n-i}{n}}(m_{{\widetilde{\tau}}_1})^{\frac{i}{n}}
\end{equation}
of $n-1$ intermediate sparticles $i$. This extremum is usually
a minimum. The only other extremum
is represented by the mass degenerate configurations
\begin{equation}
\label{eq:massdeg}
m_i=m_\textnormal{\tiny LCP}\,,\quad m_j=m_{{\widetilde{\tau}}_1} \quad
\forall\; i<k,\; j\geq k\,, \quad k = 1,\dots,n \,.
\end{equation}
which is kinematically equivalent to a direct decay of the LCP into
the stau.
These observed extrema motivate the following simplified
models.
Model $\mathcal{A}$: The direct decay
over a nearly degenerate neutralino (${m_{\widetilde{\chi}^0}}\simeqm_{{\widetilde{\tau}}_1}$).
Model $\mathcal{B}$:
The mass pattern (\ref{eq:masspatt}) for $n=2$,
LCP$\to{\widetilde{\chi}^0_1}\to{\widetilde{\tau}}_1$, i.e.
${m_{\widetilde{\chi}^0}}=\sqrt{\mlcpm_{{\widetilde{\tau}}_1}}$.
Model $\mathcal{C}$: The mass pattern (\ref{eq:masspatt})
for $n=4$, LCP$\to\widetilde{\chi}^0_2\to\widetilde{\ell}\to{\widetilde{\chi}^0_1}\to{\widetilde{\tau}}_1$.
In fact, these models represent limiting cases even for
decay chains that radiate massive SM particles, like
 LCP$\to\widetilde{\chi}^0_2\to{\widetilde{\chi}^0_1}\to{\widetilde{\tau}}_1$. Decay chains longer than
 $n=4$ are rather unlikely and will not make up the dominant
 decay mode within the MSSM.\,\cite{Heisig:2012zq}

Another important difference between the neutralino LSP
and the considered scenario is that the latter provides a
perfectly well observable direct production mechanism.
The production via the Drell-Yan (DY) process depends only
on $m_{{\widetilde{\tau}}_1}$ and (less strongly) on the stau mixing angle $\theta_{{\widetilde{\tau}}}$.
The stau mixing angle that provides the lowest cross section,
$\theta_{{\widetilde{\tau}}}^{\min}$,
is practically independent of the mass and kinematics of
the stau (at least for stau masses above roughly $150\,\textnormal{GeV}$).
Thus the DY production provides an important lower limit
on the stau mass, which is completely
model-independent.\,\cite{Heisig:2011dr}

\begin{figure}[tbhp]
\begin{center}
\setlength{\unitlength}{1\textwidth}
\begin{picture}(0.85,0.74)
\put(0.0,0.4){
  \put(0.04,0.035){\includegraphics[scale=1.2]{mgo_mstau_R05_8TeV}}
  \put(0.14,0.0){\footnotesize $m_{\widetilde{g}}=1/2 \,{m_{\widetilde{q}}}\; [ \,\textnormal{GeV} \, ] $}
  \put(0.0,0.14){\rotatebox{90}{\footnotesize $m_{{\widetilde{\tau}}_1}\;[\,\textnormal{GeV}\,]$}}
  \put(0.09,0.316){\footnotesize $\,pp\,(8\,\textnormal{TeV})$}
  \put(0.09,0.286){\footnotesize $\int\!\mathcal{L}=16\,\textnormal{fb}^{-1}$}
  \put(0.136,0.073){\tiny minimal DY production limit}
  }
 \put(0.43,0.4){
  \put(0.04,0.035){\includegraphics[scale=1.2]{mgo_mstau_R1_8TeV.eps}}
  \put(0.15,0.0){\footnotesize $m_{\widetilde{g}}={m_{\widetilde{q}}}\; [ \,\textnormal{GeV} \, ] $}
  \put(0.0,0.14){\rotatebox{90}{\footnotesize $m_{{\widetilde{\tau}}_1}\;[\,\textnormal{GeV}\,]$}}
  \put(0.09,0.316){\footnotesize $\,pp\,(8\,\textnormal{TeV})$}
  \put(0.09,0.286){\footnotesize $\int\!\mathcal{L}=16\,\textnormal{fb}^{-1}$}
  \put(0.136,0.073){\tiny minimal DY production limit}
  }
\put(0.0,0.0){
  \put(0.04,0.035){\includegraphics[scale=1.2]{mgo_mstau_R2_8TeV.eps}}
  \put(0.14,0.0){\footnotesize ${m_{\widetilde{q}}}=1/2\, m_{\widetilde{g}} \; [ \,\textnormal{GeV} \, ] $}
  \put(0.0,0.14){\rotatebox{90}{\footnotesize $m_{{\widetilde{\tau}}_1}\;[\,\textnormal{GeV}\,]$}}
  \put(0.09,0.316){\footnotesize $\,pp\,(8\,\textnormal{TeV})$}
  \put(0.09,0.286){\footnotesize $\int\!\mathcal{L}=16\,\textnormal{fb}^{-1}$}
  \put(0.136,0.073){\tiny minimal DY production limit}
  }
 \put(0.43,0){
  \put(0.04,0.035){\includegraphics[scale=1.2]{mgo_msq_8TeV.eps}}
  \put(0.18,0.0){\footnotesize $m_{\widetilde{g}}\; [ \,\textnormal{GeV} \, ] $}
  \put(0.0,0.14){\rotatebox{90}{\footnotesize ${m_{\widetilde{q}}}\;[\,\textnormal{GeV}\,]$}}
  \put(0.215,0.316){\footnotesize $\,pp\,(8\,\textnormal{TeV})$}
  \put(0.215,0.286){\footnotesize $\int\!\mathcal{L}=16\,\textnormal{fb}^{-1}$}
 \put(0.215,0.26){\tiny along minima in}
 \put(0.215,0.244){\tiny $m_{{\widetilde{\tau}}_1}$-variation}
  }
\end{picture}
\end{center}
\caption{Projected LHC sensitivity ($95\%\,\textnormal{CL}_\textnormal{s}$
exclusion, see footnote \ref{footn}) for the models $\mathcal{A}$ (blue dashed),
$\mathcal{B}$ (green solid) and $\mathcal{C}$ (red dot-dashed) as well as
$\mathcal{A}$ (cyan dotted) and $\mathcal{C}$ (purple dot-dot-dashed) for
a reduced set of selection criteria (see text for details). In the lower right panel
the curves represent the minima in the sensitivity with respect to the variation
of $m_{{\widetilde{\tau}}_1}$. We assume a common squark mass for all flavors. Taken from
reference 9.\,\,
}
\label{fig:8TeV}
\end{figure}

\begin{figure}[tbhp]
\begin{center}
\setlength{\unitlength}{1\textwidth}
\begin{picture}(0.85,0.36)
\put(0.0,0.0){
  \put(0.04,0.035){\includegraphics[scale=1.2]{mgo_msq_7TeV.eps}}
  \put(0.18,0.0){\footnotesize $m_{\widetilde{g}}\; [ \,\textnormal{GeV} \, ] $}
  \put(0.0,0.14){\rotatebox{90}{\footnotesize ${m_{\widetilde{q}}}\;[\,\textnormal{GeV}\,]$}}
  \put(0.21,0.316){\footnotesize $\,pp\,(7\,\textnormal{TeV})$}
  \put(0.21,0.286){\footnotesize $\int\!\mathcal{L}=5\,\textnormal{fb}^{-1}$}
 \put(0.215,0.26){\tiny along minima in}
 \put(0.215,0.244){\tiny $m_{{\widetilde{\tau}}_1}$-variation}
  }
 \put(0.43,0){
  \put(0.04,0.035){\includegraphics[scale=1.2]{DY_exclusion.eps}}
  \put(0.18,0.0){\footnotesize $m_{{\widetilde{\tau}}_1}\;[\,\textnormal{GeV}\,]$}
  \put(0.0,0.17){\rotatebox{90}{\footnotesize $\sigma\; [ \,\textnormal{pb} \, ] $}}
  \put(0.215,0.316){\footnotesize $\,pp\,(7\,\textnormal{TeV})$}
  \put(0.215,0.286){\footnotesize $\int\!\mathcal{L}=5\,\textnormal{fb}^{-1}$}
 \put(0.215,0.256){\footnotesize minimal direct}
 \put(0.215,0.231){\footnotesize DY production}
 \put(0.1,0.324){\tiny theoretical}
  \put(0.11,0.308){\tiny prediction}
 \put(0.1,0.125){\tiny CMS limit}
 \put(0.1,0.108){\tiny Tracker+ToF}
  }
\end{picture}
\end{center}
\caption{Left panel: Projected LHC sensitivity
($95\%\,\textnormal{CL}_\textnormal{s}$ exclusion,
see footnote \ref{footn}) for the models $\mathcal{A}$ (blue dashed),
$\mathcal{B}$ (green solid) and $\mathcal{C}$ (red dot-dashed).
The curves represent the minima in the sensitivity w.r.t. the variation
of $m_{{\widetilde{\tau}}_1}$. We assume a common
squark mass for all flavors. Right panel:
The cross section limit on stau pair production from the
$7\,\textnormal{TeV}$ CMS data and the NLO
cross section prediction for the DY production with
$\theta_{{\widetilde{\tau}}}^{\min}$.
}
\label{fig:7TeV}
\end{figure}

\section{LHC reach}

To estimate the projected LHC reach of the scenario we
performed a full-fledged Monte Carlo study. We imposed a set of
selection criteria which was optimized in order to provide
high efficiencies throughout the considered parameter space.
This was achieved by only taking into account the staus and
up to two jets. The most important background is the DY production of
muons, even when two hard jets ($p_T>200\,\textnormal{GeV}$) are required.

Figure \ref{fig:8TeV} shows the projected $95\%\,\textnormal{CL}_\textnormal{s}$
exclusion limits\,\footnote{Due to the special situation of a very clean
signal region and thus only few required events these curves also represent an
(even conservative) estimation on the $5\sigma$ discovery potential.\label{footn}}
at the $8\,\textnormal{TeV}$ LHC for the three simplified models
defined in section \ref{sec:models}.
In model $\mathcal{A}$ (blue dashed lines) the decay LCP$\to{\widetilde{\chi}^0_1}$
is the dominant decay leading to hard jets and potentially highly boosted staus.
In an intermediate region of $m_\textnormal{\tiny LCP}-m_{{\widetilde{\tau}}_1}$ this leads to an
enhancement of the significance whereas for large mass gaps
$m_\textnormal{\tiny LCP}-m_{{\widetilde{\tau}}_1}$ the significance drops sharply, despite the fact that the
jets become harder. (The cyan dotted curve shows the sensitivity
without taking a jet-optimized selection criterium into account.)
The scenario would hide very effectively from observation if it
were not for the direct production which increases the sensitivity for
lower $m_{{\widetilde{\tau}}_1}$. This results in a minimum in the sensitivity
just above the limit provided by the minimal DY production.
For small mass gaps $m_\textnormal{\tiny LCP}-m_{{\widetilde{\tau}}_1}$ the lower limit on the
velocity becomes more important. Such a lower limit might be
imposed by trigger restrictions. The purple dot-dot-dashed
curves show the sensitivity for $\mathcal{C}$ after dropping
a selection criterium that allows staus to fire the muon trigger in delay,
i.e. requiring a buffering of the tracker data up to three bunch crossings
in delay.

In the lower right panel of figure \ref{fig:8TeV} and the left panel of figure
\ref{fig:7TeV} the LHC sensitivity in the $m_{\widetilde{g}}$-${m_{\widetilde{q}}}$-plane is shown.
The curves show the $95\%\,\textnormal{CL}_\textnormal{s}$ taking the
respective $m_{{\widetilde{\tau}}_1}$ that yields the lowest cross section at each point.
The simplified models $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$
span a relatively narrow band in the $m_{\widetilde{g}}$-${m_{\widetilde{q}}}$-plane.
From the non-observation of heavy stable charged particles
this allows for an estimation of model-independent limits on
$m_{\widetilde{g}}$ and ${m_{\widetilde{q}}}$.
In the right panel of figure \ref{fig:7TeV} the cross section for direct DY
production of staus is displayed as a function of the stau mass. We
chose here the stau mixing angle which yields the smallest cross section.
Additionally, we display the cross section limit curve from
CMS\,\cite{Chatrchyan1205.0272} obtained for another choice of
$\theta_{{\widetilde{\tau}}}$. Since the kinematics do not depend on $\theta_{{\widetilde{\tau}}}$,
we do not expect any changes in the cut efficiencies.

\section{Conclusion}
The long-lived stau scenario can be very well described by simplified
models. The direct detectable stau does not give rise to regions
in parameter space where the scenario hides very effectively from
observation, as it is e.g. the case in the neutralino LSP scenario where
compressed or widely spread spectra are difficult to explore.
From the minimal direct production, long-lived staus with a mass below
$m_{{\widetilde{\tau}}_1}=216\,\textnormal{GeV}$ can be excluded. This limit is completely
model-independent. In the described framework, conservative limits
on the gluino mass and a common squark mass of about $m_{\widetilde{g}}>1100\,\textnormal{GeV}$
and ${m_{\widetilde{q}}}>1250\,\textnormal{GeV}$ can be estimated.

\section*{References}
\providecommand{\href}[2]{#2}\begingroup\raggedright