General Relativistic Hydrodynamics on a Moving-mesh I: Static Spacetimes

Multimessenger observations are beginning to address some of the most important unsolved problems in gravity and astrophysics. These include testing strong field general relativity (GR), constraining the nuclear equation of state (EOS), elucidating the origin of the elements, and uncovering plausible formation scenarios for black hole and neutron star binaries. However, this promise is predicated on having detailed theoretical models of multimessenger sources, in particular of compact binary merger events.

Construction of self-consistent theoretical models remains a central problem in the theory of compact object mergers, and their complexity cannot be overstated. There are several reasons for this. First, the fluid equations in dynamical spacetimes are considerably more complex than in the Newtonian limit; both the Einstein field equations and the those of general relativistic (magneto)hydrodynamics (GR(M)HD) must be solved. Second, the equations must be evolved over many orders of magnitude in length and time, spanning from the core of the neutron star (NS) out to the far-field wave zone. Third, this is a multiphysics problem, spanning the range from hyperaccretion physics at a black hole horizon to the core of a neutron star to interstellar space and includes the physics of radiation, neutrino radiation, nuclear reactions, nuclear equations of state, gravitational waves, shocks, and accretion, just to name a few.

Fully three-dimensional numerical (magneto-)hydrodynamic simulations are required to accurately model and understand these systems. In the context of binary neutron star GR(M)HD simulations, for example, a number of numerical codes have been developed to solve the GRMHD equations. Most codes solve the Einstein field equations without approximation either on an adaptive-mesh refined (AMR) grid or (pseudo-)spectrally and solve the equations of GR(M)HD on an AMR grid (e.g., Duez et al. 2005; Anderson et al. 2006; Giacomazzo & Rezzolla 2007; Cerdá-Durán et al. 2008; Yamamoto et al. 2008; Duez et al. 2008; Dionysopoulou et al. 2013; Palenzuela 2013; Mösta et al. 2014; Etienne et al. 2015; Kidder et al. 2017). Others solve the Einstein equations in the conformally-flat approximation and solve the equations of GRHD with smoothed-particle-hydrodynamics (SPH) methods (Oechslin et al., 2002).

These two methodologies for solving the (M)HD equations—smooth particle hydrodynamics (SPH) and grid-based solvers—have their respective advantages and disadvantages¹¹1As this paper focuses on a new GRHD solver for static spacetimes, we do not contrast the dynamical GR solvers. However, future plans do include implementation of an Einstein field equation solver (without approximation) based on SENR/NRPy+, as described in the Conclusions and Future Work section.. Smooth particle hydrodynamics (SPH) is based upon the Lagrangian view of the Euler equations where the sampling of a fluid is determined from a finite number of particles, and fluid quantities like density and pressure are determined by computing a smoothing kernel over a number of neighbors. The Lagrangian nature of SPH allows it to conserve linear and angular momentum, but comes at the expense of comparatively poor resolution of shocks due to its smoothing nature. On the other hand, grid based methods have superior shock capturing abilities due to the use of Godonov schemes, but suffer from grid effects, e.g., the presence of grid direction can affect the conservation of angular momentum.

Springel (2010) developed an arbitrary Lagrangian-Eulerian (ALE)/moving-mesh (MM) scheme in an effort to capture the best characteristics of both approaches. This scheme, which is implemented in AREPO (Springel, 2010; Weinberger et al., 2019), relies on a Voronoi tessellation to generate well-defined and unique meshes for an arbitrary distribution of points that deform continuously under the movement of the mesh generating points. Springel (2010) has argued that the use of ALE schemes is important in maintaining the Galilean invariance of Eulerian schemes in the Newtonian case. It has also been argued that these schemes are superior to SPH and Eulerian grid schemes at capturing boundary layer instabilities such as Kelvin-Helmholtz instabilities (Springel 2010 but also see Lecoanet et al. 2016). In any case they do seem ideal for modeling colliding galaxies or stars. In particular, AREPO, has been used to study a number of different astrophysical problems including cosmological galaxy formation (see for instance Vogelsberger et al., 2014), disks, and stellar mergers (Zhu et al., 2015; Ohlmann et al., 2016).

