metric_and_scalar_equations_cubic_Galileon.ipynb

Kernel: SageMath (default)

Metric and scalar equations for the cubic Galileon in quasi-isotropic coordinates

version()

'SageMath version 8.9, Release Date: 2019-09-29'

Set up the notebook to display mathematical objects using LaTeX formatting:

%display latex

Setup

Declare a 4-dimensional manifold $M$:

M = Manifold(4, 'M')

Declare a quasi-isotropic chart over $M$:

QI.<t,r,theta,p> = M.chart()

Declare the four metric functions over $M$:

N = function('N')(r,theta)
A = function('A')(r,theta)
B = function('B')(r,theta)
w = function('omega')(r,theta)

Declare the metric $g$ on $M$:

g = M.metric('g')

First set all its components to zero:

g.set_comp()

Then set the non zero components according to relation (8):

g[0,0] = - N^2 + B^2 * r^2 * sin(theta)^2 * w^2
g[3,0] = - B^2 * r^2 * sin(theta)^2 * w
g[1,1] = A^2
g[2,2] = A^2 * r^2
g[3,3] = B^2 * r^2 * sin(theta)^2

Compute the inverse metric for later use:

g_inv = g.inverse()

Compute the associated Levi-Civita connection for later use:

NAB = g.connection()

To define the scalar field $\phi$ over $M$ according to the ansatz (11), first declare the parameter $q$:

q = var('q')

Then declare the spatial dependece $\Psi$:

Psi = function('Psi')(r,theta)

Finally, declare the scalar field $\phi$ over $M$ according to the ansatz (11):

phi = M.scalar_field(q*t + Psi)

Compute the terms of tensor $T^{(\phi)}$ given in relation (3)

First declare the coupling constants $\zeta$, $\eta$ and $\gamma$:

[zeta, eta, gamma] = var('zeta eta gamma')

Term proportional to $\eta$

To compute the term proportional to $\eta$, first compute the exterior derivative of $\phi$:

d_phi = phi.differential()

Then compute its squared norm $(\partial\phi)^{2} = \nabla_{\mu}\phi\nabla^{\mu}\phi$ using the inverse metric:

NAB_phi_SQR = ( g_inv['^u.'] * d_phi['_u'] )['^v'] * d_phi['_v']

Finally, compute the term proportional to $\eta$:

eta_term = eta/zeta * ( d_phi * d_phi - 1/2 * g * NAB_phi_SQR )

Term proportional to $\gamma$

To compute the term proportional to $\gamma$, first compute the exterior derivative of NAB_phi_SQR:

d_NAB_phi_SQR = NAB_phi_SQR.differential()

Then compute the first term as the symmetrized tensor product of d_NAB_phi_SQR with $d\phi$:

gamma_term_1 = 1/2 * ( d_phi * d_NAB_phi_SQR + d_NAB_phi_SQR * d_phi )

To compute the second term, first compute $\Box\phi$ using the Hessian matrix of $\phi$:

HESS_phi = NAB(d_phi)

Extract $\Box\phi$ as the trace of HESS_phi:

BOX_phi = g_inv['^uv'] * HESS_phi['_uv']

Then compute the second term:

gamma_term_2 = - BOX_phi * d_phi * d_phi

And the third term:

gamma_term_3 = - 1/2 * g * ( d_phi.up(g)['^u'] * d_NAB_phi_SQR['_u'] )

Finally, add all three terms:

gamma_term = gamma/zeta * ( gamma_term_1 + gamma_term_2 + gamma_term_3 )

3+1 decomposition of the effective stress-energy-momentum tensor $T^{(eff)}$

Subtract the cosmological term to $T^{(\phi)}$ to form $T^{(eff)}$ in order to write the metric equations as $G_{\mu\nu} = 8\pi T^{(eff)}_{\mu\nu}$. First declare a cosmological constant:

Lambda = var('Lambda')

Then declare $T^{(eff)}$:

T_eff = 1/(8*pi)* (eta_term + gamma_term - Lambda * g)

It is shown in ref. [87] (E. Gourgoulhon, "An introduction to the theory of rotating relativistic stars") that algebraic manipulations of equations of the form $G_{\mu\nu} = 8\pi T^{(eff)}_{\mu\nu}$ written in quasi-isotropic coordinates may result in equations (3.14) to (3.17) of the same reference. These equations rewrite as $$ N^{2} \Delta_{3} N = \frac{N \left( Br\sin\theta \right)^{2}}{2} \partial\omega \partial\omega

