BSSN Quantities¶

Author: Zach Etienne¶

Formatting improvements courtesy Brandon Clark¶

This module documents and constructs a number of quantities useful for building symbolic (SymPy) expressions in terms of the core BSSN quantities $\left\{h_{i j},a_{i j},\phi, K, \lambda^{i}, \alpha, \mathcal{V}^i, \mathcal{B}^i\right\}$, as defined in Ruchlin, Etienne, and Baumgarte (2018) (see also Baumgarte, Montero, Cordero-Carrión, and Müller (2012)).¶

Notebook Status: Self-Validated

Validation Notes: This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented below. Additional validation tests may have been performed, but are as yet, undocumented. (TODO)

A Note on Notation:¶

As is standard in NRPy+,

• Greek indices refer to four-dimensional quantities where the zeroth component indicates temporal (time) component.
• Latin indices refer to three-dimensional quantities. This is somewhat counterintuitive since Python always indexes its lists starting from 0. As a result, the zeroth component of three-dimensional quantities will necessarily indicate the first spatial direction.

As a corollary, any expressions involving mixed Greek and Latin indices will need to offset one set of indices by one: A Latin index in a four-vector will be incremented and a Greek index in a three-vector will be decremented (however, the latter case does not occur in this tutorial notebook).

Table of Contents¶

$$\label{toc}$$

Each family of quantities is constructed within a given function (boldfaced below). This notebook is organized as follows

1. Step 1: Initialize needed Python/NRPy+ modules
2. Step 2: declare_BSSN_gridfunctions_if_not_declared_already(): Declare all of the core BSSN variables $\left\{h_{i j},a_{i j},\text{cf}, K, \lambda^{i}, \alpha, \mathcal{V}^i, \mathcal{B}^i\right\}$ and register them as gridfunctions
3. Step 3 Rescaling tensors to avoid coordinate singularities
   1. Step 3.a BSSN_basic_tensors(): Define all basic conformal BSSN tensors $\left\{\bar{\gamma}_{i j},\bar{A}_{i j},\bar{\Lambda}^{i}, \beta^i, B^i\right\}$ in terms of BSSN gridfunctions
4. Step 4: gammabar__inverse_and_derivs(): $\bar{\gamma}^{ij}$, and spatial derivatives of $\bar{\gamma}_{ij}$ including $\bar{\Gamma}^{i}_{jk}$
   1. Step 4.a: Inverse conformal 3-metric: $\bar{\gamma^{ij}}$
   2. Step 4.b: Derivatives of the conformal 3-metric $\bar{\gamma}_{ij,k}$ and $\bar{\gamma}_{ij,kl}$, and associated "barred" Christoffel symbols $\bar{\Gamma}^{i}_{jk}$
5. Step 5: detgammabar_and_derivs(): $\det \bar{\gamma}_{ij}$ and its derivatives
6. Step 6: AbarUU_AbarUD_trAbar(): Quantities related to conformal traceless extrinsic curvature $\bar{A}_{ij}$: $\bar{A}^{ij}$, $\bar{A}^i_j$, and $\bar{A}^k_k$
7. Step 7: RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU(): The conformal ("barred") Ricci tensor $\bar{R}_{ij}$ and associated quantities
   1. Step 7.a: Conformal Ricci tensor, part 1: The $\hat{D}_{k} \hat{D}_{l} \bar{\gamma}_{i j}$ term
   2. Step 7.b: Conformal Ricci tensor, part 2: The $\bar{\gamma}_{k(i} \hat{D}_{j)} \bar{\Lambda}^{k}$ term
   3. Step 7.c: Conformal Ricci tensor, part 3: The $\Delta^{k} \Delta_{(i j) k} + \bar{\gamma}^{k l} \left (2 \Delta_{k(i}^{m} \Delta_{j) m l} + \Delta_{i k}^{m} \Delta_{m j l} \right )$ terms
   4. Step 7.d: Summing the terms and defining $\bar{R}_{ij}$
8. Step 8: betaU_derivs(): Unrescaled shift vector $\beta^i$ and spatial derivatives $\beta^i_{,j}$ and $\beta^i_{,jk}$
9. Step 9: phi_and_derivs(): Standard BSSN conformal factor $\phi$, and its derivatives $\phi_{,i}$, $\phi_{,ij}$, $\bar{D}_j \phi$, and $\bar{D}_j\bar{D}_k \phi$
   1. Step 9.a: $\phi$ in terms of the chosen (possibly non-standard) conformal factor variable cf (e.g., cf$=W=e^{-4\phi}$)
   2. Step 9.b: Partial and covariant derivatives of $\phi$
