GiRaFFE: An Open-Source General Relativistic Force-Free Electrodynamics Code

TL;DR

GiRaFFE is the first open-source code for simulating force-free electrodynamics in dynamical, relativistic spacetimes, enabling large-scale astrophysical plasma simulations.

Contribution

It introduces GiRaFFE, a novel open-source GRFFE code based on IllinoisGRMHD, compatible with the Einstein Toolkit for scalable, dynamical spacetime simulations.

Findings

Successfully passes flat and curved spacetime tests

Compatible with adaptive-mesh refinement grids

Ready for large-scale astrophysical simulations

Abstract

We present GiRaFFE, the first open-source general relativistic force-free electrodynamics (GRFFE) code for dynamical, numerical-relativity generated spacetimes. GiRaFFE adopts the strategy pioneered by McKinney and modified by Paschalidis and Shapiro to convert a GR magnetohydrodynamic (GRMHD) code into a GRFFE code. In short, GiRaFFE exists as a modification of IllinoisGRMHD, a user-friendly, open-source, dynamical-spacetime GRMHD code. Both GiRaFFE and IllinoisGRMHD leverage the Einstein Toolkit's highly-scalable infrastructure to make possible large-scale simulations of magnetized plasmas in strong, dynamical spacetimes on adaptive-mesh refinement (AMR) grids. We demonstrate that GiRaFFE passes a large suite of both flat and curved-spacetime code tests passed by a number of other state-of-the-art GRFFE codes, and is thus ready for production-scale simulations of GRFFE phenomena of key interest to relativistic astrophysics.

Tables

Table 1. Table 1: Initial setup for 1D tests. The unprimed and primed fields are in the grid and wave frames respectively, evaluated as functions of the coordinate x 𝑥 x . μ 𝜇 \mu is the wave speed relative to the grid frame (in units where c = 1 𝑐 1 c=1 ), and γ μ subscript 𝛾 𝜇 \gamma_{\mu} is the corresponding relativistic Lorentz factor. The fields in the wave frame can be obtained from that in the grid frame via a simple Lorentz boost (see Eqs. 104 and 105 in Ref.  [ 45 ] ). From the A i subscript 𝐴 𝑖 A_{i} and E i subscript 𝐸 𝑖 E_{i} given above, all other quantities v i superscript 𝑣 𝑖 v^{i} , B i superscript 𝐵 𝑖 B^{i} and S ~ i subscript ~ 𝑆 𝑖 \tilde{S}_{i} needed for a GRFFE evolution can be computed from Eqs.  14 , 16 , and 18 respectively.

Test name
(wave speed)	Vector potential	Electric field
Fast wave
($\mu =1$)
Alfvén wave
($\mu =-0.5$)
Degenerate Alfvén wave
($\mu =0.5$)
Three waves
(stationary Alfvén: $\mu =0$,
fast right-going:
$\mu =1$,
fast left-going:
$\mu =-1$)
FFE breakdown
(none)

Table 2. Table 2: Grid setup for 1D tests.

Test name	$(N_x,N_y,N_z)$	$[x_{min},x_{max}]$	$[y_{min},y_{max}]$ & $[z_{min},z_{max}]$	CFL factor
Fast wave	$\begin{array}{cc}(1280,32,32)& \end{array}$	$[-4.0,4.0]$	$[-0.025,0.025]$
Alfvén wave
Degen. Alfvén
Three waves				$0.5$
FFE breakdown	$\begin{array}{cc}(200,8,8)& \end{array}$	$[-0.4,0.6]$	$\begin{array}{cc}[-0.02,0.02]& \end{array}$

Table 3. Table 3: Initial setup for 3D BH spacetime tests. a 𝑎 a is the dimensionless spin parameter for the BH with M 𝑀 M representing the BH mass. C 𝐶 C and C 0 subscript 𝐶 0 C_{0} are constants (both set to 1 in our simulations), and L ​ ( x ) 𝐿 𝑥 L(x) is the dilogarithm function. To avoid the logarithmic divergence as r → ∞ → 𝑟 r\to\infty , we neglect f ​ ( r ) 𝑓 𝑟 f(r) and f ′ ​ ( r ) superscript 𝑓 ′ 𝑟 f^{\prime}(r) terms in our tests, setting them to zero. From the A i subscript 𝐴 𝑖 A_{i} and E i subscript 𝐸 𝑖 E_{i} given above, all other quantities v i superscript 𝑣 𝑖 v^{i} , B i superscript 𝐵 𝑖 B^{i} and S ~ i subscript ~ 𝑆 𝑖 \tilde{S}_{i} needed for a GRFFE evolution can be computed from Eqs.  14 , 16 , and 18 respectively.