• \frac{N^{2}}{B} \partial N \partial B

• 4 \pi A^{2} N^{3} (E + S), \ N^{3} \Delta_{2} [NA] = \frac{N^4}{A}\partial A \partial A

• 2 N^3 \partial A \partial N

• \frac{3 A (N Br\sin\theta)^{2}}{4} \partial\omega \partial\omega

• 8 \pi A^{3} N^{4} S^{\varphi}{\phantom{\varphi}\varphi}, \ N^{2} \Delta{2} [NBr\sin\theta] = 8 \pi A^{2} B N^{3} r \sin\theta (S^{r}{\phantom{r}r} + S^{\theta}{\phantom{\theta}\theta}), \ N \Delta_{3} [\omega r \sin \theta] = \frac{N \omega}{r \sin\theta}

• r \sin\theta \left( \partial\omega \partial N - \frac{3 N}{B} \partial \omega \partial B \right), $$ where $E$and the$S^{i}_{\phantom{i}j}$are defined from the 3+1 decomposition of$T^{(eff)}$ performed below.

Declare the time function over $M$:

time = M.scalar_field(t)

Declare the future-oriented unit vector normal to the constant-$t$ hypersurfaces, which is constructed as the metric dual of the exterior derivative of the time function, normalized by the metric function $N$:

unit_normal = -N * time.differential().up(g)

Effective energy density $E$

Compute the quantity $E$ appearing in the equations above as the full contraction of $T^{(eff)}$ with the unit normal vector:

E = (T_eff['_u.'] * unit_normal['^u'])['_v'] * unit_normal['^v']

Extract its expression, automatically simplified for latter use:

E = E.expr().simplify_full()

Effective stress tensor $S^{(eff)}$

To construct the projector onto the tangent spaces to the constant-$t$ hypersurfaces, first declare the identity map on the full tangent spaces to $M$:

identity = M.tensor_field(1,1)

First set all its components to zero:

identity.set_comp()

And set to $1$ the diagonal components:

for mu in range(M.dim()):
    identity[mu,mu] = 1

Declare the projector from the identity and the tensor product of the unit normal with its metric dual:

projector = identity + unit_normal * unit_normal.down(g)

Declare the effective stress tensor $S^{(eff)}$ from two contractions of $T^{(eff)}$ with the projector:

S_eff = T_eff['_uv'] * (projector * projector)['^uv_..']

Extract the simplified expressions of the quantities $S^{r}_{\phantom{r}r}$, $S^{\theta}_{\phantom{\theta}\theta}$ and $S^{\varphi}_{\phantom{\varphi}\varphi}$:

S_r_up_r_down = S_eff.up(g,pos=0)[1,1].expr().simplify_full()
S_theta_up_theta_down = S_eff.up(g,pos=0)[2,2].expr().simplify_full()
S_phi_up_phi_down = S_eff.up(g,pos=0)[3,3].expr().simplify_full()

Extract the trace $S$ of $S^{(eff)}$:

S = (S_r_up_r_down + S_theta_up_theta_down + S_phi_up_phi_down).simplify_full()

Metric equations (21) - (24) with sources (B.5) - (B.8)

To check the metric equations, first define the various differential operators involved. First is the 2-dimensional Laplacian denoted $\Delta_{2}$ in the equations above:

def Lap2(f,r,theta):
    return diff(diff(f,r),r) + 1/r*diff(f,r) + 1/r^2*diff(diff(f,theta),theta)

Then the 3-dimensional Laplacian denoted $\Delta_{3}$:

def Lap3(f,r,theta):
    return diff(diff(f,r),r) + 2/r*diff(f,r) + 1/r^2*diff(diff(f,theta),theta) + 1/(r^2*tan(theta))*diff(f,theta)

Then the operator acting on any two functions $f$ and $g$ of the coordinates $r$ and $\theta$ and denoted $\partial f \partial g$:

def DD(f,g,r,theta):
    return diff(f,r)*diff(g,r) + 1/r^2*diff(f,theta)*diff(g,theta)

Then the operator acting on any three functions $f$, $g$ and $h$ of the coordinates $r$ and $\theta$ and denoted $\mathcal{H}^{(0)}_{f}[g,h]$:

def Hf0gh(f,g,h,r,theta):
    left = matrix(SR, 2, 1)
    right = matrix(SR, 2, 1)
    mat = matrix(SR,2,2)
    left[0,0] = diff(g,r)
    left[1,0] = 1/r*diff(g,theta)
    right[0,0] = diff(h,r)
    right[1,0] = 1/r*diff(h,theta)
    mat[0,0] = diff(diff(f,r),r)
    mat[0,1] = 1/r*diff(diff(f,r),theta)
    mat[1,0] = mat[0,1]
    mat[1,1] = 1/r^2 * diff(diff(f,theta),theta)
    return (left.transpose()*mat*right)[0,0].simplify_full()

Metric equation (21) with source term (B.5)

Check that the last term appearing in the right-hand side of the equation on $\Delta_{3} N$ above, i.e. $4 \pi A^{2} N^{3} (E + S)$, explicitly rewrites as

$$4 \pi A^{2} N^{3} (E + S) = N A^{2} \left( \frac{\eta q^{2}}{\zeta} - \Lambda N^{2} \right) - \frac{\gamma}{2 \zeta} \left(q^{2} + \frac{N^{2} \partial\Psi\partial\Psi}{A^{2}} \right) \left( N \Delta_{3}\Psi + \partial\Psi\partial N + \frac{N \partial\Psi\partial B}{B} \right).$$

To do so, declare the right-hand side of this relation:

Source_N_last_terms_1 = N*A^2 * ( eta*q^2/zeta - Lambda*N^2 )
Source_N_last_terms_2a = -gamma/(2*zeta) * ( q^2 + N^2*DD(Psi,Psi,r,theta)/A^2 )
Source_N_last_terms_2b = N*Lap3(Psi,r,theta) + DD(Psi,N,r,theta) + N*DD(Psi,B,r,theta)/B
Source_N_last_terms = Source_N_last_terms_1 + Source_N_last_terms_2a * Source_N_last_terms_2b  

Check that it coincides with $4 \pi A^{2} N^{3} (E + S)$:

( Source_N_last_terms - 4*pi*A^2*N^3*(E + S) ).simplify_full()

Therefore the equation on $\Delta_{3} N$ explicitly rewrites as

$$N^{2} \Delta_{3} N = \frac{N \left( Br\sin\theta \right)^{2}}{2} \partial\omega \partial\omega - \frac{N^{2}}{B} \partial N \partial B + N A^{2} \left( \frac{\eta q^{2}}{\zeta} - \Lambda N^{2} \right) - \frac{\gamma}{2 \zeta} \left(q^{2} + \frac{N^{2} \partial\Psi\partial\Psi}{A^{2}} \right) \left( N \Delta_{3}\Psi + \partial\Psi\partial N + \frac{N \partial\Psi\partial B}{B} \right).$$

Using the relations (12)-(13) to reformulate everything in terms of dimensionless quantities eventually leads to equation (21) with source (B.5).

Metric equation (22) with source term (B.6)

Check that the last term appearing in the right-hand side of the equation on $\Delta_{2} [NA]$ above, i.e. $8 \pi A^{3} N^{4} S^{\varphi}_{\phantom{\varphi}\varphi}$, explicitly rewrites as

$$8 \pi A^{3} N^{4} S^{\varphi}_{\phantom{\varphi}\varphi} = \frac{\eta N^{2} A}{2\zeta} \left( q^{2} A^{2} - N^{2} \partial\Psi\partial\Psi \right) - \Lambda A^{3} N^{4} -\frac{\gamma}{\zeta} \left( q^{2} N A \partial\Psi\partial N - \frac{N^{4} \partial\Psi\partial\Psi \partial\Psi\partial A}{A^{2}} + \frac{1}{A} \left[ N^{4} \mathcal{H}^{(0)}_{\Psi}[\Psi,\Psi] - \frac{N^{4} \partial_{r}\Psi}{r^{3}} (\partial_{\theta}\Psi)^{2} \right] \right).$$

To do so, declare the right-hand side of this relation:

Source_A_last_terms_1 = eta*N^2*A/(2*zeta)*(q^2*A^2 - N^2*DD(Psi,Psi,r,theta))
Source_A_last_terms_2 = - Lambda*A^3*N^4
Source_A_last_terms_3a = q^2*N*A*DD(Psi,N,r,theta) - N^4*DD(Psi,Psi,r,theta)*DD(Psi,A,r,theta)/A^2
Source_A_last_terms_3b = 1/A*(N^4*Hf0gh(Psi,Psi,Psi,r,theta) - N^4*diff(Psi,r)/r^3*diff(Psi,theta)^2)
Source_A_last_terms = Source_A_last_terms_1 + Source_A_last_terms_2 - gamma/zeta*(Source_A_last_terms_3a + Source_A_last_terms_3b)