10. Step 10: Code Validation against BSSN.BSSN_quantities NRPy+ module
11. Step 11: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: Initialize needed Python/NRPy+ modules [Back to top]¶

$$\label{initializenrpy}$$

Step 2: declare_BSSN_gridfunctions_if_not_declared_already(): Declare all of the core BSSN variables $\left\{h_{i j},a_{i j},\text{cf}, K, \lambda^{i}, \alpha, \mathcal{V}^i, \mathcal{B}^i\right\}$ and register them as gridfunctions [Back to top]¶

$$\label{declare_bssn_gfs}$$

Step 3: Rescaling tensors to avoid coordinate singularities [Back to top]¶

$$\label{rescaling_tensors}$$

While the covariant form of the BSSN evolution equations are properly covariant (with the potential exception of the shift evolution equation, since the shift is a freely specifiable gauge quantity), components of the rank-1 and rank-2 tensors $\varepsilon_{i j}$, $\bar{A}_{i j}$, and $\bar{\Lambda}^{i}$ will drop to zero (destroying information) or diverge (to $\infty$) at coordinate singularities.

The good news is, this singular behavior is well-understood in terms of the scale factors of the reference metric, enabling us to define rescaled version of these quantities that are well behaved (so that, e.g., they can be finite differenced).

For example, given a smooth vector in a 3D Cartesian basis $\bar{\Lambda}^{i}$, all components $\bar{\Lambda}^{x}$, $\bar{\Lambda}^{y}$, and $\bar{\Lambda}^{z}$ will be smooth (by assumption). When changing the basis to spherical coordinates (applying the appropriate Jacobian matrix transformation), we will find that since $\phi = \arctan(y/x)$, $\bar{\Lambda}^{\phi}$ is given by

\begin{align} \bar{\Lambda}^{\phi} &= \frac{\partial \phi}{\partial x} \bar{\Lambda}^{x} + \frac{\partial \phi}{\partial y} \bar{\Lambda}^{y} + \frac{\partial \phi}{\partial z} \bar{\Lambda}^{z} \\ &= -\frac{y}{x^2+y^2} \bar{\Lambda}^{x} + \frac{x}{x^2+y^2} \bar{\Lambda}^{y} \\ &= -\frac{y}{(r \sin\theta)^2} \bar{\Lambda}^{x} + \frac{x}{(r \sin\theta)^2} \bar{\Lambda}^{y}. \end{align}

Thus $\bar{\Lambda}^{\phi}$ diverges at all points where $r\sin\theta=0$ (or equivalently where $x=y=0$; i.e., the $z$-axis) due to the $\frac{1}{(r\sin\theta)^2}$ that appear in the Jacobian transformation.

This divergence might pose no problem on cell-centered grids that avoid $r \sin\theta=0$, except that the BSSN equations require that first and second derivatives of these quantities be taken. Usual strategies for numerical approximation of these derivatives (e.g., finite difference methods) will "see" these divergences and errors generally will not drop to zero with increased numerical sampling of the functions at points near where the functions diverge.

However, notice that if we define $\lambda^{\phi}$ such that

$$\bar{\Lambda}^{\phi} = \frac{1}{r\sin\theta} \lambda^{\phi},$$

then $\lambda^{\phi}$ will be smooth as well.

Avoiding such singularities can be generalized to other coordinate systems, so long as $\lambda^i$ is defined as:

$$\bar{\Lambda}^{i} = \frac{\lambda^i}{\text{scalefactor[i]}} ,$$

where scalefactor[i] is the $i$th scale factor in the given coordinate system. In an identical fashion, we define the smooth versions of $\beta^i$ and $B^i$ to be $\mathcal{V}^i$ and $\mathcal{B}^i$, respectively. We refer to $\mathcal{V}^i$ and $\mathcal{B}^i$ as vet[i] and bet[i] respectively in the code after the Hebrew letters that bear some resemblance.

Similarly, we define the smooth versions of $\bar{A}_{ij}$ and $\varepsilon_{ij}$ ($a_{ij}$ and $h_{ij}$, respectively) via

\begin{align} \bar{A}_{ij} &= \text{scalefactor[i]}\ \text{scalefactor[j]}\ a_{ij} \\ \varepsilon_{ij} &= \text{scalefactor[i]}\ \text{scalefactor[j]}\ h_{ij}, \end{align}