Test name	Vector potential	Electric field
Split Monopole
Exact Wald	$A_{\varphi }=\frac{C_0}{2}{r}^2{\sin }^2\theta .$
Magnetospheric Wald	$A_i=\frac{C_0}{2}\left(g_{i\varphi }+2a{g}_{ti}\right).$	$𝐄=0.$

Table 4. Table 4: Grid setup for 3D tests. n = 1 𝑛 1 n=1 and Δ ​ x min Δ subscript 𝑥 min \Delta x_{\rm min} represent the coarsest refinement level, and the grid spacing of the finest level respectively. L, M, H denote low, medium, and high resolution respectively. The BH mass M 𝑀 M is set to 1 in all tests for numerical convenience (all results are scale-free and can be trivially rescaled to any desired mass), and R NS = 1 subscript 𝑅 NS 1 R_{\rm NS}=1 for the Aligned Rotator test, also for numerical convenience due to the case being scale-free according to R NS subscript 𝑅 NS R_{\rm NS} .

Test name		CFL factor
Split Monopole
Exact Wald
Magnetospheric Wald
Aligned Rotator

Table 5. Table 5: L 2 subscript 𝐿 2 L_{2} norms of the difference between our numerical results and the initial data (Eq.  28 ) at t = 5 ​ M 𝑡 5 𝑀 t=5M for the 4-vector potential (top), and the velocity and magnetic field (bottom) at 3 resolutions. L, M, and H represent choice of numerical resolution, as defined in Table  4 .

Res.	$\Delta A^x\cdot {10}^4$	$\Delta A^y\cdot {10}^4$	$\Delta A^z\cdot {10}^5$	$\Delta \left({\psi }^6\Phi \right)\cdot {10}^4$
H	$4.855$	$4.857$	$4.394$	$1.053$
M	$3.9\cdot$H	$3.9\cdot$H	$4.3\cdot$H	$4.0\cdot$H
L	$15.3\cdot$H	$15.3\cdot$H	$19.3\cdot$H	$15.6\cdot$H

Equations

$$d{s}^2=-{\alpha }^{2}d{t}^2+{\gamma }_{ij}(d{x}^i+{\beta }^{i}dt)(d{x}^j+{\beta }^{j}dt),$$

$$T_{\mathrm{EM}}^{\mu \nu }=F_{\sigma }^{\mu }{F}^{\nu \sigma }-\frac{1}{4}{g}^{\mu \nu }{F}_{\sigma \eta }{F}^{\sigma \eta }.$$

$$F^{\mu \nu }=n^{\mu }{E}^{\nu }-n^{\nu }{E}^{\mu }-{ϵ}^{\mu \nu \sigma \eta }{B}_{\sigma }{n}_{\eta },$$

$$E^{\nu }$$

$$B^{\nu }$$

$$F^{\mu \nu }{u}_{\nu }=-\sqrt{4\pi }{E}_{(u)}^{\mu }=0.$$

$$F^{\mu \nu }{F}_{\mu \nu }=2(B^2-E^2)$$

$${}^{\ast }{F}^{\mu \nu }{F}_{\mu \nu }=4E_{\mu }{B}^{\mu }$$

$${\nabla }_{\nu }{T}^{\mu \nu }\approx {\nabla }_{\nu }{T}_{\mathrm{EM}}^{\mu \nu }=-F^{\mu \nu }{J}_{\nu }=0$$

$$\to \rho E^i+n_{\nu }{ϵ}^{\nu ijk}{J}_j^{(3)}{B}_k=0.$$

$$B_i{E}^i=0,\mathrm{and}$$

$$B^2>E^2.$$

$${\partial }_t{\tilde{S}}_i+{\partial }_j\left(\alpha \sqrt{\gamma }{T}_{\mathrm{EM}}^j{}_i\right)=\frac{1}{2}\alpha \sqrt{\gamma }{T}_{\mathrm{EM}}^{\mu \nu }{\partial }_i{g}_{\mu \nu },$$

$$S_{\mu }=-n_{\nu }{T}_{\mathrm{EM}\mu }^{\nu }.$$

