Notebook

Tutorial: The right-hand sides of the time-evolution equation for a massless scalar field¶

Author(s): Leonardo Werneck and Zach Etienne¶

This tutorial notebook documents and constructs the time evolution equations of the Klein-Gordon equations for a massless scalar field written in terms of the BSSN quantities.¶

Notebook Status: Validated

Validation Notes: The expressions generated by the NRPy+ module corresponding to this tutorial notebook were used to generate the results shown in Werneck et al. (in preparation).

Python module containing the final expressions constructed here: ScalarField/ScalarField_RHSs.py ¶

Table of Contents¶

$$\label{toc}$$

The module is organized as follows

Step 0: The Klein-Gordon equation in curvilinear coordinates
Step 1: Initialize core NRPy+ modules
Step 2: Right-hand side of $\partial_{t} \varphi$
Step 3: Right-hand side of $\partial_{t} \Pi$
Step 4: Code Validation against ScalarField.ScalarField_RHSs.py NRPy+ module
Step 5: Output this module to $\LaTeX$-formatted PDF

Step 0: The Klein-Gordon equation in curvilinear coordinates [Back to top]¶

$$\label{klein_gordon_eq}$$

We will begin by considering the Klein-Gordon equation for a massless scalar field, $\varphi$,

$$ \nabla_{\mu}\left(\nabla^{\mu}\varphi\right) = 0\ . $$

We then define the auxiliary field

$$ \Pi \equiv -\frac{1}{\alpha}\left(\partial_{t}\varphi - \beta^{i}\partial_{i}\varphi\right)\ , $$

so that the Klein-Gordon equation is decomposed into two first order equations (cf. eqs. (5.252) and (5.253) in B&S)

\begin{align} \partial_{t}\varphi &= \beta^{i}\partial_{i}\varphi - \alpha\,\Pi\ ,\nonumber\\ \partial_{t}\Pi &= \beta^{i}\partial_{i}\Pi + \alpha K\,\Pi -\gamma^{ij}\left(\partial_{j}\varphi\,\partial_{i}\alpha -\alpha\,\Gamma^{k}_{\ ij}\, \partial_{k}\varphi+\alpha\,\partial_{i}\partial_{j}\varphi\right)\ , \end{align}

where $K=\gamma^{ij}K_{ij}$ is the trace of the extrinsic curvature $K_{ij}$ and $\gamma^{ij}$ the inverse of the physical spatial metric $\gamma_{ij}$. We will choose not to define the auxiliary variables $\varphi_{i}\equiv\partial_{i}\varphi$ in our code and, instead, leave the equations in terms of second derivatives of $\varphi$.

To write the equations in terms of BSSN variables (see this tutorial notebook for a review) we start by considering the conformal metric, $\bar{\gamma}_{ij}$, related to the physical metric by

$$ \gamma_{ij} = e^{4\phi}\bar{\gamma}_{ij}\ , $$

and its inverse

$$ \gamma^{ij} = e^{-4\phi}\bar{\gamma}^{ij}\ . $$

Let us also look at equation (3.7) of B&S (with $i\leftrightarrow k$ and $\ln\psi = \phi$, for convenience)

$$ \Gamma^{k}_{\ ij} = \bar{\Gamma}^{k}_{\ ij} + 2\left(\delta^{k}_{\ i}\bar{D}_{j}\phi + \delta^{k}_{\ j}\bar{D}_{i}\phi - \bar{\gamma}_{ij}\bar{\gamma}^{k\ell}\bar{D}_{\ell}\phi\right)\ . $$

Then let us consider the term that contains $\Gamma^{k}_{\ ij}$ on the right-hand side of $\partial_{t}\Pi$:

\begin{align} \alpha\,\gamma^{ij}\Gamma^{k}_{\ ij}\, \partial_{k}\varphi &= \alpha e^{-4\phi}\bar{\gamma}^{ij}\left[\bar{\Gamma}^{k}_{\ ij} + 2\left(\delta^{k}_{\ i}\bar{D}_{j}\phi + \delta^{k}_{\ j}\bar{D}_{i}\phi - \bar{\gamma}_{ij}\bar{\gamma}^{k\ell}\bar{D}_{\ell}\phi\right)\right]\partial_{k}\varphi\ . \end{align}

Focusing on the term in parenthesis, we have (ignoring, for now, a few non-essential multiplicative terms and replacing $\bar{D}_{i}\phi = \partial_{i}\phi$)

