<jupyter_start><jupyter_text>Pseudo-Riemannian manifolds in SageMath The Schwarzschild spacetime example This notebook demonstrates some SageMath tools for pseudo-Riemannian geometry, developed through the [SageManifolds project](https://sagemanifolds.obspm.fr/), by these [authors](https://sagemanifolds.obspm.fr/authors.html).  This notebook requires a version of SageMath at least equal to 9.0:<jupyter_code>version()<jupyter_output><empty_output><jupyter_text>First we set up the notebook to display outputs via LaTeX rendering:<jupyter_code>%display latex<jupyter_output><empty_output><jupyter_text>Since some computations are quite heavy, we ask for running them in parallel on 8 threads:<jupyter_code>Parallelism().set(nproc=8)<jupyter_output><empty_output><jupyter_text>We introduce the **Schwarzschild spacetime**, which is the spacetime of a **static black hole** in general relativity, as a **4-dimensional Lorentzian manifold** $M$:<jupyter_code>M = Manifold(4, 'M', structure='Lorentzian')
M
print(M)<jupyter_output>4-dimensional Lorentzian manifold M<jupyter_text>$M$ is in the category of smooth manifolds over the real field:<jupyter_code>M.category()
M.base_field()<jupyter_output><empty_output><jupyter_text>At the moment, the real field is modeled by 53-bit floating-point approximations, but this plays no role in the manifold implementation:<jupyter_code>print(M.base_field())<jupyter_output>Real Field with 53 bits of precision<jupyter_text>Coordinate charts The function `Manifold` generates a manifold with no-predefined coordinate chart, so that the manifold (user) atlas is empty:<jupyter_code>M.atlas()<jupyter_output><empty_output><jupyter_text>We introduce the standard **Schwarzchild-Droste coordinates** $(t,r,\theta,\phi)$ on $M$, via the method `chart`:<jupyter_code>SD.<t, r, th, ph> = M.chart(r"t r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi:periodic")<jupyter_output><empty_output><jupyter_text>Note that the argument of `chart()` is a raw string (hence the prefix `r` in front of it), which defines the range of each coordinate, if different from $(-\infty, +\infty)$, as well as its LaTeX symbol, if different from the Python symbol to denote the coordinate. The Python variables for each coordinate are declared within the `` operator on the left-hand side, `SD` denoting the Python variable chosen for the coordinate chart.<jupyter_code>SD
print(SD)
SD.coord_range()<jupyter_output><empty_output><jupyter_text>Thanks to the SageMath operator `` used in the chart declaration, the coordinates are immediately available:<jupyter_code>th<jupyter_output><empty_output><jupyter_text>They are SageMath's symbolic expressions:<jupyter_code>th.parent()<jupyter_output><empty_output><jupyter_text>They are also accessible as items of the chart:<jupyter_code>SD[0], SD[3]
SD[:]<jupyter_output><empty_output><jupyter_text>The manifold (user) atlas is no longer empty:<jupyter_code>M.atlas()<jupyter_output><empty_output><jupyter_text>Let us introduce a second chart on the manifold, that of **Eddington-Finkelstein coordinates**  $(T,r,\theta,\phi)$:<jupyter_code>EF.<T, r, th, ph> = M.chart(r"T r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi:periodic")
EF<jupyter_output><empty_output><jupyter_text>The transition map from Schwarzschild-Droste coordinates (chart `SD`) to Eddington-Finkelstein ones (chart `EF`) depends on a parameter $m$, the mass of the Schwarzschild black hole:<jupyter_code>m = var('m')
assume(m > 0)<jupyter_output><empty_output><jupyter_text>We provide the explicit coordinate transformation via the method `transition_map`:<jupyter_code>SD_to_EF = SD.transition_map(EF, [t +2*m*ln(abs(r/(2*m)-1)), r, th, ph])
SD_to_EF.display()
SD_to_EF.inverse().display()<jupyter_output><empty_output><jupyter_text>There are now two charts in the manifold atlas:<jupyter_code>M.atlas()<jupyter_output><empty_output><jupyter_text>One of them is the so-called *default chart*: it is the chart used by any function that requires a chart as argument and none is provided by the user. At this stage, the default chart is the first chart defined on the manifold, but this can be changed by the manifold method `set_default_chart`.<jupyter_code>M.default_chart()<jupyter_output><empty_output><jupyter_text>One can plot the `SD` chart in terms of `EF` one:<jupyter_code>plot1 = SD.plot(EF, ranges={t:(0, 8), r:(2.1, 10)}, fixed_coords={th:pi/2, ph:0}, 
                ambient_coords=(r,T), style={t:'--', r:'-'}, parameters={m: 1}) \
        + SD.plot(EF, ranges={t:(0, 8), r:(0.1, 1.9)}, fixed_coords={th:pi/2, ph:0}, 
                  ambient_coords=(r,T), number_values={t: 9, r: 3},
                  style={t:'--', r:'-'}, parameters={m: 1})
plot1<jupyter_output><empty_output><jupyter_text>Manifold points To create a point on $M$, we use SageMath's *parent/element* syntax, i.e. the call operator `M(...)` acting on the parent `M`, with the point's coordinates in some chart as argument:<jupyter_code>p = M((m, 8*m, pi/2, 0), name='p')
print(p)<jupyter_output>Point p on the 4-dimensional Lorentzian manifold M<jupyter_text>Since the chart has not been specified, the default chart (i.e. `SD`) is meant:<jupyter_code>SD(p)<jupyter_output><empty_output><jupyter_text>Thanks to the transition map declared above, the coordinates of $p$ in the Eddington-Finkelstein chart can computed:<jupyter_code>EF(p)<jupyter_output><empty_output><jupyter_text>Manifold points have a `plot` method:<jupyter_code>plot1 += p.plot(EF, color='blue', ambient_coords=(r,T), 
                parameters={m: 1}, label_offset=0.4, fontsize=14)
plot1<jupyter_output><empty_output><jupyter_text>Vector fields When a chart is declared, the manifold is automatically endowed with some vector fields, those of the **coordinate vector frame**:<jupyter_code>SD.frame()
EF.frame()
M.frames()<jupyter_output><empty_output><jupyter_text>As for charts, there is a *default frame*, i.e. a vector frame that is used by default in functions having a vector frame in their arguments. The default frame can be changed by the method `set_default_frame`.<jupyter_code>M.default_frame()<jupyter_output><empty_output><jupyter_text>The first vector of the Schwarzschild-Droste coordinate frame:<jupyter_code>vt = SD.frame()[0]
vt
vt.display()
vt.display(EF.frame())<jupyter_output><empty_output><jupyter_text>The second vector of the Schwarzschild-Droste coordinate frame:<jupyter_code>vr = SD.frame()[1]
vr
vr.display(EF.frame())<jupyter_output><empty_output><jupyter_text>Creating a vector field from scratch, by providing its components with rest to a given vector frame:<jupyter_code>k = M.vector_field(1, -1, 0, 0, frame=EF.frame(), name='k')
k.display()
k.display(EF.frame())<jupyter_output><empty_output><jupyter_text>Plot with respect to Schwarzschild-Droste coordinates (default chart):<jupyter_code>k.plot(ambient_coords=(r,t), fixed_coords={th: pi/2, ph: 0}, 
       ranges={t: (0, 6), r: (0.1, 6)}, number_values=9,
       parameters={m: 1}, color='green', scale=0.3)<jupyter_output><empty_output><jupyter_text>Plot with respect to Eddington-Finkelstein coordinates:<jupyter_code>k.plot(chart=EF, ambient_coords=(r,T), chart_domain=EF,
       fixed_coords={th: pi/2, ph: 0}, ranges={T: (0, 8), r: (0.1, 10)}, 
       number_values=5, parameters={m: 1}, color='green')<jupyter_output><empty_output><jupyter_text>Vector fields as sections of the tangent bundle<jupyter_code>TM = M.tangent_bundle()
print(TM)
k1 = TM.section({EF.frame(): [1, -1, 0, 0]})
print(k1)
k1 == k<jupyter_output><empty_output><jupyter_text>The set of all vector fields on $M$ as a $C^\infty(M)$-module:<jupyter_code>XM = M.vector_field_module()
print(XM)
XM
XM.base_ring()
print(XM.base_ring())
XM.base_ring().an_element().display()<jupyter_output><empty_output><jupyter_text>Vectors at a point The value of a vector field, and more generally of any tensor field, at a given point is obtained by the method `at`:<jupyter_code>kp = k.at(p)
print(kp)
kp.display()<jupyter_output><empty_output><jupyter_text>The parent of `kp` is the tangent space at `p`:<jupyter_code>kp.parent()
print(kp.parent())<jupyter_output>Tangent space at Point p on the 4-dimensional Lorentzian manifold M<jupyter_text>It is accessible from the manifold via the method `tangent_space`:<jupyter_code>kp.parent() is M.tangent_space(p)<jupyter_output><empty_output><jupyter_text>$T_p M$ is the fiber over $p$ in the tangent bundle $TM$:<jupyter_code>kp.parent() is TM.fiber(p)<jupyter_output><empty_output><jupyter_text>Tangent vectors have a method `plot`:<jupyter_code>plot1 += kp.plot(EF, color='green', ambient_coords=(r,T), 
                 parameters={m: 1}, scale=2, label_offset=0.5, 
                 fontsize=16)
plot1<jupyter_output><empty_output><jupyter_text>Vector field defined on an open subset Let us introduce the exterior $E$ of the black hole as an open subset of $M$:<jupyter_code>E = M.open_subset('E', coord_def = {SD: r>2*m})
SD.restrict(E).coord_range()
p in E
u = E.vector_field(name='u')
u[0] = 1/sqrt(1-2*m/r)
u.display()
SD_to_EF.restrict(E)
E.atlas()
u.display(EF.frame().restrict(E))
u.plot(ambient_coords=(r, t), fixed_coords={th: pi/2, ph: 0}, 
       ranges={t: (0, 6), r: (2.1, 6)}, number_values=5, 
       parameters={m: 1}, scale=0.3)
u.at(p).display()<jupyter_output><empty_output><jupyter_text>Metric tensorWe define next the **metric tensor** $g$ from its non-vanishing components in the manifold's default frame, namely the coordinate frame associated to Schwarzschild-Droste coordinate:<jupyter_code>g = M.metric()
g[0, 0] = - (1 - 2*m/r)
g[1, 1] = 1/(1 - 2*m/r)
g[2, 2] = r^2
g[3, 3] = r^2*sin(th)^2
g.display()
SD.coframe()
g[:]
g[1,1]<jupyter_output><empty_output><jupyter_text>$g_{rr}$ is diverging at $r=2m$: this is a singularity of the Schwarszchild-Droste coordinates.<jupyter_code>g.display_comp()<jupyter_output><empty_output><jupyter_text>The components of $g$ with respect to the Eddington-Finkelstein frame are evaluated via the methods `comp` or `display`:<jupyter_code>g.display(EF.frame())
g[EF.frame(),:]<jupyter_output><empty_output><jupyter_text>Note that these components are regular at $r=2m$, contrary to the components in Schwarzschild-Droste coordinates. The metric tensor is a twice-covariant tensor (actually a field of symmetric bilinear forms):<jupyter_code>g.tensor_type()<jupyter_output><empty_output><jupyter_text>It can be thus be applied to a pair of vector fields, for instance $u$ and $k$:<jupyter_code>s = g(u, k)
s
print(s)
s.display()
s.expr()
g(k, k).display()<jupyter_output><empty_output><jupyter_text>$u$ is a unit timelike vector:<jupyter_code>g(u, u).display()<jupyter_output><empty_output><jupyter_text>Scalar products returned by the method `dot` are actually those formed with $g$:<jupyter_code>u.dot(u)
u.dot(u).display()<jupyter_output><empty_output><jupyter_text>Levi-Civita connection<jupyter_code>nabla = g.connection()
print(nabla)
nabla.display()
Parallelism().set(nproc=1)
Dk = nabla(k)
Dk
print(Dk)
Dk.display()
Dk.display_comp()<jupyter_output><empty_output><jupyter_text>The acceleration $\nabla_k k$ of the vector field $k$:<jupyter_code>Ak = Dk.contract(k)
print(Ak)
Ak.display()<jupyter_output><empty_output><jupyter_text>$k$ is thus geodesic vector. The field lines of $k$ are actually the ingoing radial null geodesics of Schwarzschild spacetime. Instead of the method contract, one can use index notations to compute $\nabla_k k$ as $k^b \nabla_b k^a$:<jupyter_code>Ak == Dk['^a_b']*k['^b']
Parallelism().set(nproc=8)<jupyter_output><empty_output><jupyter_text>The acceleration $\nabla_u u$ of the vector field $u$:<jupyter_code>Au = nabla(u).contract(u)
Au.display()<jupyter_output><empty_output><jupyter_text>$\nabla g$ is identically zero, since the connection $\nabla$ is compatible with $g$:<jupyter_code>Dg = nabla(g)
print(Dg)
Dg.display()<jupyter_output><empty_output><jupyter_text>Curvature The **Riemann curvature tensor** is computed as<jupyter_code>Riem = g.riemann()
print(Riem)
Riem.display_comp(only_nonredundant=True)<jupyter_output><empty_output><jupyter_text>The component $\mathrm{Riem}(g)^t_{\ \, rtr} = \mathrm{Riem}(g)^0_{\ \, 101}$ is returned by<jupyter_code>Riem[0,1,0,1]
Riem.display_comp(EF.frame(), only_nonredundant=True)<jupyter_output><empty_output><jupyter_text>The **Kretschmann scalar** is the "square" of the Riemann tensor defined by $$K = R_{abcd} \, R^{abcd}, \qquad R := \mathrm{Riem}(g)$$To compute it, we must first form the tensor fields whose components are $R_{abcd}$ and $R^{abcd}$. They are obtained by respectively lowering and raising the indices of the components $R^a_{\ \, bcd}$ of the Riemann tensor, via the metric $g$. These two operations are performed by the methods `down()` and `up()`. The contraction is performed by summation on repeated indices, using LaTeX notations:<jupyter_code>K = Riem.down(g)['_{abcd}'] * Riem.up(g)['^{abcd}']
print(K)
K.display()
K.expr()<jupyter_output><empty_output><jupyter_text>Since $\lim_{r\to 0} K = +\infty$, we may say that $r=0$ is a **curvature singularity** of Schwarzschild spacetime.  Ricci tensor<jupyter_code>Ric = g.ricci()
print(Ric)<jupyter_output>Field of symmetric bilinear forms Ric(g) on the 4-dimensional Lorentzian manifold M<jupyter_text>We check that the Schwarzschild metric is a solution of the vacuum **Einstein equation**:<jupyter_code>Ric.display()<jupyter_output><empty_output><jupyter_text>Geodesics First, for graphical purposes, we introduce the **Euclidean space** $\mathbb{E}^3$ and some map $M\to \mathbb{E}^3$:<jupyter_code>E3.<x,y,z> = EuclideanSpace()
X3 = E3.cartesian_coordinates()
to_E3 = M.diff_map(E3, {(SD, X3): 
                        [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)]})