$${\partial }_t{\tilde{B}}^i+{\partial }_j\left(v^j{\tilde{B}}^i-v^i{\tilde{B}}^j\right)=0,$$

$$v^i=4\pi \alpha \frac{{\gamma }^{ij}{\tilde{S}}_j}{\sqrt{\gamma }{B}^2}-{\beta }^i.$$

$${\partial }_t{A}_i={ϵ}_{ijk}{v}^j{B}^k-{\partial }_i(\alpha \Phi -{\beta }^j{A}_j),$$

$$B^i={ϵ}^{ijk}{\partial }_j{A}_k=\frac{[ijk]}{\sqrt{\gamma }}{\partial }_j{A}_k.$$

$${\partial }_t\left[\sqrt{\gamma }\Phi \right]+{\partial }_j\left(\alpha \sqrt{\gamma }{A}^j-{\beta }^j\left[\sqrt{\gamma }\Phi \right]\right)=-\xi \alpha \left[\sqrt{\gamma }\Phi \right].$$

$$v^i=\frac{u^i}{u^0}=\alpha {\bar{v}}^i-{\beta }^i.$$

$${\tilde{S}}_i=\frac{{\gamma }_{ij}(v^j+{\beta }^j)\sqrt{\gamma }{B}^2}{4\pi \alpha }.$$

$${\tilde{S}}_i\to {\tilde{S}}_i-\frac{({\tilde{S}}_j{\tilde{B}}^j){\tilde{B}}_i}{{\tilde{B}}^2}$$

$$v^z=-\frac{{\gamma }_{xz}{v}^x+{\gamma }_{yz}{v}^y}{{\gamma }_{zz}}$$

$$\mu =\frac{B_z{E}_y-B_y{E}_z}{B^2},$$

$$A_{\varphi }=\frac{\mu {\varpi }^2}{r^3};A_r=A_{\theta }=0,$$

$$\mathbf{v}=\Omega {\mathbf{e}}_z\times \mathbf{r},$$

$$B^2-E^2=B_0^2\left(1-\frac{2M{\sin }^2(\theta )}{r}\right).$$

$$\Delta Q=\sqrt{{\int }_{\mathcal{V}}{\left(Q-Q_0\right)}^2{d}^{3}x}.$$

$$\alpha$$

$${\beta }^r$$

$${\beta }^{\theta }$$

$${\gamma }_{rr}$$

$${\gamma }_{r\varphi }$$

$${\gamma }_{\theta \theta }$$

$${\gamma }_{\varphi \varphi }$$

Full text

1 Introduction

With LIGO’s first detections of gravitational waves [1, 2], the era of gravitational wave (GW) astronomy has arrived, and an exciting new window on the Universe has been opened. The prospect of multimessenger observations stemming from an observed gravitational wave source with coincident electromagnetic (EM) and/or neutrino signals promises to provide far deeper insights into the extremely violent processes that generate these GW signals.

0.4 seconds after the first GW detection, a weak gamma-ray burst (GRB) lasting 1 second was observed in the hard X-ray band ($>50$keV) by the Fermi Gamma-ray Burst Monitor, in a region of the sky overlapping the very large LIGO localization area [14]. Stellar-mass black hole binaries (BHBs) are not generally expected to exist in matter-rich environments assumed necessary to power such a GRB, suggesting that the coincidence is most likely due to a statistical/systematic fluke in the Fermi data, but would indicate the unexpected presence of a tenuous plasma surrounding the BHB if the coincidence proved genuine [6]. Regardless, if Nature does provide an electromagnetic (EM) counterpart to the GW emission from BHBs—whether from a stellar-mass binary observable by Advanced LIGO, intermediate-mass or massive BHBs observable by the Laser Interferometer Space Antenna (LISA), or supermassive BHBs observable by pulsar timing arrays—it will most likely be driven by the interaction of the BHB with a tenuous plasma.

BHBs are not the only systems in which strong-field gravitation could affect the dynamics of a tenuous, magnetically-dominated plasma. Neutron star (NS) and pulsar magnetospheres are other prime examples in which magnetic field lines originating from the intense gravitational fields at the stellar surface continue into the magnetosphere just outside the NS. So just as in the case of BHBs, fully self-consistent modeling of these systems requires that the plasma dynamics be followed in the context of nonperturbative solutions to the general relativistic field equations.