where in this case we multiply due to the fact that these tensors are purely covariant (as opposed to contravariant). To slightly simplify the notation, in NRPy+ we define the rescaling matrices ReU[i] and ReDD[i][j], such that

\begin{align} \text{ReU[i]} &= 1 / \text{scalefactor[i]} \\ \text{ReDD[i][j]} &= \text{scalefactor[i] scalefactor[j]}. \end{align}

Thus, for example, $\bar{A}_{ij}$ and $\bar{\Lambda}^i$ can be expressed as the Hadamard product of matrices :

\begin{align} \bar{A}_{ij} &= \mathbf{ReDD}\circ\mathbf{a} = \text{ReDD[i][j]} a_{ij} \\ \bar{\Lambda}^{i} &= \mathbf{ReU}\circ\mathbf{\lambda} = \text{ReU[i]} \lambda^i, \end{align}

where no sums are implied by the repeated indices.

Further, since the scale factors are time independent,

\begin{align} \partial_t \bar{A}_{ij} &= \text{ReDD[i][j]}\ \partial_t a_{ij} \\ \partial_t \bar{\gamma}_{ij} &= \partial_t \left(\varepsilon_{ij} + \hat{\gamma}_{ij}\right)\\ &= \partial_t \varepsilon_{ij} \\ &= \text{scalefactor[i]}\ \text{scalefactor[j]}\ \partial_t h_{ij}. \end{align}

Thus instead of taking space or time derivatives of BSSN quantities

$$\left\{\bar{\gamma}_{i j},\bar{A}_{i j},\phi, K, \bar{\Lambda}^{i}, \alpha, \beta^i, B^i\right\},$$

across coordinate singularities, we instead factor out the singular scale factors according to this prescription so that space or time derivatives of BSSN quantities are written in terms of finite-difference derivatives of the rescaled variables

$$\left\{h_{i j},a_{i j},\text{cf}, K, \lambda^{i}, \alpha, \mathcal{V}^i, \mathcal{B}^i\right\},$$

and exact expressions for (spatial) derivatives of scale factors. Note that cf is the chosen conformal factor (supported choices for cf are discussed in Step 6.a).

As an example, let's evaluate $\bar{\Lambda}^{i}_{\, ,\, j}$ according to this prescription:

\begin{align} \bar{\Lambda}^{i}_{\, ,\, j} &= -\frac{\lambda^i}{(\text{ReU[i]})^2}\ \partial_j \left(\text{ReU[i]}\right) + \frac{\partial_j \lambda^i}{\text{ReU[i]}} \\ &= -\frac{\lambda^i}{(\text{ReU[i]})^2}\ \text{ReUdD[i][j]} + \frac{\partial_j \lambda^i}{\text{ReU[i]}}. \end{align}

Here, the derivative ReUdD[i][j] is computed symbolically and exactly using SymPy, and the derivative $\partial_j \lambda^i$ represents a derivative of a smooth quantity (so long as $\bar{\Lambda}^{i}$ is smooth in the Cartesian basis).

Step 3.a: BSSN_basic_tensors(): Define all basic conformal BSSN tensors $\left\{\bar{\gamma}_{i j},\bar{A}_{i j},\bar{\Lambda}^{i}, \beta^i, B^i\right\}$ in terms of BSSN gridfunctions [Back to top]¶

$$\label{bssn_basic_tensors}$$

The BSSN_vars__tensors() function defines the tensorial BSSN quantities $\left\{\bar{\gamma}_{i j},\bar{A}_{i j},\bar{\Lambda}^{i}, \beta^i, B^i\right\}$, in terms of the rescaled "base" tensorial quantities $\left\{h_{i j},a_{i j}, \lambda^{i}, \mathcal{V}^i, \mathcal{B}^i\right\},$ respectively:

\begin{align} \bar{\gamma}_{i j} &= \hat{\gamma}_{ij} + \varepsilon_{ij}, \text{ where } \varepsilon_{ij} = h_{ij} \circ \text{ReDD[i][j]} \\ \bar{A}_{i j} &= a_{ij} \circ \text{ReDD[i][j]} \\ \bar{\Lambda}^{i} &= \lambda^i \circ \text{ReU[i]} \\ \beta^{i} &= \mathcal{V}^i \circ \text{ReU[i]} \\ B^{i} &= \mathcal{B}^i \circ \text{ReU[i]} \end{align}

Rescaling vectors and tensors are built upon the scale factors for the chosen (in general, singular) coordinate system, which are defined in NRPy+'s reference_metric.py (Tutorial), and the rescaled variables are defined in the stub function BSSN/BSSN_rescaled_vars.py.

