IllinoisGRMHD: Basic equations and modules¶

Authors: Leo Werneck & Zach Etienne¶

This module is currently under development

This module introduces the basic equations solved by IllinoisGRMHD and provides an introduction to the modules inside the code.¶

Introduction:¶

IllinoisGRMHD solves the equations of General Relativistic MagnetoHydroDynamics (GRMHD) using a high-resolution shock capturing scheme. It is a rewrite of the Illinois Numerical Relativity (ILNR) group's GRMHD code, and generates results that agree to roundoff error with that original code. Its feature set coincides with the features of the ILNR group's recent code (ca. 2009-2014), which was used in their modeling of the following systems:

1. Magnetized circumbinary disk accretion onto binary black holes
2. Magnetized black hole-neutron star mergers
3. Magnetized Bondi flow, Bondi-Hoyle-Littleton accretion
4. White dwarf-neutron star mergers

IllinoisGRMHD is particularly good at modeling GRMHD flows into black holes without the need for excision. Its HARM-based conservative-to-primitive solver has also been modified to check the physicality of conservative variables prior to primitive inversion, and move them into the physical range if they become unphysical.

Currently IllinoisGRMHD consists of

1. the Piecewise Parabolic Method (PPM) for reconstruction,
2. the Harten, Lax, van Leer (HLL/HLLE) approximate Riemann solver, and
3. a modified HARM Conservative-to-Primitive solver.

IllinoisGRMHD evolves the vector potential $A_{\mu}$ (on staggered grids) instead of the magnetic fields ($B^i$) directly, to guarantee that the magnetic fields will remain divergenceless even at AMR boundaries. On uniform resolution grids, this vector potential formulation produces results equivalent to those generated using the standard, staggered flux-CT scheme. This scheme is based on that of Del Zanna et al. (2003).

Required and recommended citations:¶

• (Required) Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. (arxiv:1501.07276).
• (Required) Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) (astro-ph/0512420).
• (Recommended) Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 (astro-ph/0210618).

Table of Contents¶

$$\label{toc}$$

This module is organized as follows

1. Step 1: The equations of GRMHD
2. Step 2: Recasting the equations of GRMHD into a more useful form
   1. Step 2.1: Einstein's equations
   2. Step 2.2: Conservation of baryon number
   3. Step 2.3: Conservation of energy-momentum
   4. Step 2.4: The evolution equations for the electromagnetic field
   5. Step 2.5: The equation of state (EOS)
3. Step 3: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: The equations of GRMHD [Back to top]¶

$$\label{basic_equations}$$

The basic equations solved by IllinoisGRMHD are

1. Einstein's equations: $R^{\mu\nu} + \frac{1}{2}g^{\mu\nu}R = 8\pi T^{\mu\nu} = 8\pi \left(T^{\mu\nu}_{\rm matter} + T^{\mu\nu}_{\rm EM}\right)$ ;
2. Conservation of baryon number: $\nabla_{\mu}\left(\rho_{0} u^{\mu}\right) = 0$ ;
3. Conservation of energy-momentum: $\nabla_{\mu}T^{\mu\nu} = \nabla_{\mu}\left(T^{\mu\nu}_{\rm matter} + T^{\mu\nu}_{\rm EM}\right) = 0$ ;
4. Dual of homogeneous Maxwell's equations: $\nabla_{\nu} F^{*\mu\nu}=\frac{1}{\sqrt{-g}}\partial_{\nu}\left(\sqrt{-g}F^{*\mu\nu}\right)=0$ ,

where $g_{\mu\nu}$ is the ADM 4-metric, $R_{\mu\nu}$ and $R$ are the Ricci tensor and scalar, respectively, constructed from $g_{\mu\nu}$, $\nabla_{\mu}$ is the covariant derivative compatible with $g_{\mu\nu}$, $g\equiv\det\left(g_{\mu\nu}\right)$, $F^{\mu\nu}$ is the Faraday tensor and $F^{*\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}$ its dual ($\epsilon^{\mu\nu\rho\sigma}$ is the Levi-Civita symbol).

The final equation that must be solved is the equation of state (EOS) for the matter. Currently IllinoisGRMHD implements a hybrid EOS of the form

5. Equation of State (EOS): $P(\rho_{0},\epsilon) = P_{\rm cold}(\rho_{0}) + \left(\Gamma_{\rm th}-1\right) \rho_{0} \left[\epsilon-\epsilon_{\rm cold}(\rho_{0})\right]$ ,

where P is the pressure, $\epsilon$ the specific internal energy (the subscript cold indicate the cold components of these quantities), and $\Gamma_{\rm th}$ is a constant parameter which determines the conversion efficiency of kinetic to thermal energy at shocks.

