# Geodesics in Kerr spacetime¶

Geodesics of Kerr spacetime are implemented via the class KerrGeodesic.

class kerrgeodesic_gw.kerr_geodesic.KerrGeodesic(parent, initial_point, pt0=None, pr0=None, pth0=None, pph0=None, mu=None, E=None, L=None, Q=None, r_increase=True, th_increase=True, chart=None, name=None, latex_name=None, a_num=None, m_num=None, verbose=False)

Bases: sage.manifolds.differentiable.integrated_curve.IntegratedGeodesic

Geodesic of Kerr spacetime.

Geodesics are computed by solving the geodesic equation via the generic SageMath geodesic integrator.

INPUT:

• parentIntegratedGeodesicSet, the set of curves $$\mathrm{Hom_{geodesic}}(I, M)$$ to which the geodesic belongs

• initial_point – point of Kerr spacetime from which the geodesic is to be integrated

• pt0 – (default: None) Boyer-Lindquist component $$p^t$$ of the initial 4-momentum vector

• pr0 – (default: None) Boyer-Lindquist component $$p^r$$ of the initial 4-momentum vector

• pth0 – (default: None) Boyer-Lindquist component $$p^\theta$$ of the initial 4-momentum vector

• pph0 – (default: None) Boyer-Lindquist component $$p^\phi$$ of the initial 4-momentum vector

• mu – (default: None) mass $$\mu$$ of the particle

• E – (default: None) conserved energy $$E$$ of the particle

• L – (default: None) conserved angular momemtum $$L$$ of the particle

• Q – (default: None) Carter constant $$Q$$ of the particle

• r_increase – (default: True) boolean; if True, the initial value of $$p^r=\mathrm{d}r/\mathrm{d}\lambda$$ determined from the integral of motions is positive or zero, otherwise, $$p^r$$ is negative

• th_increase – (default: True) boolean; if True, the initial value of $$p^\theta=\mathrm{d}\theta/\mathrm{d}\lambda$$ determined from the integral of motions is positive or zero, otherwise, $$p^\theta$$ is negative

• chart – (default: None) chart on the spacetime manifold in terms of which the geodesic equations are expressed; if None the default chart (Boyer-Lindquist coordinates) is assumed

• name – (default: None) string; symbol given to the geodesic

• latex_name – (default: None) string; LaTeX symbol to denote the geodesic; if none is provided, name will be used

• a_num – (default: None) numerical value of the Kerr spin parameter $$a$$ (required for a numerical integration)

• m_num – (default: None) numerical value of the Kerr mass parameter $$m$$ (required for a numerical integration)

• verbose – (default: False) boolean; determines whether some information is printed during the construction of the geodesic

EXAMPLES:

We construct first the Kerr spacetime:

sage: from kerrgeodesic_gw import KerrBH
sage: a = var('a')
sage: M = KerrBH(a); M
Kerr spacetime M
sage: BLchart = M.boyer_lindquist_coordinates(); BLchart
Chart (M, (t, r, th, ph))


We pick an initial spacetime point for the geodesic:

sage: init_point = M((0, 6, pi/2, 0), name='P')


A geodesic is constructed by providing the range of the affine parameter, the initial point and either (i) the Boyer-Lindquist components $$(p^t_0, p^r_0, p^\theta_0, p^\phi_0)$$ of the initial 4-momentum vector $$p_0 = \left. \frac{\mathrm{d}x}{\mathrm{d}\lambda}\right| _{\lambda=0}$$, (ii) the four integral of motions $$(\mu, E, L, Q)$$ or (iii) some of the components of $$p_0$$ along with with some integrals of motion. We shall also specify some numerical value for the Kerr spin parameter $$a$$. Examples of (i) and (iii) are provided below. Here, we choose $$\lambda\in[0, 300m]$$, the option (ii) and $$a=0.998 m$$, where $$m$$ in the black hole mass:

sage: geod = M.geodesic([0, 300], init_point, mu=1, E=0.883,
....:                   L=1.982, Q=0.467, a_num=0.998)
sage: geod
Geodesic of the Kerr spacetime M


The numerical integration of the geodesic equation is performed via integrate(), by providing the step in $$\delta\lambda$$ in units of $$m$$:

sage: geod.integrate(step=0.005)


We can then plot the geodesic:

sage: geod.plot()
Graphics3d Object


Actually, many options can be passed to plot(). For instance to a get a 3D spacetime diagram:

