Relativistic radiation hydrodynamics in a reference-metric formulation

We adopt a two-moment formalism, together with a reference-metric approach, to express the equations of relativistic radiation hydrodynamics in a form that is well suited for numerical implementations in curvilinear coordinates. We illustrate the approach by employing a gray opacity in an optically thick medium. As numerical demonstrations we present results for two test problems, namely stationary, slab-symmetric solutions in flat spacetimes, including shocks, and heated Oppenheimer-Snyder collapse to a black hole. For the latter, we carefully analyze the transition from an initial transient to a post-transient phase that is well described by an analytically known diffusion solution. We discuss the properties of the numerical solution when rendered in moving-puncture coordinates.

Figure 1

A continuous flat-spacetime solution, showing the fluid rest-mass density ${\rho }_0$ in the top panel and the radiation energy density $E$ in the bottom panel. The black rectangular grid represents the semianalytical solution for the data in the top row of Table 1, while the colored surface shows the numerical solution at time $t=10.053$, obtained with $N_r=320$ radial and $N_{\theta }=120$ angular grid points, with the grid extending to $r_{\mathrm{out}}=24$. The white lines represent our spherical polar coordinate system, showing every 12-th radial and every 4-th angular grid line.

Figure 2

Convergence test for the continuous solution shown in Fig. 1, except that we have also boosted the solution with a speed $\beta =0.1$ in the $z$-direction for this test. The different lines show differences $N^2\mathrm{\Delta }E$, rescaled assuming second-order convergence, between the numerical and semianalytical solutions. The different lines show interpolations to the $z$-axis, for grids with $N_r=32N$ radial and $N_{\theta }=12N$ grid points, and with the numerical grid extending to $r_{\mathrm{out}}=24$. The top left inset shows the behavior in the vicinity of the center, demonstrating second-order convergence even in the presence of the coordinate singularities at the origin. The bottom right inset shows results for the norm $||\mathrm{\Delta }E||\equiv \int |\mathrm{\Delta }E|dV$, integrated to a radius of $r=12$. The solid line in this inset represents a power-law $(1/N{)}^2$.

Figure 3

Same as Fig. 1, except for a solution featuring a shock discontinuity (see the bottom row in Table 1). For this test, performed with $N_r=256$ radial and $N_{\theta }=48$ angular grid points, and the outer boundary at $r_{\mathrm{out}}=16$, we placed the shock discontinuity at $z=2$ rather than at $z=0$, so that the shock front does not coincide with the symmetry plane of the coordinate system.

Figure 4

A schematic spacetime diagram for Oppenheimer-Snyder collapse. The (red) lines starting out vertically at $\tau =0$ trace the worldlines of selected dust particles, with the thick line marking the surface. The (black) dotted horizontal line shows a surface of constant proper time $\tau$ (where $\tau$ is measured by observers comoving with the dust), while the (black) dashed line sketches a surface of constant coordinate time $t$. Also included are two characteristics originating at the stellar surface and traveling inwards; one for the radiation field, and one for gauge perturbations (see text for details).

Figure 5

A comparison of numerical results and analytical expressions for the rest-mass density ${\rho }_0$ (top panel), the radiation energy density $E$ (middle panel), and the magnitude of the flux $F\equiv (F_a{F}^a{)}^{1/2}$ (bottom panel) for a heated Oppenheimer with $R_0=10M$ (see text for details). The markers represent individual Lagrangian fluid tracers (rather than grid points), while the solid lines represent the analytical solution in the diffusion approximation, computed from the proper times and areal radii recorded by the Lagrangian tracers.

Figure 6

The lapse $\alpha$ and the dust velocity $v^r$ at coordinate time $t=5.52M$ [compare the (green) circled data in Fig. 5]. Regions I, II, and III are labeled as in the schematic spacetime diagram of Fig. 4. Note that $\alpha$ does not depend on $R$, and $v^r=0$, in Region I.

Figure 7

Same as Fig. 6 but for the primitive radiation variables $E$ (top panel) and ${\mathcal{F}}^r$ (bottom panel). The inset shows the small increase in $E$ towards larger radii in Region II.

Figure 8

Same as Fig. 6 but for the conserved radiation variables $\overline{\tau }$ (top panel) and ${\overline{S}}_r$ (bottom panel) in the stellar interior.