Here we implement BSSN_vars__tensors():

Step 4: gammabar__inverse_and_derivs(): $\bar{\gamma}^{ij}$, and spatial derivatives of $\bar{\gamma}_{ij}$ including $\bar{\Gamma}^{i}_{jk}$ [Back to top]¶

$$\label{bssn_barred_metric__inverse_and_derivs}$$

Step 4.a: Inverse conformal 3-metric: $\bar{\gamma^{ij}}$ [Back to top]¶

$$\label{bssn_barred_metric__inverse}$$

Since $\bar{\gamma}^{ij}$ is the inverse of $\bar{\gamma}_{ij}$, we apply a $3\times 3$ symmetric matrix inversion to compute $\bar{\gamma}^{ij}$.

Step 4.b: Derivatives of the conformal 3-metric $\bar{\gamma}_{ij,k}$ and $\bar{\gamma}_{ij,kl}$, and associated "barred" Christoffel symbols $\bar{\Gamma}^{i}_{jk}$ [Back to top]¶

$$\label{bssn_barred_metric__derivs}$$

In the BSSN-in-curvilinear coordinates formulation, all quantities must be defined in terms of rescaled quantities $h_{ij}$ and their derivatives (evaluated using finite differences), as well as reference-metric quantities and their derivatives (evaluated exactly using SymPy).

For example, $\bar{\gamma}_{ij,k}$ is given by: \begin{align} \bar{\gamma}_{ij,k} &= \partial_k \bar{\gamma}_{ij} \\ &= \partial_k \left(\hat{\gamma}_{ij} + \varepsilon_{ij}\right) \\ &= \partial_k \left(\hat{\gamma}_{ij} + h_{ij} \text{ReDD[i][j]}\right) \\ &= \hat{\gamma}_{ij,k} + h_{ij,k} \text{ReDD[i][j]} + h_{ij} \text{ReDDdD[i][j][k]}, \end{align} where ReDDdD[i][j][k] is computed within rfm.reference_metric().

By extension, the second derivative $\bar{\gamma}_{ij,kl}$ is given by \begin{align} \bar{\gamma}_{ij,kl} &= \partial_l \left(\hat{\gamma}_{ij,k} + h_{ij,k} \text{ReDD[i][j]} + h_{ij} \text{ReDDdD[i][j][k]}\right)\\ &= \hat{\gamma}_{ij,kl} + h_{ij,kl} \text{ReDD[i][j]} + h_{ij,k} \text{ReDDdD[i][j][l]} + h_{ij,l} \text{ReDDdD[i][j][k]} + h_{ij} \text{ReDDdDD[i][j][k][l]} \end{align}

Finally, we compute the Christoffel symbol associated with the barred 3-metric: $\bar{\Gamma}^{i}_{kl}$: $$ \bar{\Gamma}^{i}_{kl} = \frac{1}{2} \bar{\gamma}^{im} \left(\bar{\gamma}_{mk,l} + \bar{\gamma}_{ml,k} - \bar{\gamma}_{kl,m} \right) $$

Step 5: detgammabar_and_derivs(): $\det \bar{\gamma}_{ij}$ and its derivatives [Back to top]¶

$$\label{detgammabar_and_derivs}$$

As described just before Section III of Baumgarte et al (2012), we are free to choose $\det \bar{\gamma}_{ij}$, which should remain fixed in time.

As in Baumgarte et al (2012) generally we make the choice $\det \bar{\gamma}_{ij} = \det \hat{\gamma}_{ij}$, but this need not be the case; we could choose to set $\det \bar{\gamma}_{ij}$ to another expression.

In case we do not choose to set $\det \bar{\gamma}_{ij}/\det \hat{\gamma}_{ij}=1$, below we begin the implementation of a gridfunction, detgbarOverdetghat, which defines an alternative expression in its place.

$\det \bar{\gamma}_{ij}/\det \hat{\gamma}_{ij}$=detgbarOverdetghat$\ne 1$ is not yet implemented. However, we can define detgammabar and its derivatives in terms of a generic detgbarOverdetghat and $\det \hat{\gamma}_{ij}$ and their derivatives:

\begin{align} \text{detgammabar} &= \det \bar{\gamma}_{ij} = \text{detgbarOverdetghat} \cdot \left(\det \hat{\gamma}_{ij}\right) \\ \text{detgammabar}\_\text{dD[k]} &= \left(\det \bar{\gamma}_{ij}\right)_{,k} = \text{detgbarOverdetghat}\_\text{dD[k]} \det \hat{\gamma}_{ij} + \text{detgbarOverdetghat} \left(\det \hat{\gamma}_{ij}\right)_{,k} \\ \end{align}