sage: geod.plot(coordinates='txy')
Graphics3d Object


or to get the trace of the geodesic in the $$(x,y)$$ plane:

sage: geod.plot(coordinates='xy', plot_points=2000)
Graphics object consisting of 2 graphics primitives


or else to get the trace in the $$(x,z)$$ plane:

sage: geod.plot(coordinates='xz')
Graphics object consisting of 2 graphics primitives


As a curve, the geodesic is a map from an interval of $$\mathbb{R}$$ to the spacetime $$M$$:

sage: geod.display()
(0, 300) → M
sage: geod.domain()
Real interval (0, 300)
sage: geod.codomain()
Kerr spacetime M


It maps values of $$\lambda$$ to spacetime points:

sage: geod(0)
Point on the Kerr spacetime M
sage: geod(0).coordinates()  # coordinates in the default chart  # tol 1.0e-13
(0.0, 6.0, 1.5707963267948966, 0.0)
sage: BLchart(geod(0))       # equivalent to above   # tol 1.0e-13
(0.0, 6.0, 1.5707963267948966, 0.0)
sage: geod(300).coordinates()   # tol 1.0e-13
(553.4637326813786, 3.703552505462962, 1.6613834863942039, 84.62814710987239)


The initial 4-momentum vector $$p_0$$ is returned by the method initial_tangent_vector():

sage: p0 = geod.initial_tangent_vector(); p0
Tangent vector p at Point P on the Kerr spacetime M
sage: p0 in M.tangent_space(init_point)
True
sage: p0.display()  # tol 1.0e-13
p = 1.29225788954106 ∂/∂t + 0.00438084990626460 ∂/∂r
+ 0.0189826106258554 ∂/∂th + 0.0646134478134985 ∂/∂ph
sage: p0[:]  # tol 1.0e-13
[1.29225788954106, 0.00438084990626460, 0.0189826106258554, 0.0646134478134985]


For instance, the components $$p^t_0$$ and $$p^\phi_0$$ are recovered by:

sage: p0[0], p0[3]  # tol 1.0e-13
(1.29225788954106, 0.0646134478134985)


Let us check that the scalar square of $$p_0$$ is $$-1$$, i.e. is consistent with the mass parameter $$\mu = 1$$ used in the construction of the geodesic:

sage: g = M.metric()
sage: g.at(init_point)(p0, p0).subs(a=0.998)  # tol 1.0e-13
-1.00000000000000


The 4-momentum vector $$p$$ at any value of the affine parameter $$\lambda$$, e.g. $$\lambda=200m$$, is obtained by:

sage: p = geod.evaluate_tangent_vector(200); p
Tangent vector at Point on the Kerr spacetime M
sage: p in M.tangent_space(geod(200))
True
sage: p.display()  # tol 1.0e-13
1.316592599498746 ∂/∂t - 0.07370434215844164 ∂/∂r
- 0.01091195426423706 ∂/∂th + 0.07600209768075264 ∂/∂ph


The particle mass $$\mu$$ computed at a given value of $$\lambda$$ is returned by the method evaluate_mu():

sage: geod.evaluate_mu(0)  # tol 1.0e-13
1.00000000000000


Of course, it should be conserved along the geodesic; actually it is, up to the numerical accuracy:

sage: geod.evaluate_mu(300)  # tol 1.0e-13
1.0000117978600134


Similarly, the conserved energy $$E$$, conserved angular momentum $$L$$ and Carter constant $$Q$$ are computed at any value of $$\lambda$$ by respectively evaluate_E(), evaluate_L() and evaluate_Q():

sage: geod.evaluate_E(0)  # tol 1.0e-13
0.883000000000000
sage: geod.evaluate_L(0)  # tol 1.0e-13
1.98200000000000
sage: geod.evaluate_Q(0)  # tol 1.0e-13
0.467000000000000


Let us check that the values of $$\mu$$, $$E$$, $$L$$ and $$Q$$ evaluated at $$\lambda=300 m$$ are equal to those at $$\lambda=0$$ up to the numerical accuracy of the integration scheme:

sage: geod.check_integrals_of_motion(300)  # tol 1.0e-13
quantity         value            initial value       diff.      relative diff.
$\mu^2$    1.0000235958592163   1.00000000000000    0.00002360     0.00002360
$E$      0.883067996080701    0.883000000000000   0.00006800     0.00007701
$L$       1.98248080818931    1.98200000000000    0.0004808      0.0002426
$Q$      0.467214137649741    0.467000000000000   0.0002141      0.0004585


