EXPLICIT–IMPLICIT SCHEME FOR RELATIVISTIC RADIATION HYDRODYNAMICS

Relativistic flows appear in many high-energy astrophysical phenomena, such as jets from microquasars and active galactic nuclei, pulsar winds, magnetar flares, core-collapse supernovae, and gamma-ray bursts (GRBs). In many of these systems, the magnetic field has a crucial role in dynamics. For example, magnetic fields connecting an accretion disk with a central star, or different points of accretion disks, are twisted and amplified due to the differential rotation, launching the jets (Lovelace 1976; Blandford & Payne 1982; Uchida & Shibata 1985). The production of astronomical jets has been observed in magnetohydrodynamic (MHD) simulations (Hayashi et al. 1996; Meier et al. 2001; Kato et al. 2004). The jets are also powered by a subtraction of the rotational energy of a central black hole via the magnetic field (Blandford & Znajek 1977; Koide et al. 2002; McKinney & Blandford 2009). Also the magnetic fields play an important role in accretion disks in transporting the angular momentum outward, leading to the mass accretion (Balbus & Hawley 1991; Hawley et al. 1995).

In addition to the observational point of view, the radiation field is also an important ingredient in the dynamics of relativistic phenomena. The radiation pressure force plays a key role in jet acceleration (Bisnovatyi-Kogan & Blinnikov 1977; Sikora et al. 1996; Ohsuga et al. 2005; Okuda et al. 2009). Recently, Takeuchi et al. (2010) showed a formation of radiatively accelerated and magnetically collimated jets using non-relativistic radiation MHD (RMHD) simulations (see also Ohsuga et al. 2009; Ohsuga & Mineshige 2011). The radiative effect also plays important roles in non-relativistic phenomena. For example, the dominance of the radiation pressure in optically thick, geometrically thin (standard) accretion disks could be unstable to thermal and viscous instabilities (Lightman & Eardley 1974; Shibazaki & Hōshi 1975; Shakura & Sunyaev 1976; Takahashi & Masada 2011). Here we note that the global radiation hydrodynamic (RHD) simulations showed the growth of instabilities (Ohsuga 2006), but Hirose et al. (2009) concluded that the radiation-pressure-dominated region is thermally stable using RMHD simulations of a local patch of the disk. Thus the high-energy phenomena would be understood only after we take into account the magnetic field and the radiation consistently.

Many methods have been proposed to include radiative effects in non-relativistic plasmas. Solving a full radiation transfer equation is the most challenging task for directly connecting results between numerical simulations and observations (Stone et al. 1992). But it is a hard task to couple with the hydrodynamical code due to its complexity and high computational costs. Another class of simple method is based on the moment method. In the flux-limited diffusion (FLD) approximation, a zeroth moment equation of the radiation transfer equation is solved to evaluate the radiation energy density, E'_r. The radiative flux, F'_r, is given based on the gradient of the radiation energy density without solving the first moment equation. This is quite a useful technique and gives appropriate radiation fields within the optically thick regime, but we should keep in mind that it does not always give precise radiation fields in the regime where the optical depth is around unity or less (see Ohsuga & Mineshige 2011). Thus, it is better to solve the zeroth and first moment equations.

When we solve two moment equations, the closure relation, P'^ij_r = D^ijE'_r, is used to evaluate the radiation stress, P'^ij_r, where D^ij is the Eddington tensor. In the Eddington approximation, the tensor is simply given by D^ij = δ^ij/3, where δ^ij is the Kronecker delta. This implies that the radiation field is isotropic in the comoving frame, which is a reasonable assumption in the optically thick media. Another class of the Eddington tensor is obtained by approximately taking into account the anisotropy of radiation fields (so-called M-1 closure; Kershaw 1976; Minerbo 1978; Levermore 1984). González et al. (2007) implemented the M-1 closure in their non-relativistic hydrodynamic code and demonstrated that the anisotropic radiation fields around an obstacle are accurately solved.

In the framework of relativistic radiation hydrodynamics (RRMHD), Farris et al. (2008) first proposed a numerical scheme for solving general relativistic radiation magnetohydrodynamic equations assuming an optically thick medium (i.e., D^ij = δ^ij/3). Zanotti et al. (2011) implemented the general relativistic radiation hydrodynamics code and adopted it for the Bondi–Hoyle accretion on to the black hole. Recently, Shibata et al. (2011) proposed another truncated moment formalism of radiation fields in optically thick and thin limits.

In the RR(M)HD, there are three independent timescales. The first one is the dynamical timescale t_dyn ≡ L/v, where L and v are the typical size of the system and the fluid velocity. The second one is the timescale at which the characteristic wave passes the system, t_w ≡ L/λ, where λ is the characteristic wave speed. These two timescales are comparable and close to the light crossing time, L/c (where c is the speed of light) in relativistic phenomena (i.e., v ∼ λ ∼ c). The third one is the heating/cooling and scattering timescales, t_ab and t_sc. These timescales become much shorter than the dynamical one in the optically thick medium, since the gas and the radiation frequently interact. This indicates that equations for the radiation fields become stiff. In the explicit scheme, a numerical time step Δt should be restricted to being shorter than these typical timescales (Δt < min [t_dyn, t_w, t_ab, t_sc]) to ensure the numerical stability. Therefore, high computational costs prevent us from studying long-term evolutions when the gas is optically thick, if we are interested in the phenomena of the dynamical timescale.

In this paper, we propose an explicit–implicit numerical scheme to overcome this problem. As the first step, we consider the relativistic radiation hydrodynamics by neglecting the magnetic field (but, see Section 5 for the discussion how to extend our numerical scheme to RRMHD). The proposed scheme ensures a conservation of total energy and momentum. Governed equations are integrated in time using both explicit and implicit schemes. The former solves an hyperbolic term (timescale ∼t_w ∼ t_dyn). The latter treats the exchange of energy and momentum between the radiation and the matter via the gas–radiation interaction, whose timescales are ∼t_ab and t_sc. This method allows us to take a larger time step Δt > t_ab, t_sc than that of the explicit scheme.

This paper is organized as follows. In Section 2, we introduce argument equations for RRHD and the numerical scheme is shown in Section 3. Numerical results of one- and two-dimensional tests are shown in Section 4. Discussion and summary are appeared in Sections 5 and 6.

In the following, we take the light speed as unity. The special relativistic radiation magnetohydrodynamic equations of ideal gas consist of conservation of mass,

$$\begin{equation} (\rho u^\nu)_{,\nu } = 0, \end{equation} \tag{ 1 }$$

conservation of energy momentum,

$$\begin{equation} \left(T^{\mu \nu }_\mathrm{HD} +T^{\mu \nu }_\mathrm{rad}\right)_{,\nu } =0, \end{equation} \tag{ 2 }$$

and equations of radiation energy momentum

$$\begin{equation} T^{\mu \nu }_\mathrm{rad,\nu } = -G^\mathrm{\mu }, \end{equation} \tag{ 3 }$$

where ρ is the proper mass density, u^μ = γ(1, vⁱ) is the fluid four velocity, and T^μν_HD and T^μν_rad are energy momentum tensors of fluid and radiation. Here, $\gamma = \sqrt{1+u_i u^i}$ is the bulk Lorentz factor and vⁱ is the fluid three velocity. Greek indices range over 0, 1, 2, 3 and Latin ranges over 1, 2, 3, where 0 indicates the time component and 1, 2, 3 do space components.

The energy momentum tensor of fluids is written as

$$\begin{equation} T^{\mu \nu }_\mathrm{HD} = \rho \xi u^\mu u^\nu + p_g \eta ^{\mu \nu }, \end{equation} \tag{ 4 }$$

where p_g is the gas pressure and diagη^μν = (− 1, 1, 1, 1) is the Minkowski metric. The specific enthalpy of relativistic ideal gas ξ is given by

$$\begin{equation} \xi = 1 + \frac{\Gamma }{\Gamma -1}\frac{p_g}{\rho }, \end{equation} \tag{ 5 }$$

where Γ is the specific heat ratio.

The energy momentum tensor of radiation is written as

$$\begin{equation} T^{\mu \nu }_\mathrm{rad} =\left(\begin{array}{@{}cc@{}}E_r, & F_r^j \\ F_r^{i}, & P^{ij}_r \end{array}\right) \end{equation} \tag{ 6 }$$

where E_r, Fⁱ_r, and P^ij_r are the radiation energy density, flux, and stress measured in the laboratory frame.

The radiation exchanges its energy and momentum with fluids through absorption/emission and scattering processes. The radiation four force G^μ is explicitly given by

