kerrgeodesic_gw
¶This Jupyter notebook requires SageMath (version $\geq$ 8.2), with the package kerrgeodesic_gw (version $\geq$ 0.3). To install the latter, simply run
sage -pip install kerrgeodesic_gw
in a terminal.
version()
'SageMath version 9.4, Release Date: 2021-08-22'
First, we set up the notebook to use LaTeX-formatted display:
%display latex
and we ask for CPU demanding computations to be performed in parallel on 8 processes:
Parallelism().set(nproc=8)
A Kerr black bole is entirely defined by two parameters $(m, a)$, where $m$ is the black hole mass and $a$ is the black hole angular momentum divided by $m$.
In this notebook, we shall set $m=1$ and we denote the angular momentum parameter $a$ by the symbolic variable a
, using a0
for a specific numerical value:
a = var('a')
a0 = 0.998
The spacetime object is created as an instance of the class KerrBH
:
from kerrgeodesic_gw import KerrBH
M = KerrBH(a)
print(M)
Kerr spacetime M
The object M
is endowed with many methods, which can be discovered via the TAB key:
# M.<TAB>
One of them returns the Boyer-Lindquist coordinate $r$ of the event horizon:
rH = M.event_horizon_radius()
rH
rH0 = rH.subs({a: a0})
rH0
Another one returns the chart of Boyer-Lindquist coordinates and allows the user to instanciate the Python variables (t, r, th, ph)
to the coordinates $(t,r,\theta,\phi)$:
BL.<t, r, th, ph> = M.boyer_lindquist_coordinates()
BL
The metric tensor is naturally returned by the method metric()
:
g = M.metric()
g.display()
Again many methods are available on metric objects and one can discover them via the TAB key:
# g.<TAB>
For instance christoffel_symbols_display()
computes the Christofell symbols with respect to the default chart (Boyer-Lindquist) and displays them:
g.christoffel_symbols_display()
The Riemann curvature tensor is naturally returned by the method riemann()
Riem = g.riemann()
print(Riem)
Tensor field Riem(g) of type (1,3) on the Kerr spacetime M
The $R^t_{rtr}$ and $R^t_{r\theta\phi}$ components:
Riem[0,1,0,1]
Riem[0,1,2,3]
Of course, since the Kerr metric is a solution of the vacuum Einstein equation, the Ricci tensor identically vanishes:
g.ricci().display()
Let us choose the initial point $P$ for the geodesic:
P = M.point((0, 6, pi/2, 0), name='P')
print(P)
Point P on the Kerr spacetime M
A geodesic is constructed by providing the range $[\lambda_{\rm min},\lambda_{\rm max}]$ of the affine parameter $\lambda$, the initial point and either
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, 300\, m]$, the option (ii) and $a=0.998 \,m$, where $m$ in the black hole mass::
Li = M.geodesic([0, 300], P, mu=1, E=0.883, L=1.982, Q=0.467, a_num=0.998,
name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector:
The curve was correctly set. Parameters appearing in the differential system defining the curve are [a].
print(Li)
Geodesic Li of the Kerr spacetime M
The numerical integration of the geodesic equation is performed via integrate()
, by providing the integration step $\delta\lambda$ in units of $m$:
Li.integrate(step=0.005)
We can then plot the geodesic:
Li.plot()
Actually, many options can be passed to the method plot()
. For instance to a get a 3D spacetime diagram:
Li.plot(coordinates='txy')
or to get the trace of the geodesic in the $(x,y)$ plane:
Li.plot(coordinates='xy', plot_points=2000)
or else to get the trace in the $(x,z)$ plane:
Li.plot(coordinates='xz')
As a curve, the geodesic $\mathcal{L}$ is a map from an interval of $\mathbb{R}$ to the spacetime $M$:
Li.display()
Li.domain()
Li.codomain()
It maps values of $\lambda$ to spacetime points:
Li(0)
Li(0).coordinates() # coordinates in the default chart (BL)
BL(Li(0)) # equivalent to above
Li(300).coordinates()
The initial 4-momentum vector $p_0$ is returned by the method initial_tangent_vector()
:
p0 = Li.initial_tangent_vector()
print(p0)
Tangent vector p at Point P on the Kerr spacetime M
p0 in M.tangent_space(P)
p0.display()
p0[:]
For instance, the components $p^t_0$ and $p^\phi_0$ are recovered by
p0[0], p0[3]
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:
g = M.metric()
g.at(P)(p0, p0).subs(a=0.998)
The 4-momentum vector $p$ at any value of the affine parameter $\lambda$ is obtained by
the method evaluate_tangent_vector()
; for instance, for $\lambda=200\,m$:
p = Li.evaluate_tangent_vector(200)
p.display()
p in M.tangent_space(Li(200))
The particle mass $\mu$ computed at a given value of $\lambda$ is returned by the method evaluate_mu()
:
Li.evaluate_mu(0)
Of course, it should be conserved along $\mathcal{L}$; actually it is, up to the numerical accuracy::
Li.evaluate_mu(300)
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()
:
Li.evaluate_E(0)
Li.evaluate_L(0)
Li.evaluate_Q(0)
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:
Li.check_integrals_of_motion(300)
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:
Li.integrate(step=0.001)
Li.check_integrals_of_motion(300)
quantity | value | initial value | diff. | relative diff. |
\(\mu^2\) | \(1.0000047183936422\) | \(1.00000000000000\) | \(4.718 \times 10^{-6}\) | \(4.718 \times 10^{-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\) |
Instead of providing the integral of motions, as for Li
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:
print(p0)
Tangent vector p at Point P on the Kerr spacetime M
p0[:]
Li2 = M.geodesic([0, 300], P, pt0=p0[0], pr0=p0[1], pth0=p0[2], pph0=p0[3],
a_num=0.998)
Li2.initial_tangent_vector() == p0
As a check, we recover the same values of $(\mu, E, L, Q)$ as those that were used to initialize Li
:
Li2.evaluate_mu(0)
Li2.evaluate_E(0)
Li2.evaluate_L(0)
Li2.evaluate_Q(0)
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:
Li3 = M.geodesic([0, 300], P, mu=1, pr0=p0[1], pth0=p0[2], pph0=p0[3],
a_num=0.998)
The component $p^t_0$ is then automatically computed:
Li3.initial_tangent_vector()[:]
and we check the identity with the initial vector of Li
, up to numerical errors:
(Li3.initial_tangent_vector() - p0)[:]
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:
Li4 = M.geodesic([0, 300], P, E=0.8830, L=1.982, pr0=p0[1], pth0=p0[2],
a_num=0.998)
Li4.initial_tangent_vector()[:]
Again, we get a geodesic equivalent to Li
:
(Li4.initial_tangent_vector() - p0)[:]
We choose a ingoing null geodesic in the equatorial plane with $L = -6 E < 0$, starting at the point of Boyer-Lindquist coordinates $(t,r,\theta,\phi) = (0, 12, \pi/2, 0)$:
lambda_max = 13.063
Li = M.geodesic([0, lambda_max], M((0,12,pi/2,0)), mu=0, E=1, L=-6, Q=0,
r_increase=False, a_num=a0, name='Li', latex_name=r'\mathcal{L}',
verbose=True)
Initial tangent vector:
The curve was correctly set. Parameters appearing in the differential system defining the curve are [a].
Li.integrate(step=0.00002)
We plot the trajectory of the geodesic in the equatorial plane with the tangent vector at 6 points:
Li.plot(coordinates='xy', color='green', prange=[0, lambda_max], plot_points=40000,
display_tangent=True, scale=0.002, color_tangent='brown',
plot_points_tangent=6, axes_labels=[r'$x/m$', r'$y/m$'])
Note the turning point in $\phi$ and the final winding in the direction of the black rotation (counterclockwise in the figure).
Li.check_integrals_of_motion(0.99999*lambda_max)
quantity | value | initial value | diff. | relative diff. |
\(\mu^2\) | \(0.000211815000511706\) | \(0.000000000000000\) | \(0.0002118\) | - |
\(E\) | \(1.00000186057173\) | \(1.00000000000000\) | \(1.861 \times 10^{-6}\) | \(1.861 \times 10^{-6}\) |
\(L\) | \(-5.99999596865291\) | \(-6.00000000000000\) | \(4.031 \times 10^{-6}\) | \(-6.719 \times 10^{-7}\) |
\(Q\) | \(1.31244559310830 \times 10^{-31}\) | \(0\) | \(1.312 \times 10^{-31}\) | - |