Decreasing the integration step leads to smaller errors:

sage: geod.integrate(step=0.001)
sage: geod.check_integrals_of_motion(300)  # tol 1.0e-13
quantity         value            initial value       diff.      relative diff.
$\mu^2$    1.0000047183936422   1.00000000000000     4.718e-6       4.718e-6
$E$      0.883013604456676    0.883000000000000   0.00001360     0.00001541
$L$      1.98209626120918     1.98200000000000    0.00009626     0.00004857
$Q$      0.467042771975860    0.467000000000000   0.00004277     0.00009159


Various ways to initialize a geodesic

Instead of providing the integral of motions, as for geod above, one can initialize a geodesic by providing the Boyer-Lindquist components $$(p^t_0, p^r_0, p^\theta_0, p^\phi_0)$$ of the initial 4-momentum vector $$p_0$$. For instance:

sage: p0
Tangent vector p at Point P on the Kerr spacetime M
sage: p0[:]  # tol 1.0e-13
[1.29225788954106, 0.00438084990626460, 0.0189826106258554, 0.0646134478134985]
sage: geod2 = M.geodesic([0, 300], init_point, pt0=p0[0], pr0=p0[1],
....:                    pth0=p0[2], pph0=p0[3], a_num=0.998)
sage: geod2.initial_tangent_vector() == p0
True


As a check, we recover the same values of $$(\mu, E, L, Q)$$ as those that were used to initialize geod:

sage: geod2.evaluate_mu(0)
1.00000000000000
sage: geod2.evaluate_E(0)
0.883000000000000
sage: geod2.evaluate_L(0)
1.98200000000000
sage: geod2.evaluate_Q(0)
0.467000000000000


We may also initialize a geodesic by providing the mass $$\mu$$ and the three spatial components $$(p^r_0, p^\theta_0, p^\phi_0)$$ of the initial 4-momentum vector:

sage: geod3 = M.geodesic([0, 300], init_point, mu=1, pr0=p0[1],
....:                    pth0=p0[2], pph0=p0[3], a_num=0.998)


The component $$p^t_0$$ is then automatically computed:

sage: geod3.initial_tangent_vector()[:]  # tol 1.0e-13
[1.29225788954106, 0.00438084990626460, 0.0189826106258554, 0.0646134478134985]


and we check the identity with the initial vector of geod, up to numerical errors:

sage: (geod3.initial_tangent_vector() - p0)[:]  # tol 1.0e-13
[2.22044604925031e-16, 0, 0, 0]


Another way to initialize a geodesic is to provide the conserved energy $$E$$, the conserved angular momentum $$L$$ and the two components $$(p^r_0, p^\theta_0)$$ of the initial 4-momentum vector:

sage: geod4 = M.geodesic([0, 300], init_point, E=0.8830, L=1.982,
....:                     pr0=p0[1], pth0=p0[2], a_num=0.998)
sage: geod4.initial_tangent_vector()[:]
[1.29225788954106, 0.00438084990626460, 0.0189826106258554, 0.0646134478134985]


Again, we get a geodesic equivalent to geod:

sage: (geod4.initial_tangent_vector() - p0)[:]  # tol 1.0e-13
[0, 0, 0, 0]

check_integrals_of_motion(affine_parameter, solution_key=None)

Check the constancy of the four integrals of motion

INPUT:

• affine_parameter – value of the affine parameter $$\lambda$$

• solution_key – (default: None) string denoting the numerical solution to use for evaluating the various integrals of motion; if None, the latest solution computed by integrate() is used.

OUTPUT:

• a SageMath table with the absolute and relative differences with respect to the initial values.

evaluate_E(affine_parameter, solution_key=None)

Compute the conserved energy $$E$$ at a given value of the affine parameter $$\lambda$$.

INPUT:

• affine_parameter – value of the affine parameter $$\lambda$$

• solution_key – (default: None) string denoting the numerical solution to use for the evaluation; if None, the latest solution computed by integrate() is used.

OUTPUT:

• value of $$E$$

evaluate_L(affine_parameter, solution_key=None)

Compute the conserved angular momentum about the rotation axis $$L$$ at a given value of the affine parameter $$\lambda$$.

INPUT:

• affine_parameter – value of the affine parameter $$\lambda$$

• solution_key – (default: None) string denoting the numerical solution to use for the evaluation; if None, the latest solution computed by integrate() is used.