$$\begin{eqnarray} G^0 &=& {-}\rho \kappa \left(4\pi \mathrm{B} \gamma - \gamma E_r + u_i F_r^i\right)\nonumber \\ && - \rho \sigma _s\left[\gamma u^2 E_r + \gamma u_i u_j P^{ij}_r-(\gamma ^2+u^2) u_i F_r^i\right] \end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray} G^{i} &=& - 4\pi \rho \kappa \mathrm{B} u^i + \rho (\kappa + \sigma _s)\big(\gamma F_r^i-u_jP^{ij}_r\big) \nonumber \\ && -\rho \sigma _s u^i \left(\gamma ^2 E_r - 2\gamma u_j F_r^j + u_j u_k P_r^{jk}\right), \end{eqnarray} \tag{ 8 }$$

where $u = \sqrt{u_i u^i}$. κ and σ_s are absorption and scattering coefficients measured in the comoving frame. Thus, Equation (3) is a mixed-frame radiation energy–momentum equation, with the radiation field defined in the observer frame, while the absorption and scattering coefficients are defined in the comoving frame.

In Equations (7) and (8), we assume the Kirchhoff–Planck relation so that the emissivity η is replaced by the blackbody intensity B as η = κB in the comoving frame.

The blackbody intensity B is a function of gas temperature T as B = a_RT⁴/(4π), where a_R is related to the Stefan–Boltzmann constant σ_SB = a_R/4. The plasma temperature is determined by

$$\begin{equation} p_g = \frac{\rho k_B T}{\mu m_p}, \end{equation} \tag{ 9 }$$

where k_B and m_p are the Boltzmann constant and the proton mass, respectively, and μ is the mean molecular weight.

Since we consider the moment equation of radiation fields, we need another relation between E_r, Fⁱ_r, and P^ij_r to close the system, i.e., the closure relation. Here, we leave it as a general form described by

$$\begin{equation} P^{\prime ij}_r = P^{\prime ij}\big(E_r^{\prime }, F_r^{\prime k}\big), \end{equation} \tag{ 10 }$$

where dash denotes the quantity defined in the comoving frame. The explicit form of closure relation is introduced in Section 3.3 after the formulation.

In the following, we deal with the one-dimensional conservation law in the x-direction:

$$\begin{equation} \frac{\partial \mathcal {U}}{\partial t} + \frac{\partial \mathcal {F}}{\partial x} = \mathcal {S}, \end{equation} \tag{ 11 }$$

which follows Equations (1), (2), and (3). Primitive variables $\mathcal {P}$ are defined as

$$\begin{equation} \mathcal {P} = \left(\begin{array}{c}\rho \\ u^k \\ p_g \\ E_r \\ F_r^k\\ \end{array}\right). \end{equation} \tag{ 12 }$$

Conserved variables $\mathcal {U}$, fluxes $\mathcal {F}$, and source terms $\mathcal {S}$ have forms of

$$\begin{equation} \mathcal {U} = \left(\begin{array}{c}D \\ m_t^k\\ E_t\\ E_r\\ F_r^k\\ \end{array}\right), \end{equation} \tag{ 13 }$$

$$\begin{equation} \mathcal {F} = \left(\begin{array}{c}D v^x \\ \rho \xi u^x u^k + p_g\delta ^{xk} + P_r^{xk} \\ m_t^x \\ F_r^x \\ P_r^{xk} \end{array}\right), \end{equation} \tag{ 14 }$$

and

$$\begin{equation} \mathcal {S} \equiv \left(\begin{array}{c}0 \\ 0\\ 0\\ S_E\\ S_F^k\\ \end{array}\right) = \left(\begin{array}{c}0 \\ 0\\ 0\\ -G^0\\ -G^k\\ \end{array}\right), \end{equation} \tag{ 15 }$$

where D = ργ is the mass density measured in the laboratory frame. The total energy density E_t and momentum density m_t are given by

$$\begin{equation} E_t = E_\mathrm{HD} + E_r = \rho \xi \gamma ^2 - p_g + E_r, \end{equation} \tag{ 16 }$$

and

$$\begin{equation} m^k_t = m^k_\mathrm{HD} + F_r^k = \rho \xi \gamma u^k + F_r^k, \end{equation} \tag{ 17 }$$

where E_HD and m^k_HD denote the energy and momentum density of the fluids, respectively.

It should be noted that E_r and F^k_r are not only primitive variables, but are also conserved variables. Thus, it is straightforward for the radiation fields to convert the primitive variables from the conserved variables, as we will see in Section 3.2.

In this section, we propose a new numerical scheme for solving RRHD equations. A conservative discretization of one-dimensional equation (11) over a time step Δt from t = nΔt to t = (n + 1)Δt with grid spacing Δx is written as

$$\begin{equation} \mathcal {U}^{n+1}_i= \mathcal {U}^n_i - \frac{\Delta t}{\Delta x} \left(f_{i+\frac{1}{2}} - f_{i-\frac{1}{2}}\right) + \mathcal {S}_i\Delta t, \end{equation} \tag{ 18 }$$

where $\mathcal {U}^n_i$ is the conservative variable at x = x_i and t = nΔt, and f_{i ± (1/2)} is the upwind numerical flux at the cell surfaces x = x_{i ± (1/2)}. In the numerical procedure, we divide Equation (18) into two parts, the hyperbolic term

$$\begin{equation} \mathcal {U}^{*}_i = \mathcal {U}^{n}_i - \frac{\Delta t}{\Delta x} \left(f_{i+\frac{1}{2}} - f_{i-\frac{1}{2}}\right), \end{equation} \tag{ 19 }$$

and the source term

$$\begin{equation} \mathcal {U}^{n+1}_i = \mathcal {U}^{*}_i + \mathcal {S}_i \Delta t, \end{equation} \tag{ 20 }$$

where $\mathcal {U}^{*}$ is the conservative variable at the auxiliary step. In the following subsection, we describe how to solve these two equations.

3.1. Hyperbolic Term

For the hyperbolic term, Equation (19) is solved as an initial value problem with the initial condition

$$\begin{equation} \mathcal {U}(t^n, x) = \left\lbrace \begin{array}{@{}lll}\mathcal {U}_L & \mathrm{for} & x<x_{i+\frac{1}{2}}\\ \ms \mathcal {U}_R & \mathrm{for} & x>x_{i+\frac{1}{2}} \end{array}, \right. \end{equation} \tag{ 21 }$$

where the subscript L and R denote left and right constant states on the cell interface. When we take $\mathcal {U}_L = \mathcal {U}_i$ and $\mathcal {U}_R = \mathcal {U}_{i+1}$, the numerical solution is reliable with a first-order accuracy in space.

In order to improve the spatial accuracy, primitive variables at the zone surface $\mathcal {P}_{i\pm 1/2}$ are calculated by interpolating them from a cell center to a cell surface (the so-called reconstruction step). The primitive variables at the left and right states are written as

$$\begin{eqnarray} &&\mathcal {P}_{i\pm \frac{1}{2},S} = \mathcal {P}_i \pm \frac{\delta _x \mathcal {P}}{2}, \end{eqnarray} \tag{ 22 }$$

where we take S = L(R) with the plus (minus) sign. The slope δ_xP should be determined so as to preserve the monotonicity. In our numerical code, we adopt the harmonic mean proposed by van Leer (1977), which has a second-order accuracy in space, as

$$\begin{equation} \delta _x \mathcal {P}= \frac{2 \mathrm{max(0,\Delta \mathcal {P}_+ \Delta \mathcal {P}_-)}}{\Delta \mathcal {P}_+ + \Delta \mathcal {P}_-}, \end{equation} \tag{ 23 }$$

where

$$\begin{equation} \Delta \mathcal {P}_\pm = \pm (\mathcal {P}_{i\pm 1} - \mathcal {P}_i). \end{equation} \tag{ 24 }$$

Here we adopted a second-order accurate scheme, but a higher order scheme such as piecewise constant method (Komissarov 1999), piecewise parabolic method (Colella & Woodward 1984; Martí & Mueller 1996), or other scheme, such as the ENO scheme (Del Zanna & Bucciantini 2002), can be applicable to both the hydro part and the radiation part.

