#!/usr/bin/env python # coding: utf-8 # # Spherical null geodesics in Kerr spacetime # ## Computation with `kerrgeodesic_gw` # # This Jupyter/SageMath notebook is relative to the lectures # [Geometry and physics of black holes](https://relativite.obspm.fr/blackholes/). # # It requires [SageMath](http://www.sagemath.org/) (version $\geq$ 8.2), with the package [kerrgeodesic_gw](https://github.com/BlackHolePerturbationToolkit/kerrgeodesic_gw) (version $\geq$ 0.3.2). To install the latter, simply run # ``` # sage -pip install kerrgeodesic_gw # ``` # in a terminal. # In[1]: version() # First, we set up the notebook to use LaTeX-formatted display: # In[2]: get_ipython().run_line_magic('display', 'latex') # and we ask for CPU demanding computations to be performed in parallel on 8 processes: # In[3]: 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: # In[4]: a = var('a') a0 = 0.95 # The spacetime object is created as an instance of the class `KerrBH`: # In[5]: from kerrgeodesic_gw import KerrBH M = KerrBH(a) print(M) # The Boyer-Lindquist coordinate $r$ of the event horizon: # In[6]: rH = M.event_horizon_radius() rH # In[7]: rH0 = rH.subs({a: a0}) show(LatexExpr(r'r_+ = '), rH0) # In[8]: show(LatexExpr(r'r_- = '), M.inner_horizon_radius().subs({a: a0})) # The method `boyer_lindquist_coordinates()` returns the chart of Boyer-Lindquist coordinates `BL` and allows the user to instanciate the Python variables `(t, r, th, ph)` to the coordinates $(t,r,\theta,\phi)$: # In[9]: BL. = M.boyer_lindquist_coordinates() BL # The metric tensor is naturally returned by the method `metric()`: # In[10]: g = M.metric() g.display() # ## Spherical photon orbits # Functions $\ell(r_0)$ and $q(r_0)$ for spherical photon orbits: # In[11]: r = var('r') lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1)) lsph # In[12]: qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2) qsph # $\theta$-turning points: # In[13]: theta0(a, l, q) = arccos(sqrt(1/2*(1 - (l^2+q)/a^2 + sqrt((1 - (l^2+q)/a^2)^2 + 4*q/a^2)))) theta0 # In[14]: theta1(a, l, q) = arccos(sqrt(1/2*(1 - (l^2+q)/a^2 - sqrt((1 - (l^2+q)/a^2)^2 + 4*q/a^2)))) theta1 # ### Spherical photon orbit at $r_0 = 3 m$ ($q = q_{\rm max} = 27 m^2$) # In[15]: r0 = 3. E = 1 L = lsph(a0, r0) Q = qsph(a0, r0) L, Q # In[16]: theta0(a0, L, Q) # In[17]: P = M.point((0, r0, pi/2, 0), name='P') print(P) # A geodesic is constructed by providing the range $[\lambda_{\rm min},\lambda_{\rm max}]$ of the affine parameter $\lambda$, 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_{\rm min}}$, # - (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. # In[18]: lmax = 100 # lambda_max # In[19]: Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[20]: v0 = Li.initial_tangent_vector() Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3], a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # The numerical integration of the geodesic equation is performed via `integrate()`, by providing the integration step $\delta\lambda$: # In[21]: Li.integrate(step=0.001, method='dopri5') Li.check_integrals_of_motion(0.999*lmax) # In[22]: print("Final point: ") show(BL[:], "=", BL(Li(0.999*lmax))) # In[23]: lplot = 0.32*lmax print("max lambda (plot): ", lplot) Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.15, width_tangent=2, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey') \ + line([(0,0,-3), (0,0,3)], color='black', thickness=2) # In[24]: graph = Li.plot(coordinates='xy', prange=(0, lplot), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.15, width_tangent=2, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey', axes_labels=[r'$x/m$', r'$y/m$']) graph.save("gik_spher_3d_r_30_xy.pdf") graph # In[25]: Li.plot(coordinates='xz', prange=(0, lplot), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.15, width_tangent=2, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey', axes_labels=[r'$x/m$', r'$z/m$']) # ### Prograde spherical photon orbit at $r_0 = 1.6m$ # In[26]: r0 = 1.6 L = lsph(a0, r0) Q = qsph(a0, r0) L, Q # In[27]: theta0(a0, L, Q) # In[28]: P = M.point((0, r0, pi/2, 0), name='P') lmax = 70 Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[29]: v0 = Li.initial_tangent_vector() Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3], a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[30]: Li.integrate(step=0.0004, method='dopri5') Li.check_integrals_of_motion(0.999*lmax) # In[31]: print("Final point: ") show(BL[:], "=", BL(Li(0.999*lmax))) # In[32]: lplot = 0.11*lmax print("max lambda (plot): ", lplot) Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.03, width_tangent=1, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey') \ + line([(0,0,-1.5), (0,0,1.5)], color='black', thickness=2) # In[33]: Li.plot(prange=(0, lmax), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.03, width_tangent=1, color_tangent='green', plot_points_tangent=40, horizon_color='lightgrey') \ + line([(0,0,-1.5), (0,0,1.5)], color='black', thickness=2) # ### Spherical photon orbit at $r_0=2.8m$ # In[34]: r0 = 2.8 L = lsph(a0, r0) Q = qsph(a0, r0) L, Q # In[35]: theta0(a0, L, Q) # In[36]: P = M.point((0, r0, pi/2, 0), name='P') lmax = 100 Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[37]: v0 = Li.initial_tangent_vector() Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3], a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[38]: Li.integrate(step=0.001, method='dopri5') Li.check_integrals_of_motion(0.999*lmax) # In[39]: print("Final point: ") show(BL[:], "=", BL(Li(0.999*lmax))) # In[40]: lplot = 0.38*lmax print("max lambda (plot): ", lplot) Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.12, width_tangent=2, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey') # In[41]: graph = Li.plot(coordinates='xy', prange=(0, lplot), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.01, width_tangent=1, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey', axes_labels=[r'$x/m$', r'$y/m$']) graph.save("gik_spher_3d_r_28_xy.pdf") graph # In[42]: Li.plot(coordinates='txy', prange=(0, lplot), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.2, width_tangent=2, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey') # ### Retrograde spherical photon orbit at $r_0 = 3.9 m$ # In[43]: r0 = 3.9 L = lsph(a0, r0) Q = qsph(a0, r0) L, Q # In[44]: theta0(a0, L, Q) # In[45]: P = M.point((0, r0, pi/2, 0), name='P') # In[46]: lmax = 100 Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[47]: v0 = Li.initial_tangent_vector() Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3], a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[48]: Li.integrate(step=0.001, method='dopri5') Li.check_integrals_of_motion(0.999*lmax) # In[49]: print("Final point: ") show(BL[:], "=", BL(Li(0.999*lmax))) # In[50]: lplot = 0.55*lmax print("max lambda (plot): ", lplot) Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.15, width_tangent=2, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey') # In[51]: graph = Li.plot(coordinates='xy', prange=(0, lplot), plot_points=500, thickness=3, color='green', display_tangent=True, scale=0.01, width_tangent=1, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey', axes_labels=[r'$x/m$', r'$y/m$']) graph.save("gik_spher_3d_r_39_xy.pdf") graph # ### A polar spherical photon orbit # In[52]: rph_pol(a) = 1 + 2*sqrt(1 - a^2/3)*cos(1/3*arccos((1 - a^2)/(1 - a^2/3)^(3/2))) rph_pol # In[53]: r0 = rph_pol(a0) L = lsph(a0, r0) Q = qsph(a0, r0) r0, L, Q # In[54]: theta0(a0, L, Q) # In[55]: lmax = 30 P = M.point((0, r0, pi/2, 0), name='P') Li = M.geodesic([0, lmax], P, mu=0, E=E, L=0, Q=Q, a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[56]: v0 = Li.initial_tangent_vector() Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3], a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[57]: Li.integrate(step=0.002, method='dopri5') #Li.check_integrals_of_motion(0.99*lmax) # In[58]: lplot = 0.5*lmax print("max lambda (plot): ", lplot) Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.2, width_tangent=2, color_tangent='green', plot_points_tangent=20, horizon_color='lightgrey') # ## Inner spherical orbits # To plot the inner orbits, and in particular orbits with $r<0$, we use a map from Kerr spacetime to the Euclidean space based on the radial coordinate # $$ # \hat{r} := \frac{1}{2} \left( r + \sqrt{r^2 + 4m^2} \right) # $$ # instead of the Boyer-Lindquist $r$. # In[59]: E4 = M.map_to_Euclidean().codomain() X4 = E4.cartesian_coordinates() t, x, y, z = X4[:] X4 # In[60]: rhat = 1/2*(r + sqrt(r^2 + 4)) rhat # In[61]: F = M.diff_map(E4, {(BL, X4): [t, rhat*sin(th)*cos(ph), rhat*sin(th)*sin(ph), rhat*cos(th)]}, name='F') F.display() # ### Marginally stable spherical photon orbit # In[62]: r0 = 1 - (1 - a0^2)^(1/3) r0 # In[63]: L = lsph(a0, r0) Q = qsph(a0, r0) L, Q # In[64]: theta0(a0, L, Q) # Since all orbits with $0< r_0 \leq r_{\rm ph}^*$ (such as the marginally stable spherical orbit) have $E<0$, we set `E = -1` and `L = -L`: # In[65]: lmax = 50 P = M.point((0, r0, pi/2, 0), name='P') Li = M.geodesic([0, lmax], P, mu=0, E=-1, L=-L, Q=Q, a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[66]: v0 = Li.initial_tangent_vector() Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3], a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[67]: Li.integrate(step=0.001, method='dopri5') Li.check_integrals_of_motion(0.999*lmax) # In[68]: print("Final point: ") show(BL[:], "=", BL(Li(0.999*lmax))) # Plot with respect to $\hat{r}$ coordinate: # In[69]: sing_ring = circle((0,0,0), 1, thickness=3, color='orangered') r_minf = sphere((0,0,0), 0.03, color='black') \ + text3d("r=\u2212\u221E", (0,0,0.07), fontfamily='serif', fontstyle='italic', fontsize=20) # In[70]: lplot = 0.06*lmax print("max lambda (plot): ", lplot) graph = Li.plot_integrated(chart=X4, mapping=F, ambient_coords=(x,y,z), prange=(0, lplot), thickness=3, plot_points=2000, color='green', display_tangent=True, scale=0.012, width_tangent=1, color_tangent='green', plot_points_tangent=20, label_axes=False) graph += sing_ring + r_minf \ + line([(0,0,-0.53), (0,0,0.53)], color='black', thickness=2) graph._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z'] graph # In[71]: lplot = 0.4*lmax print("max lambda (plot): ", lplot) graph = Li.plot_integrated(chart=X4, mapping=F, ambient_coords=(x,y,z), prange=(0, lplot), thickness=3, plot_points=2000, color='green', display_tangent=True, scale=0.012, width_tangent=1, color_tangent='green', plot_points_tangent=20, label_axes=False) graph += sing_ring + r_minf \ + line([(0,0,-0.53), (0,0,0.53)], color='black', thickness=2) graph._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z'] graph # ### Stable inner spherical photon orbit at $r_0 = 0.5 m$ # In[72]: r0 = 0.5 L = lsph(a0, r0) Q = qsph(a0, r0) L, Q # In[73]: theta0(a0, L, Q) # In[74]: lmax = 30 P = M.point((0, r0, pi/2, 0), name='P') Li = M.geodesic([0, lmax], P, mu=0, E=-1, L=-L, Q=Q, a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[75]: v0 = Li.initial_tangent_vector() Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3], a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[76]: Li.integrate(step=0.001, method='dopri5') Li.check_integrals_of_motion(0.999*lmax) # In[77]: print("Final point: ") show(BL[:], "=", BL(Li(0.999*lmax))) # Plot in (Cartesian) Boyer-Lindquist coordinates: # In[78]: Li.plot(prange=(0, lmax), plot_points=2000, thickness=3, color='green', display_tangent=True, scale=0.01, width_tangent=1, color_tangent='green', plot_points_tangent=20, plot_horizon=False) # Plot with respect to $\hat{r}$ coordinate: # In[79]: lplot = 0.25*lmax print("max lambda (plot): ", lplot) graph1 = Li.plot_integrated(chart=X4, mapping=F, ambient_coords=(x,y,z), prange=(0, lplot), thickness=3, plot_points=2000, color='green', display_tangent=True, scale=0.015, width_tangent=1, color_tangent='green', plot_points_tangent=20, label_axes=False) graph1 += sing_ring + r_minf graph1._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z'] graph1 # ### Vortical circular photon orbit at $r_0 = r_{\rm ph}^{**}$ # In[80]: rph_ss = n(1/2 + cos(2/3*asin(a0) + 2*pi/3)) rph_ss # In[81]: L = lsph(a0, rph_ss) Q = qsph(a0, rph_ss) L, Q # In[82]: th0 = theta0(a0, L, Q) th1 = theta1(a0, L, Q) th0, th1 # In[83]: th_ss = arcsin(sqrt(abs(L)/a0)) th_ss # In[84]: rh_ss = rhat.subs({r: rph_ss}) rh_ss # In[85]: lmax = 10 P = M.point((0, rph_ss, th_ss, 0), name='P') Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[86]: v0 = Li.initial_tangent_vector() Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=0, pph0=v0[3], a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[87]: Li.integrate(step=0.001, method='dopri5') Li.check_integrals_of_motion(0.999*lmax) # In[88]: show("Final point: ", BL[:], "=", BL(Li(0.999*lmax))) # In[89]: lplot = 0.444*lmax print("max lambda (plot): ", lplot) vort_circ = Li.plot_integrated(chart=X4, mapping=F, ambient_coords=(x,y,z), prange=(0, lplot), thickness=3, plot_points=2000, color='lightgreen', display_tangent=True, scale=0.1, width_tangent=1, color_tangent='lightgreen', plot_points_tangent=5, label_axes=False) graph = vort_circ + sing_ring + r_minf \ + line([(0,0,0), (0,0,1.3*rh_ss)], color='black', thickness=2) graph._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z'] graph # ### Vortical inner spherical photon orbit at $r_0 = -0.46 m$ # In[90]: r0 = -0.46 L = lsph(a0, r0) Q = qsph(a0, r0) L, Q # In[91]: th0 = theta0(a0, L, Q) th1 = theta1(a0, L, Q) th0, th1 # In[92]: lmax = 20 P = M.point((0, r0, (th0+th1)/2, 0), name='P') Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[93]: v0 = Li.initial_tangent_vector() Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3], a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True) # In[94]: Li.integrate(step=0.001, method='dopri5') Li.check_integrals_of_motion(0.999*lmax) # In[95]: print("Final point: ") show(BL[:], "=", BL(Li(0.999*lmax))) # In[96]: lplot = lmax print("max lambda (plot): ", lplot) graph = Li.plot_integrated(chart=X4, mapping=F, ambient_coords=(x,y,z), prange=(0, lplot), thickness=3, plot_points=2000, color='green', display_tangent=True, scale=0.1, width_tangent=1, color_tangent='green', plot_points_tangent=20, label_axes=False) # In[97]: graph += sing_ring + r_minf + vort_circ \ + line([(0,0,0), (0,0,0.9)], color='black', thickness=2) graph._extra_kwds['axes_labels'] = ['\u0302x', '\u0302y', '\u0302z'] graph # In[98]: graph + graph1