https://en.wikipedia.org/wiki/Determinant#Properties_of_the_determinant

Step 6: AbarUU_AbarUD_trAbar_AbarDD_dD(): Quantities related to conformal traceless extrinsic curvature $\bar{A}_{ij}$: $\bar{A}^{ij}$, $\bar{A}^i_j$, and $\bar{A}^k_k$ [Back to top]¶

$$\label{abar_quantities}$$

$\bar{A}^{ij}$ is given by application of the raising operators (a.k.a., the inverse 3-metric) $\bar{\gamma}^{jk}$ on both of the covariant ("down") components: $$ \bar{A}^{ij} = \bar{\gamma}^{ik}\bar{\gamma}^{jl} \bar{A}_{kl}. $$

$\bar{A}^i_j$ is given by a single application of the raising operator (a.k.a., the inverse 3-metric) $\bar{\gamma}^{ik}$ on $\bar{A}_{kj}$: $$ \bar{A}^i_j = \bar{\gamma}^{ik}\bar{A}_{kj}. $$

The trace of $\bar{A}_{ij}$, $\bar{A}^k_k$, is given by a contraction with the barred 3-metric: $$ \text{Tr}(\bar{A}_{ij}) = \bar{A}^k_k = \bar{\gamma}^{kj}\bar{A}_{jk}. $$

Note that while $\bar{A}_{ij}$ is defined as the traceless conformal extrinsic curvature, it may acquire a nonzero trace (assuming the initial data impose tracelessness) due to numerical error. $\text{Tr}(\bar{A}_{ij})$ is included in the BSSN equations to drive $\text{Tr}(\bar{A}_{ij})$ to zero.

In terms of rescaled BSSN quantities, $\bar{A}_{ij}$ is given by $$ \bar{A}_{ij} = \text{ReDD[i][j]} a_{ij}, $$ so in terms of the same quantities, $\bar{A}_{ij,k}$ is given by $$ \bar{A}_{ij,k} = \text{ReDDdD[i][j][k]} a_{ij} + \text{ReDD[i][j]} a_{ij,k}. $$

Step 7: RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU(): The conformal ("barred") Ricci tensor $\bar{R}_{ij}$ and associated quantities [Back to top]¶

$$\label{rbar}$$

Let's compute perhaps the most complicated expression in the BSSN evolution equations, the conformal Ricci tensor:

\begin{align} \bar{R}_{i j} {} = {} & - \frac{1}{2} \bar{\gamma}^{k l} \hat{D}_{k} \hat{D}_{l} \bar{\gamma}_{i j} + \bar{\gamma}_{k(i} \hat{D}_{j)} \bar{\Lambda}^{k} + \Delta^{k} \Delta_{(i j) k} \nonumber \\ & + \bar{\gamma}^{k l} \left (2 \Delta_{k(i}^{m} \Delta_{j) m l} + \Delta_{i k}^{m} \Delta_{m j l} \right ) \; . \end{align}

Let's tackle the $\hat{D}_{k} \hat{D}_{l} \bar{\gamma}_{i j}$ term first:

Step 7.a: Conformal Ricci tensor, part 1: The $\hat{D}_{k} \hat{D}_{l} \bar{\gamma}_{i j}$ term [Back to top]¶

$$\label{rbar_part1}$$

First note that the covariant derivative of a metric with respect to itself is zero $$\hat{D}_{l} \hat{\gamma}_{ij} = 0,$$ so $$\hat{D}_{k} \hat{D}_{l} \bar{\gamma}_{i j} = \hat{D}_{k} \hat{D}_{l} \left(\hat{\gamma}_{i j} + \varepsilon_{ij}\right) = \hat{D}_{k} \hat{D}_{l} \varepsilon_{ij}.$$

Next, the covariant derivative of a tensor is given by (from the wikipedia article on covariant differentiation): \begin{align} {(\nabla_{e_c} T)^{a_1 \ldots a_r}}_{b_1 \ldots b_s} = {} &\frac{\partial}{\partial x^c}{T^{a_1 \ldots a_r}}_{b_1 \ldots b_s} \\ &+ \,{\Gamma ^{a_1}}_{dc} {T^{d a_2 \ldots a_r}}_{b_1 \ldots b_s} + \cdots + {\Gamma^{a_r}}_{dc} {T^{a_1 \ldots a_{r-1}d}}_{b_1 \ldots b_s} \\ &-\,{\Gamma^d}_{b_1 c} {T^{a_1 \ldots a_r}}_{d b_2 \ldots b_s} - \cdots - {\Gamma^d}_{b_s c} {T^{a_1 \ldots a_r}}_{b_1 \ldots b_{s-1} d}. \end{align}