By computing the primitive variables at the cell surface, the flux $\mathcal {F}_{i \pm 1/2,S}$ is calculated directly from $\mathcal {P}_{i\pm ({1}/{2}),S}$ (but see Section 3.3 for the procedure of calculating P^xk_r, which appeared in Equation (14), from P'^ij_r). Then, numerical fluxes are calculated using an approximate Riemann solver. We adopted the HLL method (Harten et al. 1983), which can capture the propagation of fastest waves. The numerical flux is then computed as

$$\begin{eqnarray} &&f_{i+\frac{1}{2}}= \frac{\lambda _+\mathcal {F}_{i+\frac{1}{2},L} - \lambda _-\mathcal {F}_{i+\frac{1}{2},R} +\lambda _+\lambda _- \left(\mathcal {U}_{i+\frac{1}{2},R} - \mathcal {U}_{i+\frac{1}{2},L}\right)}{\lambda _+ - \lambda _-}. \nonumber \\ \end{eqnarray} \tag{ 25 }$$

Here λ₋ and λ₊ are

$$\begin{equation} \lambda _-= \mathrm{min}(0, \lambda _{L-}, \lambda _{R-}) \end{equation} \tag{ 26 }$$

and

$$\begin{equation} \lambda _+= \mathrm{max}(0, \lambda _{L+}, \lambda _{R+}), \end{equation} \tag{ 27 }$$

where λ_{S +} and λ_{S −} are the right and left going wave speed of the fastest mode. For example, they correspond to the sound speeds for the relativistic hydrodynamics. When the radiation field is included, the light mode determines the fastest wave speed, which depends on the closure relation. The fastest wave speed is discussed after specifying the closure relation in Section 3.3. We note that although we adopted the HLL scheme for simplicity, the higher order approximate Riemann solvers such as the HLLC (Mignone & Bodo 2005, 2006; Honkkila & Janhunen 2007) and HLLD (Mignone et al. 2009) can be implemented in the relativistic radiation (magneto)hydrodynamics (Sekora & Stone 2009).

Using the numerical flux $f_{\pm \frac{1}{2}}$, the conserved variables U* is obtained from Equation (19). It should be noted that the total energy density E^{n + 1}_t and momentum density m^{n + 1, k}_t at t = (n + 1)Δt are already solved, although we only integrate the hyperbolic term. Thus, when the radiation energy density and flux are obtained by solving Equation (20), the fluid energy density and momentum density are immediately computed (see Equations (43) and (44)).

3.2. Source Term

Next, we show how to solve Equation (20). The source term appears in Equation (3), which treats the interaction between the radiation and the matter. As discussed in Section 1, the heating/cooling and scattering timescales can be shorter than the dynamical timescale in optically thick media. It prevents us from studying the long time evolution (∼t_dyn) when the equation is explicitly integrated. To overcome this difficulty, we want to construct an implicit scheme as

$$\begin{eqnarray} &&E_r^{n+1} = E_r^{*} + \Delta t \mathcal {S}_E\big(E_r^{n+1}, {F}_r^{n+1,j}, P_r^{n+1,jk}, \mathcal {P}_h^{n+1}\big), {} \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} &&F_r^{n+1,i} = F_r^{*,i} + \Delta t \mathcal {S}_F^i\big(E_r^{n+1}, {F}_r^{n+1,j}, P_r^{n+1,jk}, \mathcal {P}_h^{n+1}\big),{} \end{eqnarray} \tag{ 29 }$$

where $\mathcal {P}_h$ is primitive variables of fluids (i.e., ρ, uⁱ, p_g). We confront two problems to solve Equations (28) and (29). The first problem comes from the appearance of $\mathcal {P}_h^{n+1}$. Since $\mathcal {P}_h^{n+1}$ should be obtained after computing $\mathcal {U}^{n+1}$, Equations (28) and (29) become nonlinear equations for $\mathcal {U}^{n+1}$. Another difficulty comes from the closure relation P^ij_r = D^ijE_r. The Eddington tensor D^ij = D^ij(E_r, Fⁱ_r) is generally a nonlinear function of E_r and Fⁱ_r, so we cannot adopt the simple implicit method by replacing E_r → E^{n + 1}_r and Fⁱ_r → F_r^{n + 1, i} in the Eddington tensor.

In this paper, we propose an iterative method, which solves the following equation:

$$\begin{eqnarray} &&{E}_r^{(m+1)} = E_r^{*} + \Delta t \mathcal {S}_E\big(E_r^{(m+1)}, {F}_r^{(m+1),j}, P_r^{(m+1),jk}, \mathcal {P}_h^{(m)}\big), \nonumber \\ \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} &&{F}_r^{(m+1),i} = F_r^{*,i} + \Delta t \mathcal {S}_F^i\big(E_r^{(m+1)}, {F}_r^{(m+1),j}, P_r^{(m+1),jk}, \mathcal {P}_h^{(m)}\big), \nonumber \\ \end{eqnarray} \tag{ 31 }$$

where m = 0, 1, 2... indicates the iteration step. We take E⁽⁰⁾_r = E_rⁿ, F^{(0), j}_r = F_r^{n, j}, and $\mathcal {P}_h^{(0)}=\mathcal {P}_h^n$ for the initial guess (also E*_r and F*, ^j_r can be candidates for the initial guess). We note that the hydrodynamic quantity is explicitly evaluated at mth step due to the complexity discussed above.

Next we introduce the following two variables:

$$\begin{eqnarray} \delta E_r^{(m+1)}&\equiv & E_r^{(m+1)} - E_r^{(m)}, \hspace{2pt} \end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray} \delta F_r^{(m+1),i}&\equiv & F_r^{(m+1),i} - F_r^{(m),i}. \end{eqnarray} \tag{ 33 }$$

By substituting Equations (32) and (33) into Equations (30) and (31), and taking the first-order Taylor series in δE^{(m + 1)}_r and $\delta \mbox{\boldsymbol $F$}_r^{(m+1)}$, we obtain

$$\begin{eqnarray} &&\mbox{\boldsymbol $C$}^{(m)} \left(\begin{array}{@{}c@{}}\delta E_r^{(m+1)}\\ \ms \delta F_r^{(m+1),j}\\ \end{array}\right)= \left(\begin{array}{@{}c@{}}\Delta t S_E^{(m)} +E_r^* - E_r^{(m)} \\ \ms \Delta t S_F^{(m),i} +F_r^{*,i} - F_r^{(m),i} \end{array}\right),\nonumber \\ \end{eqnarray} \tag{ 34 }$$

where $S_E^{(m)}=S_E[E_r^{(m)}, F_r^{(m),j}, P_r^{(m),jk}, \mathcal {P}^{(m)}]$ and $S_F^{(m),i}= S_F^{i}[E_r^{(m)}, F_r^{(m),j}, P_r^{(m),jk}, \mathcal {P}^{(m)}]$. Here C is the 4 × 4 matrix given by

$$\begin{eqnarray} &&\mbox{\boldsymbol $C$} \equiv \mbox{\boldsymbol $1$} - \Delta t \left(\begin{array}{@{}cc@{}}\frac{\partial S_E}{\partial E_r} +\frac{\partial P^{kl}}{\partial E_r} \frac{\partial S_E}{\partial P^{kl}}, \frac{\partial S_E}{\partial F_r^j} +\frac{\partial P^{kl}}{\partial F_r^j} \frac{\partial S_E}{\partial P^{kl}} \\ \ms \frac{\partial S_F^i}{\partial E_r} +\frac{\partial P^{kl}}{\partial E_r} \frac{\partial S_F^i}{\partial P^{kl}}, \frac{\partial S_F^i}{\partial F_r^j} +\frac{\partial P^{kl}}{\partial F_r^j} \frac{\partial S_F^i}{\partial P^{kl}} \end{array}\right), \end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray} \frac{\partial S_E}{\partial E_r}&=& {-}\kappa \rho \gamma + \sigma _s \rho \gamma u^2, \hspace{56pt} \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} \frac{\partial S_E}{\partial F^j_r}&=& \kappa \rho u^j -\sigma _s \rho (\gamma ^2 +u^2) u^j,\hspace{32pt} \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} \frac{\partial S_E}{\partial P_r^{kl}}&=& \sigma _s \rho \gamma u_k u_l,\hspace{88pt} \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} \frac{\partial S_F^i}{\partial E_r}&=& \sigma _s \rho \gamma ^2 u^i,\hspace{92pt} \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} \frac{\partial S_F^i}{\partial F_r^j}&=& {-}\kappa \rho \gamma \delta ^i_j -\sigma _s \rho \gamma \left(\delta ^i_j + 2 u^i u_j\right), \end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray} \frac{\partial S_F^i}{\partial P_r^{kl}}&=& \kappa \rho \delta ^i_k u_l +\sigma _s\rho \left(\delta ^i_k u_l+u^i u_k u_l\right).\hspace{12pt} \end{eqnarray} \tag{ 41 }$$

Here we use Equations (7), (8), and (15). Also (∂P^ij_r)/(∂E_r) and (∂P^ij_r)/(∂F_r^k) are required to complete matrix elements. These quantities depend on the closure relation given in Equation (10). We do not specify the closure relation here, but the explicit form of these quantities is shown in Section 3.3 and the Appendix.

When P^ij_r is the linear function of E_r and Fⁱ_r (e.g., the Eddington approximation), Equations (30) and (31) reduce to

$$\begin{eqnarray} &&\mbox{\boldsymbol $C$}^{(m)} \left(\begin{array}{@{}c@{}}E_r^{(m+1)}\\ \ms F_r^{(m+1),i}\\ \end{array}\right)= \left[\begin{array}{@{}c@{}}E^n_r + S_E^{(m)}\big(E_r=0, F_r^i=0\big)\Delta t\\ \ms F^{n,i}_r +S_F^{(m),i}\big(E_r=0, F_r^i=0\big)\Delta t \end{array}\right]. \nonumber \\ \end{eqnarray} \tag{ 42 }$$