Future gravitational wave detections by LIGO may originate from binary NSs or black hole–neutron star (BHNS) binary mergers. Observing an electromagnetic counterpart to these incredibly violent processes may provide deep insights not only into the behavior of matter at extreme densities, but also the mechanism behind short GRBs. The behavior of binary NS magnetospheres during their late-inspiral phase may launch EM counterparts prior to and during a GW detection. In the case of BHNS binaries, the motion of the BH through the magnetic field of its NS companion during the inspiral can form a unipolar inductor, beaming radiation over the full sky each orbit [18, 23, 35, 43]. At the end of the inspiral, either for a BHNS with a sufficiently light, highly-spinning BH or a BNS with sufficiently high mass, NS disruption and subsequent accretion onto the BH (formed from hypermassive NS collapse in the case of BNS) can launch a classic short GRB jet [7, 36, 44, 42, 48].

In all of the aforementioned systems, the dynamics driving the EM signal involve a strongly magnetized tenuous plasma accelerating within a curved, often highly-dynamical spacetime. Properly modeling such dynamics will require that the appropriate general relativistic versions of the equations governing plasma dynamics be chosen. For regions in which the magnetic-to-gas-pressure ratios are less than or are of order unity, such as the interior of a NS or an accretion disk surrounding a BH or BHB, the plasma dynamics are well-described by the equations of ideal general relativistic magnetohydrodynamics (GRMHD). But as these ratios far exceed unity—as is the case for tenuous plasma outside BHs, BHBs, and NSs—the GRMHD equations become stiff, and solving instead the equations of general relativistic force-free electrodynamics (GRFFE) is far more numerically tractable. The equations of GRFFE are a limiting case of ideal GRMHD, in which the magnetic fields (as opposed to hydrodynamic densities and pressures) entirely drive the plasma dynamics.

While technically it would be best to employ an algorithm that can accommodate arbitrary magnetic-to-gas-pressure ratios, in practical terms, the dynamics of many systems of interest are well approximated by either the GRMHD or GRFFE limits. In exceptional cases such as the boundary between a NS and its magnetosphere, an optimal method would allow for feedback from the ideal GRMHD region into the GRFFE region, and vice-versa. To this end, methods have been developed for feeding information from the ideal GRMHD region of the NS interior into the magnetosphere, which is well-described by the equations of GRFFE (see [31, 43, 45]). These methods effectively impose a dynamical boundary condition for the surface of the magnetosphere (i.e., the surface outside of which the equations of GRFFE are solved). While GRMHD flows can inform regions well-described by the GRFFE equations, the reverse is not yet possible with current GRFFE prescriptions, as the GRFFE equations destroy information about hydrodynamic pressures and densities—two essential quantities that must be specified within the ideal GRMHD framework.

The GRFFE equations and their application in modeling important astrophysical scenarios have a long history in the literature and continue to be an active area of study by many independent groups. Komissarov [26] was one of the first to develop a conservative GRFFE formalism and numerical code in the context of strong-field spacetimes to study the properties of the plasma-filled magnetospheres of black holes [27, 28, 29]. Later, McKinney would introduce a general relativistic force-free (GRFFE) numerical formulation that, through minimal modification of a general relativistic magnetohydrodynamic (GRMHD) code, enables evolutions of the GRFFE equations of motion [32, 33, 34]. Adopting such a prescription, Palenzuela et al [38, 39, 40] studied the interaction of binary and single black holes with ambient magnetic fields within the GRFFE approximation, observing the formation of a dual jet morphology when BHBs orbited within a uniform background magnetic field, and a single jet in the context of a single, spinning black hole—a dramatic confirmation of the Blandford-Znajek mechanism [8] (see also [49]). Motivated in part by these results, as well as a desire to more deeply explore important GRFFE phenomena in this age of multimessenger astrophysics, multiple other independent groups have developed GRFFE codes [3, 10, 11, 37, 41, 46, 47, 54, 49].

One GRFFE code of particular interest to work presented here is that of Paschalidis and Shapiro [45], who extended McKinney’s prescription for GRMHD-to-GRFFE code conversion. In part this extension focused on smoothly matching general relativistic ideal MHD to its force-free limit, so that GRMHD dynamics can inform regions of the physical domain more accurately modeled by a GRFFE code. IllinoisGRMHD is an open-source rewrite of the GRMHD code on which Paschalidis and Shapiro’s GRFFE code is based, and the code we present here, GiRaFFE 111GiRaFFE and both initial data modules (GiRaFFEfood and ShiftedKerrSchild) are open source (under the 2-clause BSD or GNU General Public License, version 2.0 or greater) and can be downloaded from https://bitbucket.org/zach_etienne/wvuthorns/., is based on IllinoisGRMHD. Unlike the code of Paschalidis and Shapiro, GiRaFFE is a pure GRFFE code, not currently designed to match the equations of GRMHD to its force-free limit.