Therefore, $$\hat{D}_{l} \bar{\gamma}_{i j} = \hat{D}_{l} \varepsilon_{i j} = \varepsilon_{i j,l} - \hat{\Gamma}^m_{i l} \varepsilon_{m j} -\hat{\Gamma}^m_{j l} \varepsilon_{i m}.$$

Since the covariant first derivative is a tensor, the covariant second derivative is given by (same as Eq. 27 in Baumgarte et al (2012))

\begin{align} \hat{D}_{k} \hat{D}_{l} \bar{\gamma}_{i j} &= \hat{D}_{k} \hat{D}_{l} \varepsilon_{i j} \\ &= \partial_k \hat{D}_{l} \varepsilon_{i j} - \hat{\Gamma}^m_{lk} \left(\hat{D}_{m} \varepsilon_{i j}\right) - \hat{\Gamma}^m_{ik} \left(\hat{D}_{l} \varepsilon_{m j}\right) - \hat{\Gamma}^m_{jk} \left(\hat{D}_{l} \varepsilon_{i m}\right), \end{align}

where the first term is the partial derivative of the expression already derived for $\hat{D}_{l} \varepsilon_{i j}$:

\begin{align} \partial_k \hat{D}_{l} \varepsilon_{i j} &= \partial_k \left(\varepsilon_{ij,l} - \hat{\Gamma}^m_{i l} \varepsilon_{m j} -\hat{\Gamma}^m_{j l} \varepsilon_{i m} \right) \\ &= \varepsilon_{ij,lk} - \hat{\Gamma}^m_{i l,k} \varepsilon_{m j} - \hat{\Gamma}^m_{i l} \varepsilon_{m j,k} - \hat{\Gamma}^m_{j l,k} \varepsilon_{i m} - \hat{\Gamma}^m_{j l} \varepsilon_{i m,k}. \end{align}

In terms of the evolved quantity $h_{ij}$, the derivatives of $\varepsilon_{ij}$ are given by: \begin{align} \varepsilon_{ij,k} &= \partial_k \left(h_{ij} \text{ReDD[i][j]}\right) \\ &= h_{ij,k} \text{ReDD[i][j]} + h_{ij} \text{ReDDdD[i][j][k]}, \end{align} and \begin{align} \varepsilon_{ij,kl} &= \partial_l \left(h_{ij,k} \text{ReDD[i][j]} + h_{ij} \text{ReDDdD[i][j][k]} \right)\\ &= h_{ij,kl} \text{ReDD[i][j]} + h_{ij,k} \text{ReDDdD[i][j][l]} + h_{ij,l} \text{ReDDdD[i][j][k]} + h_{ij} \text{ReDDdDD[i][j][k][l]}. \end{align}

We next compute three quantities derived above:

• gammabarDD_DhatD[i][j][l] = $\hat{D}_{l} \bar{\gamma}_{i j} = \hat{D}_{l} \varepsilon_{i j} = \varepsilon_{i j,l} - \hat{\Gamma}^m_{i l} \varepsilon_{m j} -\hat{\Gamma}^m_{j l} \varepsilon_{i m}$,
• gammabarDD_DhatD\_dD[i][j][l][k] = $\partial_k \hat{D}_{l} \bar{\gamma}_{i j} = \partial_k \hat{D}_{l} \varepsilon_{i j} = \varepsilon_{ij,lk} - \hat{\Gamma}^m_{i l,k} \varepsilon_{m j} - \hat{\Gamma}^m_{i l} \varepsilon_{m j,k} - \hat{\Gamma}^m_{j l,k} \varepsilon_{i m} - \hat{\Gamma}^m_{j l} \varepsilon_{i m,k}$, and
• gammabarDD_DhatDD[i][j][l][k] = $\hat{D}_{k} \hat{D}_{l} \bar{\gamma}_{i j} = \partial_k \hat{D}_{l} \varepsilon_{i j} - \hat{\Gamma}^m_{lk} \left(\hat{D}_{m} \varepsilon_{i j}\right) - \hat{\Gamma}^m_{ik} \left(\hat{D}_{l} \varepsilon_{m j}\right) - \hat{\Gamma}^m_{jk} \left(\hat{D}_{l} \varepsilon_{i m}\right)$.