By inverting the 4 × 4 matrix C in Equation (34) or (42), we obtain the radiation energy and flux at (m + 1)th iteration step from Equations (32) and (33). In general, the matrix inversion is time consuming and sometimes unstable. The matrix C is, however, only a 4 × 4 matrix so that we can invert it analytically. We also tried inverting it using LU-decomposition method. We obtain the inverse matric C⁻¹ stably in both scheme, but the analytical method is faster than the LU-decomposition method. Thus, we decide to use analytical expression of C⁻¹.

Next, we calculate the primitive variables $\mathcal {P}^{(m+1)}$ from updated conservative variables $\mathcal {U}^{(m+1)}$. As pointed in Section 2, E_r and F^k_r are the both conserved and primitive variables. Thus, the recovery step is unnecessary for the radiation field. We need to compute $\mathcal {P}_h^{(m+1)}$ from $\mathcal {U}^{(m+1)}$.

Since the total energy E^{n + 1}_t and the momentum m^{n + 1, k}_t are already determined, the energy and momentum of fluids (E^{(m + 1)}_HD and m^{(m + 1), k}_HD) can be calculated as

$$\begin{eqnarray} &&E_\mathrm{HD}^{(m+1)} = E_t^{n+1} - E_r^{(m+1)}, \hspace{2pt} \end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray} &&m^{(m+1),k}_\mathrm{HD} = m_t^{n+1,k} - F^{(m+1),k}_r. \end{eqnarray} \tag{ 44 }$$

Then, three unknown variables ρ^{(m + 1)}, u^{(m + 1), k}, p^{(m + 1)}_g are computed from Dⁿ, m^{(m + 1), k}_HD, and E^{(m + 1)}_HD. Thus, the recovery step in RRHD is the same as that in relativistic HD. We adopt the recovery method developed by Zenitani et al. (2009) for solving a quartic equation. We briefly show the method. In the following discussion, we drop the superscripts n and (m + 1) for simplicity.

The gas density D, momentum mⁱ_HD, and energy E_HD are related to the primitive variables ρ, uⁱ, p_g as

$$\begin{eqnarray} &&D=\rho \gamma, \hspace{34pt} \end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray} &&m^i_\mathrm{HD} = (\rho + \Gamma _1p_g)\gamma u^i, \end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray} &&E_\mathrm{HD} = (\rho + \Gamma _1 p_g)\gamma ^2, \hspace{2pt} \end{eqnarray} \tag{ 47 }$$

where Γ₁ ≡ Γ/(Γ − 1). From Equation (46), we obtain

$$\begin{equation} p_g = \frac{1}{\gamma \Gamma _1}\left(\frac{|m_\mathrm{HD}|}{\sqrt{\gamma ^2-1}} - D\right), \end{equation} \tag{ 48 }$$

where and $m_\mathrm{HD}=\sqrt{\vphantom{A^A}\smash{\hbox{$\eta _{ij} m_\mathrm{HD}^i m_\mathrm{HD}^j$}}}$. By substituting Equation (48) into (47), we obtain a quartic equation on $u=\sqrt{u_i u^i}$

$$\begin{eqnarray} f(u)& \equiv & \Gamma _1^2\left(E_\mathrm{HD}^2 - m_\mathrm{HD}^2 \right) u^4 - 2 \Gamma _1 m_\mathrm{HD} D u^3\nonumber \\ && + \left[\Gamma _1^2 E_\mathrm{HD}^2 - D^2 - 2 \Gamma _1(\Gamma _1-1)| m_\mathrm{HD}|^2\right]u^2\nonumber \\ && -2(\Gamma _1-1)D| m_\mathrm{HD}| u - (\Gamma _1-1)^2 m_\mathrm{HD}^2=0. \nonumber \\ \end{eqnarray} \tag{ 49 }$$

By solving the above quartic equation, p_g, ρ, and uⁱ are computed from Equations (48), (45), and (46). When we solve Equation (49), we first adopt a Brown method (Nunohiro & Hirano 2003), which gives analytical solutions of a quartic equation. If reasonable solutions are not obtained, we solve Equation (49) using the Newton–Raphson method with an accuracy of f(u) ⩽ 10⁻⁸. It is, however, noted that numerical solution converges without switching to the Newton–Raphson method in the test problems shown in Section 4. Also we note that Zenitani et al. (2009) showed that the Brown method can be applicable to obtain solutions with a wide range of parameters.

When all the primitive variables are recovered, we obtain solutions $\mathcal {P}^{(m+1)}$. These solutions would satisfy Equation (42). But we note that while E_r, Fⁱ_r, and P^ij_r are evaluated at the (m + 1)th time step, we evaluate $\mathcal {P}_h$ at the (m)th time step when we solve Equation (42). The evaluation of $\mathcal {P}_h$ at the (m)th iteration might be problematic when the density or temperature jump is large (e.g., at the shock front). Thus, we again solve Equation (42) using updated primitive variables $\mathcal {P}_h (m+1)$ and check the difference between two successive variables, |δE^{(m + 1)}_r|/E_r^{(m + 1)}, |δF^{(m + 1), i}_r/F_r^{(m + 1), i}|, and $|\mathcal {P}_h^{(m+1)}-\mathcal {P}_h^{(m)}|/|\mathcal {P}^{(m+1)}_h|$. When these quantities are larger than a specified value (typically ∼10⁻⁶), the source term is integrated again using updated primitive variables. This process is continued until two successive variables fall below a specified tolerance (Palenzuela et al. 2009). A similar iteration method is also adopted in relativistic resistive magnetohydrodynamics (Palenzuela et al. 2009). Since $\mathcal {P}_h$ is evaluated at the current iteration step, solutions obtained by solving Equation (42) might not converge when the shock appears. But as we will see in Section 4, we obtain solutions correctly even when the cooling time is much shorter than Δt (see, Section 4.1.1), or a shock appears (see Section 4.2).

Now we summarize our newly developed scheme for solving RRHD equations.

• 1.  
  Calculate primitive variables at the cell surface $\mathcal {P}_{i\pm 1/2, S}$ from $\mathcal {P}_i$.
• 2.  
  Compute the numerical flux f_{i ± 1/2} using approximate Riemann solvers (e.g., HLL scheme).
• 3.  
  Integrate the hyperbolic term from numerical flux f_{i ± 1/2} (in Equation (19)). Then, we obtain the intermediate states of conservative variables $\mathcal {U}^*$.
• 4.  
  Compute matrix elements of C given in Equation (35). Then, calculate E^{(m + 1)}_r and F^{(m + 1), i}_r by inverting the 4 × 4 matrix.
• 5.  
  Calculate primitive variables $\mathcal {P}^{(m+1)}$ from E^{(m + 1)}_HD, m^{(m + 1), k}_HD and D^{n + 1}.
• 6.  
  When the difference between two successive values does not fall below a specified tolerance, repeat from step 4. The updated primitive variables of fluids $\mathcal {P}_h^{(m+1)}$ are used to evaluate matrix C.

3.3. Closure Relation

In the above discussion, we do not specify the closure relation, but take the general form given in Equation (10). When we compute the numerical flux using the HLL scheme, the wave speeds of the fastest mode λ_± are needed, which depend on the closure relation. Another term depending on the closure relation is P^ij_r, which appears in numerical flux (Equation (14)). Also we need to evaluate its derivatives, ∂P^ij/∂E_r and ∂P^ij/∂F^k_r to compute the matrix C. In this subsection, we show how to obtain these quantities by specifying the closure relation.

Many kinds of closure relations are proposed by authors, as discussed in Section 1. As the first step of developing the RRHD code, we hereafter restrict our discussion to the closure relation being provided when we assume the Eddington approximation. Then, the closure relation is described by

$$\begin{equation} P^{\prime ij}_r = \frac{\delta ^{ij}}{3}E^{\prime }_r, \end{equation} \tag{ 50 }$$

which is valid when the radiation is well coupled with the matter so that the radiation field is isotropic in the comoving frame. This corresponds to the Eddington factor being 1/3.

First, we show the wave speed of the fastest mode. By assuming the relation given in Equation (50), the characteristic wave velocity of fastest mode in the comoving frame is $1/\sqrt{3}$, which is equivalent to the sound speed in the relativistic regime. Thus the maximum wave velocity is always $1/\sqrt{3}$ in the RHDs with the Eddington approximation. The wave velocity in the observer frame λ₋ and λ₊ is obtained by boosting $\pm 1/\sqrt{3}$ with the fluid velocity.

Next, we show how to obtain radiation pressure P^ij_r measured in the observer frame. By performing Lorentz transformation on the radiation energy–momentum tensor and combining it with Equation (50), the radiation stress tensor in the observer frame obeys the following equation:

$$\begin{eqnarray} P_r^{ij} &+&\left[{-}\frac{\delta ^{ij}}{3}+\frac{u^i u^j }{(1+\gamma)^2} \right] u_k u_mP_r^{km} \nonumber \\ &+& \frac{1}{1+\gamma } \left(u^i u_k P_r^{jk} + u^j u_k P_r^{ik}\right) = R^{ij}, \end{eqnarray} \tag{ 51 }$$

where

$$\begin{eqnarray} R^{ij}&=& \frac{\delta ^{ij}}{3}\left(\gamma ^2 E_r - 2 \gamma u_k F_r^k \right) -u^i u^j E_r \nonumber \\ && + \big(u^i F_r^j + u^j F_r^i\big) + \frac{2}{1+\gamma }u^i u^j u_k F_r^k, \end{eqnarray} \tag{ 52 }$$

(e.g., Kato et al. 2008). Since the radiation stress P^ij_r is a 3 × 3 symmetric matrix, we need to solve sixth-order linear equations to compute P^ij_r as

$$\begin{equation} \mbox{\boldsymbol $A$}(u)p = r, \end{equation} \tag{ 53 }$$

where p^T = (P¹¹_r, P_r²², P³³_r, P_r¹², P¹³_r, P_r²³), r^T = (R¹¹, R²², R³³, R¹², R¹³, R²³), and $\mbox{\boldsymbol $A$}=\mbox{\boldsymbol $A$}(u)$ is the 6 × 6 matrix. Since A is a function of velocity, the derivatives of P^ij_r are described by

$$\begin{equation} \mbox{\boldsymbol $A$}(u)\frac{\partial p}{\partial E_r} = \frac{\partial r}{\partial E_r},\quad \mbox{\boldsymbol $A$}(u)\frac{\partial p}{\partial F^k_r} = \frac{\partial r}{\partial F^k_r}. \end{equation} \tag{ 54 }$$

Thus, P^ij_r, ∂P^ij_r/∂E_r, and ∂P^ij_r/∂F_r^k are computed by inverting matrix A. Since A is a 6 × 6 matrix, it is difficult to obtain $\mbox{\boldsymbol $A$}^{-1}$ using analytical expression. Thus, we use LU-decomposition method to invert matrix A. Explicit forms of $\mbox{\boldsymbol $A$}, (\partial r)/(\partial E_r)$ and (∂r)/(∂F^k_r) are shown in the Appendix.

In this section, we show results of some numerical tests for our RRHD code. This section consists of numerical tests for the radiation field (Section 4.1) and relativistic radiation hydrodynamics (Section 4.2). We assume that the closure relation is given by Equation (50) so that the wave speed of radiation field is $1/\sqrt{3}$.

4.1. Numerical Tests for Radiation Field

In this section, we show results of numerical tests for solving radiation fields given in Equations (3). We assume that the fluid is static and uniform for simplicity. We recover the light speed as c in this subsection. Equations are solved with second-order accuracy in space.

4.1.1. Radiative Heating and Cooling

This test has been proposed by Turner & Stone (2001). We evaluate the validity for the integration of the source term that appeared in Equation (20), which is implicitly and iteratively integrated in our numerical scheme. For this purpose, we start from a static and one-zone fluid (i.e., number of grid points N_x = 1), which is initially not in the thermal equilibrium with the radiation. From these assumptions, the radiation field obeys the following equation:

$$\begin{equation} \frac{d E_r}{d t} = \rho \kappa \left(4\pi {B}-c E_r\right), \end{equation} \tag{ 55 }$$

which is analytically integrated by assuming ρ, κ, and B are constant, and we obtain

$$\begin{equation} E_r(t) = \frac{4\pi }{c} {B} - \left(\frac{4\pi }{c} {B}-E_0\right)e^{-\rho \kappa c t}, \end{equation} \tag{ 56 }$$

where E₀ is the radiation energy density at t = 0. The radiation field approaches that in the local thermal equilibrium (LTE, E_r = 4πB/c).

The mass density is set to ρ = 0.025 g cm⁻³ and the opacity is κ = 0.04 cm² g⁻¹, so that the corresponding absorption timescale is t_ab ≡ 1/(ρκc) = 3.3 × 10⁻⁸ s. We examine two models of the initial radiation energy density, E_r = 10²E_LTE and 10⁻²E_LTE, where E_LTE ≡ a_RT⁴ = 10¹⁰ erg cm⁻³. The specific heat ratio and the mean molecular weight are Γ = 5/3 and μ = 1.0.

Figure 1 shows the time evolution of E_r. Squares and crosses respectively denote results with the explicit and implicit schemes, respectively, for integrating the source term with the time step of Δt = 0.1t_ab, while solid curves show analytical solutions obtained from Equation (56). The data point for the model with Δt = 0.1t_ab is reduced for plotting. The lower and upper plots correspond to E_r(t = 0) = 10⁻²E_LTE and 10²E_LTE, respectively. We can see that results with both schemes excellently agree with analytical solutions when the time step is smaller than the absorption timescale Δt < t_ab. When Δt = 10t_ab, solutions with the implicit scheme (filled circle) stably reach the thermal equilibrium state, although they deviate from the analytical solutions before being in LTE. We note that the explicit scheme does not converge to the analytical solution when we take such a large time step. We also performed simulations with a much larger time step Δt = 10⁴t_ab and found that solutions converge to the analytical solution (not shown in the figure). We also note that the number of iteration steps used in the implicit scheme is less than or equal 2, even if Δt = 10⁴t_ab. Thus, the implicit scheme has a great advantage when the cooling timescale is much smaller than the dynamical timescale since we can take the numerical time step to be t_ab < Δt < t_dyn. The computational time of the explicit and implicit schemes with an equal number of time steps is t_exp: t_imp = 1: 1.04 with Δt = 0.1t_ab.

Zoom In Zoom Out Reset image size

4.1.2. Radiation Transport

This test has been performed by Turner & Stone (2001). When the radiation is injected into uniform matter, it propagates with the characteristic wave velocity while it exchanges energies with the plasma. As a test of this effect, we assume that the fluid is static and uniform in a simulation box bounded by x = [0, L], where L = 1 cm. The radiation is injected from the boundary at x = 0. In the uniform medium, the opacity is set to κ = 0.04 cm⁻² g. The radiation energy is E_LTE = 10¹⁰ erg cm⁻³ and the thermal energy is determined from the LTE condition. We examined two models of the mass density ρ = 0.25 and $25\ \mathrm{g\ \mathrm{cm}^{-3}}$. The corresponding optical depth is τ = ρκL = 0.01 and 1, respectively. At the boundary x = 0, the radiation is injected with the energy density of E_r = 10¹⁰E_LTE. The free boundary condition is applied at x = L. Other parameters are Γ = 5/3, and μ = 1.0. The Courant–Friedrichs–Lewy (CFL) number is taken to be 0.5.

Figure 2 shows results with τ = 0.01. Thick solid curves denote the radiation energy density, and thin solid curves show the radiation flux at t = 0.5, 1.0, 1.5τ_c from left to right, where τ_c = L/c. Since the radiation energy is injected from x = 0, the wave front propagates from left to right as time goes on. The dashed curves denote the position of the light head $l_c=ct/\sqrt{3}$. Since we assume the Eddington approximation on the closure relation, the wave front propagates with the velocity $c/\sqrt{3}$. We can see that the wave front is sharply captured in our simulation code. In FLD approximations, the radiation field evolves obeying the diffusion equation, so that the wave front has a smooth profile (see Figure 7 in Turner & Stone 2001). In our code, although we apply the Eddington approximation, the first-order moment equation is solved. Then equations of the radiation field have hyperbolic form, so that the wave front can be captured using the HLL scheme.

Zoom In Zoom Out Reset image size

Figure 3 shows results of the model of larger density (τ = 1). The radiation propagates with the velocity $c/\sqrt{3}$. Behind the wave front, the radiation field becomes steady and its energy exponentially decreases with x (equivalently, τ = κρx) due to the absorption. When we assume that the steady state and the radiation energy is much larger than that in LTE (i.e., E_r ≫ 4πB/c), Equation (3) can be solved and we obtain

$$\begin{equation} E_r = E_0 \exp (-\sqrt{3}\rho \kappa x),\quad F_r = \frac{cE_0}{\sqrt{3}} \exp (-\sqrt{3}\rho \kappa x), \end{equation} \tag{ 57 }$$

(Mihalas & Mihalas 1984). Dotted curves in Figures 2 and 3 show solutions obtained from Equation (57). We can see that numerical results excellently recover analytical ones.

Zoom In Zoom Out Reset image size

4.1.3. Radiation Transport in Two Dimensions

This test from Turner & Stone (2001) shows how the radiation field propagates in optically thin media in the x−y plane. The simulation domain is x = [ − L, L] and y = [ − L, L], where L = 1 cm. Numerical grid points are (N_x, N_y) = (400, 400). The mass density of uniform fluid is ρ = 0.25 g cm⁻³. The other quantities are the same as those in Section 4.1.2. At the initial state, a larger radiation energy is given inside $r=\sqrt{x^2+y^2}<0.1\mathrm{ cm}$ and its energy is E_r = 10¹⁰E_LTE. The free (Neumann) boundary condition is applied on each side of the box. The CFL condition is taken as 0.5.

