Relativistic magnetohydrodynamics in dynamical spacetimes: Numerical methods and tests

Many problems at the forefront of theoretical astrophysics require the treatment of magnetized fluids in dynamical, strongly curved spacetimes. Such problems include the origin of gamma-ray bursts, magnetic braking of differential rotation in nascent neutron stars arising from stellar core collapse or binary neutron star merger, the formation of jets and magnetized disks around newborn black holes, etc. To model these phenomena, all of which involve both general relativity (GR) and magnetohydrodynamics (MHD), we have developed a GRMHD code capable of evolving MHD fluids in dynamical spacetimes. Our code solves the Einstein-Maxwell-MHD system of coupled equations in axisymmetry and in full $3+1$ dimensions. We evolve the metric by integrating the Baumgarte-Shapiro-Shibata-Nakamura equations, and use a conservative, shock-capturing scheme to evolve the MHD equations. Our code gives accurate results in standard MHD code-test problems, including magnetized shocks and magnetized Bondi flow. To test our code’s ability to evolve the MHD equations in a dynamical spacetime, we study the perturbations of a homogeneous, magnetized fluid excited by a gravitational plane wave, and we find good agreement between the analytic and numerical solutions.

Figure 1

Relative error in the central rest-mass density for uniformly rotating star A. The error is plotted for three axisymmetric runs, adopting equatorial symmetry (with resolutions ${32}^2$, ${64}^2$, and ${128}^2$), and the curves are scaled for second-order convergence. All runs are performed with the standard hydrodynamic scheme (HLL with ${\mathrm{PPM}}^{+}$). Outer boundaries are placed at $7.1M$.

Figure 2

Normalized error in the central rest-mass density for star A with different reconstruction methods. All runs are axisymmetric and equatorially symmetric with a resolution of ${64}^2$ and outer boundaries at $7.1M$. With PPM, the drift in the central density is strongly reduced with respect to the MC and CENO results. With ${\mathrm{PPM}}^{+}$, the central density drift disappears almost entirely.

Figure 3

Axisymmetric evolution of uniformly rotating stars. Star A (solid lines) is stable, while star B (dashed lines) is unstable to collapse. The upper window shows the central density normalized to its initial value, while the lower gives the central lapse. The solid dot indicates the first appearance of an apparent horizon during the collapse of star B.

Figure 4

Snapshots of the rest density contours and the velocity field $(v^x,v^z)$ in the meridional plane during the axisymmetric collapse of uniformly rotating star B to a Kerr black hole. The contour lines are drawn for ${\rho }_0={10}^{-(0.3j+0.09)}{\rho }_0^{\max }$ for $j=0,1,\dots ,12$, where ${\rho }_0^{\max }$ is the maximum of ${\rho }_0$ at the time of excision. The thick dashed curve marks the excision zone. The thick solid curve is the apparent horizon. We show the system at the time of excision (upper panel), at the time at which the last of the matter falls into the excision region (middle panel), and at a late time just prior to the termination of the integration (lower panel).

Figure 5

Snapshots of the angular velocity profile for an evolution of the differentially rotating star C in axisymmetry (${48}^2$) with outer boundaries at $7.1M$. The profile is well maintained in the bulk of the star for over 15 central rotation periods ($P_{\mathrm{rot}}$).

Figure 6

Fractional error in the central density versus time for evolutions of star C in 2D and 3D. The 3D run was performed in equatorial symmetry with resolution ${96}^2\times 48$ with outer boundaries at $7.1M$. The axisymmetric run has the same spatial resolution (grid cell size) as the 3D case, and thus employs a ${48}^2$ grid with boundaries at $7.1M$. The errors in the central density are of the same order for both runs.

Figure 7

Density profiles for the nonlinear wave tests at $t=t_{\mathrm{final}}$ (see Table ). Symbols denote data from numerical simulations with resolution $\Delta x=0.01$. Solid lines in the upper 6 panels denote the exact solutions [69]. The solid line in the last panel denotes a numerical simulation with higher resolution, $\Delta x={10}^{-3}$.

Figure 8

Velocity profiles ($u^x$) for the nonlinear wave tests at $t=t_{\mathrm{final}}$. Symbols denote data from numerical simulations with resolution $\Delta x=0.01$. Solid lines in the upper 6 panels denote the exact solutions [69]. The solid line in the last panel denotes a numerical simulation with higher resolution, $\Delta x={10}^{-3}$.

Figure 9

Nonlinear Alfvén wave test. Symbols are simulation results with resolution $\Delta x=0.01$ and solid lines are the exact solution. The profiles are shown at time $t=t_{\mathrm{final}}=2.0$. Our computational domain is $x\in (-2,2)$. We only show the region $0.3\le x\le 2.0$ in this graph.