OUTPUT:

• value of $$L$$

evaluate_Q(affine_parameter, solution_key=None)

Compute the Carter constant $$Q$$ at a given value of the affine parameter $$\lambda$$.

INPUT:

• affine_parameter – value of the affine parameter $$\lambda$$

• solution_key – (default: None) string denoting the numerical solution to use for the evaluation; if None, the latest solution computed by integrate() is used.

OUTPUT:

• value of $$Q$$

evaluate_mu(affine_parameter, solution_key=None)

Compute the mass parameter $$\mu$$ at a given value of the affine parameter $$\lambda$$.

INPUT:

• affine_parameter – value of the affine parameter $$\lambda$$

• solution_key – (default: None) string denoting the numerical solution to use for the evaluation; if None, the latest solution computed by integrate() is used.

OUTPUT:

• value of $$\mu$$

evaluate_mu2(affine_parameter, solution_key=None)

Compute the square of the mass parameter $$\mu^2$$ at a given value of the affine parameter $$\lambda$$.

INPUT:

• affine_parameter – value of the affine parameter $$\lambda$$

• solution_key – (default: None) string denoting the numerical solution to use for the evaluation; if None, the latest solution computed by integrate() is used.

OUTPUT:

• value of $$\mu^2$$

evaluate_tangent_vector(affine_parameter, solution_key=None)

Return the tangent vector (4-momentum) at a given value of the affine parameter.

INPUT:

• affine_parameter – value of the affine parameter $$\lambda$$

• solution_key – (default: None) string denoting the numerical solution to use for the evaluation; if None, the latest solution computed by integrate() is used.

OUTPUT:

• instance of TangentVector representing the 4-momentum vector $$p$$ at $$\lambda$$

initial_tangent_vector()

Return the initial 4-momentum vector.

OUTPUT:

• instance of TangentVector representing the initial 4-momentum $$p_0$$

EXAMPLES:

Initial 4-momentum vector of a null geodesic in Schwarzschild spacetime:

sage: from kerrgeodesic_gw import KerrBH
sage: M = KerrBH(0)
sage: BLchart = M.boyer_lindquist_coordinates()
sage: init_point = M((0, 6, pi/2, 0), name='P')
sage: geod = M.geodesic([0, 100], init_point, mu=0, E=1,
....:                   L=3, Q=0)
sage: p0 = geod.initial_tangent_vector(); p0
Tangent vector p at Point P on the Schwarzschild spacetime M
sage: p0.display()
p = 3/2 ∂/∂t + 1/6*sqrt(30) ∂/∂r + 1/12 ∂/∂ph

integrate(step=None, method='odeint', solution_key=None, parameters_values=None, verbose=False, **control_param)

Solve numerically the geodesic equation.

INPUT:

• step – (default: None) step $$\delta\lambda$$ for the integration, where $$\lambda$$ is the affine parameter along the geodesic; default value is a hundredth of the range of $$\lambda$$ declared when constructing the geodesic

• method – (default: 'odeint') numerical scheme to use for the integration; available algorithms are:

• 'odeint' - makes use of scipy.integrate.odeint via Sage solver desolve_odeint(); odeint invokes the LSODA algorithm of the ODEPACK suite, which automatically selects between implicit Adams method (for non-stiff problems) and a method based on backward differentiation formulas (BDF) (for stiff problems).

• 'rk4_maxima' - 4th order classical Runge-Kutta, which makes use of Maxima’s dynamics package via Sage solver desolve_system_rk4() (quite slow)

• 'dopri5' - Dormand-Prince Runge-Kutta of order (4)5 provided by scipy.integrate.ode

• 'dop853' - Dormand-Prince Runge-Kutta of order 8(5,3) provided by scipy.integrate.ode

and those provided by GSL via Sage class ode_solver:

• 'rk2' - embedded Runge-Kutta (2,3)

• 'rk4' - 4th order classical Runge-Kutta

• 'rkf45' - Runge-Kutta-Felhberg (4,5)

• 'rkck' - embedded Runge-Kutta-Cash-Karp (4,5)

• 'rk8pd' - Runge-Kutta Prince-Dormand (8,9)

• 'rk2imp' - implicit 2nd order Runge-Kutta at Gaussian points

• 'rk4imp' - implicit 4th order Runge-Kutta at Gaussian points

• 'gear1' - $$M=1$$ implicit Gear

• 'gear2' - $$M=2$$ implicit Gear