Aside from AREPO, a number of MM codes have been developed based on this scheme. These include TESS (Duffell & MacFadyen, 2011), FVMHD3D (Gaburov et al., 2012), ShadowFax (Vandenbroucke & De Rijcke, 2016), RICH (Yalinewich et al., 2015), DISCO (Duffell, 2016), and MANGA (Chang et al., 2017). These MM schemes have also been extended to include magnetic fields (Pakmor et al. 2011; Mocz et al. 2014; Mocz et al. 2016), higher-order convergence (Pakmor et al., 2016a; Mocz et al., 2015), and new physics, such as cosmic rays (Pakmor et al., 2016b; Pfrommer et al., 2017). In addition, the general scheme of determining the geometry from an arbitrary collection of points has also led to derivative methods such as GIZMO (Hopkins, 2015).

We have recently developed MANGA, a MM hydrodynamic solver for ChaNGa (Chang et al., 2017), which is largely based on the Springel (2010) scheme. MANGA is geared toward the study of dynamical stellar problems, such as common envelope evolution (Prust & Chang, 2019; Prust, 2020) and tidal disruption events (Spaulding & Chang, in preparation). We have also been steadily adding new physics including radiation hydrodynamics (Chang et al., 2020), magnetic fields, and various equations of state (EOSs) such as the HELMHOLTZ EOS (Timmes & Swesty, 2000), the MESA EOS (Paxton et al., 2011; Paxton et al., 2013, 2015), and a nuclear EOS (O’Connor & Ott, 2010; Schneider et al., 2017). In addition, we have also developed moving, reactive boundary conditions (Prust, 2020) and a multi-stepping scheme (Prust & Chang, 2019).

The algorithmic advantages of the ALE scheme in its ability to capture both shocks and interface instabilities make it exceptionally well-suited for an application to GR(M)HD for the problem of compact object mergers. In this paper, we describe ongoing work toward this goal. Here, we describe our extensions to MANGA which enables it to solve the equations of GRHD on static spacetimes. Future work will tackle the problem of dynamical spacetimes. We note that this work is similar to Ryan & MacFadyen (2017), who implemented GRHD in DISCO. The difference is that the work described here is for a completely arbitary unstructured moving mesh.

This paper is organized as follows. We write the equations of GRHD in a flux-conservative form that can be solved on a moving-mesh and pictorially describe such a scheme in §2. We provide an overview of the algorithmic steps in §3 and highlight specific technical details such as the conservative-to-primitive solver, the time integrations, the Riemann solver, the metric map, and the variable validation. We describe a number of test problems in §4 including a Tolman-Oppenheimer-Volkoff (TOV) (Tolman, 1934; Oppenheimer & Volkoff, 1939) star and stellar oscillations. We conclude in §5 and close with a discussion of the road forward.

In the following equations we adopt $G=c=1$ units and the Einstein summation convention. For expressions involving tensors, Latin indices denote the spatial components and Greek indices space and time components.

The equation of hydrodynamics in arbitrary spacetimes can be written in conservative form (see for instance Duez et al. 2005, who adopt the same formulation and variable conventions) by introducing a state vector $\boldsymbol{\mathcal{U}}=(\rho_{*},\tilde{\boldsymbol{S}},\tilde{\tau})$:

where $\rho_{*}=\alpha\rho\sqrt{\gamma}u^{0}$, $\tilde{\boldsymbol{S}}=\rho_{*}h\boldsymbol{u}$, $\tilde{\tau}=\alpha^{2}\sqrt{\gamma}T^{00}-\rho_{*}$. The flux, $\boldsymbol{\mathcal{F}}$ is given by

where $v^{j}$ are the components of the 3-velocity ($\boldsymbol{v}$ is its vector form). The source, $\boldsymbol{\mathcal{S}}$, is

where

The stress energy tensor for a perfect fluid is

where $h=1+\epsilon+P/\rho$ is the specific enthalpy and $u^{\alpha}$ are the respective components of the four velocity.

The associated set of primitive variables are $(\rho,\boldsymbol{v},\epsilon)$, which are the rest mass density, fluid 3-velocity $u^{i}/u^{0}$, and internal energy (measured in the rest frame). We close the set of equations with a simple $\Gamma$-law EOS: $P=(\Gamma-1)\rho\epsilon$. In this work we pick the adiabatic index, $\Gamma$, to be equal to the polytropic index that we select for the neutron star ($\Gamma=2$) discussed in § 4. In doing so, we set the initial internal energy, $\epsilon$, using the polytopic EOS.