Step 2: Recasting the equations of GRMHD into a more useful form [Back to top]¶

$$\label{suitable_equations}$$

In this step we will write down the equations used by IllinoisGRMHD in the form they are implemented. To give an example, the GRMHD equations are written in the conservative form

$$\partial_{t}\boldsymbol{C} + \boldsymbol{\nabla}\cdot\boldsymbol{F} = \boldsymbol{S}\ ,$$

where $\boldsymbol{C} = \left\{\rho_{\star},\tilde\tau,\tilde{S}_{i},\tilde{B}^{i}\right\}$ is the vector of conservative variables, $\boldsymbol{F}$ is the flux vector, and $\boldsymbol{S}$ the vector of source terms (explict expressions for the components of $\boldsymbol{C}$, $\boldsymbol{F}$, and $\boldsymbol{S}$ can be found below).

Step 2.1: Einstein's equations [Back to top]¶

$$\label{einstein_equations}$$

IllinoisGRMHD solves Einstein's field equations (with $G_{\rm N}=1=c$) in the presence of matter sources,

$$G^{\mu\nu} = 8\pi T^{\mu\nu}\ ,$$

where $G^{\mu\nu} \equiv R^{\mu\nu} + \frac{1}{2}g^{\mu\nu}R$ is the Einstein tensor and the total energy-momentum tensor, $T^{\mu\nu}$, is the sum of the matter and electromagnetic energy-momentum tensors:

$$T^{\mu\nu} = T^{\mu\nu}_{\rm matter} + T^{\mu\nu}_{\rm EM}\ .$$

Einstein's field equations are solved using the BSSN formalism (see this tutorial module for an overview).

Step 2.2: Conservation of baryon number [Back to top]¶

$$\label{conserv_of_baryon_number}$$

The conservation of baryon number equation is written as (cf. eqs. (6), (7), and (18) in Etienne et al.)

$$\boxed{\partial_{t}\rho_{\star} + \partial_{j}\left(\rho_{\star}v^{j}\right) = 0}\ ,$$

where $\rho_{\star}\equiv \alpha \sqrt{\gamma}\rho_{0}u^{0}$ and $v^{i}\equiv u^{i}/u^{0}$.

Step 2.3: Conservation of energy-momentum [Back to top]¶

$$\label{conserv_of_em_tensor}$$

In the ideal MHD limit, we can write down the total energy momentum tensor as (cf. eq. (8) in Etienne et al. and the discussion before it)

$$T^{\mu \nu} = (\rho_0 h +b^2) u^{\mu} u^{\nu} + \left( P + \frac{b^2}{2}\right) g^{\mu \nu} - b^{\mu} b^{\nu}$$

The spatial components of the energy-momentum conservation equation give (cf. eq. (18) of Etienne et al. and eqs. (35) and (36) of Duez et al.)

$$\boxed{\partial_{t}\tilde{S}_{i} + \partial_{j}\left(\alpha\sqrt{\gamma} T^{j}_{\ i}\right) = \frac{1}{2}\sqrt{\gamma} T^{\alpha\beta}\partial_{i}g_{\alpha\beta}}\ ,$$

where

$$\tilde{S}_{i} = \sqrt{\gamma}S_{i} = \alpha \sqrt{\gamma} T^{0}_{\ i} = \left(\rho_{\star}h + \alpha u^{0}\sqrt{\gamma}b^{2}\right)u_{i} - \alpha \sqrt{\gamma} b^{0} b^{i}\ .$$

The time component of the energy-momentum conservation equation gives (cf. eq. (18) of Etienne et al.)

$$\boxed{\partial_{t}\tilde{\tau} + \partial_{i}\left(\alpha^{2}\sqrt{\gamma}T^{0i} - \rho_{\star}v^{i}\right) = s}\ ,$$

where

$$\begin{align} \tilde{\tau} &= \sqrt{\gamma}n_{\mu}n_{\nu} - \rho_{\star} = \alpha^{2}\sqrt{\gamma}T^{00} - \rho_{\star}\ ,\\ s &= -\alpha \sqrt{\gamma} T^{\mu\nu} \nabla_{\nu} n_{\mu} = \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]\ , \end{align}$$

with $n_{\mu} = \left(\alpha,0,0,0\right)$ being the normal vector and $K_{ij}$ the extrinsic curvature.

Step 2.4: The evolution equations for the electromagnetic field [Back to top]¶

$$\label{em_evo_equations}$$

From the spatial components of the dual of Maxwell's equations,