Figure 4 shows contours of radiation energy density E_r at t = 1.5τ_c. The solid, dashed, and dotted curves denote contours at 10%, 0.1%, 10⁻³% of its maximum, respectively. Due to the enhancement of radiation energy around the origin at the initial state, the radiation field propagates in a circle, forming a caldera-shaped profile of the radiation energy density. The wave front propagates with the speed of $c/\sqrt{3}$. Figure 5 shows one-dimensional cut of results at t = 1.5τ_c. Asterisks, squares, and filled circles show results along y = 0, x = 0, and y = x, respectively. We can see that the wave packet has a sharp structure. The radiation energy is mainly accumulated around the wave front, while it is very small inside the wave front (caldera floor). Such the caldera structure is not successfully reproduced when we apply FLD approximation (Turner & Stone 2001) since the FLD is formulated based on the diffusion approximation. Although the wave speed reduces to $c/\sqrt{3}$, solving the first-order moment in radiation field (radiation flux) with the Eddington approximation has an advantage for studying the propagation of the radiation pulse in optically thin media.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the problem of Sections 4.1.2 and 4.1.3, the radiation field passes the boundary x = L at $t=\sqrt{3}\tau _c$. When we adopt the free (Neumann) boundary conditions, most of the radiation energy passes the boundary, while some part of it is reflected and stays in the simulation box. The amplitude of the reflected wave E_ref is smaller than that passing the boundary E_pass, max [E_ref/E_pass] ∼ 0.8% for the one-dimensional test and ∼5% for the two-dimensional test. The waves are mainly reflected at the corner of the simulation box in the two-dimensional simulation (i.e., around [x, y] = [ ± L, ±L]). The simple free boundary condition can be applied in the radiation field with this accuracy.

4.2. Numerical Test for Relativistic Radiation Hydrodynamics

Next, we consider the coupling between the radiation and matter fields. We solve RRHD equations with a second-order accurate scheme.

Test problems for radiative shock shown in this subsection are developed by Farris et al. (2008) who proposed and solved four shock tube test problems. Initial conditions of these problems are listed in Table 1. The initial state has a discontinuity at x = 0 and the system is in LTE. Such an initial discontinuity breaks up generating waves. When the waves pass away from the simulation box, the system approaches the steady state. Then we compare our numerical results with semi-analytical solutions after the system reaches a steady state. In Farris et al. (2008), the initial condition is constructed by boosting semi-analytical solutions, so that the shock moves with an appropriate velocity. In our test, the shock is assumed to rest around x = 0 (shock rest frame). Such a problem can be a more stringent test for our code to maintain stationarity (Zanotti et al. 2011).

Table 1. List of Simulation Runs

Model	κ	Γ	State	ρ	pg	ux	uy	uz	E'r
RHDST1	0.4	$\frac{5}{3}$	L	1.0	3.0 × 10−5	0.015	0	0	1.0 × 10−8
			R	2.4	1.61 × 10−4	6.25 × 10−3	0	0	2.50 × 10−7
RHDST2	0.2	$\frac{5}{3}$	L	1.0	4.0 × 10−3	0.25	0	0	2.0 × 10−5
			R	3.11	0.04512	0.0804	0	0	3.46 × 10−3
RHDST3	0.3	2	L	1.0	60	10	0	0	2
			R	8	2.34 × 103	1.25	0	0	1.13 × 103
RHDST4	0.08	$\frac{5}{3}$	L	1.0	6.0 × 10−3	0.69	0	0	0.18
			R	3.65	3.59 × 10−2	0.189	0	0	1.30

Notes. Parameter sets of numerical tests. Scattering coefficient σ_s is taken to be zero in all models.

Download table as:  ASCIITypeset image

The simulation box is bounded by x = [ − L, L], where L = 20 in the normalized unit, and a number of numerical grid points is N_x = 3200. The free boundary condition is applied at x = ±L. We give the physical quantities at the left and right state of ρ, p, u, and E'_r and flux is assumed to be zero. Following Farris et al. (2008), the Stefan–Boltzmann constant is normalized such that 4πa_R = E'_{r, L}/T⁴_L = E'_{r, R}/T⁴_R, where the subscripts L and R denote the quantities at the left (x ⩽ 0) and right (x > 0).

We note that although we adopt an iteration method to integrate stiff source terms (Equation (34)), solutions in the following tests converge within the relative error of 10⁻⁶ without iterations.

4.2.1. Non-relativistic Shock