Check that it coincides with $8 \pi A^{3} N^{4} S^{\varphi}_{\phantom{\varphi}\varphi}$:

( Source_A_last_terms - 8*pi*A^3*N^4*S_phi_up_phi_down ).simplify_full()

Therefore the equation on $\Delta_{2} [NA]$ explicitly rewrites as

$$N^{3} \Delta_{2} [NA] = \frac{N^4}{A}\partial A \partial A + 2 N^3 \partial A \partial N + \frac{3 A (N Br\sin\theta)^{2}}{4} \partial\omega \partial\omega + \frac{\eta N^{2} A}{2\zeta} \left( q^{2} A^{2} - N^{2} \partial\Psi\partial\Psi \right) - \Lambda A^{3} N^{4} -\frac{\gamma}{\zeta} \left( q^{2} N A \partial\Psi\partial N - \frac{N^{4} \partial\Psi\partial\Psi \partial\Psi\partial A}{A^{2}} + \frac{1}{A} \left[ N^{4} \mathcal{H}^{(0)}_{\Psi}[\Psi,\Psi] - \frac{N^{4} \partial_{r}\Psi}{r^{3}} (\partial_{\theta}\Psi)^{2} \right] \right).$$

Using the relations (12)-(13) to reformulate everything in terms of dimensionless quantities eventually leads to equation (22) with source (B.6).

Metric equation (23) with source term (B.7)

Check that the last term appearing in the right-hand side of the equation on $\Delta_{2} [NBr\sin\theta]$ above, i.e. $8 \pi A^{2} B N^{3} r \sin\theta (S^{r}_{\phantom{r}r} + S^{\theta}_{\phantom{\theta}\theta})$, explicitly rewrites as

$$8 \pi A^{2} B N^{3} r \sin\theta (S^{r}_{\phantom{r}r} + S^{\theta}_{\phantom{\theta}\theta}) = - Br\sin\theta \left[ N A^{2} \left( 2 \Lambda N^{2} - \frac{q^{2} \eta}{\zeta} \right) + \frac{\gamma N^{2} \partial\Psi\partial\Psi}{\zeta A^{2}} \left( N \Delta_{3}\Psi + \partial\Psi\partial N + \frac{N \partial\Psi\partial B}{B} \right) \right].$$

To do so, declare the right-hand side of this relation:

Source_B_last_terms_1 = N*A^2*(2*Lambda*N^2 - q^2*eta/zeta)
Source_B_last_terms_2a = gamma*N^2*DD(Psi,Psi,r,theta)/(zeta*A^2)
Source_B_last_terms_2b = (N*Lap3(Psi,r,theta) + DD(Psi,N,r,theta) + N*DD(Psi,B,r,theta)/B)
Source_B_last_terms = -B*r*sin(theta)*( Source_B_last_terms_1 + Source_B_last_terms_2a*Source_B_last_terms_2b )

Check that it coincides with $8 \pi A^{2} B N^{3} r \sin\theta (S^{r}_{\phantom{r}r} + S^{\theta}_{\phantom{\theta}\theta})$:

( Source_B_last_terms - 8*pi*A^2*B*N^3*r*sin(theta)*(S_r_up_r_down + S_theta_up_theta_down) ).simplify_full()

Therefore the equation on $\Delta_{2} [NBr\sin\theta]$ explicitly rewrites as

$$N^{2} \Delta_{2} [NBr\sin\theta] = - Br\sin\theta \left[ N A^{2} \left( 2 \Lambda N^{2} - \frac{q^{2} \eta}{\zeta} \right) + \frac{\gamma N^{2} \partial\Psi\partial\Psi}{\zeta A^{2}} \left( N \Delta_{3}\Psi + \partial\Psi\partial N + \frac{N \partial\Psi\partial B}{B} \right) \right].$$

Using the relations (12)-(13) to reformulate everything in terms of dimensionless quantities eventually leads to equation (23) with source (B.7).

Finally, it is also only needed to use the relations (12)-(13) in the equation on $N \Delta_{3} [\omega r \sin \theta]$ above to obtain (24) with source (B.8).

Scalar equation (25) with explicit expression (B.9)

Declare the covariant expression of the vector current $J$ given by relation (5):

one_form_J = ( gamma * BOX_phi - eta ) * d_phi - gamma/2 * d_NAB_phi_SQR

Declare the Jacobian matrix of $J$:

JAC_J = NAB(one_form_J.up(g))

Extract the divergence of $J$ as the trace of the Jacobian matrix:

scalar_eq = JAC_J.trace()

Then define the various differential operators involved in the explicit expression of the scalar equation:

def Hf1gh(f,g,h,r,theta):
    left = matrix(SR, 2, 1)
    right = matrix(SR, 2, 1)
    mat = matrix(SR,2,2)
    left[0,0] = 1/r*diff(g,theta)
    left[1,0] = -diff(g,r)
    right[0,0] = 1/r*diff(h,theta)
    right[1,0] = -diff(h,r)
    mat[0,0] = diff(diff(f,r),r)
    mat[0,1] = 1/r*diff(diff(f,r),theta)
    mat[1,0] = mat[0,1]
    mat[1,1] = 1/r^2 * diff(diff(f,theta),theta)
    return (left.transpose()*mat*right)[0,0].simplify_full()

def Hf2gh(f,g,h,r,theta):
    left = matrix(SR, 2, 1)
    right = matrix(SR, 2, 1)
    mat = matrix(SR,2,2)
    left[0,0] = diff(g,r)
    left[1,0] = 1/r*diff(g,theta)
    right[0,0] = diff(h,r)
    right[1,0] = 1/r*diff(h,theta)
    mat[0,0] = diff(diff(f,r),r)
    mat[0,1] = 2/r*diff(diff(f,r),theta)
    mat[1,0] = mat[0,1]
    mat[1,1] = 1/r^2 * diff(diff(f,theta),theta)
    return (left.transpose()*mat*right)[0,0].simplify_full()

Declare all the terms involved in the explicit expression of the scalar equation:

term_1 = Lap3(N,r,theta) + DD(N, B,r,theta)/B - 2*DD(N, N,r,theta)/N
term_1 *= q^2 * A^2 * N^(-3) 

term_2 = Lap3(Psi,r,theta) + DD(Psi, B,r,theta)/B + DD(Psi, N,r,theta)/N
term_2 *= gamma^(-1) * eta * A^2 + 2 * DD(Psi, A,r,theta)/A - 2 * DD(Psi, N,r,theta)/N

term_3 = - 2 * (   Lap2(Psi,r,theta) + DD(Psi,B,r,theta)/B   )
term_3 *= Lap3(Psi,r,theta) - 1/r*diff(Psi,r)

term_4 = 2/r^2*diff(diff(Psi,theta),theta) 
term_4 *= Lap2(Psi,r,theta) - DD(Psi,A,r,theta)/A

term_5 = A^(-1) * (Lap3(A,r,theta) - 4/r*diff(A,r))
term_5 += B^(-1) * Lap2(B,r,theta)
term_5 += DD(A, B,r,theta)/(A*B)
term_5 -= 3 * DD(A, A,r,theta)/A^2
term_5 += DD(N, A,r,theta)/(N*A)
term_5 *= - DD(Psi,Psi,r,theta)

term_6 = 2 * (DD(Psi,N,r,theta)/N)^2

term_7 = - N^(-1) * Hf0gh(N,Psi,Psi,r,theta)

term_8 = B^(-1) * Hf1gh(B,Psi,Psi,r,theta)

term_9 = - 2 * A^(-1) * Hf2gh(Psi,Psi,A,r,theta)

term_10 = 2 * (   diff(diff(Psi,r),r)^2 + ( 1/r^2*diff(Psi,theta) - 1/r*diff(diff(Psi,r),theta) )^2   )

term_11 = - 2 * diff(Psi,r) * diff(A,r)/A
term_11 *= diff(diff(Psi,r),r) + 2/r*diff(Psi,r) - 1/r^2*diff(diff(Psi,theta),theta)

term_12 = - 1/r * diff(Psi,r)^2 * diff(B,r)/B

term_13 = 1/r^3 * diff(Psi,theta) * (   2*diff(Psi,r)*diff(N,theta)/N - diff(Psi,theta)*diff(N,r)/N   )

explicit_scalar_eq = term_1 + term_2 + term_3 + term_4 + term_5 + term_6 + term_7 + term_8 + term_9 + term_10 + term_11 + term_12 + term_13
explicit_scalar_eq *= -gamma * A^(-4)

Check that both expressions of the scalar equation coincide:

(scalar_eq.expr() - explicit_scalar_eq).simplify_full()

The relations (12)-(13) finally allow to rewrite the explicit expression as (B.9).