GiRaFFE represents the first open-source GRFFE code designed for dynamical spacetime simulations. Based on IllinoisGRMHD, GiRaFFE employs the prescription devised by McKinney [32] and refined by Paschalidis and Shapiro [45] for converting a GRMHD code into a GRFFE code. Like the code on which it was based, GiRaFFE is fully compatible with the Einstein Toolkit adaptive-mesh refinement (AMR) Cactus/Carpet infrastructure [19, 12], and is therefore highly scalable, to tens of thousands of cores. Also, GiRaFFE is quite compact, existing in only about 1,600 lines of code. For dynamical-spacetime GRFFE simulations, GiRaFFE may be immediately coupled to McLachlan [9], a state-of-the-art, open-source, Kranc-generated [25, 30] spacetime evolution module within the Einstein Toolkit.

The rest of the paper is structured as follows. In Sec. 2, we review our adopted GRFFE formalism and in Sec. 3 how these equations are solved within the GiRaFFE code. Then in Sec. 4, we present a large suite of GRFFE code tests, and finally in Sec. 5 we summarize the paper and present plans for future work.

2 Adopted GRFFE Formalism

All equations are written in geometrized units, such that $G=c=1$, and Einstein notation is chosen for implied summation. Greek indices span all 4 spacetime dimensions, and Latin indices span only the 3 spatial dimensions.

The line element for our spacetime in standard 3+1 form is given by

[TABLE]

where $\alpha dt$ denotes the proper time interval between adjacent spatial hypersurfaces separated by coordinate times $t=t_{0}$ and $t=t_{0}+dt$, $\beta^{i}$ the magnitude of the spatial coordinate shift between adjacent hypersurfaces, and $\gamma_{ij}$ is the three-metric within a given hypersurface at coordinate time $t$.

Our 3+1 GRFFE formalism is written in terms of electric and magnetic fields as measured by an observer co-moving and normal to the spatial hypersurface, with 4-velocity $n^{\mu}=(1/\alpha,-\beta^{i}/\alpha)$. The stress-energy tensor of the electromagnetic field is defined as:

[TABLE]

where $F^{\mu\nu}$ is the electromagnetic (or Faraday) tensor:

[TABLE]

In terms of the Faraday tensor, the electric $E^{\nu}$ and magnetic $B^{\nu}$ fields in this frame are given by

[TABLE]

where ${}^{*}F^{\mu\nu}$ is the dual of the Faraday tensor, and $\epsilon^{\nu\mu\sigma\eta}=[\nu\mu\sigma\eta]/\sqrt{|g|}$ is the rank-4 Levi-Civita tensor, with $[\nu\mu\sigma\eta]$ the regular permutation symbol. In ideal MHD, the electric field vanishes for observers moving parallel to the plasma with 4-velocity $u^{\nu}$:

[TABLE]

I.e., when an observer moves in a direction parallel to the magnetic field lines, the electric field vanishes. This is simply a statement of Ohm’s law for the case of perfect conductivity. One implication of perfect conductivity is that “magnetic field lines remain attached to the fluid elements they connect”—the so-called “frozen-in condition of MHD”. In addition, so long as $F^{\mu\nu}\neq 0$, the ideal MHD condition implies that

[TABLE]

See, e.g., [26, 45] for further discussion.

In addition to the ideal MHD condition $E^{\mu}_{(u)}=0$, force-free electrodynamics also assumes that the plasma dynamics are completely driven by the electromagnetic fields (as opposed to, e.g., hydrodynamic pressure gradients). This implies that the stress-energy of the plasma: $T^{\mu\nu}=T_{\rm matter}^{\mu\nu}+T_{\rm EM}^{\mu\nu}$, is completely dominated by the electromagnetic terms, which yields the conservation equation [39, 41]:

[TABLE]

