Magnetohydrodynamics in full general relativity: Formulation and tests

A new implementation for magnetohydrodynamics (MHD) simulations in full general relativity (involving dynamical spacetimes) is presented. In our implementation, Einstein’s evolution equations are evolved by a Baumgarte-Shapiro-Shibata-Nakamura formalism, MHD equations by a high-resolution central scheme, and induction equation by a constraint transport method. We perform numerical simulations for standard test problems in relativistic MHD, including special relativistic magnetized shocks, general relativistic magnetized Bondi flow in stationary spacetime, and a long-term evolution for self-gravitating system composed of a neutron star and a magnetized disk in full general relativity. In the final test, we illustrate that our implementation can follow winding-up of the magnetic field lines of magnetized and differentially rotating accretion disks around a compact object until saturation, after which magnetically driven wind and angular momentum transport inside the disk turn on.

Figure 1

(a) Propagation of a fast shock. The snapshot at $t=2.5$ is shown. Numerical simulation is performed with $N=100$ and $\Delta x=0.02$. (b) Propagation of a slow shock. The snapshot at $t=2$ is shown. The numerical simulation is performed with $N=400$ and $\Delta x=0.005$. Only one fourth of data points are plotted except for [0.92, 1.08] in which all the data points are plotted. For both cases, the initial discontinuities were located at $x=0$, and the shock fronts move with a constant velocity $\mu$ where $\mu =0.2$ and 0.5 for the fast and slow shocks, respectively.

Figure 2

(a) Propagation of a fast switch-off rarefaction wave. The snapshot at $t=1.0$ is shown. (b) Propagation of a slow switch-on rarefaction wave. The snapshot at $t=2$ is shown. For both cases, the initial discontinuities were located at $x=0$, and the numerical simulations are performed with $N=400$ and $\Delta x=0.005$. Only half of all the data are plotted for both figures.

Figure 3

(a) Propagation of a strong continuous Alfvén wave. The snapshot at $t=2$ is shown. The initial configurations of the analytic solution at $t=2$ are shown by the dashed curves. The numerical simulation is performed with $N=500$ and $\Delta x=0.005$. Only $1/4$ of all the data are plotted. (b) L1 norm of the error for $\rho$ (circles) and $P$ (triangles) as a function of $\Delta x$. These variables should be constant in this problem, and the deviation from the stationary values is due to a numerical error. The dotted line denotes the slope expected for second-order convergence.

Figure 4

(a) Shock-tube problem with the initial discontinuities at $x=0$ normal to the magnetic field. The numerical simulation is performed with $N=500$ and $\Delta x=0.005$. We note that a large bump at $x\sim 0.9$ for $u^x$ is due to a numerical error associated with the limiter ($b=2$) of a weak dissipation. The height of this bump is decreased if we use a more dissipative limiter (e.g., the Minmod Limiter with $b=1$). (b) Shock-tube problem with the initial discontinuities at $x=0$ parallel to the magnetic field. The numerical simulation is performed with $N=500$ and $\Delta x=0.005$. For both cases, the snapshot at $t=1.0$ is shown.

Figure 5

Collision of two flows with opposite directions of the tangential component of the magnetic field. The snapshot at $t=1.22$ is shown. The numerical simulation is performed with $N=400$ and $\Delta x=0.005$.

Figure 6

Snapshot of cylindrical blast explosion [multidimensional test (i)] at $t=0.4$. The contour curves for the density, pressure, and magnetic pressure $P_{\mathrm{mag}}$, as well as the magnetic field lines, are shown. The contour curves for each quantity (denoted by $Q$) are drawn for $Q=Q_{\max }\times {10}^{-0.4\times i}$ $(i=1–8)$ (solid curves) and $Q=Q_{\max }\times {10}^{-0.01}$ (thick solid curves). Here $Q_{\max }$ denotes the maximum value. The results with $\Delta x=0.004$ are presented.

Figure 7

Configuration of various quantities in cylindrical blast explosion along $x$ and $y$ axes at $t=0.4$ with different grid resolutions ($\Delta$ is the grid spacing).

Figure 8

The same as Fig. 6 but for an explosion of a rotating cylinder [multidimensional test (ii)] at $t=0.4$. The grid spacing for the corresponding simulation is 0.0025.

Figure 9

Configuration of various quantities in an explosion of a rotating cylinder at $t=0.4$ with different grid resolutions ($\Delta$ denotes the grid spacing).

Figure 10

The density contour curves and velocity vectors (left panel) and the magnetic field lines (right panel) at $t=15$ and $t=30$ for an axisymmetric jet. The grid spacing for the corresponding simulation is 0.06. The density contour curves are drawn by the same method as in Fig. 6.

Figure 11

L1 norm for the error of density ${\rho }_{*}$ as a function of the grid spacing in units of $M$ for the magnetized Bondi flow. The numerical numbers attached for each curve denote the values of $\widehat{\beta }$ [see Eq. (110) for its definition].

Figure 12

(a) Evolution of the central density of a neutron star and mass and angular momentum of a disk around the neutron star. (b) Evolution of the ADM mass, angular momentum, and averaged violation of the Hamiltonian constraint. The ADM mass and angular momentum are shown in units of their initial values. The solid, dashed, and long-dashed curves show the results with (241, 241), (193, 193), and (161, 161) grid sizes, respectively.

Figure 13

Evolution of $U_{\mathrm{mag}}$ in units of initial internal energy of $U_{\mathrm{disk}0}$ for $R_B=1.3\times {10}^{-5}$ (long-dashed curve), $5\times {10}^{-5}$ (solid curve), and $2\times {10}^{-4}$ (dashed curve). In the smaller panel, we display evolution of $U_{\mathrm{mag}}/U_{\mathrm{disk}0}$ for $R_B=5\times {10}^{-5}$ with three levels of grid resolutions with sizes (301, 301) (solid curve), (241, 241) (dashed curve), and (201, 201) (long-dashed curve). Three curves approximately agree.

Figure 14

Snapshots of the density profile of a neutron and a disk in $x-z$ plane at $t=0,9.78,\mathrm{and}19.40\mathrm{ms}$ for $R_B=2\times {10}^{-4}$. For $\rho >{10}^{10}\mathrm{g}/\mathrm{c}{\mathrm{m}}^3$ (and for atmosphere), the density is denoted by the same dark color, and for ${10}^{10}\mathrm{g}/\mathrm{c}{\mathrm{m}}^3\ge \rho \ge {10}^7\mathrm{g}/\mathrm{c}{\mathrm{m}}^3$, the color is changed (from dark to light).

Figure 15

Evolution of $M_{\mathrm{disk}}$ and $J_{\mathrm{disk}}$ for $R_B=0$ (dotted curves), $1.3\times {10}^{-5}$ (long-dashed curves), $5\times {10}^{-5}$ (solid curves), and $2\times {10}^{-4}$ (dashed curves).