Springel (2010) showed that any generic flux-conservative equation (1) can be solved using a finite volume strategy on a moving unstructured mesh. For instance, the MHD equations can also be cast in this form (Pakmor et al., 2011; Duffell & MacFadyen, 2011; Gaburov et al., 2012; Mocz et al., 2014; Mocz et al., 2016). We refer the interested reader to Chang et al. (2017), Prust & Chang (2019), and Chang et al. (2020) for a more detailed discussion of the scheme in MANGA. Here, we provide a summary:

For each cell, the integral over the volume of the $i$th cell defines the charge of the $i$th cell, $\boldsymbol{U}_{i}$, to be

where $V_{i}$ is the volume of the cell. We then use Gauss’ theorem to convert the volume integral over the divergence of the flux in equation (1) to a surface integral:

We now take advantage of the fact that the volumes are Voronoi cells with a finite number of neighbors to define a integrated flux

where $\boldsymbol{F}_{ij}$ and $A_{ij}$ are the average flux and area of the common face between cells $i$ and $j$. The discrete time evolution of the charges in the system is given by:

where $\hat{\boldsymbol{F}}^{n+1/2}_{ij}$ is an estimate of the half-timestep flux between the initial, $\boldsymbol{U}_{i}^{n}$, and final states $\boldsymbol{U}^{n+1}_{i}$; $A^{n+1/2}_{ij}$ is the time-averaged area of the face between $i$ and $j$; and $\boldsymbol{S}_{i}^{n+1/2}=\int_{i}\boldsymbol{\mathcal{S}}dV$ is the time-averaged integrated source function. ²²2We note that $\hat{\boldsymbol{F}}^{n+1/2}_{ij}$ given by the Riemann flux solved in the “rest” frame of the face and boosted back into the “lab” frame.

Because the GR(M)HD equations can be written in flux-conservative form as well, they can also be solved on a moving unstructured mesh. Our algorithm is outlined below:

1. 1.

   At the beginning of a timestep, the Voronoi cells are built and the volume integrated charges, $\boldsymbol{U}$, are mapped to the state vector, $\boldsymbol{\mathcal{U}}$.

2. 2.

   A predictor step is applied to obtain the half-timestep state vector. The conservative variables in the state vector are mapped to half-timestep primitive variables, via the conservatives-to-primitives solver.

3. 3.

   The mesh generating point is drifted a half-timestep forward and the Voronoi mesh is rebuilt on this half-timestep.

4. 4.

   The primitive variables on the faces are reconstructed via slope-limited linear interpolation. This is combined with the metric to produce the state vector on the faces.

5. 5.

   The state vector flux across the moving faces is estimated at the half-timestep, using the relativistic version of the Harten-Lax-van Leer (HLL) approximate Riemann solver (Harten et al., 1983).

6. 6.

   The cells are drifted another half timestep.

7. 7.

   All the fluxes are summed and the source terms are included to update the $\boldsymbol{U}$of the cell to the full timestep.

8. 8.

   The updated charges are mapped to a state vector, and the conservative variables in the state vector are mapped to the primitive variables (via the conservatives-to-primitives solver). The physicality of the primitives is checked and prescribed fixes are applied as needed. This marks the end of a timestep; return to the top as needed until the chosen final time is reached.

We can also use the multi-stepping scheme as described in Prust & Chang (2019), which we do by default. Here the key difference is that each cell is associated with the largest timestep possible for that cell from a factor-of-two hierarchy. Each face is then integrated on the smallest timestep of two neighbors that define it. The changes to the charge of each cell over the timestep are accumulated and then applied at the end of the cell’s timestep.

Most of these steps are fairly straightforward moving-mesh methods. We discuss below a few important technical details of the algorithm that are significantly different than standard moving-mesh techniques described in the literature, including the conservative-to-primitive variable solver, details of the predictor-corrector time stepping, the solution of the Riemann problem on moving faces in GRHD, the validation of the primitive variables to ensure physically relevant states, and the mapping of the metric to arbitrary points in coordinate space.