where $J^{(3)}_{i}$ is the 3-current, and $\rho$ the charge density. This formalism is valid in the tenuous plasma of the stellar magnetosphere, where the rest-mass density is vanishingly small and assumed to be zero. From these, the 4-current can be expressed $J^{\nu}=\rho n^{\nu}+\gamma^{i}_{j}J^{j}=\gamma^{\nu}_{\mu}J^{\mu}$

The left-hand side of Eq. 10 is simply the general relativistic expression for the Lorentz force, indicating that indeed the Lorentz force is zero in GRFFE. Notice that if we assume we are not in electrovacuum ($J^{\mu}\neq 0$), multiplying Eq. 10 by $B_{i}$ yields the familiar $B^{i}E_{i}=0$ constraint of force-free electrodynamics (see [26, 45] for full derivation).

In summary, the force-free constraints can be written

[TABLE]

Under these constraints, the GRFFE evolution equations consist of the Cauchy momentum equation and the induction equation (see [26, 32, 45] for derivation):

The Cauchy momentum equation follows directly from the spatial components of $\nabla_{\mu}T^{\mu\nu}=\nabla_{\mu}T^{\mu\nu}_{\rm EM}=0$ (the time component yields the energy equation, which in GRFFE is redundant). We choose to write the momentum equation in conservative form and in terms of the densitized spatial Poynting flux one-form $\tilde{S}_{i}=\sqrt{\gamma}S_{i}$,

[TABLE]

where $S_{i}$ can be derived from the expression of the Poynting one-form,

[TABLE] 2. 2.

The induction equation in the force-free limit is derived from the spatial components of $\nabla_{\mu}{}^{*}F^{\mu\nu}=0$ (the time components yield the “no-monopoles” constraint), and can be written in terms of the densitized magnetic field $\tilde{B^{i}}=\sqrt{\gamma}B^{i}$ as

[TABLE]

where $v^{j}=u^{j}/u^{0}$. As detailed in Appendix A of[45], the force-free conditions do not uniquely constrain $u^{\mu}$, allowing for the freedom to choose from a one-parameter family. As in [32, 45], we choose $u^{\mu}$ to correspond to the minimum plasma 3-velocity that satisfies $F^{\mu\nu}u_{\nu}=0$. This choice of $v^{j}$ is often referred to as the drift velocity, which can be defined in terms of known variables as

[TABLE]

3 Numerical Algorithms

We briefly review the numerical algorithms employed in GiRaFFE to solve the equations of GRFFE as outlined in Sec. 2.

GiRaFFE fully supports Cartesian adaptive mesh refinement (AMR) grids via the Cactus/Carpet [12] infrastructure within the Einstein Toolkit [19].

As in IllinoisGRMHD, GiRaFFE guarantees that the magnetic fields remain divergenceless to roundoff error even on AMR grids by evolving the vector potential $\mathcal{A}_{\mu}=\Phi n_{\mu}+A_{\mu}$, where $A_{\mu}$ is spatial ($A_{\mu}n^{\mu}=0$), instead of the magnetic fields directly. The vector potential fields exist on a staggered grid (as defined in Table 1 of Ref. [20]) such that our magnetic fields are evolved according to the flux constrained transport (FluxCT) algorithm of Refs. [5, 52].

Our choice of vector potential requires that we solve the vector potential version of the induction equation

[TABLE]

where $\epsilon_{ijk}=[ijk]\sqrt{\gamma}$ is the anti-symmetric Levi-Civita tensor and $\gamma$ is the 3-metric determinant, which in a flat spacetime in Cartesian coordinates reduces to 1. $B^{k}$ in Eq. 17 is computed from the vector potential via

[TABLE]

$\Phi$ is evolved via an additional electromagnetic gauge evolution equation, which was devised specifically to avoid the buildup of numerical errors due to zero-speed characteristic modes [21] on AMR grids. Our electromagnetic gauge is identical to the Lorenz gauge, but with an exponential damping term with damping constant $\xi$ [22]:

[TABLE]

Step 0: Initial data: In addition to 3+1 metric quantities in the Arnowitt-Deser-Misner (ADM) formalism [4], GiRaFFE requires that the “Valencia” 3-velocity $\bar{v}^{i}$ and vector potential one-form $A_{\mu}$ be set initially. Regarding the former, the “Valencia” 3-velocity $\bar{v}^{i}$ is related to the 3-velocity appearing in the induction equation $v^{i}$ via

[TABLE]