In this test, the non-relativistic strong shock exists at x = 0 (model RHDST1). The ratio of the radiation to thermal energy in the upstream (x < 0) is 2.2 × 10⁻⁴. We take the CFL condition to be 0.9. Figure 6 shows profile at t = 5000. The physical quantities shown in this figure are ρ, p_g, v^x, E'_r, F'^x_r from top to bottom. (We note that F'^x_r, which is the radiation flux measured in the comoving frame, is different from F^x defined in Farris et al. 2008 by factor γ). Dots and curves denote for the numerical and semi-analytical solutions, respectively.

Zoom In Zoom Out Reset image size

Since we give an initial condition by a step function at x = 0, waves arisen at the discontinuity propagate in the ±x-direction (mainly in the +x-direction since the upstream is supersonic). After the waves pass the boundary at x = ±L, the system reaches the steady state.

Since the radiation energy is negligible, fluid quantities have jumps around x = 0, similar to a pure hydrodynamical shock. The radiation field, on the other hand, has a smooth profile in which the radiation energy is transported in front of the shock. The radiation energy and flux have no discontinuities. The simulation results are well consistent with analytical solutions.

4.2.2. Mildly Relativistic Shock

In this test, the mildly relativistic strong shock exists at x = 0 (model RHDST2). The ratio of the radiation to thermal energy in the upstream is 3.3 × 10⁻³. We take the CFL condition to be 0.9. Figure 7 shows profiles at t = 5000. The physical quantities shown in this figure are ρ, p_g, v^x, E'_r, and F'^x_r from top to bottom. We can see that the radiation energy density is no longer continuous, but jumps at x = 0. In the non-relativistic limit, the radiation energy density and its flux are continuous, but it is not the case in general. They are no longer conserved variables at the shock in the relativistic flow (Farris et al. 2008). The solutions obtained in our numerical simulation (dots) are qualitatively and quantitatively consistent with semi-analytical solutions (solid curves), respectively.

Zoom In Zoom Out Reset image size

4.2.3. Relativistic Shock

In this test, the highly relativistic strong shock exists at x ≃ 0 (model RHDST3). The upstream Lorentz factor is ∼10 and the ratio of the radiation to thermal energy is 3.3 × 10⁻². We take the CFL condition to be 0.9. Figure 8 shows profile at t = 5000. The physical quantities shown in this figure are ρ, p_g, v^x, E'_r, and F'^x_r from top to bottom.

Zoom In Zoom Out Reset image size

After waves generated at the jump (x = 0) pass the boundary, the system approaches the steady state. In this case, all the physical quantities are continuous.

The solutions obtained in our numerical simulation excellently recover analytical solutions even when the flow speed is highly relativistic. To validate our numerical code, we compute the L₁ norms from

$$\begin{equation} L_1[f] = \Delta x \sum _{i=1}{N_x}|f_i - f_\mathrm{a}(x_i)|, \end{equation} \tag{ 58 }$$

where f_i is the physical quantity at the ith grid and f_a is the semi-analytic solution.

We perform simulations with number of grid points N_x = 400, 800, 1600, 3200. In these tests, we use the semi-analytic solution as the initial condition. Figure 9 shows the L₁ norms of errors in ρ, u^x, p_g, E'_r, and F'^x_r. We can see that all errors converge at second order in Δx.

Zoom In Zoom Out Reset image size

4.2.4. Radiation-dominated Shock

In this test, the mildly relativistic, radiation-dominated shock exists at x ≃ 0 (model RHDST4). The ratio of the radiation to thermal energy in the upstream is 20. We take the CFL condition to be 0.3. Figure 10 shows profile at t = 5000. The physical quantities shown in this figure are ρ, p_g, v^x, E'_r, and F'^x_r from top to bottom.

Zoom In Zoom Out Reset image size

After waves generated at the jump (x = 0) pass the boundary, the system approaches the steady state. In this case, all the physical quantities are continuous and their profiles are very smooth since a precursor generated from the shock strongly affects the plasma.

We again find very good agreement with semi-analytical solutions even when the radiation energy dominates the thermal energy. Figure 11 shows the L₁ norms of errors in ρ, u^x, p_g, E'_r, and F'^x_r. In this test, we adopt the analytical solution as an initial condition following Farris et al. (2008). We can see that all errors converge at second order in Δx.

Zoom In Zoom Out Reset image size

In the current study, we construct the relativistic hydrodynamic simulation code including the radiation field. The magnetic field, which would play a crucial role in the relativistic phenomena, is neglected for simplicity. It would be quiet simple to include the magnetic field using a well-developed numerical scheme. When we extend our radiation hydrodynamic code to the radiation magnetohydrodynamic one, the energy–momentum tensor T^{ν, μ}_HD is replaced by that of the magnetofluids,

$$\begin{equation} T^{\mu \nu }_\mathrm{MHD} = (\rho \xi +b^2)u^\mu u^\nu - b^\mu b^\nu +\left(p_g + \frac{b^2}{2}\right)\eta ^{\mu \nu }, \end{equation} \tag{ 59 }$$

where $b^\nu = \left\lbrace \mbox{\boldsymbol $u$} \cdot \mbox{\boldsymbol $B$}, \left[\mbox{\boldsymbol $B$} + (\mbox{\boldsymbol $u$} \cdot \mbox{\boldsymbol $B$})\mbox{\boldsymbol $u$}\right]/\gamma \right\rbrace$ is the covariant form of the magnetic fields. The magnetic field evolves according to the induction equation for the ideal MHD,

$$\begin{equation} \partial _\nu (u^\nu b^\mu - u^\mu b^\nu)=0. \end{equation} \tag{ 60 }$$

Then the energy–momentum conservation equation and the induction equation are explicitly integrated using the HLL or higher order approximate Riemann solvers (Mignone & Bodo 2005, 2006; Honkkila & Janhunen 2007; Mignone et al. 2009). The source term describing the interaction between the matter and the radiation is integrated by applying our proposed scheme.

Another simplification made in this paper is adopting the Eddington approximation that the radiation field is assumed to be isotropic in the comoving frame. The wave velocity of radiation field has a propagation speed of $c/\sqrt{3}$. This would become an issue when the flow speed is relativistic (v ≃ c). The flow speed, in principle, exceeds the phase speed of light mode so that the radiation energy might accumulate in front of the plasma flow. Also the problem appears when we consider the magnetofluids. The fast magnetosonic wave speed can exceeds the reduced light speed $c/\sqrt{3}$. Another problem is that the radiation does not propagate in a straight line since we assume that the radiation field is isotropic. When the optically thick medium with a finite volume is irradiated from one side, there should be a shadow on the other side (Hayes & Norman 2003). Such a shadow is no longer formed when we utilize the Eddington approximation (González et al. 2007). To overcome these problems, we should admit an anisotropy of radiation flux in a comoving frame (Kershaw 1976; Minerbo 1978; Levermore 1984). In our formulation, we do not specify the closure relation until Section 3.3. Thus we can employ M-1 closure by replacing matrix A(u) without the other modification. The explicit–implicit scheme of the relativistic radiation magnetohydrodynamics with the M-1 closure will be reported in the near future.

We have developed a numerical scheme for solving special relativistic radiation hydrodynamics, which ensures the conservation of total energy and flux. The hyperbolic term is explicitly solved using an approximate Riemann solver, while the source term describing the interaction between the matter and the radiation is implicitly integrated using an iteration method (a similar approach is found in Roedig et al. 2012). The advantage of the implicit scheme is that we can take a numerical time step larger than the absorption and scattering timescale. This allows us to study the long-term evolution of the system (typically, the dynamical timescale).

When integrating the source term, we need to invert matrices C and A in our proposed scheme. We have to note that the rank of these matrices used in the implicit scheme is very small (4 × 4 for C and 6 × 6 for A). This is because the interaction between the radiation and the gas is described by a local nature (source term) when we solve zeroth and first moments for radiation. This means that the numerical code can be easily parallelized and a high parallelization efficiency is expected.

We note that since the matrix C is only 4 × 4 matrix, we can invert it analytically. For the matrix A, it is relatively small matrix but it is difficult to invert it analytically. Thus we decide to use LU-decomposition method. Since the LU-decomposition method can invert the matrix directly without iterations, we obtain A⁻¹ stably.

We find that the wave front of radiation field can be sharply captured using the HLL scheme even when we adopt the simple closure relation P'^ij_r = δ^ij/3 (i.e., the Eddington approximation). When we adopt the FLD approximations, such a sharp wave front cannot be captured due to the diffusion approximations. Thus, solving the first-order moment of radiation has an advantage when we consider the optically thin medium although the wave speed reduces to $c/\sqrt{3}$ and the radiation field is isotropic.

We adopted the iteration scheme to integrate stiff source terms. If we need many counts in iteration scheme, it causes a load imbalance between each core when the code is parallelized. But we have to stress that solutions converged only within two steps even if Δt = 10⁴t_ab for the problem of radiative heating and cooling given in Section 4.1.1. For the other tests, solution converges without iteration. Thus, we can expect that the load imbalance due to the iteration scheme is not severe even in the more realistic problems.

In previous papers (Farris et al. 2008; Zanotti et al. 2011), radiation fields are defined in different frames for primitive and conservative variables (mixed frame). In such a method, the radiation moment equations are simply described. However, this method is not suitable for implicit integration. This is because a Lorentz transformation between two frames is needed. By the transformation, expressions of the radiation fields become quite complex even if the Eddington approximation is adopted. Moreover, the radiation stress tensor, P_r is generally nonlinear function of the radiation energy density and the radiation flux through the Eddington tensor. Thus, extension of explicit method to implicit method is not straightforward for the relativistic radiation hydrodynamics.

In our method, we treat radiation fields in the laboratory frame, only. Such a treatment simplifies implicit integration. Although we need to invert 6 × 6 matrix and 4 × 4 matrix, our new method is easier than the method with the mixed frame. Two matrices of C (which is needed for implicit integration) and A (which is needed to compute P_r) are directly inverted using analytical solution and LU-decomposition method. Since we do not need to use some iteration methods to invert them, our method is stable and comparatively easy. Our method would be quite useful for numerical simulations of relativistic astrophysical phenomena, e.g., black hole accretion disks, relativistic jets, GRBs, and so on, since both the high density region and high velocity region exists, and since the radiation processes play important role for dynamics as well as energetics.

We thank Tomoya Takiwaki for fruitful discussions. Numerical computations were carried out on Cray XT4 at the Center for Computational Astrophysics, CfCA, at the National Astronomical Observatory of Japan, on Fujitsu FX-1 at the JAXA Supercomputer System (JSS) at the Japan Aerospace Exploration Agency (JAXA), and on T2K at the University of Tokyo. This work is supported in part by Ministry of Education, Culture, Sports, Science, and Technology (MEXT) for Research Activity Start-up (H.R.T.) 23840045, and for Young Scientist (K.O.) 20740115 (Y.S.) 21018008, 21105511, 23740160 (T.I.) 22·3369, 23740154. K.T. is supported by the Research Fellowship from the Japan Society for the Promotion of Science (JSPS) for Young Scientists. A part of this research has been funded by MEXT HPCI STRATEGIC PROGRAM.

In this Appendix, we show how to obtain the radiation stress in the observer frame from E_r and Fⁱ_r by assuming the Eddington approximation. P^ij_r is a 3 × 3 symmetric matrix and there are six unknowns. Since Equation (51) is a linear with E_r and Fⁱ_r, we need to solve sixth-order linear equations,

$$\begin{equation} \mbox{\boldsymbol $A$} p = r, \end{equation} \tag{ A1 }$$

where $\mbox{\boldsymbol $A$}=\lbrace a_{ij}\rbrace$, p^T = (P¹¹_r, P_r²², P³³_r, P_r¹², P¹³_r, P_r²³), and r^T = (R¹¹, R²², R³³, R¹², R¹³, R²³). The explicit form of a_ij is written as

$$\begin{eqnarray} a_{11}&=& 1+(u_1)(u_1)\left[-\frac{1}{3}+\frac{(u^1)^2 }{(1+\gamma)^2}\right]+\frac{2u^1u_1}{1+\gamma }, \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} a_{12}&=& (u_2)(u_2)\left[-\frac{1}{3}+\frac{(u^1)^2}{(1+\gamma)^2}\right], \hspace{48pt} \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} a_{13}&=& (u_3)(u_3)\left[-\frac{1}{3}+\frac{(u^1)^2}{(1+\gamma)^2}\right],\hspace{48pt} \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} a_{14}&=& (u_1)(u_2)\left[-\frac{1}{3}+\frac{(u^1)^2}{(1+\gamma)^2}\right]+\frac{2u^1u_2}{1+\gamma },\hspace{14pt} \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} a_{15}&=& (u_1)(u_3)\left[-\frac{1}{3}+\frac{(u^1)^2}{(1+\gamma)^2}\right]+\frac{2u^1u_3}{1+\gamma },\hspace{14pt} \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} a_{16}&=& (u_2)(u_3)\left[-\frac{1}{3}+\frac{(u^1)^2}{(1+\gamma)^2}\right]\hspace{54pt} \end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray} a_{21}&=& (u_1)(u_1)\left[-\frac{1}{3}+\frac{(u^2)^2}{(1+\gamma)^2}\right],\hspace{50pt} \end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray} a_{22}&=& 1+(u_2)(u_2)\left[-\frac{1}{3}+\frac{(u^2)^2}{(1+\gamma)^2}\right] +\frac{2u^2u_2}{1+\gamma }, \end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray} a_{23}&=& (u_3)(u_3)\left[-\frac{1}{3}+\frac{(u^2)^2}{(1+\gamma)^2}\right],\hspace{48pt} \end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray} a_{24}&=& (u_1)(u_2)\left[-\frac{1}{3}+\frac{(u^2)^2}{(1+\gamma)^2}\right] +\frac{2u^2u_1}{1+\gamma }, \hspace{32pt} \end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray} a_{25}&=& (u_1)(u_3)\left[-\frac{1}{3}+\frac{(u^2)^2}{(1+\gamma)^2}\right], \hspace{66pt} \end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray} a_{26}&=& (u_2)(u_3)\left[-\frac{1}{3}+\frac{(u^2)^2}{(1+\gamma)^2}\right] +\frac{2u^2u_3}{1+\gamma }, \hspace{32pt} \end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray} a_{31}&=& (u_1)(u_1)\left[-\frac{1}{3}+\frac{(u^3)^2}{(1+\gamma)^2}\right], \hspace{66pt} \end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray} a_{32}&=& (u_2)(u_2)\left[-\frac{1}{3}+\frac{(u^3)^2}{(1+\gamma)^2}\right], \hspace{68pt} \end{eqnarray} \tag{ A15 }$$