The GRHD flux and source terms themselves are most cleanly written in Einstein notation. Expanding these equations in full by hand in, e.g., the C language directly would be both time consuming and error prone. NRPy+ ³³3http://nrpyplus.net (Ruchlin et al., 2018), “Python-based code generation for Numerical Relativity… and Beyond!”, converts equations written in Einstein notation into highly optimized C-code kernels. The GRHD flux and source terms, as well as equations needed for the primitives-to-conservatives variable conversion, were written within the NRPy+ framework as part of this work. The GRHD flux and source terms were validated against the hand-coded implementations within IllinoisGRMHD, and the primitive to/from conservative equations were validated by converting many physically valid sets of primitive variables to their conservative form and back (using the Newton-Raphson-based root-finder of Noble et al. 2006 for the conservatives-to-primitives step, as described below). The detailed Jupyter notebook used to generate needed GRHD equations in MANGA may be found in the NRPy+ github repo⁴⁴4https://github.com/zachetienne/nrpytutorial or be viewed directly via nbviewer⁵⁵5https://nbviewer.jupyter.org/github/zachetienne/nrpytutorial/blob/master/Tutorial-GRHD_Equations-Cartesian-c-code.ipynb.

Before proceeding, we note two important points. Previously, MANGA used only cgs units, but for GR, $G=c=1$ units are superior so the code was extended to support them. Second, although our moving mesh code is limited to periodic boundary conditions, no boundary conditions are applied to the metric (as it is static and known) so there is no “periodic gravity.” We simply ensure that the outer boundaries of our numerical domain are sufficiently far away from regions of interest, so that as outer (periodic) boundary conditions are applied to the hydrodynamic variables, they have little to no effect on the simulation.

In future work, we may as needed adopt radiation inflow and outflow boundary conditions as described in Chang et al. (2020) to the hydrodynamic variables on the boundary, to produce arbitrary boundary conditions. We note that this implementation would be similar to the boundary conditions for a “sphered cube” (Burns et al., 2019).

We validate our GRHD extensions to MANGA with a set of two very challenging code tests in full 3D. In the first test, we evolve TOV initial data (§ 4.1) with a fixed background spacetime (§ 4.2). This is a useful code test, as the exact solution is stationary; thus any dynamics in time are purely a result of numerical errors. We apply this fact to directly measure the rate at which our numerical errors converge to zero with increased numerical resolution. In the second test, we evolve the same initial data but with the initial pressure in the TOV star depleted by 10% (§ 4.3), and compare the oscillations induced with those observed in a trusted code IllinoisGRMHD at very high resolution.

We performed these computations on the Stampede 2 supercomputer at the Texas Advanced Computing Center, the Niagara supercomputer at the University of Toronto (Loken et al., 2010; Ponce et al., 2019), and Thorny Flat HPC System at West Virginia University.

We have reviewed the structure and implementation of a moving-mesh general relativistic hydrodynamics solver for static spacetimes in MANGA. Taking advantage that any flux-conservative equation can be solved on a moving, unstructured Voronoi mesh, we briefly describe the general algorithm to solve the GRHD equations. We then elaborate on a number of technical points including implementation of the van Leer method of second-order time integration, strategy for solving the Riemann problem across moving faces, transformation of primitives to/from conservatives variables, validation and correction of hydrodynamic variables, and mapping a spherically symmetric metric to a Voronoi grid.

We then apply our code to the analytic solution of a TOV star and show that our algorithm integrates the star stably as expected. In particular, we show that at high resolution ($10^{6}$ mesh generating points inside the star) that the star is stable and the central density varies by at most 2% over 24 dynamical times. At lower resolution ($10^{5}$ mesh generating points inside the star), the central density suffers a systematic drift toward lower density similar to that observed by Duez et al. (2005) on Eulerian (Cartesian) grids at low resolution with second order spatial reconstruction. Because our current MM method is limited to second-order reconstruction, we find that our results are consistent with that of Duez et al. (2005). Moreover, we also find that we can mitigate this effect in large part by going to higher resolution. Finally we demonstrate that when evolving the same star, but with its initial pressure depleted by 10%, we recover the same fundamental frequency of ensuing oscillations.

This paper is the first in a series that will describe the eventual implementation of a moving-mesh GRHD code on dynamical spacetimes. While the code can already be used to study flows around black holes and neutron stars, this is not our primary interest. Instead we aim to study the mergers of compact binaries with matter. Toward that end a number of future improvements are planned and will be described in future work. These are listed in order of priority:

1. 1.

   Incorporation of a full dynamical spacetime solver. While there are many applications of a moving-mesh GRHD solver in static spacetimes, our main interest is adapting our moving-mesh GRHD to NS-NS and NS-BH mergers. As such, we are incorporating the dynamical spacetime solver SENR/NRPy+ in MANGA. Here the major challenge is to couple the MANGA and SENR/NRPy+ codes, and get them to pass information between one another. We note that although SENR/NRPy+ is currently geared toward solving the binary BH problem (for which $T^{\mu\nu}=0$), all $T^{\mu\nu}$ source terms were recently added to Einstein’s equations within SENR/NRPy+ and validated⁶⁶6https://nbviewer.jupyter.org/github/zachetienne/nrpytutorial/blob/master/Tutorial-Start_to_Finish-BSSNCurvilinear-Neutron_Star-Hydro_without_Hydro.ipynb; download here: https://github.com/zachetienne/nrpytutorial for the case of a TOV star with a fixed stress-energy tensor but an evolving spacetime metric (the so-called “hydro without hydro” test of numerical relativity; Baumgarte et al. 1999).

2. 2.

   Implementation of piecewise-polytropic and tabulated equations of state. As part of our strategy to create realistic simulations of NS-NS and NS-BH mergers, we plan to implement realistic equations of state into MANGA with the GRHD solver. In this regard, we plan to both implement piecewise-polytropic and tabulated equations of state in MANGA. In fact, a nuclear equation of state has already been implemented (Schneider et al., 2017) in the Newtonian version of MANGA. The framework for implementation of generic equations of state is general and can be used (in the Newtonian case) for both the moving-mesh solver and the SPH solver (if required) of MANGA without modification. Here we plan to modify the framework to support the GRHD solver.

3. 3.

   Implementation of GRMHD. Magnetic fields play an important dynamical role in compact binary mergers with matter. Toward that end, we plan on implementing GRMHD in MANGA. To do this, we will need to modify the equation of hydrodynamics to include magnetic fields and include an evolution equation for the magnetic field. In moving, unstructured meshes, two different schemes have emerged for evolving the magnetic field while maintaining the divergence-free condition, $\partial_{i}B^{i}=0$, where $B$ is the magnetic field. These schemes are the divergence cleaning/diffusion methods, e.g., Dedner scheme (Dedner et al., 2002; Pakmor et al., 2011) and vector potential methods (Mocz et al., 2016). We have implemented an MHD scheme in MANGA based on the vector potential scheme of Mocz et al. (2016) and plan to adapt the scheme to the general relativistic case. As this method is similar to the currently implemented method in IllinoisGRMHD, we do not anticipate major technical challenges.

4. 4.

   Implementation of neutrino physics/radiation. Neutrinos affect the outflows and eventual r-process yields of NS-NS and NS-BH mergers, and can act as energizers of the outflows as well as changing the $Y_{e}$ of the resulting outflowing material to change both the total mass of the outflow and its composition. To capture this physics, we plan to adapt our recent time-dependent radiative transfer algorithm (Chang et al., 2020) to include neutrino radiation. In particular, if we ignore metric effects on the radiation, e.g., straight-line propagation, the adaption of our current radiation methods would require the expansion from one photon species to the number of neutrino species and the addition of an evolution equation and source terms for the electron fraction (see for instance Müller et al., 2010; Skinner et al., 2019).

We acknowledge Charles Gammie, Andrew MacFadyen, Philipp Moesta, Kenta Hotokezaka, and Chris White for useful discussions. We thank Paul Duffell and Francois Foucart for crucial insights on the computation of the Riemann flux on unstructured grids. We also thank Scott Noble for providing a thread-safe version of his conservative-to-primitive solver. We thank the anonymous reviewer for helpful comments and criticisms. PC is supported by the NASA ATP program through NASA grant NNH17ZDA001N-ATP. PC acknowledges the support and hospitality of the Center for Computational Astrophysics at the Flatiron Institute where much of this work was carried out. ZBE gratefully acknowledges the NSF for financial support from awards OIA-1458952 and PHY-1806596; and NASA for financial support from awards ISFM-80NSSC18K0538 and TCAN-80NSSC18K1488. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper, http://www.tacc.utexas.edu. Computations were performed on the Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto. IllinoisGRMHD simulations were completed on the Thorny Flat HPC System at West Virginia University (WVU), which was funded in part by the National Science Foundation (NSF) Major Research Instrumentation Program (MRI) Award #1726534, and WVU. We also use the yt software platform for the analysis of the data and generation of plots in this work (Turk et al., 2011). The Flatiron Institute is supported by the Simons Foundation.