• 'bsimp' - implicit Bulirsch-Stoer (requires Jacobian)

• solution_key – (default: None) string to tag the numerical solution; if None, the string method is used.

• parameters_values – (default: None) list of numerical values of the parameters present in the system defining the geodesic, to be substituted in the equations before integration

• verbose – (default: False) prints information about the computation in progress

• **control_param – extra control parameters to be passed to the chosen solver

EXAMPLES:

Bound timelike geodesic in Schwarzschild spacetime:

sage: from kerrgeodesic_gw import KerrBH
sage: M = KerrBH(0)
sage: BLchart = M.boyer_lindquist_coordinates()
sage: init_point = M((0, 10, pi/2, 0), name='P')
sage: lmax = 1500.
sage: geod = M.geodesic([0, lmax], init_point, mu=1, E=0.973,
....:                   L=4.2, Q=0)
sage: geod.integrate()
sage: geod.plot(coordinates='xy')
Graphics object consisting of 2 graphics primitives


With the default integration step, the accuracy is not very good:

sage: geod.check_integrals_of_motion(lmax)  # tol 1.0e-13
quantity           value             initial value       diff.     relative diff.
$\mu^2$      1.000761704316941     1.00000000000000    0.0007617     0.0007617
$E$       0.9645485805304451     0.973000000000000   -0.008451     -0.008686
$L$       3.8897905637080923     4.20000000000000     -0.3102       -0.07386
$Q$      5.673017835722329e-32           0           5.673e-32         -


Let us improve it by specifying a smaller integration step:

sage: geod.integrate(step=0.1)
sage: geod.check_integrals_of_motion(lmax)  # tol 1.0e-13
quantity           value             initial value        diff.      relative diff.
$\mu^2$     1.0000101879128696     1.00000000000000    0.00001019      0.00001019
$E$       0.9729448260004574     0.973000000000000   -0.00005517    -0.00005671
$L$        4.197973829219027     4.20000000000000     -0.002026      -0.0004824
$Q$      6.607560764960032e-32           0            6.608e-32          -


We may set the parameter solution_key to keep track of various numerical solutions:

sage: geod.integrate(step=0.1, solution_key='step_0.1')
sage: geod.integrate(step=0.02, solution_key='step_0.02')


and use it in the various evaluation functions:

sage: geod.evaluate_mu(lmax, solution_key='step_0.1')  # tol 1.0e-13
1.0000050939434606
sage: geod.evaluate_mu(lmax, solution_key='step_0.02')  # tol 1.0e-13
1.0000010212056811

plot(coordinates='xyz', prange=None, solution_key=None, style='-', thickness=1, plot_points=1000, color='red', include_end_point=(True, True), end_point_offset=(0.001, 0.001), verbose=False, label_axes=None, plot_horizon=True, horizon_color='black', fill_BH_region=True, BH_region_color='grey', display_tangent=False, color_tangent='blue', plot_points_tangent=10, width_tangent=1, scale=1, aspect_ratio='automatic', **kwds)

Plot the geodesic in terms of the coordinates $$(t,x,y,z)$$ deduced from the Boyer-Lindquist coordinates $$(t,r,\theta,\phi)$$ via the standard polar to Cartesian transformations.

INPUT:

• coordinates – (default: 'xyz') string indicating which of the coordinates $$(t,x,y,z)$$ to use for the plot

• prange – (default: None) range of the affine parameter $$\lambda$$ for the plot; if None, the entire range declared during the construction of the geodesic is considered

• solution_key – (default: None) string denoting the numerical solution to use for the plot; if None, the latest solution computed by integrate() is used.

• verbose – (default: False) determines whether information is printed about the plotting in progress

• plot_horizon – (default: True) determines whether the black hole event horizon is drawn

• horizon_color – (default: 'black') color of the event horizon

• fill_BH_region – (default: True) determines whether the black hole region is colored (for 2D plots only)

• BH_region_color – (default: 'grey') color of the event horizon

• display_tangent – (default: False) determines whether some tangent vectors should also be plotted

• color_tangent – (default: blue) color of the tangent vectors when these are plotted

• plot_points_tangent – (default: 10) number of tangent vectors to display when these are plotted

• width_tangent – (default: 1) sets the width of the arrows representing the tangent vectors

• scale – (default: 1) scale applied to the tangent vectors before displaying them

• either a 2D graphic oject (2 coordinates specified in the parameter coordinates) or a 3D graphic object (3 coordinates in coordinates)