On Integrability of the Geodesic Deviation Equation

The Jacobi equation for geodesic deviation describes finite size effects due to the gravitational tidal forces. In this paper we show how one can integrate the Jacobi equation in any spacetime admitting completely integrable geodesics. Namely, by linearizing the geodesic equation and its conserved charges, we arrive at the invariant Wronskians for the Jacobi system that are linear in the `deviation momenta' and thus yield a system of first-order differential equations that can be integrated. The procedure is illustrated on an example of a rotating black hole spacetime described by the Kerr geometry and its higher-dimensional generalizations. A number of related topics, including the phase space formulation of the theory and the derivation of the covariant Hamiltonian for the Jacobi system are also discussed.


I. INTRODUCTION
The Jacobi equation of geodesic deviation, e.g. [1], has an important place in General Relativity. The motion of a small enough particle with no internal structure such as charge or spin is described by a geometrical object, a geodesic, so that by the principle of equivalence an observer freely falling with the particle would not be able to feel the effects of gravity in a small neighbourhood of spacetime. Gravitational effects will become apparent considering objects of a finite size, whose evolution can be thought of as that of a bundle of nearby geodesics: at first order in the separation parameter between the geodesics these effects are described by the Jacobi equation. This happens for example when studying the effects of a passing gravitational wave, as in the gravitationalwave memory effect [2][3][4][5][6]. The Jacobi equation can be generalized to the case of particles with electric charge [7], particles with spin [8], or to the non-linear case where the dependence on relative velocities is not linearized [9].
The geodesic deviation has been used to give a geometrical and physical interpretation of spacetimes in ordinary four dimensions and higher [10][11][12], while using first and higher order it has been used to construct approximations to generic geodesics starting from simple * marco.cariglia@ufop.edu.br † t.houri@maizuru-ct.ac.jp ‡ Pavel.Krtous@utf.mff.cuni.cz § dkubiznak@perimeterinstitute.ca ones [13,14], that can be used to model extreme massratio systems [15]. In another approach, it has been used to study geodesic (in)stability for dynamical systems, and Lyapunov exponents [16,17].
In this work we analyze the Jacobi equations from the point of view of Hamiltonian dynamics and symmetries of the dynamics. The Jacobi dynamical system is nontrivial, even though it is linear, because of its time dependence and the complexity of its differential equations. Time dependent Hamiltonian systems have been recently discussed in [18] in the light of the Eisenhart lift technique [19][20][21][22]. On the other hand, due to their complexity, few explicit solutions of the equations are given in the literature [7,13,23,24].
From this point of view it is interesting that hidden symmetries of the related geodesic equations generate non-trivial solutions of the Jacobi equation [25][26][27], displaying the relationship between the Jacobi equations and the concept of hidden symmetries.
In fact the equations of geodesic deviation are naturally linked to those of geodesic motion. In [28] a method is presented that generalizes the Hamilton-Jacobi equation, allowing to obtain at the same time solutions of the geodesic and the geodesic deviation equations. A Lagrangian formulation of the geodesic deviation equations, including an electromagnetic field and spin, and a treatment of higher order deviation equations can be found in [7,13,29], and is obtained by an expansion of that of geodesic motion. Given this natural connection, it is reasonable to expect that symmetries of dynamics of the geodesic motion, when present, descend to symmetries of the Jacobi equation. This is the focus of the present work, where we show that integrals of geodesic motion give rise, through linearization, to integrals of the Jacobi equation that are expressed in the form of invariant Wronskians. Of particular interest are the geodesic integrals built from Killing vectors and Killing tensors, which are associated to hidden symmetries of the geodesic dynamics and that generate integrals for the Jacobi equation. The latter integrals inherit an algebraic structure via Poisson brackets that is isomorphic to the structure of the geodesic integrals. In particular, if in n dimensions there are n functionally independent, mutually Poisson commuting integrals of the geodesic motion, then these induce via linearization a set of n independent, mutually commuting integrals for the Jacobi equation, thus showing that integrability of the geodesic equations implies integrability of the Jacobi equations.
The structure of the work is as follows. We begin in Sec. II by reviewing the concept of geodesic motion and setting some of the notation. In Sec. III we present the Jacobi equation and discuss its Lagrangian and Hamiltonian formulation. We first discuss a coordinate approach where the Hamiltonian is not a scalar with respect to changes of coordinates in the base manifold: in this case there is a coordinate Hamiltonian for each coordinate chart, different such Hamiltonians being related by canonical transformations in overlapping regions, and global motion is obtained sewing solutions from different charts. We then introduce a covariant Hamiltonian approach where the Hamiltonian is a globally defined scalar. In Sec. IV we discuss how geodesic integrals of motion descend to integrals of the Jacobi equation. We present the notion of invariant Wronskians, followed by the conserved quantities generated by Killing tensors and then discuss integrability. In Sec. V we apply our results to rotating black holes in four dimensions, and mention how these generalize to higher dimensional Kerr-NUT-(A)dS black holes. The Jacobi equation is integrable in these geometries and from our results it is possible to build a complete set of mutually commuting conserved charges. Sec. VI presents concluding remarks and possible future lines of research. Appendix A presents the Jacobi equation, the linearization procedure and the Wronskians from the point of view of a general phase space and a general set of equations of motion, before the introduction of a cotangent bundle or a specific Hamiltonian, and Appendix B details the construction of a covariant Lagrangian and Hamiltonian.

II. GEODESIC MOTION
Let M be a spacetime manifold of dimension n, equipped with metric g ab , and Riemann tensor 1 R abcd .
Having in mind applications to relativity we take the metric to be of almost plus type, though with minor adjustments everything would readily generalize to a generic metric.
There are several approaches to the Lagrangian and Hamiltonian descriptions of the relativistic particle. We follow the approach used by Carter [30] and others , where one identifies the spacetime M with a configuration space of the system. The particle is described by its trajectory x a (λ) parametrized by an external time λ. The phase space is then described by the cotangent space T * M. The Lagrangian for geodesic motion is given by The momentum and the Hamiltonian read This gives, of course, the geodesic equation. This formalism describes a free particle of an arbitrary mass. (In this paper we concentrate on the case of massive particles for which dx a /dλ is not a null vector.) The mass is fixed by the value of the Hamiltonian, i.e., by the normalization of the momentum Since the Hamiltonian does not depend on the external time explicitly, it is a conserved quantity and the mass is for a given trajectory fixed. We can always rescale the time variable and introduce the proper time τ τ = mλ , (2.6) with the normalization of the velocity u a Of course, free particles of different (non-zero) masses follow the same geometric trajectories-they differ only in time parametrization. We can thus ignore a particular value of the mass and describe the geodesic in proper time parametrization. In what follows we consider a situation where the geodesic motion admits a nontrivial number of integrals of motion. In particular, the conserved quantities that are homogeneous in particle's momentum, (2.8) are in one to one correspondence with Killing tensors [31]. A Killing tensor of rank r is a symmetric tensor K a1...ar = K (a1...ar ) , such that For r = 1 it reduces to the Killing vector. It is the purpose of the present paper to linearize these conserved quantities and show that they give rise to simple integrals of motion for the Jacobi geodesic deviation equation described in the next section. In particular, it implies that when the original geodesic motion is completely integrable, so will be the Jacobi equation.

A. Geodesic deviation equation
Given a geodesicx a (τ ), we want to study its nearby trajectories.
To that purpose we consider a oneparameter family of curves x a (σ, τ ), with σ being the parameter labeling different curves and τ the time parameter along each curve. We call the geodesicx a (τ ) = x a (0, τ ) a central geodesic. We assume τ to be the proper time along the central geodesic, however, it does not have to be the proper time for trajectories with nonvanishing σ. Nevertheless, we still denote the velocity with respect to τ as u a , We call by n a a vector that links the nearby curves, For σ = 0,ū a (τ ) ≡ u a (0, τ ) reduces to the normalized velocity of the central geodesic and n a (τ ) ≡ n a (0, τ ) describes the deviation from the central geodesic. One can think of this vector as a linear approximation for the trajectories close to the central geodesic.
A special degenerate case occurs when the whole family x a (τ, σ) lies entirely on the central geodesic. Clearly, n a is then proportional to u a .
Specifically, we are interested in the trajectories in the vicinity of the central geodesic which are geodesics as well, i.e., curves that (for all values of the parameter σ) satisfy In this case it is well known that the deviation n a (τ ) connecting nearby geodesics must satisfy the geodesic deviation Jacobi equation, Here, the bar above the Riemann tensor (and, similarly, above other quantities) indicates that it is evaluated at the central geodesic,R abcd = R abcd (x). In order to derive this equation, the key observation is to realize that since u a and n a are essentially coordinate vectors, u ≡ ∂ τ , n ≡ ∂ σ , on the 2-surface x a (σ, τ ), they Lie-commute, [u, n] = 0, which when expressed in terms of the metric covariant derivative yields see e.g. [1] for more details.

B. Lagrangian for the Jacobi equation
The Jacobi equation admits a Lagrangian formulation, with the Lagrangian l given by [7] This is the Lagrangian for independent 'deviation' variable n a which represents a general curve close to the central geodesic. As we will discuss below, it can be understood as a function l(n a ,ṅ a ) of coordinate velocitẏ n a = dn a dτ , or it can be treated 'covariantly', as a function l(n a , Dn a dτ ) of covariant velocity Dn a dτ . In either case the corresponding Euler-Lagrange equations pick up a nearby geodesic specified by the Jacobi equation.
The Lagrangian (3.6) can be derived by the linearization process, starting from the Lagrangian of the geodesic motion. 2 We refer to the Appendix B for the derivation and further technical details.
The Lagrangian (3.6) is obviously time dependent-it depends on the time-dependent positionx(τ ) and velocityū(τ ) of the central geodesic. Therefore, there is one such Lagrangian for each central geodesic. While the dependence on the velocityū(τ ) is explicit, the dependence onx(τ ) enters through the spacetime dependence of the metric, Christoffel symbols, and the Riemann tensor, and to stress this fact we write these objects with bar. All these should be considered as given functions of the time variable τ and provide apriori data for the Jacobi equation. Of course, for a concrete spacetime, it may be difficult to obtain the expression forx(τ ) andū(τ ) in an explicit and closed form. We return to the question of time-dependency of the Lagrangian again below, when we discuss it from a covariant perspective.

Transverse and tangent splitting
Before we proceed to the Hamiltonian formulation, let us first discuss decoupling of the transverse and parallel degrees of freedom in the deviation variable. The Jacobi equation (3.4) always admits the following two solutions: which arise from the degenerate case when the entire family x a (σ, τ ) lies on the central geodesic x a (τ ): . The solutions of (3.10) are precisely those in Eq. (3.7). Correspondingly, for timelikeū a , the Lagrangian separates as l = l + l ⊥ , with as follows from the fact that (3.14) When the geodesics are timelike, the metric can be decomposedḡ ab = −ū aūb +ḡ ⊥ab , (3.15) whereḡ ⊥ab is the projector to the transverse space, which is positive definite, and (3.13) can be written as

Hamiltonian formalism: coordinate approach
In what follows we shall study the Jacobi equation from the point of the Hamiltonian dynamics using both, coordinate (this subsection) and covariant (next subsection) approaches.
In the coordinate approach, the Lagrangian (3.6) is understood as a function of the 'positions' n a and the 'coordinate velocities'ṅ a , kalbū kūl n a n b .
(3.17) To write down the corresponding Hamiltonian formulation, we define the momentum canonically conjugated to n a π a = ∂ l ∂ṅ a = mḡ ab ṅ b +ū kΓb kc n c = mḡ ab and introduce what we call the coordinate Hamiltonian, h c = π aṅ a − l, (3.19) An explicit calculation shows that the equations of motion obtain from h c imply the Jacobi equation (3.4).
It is obvious that the Hamiltonian (3.19) is not covariant: it depends in a non-tensorial way on a particular choice of coordinates through the Christoffel symbols Γ a bc , reflecting the fact that we used a non-covariant form of the velocityṅ a . Such a velocity does not transform as a vector under a coordinate transformation and therefore the related Hamiltonian is also non-covariant. The Hamiltonian (3.19) generates the coordinate timeevolution of vector quantities.
As we will see in the next subsection, it is possible to proceed in a more covariant way, starting with the covariant velocity Dn a Dτ , and arrive at a simpler covariant Hamiltonian (3.27) below. This, however, requires an additional care about technical details and before we explore such approach, let us first discuss the behavior of the coordinate Hamiltonian under a change of coordinates.
Evolving the parameter τ , there can come a moment when the geodesic leaves the given coordinate chart {x a }. Then it is necessary to use a different set of coordinates From the point of view of Hamiltonian dynamics, such a change of coordinates induces a time-dependent canonical transformation of (n a , π a ) variables, and the Hamiltonian in general changes. The non-invariance property is expressed by the non-tensorial term in h c . To be explicit, the change of coordinate chart (3.20) induces the following transformation: Such a transformation is canonical and amounts to the following generating function: By this we mean that Eqs. (3.21) are obtained via solving for n ′ and π ′ . The standard theory then gives that the Hamiltonian transforms to where the term ∂G ∂τ provides the non-tensorial term in the transformation rule of the connection. The Hamiltonian is thus locally invariant in form, although it is not build from pure tensorial expressions. However, there is no such thing as a 'global coordinate Hamiltonian'. For each coordinate chart there is a different Hamiltonian which governs time evolution in that chart. These Hamiltonians differ by ∂G ∂τ terms, and so they cannot be understood as just different coordinate expressions for one Hamiltonian function. The global motion has to be sewed from solutions in different charts. Nevertheless, the global covariant Hamiltonian can be defined in the covariant approach as we will show next.

Hamiltonian formalism: covariant approach
The linearized configuration space of vectors n a is time dependent: it is a tangent space atx(τ ). In the discussion of the coordinate Hamiltonian above, we have implicitly identified such tangent spaces by choosing particular coordinates. Namely, we have naturally identified vectors with the same components with respect to the coordinate frame. More precisely, we regarded such vectors as 'not changing'. We have used the coordinate-time derivativeṅ a = dn a dτ (τ ) to define the velocity and employed it in the construction of the Hamiltonian. Since such an identification of vectors is non-covariant, we obtained a non-covariant Hamiltonian which was dependent explicitly on the Christoffel symbols.
However, we can identify the linearized configuration spaces at different times in a more covariant way: by a parallel transport along the central geodesic. For that, we use the covariant time derivative Dn a Dτ to define the velocity. 3 Thus, in the covariant approach we must interpret the Lagrangian l as a function of a linearized position n a and of a covariant velocity v a = Dn a Dτ , A covariant version of canonically conjugate momentum π a then reads The Hamilton equations have to be written again using the covariant time derivative (3.28) Combining both equations together yields the Jacobi equation, The evolution of a general phase space observable A expressed in new variables n, π is given by where the Poisson bracket assumes the standard form Let us stress that even the observables, that were independent of the time parameter in the original variables (x, p), typically become explicitly time dependent and the second term in (3.30) is non-trivial. The reason is that a 'simple' function of x and p is re-expressed in terms of the central trajectoryx,p and linearized variables n and π. The explicit time dependency then enters through the time dependent central trajectory.
In particular, this is true for the covariant Hamiltonian itself, and we have that is, the new Hamiltonian h for the linearized system is not conserved. However, for the Hamiltonian, the time dependency is not caused only by the introduction of the linearized variables with respect to the central trajectory, but it is also related to a time-dependent canonical transformation that relates the original Hamiltonian H and its (linearized) version h, see appendix B for more details.

IV. INTEGRALS OF MOTION FOR THE JACOBI SYSTEM
In this section we will describe how the integrals of motion for geodesics give rise to the particularly simple integrals of motion for the Jacobi system. We start with a discussion applicable to any linearized dynamical system and only later specify to the case of the geodesic motion.

A. Wronskian as a linear integral of motion
Consider a linearized dynamical system, with trajectories near the central trajectoryx described by a linearized trajectory n a (τ ). A general feature of linearized systems is that their dynamics is governed by a quadratic Lagrangian l(n, v), and, in the Hamiltonian picture, by a quadratic Hamiltonian h(n, π), see (3.25) and (3.27) for the specific example of linearized geodesic motion.
For two linearized trajectories n a 1 (τ ) and n a 2 (τ ), we define the Wronskian as W [n 1 |n 2 ] = n a 1 ∂ l ∂v a n 2 , (4.1) In the phase-space variables the trajectory is characterized by position n a and momentum π a and the Wronskian can be written as It is well known (see, e.g., [25,26] for the case of the Jacobi system) that for any two solutions n a 1 (τ ) and n a 2 (τ ) of the equation of motion the Wronskian is conserved in time τ . To show this, let us use the Hamiltonian picture. By employing the general quadratic Hamiltonian, we have the following Hamiltonian equations: (4.4) Taking the time derivative of (4.2) and substituting (4.4) for Dn1 Dτ , Dπ1 Dτ and Dn2 Dτ , Dπ2 Dτ , we find This means that any fixed solutionñ(τ ) generates a quantity which is conserved along any solution n(τ ). In the phasespace language, any solutionñ(τ ),π(τ ) defines a conserved quantity Wñ ,π (n, π) = W [ñ,π|n, π] =ñ a π a −π a n a . Clearly, such a conserved quantity is linear in n and π. This observation can be reversed: the most general linear conserved quantity of a linearized system is given by Wñ ,π , whereñ(τ ),π(τ ) are explicit solutions of the Hamilton equations. Indeed, let us consider a general linear observable C =ñ a π a −π a n a with yet unspecified coefficientsπ a (τ ) andñ a (τ ). Its conservation means (4.8) Substituting (4.3), re-arranging terms and using (4.3) again, we get 0 = π a Dñ a Dτ − ∂h ∂π a (ñ,π) − n a Dπ a Dτ + ∂h ∂n a (ñ,π) . (4.9) Since C should be conserved at any phase-space point (n, π), we obtain thatñ a (τ ) andπ a (τ ) must satisfy the Hamilton equations and, thus, the observable C has to be of the form C = Wñ ,π .
A similar statement can be obviously formulated in the configuration language: any conserved quantity linear in the trajectory n(τ ) and its time derivative has to have the form (4.6) for a solutionñ(τ ).
Finally, thanks to the linear structure of the linearized system, the Wronskian of two solutions n 1 (τ ), n 2 (τ ) can be related to the Poisson bracket of the corresponding conserved quantities W n1,π1 and W n2,π2 , (4.10) Let us finally return back to the geodesic motion and the corresponding linearized Jacobi system. In this case the Wronskian (4.1) takes the particular form  11) and the above formulae directly apply. In what follows we concentrate on this case.

B. Canonical observables
It is easy to see that the Poisson bracket of two observables linear in (n, π) is equal to a constant, i.e., an observable independent of (n, π). This observation can be used to construct a set of (time dependent) observables for the Jacobi system (F j , G j ), j = 1, . . . , n, which form canonical coordinates at all times, For example, choosing at time τ 0 an orthonormal frame of vectors e (i) and the dual frame of 1-forms e (i) , both at x(τ 0 ), one can define the solutions of the Jacobi system f j (τ ) and g j (τ ) with initial values at τ 0 given by (4.13) Clearly, at time τ 0 , and since the Wronskian is conserved, the relation (4.14) remains true at all times. Using this solution, and thanks to (4.10), we can thus define canonical coordinates (F j , G i ) satisfying (4.12) by the corresponding Wronskian observables In particular, the set of coordinates F j (as well as the set of G i ) forms a maximal set of commuting conserved quantities of the linearized system. However, these conserved quantities are not very useful. To find them, one has to find first the solutions f j and g j , i.e., to solve the linearized system. In the following, we want to discuss more useful conserved quantities-given by the symmetries of the spacetime. To find these observables, in addition to symmetries one only needs to know the central trajectory.

C. Conserved quantities generated by Killing tensors
As we reviewed in Sec. II, generic (homogeneous in momentum) integrals of geodesic motion are generated by Killing tensors, see (2.8). As shown by Caviglia, Zordan and Salmistraro (CZS) [25,32] these tensors also generate the following linearized solutionsñ a (τ ),π a (τ ) for the Jacobi system: Here, K is the conserved quantity of geodesic motion generated by the Killing tensor K a1...ar , (2.8), ∂K ∂p denotes the derivative with respect to momentum p with x fixed, and n a ∇aK ∂x (x, p) is the covariant derivative in direction n a with p parallelly transported, cf. (A21) in Appendix A. The latter derivative acts only on x-dependent terms in K and essentially ignores momentum p.
Let us verify the solution (4.16) by checking the Hamilton equations (3.28). Taking advantage of the fact that the central trajectory is geodesic, Dpa Dτ = 0, and p a = mū a , we obtain Using the identity which follows from the Killing condition (2.9), yields For momentumπ a we get (4.20) Here, we have used the Ricci identity and the fact that all r terms with the Riemann tensor give the same contribution. Now, the first term vanishes thanks to (2.9) and in the second term we can identifyñ a , which concludes the proof. The CZS solution (4.16) is of geometrical nature. It is obtained from the canonical transformation associated with the hidden symmetry of the geodesic equation. Let {·, ·} g be the Poisson bracket of the full geodesic theory, cf. (A22). Thenñ a = {x a , K} g andπ a = {p a , K} g . This is the same as the infinitesimal transformation δx a , δp a of the central geodesic generated by the canonical transformation induced by K.
The solution (4.16) can be related to the linearization k of the conserved quantity K. The expansion of any phase-space observable K(x, p) to the first order can be written as Obviously, using (4.16) we get k ≡ K −K =ñ a π a −π a n a = Wñ ,π . (4.23) The linearized observable k = Wñ ,π is thus again a conserved quantity that is linear in position and momentum and generated by the CZS solution (ñ,π). Let us finally mention that the CZS solution forñ a , (4.16), need not be generated from a Killing tensor. As noted in [25,32] (see also [33,34]), its existence is in oneto-one correspondence with a new object, called the affine tensor. An affine tensor of rank r, K a1...ar = K (a1...ar) , is an object that satisfies That is, the definition of an affine tensor 'generalizes' that of a Killing tensor by requiring that its symmetrized derivative need not vanish but can be a covariantly constant tensor. Of course, the above presented construction of conserved quantities through Wronskians immediately generalizes to the CZS solutions generated by affine tensors. Let us stress, however, that although the Killing tensors are formally a subfamily of affine tensors, the requirement on the existence of non-trivial h aa1...ar is very strong and at the moment there are no known physical spacetimes admitting affine tensors that are not at the same time Killing tensors. For this reason we shall not probe this possibility in this paper any further.

D. Integrability of the linearized system
Let us now assume that we have at least two Killing tensors K a...  corresponding to the conserved quantities K 1 and K 2 of the full geodesic motion. In general, such integrals of motion do not Poisson-commute. Their Poisson bracket generates a new conserved quantity K, which corresponds also to a Killing tensor K a... , given by the (symmetric) Schouten-Nijenhuis bracket [35-37] cf., e.g., [38]. For the linearized quantities k 1 , k 2 we have (4.27) The Wronskian can be expressed using the quantities related to the central trajectory. Substituting (4.2) and (4.16), we get the result already shown in [32]. It follows that if the original integrals of motion K 1 , K 2 Poisson-commute, K = 0, the linearized conserved quantities k 1 , k 2 also Poisson-commute. In particular, if the spacetime geometry possesses a full set of commuting integrals of motion K j generated by Killing tensors, the linearized system has also a full set of mutually commuting integrals of motion k j ≡ Wñ j ,πj . The complete integrability of the geodesic motion thus naturally implies the complete integrability of the Jacobi system.

V. INTEGRABILITY OF JACOBI EQUATION IN ROTATING BLACK HOLE SPACETIMES
Let us now apply the above developed formalism to explicitly demonstrate the integrability of the Jacobi equation in the Kerr black hole spacetime [39] and its higher-dimensional generalizations. Such an integrability stems from the existence of hidden symmetries in these spacetimes and derives from the integrability of the full geodesic motion. The same results remain also true for charged black holes and will be discussed elsewhere.
The Kerr metric represents a unique rotating black hole solution of vacuum Einstein equations. In the Boyer-Lindquist coordinates it reads where Σ = r 2 + a 2 cos 2 θ, and ∆ = r 2 − 2M r + a 2 . The metric admits two Killing vectors ∂ t and ∂ ϕ . For a geodesic motion, these imply the following conserved charges: In addition, there is a hidden symmetry encoded in the Killing tensor K ab that gives rise to Carter's constant [30,31,40] Lastly there is a conserved quantity generated by the metric g ab , seen as a (covariantly constant) Killing tensor, The four integrals of motion, {K E , K L , K C , K m 2 }, are functionally independent and mutually Poisson commute, yielding the geodesic motion completely integrable [40]. The explicit solution in terms of special functions can be for example found in [41,42].
We can now pick our favourite geodesic and turn to the corresponding linearized Jacobi system. The CZS solution, (4.16), yields the following independent solutions: From these we construct the independent conserved quantities for the Jacobi equation in Kerr k i =ñ a i π a − n a Dñ ia Dτ , i = 1, . . . , 4 . (5.10) We set the constant values of the Wronskians to w i , i = 1, . . . , 4, and introduce the abbreviated notation g i = n a Dñia Dτ , where the g i do not depend on the momenta π. From (5.10) it is easy to extract the value of the momenta π t and π ϕ , which are given by π t = π t (x a , n a ) = −g 1 − w 1 , (5.11) π ϕ = π ϕ (t,r,θ,φ) = g 2 + w 2 . (5.12) For geodesics with p r = 0, p θ = 0 it is possible to invert (5.10) with respect to π r , π θ . The result is π r = α r + β r π ϕ + γ r π t , (5.13) π θ = α θ + β θ π ϕ + γ θ π t , (5.14) with Here we have set f 3 = w 3 + g 3 , f 4 = w 4 + g 4 , these are functions of n a and not of the momenta. These expressions are involved, although in a closed form: it is a good example of the fact that the Jacobi equation is complicated even if it is linear. Let us finally mention that the procedure described in this section directly generalizes to higher-dimensional Kerr-NUT-(A)dS black hole spacetimes [43]. Such spacetimes are known to admit a number of hidden symmetries that yield the geodesic motion completely integrable [44,45]. It follows that the corresponding Jacobi system is also integrable and in principle can be solved by the same steps described in this section. Let us, however, stress that in higher dimensions, the generic geodesic is given only in terms of complicated integrals [38], see also [46][47][48] for special cases, and the solution of the Jacobi system thus becomes far from explicit.

VI. CONCLUSIONS
In this paper we have analyzed the Jacobi geodesic deviation equation from a point of view of Hamiltonian dynamics. It represents a dynamical system that is (although linear) explicitly time dependent. Consequently, the coordinate Hamiltonian is not covariant and varies from chart to chart. Nevertheless, we have shown that a covariant Hamiltonian can be constructed and shown (see Appendix B) how it can be obtained by the canonical transformation (accompanied by a due linearization) of the geodesic Hamiltonian. Although the geodesic Hamiltonian is a constant of motion, the linearized Hamiltonian for the geodesic deviation depends explicitly on time.
The main result of our paper regards the observation that the integrals of geodesic motion give rise to the corresponding integrals for the Jacobi system that are linear and given by the invariant Wronskians. In particular, this is true for the integrals generated by hidden symmetries of the spacetime. We have shown that if the geodesic motion is completely integrable, so will be the corresponding linearized motion described by the Jacobi equation. This has been further illustrated on an example of rotating black hole spacetimes in four and higher dimensions.
There is a number of topics that we have not discussed and that we point out as suitable for future research. One of these is the inclusion of spin. For example, a Lagrangian for the Jacobi deviation in the presence of Grassmannian spin variables can be found in [8], and a discussion of conserved quantities for geodesics in the presence of spin in Kerr-NUT-(A)dS spaces in [49]. Another one is the possibility of extending our results to higher order geodesic perturbations. These have been used to build analytic approximations of generic geodesics from simple exact solutions [13,14], and have been used to model extreme mass-ratio systems [15]. It would be interesting to find out if for example the standard conserved charges of Kerr can be used to build conserved charges for higher order geodesic perturbations. Lastly our results can be used to discuss the issue of (in)stability of dynamical systems [16,17]. In this appendix we want to elucidate the origin of the Jacobi equation from the point of view of the phase space formalism. We start by discussing the linearized evolution of a completely general system and re-derive the material from section IV for a general phase space, without the cotangent bundle structure. Specifying next to a phase space with the cotangent bundle structure and the Hamiltonian given by the configuration space metric, we show how to return back to the previous formalism.

Linearized phase space and Wronskian
In general, the phase space is a symplectic space with the symplectic structure Ω. It allows to define the Poisson brackets of two observables 4 and the Hamiltonian vector flow X F associated with an observable F Given a Hamiltonian H, the time evolution (classical trajectories) are given by orbits of the Hamiltonian flow X H . Assuming such a congruence of classical phase-space trajectories generated by X H , let us pick up one particular trajectoryX (τ ), which we call the central trajectory. We want to study a one-parametric subfamily X (σ, τ ) of these classical trajectories, which is close to the central trajectory, X (0, τ ) =X (τ ). In a linear approximation, such a family is generated by a phase-space vector field N (τ ) along the central trajectorȳ X (τ ), We call N (τ ) a linearized trajectory based on the central trajectoryX (τ ). We can extend the definition of N , (A4), also to nonzero values of σ. Since the vector fields X H and N can be viewed as coordinate fields X H ≡ ∂ τ and N ≡ ∂ σ on a two-dimensional sheet X (σ, τ ) with coordinates σ and τ , they must Lie-commute, This equation can be also rephrased as a condition on N (τ ) along the trajectoryX (τ ), Either of the last two equations can be viewed as the equation of motion for the linearized trajectory N (τ ). The space of phase-space vectors atX (τ ) represents the linearized phase space at time τ . It is a linear symplectic space with the symplectic structure given by the original symplectic structure Ω AB evaluated at X (τ ). Similar to the previous discussion, where the phase space has been represented by pairs (n, π), the linearized phase spaces at different times are geometrically different spaces. Therefore, we cannot directly apply the standard Hamiltonian formalism. We have seen that for that we would need to identify somehow the spaces at different times. However, let us discuss general phase-space situation without such an identification first. To this purpose we use the fact that the linearized evolution is given by the condition (A6) which naturally relates the phase spaces at different times.
As we have seen, the linear structure of linearized phase spaces allows us to define the Wronskian of two phase-space vectors N 1 , N 2 . It is given by the symplectic structure as cf. (4.2). It immediately follows that it is conserved along the time evolution Here, we have used (A6) and the fact that the symplectic structure is conserved along any Hamiltonian flow, L X F Ω = 0. Any chosen linearized solutionÑ (τ ) thus defines a conserved quantity WÑ on linearized solutions,

Conserved quantities and integrability
Let us now assume that the original system admits an integral of motion K, i.e., that there exists an observable K (not explicitly dependent on time parameter τ ) which Poisson-commutes with the Hamiltonian, A simple manipulation yields that Any conserved quantity K thus induces a linearized solution given by X K evaluated alongX (τ ). Clearly, the Wronskian of such two linearized trajectories is given by the Poisson bracket of the original conserved quantities, Following the previous discussion of the Jacobi system, a linearized solution X K associated with a conserved quantity K defines a conserved quantity k on the linearized phase space by the relation (A9) above, It is straightforward to show that cf. eqs. (A7) and (A2). It means, that k is the linearization of the original conserved quantity K, K(X (σ, τ )) = K(X (τ )) + σ k(N (τ )) + O(σ 2 ) , (A15) i.e., k(σN ) = K −K, cf. (4.23). Till now in this section, we have not assumed anything particular about the phase space and the Hamiltonian. We have just observed that any linearized solution defines through the Wronskian a conserved quantity (A9) and that any conserved quantity of the original system defines the linearized solution X K and the linearized conserved quantity k which is explicitly linear in N , as seen from (A14). We see that such a construction is completely general.
For a completely integrable system we have n mutually Poisson-commuting integrals of motion K i . They generate linearized solutions X K i . Thanks to the linearity of (A6) any linear (with constant coefficients) combination N of these solutions is again a linearized solution. Moreover, each of the conserved quantities K i induces the linearized quantity k i . Evaluating this linearized quantity on solutions X K j , one finds The same is true for any linear combination N of X K j , i.e., k i (N ) = 0. The conserved quantities K j thus directly generate a family of linearized solutions which all have vanishing values of linearized quantities k i . In other words, all these solutions have the same values of the conserved quantities K i as the central trajectory, K j =K j . It is not surprising since for the integrable system X K j generate symmetries of the evolution. These vector fields are tangent to the Lagrangian submanifolds given by K j = const. The linearized solutions X K j thus corresponds to trajectories which remain in the same Lagrangian submanifold as the central trajectory.
To conclude, we have just seen that the linearized conserved quantities constructed by using Wronskian and the linearized solutions of the equations of motion generated by the integrals of motion of the full system are general features of any Hamiltonian system admitting integrals of motion.

Cotangent bundle structure of the phase space
Let us now relate this general formulation to the configuration space description presented in sections III and IV. To that purpose we consider the phase space to be built from a configuration space M. While the configuration space is a space of "positions" x, the phase space is a space of "positions and momenta" (x, p). It is well known that such a phase space can be represented as a cotangent bundle T * M. The cotangent bundle has a natural symplectic structure Ω. If one chooses the configurationspace coordinates x a and the components p a of the momentum with respect to the frame dx a as coordinates in the phase space, X A ≡ (x a , p a ), the symplectic structure takes the canonical form, where the sum over spacetime index is naturally assumed. It means that (x a , p a ) are canonical coordinates in which the Poisson bracket takes the following form: A phase-space tangent vector N can be written with respect to coordinate frame ∂ ∂x a and ∂ ∂pa as Components n a can be understood as components of a configuration-space vector n, which is independent of the choice of coordinates x a . On other hand, componentsπ a cannot be combined to a 1-form π which would be independent of the choice of coordinates x a . Splitting (A19) to the configuration and momentum part is coordinate dependent. However, one can formulate a similar decomposition in a covariant way. As described in Appendix of [50], having a (torsion-free) covariant derivative ∇ on the configuration space, one can introduce the covariant splitting Here, the first term corresponds to a direction in the phase space given by the changing position x → x + εn with momentum p covariantly fixed (i.e., parallelly transported along n using the covariant derivative ∇). The second term corresponds to a direction in the phase space given by the changing momentum p → p + επ with x fixed. Such a splitting uniquely relates a phase-space vector N with a pair of the configuration-space vector n and the configuration-space 1-form π. Phase-space vectors ∇a ∂x and ∂ ∂pa thus describe horizontal and vertical directions in the cotangent bundle, where the 'horizontality' is given by the covariant derivative ∇.
Using this decomposition one can define in a covariant way derivatives of an observable F with respect to the position and the momentum, In terms of these, the Poisson bracket reads The linearized equations of motion can be written in the form (A5). It will be useful to write the Lie bracket of two phase-space vector fields N 1 , N 2 in terms of the covariant decomposition. In the splitting (A20) each vector field N i (x, p) corresponds to a pair n i (x, p), π i (x, p) of the configuration-space vector and 1-form, which both depend on x and p. Acting with the Lie bracket on a phase-space scalar F , one gets (A23) It can be shown that the first term is non-trivial while the next two terms vanish, see the end of this Appendix. The covariant splitting of the Lie bracket thus reads Let us now turn to the study of a geodesic motion in the configuration space. We assume the existence of a metric g ab and naturally choose ∇ to be the metric (torsion-free) covariant derivative, ∇ c g ab = 0. The geodesic motion in the configuration space is given by a simple quadratic Hamiltonian Since the metric is covariantly constant, we have ∇aH ∂x = 0 and ∂H ∂pa = g ab p b . Therefore, the covariant splitting of the Hamiltonian flow is Now, let us return to the linearized equations of motion (A6) near the central geodesicX (τ ), which is given in the configuration-space language byx(τ ) andp(τ ). The linearized trajectory N is given by (A20) and X H by (A27). Substituting these into (A6) and employing (A25), the linearized equations of motion yield Here, as before, the bar indicates quantities evaluated at the central geodesic. For the Hamiltonian (A26), the momentump of the central geodesic is proportional to the velocity,p = mū, and therefore 1 mp a ∇a ∂x ≡ ∇ dτ is the covariant time derivative along the central geodesic. We finally obtain (3.28), which gives the Jacobi equation for the linearized geodesics,

Commutation relations for derivatives on the cotangent bundle
To prove the relations (A24), it is useful to consider an observable F monomial in momentum, where in the first step we have used the fact that all r contributions from the Ricci identities are the same and in the second step employed the relation (A32). The remaining two commutators in (A24) are trivial. They reflect the fact that there is no curvature in p directions.
Appendix B: Covariant Lagrangian and Hamiltonian for the Jacobi system

Covariant expansion of the Lagrangian
The dynamics of a linearized system, i.e., the expansion of the equations of motion to the first order in the deviation variable, can be derived from the approximation of the Lagrangian to the second order. The first order contribution to the Lagrangian is trivial (given by a total derivative) since we are expanding near a solution of the equations of motion, i.e., near an extremal trajectory.
In order to expand the Lagrangian to the second order, we have to first generalize the notion of a deviation vector. In the main text, we have introduced the deviation vector n a as a tangent vector in σ-direction of the family of trajectories x(σ, τ ). This approach is sufficient for the first order approximation, e.g., for the derivation of the Jacobi equation (3.4). However, for an approximation to a higher order, and in particular to derive the Lagrangian (3.6), one has to first define the deviation vector more carefully.
We start again by choosing the central geodesicx(τ ) parameterized by its proper time τ and use this time as the external time for general trajectory x(τ ) nearby. Next we perform a time dependent change of the configuration variables, where instead of a configuration point x at time τ we use the deviation vector n a atx(τ ), n a = n a (x|x) . (B1) That is, we define n a (x|x) as a vector atx tangent to the geodesic joining pointsx and x, assuming that x is in a normal neighborhood ofx. To specify this vector uniquely, we normalize it to the geodesic distance of the corresponding geodesic segment. Intuitively, this vector plays a role of a difference between positions, n = x −x, generalized to the curved spacetime. We stress that the transformation (B1) is time dependent through the dependence on the pointx on the central geodesic. (For brevity, in what follows we will not write this τ -dependence explicitly.) This has 'unexpected consequences' for the linearized observables. In particular, as we shall see below, the covariant Hamiltonian for the Jacobi system will explicitly depend on time.
An advantage of the new variable n a is that it in oneto-one correspondence with the position x of the original system (at least, in the normal neighborhood ofx) while at the same time it is linear and belongs to a vector space. Therefore, one can perform an expansion in small n a and solve the system pertubatively.
The deviation vector n a , (B1), can be expressed as a derivative of the Synge world function 5 σ(x|y), [51]. The world function is defined as a half of the square of geodesic distance, with the sign given by the causal character of the geodesic, σ(x|y) = ± 1 2 (geodesic distance between x and y) 2 .
(B2) The deviation vector can be written as n a = n a (x|x) ≡ −∇ a σ(x|x) , see, e.g., [51][52][53]. Here we use the convention that ∇σ(x|x) denotes the derivative in the first argumentx and ∇σ(x|x) denotes the derivative in the second argument x. The normalization of n a is encoded in the relation When dealing with the Lagrangian, we also need a relation between velocities associated with variables x and n. For that we assume that all variables in (B3) are time dependent, x(τ ),x(τ ), and n a (τ ). Taking the covariant time derivative gives v a ≡ Dn a Dτ = −∇ a ∇ b σ(x|x) u b −ū b∇ b∇ a σ(x|x) . (B5) Let us notice, that in the coincidence limit x =x, u =ū, i.e., setting x and u to the values of the central trajectory, we have n = 0 and v = 0 (cf. (B7) below). The covariant expansion of the Lagrangian assumes that L is written as a series in n and v, and we ignore all terms of higher than second order, L(x, u) = L 0,0 + L 1,0 a n a + L 0,1 a v a + 1 2 L 2,0 ab n a n b + L 1,1 ab n a v b + 1 2 L 0,2 ab v a v b + . . . , with coefficients L k,l to be determined. For that, we need to take derivatives ∇ ∂x and ∂ ∂u of this relation, followed by the coincidence limit x =x, u =ū.
First, we calculate the coincidence limit x =x, u =ū of derivatives of relations for n(x) and v(x, u) given by (B3) and (B5). In this process we employ the following relations [51][52][53]: σ| x=x = 0 , ∇ a σ| x=x = ∇ a σ| x=x = 0 , This yields the following: ∇ a n k ∂x x=x = δ k a , ∂v k ∂u a x=x u=ū = δ k a ,  If we use that the central trajectory satisfies the Euler-Lagrange equations the definition of the canonical momentump, and express the velocity as a derivative we get L −L = Dp a Dτ n a +p a Dn a Dτ = d dτ p a n a .
Thus, for an arbitrary Lagrangian the absolute and linear terms can be written as a total derivative term d dτ f (n) with f (n) =S + n ap a , Finally, the Hamiltonian h, (B20), for the linearized system gives h = H + ∂G ∂τ = 1 2mḡ ab π a π b + m 2ū kūlR kalb n a n b + . . . , (B27) which is the Hamiltonian (3.27).
To summarize, due to the time dependence of the canonical transformation (B1), and due to time-dependent shift (B15) of the Lagrangian, the new Hamiltonian h for the linearized system is not conserved, despite the fact that the original Hamiltonian H is a conserved quantity along the geodesic. Clearly,