Step 7.b: Conformal Ricci tensor, part 2: The $\bar{\gamma}_{k(i} \hat{D}_{j)} \bar{\Lambda}^{k}$ term [Back to top]¶

$$\label{rbar_part2}$$

By definition, the index symmetrization operation is given by: $$\bar{\gamma}_{k(i} \hat{D}_{j)} \bar{\Lambda}^{k} = \frac{1}{2} \left( \bar{\gamma}_{ki} \hat{D}_{j} \bar{\Lambda}^{k} + \bar{\gamma}_{kj} \hat{D}_{i} \bar{\Lambda}^{k} \right),$$

and $\bar{\gamma}_{ij}$ is trivially computed ($=\varepsilon_{ij} + \hat{\gamma}_{ij}$) so the only nontrival part to computing this term is in evaluating $\hat{D}_{j} \bar{\Lambda}^{k}$.

The covariant derivative is with respect to the hatted metric (i.e. the reference metric), so $$\hat{D}_{j} \bar{\Lambda}^{k} = \partial_j \bar{\Lambda}^{k} + \hat{\Gamma}^{k}_{mj} \bar{\Lambda}^m,$$ except we cannot take derivatives of $\bar{\Lambda}^{k}$ directly due to potential issues with coordinate singularities. Instead we write it in terms of the rescaled quantity $\lambda^k$ via $$\bar{\Lambda}^{k} = \lambda^k \text{ReU[k]}.$$

Then the expression for $\hat{D}_{j} \bar{\Lambda}^{k}$ becomes $$ \hat{D}_{j} \bar{\Lambda}^{k} = \lambda^{k}_{,j} \text{ReU[k]} + \lambda^{k} \text{ReUdD[k][j]} + \hat{\Gamma}^{k}_{mj} \lambda^{m} \text{ReU[m]}, $$ and the NRPy+ code for this expression is written

Step 7.c: Conformal Ricci tensor, part 3: The $\Delta^{k} \Delta_{(i j) k} + \bar{\gamma}^{k l} \left (2 \Delta_{k(i}^{m} \Delta_{j) m l} + \Delta_{i k}^{m} \Delta_{m j l} \right )$ terms [Back to top]¶

$$\label{rbar_part3}$$

Our goal here is to compute the quantities appearing as the final terms of the conformal Ricci tensor: $$ \Delta^{k} \Delta_{(i j) k} + \bar{\gamma}^{k l} \left (2 \Delta_{k(i}^{m} \Delta_{j) m l} + \Delta_{i k}^{m} \Delta_{m j l} \right). $$

• DGammaUDD[k][i][j]$= \Delta^k_{ij}$ is simply the difference in Christoffel symbols: $\Delta^{k}_{ij} = \bar{\Gamma}^i_{jk} - \hat{\Gamma}^i_{jk}$, and
• DGammaU[k]$= \Delta^k$ is the contraction: $\bar{\gamma}^{ij} \Delta^{k}_{ij}$

Adding these expressions to Ricci is straightforward, since $\bar{\Gamma}^i_{jk}$ and $\bar{\gamma}^{ij}$ were defined above in Step 4, and $\hat{\Gamma}^i_{jk}$ was computed within NRPy+'s reference_metric() function:

Next we define $\Delta_{ijk}=\bar{\gamma}_{im}\Delta^m_{jk}$:

Step 7.d: Summing the terms and defining $\bar{R}_{ij}$ [Back to top]¶

$$\label{summing_rbar_terms}$$

We have now constructed all of the terms going into $\bar{R}_{ij}$:

\begin{align} \bar{R}_{i j} {} = {} & - \frac{1}{2} \bar{\gamma}^{k l} \hat{D}_{k} \hat{D}_{l} \bar{\gamma}_{i j} + \bar{\gamma}_{k(i} \hat{D}_{j)} \bar{\Lambda}^{k} + \Delta^{k} \Delta_{(i j) k} \nonumber \\ & + \bar{\gamma}^{k l} \left (2 \Delta_{k(i}^{m} \Delta_{j) m l} + \Delta_{i k}^{m} \Delta_{m j l} \right ) \; . \end{align}

Step 8: betaU_derivs(): The unrescaled shift vector $\beta^i$ spatial derivatives: $\beta^i_{,j}$ & $\beta^i_{,jk}$, written in terms of the rescaled shift vector $\mathcal{V}^i$ [Back to top]¶

$$\label{beta_derivs}$$

