Tutorial-BSSN_formulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# The BSSN Formulation of General Relativity in Generic Curvilinear Coordinates: An Overview\n",
    "\n",
    "## Author: Zach Etienne\n",
    "\n",
    "## This tutorial notebook demonstrates how Einstein's equations of general relativity in this formulation are constructed and output within NRPy+.\n",
    "\n",
    "\n",
    "### As Einstein's equations in this formalism take the form of highly nonlinear, coupled *wave equations*, the [tutorial notebook on the scalar wave equation in curvilinear coordinates](Tutorial-ScalarWaveCurvilinear.ipynb) is *required* reading before beginning this module. That module, as well as its own prerequisite [module on reference metrics within NRPy+](Tutorial-Reference_Metric.ipynb) provides the needed overview of how NRPy+ handles reference metrics.\n",
    "\n",
    "## Introduction:\n",
    "NRPy+'s original purpose was to be an easy-to-use code capable of generating Einstein's equations in a broad class of [singular](https://en.wikipedia.org/wiki/Coordinate_singularity), curvilinear coordinate systems, where the user need only input the scale factors of the underlying reference metric. Upon generating these equations, NRPy+ would then leverage SymPy's [common-expression-elimination (CSE)](https://en.wikipedia.org/wiki/Common_subexpression_elimination) and C code generation routines, coupled to its own [single-instruction, multiple-data (SIMD)](https://en.wikipedia.org/wiki/SIMD) functions, to generate highly-optimized C code.\n",
    "\n",
    "### Background Reading/Lectures:\n",
    "\n",
    "* Mathematical foundations of BSSN and 3+1 initial value problem decompositions of Einstein's equations:\n",
    "    * [Thomas Baumgarte's lectures on mathematical formulation of numerical relativity](https://www.youtube.com/watch?v=t3uo2R-yu4o&list=PLRVOWML3TL_djTd_nsTlq5aJjJET42Qke)\n",
    "    * [Yuichiro Sekiguchi's introduction to BSSN](http://www2.yukawa.kyoto-u.ac.jp/~yuichiro.sekiguchi/3+1.pdf) \n",
    "* Extensions to the standard BSSN approach used in NRPy+\n",
    "    * [Brown's covariant \"Lagrangian\" formalism of BSSN](https://arxiv.org/abs/0902.3652)\n",
    "    * [BSSN in spherical coordinates, using the reference-metric approach of Baumgarte, Montero, Cordero-Carrión, and Müller (2012)](https://arxiv.org/abs/1211.6632)\n",
    "    * [BSSN in generic curvilinear coordinates, using the extended reference-metric approach of Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658)\n",
    "\n",
    "### A Note on Notation:\n",
    "\n",
    "As is standard in NRPy+, \n",
    "\n",
    "* Greek indices refer to four-dimensional quantities where the zeroth component indicates temporal (time) component.\n",
    "* 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.\n",
    "\n",
    "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)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This module lays out the mathematical foundation for the BSSN formulation of Einstein's equations, as detailed in the references in the above Background Reading/Lectures section. It is meant to provide an overview of the basic equations and point of reference for **full tutorial notebooks** linked below:\n",
    "\n",
    "1. [Step 1](#brownslagrangebssn): [Brown](https://arxiv.org/abs/0902.3652)'s covariant formulation of the BSSN time-evolution equations (see next section for gauge conditions)\n",
    "    1. [Step 1.a](#fullequations): Numerical implementation of BSSN time-evolution equations\n",
    "        1. [Step 1.a.i](#liederivs) ([**BSSN quantities module [start here]**](Tutorial-BSSN_quantities.ipynb); [**BSSN time-evolution module**](Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb)): Expanding the Lie derivatives; the BSSN time-evolution equations in their final form\n",
    "1. [Step 2](#gaugeconditions) ([**full tutorial notebook**](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb)): Time-evolution equations for the BSSN |gauge quantities $\\alpha$ and $\\beta^i$\n",
    "1. [Step 3](#constraintequations) ([**full tutorial notebook**](Tutorial-BSSN_constraints.ipynb)): The BSSN constraint equations\n",
    "    1. [Step 3.a](#hamiltonianconstraint): The Hamiltonian constraint\n",
    "    1. [Step 3.b](#momentumconstraint): The momentum constraint\n",
    "1. [Step 4](#gammaconstraint) ([**full tutorial notebook**](Tutorial-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.ipynb)): The BSSN algebraic constraint $\\hat{\\gamma}=\\bar{\\gamma}$\n",
    "1. [Step 5](#latex_pdf_output) Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='brownslagrangebssn'></a>\n",
    "\n",
    "# Step 1: [Brown](https://arxiv.org/abs/0902.3652)'s covariant formulation of BSSN \\[Back to [top](#toc)\\]\n",
    "$$\\label{brownslagrangebssn}$$\n",
    "\n",
    "The covariant \"Lagrangian\" BSSN formulation of [Brown (2009)](https://arxiv.org/abs/0902.3652), which requires\n",
    "\n",
    "$$\n",
    "\\partial_t \\bar{\\gamma} = 0,\n",
    "$$\n",
    "\n",
    "results in the BSSN equations taking the following form (Eqs. 11 and 12 in [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658)):\n",
    "\n",
    "\\begin{align}\n",
    "  \\partial_{\\perp} \\bar{\\gamma}_{i j} {} = {} & \\frac{2}{3} \\bar{\\gamma}_{i j} \\left (\\alpha \\bar{A}_{k}^{k} - \\bar{D}_{k} \\beta^{k}\\right ) - 2 \\alpha \\bar{A}_{i j} \\; , \\\\\n",
    "  \\partial_{\\perp} \\bar{A}_{i j} {} = {} & -\\frac{2}{3} \\bar{A}_{i j} \\bar{D}_{k} \\beta^{k} - 2 \\alpha \\bar{A}_{i k} {\\bar{A}^{k}}_{j} + \\alpha \\bar{A}_{i j} K \\nonumber \\\\\n",
    "  & + e^{-4 \\phi} \\left \\{-2 \\alpha \\bar{D}_{i} \\bar{D}_{j} \\phi + 4 \\alpha \\bar{D}_{i} \\phi \\bar{D}_{j} \\phi \\right . \\nonumber \\\\\n",
    "    & \\left . + 4 \\bar{D}_{(i} \\alpha \\bar{D}_{j)} \\phi - \\bar{D}_{i} \\bar{D}_{j} \\alpha + \\alpha \\bar{R}_{i j} \\right \\}^{\\text{TF}} \\; , \\\\\n",
    "  \\partial_{\\perp} \\phi {} = {} & \\frac{1}{6} \\left (\\bar{D}_{k} \\beta^{k} - \\alpha K \\right ) \\; , \\\\\n",
    "  \\partial_{\\perp} K {} = {} & \\frac{1}{3} \\alpha K^{2} + \\alpha \\bar{A}_{i j} \\bar{A}^{i j} \\nonumber \\\\\n",
    "  & - e^{-4 \\phi} \\left (\\bar{D}_{i} \\bar{D}^{i} \\alpha + 2 \\bar{D}^{i} \\alpha \\bar{D}_{i} \\phi \\right ) \\; , \\\\\n",
    "  \\partial_{\\perp} \\bar{\\Lambda}^{i} {} = {} & \\bar{\\gamma}^{j k} \\hat{D}_{j} \\hat{D}_{k} \\beta^{i} + \\frac{2}{3} \\Delta^{i} \\bar{D}_{j} \\beta^{j} + \\frac{1}{3} \\bar{D}^{i} \\bar{D}_{j} \\beta^{j} \\nonumber \\\\\n",
    "  & - 2 \\bar{A}^{i j} \\left (\\partial_{j} \\alpha - 6 \\partial_{j} \\phi \\right ) + 2 \\bar{A}^{j k} \\Delta_{j k}^{i} \\nonumber \\\\\n",
    "  & -\\frac{4}{3} \\alpha \\bar{\\gamma}^{i j} \\partial_{j} K \\\\\n",
    "\\end{align}\n",
    "where \n",
    "* the $\\text{TF}$ superscript denotes the trace-free part.\n",
    "* $\\bar{\\gamma}_{ij} = \\varepsilon_{i j} + \\hat{\\gamma}_{ij}$, where $\\bar{\\gamma}_{ij} = e^{-4\\phi} \\gamma_{ij}$ is the conformal metric, $\\gamma_{ij}$ is the physical metric (see below), and $\\varepsilon_{i j}$ encodes information about the non-hatted metric.\n",
    "* $\\gamma_{ij}$, $\\beta^i$, and $\\alpha$ are the physical (as opposed to conformal) spatial 3-metric, shift vector, and lapse, respectively, which may be defined via the 3+1 decomposition line element (in [$G=c=1$ units](https://en.wikipedia.org/wiki/Planck_units)):\n",
    "$$ds^2 = -\\alpha^2 dt^2 + \\gamma_{ij}\\left(dx^i + \\beta^i dt\\right)\\left(dx^j + \\beta^j dt\\right).$$\n",
    "* $\\bar{R}_{ij}$ is the conformal Ricci tensor, computed via\n",
    "\\begin{align}\n",
    "  \\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 \\\\\n",
    "  & + \\bar{\\gamma}^{k l} \\left (2 \\Delta_{k(i}^{m} \\Delta_{j) m l} + \\Delta_{i k}^{m} \\Delta_{m j l} \\right ) \\; .\n",
    "\\end{align}\n",
    "* $\\partial_{\\perp} = \\partial_t - \\mathcal{L}_\\beta$; $\\mathcal{L}_\\beta$ is the [Lie derivative](https://en.wikipedia.org/wiki/Lie_derivative) along the shift vector $\\beta^i$.\n",
    "* $\\partial_0 = \\partial_t - \\beta^i \\partial_i$ is an advective time derivative.\n",
    "* $\\hat{D}_j$ is the [covariant derivative](https://en.wikipedia.org/wiki/Covariant_derivative) with respect to the reference metric $\\hat{\\gamma}_{ij}$.\n",
    "* $\\bar{D}_j$ is the [covariant derivative](https://en.wikipedia.org/wiki/Covariant_derivative) with respect to the barred spatial 3-metric $\\bar{\\gamma}_{ij}$\n",
    "* $\\Delta^i_{jk}$ is the tensor constructed from the difference of barred and hatted Christoffel symbols:\n",
    "$$\\Delta^i_{jk} = \\bar{\\Gamma}^i_{jk} - \\hat{\\Gamma}^i_{jk}$$\n",
    "    * The related quantity $\\Delta^i$ is defined $\\Delta^i \\equiv \\bar{\\gamma}^{jk} \\Delta^i_{jk}$.\n",
    "* $\\bar{A}_{ij}$ is the conformal, trace-free extrinsic curvature: \n",
    "$$\\bar{A}_{ij} = e^{-4\\phi} \\left(K_{ij} - \\frac{1}{3}\\gamma_{ij} K\\right),$$\n",
    "where $K$ is the trace of the extrinsic curvature $K_{ij}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='fullequations'></a>\n",
    "\n",
    "## Step 1.a: Numerical implementation of BSSN time-evolution equations \\[Back to [top](#toc)\\]\n",
    "$$\\label{fullequations}$$\n",
    "\n",
    "Regarding the numerical implementation of the above equations, first notice the left-hand sides of the equations include the time derivatives. Numerically, these equations are solved using as an [initial value problem](https://en.wikipedia.org/wiki/Initial_value_problem), where data are specified at a given time $t$, so that the solution at any later time can be obtained using the [Method of Lines (MoL)](https://en.wikipedia.org/wiki/Method_of_lines). MoL requires that the equations be written in the form:\n",
    "\n",
    "$$\\partial_t \\vec{U} = \\vec{f}\\left(\\vec{U},\\partial_i \\vec{U}, \\partial_i \\partial_j \\vec{U},...\\right),$$\n",
    "\n",
    "for the vector of \"evolved quantities\" $\\vec{U}$, where the right-hand side vector $\\vec{f}$ *does not* contain *explicit* time derivatives of $\\vec{U}$.\n",
    "\n",
    "Thus we must first rewrite the above equations so that *only* partial derivatives of time appear on the left-hand sides of the equations, meaning that the Lie derivative terms must be moved to the right-hand sides of the equations.\n",
    "\n",
    "<a id='liederivs'></a>\n",
    "\n",
    "### Step 1.a.i: Expanding the Lie derivatives; BSSN equations in their final form \\[Back to [top](#toc)\\]\n",
    "$$\\label{liederivs}$$\n",
    "\n",
    "In this Step, we provide explicit expressions for the [Lie derivatives](https://en.wikipedia.org/wiki/Lie_derivative) $\\mathcal{L}_\\beta$ appearing inside the $\\partial_\\perp = \\partial_t - \\mathcal{L}_\\beta$ operators for $\\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},W, K, \\bar{\\Lambda}^{i}\\right\\}$.\n",
    "\n",
    "In short, the Lie derivative of tensor weight $w$ is given by (from [the wikipedia article on Lie derivatives](https://en.wikipedia.org/wiki/Lie_derivative))\n",
    "\\begin{align}\n",
    "(\\mathcal {L}_X T) ^{a_1 \\ldots a_r}{}_{b_1 \\ldots b_s} &= X^c(\\partial_c T^{a_1 \\ldots a_r}{}_{b_1 \\ldots b_s}) \\\\\n",
    "&\\quad - (\\partial_c X ^{a_1}) T ^{c a_2 \\ldots a_r}{}_{b_1 \\ldots b_s} - \\ldots - (\\partial_c X^{a_r}) T ^{a_1 \\ldots a_{r-1}c}{}_{b_1 \\ldots b_s} \\\\\n",
    " &\\quad  +  (\\partial_{b_1} X^c) T ^{a_1 \\ldots a_r}{}_{c b_2 \\ldots b_s} + \\ldots + (\\partial_{b_s} X^c) T ^{a_1 \\ldots a_r}{}_{b_1 \\ldots b_{s-1} c} + w (\\partial_{c} X^c) T ^{a_1 \\ldots a_r}{}_{b_1 \\ldots b_{s}}\n",
    "\\end{align}\n",
    "\n",
    "Thus to evaluate the Lie derivative, one must first know the tensor density weight $w$ for each tensor. In this formulation of Einstein's equations, **all evolved quantities have density weight $w=0$**, so according to the definition of Lie derivative above,\n",
    "\\begin{align}\n",
    "\\mathcal{L}_\\beta \\bar{\\gamma}_{ij} &= \\beta^k \\partial_k \\bar{\\gamma}_{ij} + \\partial_i \\beta^k \\bar{\\gamma}_{kj} + \\partial_j \\beta^k \\bar{\\gamma}_{ik}, \\\\\n",
    "\\mathcal{L}_\\beta \\bar{A}_{ij} &= \\beta^k \\partial_k \\bar{A}_{ij} + \\partial_i \\beta^k \\bar{A}_{kj} + \\partial_j \\beta^k \\bar{A}_{ik}, \\\\\n",
    "\\mathcal{L}_\\beta \\phi &= \\beta^k \\partial_k \\phi, \\\\\n",
    "\\mathcal{L}_\\beta K &= \\beta^k \\partial_k K, \\\\\n",
    "\\mathcal{L}_\\beta \\bar{\\Lambda}^i &= \\beta^k \\partial_k \\bar{\\Lambda}^i - \\partial_k \\beta^i \\bar{\\Lambda}^k\n",
    "\\end{align}\n",
    "\n",
    "With these definitions, the BSSN equations for the un-rescaled evolved variables in the form $\\partial_t \\vec{U} = f\\left(\\vec{U},\\partial_i \\vec{U}, \\partial_i \\partial_j \\vec{U},...\\right)$ become\n",
    "\n",
    "\\begin{align}\n",
    "  \\partial_t \\bar{\\gamma}_{i j} {} = {} & \\left[\\beta^k \\partial_k \\bar{\\gamma}_{ij} + \\partial_i \\beta^k \\bar{\\gamma}_{kj} + \\partial_j \\beta^k \\bar{\\gamma}_{ik} \\right] + \\frac{2}{3} \\bar{\\gamma}_{i j} \\left (\\alpha \\bar{A}_{k}^{k} - \\bar{D}_{k} \\beta^{k}\\right ) - 2 \\alpha \\bar{A}_{i j} \\; , \\\\\n",
    "  \\partial_t \\bar{A}_{i j} {} = {} & \\left[\\beta^k \\partial_k \\bar{A}_{ij} + \\partial_i \\beta^k \\bar{A}_{kj} + \\partial_j \\beta^k \\bar{A}_{ik} \\right] - \\frac{2}{3} \\bar{A}_{i j} \\bar{D}_{k} \\beta^{k} - 2 \\alpha \\bar{A}_{i k} {\\bar{A}^{k}}_{j} + \\alpha \\bar{A}_{i j} K \\nonumber \\\\\n",
    "  & + e^{-4 \\phi} \\left \\{-2 \\alpha \\bar{D}_{i} \\bar{D}_{j} \\phi + 4 \\alpha \\bar{D}_{i} \\phi \\bar{D}_{j} \\phi  + 4 \\bar{D}_{(i} \\alpha \\bar{D}_{j)} \\phi - \\bar{D}_{i} \\bar{D}_{j} \\alpha + \\alpha \\bar{R}_{i j} \\right \\}^{\\text{TF}} \\; , \\\\\n",
    "  \\partial_t \\phi {} = {} & \\left[\\beta^k \\partial_k \\phi \\right] + \\frac{1}{6} \\left (\\bar{D}_{k} \\beta^{k} - \\alpha K \\right ) \\; , \\\\\n",
    "  \\partial_{t} K {} = {} & \\left[\\beta^k \\partial_k K \\right] + \\frac{1}{3} \\alpha K^{2} + \\alpha \\bar{A}_{i j} \\bar{A}^{i j} - e^{-4 \\phi} \\left (\\bar{D}_{i} \\bar{D}^{i} \\alpha + 2 \\bar{D}^{i} \\alpha \\bar{D}_{i} \\phi \\right ) \\; , \\\\\n",
    "  \\partial_t \\bar{\\Lambda}^{i} {} = {} & \\left[\\beta^k \\partial_k \\bar{\\Lambda}^i - \\partial_k \\beta^i \\bar{\\Lambda}^k \\right] + \\bar{\\gamma}^{j k} \\hat{D}_{j} \\hat{D}_{k} \\beta^{i} + \\frac{2}{3} \\Delta^{i} \\bar{D}_{j} \\beta^{j} + \\frac{1}{3} \\bar{D}^{i} \\bar{D}_{j} \\beta^{j} \\nonumber \\\\\n",
    "  & - 2 \\bar{A}^{i j} \\left (\\partial_{j} \\alpha - 6 \\partial_{j} \\phi \\right ) + 2 \\alpha \\bar{A}^{j k} \\Delta_{j k}^{i}  -\\frac{4}{3} \\alpha \\bar{\\gamma}^{i j} \\partial_{j} K\n",
    "\\end{align}\n",
    "\n",
    "where the terms moved from the right-hand sides to the left-hand sides are enclosed in square braces. \n",
    "\n",
    "Notice that the shift advection operator $\\beta^k \\partial_k \\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},\\phi, K, \\bar{\\Lambda}^{i}, \\alpha, \\beta^i, B^i\\right\\}$ appears on the right-hand side of *every* expression. As the shift determines how the spatial coordinates $x^i$ move on the next 3D slice of our 4D manifold, we find that representing $\\partial_k$ in these shift advection terms via an *upwinded* finite difference stencil results in far lower numerical errors. This trick is implemented below in all shift advection terms.\n",
    "\n",
    "As discussed in the [NRPy+ tutorial notebook on BSSN quantities](Tutorial-BSSN_quantities.ipynb), tensorial expressions can diverge at coordinate singularities, so each tensor in the set of BSSN variables\n",
    "\n",
    "$$\\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},\\phi, K, \\bar{\\Lambda}^{i}, \\alpha, \\beta^i, B^i\\right\\},$$\n",
    "\n",
    "is written in terms of the corresponding rescaled quantity in the set\n",
    "\n",
    "$$\\left\\{h_{i j},a_{i j},\\text{cf}, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\},$$ \n",
    "\n",
    "respectively, as defined in the [BSSN quantities tutorial](Tutorial-BSSN_quantities.ipynb). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='gaugeconditions'></a>\n",
    "\n",
    "# Step 2: Time-evolution equations for the BSSN gauge quantities $\\alpha$ and $\\beta^i$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{gaugeconditions}$$\n",
    "\n",
    "As described in the **Background Reading/Lectures** linked to above, the gauge quantities $\\alpha$ and $\\beta^i$ specify how coordinate time and spatial points adjust from one spatial hypersurface to the next, in our 3+1 decomposition of Einstein's equations.\n",
    "\n",
    "As choosing $\\alpha$ and $\\beta^i$ is equivalent to choosing coordinates for where we sample our solution to Einstein's equations, we are completely free to choose $\\alpha$ and $\\beta^i$ on any given spatial hypersuface. It has been found that fixing $\\alpha$ and $\\beta^i$ to constant values in the context of dynamical spacetimes results in instabilities, so we generally need to define expressions for $\\partial_t \\alpha$ and $\\partial_t \\beta^i$ and couple these equations to the rest of the BSSN time-evolution equations.\n",
    "\n",
    "Though we are free to choose the form of the right-hand sides of the gauge time evolution equations, very few have been found robust in the presence of (puncture) black holes.\n",
    "\n",
    "The most commonly adopted gauge conditions for BSSN (i.e., time-evolution equations for the BSSN gauge quantities $\\alpha$ and $\\beta^i$) are the \n",
    "\n",
    "* $1+\\log$ lapse condition:\n",
    "$$\n",
    "\\partial_0 \\alpha = -2 \\alpha K\n",
    "$$\n",
    "* Second-order Gamma-driving shift condition:\n",
    "\\begin{align}\n",
    "\\partial_0 \\beta^i &= B^{i} \\\\\n",
    "\\partial_0 B^i &= \\frac{3}{4} \\partial_{0} \\bar{\\Lambda}^{i} - \\eta B^{i},\n",
    "\\end{align}\n",
    "\n",
    "where $\\partial_0$ is the advection operator; i.e., $\\partial_0 A^i = \\partial_t A^i - \\beta^j \\partial_j A^i$. Note that $\\partial_{0} \\bar{\\Lambda}^{i}$ in the right-hand side of the $\\partial_{0} B^{i}$ equation is computed by adding $\\beta^j \\partial_j \\bar{\\Lambda}^i$ to the right-hand side expression given for $\\partial_t \\bar{\\Lambda}^i$, so no explicit time dependence occurs in the right-hand sides of the BSSN evolution equations and the Method of Lines can be applied directly.\n",
    "\n",
    "While it is incredibly robust in Cartesian coordinates, [Brown](https://arxiv.org/abs/0902.3652) pointed out that the above time-evolution equation for the shift is not covariant. In fact, we have found this non-covariant version to result in very poor results when solving Einstein's equations in spherical coordinates for a spinning black hole with spin axis pointed in the $\\hat{x}$ direction. Therefore we adopt Brown's covariant version as described in the [**full time-evolution equations for the BSSN gauge quantities $\\alpha$ and $\\beta^i$ tutorial notebook**](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='constraintequations'></a>\n",
    "\n",
    "# Step 3: The BSSN constraint equations \\[Back to [top](#toc)\\]\n",
    "$$\\label{constraintequations}$$\n",
    "\n",
    "In a way analogous to Maxwell's equations, the BSSN decomposition of Einstein's equations are written as a set of time-evolution equations and a set of constraint equations. In this step we present the BSSN constraints\n",
    "\n",
    "\\begin{align}\n",
    "\\mathcal{H} &= 0 \\\\\n",
    "\\mathcal{M^i} &= 0,\n",
    "\\end{align}\n",
    "\n",
    "where $\\mathcal{H}=0$ is the **Hamiltonian constraint**, and $\\mathcal{M^i} = 0$ is the **momentum constraint**. When constructing our spacetime from the initial data, one spatial hypersurface at a time, to confirm that at a given time, the Hamiltonian and momentum constraint violations converge to zero as expected with increased numerical resolution.\n",
    "\n",
    "<a id='hamiltonianconstraint'></a>\n",
    "\n",
    "## Step 3.a: The Hamiltonian constraint $\\mathcal{H}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{hamiltonianconstraint}$$\n",
    "\n",
    "The Hamiltonian constraint is written (Eq. 13 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf)):\n",
    "$$\n",
    "\\mathcal{H} = \\frac{2}{3} K^2 - \\bar{A}_{ij} \\bar{A}^{ij} + e^{-4\\phi} \\left(\\bar{R} - 8 \\bar{D}^i \\phi \\bar{D}_i \\phi - 8 \\bar{D}^2 \\phi\\right)\n",
    "$$\n",
    "\n",
    "<a id='momentumconstraint'></a>\n",
    "\n",
    "## Step 3.b: The momentum constraint $\\mathcal{M}^i$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{momentumconstraint}$$\n",
    "\n",
    "The momentum constraint is written (Eq. 47 of [Ruchlin, Etienne, & Baumgarte](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "\n",
    "$$    \\mathcal{M}^i = e^{-4\\phi} \\left(\n",
    "    \\frac{1}{\\sqrt{\\bar{\\gamma}}} \\hat{D}_j\\left(\\sqrt{\\bar{\\gamma}}\\bar{A}^{ij}\\right) +\n",
    "    6 \\bar{A}^{ij}\\partial_j \\phi -\n",
    "    \\frac{2}{3} \\bar{\\gamma}^{ij}\\partial_j K +\n",
    "    \\bar{A}^{jk} \\Delta\\Gamma^i_{jk} + \\bar{A}^{ik} \\Delta\\Gamma^j_{jk}\\right)\n",
    "$$\n",
    "\n",
    "Notice the momentum constraint as written in [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf) is missing a term, as described in [Ruchlin, Etienne, & Baumgarte](https://arxiv.org/pdf/1712.07658.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='gammaconstraint'></a>\n",
    "\n",
    "# Step 4: The BSSN algebraic constraint: $\\hat{\\gamma}=\\bar{\\gamma}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{gammaconstraint}$$\n",
    "\n",
    "[Brown](https://arxiv.org/abs/0902.3652)'s covariant Lagrangian formulation of BSSN, which we adopt, requires that $\\partial_t \\bar{\\gamma} = 0$, where $\\bar{\\gamma}=\\det \\bar{\\gamma}_{ij}$. We generally choose to set $\\bar{\\gamma}=\\hat{\\gamma}$ in our initial data.\n",
    "\n",
    "Numerical errors will cause $\\bar{\\gamma}$ to deviate from a constant in time. This actually disrupts the hyperbolicity of the PDEs (causing crashes), so to cure this, we adjust $\\bar{\\gamma}_{ij}$ at the end of each Runge-Kutta timestep, so that its determinant satisfies $\\bar{\\gamma}=\\hat{\\gamma}$ at all times. We adopt the following, rather standard prescription (Eq. 53 of [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658)):\n",
    "\n",
    "$$\n",
    "\\bar{\\gamma}_{ij} \\to \\left(\\frac{\\hat{\\gamma}}{\\bar{\\gamma}}\\right)^{1/3} \\bar{\\gamma}_{ij}.\n",
    "$$\n",
    "Notice the expression on the right is guaranteed to have determinant equal to $\\hat{\\gamma}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "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\n",
    "[Tutorial-BSSN_formulation.pdf](Tutorial-BSSN_formulation.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:21.309338Z",
     "iopub.status.busy": "2021-03-07T17:21:21.308582Z",
     "iopub.status.idle": "2021-03-07T17:21:24.791779Z",
     "shell.execute_reply": "2021-03-07T17:21:24.792580Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-BSSN_formulation.tex, and compiled LaTeX file to PDF file\n",
      "    Tutorial-BSSN_formulation.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-BSSN_formulation\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}