As for $A_{\mu}$, for all cases in this paper, we set the evolution variable $\left[\sqrt{\gamma}\Phi\right]=0$ initially, and $A_{i}$ is set based on our initial physical scenario.

After $v^{i}$ and $A_{\mu}$ are set, $B^{i}$ is computed via Eq. 18, and then the evolution variable $\tilde{S}_{i}$ is given by

[TABLE]

Step 1: Evaluation of evolution equations: In tandem with the high-resolution shock-capturing (HRSC) scheme within GiRaFFE, the Runge-Kutta fourth-order (RK4) scheme is chosen to march our evolution variables $A_{i}$ and $\tilde{S}_{i}$ forward in time from their initial values, adopting precisely the same reconstruction and Riemann solver algorithms as in IllinoisGRMHD (see Ref. [20] for more details). In short, $A_{i}$ and $\tilde{S}_{i}$ are evolved forward in time using the Piecewise Parabolic Method (PPM) [13] for reconstruction and a Harten-Lax-van Leer (HLL)-based algorithm [24, 17] for approximately solving the Riemann problem. Meanwhile, spatial derivatives within $\left[\sqrt{\gamma}\Phi\right]$’s evolution equation (Eq. 19) are evaluated via finite difference techniques (as in IllinoisGRMHD).

Step 2: Boundary conditions on $A_{\mu}$: At the end of each RK4 substep, the evolved variables $A_{i}$ and $\tilde{S}_{i}$ have been updated at all points except the outer boundaries. So next the outer boundary conditions on $A_{i}$ and $\left[\sqrt{\gamma}\Phi\right]$ are applied. As no exact outer boundary conditions typically exist for systems of interest to GiRaFFE, we typically take advantage of AMR and push our outer boundary out of causal contact from the physical system of interest. However, to retain good numerical stability, we apply “reasonable” outer boundary conditions. Specifically, values of $A_{i}$ and $\left[\sqrt{\gamma}\Phi\right]$ in the interior grid are linearly extrapolated to the outer boundary.

Step 3: Computing $B^{i}$: $B^{i}$ is next computed from $A_{i}$ via Eq. 18.

Step 4: Applying GRFFE constraints & computing $v^{i}$: Truncation, roundoff, and under-sampling errors will at times push physical quantities into regions that violate the GRFFE constraints. To nudge the variables back into a physically realistic domain, we apply the same strategy as was devised in Ref. [45] to guarantee that the GRFFE constraints remain satisfied:

First, we adjust $\tilde{S}_{i}$ via

[TABLE]

to enforce $B^{i}S_{i}=0$, which as shown by Ref. [45], is equivalent to the GRFFE constraint Eq. 11.

Next, we limit the Lorentz factor of the plasma, typically to be 2,000, by rescaling $\tilde{S}_{i}$ according to Eq. 92 in Ref. [45]. After $\tilde{S}_{i}$ is rescaled the 3-velocity $v^{i}$ is recomputed via Eq. 16.

Finally, errors within our numerical scheme dissipate sharp features, so when current sheets appear, they are quickly and unphysically dissipated. This is unfortunate because current sheets lie at the heart of many GRFFE phenomena. So to remedy the situation, we apply the basic strategy of McKinney [32] (that was also adopted by Paschalidis and Shapiro [45]) and set the velocity perpendicular to the current sheet $v^{\perp}$ to zero. For example, if the current sheet exists on the $z=0$ plane, then $v^{\perp}=v^{z}$, which we set to zero via $n_{i}v^{i}=0$, where $n_{i}=\gamma_{ij}n^{j}$ is a unit normal one-form with $n^{j}=\delta^{jz}$. Specifically, in the case of a current sheet on the $z=0$ plane, we set

[TABLE]

at all gridpoints that lie within $|z|\leq 4\Delta z$ of the current sheet.

At present the code addresses numerical dissipation of current sheets only if they appear on the equatorial plane. For cases in which current sheets appear off of the equatorial plane or spontaneously, [32] suggest the development of algorithms akin to reconnection-capturing strategies [50]. We intend to explore such approaches in future work.

Step 5: Boundary conditions on $v^{i}$: $v^{i}$ is set to zero at a given face of our outermost AMR grid cube unless the velocity is outgoing. Otherwise the value for the velocity is simply copied from the interior grid to the nearest neighbor on a face-by-face basis.

After boundary conditions on $v^{i}$ are updated, all data needed for the next RK4 substep have been generated, so we return to Step 1.