# coding: utf-8 # # Timelike orbits in Schwarzschild spacetime # # This Jupyter/SageMath worksheet is relative to the lectures # [Geometry and physics of black holes](http://luth.obspm.fr/~luthier/gourgoulhon/bh16/). # # Click [here](https://raw.githubusercontent.com/egourgoulhon/BHLectures/master/sage/ges_orbits.ipynb) to download the worksheet file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, with the command sage -n jupyter # # It uses the integrated_geodesic functionality introduced by Karim Van Aelst in **SageMath 8.1**, in the framework of the [SageManifolds](http://sagemanifolds.obspm.fr) project. # A version of SageMath at least equal to 8.1 is required to run this worksheet: # In[1]: version() # In[2]: get_ipython().run_line_magic('display', 'latex # LaTeX rendering turned on') # We define first the spacetime manifold $M$ and the standard Schwarzschild-Droste coordinates on it: # In[3]: M = Manifold(4, 'M') X. = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:\phi') X # For graphical purposes, we introduce $\mathbb{R}^3$ and some coordinate map $M\to \mathbb{R}^3$: # In[4]: R3 = Manifold(3, 'R^3', latex_name=r'\mathbb{R}^3') X3. = R3.chart() to_R3 = M.diff_map(R3, {(X, X3): [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)]}) to_R3.display() # Next, we define the Schwarzschild metric: # In[5]: g = M.lorentzian_metric('g') m = var('m'); assume(m >= 0) g[0,0], g[1,1] = -(1-2*m/r), 1/(1-2*m/r) g[2,2], g[3,3] = r^2, (r*sin(th))^2 g.display() # We set the specific conserved energy and angular momentum: # In[6]: #var('E L', domain='real') E = 0.973 L = 4.2*m # In[7]: u = M.vector_field(name='u') u[0] = E/(1-2*m/r) u[1] = sqrt(E^2 - (1-2*m/r)*(1+L^2/r^2)) u[3] = L/(r^2*sin(th)^2) u.display() # In[8]: g(u,u).display() # Pericenter and apocenter: # In[9]: eq = (r^3*((1-2/r)*(1+L^2/m^2/r^2) - E^2)).simplify_full() eq # In[10]: R. = RR[] print(R) # In[11]: eqp = R(eq.subs(r=w)) eqp # In[12]: roots = sorted(eqp.real_roots()) roots # In[13]: r_per = roots[1]*m r_apo = roots[2]*m r_per, r_apo # We pick an initial point and an initial tangent vector: # In[14]: r0 = r_per p0 = M.point((0, r0, pi/2, 1e-12), name='p_0') Tp0 = M.tangent_space(p0) u0 = u.at(p0) u0.set_name('u_0') print(u0) # In[15]: u0.display() # Let us check that the scalar square of $u_0$ is $-1$, by means of the metric tensor at $p_0$: # In[16]: g0 = g.at(p0) print(g0) # The scalar square is indeed equal to $-1$: # In[17]: g0(u0,u0).simplify_full() # ### Check # # Let us compute the conserved energy and angular momentum along the geodesic by # taking scalar products of the 4-velocity $v_0$ with the Killing vectors $\xi=\frac{\partial}{\partial t}$ and $\eta=\frac{\partial}{\partial t}$: # In[18]: xi = X.frame()[0] eta = X.frame()[3] xi, eta # In[19]: xi0 = xi.at(p0) eta0 = eta.at(p0) print(xi0) print(eta0) # The specific conserved energy $\varepsilon$ is: # In[20]: - g0(xi0, u0) # while the specific conserved angular momentum $\ell$ is # In[21]: g0(eta0, u0) # We declare the geodesic through $p_0$ having $v_0$ as inital vector, denoting by $s$ the affine parameter (proper time): # In[22]: s = var('s') geod = M.integrated_geodesic(g, (s, 0, 1500), u0); geod # We ask for the numerical integration of the geodesic, providing some numerical value for the parameter $m$, and then plot it in terms of the Cartesian chart: # In[23]: sol = geod.solve(step=4, parameters_values={m: 1}) # numerical integration interp = geod.interpolate() # interpolation of the solution for the plot # In[24]: graph = geod.plot_integrated(chart=X3, mapping=to_R3, ambient_coords=(x,y), plot_points=500, thickness=2, label_axes=False) # In[25]: graph += circle((0,0), r_per/m, linestyle=':') graph += circle((0,0), r_apo/m, linestyle=':') graph += circle((0,0), 2, fill=True, edgecolor='grey', facecolor='grey') show(graph, aspect_ratio=1, axes_labels=[r'$x/m$', r'$y/m$']) # In[26]: file_name = "ges_orbit_e{:.3f}_l{:.3f}.pdf".format(float(E), float(L/m)) graph.save(file_name, aspect_ratio=1, axes_labels=[r'$x/m$', r'$y/m$']) file_name # Some details about the system solved to get the geodesic: # In[27]: geod.system(verbose=True) # We recognize in the above list the Christoffel symbols of the metric $g$: # In[28]: g.christoffel_symbols_display()