to_E3.display()<jupyter_output><empty_output><jupyter_text>A timelike geodesic Let us consider the geodesic starting at point $p$ and having the following tangent vector at $p$ (note that the tangent vector $v_0$ is constructed by means of the call operator `()` acting of the parent, which is the tangent space to $M$ at $p$):<jupyter_code>v0 = M.tangent_space(p)((1.3, 0, 0, 0.064/m), name='v_0')
v0.display()<jupyter_output><empty_output><jupyter_text>We declare the geodesic with such initial conditions, denoting by $s$ the affine parameter (proper time), with $(s_{\rm min}, s_{\rm max})=(0, 1500\,m)$:<jupyter_code>s = var('s')
geod = M.integrated_geodesic(g, (s, 0, 2000), v0)
geod<jupyter_output><empty_output><jupyter_text>Note that the initial point `p` is not explicitely passed in the argument list of `integrated_geodesic`, because this piece of information is contained in `v0`:<jupyter_code>p is v0.parent().base_point()
sol = geod.solve(parameters_values={m: 1})  # numerical integration
interp = geod.interpolate()                 # interpolation of the solution for the plot
geod.plot_integrated(chart=X3, mapping=to_E3, plot_points=1000, 
                     thickness=2, label_axes=False) \
+ p.plot(chart=X3, mapping=to_E3, size=4, parameters={m: 1}) \
+ sphere(size=2, color='grey')<jupyter_output><empty_output><jupyter_text>A 2D view by suppressing $z$ from the ambient coordinates:<jupyter_code>bh_plot = circle((0, 0), 2, edgecolor='black', fill=True, facecolor='grey', alpha=0.5)
geod.plot_integrated(chart=X3, mapping=to_E3, ambient_coords=(x,y), plot_points=1000, 
                     thickness=2) \