\begin{align} 2\bar{\gamma}^{ij}\left(\delta^{k}_{\ i}\bar{D}_{j}\phi + \delta^{k}_{\ j}\bar{D}_{i}\phi - \bar{\gamma}_{ij}\bar{\gamma}^{k\ell}\bar{D}_{\ell}\phi\right)\partial_{k}\varphi &= 2\left(\bar{\gamma}^{kj}\partial_{j}\phi + \bar{\gamma}^{ki}\partial_{i}\phi - 3\bar{\gamma}^{k\ell}\partial_{\ell}\phi\right)\partial_{k}\varphi\nonumber\\ &=2\left(\bar{\gamma}^{ij}\partial_{i}\phi + \bar{\gamma}^{ij}\partial_{i}\phi - 3 \bar{\gamma}^{ij}\partial_{i}\phi\right)\partial_{j}\varphi\nonumber\\ &= -2\bar{\gamma}^{ij}\partial_{j}\varphi\partial_{i}\phi\ , \end{align}

so that

$$ \alpha\,\gamma^{ij}\Gamma^{k}_{\ ij}\, \partial_{k}\varphi = e^{-4\phi}\bar{\gamma}^{ij}\left(\alpha \bar{\Gamma}^{k}_{\ ij}\, \partial_{k}\varphi - 2\alpha\partial_{j}\varphi\partial_{i}\phi\right) $$

For the rest of the equation, all we need to do is replace $\gamma^{ij}\to e^{-4\phi}\bar{\gamma}^{ij}$, so that the Klein-Gordon equation becomes

\begin{align} \partial_{t}\Pi = \beta^{i}\partial_{i}\Pi + \alpha K\,\Pi -e^{-4\phi}\bar{\gamma}^{ij}\left(\partial_{j}\varphi\,\partial_{i}\alpha - \alpha\,\bar{\Gamma}^{k}_{\ ij}\,\partial_{k}\varphi + \alpha\,\partial_{i}\partial_{j}\varphi + 2\alpha\,\partial_{j}\varphi\,\partial_{i}\phi\right)\ . \end{align}

Note that the evolution equation for $\varphi$ is left unchanged

$$ \partial_{t}\varphi = \beta^{i}\partial_{i}\varphi - \alpha\,\Pi\ . $$

Step 1: Initialize core NRPy+ modules [Back to top]¶

$$\label{initializenrpy}$$

Let's start by importing all the needed modules from NRPy+:

In [1]:

# Step 1.a: import all needed modules from NRPy+:
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par    # NRPy+: Parameter interface
import grid as gri                # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm    # NRPy+: Reference metric support
import BSSN.BSSN_quantities as Bq # NRPy+: BSSN quantities

# Step 1.b: Set the coordinate system for the numerical grid
coordsystem = "Spherical"
par.set_parval_from_str("reference_metric::CoordSystem",coordsystem)

# Step 1.c: Then we set the theta and phi axes to be the symmetry
#           axes; i.e., axis "1" and "2", corresponding to the
#           i1 and i2 directions.
#           This sets all spatial derivatives in the theta and
#           phi directions to zero (analytically).
# par.set_parval_from_str("indexedexp::symmetry_axes")

# Step 1.d: Given the chosen coordinate system, set up
#           corresponding reference metric and needed
#           reference metric quantities
#    The following function call sets up the reference metric
#    and related quantities, including rescaling matrices ReDD,
#    ReU, and hatted quantities.
rfm.reference_metric()

# Step 1.e: Set spatial dimension (must be 3 for BSSN, as BSSN is
#           a 3+1-dimensional decomposition of the general
#           relativistic field equations)
DIM = 3

# Step 1.f: Import all needed BSSN quantities
Bq.BSSN_basic_tensors()
trK   = Bq.trK
alpha = Bq.alpha
betaU = Bq.betaU
Bq.gammabar__inverse_and_derivs()
gammabarUU = Bq.gammabarUU

Step 2: The right-hand side of $\partial_{t}\varphi$ [Back to top]¶

$$\label{sf_rhs}$$

Let us label each of the terms in the RHS of the $\partial_{t}\varphi$ equation so that it is easier to understand their implementation:

$$ \partial_{t}\varphi = - \underbrace{\alpha \Pi}_{\text{Term 1}} + \underbrace{\beta^{i}\partial_{i}\varphi}_{\text{Term 2}} \ . $$

The first term is an advection term and therefore we need to set up the appropriate derivate. We begin by declaring the grid functions that we need to implement this equation.

A note on notation: We choose to declare the $\color{red}{\textbf{s}}$calar $\color{red}{\textbf{f}}$ield variable, $\varphi$, as $\color{red}{\rm sf}$ and the $\color{blue}{\textbf{s}}$calar $\color{blue}{\textbf{f}}$ield conjugate $\color{blue}{\textbf{M}}$omentum variable, $\Pi$, as $\color{blue}{\rm sfM}$ to avoid possible conflicts with other variables which might be commonly denoted by $\psi$ or $\Pi$.

In [2]:

# Step 2.a: Declare grid functions for varphi and Pi
sf, sfM = gri.register_gridfunctions("EVOL",["sf", "sfM"])

