The 2 parameters $m$ and $a$ of the Kerr spacetime are declared as symbolic variables:
# In[8]: var('m, a', domain='real') # We get the (yet undefined) spacetime metric: # In[9]: g = M.metric() #The metric is set by its components in the coordinate frame associated with Boyer-Lindquist coordinates, which is the current manifold's default frame:
# In[10]: rho2 = r^2 + (a*cos(th))^2 Delta = r^2 -2*m*r + a^2 g[0,0] = -(1-2*m*r/rho2) g[0,3] = -2*a*m*r*sin(th)^2/rho2 g[1,1], g[2,2] = rho2/Delta, rho2 g[3,3] = (r^2+a^2+2*m*r*(a*sin(th))^2/rho2)*sin(th)^2 g.display() #A matrix view of the components with respect to the manifold's default vector frame:
# In[11]: g[:] #The list of the non-vanishing components:
# In[12]: g.display_comp() #The Levi-Civita connection $\nabla$ associated with $g$:
# In[13]: nabla = g.connection() print(nabla) #Let us verify that the covariant derivative of $g$ with respect to $\nabla$ vanishes identically:
# In[14]: nabla(g) == 0 #Another view of the above property:
# In[15]: nabla(g).display() #The nonzero Christoffel symbols (skipping those that can be deduced by symmetry of the last two indices):
# In[16]: g.christoffel_symbols_display() #The default vector frame on the spacetime manifold is the coordinate basis associated with Boyer-Lindquist coordinates:
# In[17]: M.default_frame() is BL.frame() # In[18]: BL.frame() #Let us consider the first vector field of this frame:
# In[19]: xi = BL.frame()[0]; xi # In[20]: print(xi) #The 1-form associated to it by metric duality is
# In[21]: xi_form = xi.down(g) xi_form.display() #Its covariant derivative is
# In[22]: nab_xi = nabla(xi_form) print(nab_xi) nab_xi.display() #Let us check that the Killing equation is satisfied:
# In[23]: nab_xi.symmetrize() == 0 #Similarly, let us check that $\frac{\partial}{\partial\phi}$ is a Killing vector:
# In[24]: chi = BL.frame()[3] ; chi # In[25]: nabla(chi.down(g)).symmetrize() == 0 #The Ricci tensor associated with $g$:
# In[26]: Ric = g.ricci() print(Ric) #Let us check that the Kerr metric is a solution of the vacuum Einstein equation:
# In[27]: Ric == 0 #Another view of the above property:
# In[28]: Ric.display() #The Riemann curvature tensor associated with $g$:
# In[29]: R = g.riemann() print(R) #Contrary to the Ricci tensor, the Riemann tensor does not vanish; for instance, the component $R^0_{\ \, 123}$ is
# In[30]: R[0,1,2,3] #Let us check the Bianchi identity $\nabla_p R^i_{\ \, j kl} + \nabla_k R^i_{\ \, jlp} + \nabla_l R^i_{\ \, jpk} = 0$:
# In[31]: DR = nabla(R) #long (takes a while) print(DR) # In[32]: from __future__ import print_function # because the print below is valid only in Python3 for i in M.irange(): for j in M.irange(): for k in M.irange(): for l in M.irange(): for p in M.irange(): print(DR[i,j,k,l,p] + DR[i,j,l,p,k] + DR[i,j,p,k,l], end=' ') #If the last sign in the Bianchi identity is changed to minus, the identity does no longer hold:
# In[33]: DR[0,1,2,3,1] + DR[0,1,3,1,2] + DR[0,1,1,2,3] # should be zero (Bianchi identity) # In[34]: DR[0,1,2,3,1] + DR[0,1,3,1,2] - DR[0,1,1,2,3] # note the change of the second + to - # ### Kretschmann scalar # # The tensor $R^\flat$, of components $R_{abcd} = g_{am} R^m_{\ \, bcd}$: # In[35]: dR = R.down(g) print(dR) # The tensor $R^\sharp$, of components $R^{abcd} = g^{bp} g^{cq} g^{dr} R^a_{\ \, pqr}$: # In[36]: uR = R.up(g) print(uR) # The Kretschmann scalar $K := R^{abcd} R_{abcd}$: # In[37]: Kr_scalar = uR['^{abcd}']*dR['_{abcd}'] Kr_scalar.display() #A variant of this expression can be obtained by invoking the factor() method on the coordinate function representing the scalar field in the manifold's default chart:
# In[38]: Kr = Kr_scalar.coord_function() Kr.factor() #As a check, we can compare Kr to the formula given by R. Conn Henry, Astrophys. J. 535, 350 (2000):
# In[39]: Kr == 48*m^2*(r^6 - 15*r^4*(a*cos(th))^2 + 15*r^2*(a*cos(th))^4 - (a*cos(th))^6) / (r^2+(a*cos(th))^2)^6 #The Schwarzschild value of the Kretschmann scalar is recovered by setting $a=0$:
# In[40]: Kr.expr().subs(a=0) #Let us plot the Kretschmann scalar for $m=1$ and $a=0.9$:
# In[41]: K1 = Kr.expr().subs(m=1, a=0.9) plot3d(K1, (r,1,3), (th, 0, pi), viewer=viewer3D, online=True, axes_labels=['r', 'theta', 'Kr']) # In[ ]: