Tutorial-Derivation_of_GRMHD_Evolution_Equations.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",
    "# Derivation of the GRMHD Evolution Equations\n",
    "\n",
    "## Author: Samuel Cupp\n",
    "\n",
    "## This notebook provides an in-depth derivation of the evolution equations for GRMHD variables, based on the work of [Duez et al. (2005)](https://arxiv.org/abs/astro-ph/0503420). The notebook explains the two sets of variables often used in a general relativistic setting: the primitive variables ($\\mathbf{P} = \\left\\{\\rho_b,P,v^i,B^i\\right\\}$) and the conservative variables ($\\mathbf{C} = \\left\\{\\tilde{\\rho},\\tilde{\\tau},\\tilde{S}_i,\\tilde{B}^i\\right\\}$). Detailed steps are presented for deriving the evolution equations for these conservative variables, along with insights into the electromagnetic contribution to the stress-energy tensor, and the decision not to directly evolve the energy variable due to its potential inaccuracy.\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "When considering the evolution of magnetohydrodynamics in a general relativistic setting, there are generally two sets of variables, the primitive variables\n",
    "\n",
    "$$\n",
    "\\mathbf{P} = \\left\\{\\rho_b,P,v^i,B^i\\right\\}\n",
    "$$\n",
    "\n",
    "and the conservative variables\n",
    "\n",
    "$$\n",
    "\\mathbf{C} = \\left\\{\\tilde{\\rho},\\tilde{\\tau},\\tilde{S}_i,\\tilde{B}^i\\right\\}.\n",
    "$$\n",
    "\n",
    "The details of the individual variables in each set are discussed more in [Step 2](#suitable_equations). While the primitive variables are the quantities we associate with the usual description of physical properties, the conservative variables have far better numerical stability. As such, IllinoisGRMHD evolves $\\mathbf{C}$, and so the following evolution equations are in terms of these variables. The derivation in this document is based on the work by Duez _et al_ (citation below), and I have rederived their results. They label their primitive density as $\\rho_0$ (0 for rest mass density), while we label ours as $\\rho_b$ ($b$ for baryonic density). These are the same quantity despite the different labels.\n",
    "\n",
    "### Original Source\n",
    "* M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens. Relativistic magnetohydrodynamics in dynamical spacetimes: Numerical methods and tests. Phys. Rev. D 72, 024028 (2005). ([arxiv:astro-ph/0503420](https://arxiv.org/abs/astro-ph/0503420))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This module is organized as follows\n",
    "\n",
    "1. [Step 1](#EM_tensor): **Deriving $T_{EM}^{\\mu\\nu}$**\n",
    "    1. [Step 1.1](#F_scalar): Solving for $F_{\\alpha\\beta}F^{\\alpha\\beta}$\n",
    "    1. [Step 1.2](#F_tensor): Solving for $F^{\\mu\\lambda}F^{\\nu}_{\\hphantom{\\nu}\\lambda}$\n",
    "    1. [Step 1.3](#T_solve): Finishing Touches\n",
    "1. [Step 2](#suitable_equations): **Evolution Equations for the Conservative Variables**\n",
    "    1. [Step 2.1](#density_eqn): Evolution Equation for Fluid Density $\\tilde{\\rho}$\n",
    "    1. [Step 2.2](#momentum_eqn): Evolution Equation for Momentum Density Variable $\\tilde{S}_i$\n",
    "    1. [Step 2.3](#energy_eqn): Evolution Equation for Energy Density Variable $\\tilde{\\tau}$\n",
    "    1. [Step 2.4](#induction_eqn): Evolution Equation for the Magnetic Field $\\tilde{B}^i$\n",
    "    1. [Step 2.5](#summary): Summary of the Conservative Variable Evolution Equations\n",
    "1. [Appendix A](#levi_civita): **Levi-Civita Contractions**\n",
    "1. [Appendix B](#curv_tensor): **Extrinsic Curvature Identities**\n",
    "1. [Step 3](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF file**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='EM_tensor'></a>\n",
    "\n",
    "# Step 1: Deriving $T_{EM}^{\\mu\\nu}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{EM_tensor}$$\n",
    "\n",
    "In order to find evolution equations for the conserved quantities for the system, we start by finding a convenient form of the electromagnetic energy-momentum tensor $T_{EM}^{\\mu\\nu}$, which is defined as\n",
    "\n",
    "$$\n",
    "T_{EM}^{\\mu\\nu} = \\frac{1}{4\\pi}\\left( F^{\\mu\\lambda}F^{\\nu}_{\\hphantom{\\nu}\\lambda} - \\frac{1}{4}g^{\\mu\\nu}F_{\\alpha\\beta}F^{\\alpha\\beta} \\right)\n",
    "$$\n",
    "\n",
    "where $F^{\\mu\\nu}$ is the Faraday tensor. We want to examine this with respect to an arbitrary observer with normalized four-velocity $\\xi^\\mu$. If we consider the electromagnetic fields in the rest frame of this observer,\n",
    "\n",
    "$$\n",
    "\\xi_\\mu E^\\mu = 0 = \\xi_\\mu B^\\mu\n",
    "$$\n",
    "\n",
    "since $E$ and $B$ are purely spatial. The Faraday tensor can then be decomposed in terms of $E$ and $B$:\n",
    "\n",
    "$$\n",
    "F^{\\mu\\nu} = \\xi^\\mu E^\\nu - \\xi^\\nu E^\\mu + \\xi_\\gamma \\epsilon^{\\gamma\\mu\\nu\\delta} B_{\\delta}.\n",
    "$$\n",
    "\n",
    "##TODO: All the EM derivation stuff has a lot in common with Terrence's section. Further review and combination is probably necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='F_scalar'></a>\n",
    "\n",
    "## Step 1.1: Solving for $F_{\\alpha\\beta}F^{\\alpha\\beta}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{F_scalar}$$\n",
    "\n",
    "Now, we simply want to find an expression for of $T_{EM}^{\\mu\\nu}$ in terms of $\\xi$, $E$, and $B$. Let's first consider the relatively simpler term\n",
    "\n",
    "\\begin{align}\n",
    "F_{\\alpha\\beta}F^{\\alpha\\beta} &= \\left(\\xi_\\alpha E_\\beta - \\xi_\\beta E_\\alpha \n",
    "  + \\xi_\\gamma {\\epsilon^{\\gamma\\hphantom{\\alpha\\beta}\\delta}_{\\hphantom{\\gamma}\\alpha\\beta}} B_{\\delta}\\right)\n",
    "  \\left(\\xi^\\alpha E^\\beta - \\xi^\\beta E^\\alpha\n",
    "  + \\xi_\\sigma \\epsilon^{\\sigma\\alpha\\beta\\mu} B_{\\mu}\\right) \\\\\n",
    "%\n",
    "&= \\xi_\\alpha \\xi^\\alpha E_\\beta E^\\beta - \\xi_\\alpha E^\\alpha E_\\beta \\xi^\\beta - E_\\alpha \\xi^\\alpha \\xi_\\beta E^\\beta + \\xi_\\beta \\xi^\\beta E_\\alpha E^\\alpha \n",
    "  + \\xi_\\alpha E_\\beta \\xi_\\sigma \\epsilon^{\\sigma\\alpha\\beta\\mu} B_{\\mu} - \\xi_\\beta E_\\alpha \\xi_\\sigma \\epsilon^{\\sigma\\alpha\\beta\\mu} B_{\\mu} \\\\\n",
    "&\\ \\ \\ \\ \\ + \\xi^\\alpha E^\\beta \\xi_\\gamma {\\epsilon^{\\gamma\\hphantom{\\alpha\\beta}\\delta}_{\\hphantom{\\gamma}\\alpha\\beta}} B_{\\delta} - \\xi^\\beta E^\\alpha \\xi_\\gamma {\\epsilon^{\\gamma\\hphantom{\\alpha\\beta}\\delta}_{\\hphantom{\\gamma}\\alpha\\beta}} B_{\\delta} +  \\xi_\\sigma \\xi_\\gamma {\\epsilon^{\\gamma\\hphantom{\\alpha\\beta}\\delta}_{\\hphantom{\\gamma}\\alpha\\beta}} \\epsilon^{\\sigma\\alpha\\beta\\mu} B_{\\delta} B_{\\mu} \\\\\n",
    "%\n",
    "&= -2E_\\beta E^\\beta + 4\\xi_\\alpha \\xi_\\sigma \\epsilon^{\\sigma\\alpha\\beta\\mu} E_\\beta B_{\\mu}\n",
    "  + \\xi_\\sigma \\xi^\\gamma \\epsilon_{\\gamma\\alpha\\beta\\delta} \\epsilon^{\\sigma\\alpha\\beta\\mu} B^{\\delta} B_{\\mu}\n",
    "\\end{align}\n",
    "\n",
    "where we have used $\\xi^\\mu \\xi_\\mu=-1$ and the permutation rules of the Levi-Civita tensor. Further simplification requires abusing the properties of the Levi-Civita tensor, which is given in [Appendix A](#levi_civita):\n",
    "\n",
    "\\begin{align}\n",
    "F_{\\alpha\\beta}F^{\\alpha\\beta} &= -2E_\\beta E^\\beta + 4\\xi_\\alpha \\xi_\\sigma \\epsilon^{\\sigma\\alpha\\beta\\mu} E_\\beta B_{\\mu}\n",
    "  + 2\\xi_\\sigma \\xi^\\gamma (\\delta^\\sigma_\\delta \\delta^\\mu_\\gamma - \\delta^\\sigma_\\gamma \\delta^\\mu_\\delta) B^{\\delta} B_{\\mu} \\\\\n",
    "%\n",
    "&= -2E_\\beta E^\\beta + 4\\xi_\\alpha \\xi_\\sigma \\epsilon^{\\sigma\\alpha\\beta\\mu} E_\\beta B_{\\mu}\n",
    "  + 2\\xi_\\sigma \\xi^\\mu B^{\\sigma} B_{\\mu} - 2\\xi_\\sigma \\xi^\\sigma B^{\\mu} B_{\\mu} \\\\\n",
    "%\n",
    "&= -2E_\\beta E^\\beta + 4\\xi_\\alpha \\xi_\\sigma \\epsilon^{\\sigma\\alpha\\beta\\mu} E_\\beta B_{\\mu}\n",
    "  + 2B^{\\mu} B_{\\mu}\n",
    "\\end{align}\n",
    "\n",
    "Now, since we are in the rest frame of the observer, the spatial components of the observer's velocity are zero. Then, \n",
    "\\begin{align}\n",
    "\\xi_\\alpha \\xi_\\sigma \\epsilon^{\\sigma\\alpha\\beta\\mu} &= \\xi_0 \\xi_\\sigma \\epsilon^{\\sigma 0\\beta\\mu} \\\\\n",
    "&= \\xi_0^2 \\epsilon^{0 0\\beta\\mu} \\\\\n",
    "&= 0\n",
    "\\end{align}\n",
    "\n",
    "since any components of the Levi-Civita tensor with repeated indices are zero. Therefore,\n",
    "\n",
    "$$\n",
    "F_{\\alpha\\beta}F^{\\alpha\\beta} = -2E^\\beta E_\\beta + 2B^{\\mu} B_{\\mu}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='F_tensor'></a>\n",
    "\n",
    "## Step 1.2: Solving for $F^{\\mu\\lambda}F^{\\nu}_{\\ \\ \\lambda}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{F_tensor}$$\n",
    "\n",
    "In this section, we tackle the term $F^{\\mu\\lambda}F^{\\nu}_{\\hphantom{\\nu}\\lambda}$ in the energy-momentum tensor.\n",
    "\n",
    "\\begin{align}\n",
    "F^{\\mu\\lambda}F^{\\nu}_{\\hphantom{\\nu}\\lambda} &= \\left(\\xi^\\mu E^\\lambda - \\xi^\\lambda E^\\mu + \\xi_\\gamma \\epsilon^{\\gamma\\mu\\lambda\\delta} B_{\\delta}\\right)\n",
    "  \\left(\\xi^\\nu E_\\lambda - \\xi_\\lambda E^\\nu + \\xi_\\alpha \\epsilon^{\\alpha\\nu\\hphantom{\\lambda}\\beta}_{\\hphantom{\\alpha\\nu}\\lambda} B_{\\beta}\\right) \\\\\n",
    "%\n",
    "&= \\xi^\\mu E^\\lambda \\xi^\\nu E_\\lambda - \\xi^\\mu E^\\lambda \\xi_\\lambda E^\\nu -\\xi^\\lambda E^\\mu \\xi^\\nu E_\\lambda + \\xi^\\lambda E^\\mu \\xi_\\lambda E^\\nu \\\\\n",
    "&\\ \\ \\ \\ \\ + \\xi^\\mu E_\\lambda \\xi_\\alpha \\epsilon^{\\alpha\\nu\\lambda\\beta} B_{\\beta} - \\xi_\\lambda E^\\mu \\xi_\\alpha \\epsilon^{\\alpha\\nu\\lambda\\beta} B_{\\beta}\n",
    "+ \\xi_\\gamma \\epsilon^{\\gamma\\mu\\lambda\\delta} B_{\\delta} \\xi^\\nu E_\\lambda \\\\\n",
    "&\\ \\ \\ \\ \\ - \\xi_\\gamma \\epsilon^{\\gamma\\mu\\lambda\\delta} B_{\\delta} \\xi_\\lambda E^\\nu\n",
    "+ \\xi_\\gamma \\epsilon^{\\gamma\\mu\\lambda\\delta} B_{\\delta} \\xi_\\alpha \\epsilon^{\\alpha\\nu\\hphantom{\\lambda}\\beta}_{\\hphantom{\\alpha\\nu}\\lambda} B_{\\beta} \\\\\n",
    "%\n",
    "&= \\xi^\\mu \\xi^\\nu E^\\lambda E_\\lambda - E^\\mu E^\\nu +  \\xi_\\alpha E_\\lambda B_{\\beta} (\\xi^\\mu \\epsilon^{\\alpha\\nu\\lambda\\beta} + \\xi^\\nu \\epsilon^{\\alpha\\mu\\lambda\\beta}) \\\\\n",
    "&\\ \\ \\ \\ \\ - \\xi_\\lambda \\xi_\\alpha B_{\\beta} (\\epsilon^{\\alpha\\nu\\lambda\\beta} E^\\mu + \\epsilon^{\\alpha\\mu\\lambda\\beta} E^\\nu) + \\xi_\\alpha \\xi_\\gamma \\epsilon^{\\gamma\\mu\\lambda\\delta} \\epsilon^{\\alpha\\nu\\hphantom{\\lambda}\\beta}_{\\hphantom{\\alpha\\nu}\\lambda} B_{\\delta} B_{\\beta}\n",
    "\\end{align}\n",
    "\n",
    "where we have again used the relations $\\xi^\\mu \\xi_\\mu=-1$ and $\\xi_\\mu E^\\mu = 0 = \\xi_\\mu B^\\mu$. We've also taken the liberty to relabel dummy indices to simplify the resulting expression. To simplify the final term, we must (unfortunately) once again rely on the generalized Kronecker delta. This time we use the second (far more complicated) expression in [Appendix A](#levi_civita):\n",
    "\n",
    "\\begin{align}\n",
    "\\xi_\\alpha \\xi_\\gamma \\epsilon^{\\gamma\\mu\\lambda\\delta} \\epsilon^{\\alpha\\nu\\hphantom{\\lambda}\\beta}_{\\hphantom{\\alpha\\nu}\\lambda} B_{\\delta} B_{\\beta} &= -g^{\\nu\\phi} \\xi^\\alpha \\xi_\\gamma B_{\\delta} B^{\\beta} \\left( \\delta^\\gamma_\\alpha (\\delta^\\mu_\\phi \\delta^\\delta_\\beta - \\delta^\\mu_\\beta \\delta^\\delta_\\phi)\n",
    "  - \\delta^\\gamma_\\phi (\\delta^\\mu_\\alpha \\delta^\\delta_\\beta - \\delta^\\mu_\\beta \\delta^\\delta_\\alpha)\n",
    "  + \\delta^\\gamma_\\beta (\\delta^\\mu_\\alpha  \\delta^\\delta_\\phi - \\delta^\\mu_\\phi \\delta^\\delta_\\alpha) \\right) \\\\\n",
    "%\n",
    "&= -g^{\\nu\\phi} \\xi^\\alpha \\xi_\\alpha B_{\\delta} B^{\\beta} (\\delta^\\mu_\\phi \\delta^\\delta_\\beta - \\delta^\\mu_\\beta \\delta^\\delta_\\phi)\n",
    "  + g^{\\nu\\gamma} \\xi^\\alpha \\xi_\\gamma B_{\\delta} B^{\\beta} (\\delta^\\mu_\\alpha \\delta^\\delta_\\beta - \\delta^\\mu_\\beta \\delta^\\delta_\\alpha)\n",
    "  - g^{\\nu\\phi} \\xi^\\alpha \\xi_\\beta B_{\\delta} B^{\\beta} (\\delta^\\mu_\\alpha  \\delta^\\delta_\\phi - \\delta^\\mu_\\phi \\delta^\\delta_\\alpha) \\\\\n",
    "%\n",
    "&= -g^{\\mu\\nu} \\xi^\\alpha \\xi_\\alpha B_{\\beta} B^{\\beta}\n",
    "  + g^{\\nu\\delta} \\xi^\\alpha \\xi_\\alpha B_{\\delta} B^{\\mu}\n",
    "  + g^{\\nu\\gamma} \\xi^\\mu \\xi_\\gamma B_{\\beta} B^{\\beta}\n",
    "  - g^{\\nu\\gamma} \\xi^\\alpha \\xi_\\gamma B_{\\alpha} B^{\\mu}\n",
    "  - g^{\\nu\\delta} \\xi^\\mu \\xi_\\beta B_{\\delta} B^{\\beta}\n",
    "  + g^{\\mu\\nu} \\xi^\\alpha \\xi_\\beta B_{\\alpha} B^{\\beta} \\\\\n",
    "%\n",
    "&= g^{\\mu\\nu} B_{\\beta} B^{\\beta} - B^{\\mu} B^{\\nu} + \\xi^\\mu \\xi^\\nu B_{\\beta} B^{\\beta} \\\\\n",
    "\\end{align}\n",
    "\n",
    "which leaves us with\n",
    "\n",
    "\\begin{align}\n",
    "F^{\\mu\\lambda}F^{\\nu}_{\\hphantom{\\nu}\\lambda} &= \\xi^\\mu \\xi^\\nu E^\\lambda E_\\lambda - E^\\mu E^\\nu +  \\xi_\\alpha E_\\lambda B_{\\beta} (\\xi^\\mu \\epsilon^{\\alpha\\nu\\lambda\\beta} + \\xi^\\nu \\epsilon^{\\alpha\\mu\\lambda\\beta}) \\\\\n",
    "&\\ \\ \\ \\ \\ - \\xi_\\lambda \\xi_\\alpha B_{\\beta} (\\epsilon^{\\alpha\\nu\\lambda\\beta} E^\\mu + \\epsilon^{\\alpha\\mu\\lambda\\beta} E^\\nu) + g^{\\mu\\nu} B_{\\beta} B^{\\beta} - B^{\\mu} B^{\\nu} + \\xi^\\mu \\xi^\\nu B_{\\beta} B^{\\beta} \\\\\n",
    "%\n",
    "&= \\xi^\\mu \\xi^\\nu E^\\lambda E_\\lambda + \\left(g^{\\mu\\nu} + \\xi^\\mu \\xi^\\nu\\right) B_{\\beta} B^{\\beta} - E^\\mu E^\\nu - B^{\\mu} B^{\\nu} +  \\xi_\\alpha E_\\lambda B_{\\beta} \\left(\\xi^\\mu \\epsilon^{\\alpha\\nu\\lambda\\beta} + \\xi^\\nu \\epsilon^{\\alpha\\mu\\lambda\\beta}\\right)\n",
    "\\end{align}\n",
    "\n",
    "where we have again used the fact that $\\xi_\\alpha \\xi_\\sigma \\epsilon^{\\sigma\\alpha\\beta\\mu} = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='T_solve'></a>\n",
    "\n",
    "## Step 1.3: Finishing Touches \\[Back to [top](#toc)\\]\n",
    "$$\\label{T_solve}$$\n",
    "\n",
    "Now we just need to put all these pieces together. The final expression for the electromagnetic stress-energy tensor for an arbitrary observer is\n",
    "\n",
    "\\begin{align}\n",
    "T_{EM}^{\\mu\\nu}|_{\\xi\\rightarrow n} &= \\frac{1}{4\\pi}\\left( \\xi^\\mu \\xi^\\nu E^\\lambda E_\\lambda + \\left(g^{\\mu\\nu} + \\xi^\\mu \\xi^\\nu\\right) B_{\\beta} B^{\\beta} - E^\\mu E^\\nu - B^{\\mu} B^{\\nu} +  \\xi_\\alpha E_\\lambda B_{\\beta} \\left(\\xi^\\mu \\epsilon^{\\alpha\\nu\\lambda\\beta} + \\xi^\\nu \\epsilon^{\\alpha\\mu\\lambda\\beta}\\right) \\\\\n",
    "  - \\frac{1}{4}g^{\\mu\\nu}(-2E^\\beta E_\\beta + 2B^{\\mu} B_{\\mu}) \\right) \\\\\n",
    "&= \\frac{1}{8\\pi}\\left(g^{\\mu\\nu} + 2\\xi^\\mu \\xi^\\nu\\right)\\left(E^\\alpha E_\\alpha + B^{\\alpha} B_{\\alpha}\\right)\n",
    "  - \\frac{1}{4\\pi}\\left(E^\\mu E^\\nu + B^{\\mu} B^{\\nu}\\right)\n",
    "  + \\frac{1}{4\\pi} + \\xi_\\alpha E_\\lambda B_{\\beta} \\left(\\xi^\\mu \\epsilon^{\\alpha\\nu\\lambda\\beta} + \\xi^\\nu \\epsilon^{\\alpha\\mu\\lambda\\beta}\\right)\n",
    "\\end{align}\n",
    "\n",
    "Now, we can consider two special cases. First, let's consider the normal observer, in which case the normalized 4-velocity becomes the unit normal vector $\\xi\\rightarrow n$\n",
    "\n",
    "\\begin{align}\n",
    "n^\\nu &= \\left( \\frac{1}{\\alpha}, -\\frac{\\beta^i}{\\alpha} \\right) \\\\\n",
    "n_\\nu &= \\left( -\\alpha, 0 \\right).\n",
    "\\end{align}\n",
    "\n",
    "If we also define\n",
    "\n",
    "$$\n",
    "\\xi_\\alpha \\epsilon^{\\alpha\\nu\\lambda\\beta} \\equiv \\epsilon^{\\nu\\lambda\\beta}\n",
    "$$\n",
    "\n",
    "then the stress-energy tensor becomes\n",
    "\n",
    "\\begin{align}\n",
    "T_{EM}^{\\mu\\nu}|_{\\xi\\rightarrow n} &= \\frac{1}{8\\pi}\\left(g^{\\mu\\nu} + 2n^\\mu n^\\nu\\right)\\left(E^\\alpha E_\\alpha + B_{\\alpha} B^{\\alpha}\\right)\n",
    "  - \\frac{1}{4\\pi}\\left(E^\\mu E^\\nu + B^{\\mu} B^{\\nu}\\right)\n",
    "  + \\frac{1}{4\\pi}E_\\alpha B_{\\beta} n^{(\\mu} \\epsilon^{\\nu)\\alpha\\beta}\n",
    "\\end{align}\n",
    "\n",
    "We can also consider the observer co-moving with the fluid. We can call this co-moving velocity $u_\\mu$. In the case of ideal conditions (that is to say, a fluid with perfect conductivity), we have (from Ohm's law) that\n",
    "\n",
    "$$\n",
    "u_\\mu F^{\\mu\\nu}=0\n",
    "$$\n",
    "\n",
    "which implies that\n",
    "\n",
    "$$\n",
    "E^\\mu = u_\\nu F^{\\mu\\nu} = 0\n",
    "$$\n",
    "\n",
    "Since $E^\\mu = 0$, the stress-energy tensor simplifies significantly in this frame:\n",
    "\n",
    "$$\n",
    "T_{EM}^{\\mu\\nu}|_{\\xi\\rightarrow u} = \\frac{1}{8\\pi}\\left(g^{\\mu\\nu} + 2u^\\mu u^\\nu\\right)B_{\\beta} B^{\\beta}\n",
    "  - \\frac{1}{4\\pi}B^{\\mu} B^{\\nu}\n",
    "$$\n",
    "\n",
    "However, variable $b$ is often used instead of $B^\\mu_{(u)}$, given by\n",
    "\n",
    "$$\n",
    "b^\\mu = \\frac{B^\\mu_{(u)}}{\\sqrt{4\\pi}}\n",
    "$$\n",
    "\n",
    "In terms of $b$, the electromagnetic energy-momentum tensor is given by\n",
    "\n",
    "$$\n",
    "T_{EM}^{\\mu\\nu}|_{\\xi\\rightarrow u} = b^2\\left(\\frac{1}{2}g^{\\mu\\nu} + u^\\mu u^\\nu\\right)\n",
    "  - b^{\\mu} b^{\\nu}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='suitable_equations'></a>\n",
    "\n",
    "# Step 2: Evolution Equations for the Conservative Variables \\[Back to [top](#toc)\\]\n",
    "$$\\label{suitable_equations}$$\n",
    "\n",
    "Before discussing the evolution variables, recall the 3+1 decomposition which is used for GRMHD. The metric is given by\n",
    "$$/\n",
    "ds^2 = -\\alpha^2 dt^2 + \\gamma_{i j}\\left(dx^i + \\beta^j dt\\right)\\left(dx^j + \\beta^i dt\\right)\n",
    "$$\n",
    "where $\\alpha$, $\\beta^i$, and $\\gamma_{i j}$ are the lapse, shift, and spatial metric, respectively. These quantities are used in the definitions of the conservative variables which we will soon introduce.\n",
    "\n",
    "In GRMHD simulations, there are two sets of variables which can be used for evolution. The first set, called primitive variables, are\n",
    "\n",
    "$$\n",
    "\\mathbf{P} = \\left\\{\\rho_b,P,v^i,B^i\\right\\}\n",
    "$$\n",
    "\n",
    "where $\\rho_b$ is the fluid rest-mass density, $P$ is the pressure, $v^i$ is the rescaled fluid three-velocity $u^i/u^0$, and $B^i$ is the magnetic field. An important sidenote here is that different GRMHD codes define different $v^i$. As an example, the widely used GRHydro thorn in the Einstein Toolkit uses the Valencia formalism\n",
    "\n",
    "$$\n",
    "v^i = \\frac{u^i}{W} + \\frac{\\beta^i}{\\alpha}\n",
    "$$\n",
    "\n",
    "where $W=\\sqrt{1-v^i v_i}$ is the Lorentz factor. The quantity $v^i$ used by IllinoisGRMHD is **not** the Valencia three-velocity, so those interested in seeing more about the Valencia formulation of the evolution equations can look to the [GRHydro paper](https://arxiv.org/abs/1304.5544) and its references.\n",
    "\n",
    "The second set, the conservative variables, are\n",
    "\n",
    "$$\n",
    "\\mathbf{C} = \\left\\{\\tilde{\\rho},\\tilde{\\tau},\\tilde{S}_i,\\tilde{B}^i\\right\\}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{align}\n",
    "\\tilde{\\rho} &= \\sqrt{-g}\\rho_b u^0 = \\alpha\\sqrt{\\gamma}\\rho_b u^0 \\\\\n",
    "\\tilde{\\tau} &= \\sqrt{\\gamma} n_\\mu n_\\nu T^{\\mu\\nu} - \\tilde{\\rho} = \\alpha^2 \\sqrt{\\gamma} T^{00} - \\tilde{\\rho} \\\\\n",
    "\\tilde{S}_i &= \\sqrt{\\gamma}S_i = \\alpha \\sqrt{\\gamma}T^0_i \\\\\n",
    "\\tilde{B}^i &= \\sqrt{\\gamma}B^i\n",
    "\\end{align}\n",
    "\n",
    "One can design numerical codes to evolve either set of variables, and we choose to evolve $\\mathbf{C}$. Conservative variables are used to evolve these systems because the evolution of primitive variables leads to highly unstable simulations. This change of variables for the sake of stability has some similarity to the change from ADM variables to BSSN variables when evolving space-time, where one set of variables has substantially better numerical stability.\n",
    "\n",
    "Another question that may arise is why we do not simply evolve the energy variable $\\alpha^2 \\sqrt{\\gamma} T^{00}$ and instead subtract $\\tilde{\\rho}$. It turns out that evolving $T^{00}$ directly is inaccurate because the magnetic and internal energy density can be orders of magnitude smaller than the rest mass density (this is mentioned in e.g. the [`HARM` software paper](https://arxiv.org/pdf/astro-ph/0301509.pdf)). Because of this, we subtract the baryonic density evolution equation from the energy evolution equation to find a more accurate and stable evolution equation.\n",
    "\n",
    "All of our previous derivations, of course, only considered the electromagnetic contribution. We do not derive the hydrodynamic part of the stress-energy tensor and instead merely state that in the ideal MHD limit the full stress-energy tensor takes the form\n",
    "\n",
    "$$\n",
    "T^{\\mu\\nu} = \\rho_b h u^\\mu u^\\nu + P g^{\\mu\\nu} + T_{EM}^{\\mu\\nu}\n",
    "$$\n",
    "\n",
    "Our next step is to derive the actual evolution equations for our designated evolution variables $\\mathbf{C}$. To this end, we will use the identity\n",
    "\n",
    "$$\n",
    "\\Gamma^\\nu_{\\gamma\\nu} = \\frac{1}{\\sqrt{-g}} \\partial_\\gamma\\sqrt{-g}\n",
    "$$\n",
    "\n",
    "to simplify the equations. We will be frequently using this to condense two terms into one in the following way:\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_\\nu T_\\mu^\\nu + \\Gamma^\\nu_{\\gamma\\nu} T^\\gamma_\\mu &= \\partial_\\nu T_\\mu^\\nu + \\frac{1}{\\sqrt{-g}} T^\\gamma_\\mu \\partial_\\gamma\\sqrt{-g} \\\\\n",
    "%\n",
    "&= \\partial_\\nu T_\\mu^\\nu + \\frac{1}{\\sqrt{-g}} T^\\gamma_\\mu \\partial_\\gamma\\sqrt{-g} \\\\\n",
    "%\n",
    "&= \\frac{1}{\\sqrt{-g}} \\left(\\sqrt{-g} \\partial_\\nu T_\\mu^\\nu + T^\\gamma_\\mu \\partial_\\gamma\\sqrt{-g}\\right) \\\\\n",
    "%\n",
    "&= \\frac{1}{\\sqrt{-g}} \\partial_\\gamma \\left(\\sqrt{-g}T^\\gamma_\\mu\\right)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='density_eqn'></a>\n",
    "\n",
    "## Step 2.1: Evolution Equation for the Fluid Density $\\tilde{\\rho}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{density_eqn}$$\n",
    "\n",
    "The first evolution equation comes from the baryon number conservation equation\n",
    "\n",
    "\\begin{align}\n",
    "0 &= \\nabla_\\nu\\left( \\rho_b u^\\nu \\right) \\\\\n",
    "&= \\partial_\\nu\\left( \\rho_b u^\\nu \\right) + \\Gamma^\\nu_{\\gamma\\nu} \\rho_b u^\\gamma \\\\\n",
    "&= \\partial_\\nu\\left( \\rho_b u^\\nu \\right) + \\rho_b u^\\gamma \\frac{1}{\\sqrt{-g}} \\partial_\\gamma\\sqrt{-g} \\\\\n",
    "&= \\frac{1}{\\sqrt{-g}}\\partial_\\nu \\left(\\sqrt{-g}\\rho_b u^\\nu \\right) \\\\\n",
    "&= \\partial_t \\left(\\alpha\\sqrt{\\gamma}\\rho_b u^0 \\right) + \\partial_i \\left(\\alpha\\sqrt{\\gamma}\\rho_b u^i \\right) \\\\\n",
    "&= \\partial_t \\tilde{\\rho} + \\partial_i \\left(\\frac{\\tilde{\\rho}}{u^0} u^i \\right) \\\\\n",
    "&= \\partial_t \\tilde{\\rho} + \\partial_i \\left(\\tilde{\\rho} v^i \\right)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='momentum_eqn'></a>\n",
    "\n",
    "## Step 2.2: Evolution Equation for Momentum Density Variable $\\tilde{S}_i$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{momentum_eqn}$$\n",
    "\n",
    "The next two sections derive their equations from the energy-momentum conservation equation\n",
    "\n",
    "\\begin{align}\n",
    "\\nabla_\\nu T_\\mu^\\nu &= 0 \\\\\n",
    "&= \\partial_\\nu T_\\mu^\\nu + \\Gamma^\\nu_{\\gamma\\nu} T^\\gamma_\\mu - \\Gamma^\\gamma_{\\mu\\nu} T^\\nu_\\gamma\n",
    "\\end{align}\n",
    "\n",
    "Rearranging and using the identity for $\\Gamma^\\nu_{\\gamma\\nu}$, we get\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\sqrt{-g}}\\partial_\\nu \\left(\\sqrt{-g} T_\\mu^\\nu \\right) = \\Gamma^\\gamma_{\\mu\\nu} T^\\nu_\\gamma\n",
    "$$\n",
    "\n",
    "For the momentum density variable $\\tilde{S}_i$, we consider the indices $\\mu\\in\\{1,2,3\\}$. Replacing $\\mu$ with $i$ to find the evolution equation,\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{1}{\\sqrt{-g}}\\partial_\\nu \\left(\\sqrt{-g} T_i^\\nu \\right) &= \\Gamma^\\gamma_{i\\nu} T^\\nu_\\gamma \\\\\n",
    "%\n",
    "\\frac{1}{\\sqrt{-g}}\\partial_t \\left(\\sqrt{-g} T_i^0 \\right) + \\frac{1}{\\sqrt{-g}}\\partial_j \\left(\\sqrt{-g} T_i^j \\right) &= \\frac{1}{2}T^\\nu_\\gamma g^{\\gamma\\beta}\\left( g_{\\beta\\nu,i} + g_{i\\beta,\\nu} - g_{i\\nu,\\beta} \\right) \\\\\n",
    "%\n",
    "\\partial_t \\tilde{S}_i + \\partial_j \\left(\\alpha\\sqrt{\\gamma} T_i^j \\right) &= \\frac{\\alpha\\sqrt{\\gamma}}{2} T^{\\beta\\nu}g_{\\beta\\nu,i} + \\frac{\\alpha\\sqrt{\\gamma}}{2} T^{\\beta\\nu}\\left( g_{i\\beta,\\nu} - g_{i\\nu,\\beta} \\right) \\\\\n",
    "%\n",
    "\\partial_t \\tilde{S}_i + \\partial_j \\left(\\alpha\\sqrt{\\gamma} T_i^j \\right) &= \\frac{\\alpha\\sqrt{\\gamma}}{2} T^{\\beta\\nu}g_{\\beta\\nu,i}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='energy_eqn'></a>\n",
    "\n",
    "## Step 2.3: Evolution Equation for Energy Density Variable $\\tilde{\\tau}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{energy_eqn}$$\n",
    "\n",
    "The energy evolution equation comes from the same conservation equation as the previous section, but with $\\mu=0$. However, we choose to start with $\\mu$ raised instead of lowered. Then,\n",
    "\n",
    "\\begin{align}\n",
    "0 &= \\nabla_\\nu T^{\\mu\\nu} \\\\\n",
    "&= \\partial_\\nu T^{\\mu\\nu} + \\Gamma^\\mu_{\\sigma\\nu} T^{\\sigma\\nu} + \\Gamma^\\nu_{\\sigma\\nu} T^{\\mu\\sigma} \\\\\n",
    "&= \\frac{1}{\\sqrt{-g}}\\partial_\\nu \\left(\\sqrt{-g} T^{\\mu\\nu} \\right) + \\Gamma^\\mu_{\\sigma\\nu} T^{\\sigma\\nu}\n",
    "\\end{align}\n",
    "\n",
    "Setting $\\mu=0$,\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{1}{\\sqrt{-g}}\\partial_\\nu \\left(\\sqrt{-g} T^{0\\nu} \\right) &= -\\Gamma^0_{\\sigma\\nu} T^{\\sigma\\nu} \\\\\n",
    "\\partial_t \\left(\\alpha\\sqrt{\\gamma} T^{00} \\right) + \\partial_i \\left(\\alpha\\sqrt{\\gamma} T^{0i} \\right) &= -\\alpha\\sqrt{\\gamma}\\Gamma^0_{\\sigma\\nu} T^{\\sigma\\nu} \\\\\n",
    "&= -\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{\\sigma\\nu} g^{0\\beta}\\left( g_{\\beta\\sigma,\\nu} + g_{\\beta\\nu,\\sigma} - g_{\\sigma\\nu,\\beta} \\right) \\\\\n",
    "&= -\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{\\sigma\\nu} g^{0\\beta}\\left(2 g_{\\beta\\sigma,\\nu} - g_{\\sigma\\nu,\\beta} \\right)\n",
    "\\end{align}\n",
    "\n",
    "where the final step follows from the fact that $T^{\\sigma\\nu}$ is symmetric. Now, we will consider the right-hand side for various components of $T^{\\sigma\\nu}$. For these derivations, we will use derivations relating to the extrinsic curvature in [Appendix B](#curv_tensor). For $T^{00}$,\n",
    "\n",
    "\\begin{align}\n",
    "-\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{00} g^{0\\beta}\\left(2 g_{\\beta 0,0} - g_{00,\\beta} \\right) &= -\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{00} \\left[2 \\left(g^{00} g_{00,0} + g^{0i}g_{0i,0}\\right) - g^{00}g_{00,0} - g^{0i}g_{00,i} \\right] \\\\\n",
    "%\n",
    "&= -\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{00} \\left( g^{00} g_{00,0} + 2g^{0i}g_{0i,0}  - g^{0i}g_{00,i} \\right) \\\\\n",
    "%\n",
    "&= -\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{00} \\left[ -\\alpha^{-2} \\partial_t \\left(\\beta^2 - \\alpha^2\\right) + \\frac{2\\beta^i}{\\alpha^2} \\partial_t \\beta_i - \\frac{\\beta^i}{\\alpha^2} \\partial_i \\left(\\beta^2 - \\alpha^2\\right) \\right] \\\\\n",
    "%\n",
    "&= -\\frac{\\sqrt{\\gamma}}{2\\alpha} T^{00} \\left[ -\\partial_t \\left(\\beta^2 - \\alpha^2\\right) + 2\\beta^i \\partial_t \\beta_i - \\beta^i \\partial_i \\left(\\beta^2 - \\alpha^2\\right) \\right] \\\\\n",
    "%\n",
    "&= -\\frac{\\sqrt{\\gamma}}{\\alpha} T^{00} \\left( \\alpha\\partial_t \\alpha + \\beta^i \\alpha\\partial_i \\alpha - \\beta^i \\beta^j \\partial_i \\beta_j \\right) \\\\\n",
    "%\n",
    "&= \\sqrt{\\gamma} T^{00} \\left( \\beta^i \\beta^j K_{ij} - \\partial_t \\alpha - \\beta^i \\partial_i \\alpha \\right)\n",
    "\\end{align}\n",
    "\n",
    "where we have used the identity in [Appendix B](#curv_tensor). Next, we look at the mixed term $T^{0i} + T^{i0}$:\n",
    "\n",
    "\\begin{align}\n",
    "-\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{0i} g^{0\\beta}\\left(2 g_{\\beta 0,i} - g_{0i,\\beta} + 2 g_{\\beta i,0} - g_{i 0,\\beta} \\right) &= -\\alpha\\sqrt{\\gamma} T^{0i} \\left[ g^{00}\\left( g_{00,i} + g_{0i,0} \\right) + g^{0j}\\left( g_{0j,i} + g_{ij,0} \\right) - \\left(g^{00} g_{0i,0} + g^{0j} g_{0i,j}\\right) \\right] \\\\\n",
    "%\n",
    "&= -\\alpha\\sqrt{\\gamma} T^{0i} \\left[ g^{00} g_{00,i} + g^{0j}\\left( g_{0j,i} + g_{ij,0} - g_{0i,j} \\right) \\right] \\\\\n",
    "%\n",
    "&= -\\alpha\\sqrt{\\gamma} T^{0i} \\left[ -\\alpha^{-2} \\partial_i \\left( \\beta^2 - \\alpha^2 \\right) + \\frac{\\beta^j}{\\alpha^2} \\left( \\partial_i \\beta_j + \\partial_t \\gamma_{ij} - \\partial_j \\beta_i \\right) \\right] \\\\\n",
    "%\n",
    "&= -\\frac{\\sqrt{\\gamma}}{\\alpha} T^{0i} \\left[ 2\\alpha\\partial_i \\alpha + \\beta^j \\left(  + \\partial_t \\gamma_{ij} - \\partial_i \\beta_j - \\partial_j \\beta_i \\right) \\right] \\\\\n",
    "%\n",
    "&= -\\frac{\\sqrt{\\gamma}}{\\alpha} T^{0i} \\left[ 2\\alpha\\partial_i \\alpha + \\beta^j \\left( - 2\\alpha K_{ij} - 2\\Gamma^k_{ij}\\beta_k \\right) \\right] \\\\\n",
    "%\n",
    "&= 2\\sqrt{\\gamma} T^{0i} \\left( \\beta^j K_{ij} - \\partial_i \\alpha \\right) + \\frac{2\\sqrt{\\gamma}}{\\alpha} T^{0i} \\beta^j \\Gamma^k_{ij}\\beta_k \\\\\n",
    "%\n",
    "&= 2\\sqrt{\\gamma} T^{0i} \\left( \\beta^j K_{ij} - \\partial_i \\alpha \\right)\n",
    "\\end{align}\n",
    "\n",
    "Finally, for $T^{ij}$ we have\n",
    "\n",
    "\\begin{align}\n",
    "-\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{ij} g^{0\\beta}\\left(2 g_{\\beta i,j} - g_{ij,\\beta} \\right) &= -\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{ij} \\left[ 2\\left(g^{00} g_{0i,j} + g^{0k} g_{ki,j}\\right) - g^{0\\beta} g_{ij,\\beta} \\right] \\\\\n",
    "&= -\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{ij} \\left[ 2\\left(g^{00} \\partial_j \\beta_i + g^{0k} \\partial_j \\gamma_{ki}\\right) - g^{0\\beta} \\partial_\\beta \\gamma_{ij} \\right] \\\\\n",
    "&= -\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{ij} \\left[ g^{00}\\left( 2\\partial_j \\beta_i - \\partial_t \\gamma_{ij} \\right) + g^{0k}\\left( 2\\partial_j \\gamma_{ki} - \\partial_k \\gamma_{ij} \\right) \\right] \\\\\n",
    "&= -\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{ij} \\left[ -\\alpha^{-2}\\left( 2\\partial_j \\beta_i - \\partial_t \\gamma_{ij} \\right) + \\frac{\\beta^k}{\\alpha^2}\\left( 2\\partial_j \\gamma_{ki} - \\partial_k \\gamma_{ij} \\right) \\right] \\\\\n",
    "\\end{align}\n",
    "\n",
    "We can see that, thanks to the symmetry of $T^{ij}$,\n",
    "\n",
    "\\begin{align}\n",
    "T^{ij}\\beta^k\\left(2\\partial_j \\gamma_{ki} - \\partial_k \\gamma_{ij}\\right) &= T^{ij}\\beta_n g^{kn}\\left(\\partial_j \\gamma_{ki} + \\partial_i \\gamma_{kj} - \\partial_k \\gamma_{ij}\\right) \\\\\n",
    "&= 2T^{ij}\\Gamma^k_{ij}\\beta_k\n",
    "\\end{align}\n",
    "\n",
    "Using this and using the definition of $K_{ij}$ to replace $\\partial_t \\gamma_{ij}$,\n",
    "\n",
    "\\begin{align}\n",
    "-\\frac{\\alpha\\sqrt{\\gamma}}{2} T^{ij} g^{0\\beta}\\left(2 g_{\\beta i,j} - g_{ij,\\beta} \\right) &= \\frac{\\sqrt{\\gamma}}{2\\alpha} T^{ij} \\left( 2\\partial_j \\beta_i + 2\\alpha K_{ij} - \\partial_i\\beta_j - \\partial_j\\beta_i + 2\\Gamma^k_{ij}\\beta_k - 2\\Gamma^k_{ij}\\beta_k \\right) \\\\\n",
    "&= \\frac{\\sqrt{\\gamma}}{2\\alpha} T^{ij} \\left( 2\\alpha K_{ij} + \\partial_j \\beta_i - \\partial_i\\beta_j \\right) \\\\\n",
    "&= \\sqrt{\\gamma} T^{ij} K_{ij} \\\\\n",
    "\\end{align}\n",
    "\n",
    "where the final step again follows from the symmetry of $T^{ij}$. Putting all this together and multiplying through by $\\alpha$,\n",
    "\n",
    "\\begin{align}\n",
    "\\alpha\\partial_t \\left(\\alpha\\sqrt{\\gamma} T^{00} \\right) + \\alpha\\partial_i \\left(\\alpha\\sqrt{\\gamma} T^{0i} \\right) &= \\alpha\\sqrt{\\gamma} T^{00} \\left( \\beta^i \\beta^j K_{ij} - \\partial_t \\alpha - \\beta^i \\partial_i \\alpha \\right) \\\\\n",
    "&\\ \\ \\ \\ \\ \\ \\ + 2\\alpha\\sqrt{\\gamma} T^{0i} \\left( \\beta^j K_{ij} - \\partial_i \\alpha \\right) + \\alpha\\sqrt{\\gamma} T^{ij} K_{ij} \\\\\n",
    "%\n",
    "\\alpha\\sqrt{\\gamma} T^{00}\\partial_t \\alpha + 2\\alpha\\sqrt{\\gamma} T^{0i}\\partial_i \\alpha + \\alpha\\partial_t \\left(\\alpha\\sqrt{\\gamma} T^{00} \\right) + \\alpha\\partial_i \\left(\\alpha\\sqrt{\\gamma} T^{0i} \\right) &= \\alpha\\sqrt{\\gamma} T^{00} \\left( \\beta^i \\beta^j K_{ij} - \\beta^i \\partial_i \\alpha \\right) + 2\\alpha\\sqrt{\\gamma} T^{0i}\\beta^j K_{ij} \\\\\n",
    "&\\ \\ \\ \\ \\ \\ \\ - \\alpha\\sqrt{\\gamma} T^{0i} \\partial_i \\alpha + \\alpha\\sqrt{\\gamma} T^{ij} K_{ij} \\\\\n",
    "%\n",
    "\\partial_t \\left(\\alpha^2\\sqrt{\\gamma} T^{00} \\right) + \\partial_i \\left(\\alpha^2\\sqrt{\\gamma} T^{0i} \\right) &= \\alpha\\sqrt{\\gamma} \\left[ \\left( T^{00} \\beta^i \\beta^j + 2T^{0i}\\beta^j + T^{ij} \\right) K_{ij} - \\left( T^{00} \\beta^i + T^{0i} \\right)\\partial_i \\alpha \\right]\n",
    "\\end{align}\n",
    "\n",
    "We have finally arrived at the evolution equation. However, we still need to get it in terms of the evolution variable $\\tilde{\\tau} = \\alpha^2 \\sqrt{\\gamma} T^{00} - \\tilde{\\rho}$. To do so, we can add the fluid density evolution equation to the left-hand side. We also define the source term\n",
    "\n",
    "$$\n",
    "s = \\alpha\\sqrt{\\gamma} \\left[ \\left( T^{00} \\beta^i \\beta^j + 2T^{0i}\\beta^j + T^{ij} \\right) K_{ij} - \\left( T^{00} \\beta^i + T^{0i} \\right)\\partial_i \\alpha \\right]\n",
    "$$\n",
    "\n",
    "Then, the energy evolution equation becomes\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_t \\left(\\alpha^2\\sqrt{\\gamma} T^{00} \\right) + \\partial_i \\left(\\alpha^2\\sqrt{\\gamma} T^{0i} \\right) - \\partial_t \\tilde{\\rho} - \\partial_i \\left(\\tilde{\\rho} v^i \\right) &= s \\\\\n",
    "\\partial_t \\tilde{\\tau} + \\partial_i \\left(\\alpha^2\\sqrt{\\gamma} T^{0i} - \\tilde{\\rho} v^i \\right) &= s\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='induction_eqn'></a>\n",
    "\n",
    "## Step 2.4: Evolution Equation for the Magnetic Field $\\tilde{B}^i$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{induction_eqn}$$\n",
    "\n",
    "\n",
    "To find a convenient start for the derivation, I find the dual of Maxwell's equation. Recall that the dual of the Faraday tensor is\n",
    "\n",
    "$$\n",
    "F^{*\\mu\\nu} = \\frac{1}{2} \\epsilon^{\\mu\\nu\\alpha\\beta}F_{\\alpha\\beta}\n",
    "$$\n",
    "\n",
    "Also, since the covariant derivative of the metric is 0,\n",
    "\n",
    "$$\n",
    "\\epsilon^{\\alpha\\lambda\\mu\\nu} F_{\\mu\\nu;\\lambda} = \\nabla_\\lambda \\left(\\epsilon^{\\alpha\\lambda\\mu\\nu} F_{\\mu\\nu}\\right) = 2\\nabla_\\lambda F^{*\\alpha\\lambda}\n",
    "$$\n",
    "\n",
    "Therefore, Maxwell's equation becomes\n",
    "\n",
    "\\begin{align}\n",
    "0 &= \\epsilon^{\\alpha\\lambda\\mu\\nu} F_{[\\mu\\nu;\\lambda]} \\\\\n",
    "&= \\epsilon^{\\alpha\\lambda\\mu\\nu} \\left( F_{\\mu\\nu;\\lambda} - F_{\\lambda\\nu;\\mu} + F_{\\nu\\lambda;\\mu} - F_{\\mu\\lambda;\\nu} + F_{\\lambda\\mu;\\nu} - F_{\\nu\\mu;\\lambda} \\right) \\\\\n",
    "%\n",
    "&= \\epsilon^{\\alpha\\lambda\\mu\\nu}F_{\\mu\\nu;\\lambda} + \\epsilon^{\\alpha\\lambda\\mu\\nu} \\nabla_\\mu \\left( F_{\\nu\\lambda} - F_{\\lambda\\nu} \\right) + \\epsilon^{\\alpha\\lambda\\mu\\nu} \\nabla_\\nu \\left( F_{\\lambda\\mu} - F_{\\mu\\lambda}\\right) - \\epsilon^{\\alpha\\lambda\\mu\\nu}\\nabla_\\lambda F_{\\nu\\mu} \\\\\n",
    "%\n",
    "&= 2\\nabla_\\lambda F^{*\\alpha\\lambda} + 2\\epsilon^{\\alpha\\lambda\\mu\\nu} \\nabla_\\mu F_{\\nu\\lambda} + 2\\epsilon^{\\alpha\\lambda\\mu\\nu} \\nabla_\\nu F_{\\lambda\\mu} - \\epsilon^{\\alpha\\lambda\\mu\\nu}\\nabla_\\lambda F_{\\nu\\mu}\n",
    "\\end{align}\n",
    "\n",
    "By exploiting the symmetries of the Levi-Civita tensor,\n",
    "\n",
    "\\begin{align}\n",
    "0 &= 2\\nabla_\\lambda F^{*\\alpha\\lambda} + 2\\epsilon^{\\alpha\\mu\\nu\\lambda} \\nabla_\\mu F_{\\nu\\lambda} + 2\\epsilon^{\\alpha\\nu\\lambda\\mu} \\nabla_\\nu F_{\\lambda\\mu} + \\epsilon^{\\alpha\\lambda\\nu\\mu}\\nabla_\\lambda F_{\\nu\\mu} \\\\\n",
    "%\n",
    "&= 2\\nabla_\\lambda F^{*\\alpha\\lambda} + 4\\nabla_\\mu F^{*\\alpha\\mu} + 4\\nabla_\\nu F^{*\\alpha\\nu} + 2\\nabla_\\lambda F^{*\\alpha\\lambda} \\\\\n",
    "%\n",
    "&= 12\\nabla_\\lambda F^{*\\alpha\\lambda} \\\\\n",
    "0 &= \\nabla_\\lambda F^{*\\alpha\\lambda} = \\partial_\\lambda F^{*\\alpha\\lambda} + \\Gamma^\\alpha_{\\beta\\lambda} F^{*\\beta\\lambda} + \\Gamma^\\lambda_{\\beta\\lambda} F^{*\\alpha\\beta} \\\\\n",
    "0 &= \\frac{1}{\\sqrt{-g}}\\partial_\\lambda \\left(\\sqrt{-g} F^{*\\alpha\\lambda}\\right)\n",
    "\\end{align}\n",
    "\n",
    "where the first $\\Gamma$ term dissapears because it is summing over the multiplication of symmetric and anti-symmetric objects.\n",
    "\n",
    "We can now find the magnetic equations from this. Taking the time component simply gives us the no-monopole constraint. To see this, we need to first consider the individual components of the dual. By the definition of $B$,\n",
    "\n",
    "$$\n",
    "B^\\mu = \\frac{1}{2} \\epsilon^{\\mu\\nu\\beta\\alpha}n_\\nu F_{\\alpha\\beta} = n_\\nu F^{*\\nu\\mu}\n",
    "$$\n",
    "\n",
    "Since $B$ for the normal observer is purely spatial ($B^\\mu n_\\mu=0$), this implies that $F^{*00}=0$. Then, the time component of the magnetic equations is\n",
    "\n",
    "\\begin{align}\n",
    "0 &= \\frac{1}{\\sqrt{-g}}\\partial_\\lambda \\left(\\sqrt{-g} F^{*0\\lambda}\\right) \\\\\n",
    "&= \\partial_i \\left(\\alpha\\sqrt{\\gamma} F^{*0i}\\right) \\\\\n",
    "&= \\partial_i \\left(\\alpha\\sqrt{\\gamma} \\frac{B^i}{\\alpha}\\right) \\\\\n",
    "0 &= \\partial_i \\left(\\tilde{B}^i \\right)\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "Before examining the spatial components of Maxwell's equations, we will do some preliminary work to make our lives easier later. In the co-moving frame,\n",
    "\n",
    "\\begin{align}\n",
    "F^{\\mu\\nu} &= u^\\mu E^\\nu - u^\\nu E^\\mu + u_\\gamma \\epsilon^{\\gamma\\mu\\nu\\delta} B_{\\delta} \\\\\n",
    "&= u_\\gamma \\epsilon^{\\gamma\\mu\\nu\\delta} B_{\\delta}\n",
    "\\end{align}\n",
    "\n",
    "Then, the dual is\n",
    "\n",
    "\\begin{align}\n",
    "F^{*\\mu\\nu} &= \\frac{1}{2} \\epsilon^{\\mu\\nu\\alpha\\beta} F_{\\alpha\\beta} \\\\\n",
    "&= \\frac{1}{2} \\epsilon^{\\mu\\nu\\alpha\\beta} u^\\gamma \\epsilon_{\\gamma\\alpha\\beta\\delta} B^{\\delta} \\\\\n",
    "&= \\frac{1}{2} \\epsilon^{\\alpha\\beta\\mu\\nu}\\epsilon_{\\alpha\\beta\\gamma\\delta} u^\\gamma  B^{\\delta} \\\\\n",
    "&= \\left( \\delta^\\mu_\\delta \\delta^\\nu_\\gamma - \\delta^\\mu_\\gamma \\delta^\\nu_\\delta \\right) u^\\gamma  B^{\\delta} \\\\\n",
    "&= u^\\nu  B^{\\mu} - u^\\mu  B^{\\nu}\n",
    "\\end{align}\n",
    "\n",
    "where we have again used the Levi-Civita identity from [Appendix A](#levi_civita). Next, we need to find a relationship between the magnetic field of the normal observer and the co-moving observer. For clarity, let the co-moving magnetic field be $B^\\mu_{(u)}$ and the normal magnetic field be $B^\\mu$. Then, we can define a projection operator\n",
    "\n",
    "$$\n",
    "P^{\\mu\\nu} = g_{\\mu\\nu} + u_\\mu u_\\nu\n",
    "$$\n",
    "\n",
    "Naturally, the projection of the co-moving magnetic field should simply project back into the same field:\n",
    "\n",
    "\\begin{align}\n",
    "P^\\mu_\\nu B^\\nu_{(u)} &= \\left( \\delta^\\mu_\\nu + u^\\mu u_\\nu \\right)B^\\nu_{(u)} \\\\\n",
    "&= B^\\mu_{(u)}\n",
    "\\end{align}\n",
    "\n",
    "where we have used the orthogonality relation $u_\\nu B^\\nu_{(u)}=0$. Projecting the normal observer's magnetic field,\n",
    "\n",
    "\\begin{align}\n",
    "P^\\mu_\\nu B^\\nu &= P^\\mu_\\nu n_\\alpha F^{*\\alpha\\nu} \\\\\n",
    "&= P^\\mu_\\nu n_\\alpha \\left( u^\\nu  B^{\\alpha}_{(u)} - u^\\alpha  B^{\\nu}_{(u)} \\right) \\\\\n",
    "&= n_\\alpha \\left( \\delta^\\mu_\\nu + u^\\mu u_\\nu \\right) \\left( u^\\nu  B^{\\alpha}_{(u)} - u^\\alpha  B^{\\nu}_{(u)} \\right) \\\\\n",
    "%\n",
    "&= n_\\alpha \\left( u^\\mu  B^{\\alpha}_{(u)} +  u^\\mu u_\\nu u^\\nu  B^{\\alpha}_{(u)} - u^\\alpha  B^{\\mu}_{(u)} - u^\\mu u_\\nu u^\\alpha  B^{\\nu}_{(u)} \\right) \\\\\n",
    "%\n",
    "&= n_\\alpha u^\\mu B^{\\alpha}_{(u)} \\left( 1 + u_\\nu u^\\nu \\right) - n_\\alpha u^\\alpha \\left( B^{\\mu}_{(u)} + u^\\mu u_\\nu B^{\\nu}_{(u)} \\right) \\\\\n",
    "%\n",
    "&= - \\alpha u^0 B^{\\mu}_{(u)} \\\\\n",
    "\\Rightarrow B^{\\mu}_{(u)} &= \\frac{B^\\mu + u^\\mu u_\\nu B^\\nu}{\\alpha u^0}\n",
    "\\end{align}\n",
    "\n",
    "Since we will only need the spatial component for our purposes,\n",
    "\n",
    "\\begin{align}\n",
    "B^{i}_{(u)} &= \\frac{B^i + u^i u_\\nu B^\\nu}{\\alpha u^0} \\\\\n",
    "&= \\frac{B^i + u^i u_j B^j}{\\alpha u^0} \\\\\n",
    "&= \\frac{B^i}{\\alpha u^0} + \\frac{v^i u_j B^j}{\\alpha}\n",
    "\\end{align}\n",
    "\n",
    "where $v^i$ is the conservative variable $u^i/u^0$. Finally, the spatial components of the dual of Maxwell's equations gives the induction equation for the magnetic field:\n",
    "\n",
    "\\begin{align}\n",
    "0 &= \\frac{1}{\\sqrt{-g}}\\partial_\\nu \\left(\\sqrt{-g} F^{*i \\nu}\\right) \\\\\n",
    "0 &= \\partial_t \\left(\\alpha\\sqrt{\\gamma} F^{*i 0}\\right) + \\partial_j \\left(\\alpha\\sqrt{\\gamma} F^{*i j}\\right) \\\\\n",
    "%\n",
    "0 &= \\partial_t \\left(\\sqrt{\\gamma} B^i \\right) + \\partial_j \\left(\\alpha\\sqrt{\\gamma} \\left[ u^j B^i_{(u)} - u^i B^j_{(u)} \\right] \\right) \\\\\n",
    "%\n",
    "0 &= \\partial_t \\tilde{B}^i + \\partial_j \\left(\\alpha\\sqrt{\\gamma} \\left[ u^j \\left( \\frac{B^i}{\\alpha u^0} + \\frac{v^i u_k B^k}{\\alpha} \\right) - u^i \\left( \\frac{B^j}{\\alpha u^0} + \\frac{v^j u_k B^k}{\\alpha} \\right) \\right] \\right) \\\\\n",
    "%\n",
    "0 &= \\partial_t \\tilde{B}^i  + \\partial_j \\left(v^j \\tilde{B}^i - v^i \\tilde{B}^j + \\frac{\\sqrt{\\gamma}}{u^0}\\left[ u^j u^i u_k B^k - u^i u^j u_k B^k \\right] \\right) \\\\\n",
    "%\n",
    "0 &= \\partial_t \\tilde{B}^i  + \\partial_j \\left(v^j \\tilde{B}^i - v^i \\tilde{B}^j \\right)\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='summary'></a>\n",
    "\n",
    "## Step 2.5: Summary of the Conservative Variable Evolution Equations \\[Back to [top](#toc)\\]\n",
    "$$\\label{summary}$$\n",
    "\n",
    "In the previous sections, we have derived the evolution equations for the conservative variables \\mathbf{C}. To summarize, these are\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_t \\tilde{\\rho} + \\partial_i \\left(\\tilde{\\rho} v^i \\right) &= 0 \\\\\n",
    "\\partial_t \\tilde{\\tau} + \\partial_i \\left(\\alpha^2\\sqrt{\\gamma} T^{0i} - \\tilde{\\rho} v^i \\right) &= s \\\\\n",
    "\\partial_t \\tilde{S}_i + \\partial_j \\left(\\alpha\\sqrt{\\gamma} T_i^j \\right) &= \\frac{\\alpha\\sqrt{\\gamma}}{2} T^{\\beta\\nu}g_{\\beta\\nu,i} \\\\\n",
    "\\partial_t \\tilde{B}^i  + \\partial_j \\left(v^j \\tilde{B}^i - v^i \\tilde{B}^j \\right) &= 0\n",
    "\\end{align}\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "s = \\alpha\\sqrt{\\gamma} \\left[ \\left( T^{00} \\beta^i \\beta^j + 2T^{0i}\\beta^j + T^{ij} \\right) K_{ij} - \\left( T^{00} \\beta^i + T^{0i} \\right)\\partial_i \\alpha \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='levi_civita'></a>\n",
    "\n",
    "# Appendix A: Levi-Civita Contractions \\[Back to [top](#toc)\\]\n",
    "$$\\label{levi_civita}$$\n",
    "\n",
    "In order to simplify the various combinations of $F^{\\mu\\nu}$, contractions of the Levi-Civita tensor are required. First, summing over the first indices of two such tensors yields\n",
    "\n",
    "\\begin{align}\n",
    "\\epsilon^{\\alpha\\beta\\sigma\\mu} \\epsilon_{\\alpha\\beta\\gamma\\delta} &= -\\delta^{\\beta\\sigma\\mu}_{\\hphantom{\\beta\\sigma\\nu}\\beta\\gamma\\delta} \\\\\n",
    "%\n",
    "&= -\\delta^\\beta_\\beta \\delta^{\\sigma\\mu}_{\\gamma\\delta} \n",
    "  + \\delta^\\beta_\\gamma \\delta^{\\sigma\\mu}_{\\beta\\delta}\n",
    "  - \\delta^\\beta_\\delta \\delta^{\\sigma\\mu}_{\\beta\\gamma} \\\\\n",
    "%\n",
    "&= -4 \\left(\\delta^\\sigma_\\gamma \\delta^\\mu_\\delta - \\delta^\\sigma_\\delta \\delta^\\mu_\\gamma\\right)\n",
    "  + \\delta^\\beta_\\gamma \\left(\\delta^\\sigma_\\beta \\delta^\\mu_\\delta - \\delta^\\sigma_\\delta \\delta^\\mu_\\beta\\right)\n",
    "  - \\delta^\\beta_\\delta \\left(\\delta^\\sigma_\\beta \\delta^\\mu_\\gamma - \\delta^\\sigma_\\gamma \\delta^\\mu_\\beta\\right) \\\\\n",
    "%\n",
    "&= -4 \\left(\\delta^\\sigma_\\gamma \\delta^\\mu_\\delta - \\delta^\\sigma_\\delta \\delta^\\mu_\\gamma\\right)\n",
    "  + 2\\left(\\delta^\\sigma_\\gamma \\delta^\\mu_\\delta - \\delta^\\sigma_\\delta \\delta^\\mu_\\gamma\\right) \\\\\n",
    "%\n",
    "&= 2\\left(\\delta^\\sigma_\\delta \\delta^\\mu_\\gamma - \\delta^\\sigma_\\gamma \\delta^\\mu_\\delta\\right)\n",
    "\\end{align}\n",
    "\n",
    "Second, summing over the third index (with all other components raised) yields\n",
    "\n",
    "\\begin{align}\n",
    "\\epsilon^{\\gamma\\mu\\lambda\\delta} \\epsilon^{\\alpha\\nu\\hphantom{\\lambda}\\beta}_{\\hphantom{\\alpha\\nu}\\lambda} &= g^{\\alpha\\tau}g^{\\nu\\phi}g^{\\beta\\theta} \\epsilon^{\\gamma\\mu\\lambda\\delta} \\epsilon_{\\tau\\phi\\lambda\\theta} \\\\\n",
    "%\n",
    "&= -g^{\\alpha\\tau}g^{\\nu\\phi}g^{\\beta\\theta} \\delta^{\\gamma\\mu\\delta}_{\\tau\\phi\\theta} \\\\\n",
    "%\n",
    "&= -g^{\\alpha\\tau}g^{\\nu\\phi}g^{\\beta\\theta}\\left( \\delta^\\gamma_\\tau\\delta^{\\mu\\delta}_{\\phi\\theta}\n",
    "  - \\delta^\\gamma_\\phi \\delta^{\\mu\\delta}_{\\tau\\theta}\n",
    "  + \\delta^\\gamma_\\theta \\delta^{\\mu\\delta}_{\\tau\\phi} \\right) \\\\\n",
    "%\n",
    "&= -g^{\\alpha\\tau}g^{\\nu\\phi}g^{\\beta\\theta}\n",
    "  \\left( \\delta^\\gamma_\\tau (\\delta^\\mu_\\phi \\delta^\\delta_\\theta - \\delta^\\mu_\\theta \\delta^\\delta_\\phi)\n",
    "  - \\delta^\\gamma_\\phi (\\delta^\\mu_\\tau \\delta^\\delta_\\theta - \\delta^\\mu_\\theta \\delta^\\delta_\\tau)\n",
    "  + \\delta^\\gamma_\\theta (\\delta^\\mu_\\tau  \\delta^\\delta_\\phi - \\delta^\\mu_\\phi \\delta^\\delta_\\tau) \\right) \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='curv_tensor'></a>\n",
    "\n",
    "# Appendix B: Extrinsic Curvature Identities \\[Back to [top](#toc)\\]\n",
    "$$\\label{curv_tensor}$$\n",
    "\n",
    "The source term of the energy equation involves the extrinsic curvature. Defined in terms of the 3+1 formalism quantities, the extrinsic curvature is\n",
    "\\begin{align}\n",
    "K_{ij} &= -\\frac{1}{2\\alpha}\\left( \\partial_t \\gamma_{ij} - \\nabla_i\\beta_j - \\nabla_j\\beta_i \\right) \\\\\n",
    "&= -\\frac{1}{2\\alpha}\\left( \\partial_t \\gamma_{ij} - \\partial_i\\beta_j - \\partial_j\\beta_i + 2\\Gamma^k_{ij}\\beta_k \\right)\n",
    "\\end{align}\n",
    "\n",
    "As for how the extrinsic curvature enters into the energy equations, it involves several relations between it and the other quantities which we must show. First, consider the contraction of the Christoffel symbol with $\\beta^i \\beta^j$:\n",
    "\n",
    "\\begin{align}\n",
    "\\beta^j \\Gamma^k_{ij}\\beta_k &= \\beta^j \\beta^k \\left( \\gamma_{k i,j} + \\gamma_{k j,i} - \\gamma_{i j,k} \\right) \\\\\n",
    "&= \\beta^j \\beta^k \\gamma_{k j,i} \\\\\n",
    "&= \\alpha^2 \\left( \\gamma^{jk} - g^{jk} \\right) \\gamma_{k j,i} \\\\\n",
    "&\\propto \\gamma^{jk}\\partial_i \\gamma_{k j} - g_{jk} \\partial_i \\gamma^{k j} \\\\\n",
    "&\\propto \\gamma^{jk}\\partial_i \\gamma_{k j} - \\gamma_{jk} \\partial_i \\gamma^{k j} \\\\\n",
    "&\\propto \\gamma^{jk} \\partial_i \\gamma_{k j} - \\gamma^{jk} \\partial_i \\gamma_{k j} \\\\\n",
    "&= 0\n",
    "\\end{align}\n",
    "\n",
    "The change $g_{jk}\\rightarrow\\gamma_{jk}$ is allowed because the 3-metric is identical to the 4-metric with only spatial indices. We can use this to find a relationship between the partial derivative of $\\beta_j$ and the extrinsic curvature:\n",
    "\n",
    "\\begin{align}\n",
    "\\beta^i \\beta^j \\partial_i \\beta_j &= \\beta^i \\beta^j \\left( 2\\alpha K_{ij} + \\partial_t \\gamma_{ij} - \\partial_j\\beta_i + 2\\Gamma^k_{ij}\\beta_k \\right) \\\\\n",
    "&= \\beta^i \\beta^j \\left( 2\\alpha K_{ij} + \\partial_t \\gamma_{ij} - \\partial_i\\beta_j + 2\\Gamma^k_{ij}\\beta_k \\right) \\\\\n",
    "&= \\beta^i \\beta^j \\left( 2\\alpha K_{ij} + \\partial_t \\gamma_{ij} - \\partial_i\\beta_j \\right) \\\\\n",
    "\\Rightarrow \\beta^i \\beta^j \\partial_i \\beta_j &= \\beta^i \\beta^j \\left( \\alpha K_{ij} + \\frac{1}{2}\\partial_t \\gamma_{ij} \\right) \\\\\n",
    "\\Rightarrow \\beta^i \\beta^j \\partial_i \\beta_j &= \\beta^i \\beta^j \\alpha K_{ij}\n",
    "\\end{align}\n",
    "\n",
    "where in the second step we use the symmetry of $\\beta^i \\beta^j$. The final step uses the same trick as with $\\beta^j \\Gamma^k_{ij}\\beta_k$, just with the derivative being $\\partial_t$ instead of $\\partial_i$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 3: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Derivation_of_GRMHD_Evolution_Equations.tex, and compiled\n",
      "    LaTeX file to PDF file Tutorial-\n",
      "    Derivation_of_GRMHD_Evolution_Equations.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Derivation_of_GRMHD_Evolution_Equations\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}