# Step 2.a: Add Term 1 to sf_rhs: -alpha*Pi
sf_rhs = - alpha * sfM

# Step 2.b: Add Term 2 to sf_rhs: beta^{i}\partial_{i}\varphi
sf_dupD = ixp.declarerank1("sf_dupD")
for i in range(DIM):
    sf_rhs += betaU[i] * sf_dupD[i]

Step 3: The right-hand side of $\partial_{t}\Pi$ [Back to top]¶

$$\label{pi_rhs}$$

Now let us (slightly) rearrange the RHS of $\partial_{t}\Pi$ so that we are able to group relevant terms together as we label them

\begin{align} \partial_{t}\Pi &= \underbrace{\alpha K\,\Pi}_{\text{Term 1}} + \underbrace{\beta^{i}\partial_{i}\Pi}_{\text{Term 2}} - \underbrace{e^{-4\phi}\bar{\gamma}^{ij}\left(\partial_{i}\alpha\partial_{j}\varphi+\alpha\partial_{i}\partial_{j}\varphi + 2\alpha\partial_{j}\varphi\partial_{i}\phi\right)}_{\text{Term 3}} + \underbrace{e^{-4\phi}\alpha\bar{\gamma}^{ij}\bar{\Gamma}^{k}_{ij}\partial_{k}\varphi}_{\text{Term 4}}\ . \end{align}

Step 3a: Term 1 of $\partial_{t}\Pi$ = sfM_rhs : $\alpha K\, \Pi$¶

In [3]:

# Step 3a: Adding Term 1 to sfM_rhs: alpha * K * Pi
sfM_rhs = alpha * trK * sfM

Step 3b: Term 2 of $\partial_{t}\Pi$ = sfM_rhs : $\beta^{i}\partial_{i}\Pi$¶

In [4]:

# Step 3b: Adding Term 2 to sfM_rhs: beta^{i}\partial_{i}Pi
sfM_dupD = ixp.declarerank1("sfM_dupD")
for i in range(DIM):
    sfM_rhs += betaU[i] * sfM_dupD[i]

Step 3c: Term 3 of $\partial_{t}\Pi$ = sfM_rhs : $-e^{-4\phi}\bar{\gamma}^{ij}\left(\partial_{i}\alpha\partial_{j}\varphi+\alpha\partial_{i}\partial_{j}\varphi + 2\alpha\partial_{j}\varphi\partial_{i}\phi\right)$¶

Now let's focus on Term 3:

$$ e^{-4\phi}\left(-\underbrace{\bar{\gamma}^{ij}\partial_{i}\alpha\partial_{j}\varphi}_{\text{Term 3a}}-\underbrace{\alpha\bar{\gamma}^{ij}\partial_{i}\partial_{j}\varphi}_{\text{Term 3b}} - \underbrace{2\alpha\bar{\gamma}^{ij}\partial_{j}\varphi\partial_{i}\phi}_{\text{Term 3c}}\right) $$

In [5]:

# Step 3c: Adding Term 3 to sfM_rhs
# Step 3c.i: Term 3a: gammabar^{ij}\partial_{i}\alpha\partial_{j}\varphi
alpha_dD    = ixp.declarerank1("alpha_dD")
sf_dD       = ixp.declarerank1("sf_dD")
sfMrhsTerm3 = sp.sympify(0)
for i in range(DIM):
    for j in range(DIM):
        sfMrhsTerm3 += - gammabarUU[i][j] * alpha_dD[i] * sf_dD[j]

# Step 3c.ii: Term 3b: \alpha*gammabar^{ij}\partial_{i}\partial_{j}\varphi
sf_dDD = ixp.declarerank2("sf_dDD","sym01")
for i in range(DIM):
    for j in range(DIM):
        sfMrhsTerm3 += - alpha * gammabarUU[i][j] * sf_dDD[i][j]

# Step 3c.iii: Term 3c: 2*alpha*gammabar^{ij}\partial_{j}\varphi\partial_{i}\phi
Bq.phi_and_derivs() # sets exp^{-4phi} = exp_m4phi and \partial_{i}phi = phi_dD[i]
for i in range(DIM):
    for j in range(DIM):
        sfMrhsTerm3 += - 2 * alpha * gammabarUU[i][j] * sf_dD[j] * Bq.phi_dD[i]

# Step 3c.iv: Multiplying Term 3 by e^{-4phi} and adding it to sfM_rhs
sfMrhsTerm3 *= Bq.exp_m4phi
sfM_rhs     += sfMrhsTerm3

Step 3d: Term 4 of $\partial_{t}\Pi$ = sfM_rhs : $e^{-4\phi}\alpha\bar{\gamma}^{ij}\bar{\Gamma}^{k}_{\ ij}\partial_{k}\varphi$¶