$$\nabla_{\nu} F^{*\mu\nu}=\frac{1}{\sqrt{-g}}\partial_{\nu}\left(\sqrt{-g}F^{*\mu\nu}\right)=0\ ,$$

we get the magnetic induction equation, which in conservative form may be written as (cf. eq. (12) Etienne et al.)

$$\partial_{t}\tilde{B}^{i} + \partial_{j}\left(v^{j}\tilde{B}^{i} - v^{i}\tilde{B}^{j}\right) = 0\ ,$$

where $\tilde{B}^{i} = \sqrt{\gamma}B^{i}$. We must also guarantee that no magnetic mopoles form, via the constraint

$$\partial_{i}\tilde{B}^{i} = 0\ ,$$

which is the time component of the dual of Maxwell's equations. Satisfying this constraint equation while evolving the magnetic field forward in time via the evolution equation above turns out to be a nontrivial endeavor, particularly on AMR grids. Instead, we choose to evolve the magnetic 4-vector potential $\mathcal{A}_{\mu}$ instead of the magnetic fields directly, so that

$$\begin{align} \mathcal{A}_{\mu} &= \Phi n_{\mu} + A_{\mu}\ ,\\ \tilde{B}^{i} &= \tilde{\epsilon}^{ijk}\partial_{j}A_{k}\ , \end{align}$$

where $A_{\mu}$ is purely spatial ($A_{\mu}n^{\mu} = 0$) and $\Phi$ is the EM scalar potential.

Special finite difference operators for the vector potential are defined in IllinoisGRMHD so that the divergence of a curl is zero to roundoff error, which implies that the divergence of $\tilde{B}^{i}$ (as defined above) is zero and the condition $\partial_{i}\tilde{B}^{i}=0$ is satisfied automatically, even on AMR grids.

In terms of $A_{i}$, the induction equation becomes

$$\boxed{\partial_t A_i = \tilde\epsilon_{ijk} v^j \tilde{B}^k - \partial_i (\alpha \Phi - \beta^j A_j)}\ .$$

The final equation comes from choosing the covariant version of the "generalized Lorenz gauge condition" that was introduced by the Illinois Relativity group,

$$\nabla_{\mu}\mathcal{A}^{\mu} = \xi n_{\mu} \mathcal{A}^{\mu}\ ,$$

where $\xi$ is a parameter with dimensions 1/Length, chosen so that the CFL condition remains satisfied. This gauge choice yields the additional evolution equation

$$\boxed{\partial_t [\sqrt{\gamma} \Phi] + \partial_j (\alpha \sqrt{\gamma} A^j - \beta^j [\sqrt{\gamma}\Phi]) = -\xi\alpha \sqrt{\gamma}\Phi}\ .$$

Note that in IllinoisGRMHD the evolved variable is $\sqrt{\gamma}\Phi$, not $\Phi$.

Step 2.5: The equation of state (EOS) [Back to top]¶

$$\label{equation_of_state}$$

IllinoisGRMHD currently implements a hybrid EOS of the form

$$\boxed{P(\rho_{0},\epsilon) = P_{\rm cold}(\rho_{0}) + \left(\Gamma_{\rm th}-1\right) \rho_{0} \left[\epsilon-\epsilon_{\rm cold}(\rho_{0})\right]}\ .$$

The function $\epsilon_{\rm cold}$ is related to $P_{\rm cold}$ by the fist law of thermodynamics,

$$\epsilon_{\rm cold}(\rho_{0}) = \int \frac{P_{\rm cold}(\rho_{0})}{\rho_{0}^{2}}d\rho_{0}\ .$$

Currently, all functions within IllinoisGRMHD support piecewise-defined $P_{\rm cold}(\rho_{0})$ (the so-called "piecewise polytrope" EOS) with up to nine different polytropic indices.

In typical runs of the code, particularly on the ones presented in the release paper, the $\Gamma$-law EOS, $P = (\Gamma - 1)\rho_{0}\epsilon$, is adopted. This corresponds to setting $P_{\rm cold} = (\Gamma - 1)\rho_{0}\epsilon_{\rm cold}$ in the boxed equation above, which is equivalent to $P_{\rm cold} = \kappa\rho_{0}^{\Gamma}$ (with constant $\kappa$), and $\Gamma_{\rm th} = \Gamma$. In the absence of shocks, $\epsilon = \epsilon_{\rm cold}$ so that $P = P_{\rm cold}$.

IllinoisGRMHD has also performed successful binary neutron star mergers using the widely adopted piecewise polytropic EOSs of Read et al.. In the EOS low level functions tutorial, you will find special Python functions designed to help the user who wishes to use Read et al. EOSs or any other piecewise polytropic EOS.

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

$$\label{latex_pdf_output}$$