+ p.plot(chart=X3, mapping=to_E3, ambient_coords=(x,y), size=4, parameters={m: 1}) \
+ bh_plot<jupyter_output><empty_output><jupyter_text>Null geodesics Let us consider a null geodesique $\mathscr{L}$ in the equatorial plane ($\theta = \pi/2$). The null vector $v$ tangent to $\mathscr{L}$ and associated to some affine parameter $\lambda$ is given by the following first integrals of the geodesic equation:$$    v^t = \frac{\mathrm{d}t}{\mathrm{d}\lambda} =  \left(1 - \frac{2m}{r} \right)^{-1}$$$$    v^r = \frac{\mathrm{d}r}{\mathrm{d}\lambda} = \pm \sqrt{1 - \frac{b^2}{r^2} \left(1 - \frac{2m}{r} \right) }$$$$    v^\theta = \frac{\mathrm{d}\theta}{\mathrm{d}\lambda} = 0$$$$    v^\varphi = \frac{\mathrm{d}\varphi}{\mathrm{d}\lambda} = \frac{b}{r^2}$$where the constant $b$ is related to the conserved energy $E$ and conserved angular momentum $L$ along the geodesic by $$    b = \frac{L}{E} . $$For a geodesic arising from infinity, $b$ is the **impact parameter**. To set up the initial vector for the computation of a null geodesic, let us define a function that takes $b$ and some initial radius $r_0$ as input and returns the initial vector $v_0$. We take advantage that SageMath is built atop Python to construct this function as a pure Python function:<jupyter_code>def initial_vector(r0, b, phi0=0, inward=True):
    r"""
    Evaluate the initial tangent vector along a null geodesic. 
    
    INPUT:
    
    - r0: radial SD coordinate of the initial point
    - b: impact parameter
    - phi0: azimuthal SD coordinate of the initial point (default: 0)
    - inward: determines whether the geodesic has initially v^r < 0 (default: True)
    
    """
    vt0 = 1/(1 - 2*m/r0)
    vr0 = sqrt(1 - b^2/r0^2*(1 - 2*m/r0))
    if inward:
        vr0 = - vr0
    vth0 = 0
    vph0 = b / r0^2
    p0 = M((0, r0, pi/2, phi0), name='p_0')  # initial point
    return M.tangent_space(p0)((vt0, vr0, vth0, vph0), name='v_0')<jupyter_output><empty_output><jupyter_text>Let us use this function to construct the initial vector for $r_0 = 10m$ and $b=7m$:<jupyter_code>v0 = initial_vector(10*m, 7*m)
