This SageMath notebook implements the computation of the Ricci rotation coefficients of a given tetrad, as well as the Newman-Penrose and GHP spin coefficients.
An example is provided by the Kinnersley null tetrad in Kerr spacetime.
Author: Eric Gourgoulhon (2024)
version()
'SageMath version 10.4.rc1, Release Date: 2024-06-27'
%display latex
Since some computations are quite heavy, we ask for running them in parallel on 8 threads:
Parallelism().set(nproc=8)
M = Manifold(4, 'M', structure='Lorentzian')
We sall use the advanced Kerr coordinates $(v,r,\theta,\varphi)$:
X.<v,r,th,ph> = M.chart(r'v r th:\theta:(0,pi) ph:\varphi:(0,2*pi)')
X
Expression of the Kerr metric in terms of the advanced Kerr coordinates:
m, a = var('m a', domain='real')
Delta = r^2 - 2*m*r + a^2
Sigma = r^2 + a^2*cos(th)^2
g = M.metric()
g[0,0] = - 1 + 2*m*r/Sigma
g[0,1] = 1
g[0,3] = - 2*a*m*r*sin(th)^2/Sigma
g[1,3] = - a*sin(th)^2
g[2,2] = Sigma
g[3,3] = (r^2 + a^2 + 2*a^2*m*r*sin(th)^2/Sigma)*sin(th)^2
g.display()
g.inverse().display()
Let us first introduce some auxiliary quantities (remark: we use I
for the imaginary number $i$, keeping i
for indices):
rho = 1 / (I*a*cos(th) - r)
rho
rho_bar = rho.conjugate().simplify_full()
rho_bar
(rho*rho_bar).simplify_full()
We then introduce the vector fields of the Kinnersley tetrad $(\ell_{\rm K}, n_{\rm K}, m_{\rm K}, \bar{m}_{\rm K})$ one by one:
lK = M.vector_field(2*(r^2 + a^2)/Delta, 1, 0, 2*a/Delta,
name='lK', latex_name=r'\ell_{\rm K}')
lK.display()
nK = M.vector_field(0, -Delta/(2*Sigma), 0, 0,
name='nK', latex_name=r'n_{\rm K}')
nK.display()
mK = - rho_bar/sqrt(2) * M.vector_field(I*a*sin(th), 0, 1, I/sin(th))
mK.set_name('mK', latex_name=r'm_{\rm K}')
mK.display()
mbK = M.vector_field(name='mbK', latex_name=r'\bar{m}_{\rm K}')
for i in M.irange():
mbK[i] = mK[i].expr().conjugate()
mbK.display()
Let us check that each vector of the tetrad is a null vector:
g(lK, lK).expr()
g(nK, nK).expr()
g(mK, mK).expr()
g(mbK, mbK).expr()
Furthermore, we have the following scalar products:
g(lK, nK).expr(), g(mK,mbK).expr()
g(lK, mK).expr(), g(nK, mK).expr()
g(lK, mbK).expr(), g(nK, mbK).expr()
Let us define the Kinnersley tetrad K
as the vector frame $(\ell_{\rm K}, n_{\rm K}, m_{\rm K}, \bar{m}_{\rm K})$ on $M$, using the symbol $e_{\rm K}$ for a generic member of the tetrad:
K = M.vector_frame('eK', (lK, nK, mK, mbK), latex_symbol=r'{e_{\rm K}}')
K
We have then:
for vec in (lK, nK, mK, mbK):
show(vec.display(K))
The dual basis of the Kinnersley tetrad (each member of the dual is a 1-form and is labelled with an upper index):
for form in K.coframe():
show(form.display())
We notice a lack of simplification. To fix it, we define a simplifying function that invokes the algebraic field $\bar{\mathbb{Q}}$:
def algebraic_simplify(expr):
r"""
Simplifies a symbolic expression through the fraction field of polynomials
over QQbar
INPUT:
- ``expr`` -- a Sage's symbolic expression (element of the symbolic ring SR)
OUTPUT:
- a Sage's symbolic expression
"""
a00, m00, r00, y00, z00, s00 = SR.var('a00 m00 r00 y00 z00 s00')
# From symbolic variables to variables in QQbar:
expr = expr.subs({a: a00, m: m00, r: r00, cos(th): y00, sin(th): z00, sqrt(2): s00})
expr = expr.subs({z00^2: 1 - y00^2}) # sin(th)^2 --> 1 - cos(th)^2
F = QQbar['a00', 'm00', 'r00', 'y00', 'z00', 's00'].fraction_field()
expr = SR(F(str(expr))).simplify()
# Back to the original symbolic variables:
expr = expr.subs({a00: a, m00: m, r00: r, y00: cos(th), z00: sin(th), s00: sqrt(2)})
# Some extra simplifications:
expr = expr.subs({1 - cos(th)^2: sin(th)^2})
expr = expr.subs({cos(th)^2 - 1: -sin(th)^2})
w0 = SR.wild(0)
expr = expr.subs({w0 - w0*cos(th)^2: w0*sin(th)^2})
expr = expr.subs({w0*cos(th)^2 - w0: -w0*sin(th)^2})
expr = expr.factor().simplify_trig().factor()
return expr
Let us check it:
s1, s2 = K.coframe()[3][3].expr(), K.coframe()[1][3].expr()
s1, s2
algebraic_simplify(s1), algebraic_simplify(s2)
Let us apply the simplifying function to all the coframe elements:
for f in K.coframe():
f.apply_map(algebraic_simplify, keep_other_components=True)
show(f.display())
Let us check that K
is a null NP tetrad by asking for the components of the metric tensor w.r.t. it:
g.display(K)
A matrix view of the components:
g[K,:]
Let us evaluate the metric duals (not the frame dual!) of each of the vector of the Kinnersley frame. These 1-forms are obtained via down(g)
(index lowering with the metric $g$):
lKf = lK.down(g)
lKf.set_name('lKf', latex_name=r'\underline{\ell}_{\rm K}')
nKf = nK.down(g)
nKf.set_name('nKf', latex_name=r'\underline{n}_{\rm K}')
mKf = mK.down(g)
mKf.set_name('mKf', latex_name=r'\underline{m}_{\rm K}')
mbKf = mbK.down(g)
mbKf.set_name('mbKf', latex_name=r'\underline{\bar{m}}_{\rm K}')
for f in (lKf, nKf, mKf, mbKf):
f.apply_map(algebraic_simplify, keep_other_components=True)
show(f.display())
for f in (lKf, nKf, mKf, mbKf):
show(f.display(K))
The Levi-Civita connection $\nabla$ is returned by the method connection
:
nabla = g.connection()
The connection coefficients with respect to the Kinnersley tetrad would be returned by simply typing nabla.coef(K)
. However the Ricci rotation coefficients do not coincide exactly with those.
Indeed, the connection coefficients $\Gamma^k_{\ \, ij}$ are defined by
$$ \Gamma^i_{\ \, jk} := \langle e^i ,\nabla_{e_k} e_j\rangle $$
while the Ricci rotation coefficients are defined by
$$ \Gamma_{ijk} := e_i\cdot \nabla_{e_k} e_j \qquad\qquad \mbox{(1)}$$
where a dot stands for the scalar product with respect to the metric $g$.
If $g_{ij} := g(e_i,e_j) = e_i\cdot e_j$, then
$$ \Gamma_{ijk} = g_{il} \Gamma^{l}_{\ \, jk} . $$
We implement here a Python function that computes the Ricci rotation coefficients from (1) and returns them as an object of type CompWithSym
, in order to take into account the antisymmetry with respect to the first two indices ($\Gamma_{ijk} = - \Gamma_{jik}$):
def ricci_coef(tetrad, metric=g, chart=M.default_chart(),
simplify=algebraic_simplify, verbose=False):
r"""
Compute the Ricci rotation coefficients with respect to a given tetrad.
The Ricci rotation coefficients `Gamma_{ijk}` are stored in a CompWithSym
object ``C`` such that ``C[i,j,k]`` ` = \Gamma_{ijk}`
INPUT:
- ``tetrad`` -- the tetrad, as an instance of VectorFrame
- ``metric`` -- the metric defining the Levi-Civita connection
- ``chart`` -- the coordinate chart in which the computation is performed
- ``simplify`` -- a simplifying function
- ``verbose`` -- if ``True``, loop indices are printed out during
the computation
OUTPUT:
- instance of ``CompWithSym`` providing the Ricci rotation coefficients
"""
from sage.tensor.modules.comp import CompWithSym
from sage.manifolds.differentiable.scalarfield import DiffScalarField
ee = tetrad # shortcut notations
gg = metric #
nabla = gg.connection()
resu = CompWithSym(M.scalar_field_algebra(), # each $\Gamma_{ijk}$ is a scalar field
tetrad, # these are coefficients with respect to the given tetrad
3, # 3-index objects
start_index=M.start_index(), # same starting index as for the manifold
output_formatter=DiffScalarField.coord_function, # for display
antisym=(0,1)) # antisymmetry with respect to the first two indices
for j in M.irange(end=M.dim() - 2 - M.start_index()):
if verbose: print("j = ", j)
nab_ej = nabla(ee[j]) # \nabla e_j
for k in M.irange():
nab_ek_ej = nab_ej.contract(ee[k]) # \nabla_{e_k} e_j
nab_ek_ej.apply_map(simplify, frame=chart.frame(), chart=chart)
for i in M.irange(start=j+1):
if verbose: print('i j k : ', i, j, k)
# formula (1) above :
resu[[i,j,k]] = simplify( gg(ee[i], nab_ek_ej).expr(chart) )
return resu
%time omK = ricci_coef(K, verbose=True)
j = 0 i j k : 1 0 0 i j k : 2 0 0 i j k : 3 0 0 i j k : 1 0 1 i j k : 2 0 1 i j k : 3 0 1 i j k : 1 0 2 i j k : 2 0 2 i j k : 3 0 2 i j k : 1 0 3 i j k : 2 0 3 i j k : 3 0 3 j = 1 i j k : 2 1 0 i j k : 3 1 0 i j k : 2 1 1 i j k : 3 1 1 i j k : 2 1 2 i j k : 3 1 2 i j k : 2 1 3 i j k : 3 1 3 j = 2 i j k : 3 2 0 i j k : 3 2 1 i j k : 3 2 2 i j k : 3 2 3 CPU times: user 2min 27s, sys: 2.44 s, total: 2min 29s Wall time: 2min 41s
def spin_coef_NP(rcf, simplify=algebraic_simplify, signature='positive'):
r"""
Evaluate the Newman-Penrose spin coefficients associated to a tetrad,
from the Ricci rotation coefficients of this tetrad.
INPUT:
- ``rcf`` -- the Ricci rotation coefficients, as an instance of CompWithSym
- ``simplify`` -- a simplifying function
- ``signature`` -- the metric signature (either ``'positive'``
or ``'negative'``)
OUPUT:
- a dictionary containing the 12 Newman-Penrose spin coefficients, with
the Greek symbols of the coefficients as the dictionary keys.
"""
resu = {}
resu['κ'] = rcf[0,2,0]
resu['ρ'] = rcf[0,2,3]
resu['σ'] = rcf[0,2,2]
resu['τ'] = rcf[0,2,1]
resu['ν'] = rcf[3,1,1]
resu['μ'] = rcf[3,1,2]
resu['λ'] = rcf[3,1,3]
resu['π'] = rcf[3,1,0]
resu['ϵ'] = (rcf[0,1,0] - rcf[2,3,0])/2
resu['α'] = (rcf[0,1,3] - rcf[2,3,3])/2
resu['β'] = (rcf[0,1,2] - rcf[2,3,2])/2
resu['γ'] = (rcf[0,1,1] - rcf[2,3,1])/2
if signature == 'negative':
for name, sc in resu.items():
resu[name] = - sc
elif signature != 'positive':
raise ValueError("signature must be either 'positive' or 'negative'")
# Conversion to symbolic expressions and simplification:
for name, sc in resu.items():
resu[name] = simplify(sc.expr())
return resu
scNP = spin_coef_NP(omK)
scNP
def spin_coef_GHP(rcf, simplify=algebraic_simplify, signature='positive'):
r"""
Evaluate the Geroch-Held-Penrose spin coefficients associated to a tetrad,
from the Ricci rotation coefficients of this tetrad.
INPUT:
- ``rcf`` -- the Ricci rotation coefficients, as an instance of CompWithSym
- ``simplify`` -- a simplifying function
- ``signature`` -- the metric signature (either ``'positive'``
or ``'negative'``)
OUPUT:
- a dictionary containing the 12 Geroch-Help-Penrose spin coefficients, with
the Greek symbols of the coefficients as the dictionary keys.
"""
scNP = spin_coef_NP(rcf, simplify=simplify, signature=signature)
resu = {}
for symb in ('β', 'ϵ', 'ρ', 'τ', 'σ', 'κ'):
resu[symb] = scNP[symb]
resu["β'"] = -scNP["α"]
resu["ϵ'"] = -scNP["γ"]
resu["ρ'"] = -scNP["μ"]
resu["τ'"] = -scNP["π"]
resu["σ'"] = -scNP["λ"]
resu["κ'"] = -scNP["ν"]
return resu
sc = spin_coef_GHP(omK)
sc
sc["ρ"]
bool(sc["ρ"] == rho)
sc["ρ'"]
bool(sc["ρ'"] == -rho*Delta/(2*Sigma))
sc["τ"]
bool(sc["τ"] == -I*a*sin(th)/(sqrt(2)*Sigma))
sc["τ'"]
bool(sc["τ'"] == -I*a*rho^2*sin(th)/sqrt(2))
sc["β"]
bool(sc["β"] == -rho_bar/(2*sqrt(2)*tan(th)))
sc["β'"]
bool(sc["β'"] == sc["β"].conjugate() + sc["τ'"])
Besides, we have $\epsilon' \neq 0$:
sc["ϵ'"]
eps_prime0 = (m*(a^2*cos(th)^2 - r^2) + a^2*r*sin(th)^2
+ I*a*Delta*cos(th)) / (2*(r^2 + a^2*cos(th)^2)^2)
eps_prime0
bool(sc["ϵ'"] == eps_prime0)
All the other spin coefficients are zero:
sc["ϵ"], sc["σ"], sc["σ'"], sc["κ"], sc["κ'"]