Figure 10

L1 norms of the errors in $u^x$, $u^y$, $B^y$ and $B^z$ for the nonlinear Alfvén wave test at $t=t_{\mathrm{final}}=2.0$. This log-log plot shows that the L1 norms of the errors in $u^x$, $u^y$ and $B^y$ are proportional to $(\Delta x{)}^2$, and are thus second-order convergent. The error in $B^z$ goes as a slightly higher power of $\Delta x$.

Figure 11

Unmagnetized Bondi accretion onto a Schwarzschild black hole. Four different methods are compared: MC reconstruction, CENO reconstruction, and PPM reconstruction, both with and without Kreiss-Oliger (KO) dissipation. Each run was performed on a ${64}^2$ grid using axisymmetry.

Figure 12

$\delta {\rho }_{\star }$ after $100M$ for various initial values of $b^2/{\rho }_0{|}_{r=2M}$. PPM with Kreiss-Oliger dissipation is used in each run. Results from ${64}^2$ and ${128}^2$ grid runs are compared. The $\delta {\rho }_{\star }$ found with the ${128}^2$ grid is multiplied by 4 to allow scaling to be checked for second-order convergence.

Figure 13

Analytic and numerical solutions for the perturbations of a magnetized fluid due to the presence of a gravitational wave (see Sec. 4d1). The thick solid and thin dotted lines represent, respectively, the analytic and numerical solutions, though the two lines are not readily distinguishable in plots (a) and (b). All quantities are evaluated at $z=1/8$ and are normalized as indicated. Time is normalized by the gravitational wave period.

Figure 14

Relative error in the pressure perturbation (evaluated at $z=1/8$) for three different initial data sets. The solid lines give the numerical results while the dotted lines give the fits from Eq. (80). The three initial data sets are derived from the standard case in Eq. (79) by taking $\{{\rho }_0,P_0{\}}^{\mathrm{new}}=\xi \{{\rho }_0,P_0{\}}^{\mathrm{old}}$, $(B_0^i{)}^{\mathrm{new}}=\sqrt{\xi }(B_0^i{)}^{\mathrm{old}}$, and $\{h_{+0},h_{\times 0}{\}}^{\mathrm{new}}=\zeta \{h_{+0},h_{\times 0}{\}}^{\mathrm{old}}$, where $\xi$ and $\zeta$ are constants and “old” refers to the values in Eq. (79). (Effectively, $T_{\mu \nu }^{\mathrm{new}}=\xi T_{\mu \nu }^{\mathrm{old}}$ and $h_0^{\mathrm{new}}=\zeta h_0^{\mathrm{old}}$.) The regime of validity for the analytic solution is approached for small $h_0$ and $|T_{\mu \nu }|/h_0$ [see Eqs. (70, 71)], or equivalently, for decreasing $\zeta$ and $\xi /\zeta$. For the curves shown, these values are (a) $\zeta =3$, $\xi /\zeta =3$, (b) $\zeta =2$, $\xi /\zeta =2$, and (c) $\zeta =1$, $\xi /\zeta =1$. Moving from (a) to (c), we find that the relative error decreases as expected and that the errors are well fit by Eq. (80). Note that the normalization differs in the three cases to reflect the differing values of $P_0$ and $h_0$.

Figure 15

Velocity projected along the direction of the Alfvén mode eigenvector ($\delta v_A\equiv \delta \mathit{v}·{\tilde{\mathit{u}}}_A/|{\tilde{\mathit{u}}}_A|$), for the case discussed in Sec. 4d1 (the general case). (a) Analytic and numerical plots of the projected velocity. (The curves cannot be distinguished on this scale.) (b) Contributions to the analytic solution for $\delta v_A/h_0$. The particular solution (dashed curve) and Alfvén wave (dotted curve) contributions added together give the total perturbation (solid line). All quantities are evaluated at $z=1/8$.

Figure 16

Velocity projected along the direction of the slow magnetosonic mode eigenvector ($\delta v_{m2}\equiv \delta \mathit{v}·{\tilde{\mathit{u}}}_{m2}/|{\tilde{\mathit{u}}}_{m2}|$), for the same case as in Fig. 15. (a) Analytic and numerical plots of the projected velocity. (b) Contributions to the analytic solution for $\delta v_{m2}/h_0$. The fast magnetosonic (dotted curve), slow magnetosonic (short-dashed curve), and particular solution (long-dashed curve) contributions added together give the total perturbation (solid line). All quantities are evaluated at $z=1/8$.