v0.display()<jupyter_output><empty_output><jupyter_text>Let us check that $v_0$ is a null vector:<jupyter_code>p0 = v0.parent().base_point()
g.at(p0)(v0, v0)
geod = M.integrated_geodesic(g, (s, 0, 40), v0)
sol = geod.solve(step=0.01, parameters_values={m: 1}) 
interp = geod.interpolate()   
plot2 = geod.plot_integrated(chart=X3, mapping=to_E3, ambient_coords=(x,y), 
                             plot_points=500, color='green', thickness=1.5, display_tangent=True, 
                             color_tangent='green', plot_points_tangent=4, scale=1) 
plot2 += bh_plot
plot2<jupyter_output><empty_output><jupyter_text>A null geodesic plunging into the black hole:<jupyter_code>v0 = initial_vector(10*m, 5*m)
v0.display()
geod = M.integrated_geodesic(g, (s, 0, 13), v0)
sol = geod.solve(step=0.01, parameters_values={m: 1}) 
interp = geod.interpolate()   
plot2 += geod.plot_integrated(chart=X3, mapping=to_E3, ambient_coords=(x,y), 
                              plot_points=500, color='green', thickness=1.5, display_tangent=True, 
                              color_tangent='green', plot_points_tangent=3, scale=0.2)