$$\begin{eqnarray} a_{33}&=& 1+(u_3)(u_3)\left[-\frac{1}{3}+\frac{(u^3)^2}{(1+\gamma)^2}\right] +\frac{2u^3u_3}{1+\gamma }, \hspace{16pt} \end{eqnarray} \tag{ A16 }$$

$$\begin{eqnarray} a_{34}&=& (u_1)(u_2)\left[-\frac{1}{3}+\frac{(u^3)^2}{(1+\gamma)^2}\right], \hspace{66pt} \end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray} a_{35}&=& (u_1)(u_3)\left[-\frac{1}{3}+\frac{(u^3)^2}{(1+\gamma)^2}\right] +\frac{2u^3u_1}{1+\gamma }, \hspace{32pt} \end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray} a_{36}&=& (u_2)(u_3)\left[-\frac{1}{3}+\frac{(u^3)^2}{(1+\gamma)^2}\right] +\frac{2u^3u_2}{1+\gamma }, \hspace{30pt} \end{eqnarray} \tag{ A19 }$$

$$\begin{eqnarray} a_{41} &=& \frac{(u^1) (u^2) (u_1) (u_1)}{(1+\gamma)^2} +\frac{(u^2) (u_1)}{1+\gamma }, \hspace{62pt} \end{eqnarray} \tag{ A20 }$$

$$\begin{eqnarray} a_{42} &=& \frac{(u^1) (u^2) (u_2) (u_2)}{(1+\gamma)^2} +\frac{(u^1) (u_2)}{1+\gamma }, \hspace{60pt} \end{eqnarray} \tag{ A21 }$$

$$\begin{eqnarray} a_{43} &=& \frac{(u^1) (u^2) (u_3) (u_3)}{(1+\gamma)^2}, \hspace{104pt} \end{eqnarray} \tag{ A22 }$$

$$\begin{eqnarray} a_{44} &=& 1 + \frac{(u^1) (u^2) (u_1) (u_2)}{(1+\gamma)^2} +\frac{(u^1) (u_1)+(u^2)(u_2)}{1+\gamma },\nonumber \\ \end{eqnarray} \tag{ A23 }$$

$$\begin{eqnarray} a_{45} &=& \frac{(u^1) (u^2) (u_1) (u_3)}{(1+\gamma)^2} +\frac{(u^2) (u_3)}{1+\gamma }, \hspace{60pt} \end{eqnarray} \tag{ A24 }$$

$$\begin{eqnarray} a_{46} &=& \frac{(u^1) (u^2) (u_2) (u_3)}{(1+\gamma)^2} +\frac{(u^1) (u_3)}{1+\gamma }, \hspace{62pt} \end{eqnarray} \tag{ A25 }$$

$$\begin{eqnarray} a_{51} &=& \frac{(u^1) (u^3) (u_1) (u_1)}{(1+\gamma)^2} +\frac{(u^3) (u_1)}{1+\gamma }, \hspace{62pt} \end{eqnarray} \tag{ A26 }$$

$$\begin{eqnarray} a_{52} &=& \frac{(u^1) (u^3) (u_2) (u_2)}{(1+\gamma)^2}, \hspace{104pt} \end{eqnarray} \tag{ A27 }$$

$$\begin{eqnarray} a_{53} &=& \frac{(u^1) (u^3) (u_3) (u_3)}{(1+\gamma)^2} +\frac{(u^1) (u_3)}{1+\gamma }, \hspace{60pt} \end{eqnarray} \tag{ A28 }$$

$$\begin{eqnarray} a_{54} &=& \frac{(u^1) (u^3) (u_1) (u_2)}{(1+\gamma)^2} +\frac{(u^3) (u_2)}{1+\gamma }, \hspace{60pt} \end{eqnarray} \tag{ A29 }$$

$$\begin{eqnarray} a_{55} &=& 1+\frac{(u^1) (u^3) (u_1) (u_3)}{(1+\gamma)^2} +\frac{(u^1) (u_1)+(u^3)(u_3)}{1+\gamma },\nonumber \\ \end{eqnarray} \tag{ A30 }$$

$$\begin{eqnarray} a_{56} &=& \frac{(u^1) (u^3) (u_2) (u_3)}{(1+\gamma)^2} +\frac{(u^1) (u_2)}{1+\gamma }, \hspace{56pt} \end{eqnarray} \tag{ A31 }$$

$$\begin{eqnarray} a_{61} &=& \frac{(u^2) (u^3) (u_1) (u_1)}{(1+\gamma)^2}, \hspace{102pt} \end{eqnarray} \tag{ A32 }$$

$$\begin{eqnarray} a_{62} &=& \frac{(u^2) (u^3) (u_2) (u_2)}{(1+\gamma)^2} +\frac{(u^3) (u_2)}{1+\gamma }, \hspace{58pt} \end{eqnarray} \tag{ A33 }$$

$$\begin{eqnarray} a_{63} &=& \frac{(u^2) (u^3) (u_3) (u_3)}{(1+\gamma)^2} +\frac{(u^2) (u_3)}{1+\gamma }, \hspace{58pt} \end{eqnarray} \tag{ A34 }$$

$$\begin{eqnarray} a_{64} &=& \frac{(u^2) (u^3) (u_1) (u_2)}{(1+\gamma)^2} +\frac{(u^3) (u_1)}{1+\gamma }, \hspace{58pt} \end{eqnarray} \tag{ A35 }$$

$$\begin{eqnarray} a_{65} &=& \frac{(u^2) (u^3) (u_1) (u_3)}{(1+\gamma)^2} +\frac{(u^2)(u_1)}{1+\gamma }, \hspace{58pt} \end{eqnarray} \tag{ A36 }$$

$$\begin{eqnarray} a_{66} &=& 1+\frac{(u^2) (u^3) (u_2) (u_3)}{(1+\gamma)^2} +\frac{(u^2) (u_2)+(u^3)(u_3)}{1+\gamma }.\nonumber \\ \end{eqnarray} \tag{ A37 }$$

In our numerical code, the matrix A is inverted using LU-decomposition method.

Also (∂P^ij)/(∂E_r) and (∂P^ij)/(∂F^k_r) are obtained by solving the following equations:

$$\begin{eqnarray} &&A \frac{\partial p}{\partial E_r} = \frac{\partial r}{\partial E_r},\quad A \frac{\partial p}{\partial F^k_r} = \frac{\partial r}{\partial F^k_r}, \end{eqnarray} \tag{ A38 }$$

where

$$\begin{eqnarray} \frac{\partial R^{ij}}{\partial E_r} &=& \frac{\gamma ^2}{3}\delta ^{ij} - u^i u^j, \hspace{110pt} \end{eqnarray} \tag{ A39 }$$

$$\begin{eqnarray} \frac{\partial R^{ij}}{\partial F_r^k} &=&-\frac{2}{3} \gamma \delta ^{ij}u^k +u^i \delta ^{j}_k + u^j\delta ^{i}_k + \frac{2}{1+\gamma }u^i u^j u_k. \nonumber \\ \end{eqnarray} \tag{ A40 }$$

Since the matrix A is a function of u and independent of E_r and F^k_r, matrix A⁻¹ is computed once before calculating P^ij_r, (∂P^ij_r)/(∂E_r), and (∂P^ij_r)/(∂F_r^k).