This step, which documents the function betaUbar_and_derivs() inside the BSSN.BSSN_unrescaled_and_barred_vars module, defines three quantities:

• betaU_dD[i][j]$=\beta^i_{,j} = \left(\mathcal{V}^i \circ \text{ReU[i]}\right)_{,j} = \mathcal{V}^i_{,j} \circ \text{ReU[i]} + \mathcal{V}^i \circ \text{ReUdD[i][j]}$
• betaU_dupD[i][j]: the same as above, except using upwinded finite-difference derivatives to compute $\mathcal{V}^i_{,j}$ instead of centered finite-difference derivatives.
• betaU_dDD[i][j][k]$=\beta^i_{,jk} = \mathcal{V}^i_{,jk} \circ \text{ReU[i]} + \mathcal{V}^i_{,j} \circ \text{ReUdD[i][k]} + \mathcal{V}^i_{,k} \circ \text{ReUdD[i][j]}+\mathcal{V}^i \circ \text{ReUdDD[i][j][k]}$

Step 9: phi_and_derivs(): Standard BSSN conformal factor $\phi$, and its derivatives $\phi_{,i}$, $\phi_{,ij}$, $\bar{D}_j \phi$, and $\bar{D}_j\bar{D}_k \phi$, all written in terms of BSSN gridfunctions like $\text{cf}$ [Back to top]¶

$$\label{phi_and_derivs}$$

Step 9.a: $\phi$ in terms of the chosen (possibly non-standard) conformal factor variable $\text{cf}$ (e.g., $\text{cf}=\chi=e^{-4\phi}$) [Back to top]¶

$$\label{phi_ito_cf}$$

When solving the BSSN time evolution equations across the coordinate singularity (i.e., the "puncture") inside puncture black holes for example, the standard conformal factor $\phi$ becomes very sharp, whereas $\chi=e^{-4\phi}$ is far smoother (see, e.g., Campanelli, Lousto, Marronetti, and Zlochower (2006) for additional discussion). Thus if we choose to rewrite derivatives of $\phi$ in the BSSN equations in terms of finite-difference derivatives cf$=\chi$, numerical errors will be far smaller near the puncture.

The BSSN modules in NRPy+ support three options for the conformal factor variable cf:

1. cf$=\phi$,
2. cf$=\chi=e^{-4\phi}$, and
3. cf$=W = e^{-2\phi}$.

The BSSN equations are written in terms of $\phi$ (actually only $e^{-4\phi}$ appears) and derivatives of $\phi$, we now define $e^{-4\phi}$ and derivatives of $\phi$ in terms of the chosen cf.

First, we define the base variables needed within the BSSN equations:

Then we define $\phi_{,i}$, $\phi_{,ij}$, and $e^{-4\phi}$ for each of the choices of cf.

For cf$=\phi$, this is trivial:

For cf$=W=e^{-2\phi}$, we have

• $\phi_{,i} = -\text{cf}_{,i} / (2 \text{cf})$
• $\phi_{,ij} = (-\text{cf}_{,ij} + \text{cf}_{,i}\text{cf}_{,j}/\text{cf}) / (2 \text{cf})$
• $e^{-4\phi} = \text{cf}^2$

*Exercise to student: Prove the above relations*

For cf$=\chi=e^{-4\phi}$, we have

• $\phi_{,i} = -\text{cf}_{,i} / (4 \text{cf})$
• $\phi_{,ij} = (-\text{cf}_{,ij} + \text{cf}_{,i}\text{cf}_{,j}/\text{cf}) / (4 \text{cf})$
• $e^{-4\phi} = \text{cf}$

*Exercise to student: Prove the above relations*

Step 9.b: Covariant derivatives of $\phi$ [Back to top]¶

$$\label{phi_covariant_derivs}$$

Since $\phi$ is a scalar, $\bar{D}_i \phi = \partial_i \phi$.

Thus the second covariant derivative is given by \begin{align} \bar{D}_i \bar{D}_j \phi &= \phi_{;\bar{i}\bar{j}} = \bar{D}_i \phi_{,j}\\ &= \phi_{,ij} - \bar{\Gamma}^k_{ij} \phi_{,k}. \end{align}

Step 10: Code validation against BSSN.BSSN_quantities NRPy+ module [Back to top]¶

$$\label{code_validation}$$

As a code validation check, we verify agreement in the SymPy expressions for the RHSs of the BSSN equations between

1. this tutorial and
2. the NRPy+ BSSN.BSSN_quantities module.

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

Step 11: Output this notebook 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-BSSN_quantities.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)