plot2<jupyter_output><empty_output><jupyter_text>The photon orbit The photon orbit corresponds to $r_0=3m$ with the following critical value of the impact parameter: $$    b_{\rm c} = 3\sqrt{3}\, m$$<jupyter_code>bc = 3*sqrt(3)*m
n(bc/m)
v0 = initial_vector(3*m, bc)
v0.display()
geod = M.integrated_geodesic(g, (s, 0, 100), v0)
sol = geod.solve(step=0.01, parameters_values={m: 1}) 
interp = geod.interpolate()   
plot2 += geod.plot_integrated(chart=X3, mapping=to_E3, ambient_coords=(x,y), 
                              plot_points=500, color='green', thickness=1.5)
plot2<jupyter_output><empty_output><jupyter_text>A geodesic with $b$ close to $b_{\rm c}$ wraps around the circular orbit:<jupyter_code>v0 = initial_vector(10*m, 5.2025*m)
geod = M.integrated_geodesic(g, (s, 0, 40), v0)
sol = geod.solve(step=0.01, parameters_values={m: 1}) 
interp = geod.interpolate()   
plot2 += geod.plot_integrated(chart=X3, mapping=to_E3, ambient_coords=(x,y), 
                              plot_points=500, color='orange', thickness=1.5,
                              display_tangent=True, color_tangent='orange', 
                              plot_points_tangent=4, scale=1)
plot2<jupyter_output><empty_output><jupyter_text>Using SymPy as the symbolic backendBy default, the symbolic backend used in tensor calculus is SageMath's one (Pynac + Maxima), implemented via the symbolic ring `SR`. We can choose to use SymPy instead:<jupyter_code>M.set_calculus_method('sympy')
v = 2*k
v.display()
v[0]
v[0].expr()
type(v[0].expr())
M.set_calculus_method('SR')
type(v[0].expr())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13954990312484475, 'fraction_numerical': 0.01912663321575836, 'mean_word_length': 5.78006062647356, 'pattern_counts': {'":': 0, '<': 268, '<?xml version=': 0, '>': 269, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '67', 'dir': 'jupyter-structured-clean-dedup', 'id': '370189', 'license': 'MIT', 'n_tokens_mistral': 7547, 'n_tokens_neox': 7123, 'n_words': 1844, 'path': 'otros/demo_pseudo_Riemannian_Schwarzschild.ipynb', 'repo_name': 'ashcat2005/SchwarzschildBH', 'size': '20129'}