We are now going to rewrite this term a bit before implementation. This rewriting is useful in order to reduce the number of finite differences approximations used to evaluate this term. Let us consider the following definitions (see equations 12a and 12b and the discussion above equation 15 in Brown (2009))

\begin{align} \Delta\Gamma^{k}_{\ ij} &\equiv \bar\Gamma^{k}_{\ ij} - \hat\Gamma^{k}_{\ ij}\ ,\\ \bar\Lambda^{k} &\equiv \bar\gamma^{ij}\Delta\Gamma^{k}_{\ ij}\ . \end{align}

Then we have

\begin{align} \text{Term 4} &= e^{-4\phi}\alpha\bar{\gamma}^{ij}\bar{\Gamma}^{k}_{\ ij}\partial_{k}\varphi\nonumber\\ &=e^{-4\phi}\alpha\bar\gamma^{ij}\Delta\Gamma^{k}_{\ ij}\partial_{k}\varphi + e^{-4\phi}\alpha\bar{\gamma}^{ij}\hat{\Gamma}^{k}_{\ ij}\partial_{k}\varphi\nonumber\\ &=e^{-4\phi}\left(\underbrace{\alpha\bar\Lambda^{i}\partial_{i}\varphi}_{\text{Term 4a}} + \underbrace{\alpha\bar{\gamma}^{ij}\hat{\Gamma}^{k}_{\ ij}\partial_{k}\varphi}_{\text{Term 4b}}\right) \end{align}

In [6]:

# Step 3d: Adding Term 4 to sfM_rhs
# Step 3d.i: Term 4a: \alpha \bar\Lambda^{i}\partial_{i}\varphi
LambdabarU  = Bq.LambdabarU
sfMrhsTerm4 = sp.sympify(0)
for i in range(DIM):
    sfMrhsTerm4 += alpha * LambdabarU[i] * sf_dD[i]

# Step 3d.ii: Evaluating \bar\gamma^{ij}\hat\Gamma^{k}_{ij}
GammahatUDD = rfm.GammahatUDD
gammabarGammahatContractionU = ixp.zerorank1()
for k in range(DIM):
    for i in range(DIM):
        for j in range(DIM):
            gammabarGammahatContractionU[k] += gammabarUU[i][j] * GammahatUDD[k][i][j]

# Step 3d.iii: Term 4b: \alpha \bar\gamma^{ij}\hat\Gamma^{k}_{ij}\partial_{k}\varphi
for i in range(DIM):
    sfMrhsTerm4 += alpha * gammabarGammahatContractionU[i] * sf_dD[i]

# Step 3d.iii: Multplying Term 4 by e^{-4phi} and adding it to sfM_rhs
sfMrhsTerm4 *= Bq.exp_m4phi
sfM_rhs     += sfMrhsTerm4

Step 4: Code Validation against ScalarField.ScalarField_RHSs.py NRPy+ module [Back to top]¶

$$\label{code_validation}$$

Here we perform a code validation. We verify agreement in the SymPy expressions for the RHSs of the scalar field equations between

this tutorial notebook and
the ScalarField.ScalarField_RHSs.py NRPy+ module.

By default, we analyze the RHSs in Spherical coordinates, though other coordinate systems may be chosen.

In [7]:

# Step 4: Code validation against NRPy+ module
# Step 4.a: Import the ScalarFieldCollapse module and
#           run the ScalarFieldCollapse.scalarfield_RHSs()
#           function to evaluate the RHSs.
import ScalarField.ScalarField_RHSs as sfrhs # NRPyCritCol: Scalar field right-hand sides
sfrhs.ScalarField_RHSs()

# Step 4.b: Perform the consistency check by subtracting
#           the RHSs computed in this tutorial from the
#           ones computed in the ScalarFieldCollapse module
print("Consistency check between Scalar Field RHSs tutorial and NRPy+ module")
print(" sf_rhs - sfrhs.sf_rhs  = " + str( sf_rhs - sfrhs.sf_rhs )+" Should be zero.")
print("sfM_rhs - sfrhs.sfM_rhs = " + str(sfM_rhs - sfrhs.sfM_rhs)+" Should be zero.")

Consistency check between Scalar Field RHSs tutorial and NRPy+ module
 sf_rhs - sfrhs.sf_rhs  = 0 Should be zero.
sfM_rhs - sfrhs.sfM_rhs = 0 Should be zero.

Step 5: Output this module to $\LaTeX$-formatted PDF file [Back to top]¶

$$\label{latex_pdf_output}$$

The following code cell converts this Jupyter notebook into a proper, clickable $\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename Tutorial-ScalarField_RHSs.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

In [8]:

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-ScalarField_RHSs")

Created Tutorial-ScalarField_RHSs.tex, and compiled LaTeX file to PDF file
    Tutorial-ScalarField_RHSs.pdf