A radiative transfer module for relativistic magnetohydrodynamics in the PLUTO code

TL;DR

This paper introduces a new module within the PLUTO code for solving relativistic radiation magnetohydrodynamics equations, capable of handling complex astrophysical phenomena with various coordinate systems and adaptive grids.

Contribution

It implements a novel numerical scheme with an implicit-explicit approach and a new HLLC solver for radiation transport, enhancing the capabilities of the PLUTO code.

Findings

Successfully demonstrated through numerical benchmarks

Handles both free-streaming and diffusion limits effectively

Supports multidimensional and adaptive grid computations

Abstract

We present a numerical implementation for the solution of the relativistic radiation hydrodynamics and magnetohydrodynamics equations, designed as an independent module within the freely available code PLUTO. The radiation transfer equations are solved under the grey approximation and imposing the M1 closure, which allows the radiation transport to be handled in both the free-streaming and diffusion limits. Equations are integrated following an implicit-explicit scheme, where radiation-matter interaction terms are integrated implicitly, whereas transport and all of the remaining source terms are solved explicitly by means of the same Godunov-type solvers included in PLUTO. Among these, we introduce a new Harten-van Leer-contact (HLLC) solver for optically thin radiation transport. The code is suitable for multidimensional computations in Cartesian, spherical and cylindrical coordinates, using either a single processor or parallel architectures. Adaptive grid computations are also made possible, by means of the CHOMBO library. The algorithm performance is demonstrated through a series of numerical benchmarks by investigating various different configurations with a particular emphasis on the behavior of the solutions in the free-streaming and diffusion limits.

Figures

Tables

Table 1. Table 1: Parameters used in the shock tests, in code units. The subscripts R 𝑅 R and L 𝐿 L correspond, respectively, to the initial conditions for x > 0 𝑥 0 x>0 and x < 0 𝑥 0 x<0 .

Test	${\rho }_L$	$p_{g,L}$	$u_L^x$	${\tilde{E}}_{r,L}$	${\rho }_R$	$p_{g,R}$	$u_R^x$	${\tilde{E}}_{r,R}$	$\Gamma$	$\kappa$
1	$1.0$	$3.0\times {10}^{-5}$	$0.015$	$1.0\times {10}^{-8}$	$2.4$	$1.61\times {10}^{-4}$	$6.25\times {10}^{-3}$	$2.51\times {10}^{-7}$	$5/3$	$0.4$
2	$1.0$	$4.0\times {10}^{-3}$	$0.25$	$2.0\times {10}^{-5}$	$3.11$	$0.04512$	$0.0804$	$3.46\times {10}^{-3}$	$5/3$	$0.2$
3	$1.0$	$60.0$	$10.0$	$2.0$	$8.0$	$2.34\times {10}^3$	$1.25$	$1.14\times {10}^3$	$2$	$0.3$
4	$1.0$	$6.0\times {10}^{-3}$	$0.69$	$0.18$	$3.65$	$3.59\times {10}^{-2}$	$0.189$	$1.3$	$5/3$	$0.08$

Equations

$$\begin{array}{rl}\frac{\partial I_{\nu }(t,\mathbf{x},\mathbf{n})}{\partial t}& +\mathbf{n}\cdot \nabla I_{\nu }(t,\mathbf{x},\mathbf{n})\\ & ={\eta }_{\nu }(t,\mathbf{x},\mathbf{n})-{\chi }_{\nu }(t,\mathbf{x},\mathbf{n})I_{\nu }(t,\mathbf{x},\mathbf{n}).\end{array}$$

$${\chi }_{\nu }(t,\mathbf{x},\mathbf{n})={\kappa }_{\nu }(t,\mathbf{x},\mathbf{n})+{\sigma }_{\nu }(t,\mathbf{x},\mathbf{n}),$$

$$T^{\mu \nu }=T_g^{\mu \nu }+T_{em}^{\mu \nu }+T_r^{\mu \nu }.$$

$$T_g^{\mu \nu }=\rho h{u}^{\mu }{u}^{\nu }+p_g{\eta }^{\mu \nu },$$

$$T_{em}^{\mu \nu }=F^{\mu \alpha }{F}_{\alpha }^{\nu }-\frac{1}{4}{\eta }^{\mu \nu }{F}_{\alpha \beta }{F}^{\alpha \beta },$$

$$F^{\mu \nu }=\left(\begin{array}{cccc}0& -E_1& -E_2& -E_3\\ E_1& 0& -B_3& B_2\\ E_2& B_3& 0& -B_1\\ E_3& -B_2& B_1& 0\end{array}\right),$$

$$T_r^{\alpha \beta }={\int }_0^{\infty }\mathrm{d}\nu \oint \mathrm{d}\Omega I_{\nu }(t,\mathbf{x},\mathbf{n})n^{\alpha }{n}^{\beta },$$

$$T_{r}=\left(\begin{array}[]{cc}E_{r}&F_{r}^{i}\\ F_{r}^{j}&P_{r}^{ij}\\ \end{array}\right),$$

$$E_r$$

$$F_r^i$$

$$P_r^{ij}$$

$${\nabla }_{\mu }(\rho u^{\mu })=0$$

$${\nabla }_{\mu }{T}^{\mu \nu }=0.$$

$${\nabla }_{\mu }\left(T_g^{\mu \nu }+T_{em}^{\mu \nu }\right)=-{\nabla }_{\mu }{T}_r^{\mu \nu }\equiv G^{\nu }$$

$$G^{\mu }={\int }_0^{\infty }\mathrm{d}\nu \oint \mathrm{d}\Omega \left({\chi }_{\nu }{I}_{\nu }-{\eta }_{\nu }\right)n^{\mu }.$$

$${\eta }_{\nu }(t,\mathbf{x},\mathbf{n})={\kappa }_{\nu }{B}_{\nu }(T)+{\sigma }_{\nu }{J}_{\nu },$$

$$T=\frac{\mu m_p{p}_g}{k_B\rho },$$

$${\tilde{G}}^{\mu }=\rho \left[\kappa {\tilde{E}}_r-4\pi \kappa B(T),\chi {\tilde{\mathbf{F}}}_r\right]$$

$$\begin{array}{rl}{G}^{\mu }=& -\kappa \rho \left(T_r^{\mu \alpha }{u}_{\alpha }+4\pi B(T)u^{\mu }\right)\\ & -\sigma \rho \left(T_r^{\mu \alpha }{u}_{\alpha }+T_r^{\alpha \beta }{u}_{\alpha }{u}_{\beta }{u}^{\mu }\right).\end{array}$$

$$\frac{\partial \left(\rho \gamma \right)}{\partial t}+\nabla \cdot \left(\rho \gamma \mathbf{v}\right)$$

$$\frac{\partial \mathcal{E}}{\partial t}+\nabla \cdot \left(\mathbf{m}-\rho \gamma \mathbf{v}\right)$$

$$\frac{\partial \mathbf{m}}{\partial t}+\nabla \cdot \left(\rho h{\gamma }^2\mathbf{vv}-\mathbf{BB}-\mathbf{EE}\right)+\nabla p$$

$$\frac{\partial \mathbf{B}}{\partial t}+\nabla \times \mathbf{E}$$

$$\frac{\partial E_r}{\partial t}+\nabla \cdot {\mathbf{F}}_r$$

$$\frac{\partial {\mathbf{F}}_r}{\partial t}+\nabla \cdot P_r$$

$$p=p_g+\frac{{\mathbf{E}}^2+{\mathbf{B}}^2}{2},$$

$$\mathbf{m}=\rho h{\gamma }^2\mathbf{v}+\mathbf{E}\times \mathbf{B},$$

$$\mathcal{E}=\rho h{\gamma }^2-p_g-\rho \gamma +\frac{{\mathbf{E}}^2+{\mathbf{B}}^2}{2},$$

$$h=1+\frac{\Gamma }{\Gamma -1}\Theta ,$$

$$h=\frac{5}{2}\Theta +\sqrt{1+\frac{9}{4}{\Theta }^2},$$

$$P_r^{ij}=D^{ij}{E}_r,$$

$$D^{ij}=\frac{1-\xi }{2}{\delta }^{ij}+\frac{3\xi -1}{2}{n}^i{n}^j,$$

$$\xi =\frac{3+4f^2}{5+2\sqrt{4-3f^2}},$$

$$||{\mathbf{F}}_r||\le E_r,$$

Full text

radiative transfer — magnetohydrodynamics (MHD) — relativity — methods: numerical — adaptive mesh refinement

1 Introduction

Radiative transfer is of great relevance in many different physical systems, ocurring in a broad range of size scales. In the context of astrophysics, for instance, it is of fundamental importance in the modeling of star atmospheres (Kippenhahn et al., 2012), pulsars (Becker et al., 2009), supernovae (Fryer et al., 2007), and black hole accretion disks (Thorne, 1974). In some high-energy environments (e.g. gamma-ray bursts), matter can emit radiation while being accelerated to relativistic speeds, sometimes being simultaneously subject to strong electromagnetic (EM) fields (see, e.g., the review by Mészáros, 2006). In situations where none of these effects can be disregarded, it is convenient to rely on numerical schemes that are able to deal with them simultaneously.

Several numerical methods for radiation transport found in the literature are based on different solutions of the radiative transfer equation (see, e.g., Mihalas et al., 1999), which provides a simplified yet strong formalism to the problem of radiation transport, absorption, and emission in presence of matter. This approach neglects every wave-like behavior of photons, and focuses, instead, on energy and momentum transport. Regardless of its simplicity, solving the frequency-dependent radiative transfer equation is not a trivial task, due to the high degree of nonlinearity present in the underlying mathematical description and the number of variables involved in it. For this reason, several simplified schemes can be found in the literature. Typical examples are the post-processing of ideal hydrodynamical calculations, (sometimes used in cases where radiation back-reaction can be neglected, see, e.g., Mimica et al., 2009), and Monte Carlo methods (see e.g. Mazzali et al., 1993; Kromer et al., 2009), where radiation densities and fluxes are computed by following the evolution of a large number of effective ‘photon packets’ along selected trajectories.

An alternative to these methods, followed throughout our work, is given by the moment approach to the radiative transfer equation. This consists in taking succesive angular moments of this equation, in the same fashion as the hydrodynamics (HD) and magnetohydrodynamics (MHD) can be obtained from the collisional Boltzmann-Vlasov equation (see e.g. Goedbloed et al., 2004). The resulting scheme provides an extension to relativistic and non-relativistic MHD, that can be used to compute the evolution of the total radiation energy density and its flux, considering its interaction with a material fluid. In this work we focus on the relativistic case, to which we refer as relativistic radiation MHD (Rad-RMHD henceforth). The model involves a series of additional approximations which we now describe. First, as proposed by Levermore (1984), we close the system of equations by assuming that the radiation intensity is isotropic in a certain reference frame (the M1 closure). In addition, we consider the fluid to be a perfect conductor, and assume the validity of the equations of ideal MHD for the interaction between matter and EM fields. Lastly, we adopt an effective gray body approximation, by replacing the opacity coefficients for a set of conveniently chosen frequency-averaged values.

Our implementation has been built as a supplementary module in the multiphysics, multialgorithm, high-resolution code PLUTO, designed for time-dependent explicit computations of either classical, relativistic unmagnetized or magnetized flows (Mignone et al., 2007a). The new module is fully parallel, has been adapted to all available geometries (Cartesian, cylindrical and spherical) and supports calculations on adaptively refined grids using the standard PLUTO-CHOMBO framework (see Mignone et al., 2012; Adams et al., 2015). In addition, we have introduced a new HLLC-type Riemann solver, suitable for optically thin radiation transport. In particular, our scheme is based on the HLLC solver for relativistic HD by Mignone & Bodo (2005) and it is designed to improve the code’s ability to resolve contact discontinuities when compared to HLL (Harten-van Leer) formulations (see e.g. Toro, 2013).

To integrate the transport terms of the equations of Rad-RMHD, our implementation employs the same sort of explicit methods used in PLUTO for the non-radiative case. However, gas-radiation interaction is treated differently, since this process may occur in times that are much shorter than the dynamical times; for instance, when matter is highly opaque. Hence, direct explicit integration of the interaction terms would lead to prohibitively small time steps and inefficient calculations. For this reason, our method of choice relies on Implicit-Explicit (IMEX) Runge-Kutta methods (Pareschi & Russo, 2005) whereby spatial gradients are treated explicitly while point-local interaction terms are integrated via an implicit scheme.

Similar approaches in the context of radiation HD and MHD have been followed by González et al. (2007), Commerçon et al. (2011), Roedig et al. (2012), Sadowski et al. (2013), Skinner et al. (2013), Takahashi et al. (2013), McKinney et al. (2014) and Rivera-Paleo et al. (2016). In particular, it is our intention to include our module in the following freely distributed versions of PLUTO.

This paper is structured as follows: in Section 2, we provide a brief summary of radiative transfer and the relevant equations used in this work, while in Section 3 we give a description of the implemented algorithms. Section 4 shows the code’s performance on several selected tests, and in Section 5 we summarize the main results of our work.

2 Radiation hydrodynamics

2.1 The equation of radiative transfer

In this section we outline the basic derivation that leads to the equations of Rad-RMHD, which are described in Section 2.3. We follow the formalism shown in Mihalas et al. (1999), taking as a starting point the radiative transfer equation,

[TABLE]

In this framework, photons are treated as point-like wave packets, that can be instantly emitted or absorbed by matter particles. As outlined in the introduction, this approach rules out effects due to the wave-like nature of light such as diffraction, refraction, dispersion, and polarization, and takes care only of energy and momentum transport (see e.g. Pomraning, 1973). Macroscopic EM fields, however, do not get such treatment along this work, and are instead regarded separately as classical fields.

Equation (1) describes the evolution of the radiation specific intensity $I_{\nu}$, defined as the amount of energy per unit area transported in a time interval $dt$ through an infinitesimal solid angle around the direction given by $\mathbf{n}$, in a range of frequencies between $\nu$ and $\nu+d\nu$. The quantities on the right hand side of this equation describe the interaction of the gas with the radiation field. The function $\eta_{\nu}$, known as emissivity, accounts for the energy released by the material per unit length, while the last term, proportional to $I_{\nu}$, measures the energy removed from the radiation field, also per unit length. The total opacity $\chi_{\nu}$ comprises absorption and scattering in the medium:

[TABLE]

where $\kappa_{\nu}$ and $\sigma_{\nu}$ are, respectively, the absorption and scattering frequency-dependent opacities.

Solving Equation (1) in the presented form is not a trivial task since integration must be in general carried out considering the dependency of $I_{\nu}$ on multiple variables $(t,\mathbf{r},\nu,\mathbf{n})$, while concurrently taking into account changes in the moving material. It also requires a precise knowledge of the functions $\eta_{\nu}$ and $\chi_{\nu}$, including effects such as the anisotropy caused by the Doppler shift. Instead of attempting a full solution, we adopt a frequency-integrated moment-based approach: we integrate Equation (1) over the frequency domain and take convenient averages in the angle -the moments- that can be naturally introduced in the equations of hydrodynamics. This procedure is described in the next section.

2.2 Energy-momentum conservation and interaction terms

We now explicitly derive the set of conservation laws describing the coupled evolution of fluid, EM, and radiation fields. While MHD quantities and radiation fields are calculated in an Eulerian frame of reference, absorption and scattering coefficients are best obtained in the fluid’s comoving frame (comoving frame henceforth), following the formalism described in Mihalas et al. (1999). The convenience of this choice relies on the fact that the opacity coefficients can be averaged easily without taking into account anisotropies due to a non-null fluid’s velocity, while the hyperbolic form of the conservation equations is kept. In this formalism, we split the total energy-momentum-stress tensor $T^{\mu\nu}$ into a gas, EM, and a radiative contribution:

[TABLE]

The first of these can be written as

[TABLE]

where $u^{\mu}$ is the fluid’s four-velocity and $\eta^{\mu\nu}$ is the Minkowski tensor, while $\rho$, $h$, and $p_{g}$ are, respectively, the fluid’s matter density, specific enthalpy, and pressure, measured in the comoving frame (our units are chosen so that $c=1$).

This expression is valid as long as material particles are in local thermal equilibrium (LTE henceforth), which is one of the assumptions of the hydrodynamical treatment.

The electromagnetic contribution is given by the EM stress-energy tensor:

[TABLE]

where the components of the field tensor $F^{\mu\nu}$ are given by

[TABLE]

where $E_{i}$ and $B_{i}$ are, respectively, the components of the electric and magnetic fields.

Lastly, $T_{r}^{\mu\nu}$ can be written in terms of the specific intensity $I_{\nu}$, as

[TABLE]

where $n^{\mu}\equiv(1,\mathbf{n})$ denotes the direction of propagation, $\mathrm{d}\nu$ the differential frequency, and $\mathrm{d}\Omega$ the differential solid angle around $\mathbf{n}$. This expression, by definition covariant (see e.g. Mihalas et al., 1999), can be shortened as

[TABLE]

where

[TABLE]

are the first three moments of the radiation field, namely, the radiation energy density, the flux, and the pressure tensor. In our scheme, we follow separately the evolution of $E_{r}$ and $F_{r}^{i}$, and define the pressure tensor in terms of these fields by means of a closure relation, as it is described in Section 2.4.

Following these definitions, and imposing conservation of mass, total energy, and momentum, we have

[TABLE]

and

[TABLE]

From equations (3) and (13), we immediately obtain

[TABLE]

where $G^{\mu}$ - the radiation four-force density - is computed by integrating Eq. (1) over the frequency and the solid angle, as

[TABLE]

The equations of Rad-RMHD can then be derived from Eq. (14), where the term $G^{\mu}$ accounts for the interaction between radiation and matter.

The previous expression can be simplified in the comoving frame provided some conditions are met. Firstly, we assume coherent and isotropic scattering and calculate the total comoving emissivity as

[TABLE]

where $B_{\nu}(T)$ is the Planck’s spectral radiance at a temperature $T$, while $J_{\nu}$ is the angle-averaged value of $I_{\nu}$. The temperature can be determined from the ideal gas law

[TABLE]

where $\mu$ is the mean molecular weight, $m_{p}$ is the proton mass, and $k_{B}$ the Boltzmann constant.

We can then insert these expressions in Eq. (15) and replace the opacities by their corresponding frequency-averaged values, such as the Planck and Rosseland means (see e.g. Mihalas et al., 1999; Skinner et al., 2013). In this way, we obtain the following comoving-frame source terms

[TABLE]

where $B(T)=\sigma_{\textsc{SB}}T^{4}/\pi c$, $\sigma_{\textsc{SB}}$ is the Stefan-Boltzmann constant, and $\chi$, $\kappa$, and $\sigma$ are the mentioned frequency-averaged opacities, per unit density. In the code, these can either be set as constants, or defined by the user as functions of any set of fields (for instance, $\rho$ and $T$). From now on and except for the opacity coefficients, we label with a tilde sign quantities in the comoving frame.

Finally, $G^{\mu}$ can be obtained in the Eulerian frame by means of a Lorentz boost applied to Equation (18) (see e.g. McKinney et al., 2014):

[TABLE]

2.3 The equations of Rad-RMHD

Assuming ideal MHD for the interaction between matter and EM fields, we obtain the equations of Rad-RMHD in quasi-conservative form:

[TABLE]

where $\mathbf{v}$ is the fluid’s velocity, $\gamma$ is the Lorentz factor, $\mathbf{B}$ the mean magnetic field, $\mathbf{E}=-\mathbf{v}\times\mathbf{B}$ the electric field. In addition, we have introduced the quantities

[TABLE]

[TABLE]

[TABLE]

which account, respectively, for the total pressure, momentum density, and energy density of matter and EM fields. The system (20)-(25) is subject to the constraint $\nabla\cdot\mathbf{B}=0$, and the non-magnetic case (Rad-RHD) is recovered by taking the limit $\mathbf{B}\rightarrow\mathbf{0}$ in the previous expressions.

In our current scheme, Equations (20) to (25) can be solved in Cartesian, cylindrical or spherical coordinates.

2.4 Closure relations

An additional set of relations is required in order to close the system of Equations (20)–(25). An equation of state (EoS) provides closure between thermodynamical quantities and it can be specified as the constant-$\Gamma$ law

[TABLE]

or the Taub-Mathews equation, introduced by Mathews (1971),

[TABLE]

where $\Theta=p_{g}/\rho$. The properties of these equations are known and can be found, e.g., in Mignone et al. (2007b).

A further closure relation is needed for the radiation fields, i.e., an equation relating $P_{r}^{ij}$ to $E_{r}$ and $\mathbf{F}_{r}$. We have chosen to implement the M1 closure, proposed by Levermore (1984), which permits to handle both the optically thick and optically thin regimes. In this scheme, it is assumed that $I_{\nu}$ is isotropic in some inertial frame, where the radiation stress-energy tensor takes the form $T^{\prime\mu\nu}_{r}=\operatorname{diag}(E_{r}^{\prime},E_{r}^{\prime}/3,E_{r}^{\prime}/3,E_{r}^{\prime}/3)$. This leads to the following relations, which hold in any frame:

[TABLE]

[TABLE]

[TABLE]

where now $\bm{n}=\mathbf{F}_{r}/||\mathbf{F}_{r}||$ and $f=||\mathbf{F}_{r}||/E_{r}$, while $\delta^{ij}$ is the Kronecker delta.

These relations are well behaved, as Equations (9) and (10) provide an upper limit to the flux, namely

[TABLE]

and therefore $0\leq f\leq 1$.

In our scheme, we apply Equations (31)-(33) in the laboratory frame. In the diffusion limit, namely, if $||\mathbf{F}_{r}||\ll E_{r}$, this closure leads to $P_{r}^{ij}=\left(\delta^{ij}\middle/3\right)E_{r}$, which reproduces an isotropic specific intensity known as Eddington limit. Likewise, in the free-streaming limit given by $||\mathbf{F}_{r}||\rightarrow E_{r}$, the pressure tensor tends to $P_{r}^{ij}=E_{r}\,n^{i}n^{j}$, which corresponds to a delta-like $I_{\nu}$ pointing in the same direction and orientation as $\mathbf{F}_{r}$.

We point out that, even though both the free-streaming and the diffusion limits are reproduced correctly, the M1 closure may fail in some cases, since it implicitly assumes, for example, that the intensity $I_{\nu}$ is axially symmetric in every reference frame with respect to the direction of $\mathbf{F}_{r}$. This is not the case, for example, when two or more radiation sources are involved, in which case direct employment of the M1 closure may become inaccurate, leading to instabilities (see e.g. Sadowski et al., 2013; Skinner et al., 2013).

3 Numerical scheme

For numerical purposes we write equations (20)-(25) in conservative form as

[TABLE]

where ${\cal U}\equiv\left(\rho\gamma,\,\mathcal{E},\,\mathbf{m},\,\mathbf{B},\,E_{r},\,\mathbf{F}_{r}\right)^{\intercal}$ is an array of conserved quantities, $\mathsf{F}({\cal U})$ is the flux tensor and ${\cal S}\equiv\left(0,G^{0},\mathbf{G},\mathbf{0},-G^{0},-\mathbf{G}\right)^{\intercal}$ contains the radiation-matter interaction terms. The explicit expressions of $\mathsf{F}$ can be extracted from Equations (20)-(25).

As outlined in the introduction, gas-radiation interaction may occur in timescales that are much smaller than any dynamical characteristic time and an explicit integration of the interaction terms would lead either to instabilities or to excessively large computing times. For this reason, the time discretization of Equations (35) is achieved by means of IMEX (implicit-explicit) Runge-Kutta schemes (see e.g. Pareschi & Russo, 2005). In the presented module, the user can choose between two different IMEX schemes, as described in Section 3.1.

In our implementation of the IMEX formalism, fluxes and geometrical source terms are integrated explicitly by means of standard shock-capturing Godunov-type methods, following a finite volume approach. Fluxes are thus evaluated at cell interfaces by means of a Riemann solver between left and right states properly reconstructed from the two adjacent zones. Geometrical source terms can be obtained at the cell center or following the approach outlined in Mignone (2014). This explicit step is thoroughly described in Section 3.2. Within this stage, we have included a new Riemann solver for radiation transport, which we introduce in Section 3.3. On the other hand, the integration of $G^{\mu}$ is performed implicitly through a separate step (the implicit step), as described in Section 3.4.

3.1 Implemented IMEX schemes

A commonly used second-order scheme is the IMEX-SSP2(2,2,2) method by Pareschi & Russo (2005) which, when applied to (35), results in the following discrete scheme:

[TABLE]

Here ${\cal U}$ is an array of volume averages inside the zone $i,j,k$ (indices have been omitted to avoid cluttered notations), $n$ denotes the current step number, $\Delta t^{n}$ is the time step, $a=1-\sqrt{2}/2$, and the operator ${\cal R}$, which approximates the contribution of $(-\nabla\cdot\mathsf{F})$, is computed in an explicit fashion in terms of the conserved fields as detailed in Section 3.2. Potentially stiff terms

• i.e., those poportional to $\kappa$ and $\sigma$ - are included in the operator ${\cal S}$ which is solved implicitly during the first and second stages in Eq. (36).

An alternative scheme which we also consider in the present context is the following scheme (IMEX1 henceforth):

[TABLE]

This method is an extension to the second-order total variation diminishing Runge-Kutta scheme (RK2) by Gottlieb et al. (1996), where we have just added an implicit step after every flux integration. In the same way, we have included in the code a third-order version of this scheme that extends the third-order Runge-Kutta scheme by the same authors. Both the second- and third-order of this method are similar to those described in McKinney et al. (2014).

Using general methods for IMEX-RK schemes (see e.g. Pareschi & Russo, 2005), it can be shown that IMEX-SSP2(2,2,2) and IMEX1 are of order 2 and 1 in time and L- and A-stable respectively, which makes IMEX-SSP2(2,2,2) a seemingly better option when it comes to the schemes’ stability. However, as we have observed when testing the module, the explicit addition of previously-calculated source terms in the last step of IMEX-SSP2(2,2,2) can cause inaccuracies whenever interaction terms are stiff and there are large differences in the orders of magnitude of matter and radiation fields (see Sections 4.5.2 and 4.6). Contrarily, IMEX1 seems to have better positivity-preserving properties and a higher accuracy in those cases. In general, as it is shown in Section 4, we have obtained equivalent results with both methods in every other case. Whenever source terms can be neglected, both methods reduce to the standard RK2, which makes them second-order accurate in time for optically thin transport.

3.2 Explicit step

In order to compute the explicit operator ${\cal R}$, we implement a standard reconstruct-solve-update strategy (see e.g. Rezzolla et al., 2013). First, the zone averages ${\cal U}$ are used to compute cell-centered values of a set of primitive fields, defined as

[TABLE]

Although this is a straightforward step for the radiation fields, as in their case primitive and conserved quantities coincide, this is not the case for the remaining variables. Primitive fields are obtained from conservative ones by means of a root-finding algorithm, paying special attention to avoiding problems related to small number handling that arise when large Lorentz factors are involved. To perform this conversion, we follow the procedure detailed in Mignone et al. (2007b).

Next, primitive fields are reconstructed to zone interfaces (reconstruction step). In more than one dimensions, reconstruction is carried direction-wise. In order to avoid spurious oscillations next to discontinuities and steep gradients, reconstruction must use slope limiters in order to satisfy monotonicity constraints. During this step, some physical constraints are imposed, such as gas pressure positivity, an upper boundary for the velocity given by $||\mathbf{v}||<1$, and the upper limit to the radiation flux given by Equation (34).

The reconstruction step produces left and right discontinuous states adjacent to zone interfaces, which we denote with ${\cal V}_{L}$ and ${\cal V}_{R}$. This poses a local initial-value problem that is solved by means of an approximate Riemann solver, whose outcome is an estimation of the fluxes on each interface. In our implementation, the user can choose among three of these methods. The simplest one of these is the Lax-Friedrichs-Rusanov solver (see e.g. Toro, 2013), which yields the following flux:

[TABLE]

In this expression, ${\cal U}_{L/R}$ and ${\cal F}_{L/R}=\hat{\bm{e}}_{d}\cdot\mathsf{F}({\cal U}_{L/R})$ are the conserved fields and flux components in the coordinate direction $\hat{\bm{e}}_{d}$ (here $d=x,y,z$ in Cartesian coordinates or $d=r,\theta,\phi$ in spherical coordinates) evaluated at the left and right of the interface, while $\lambda_{\max}$ is the fastest signal speed at both sides, computed using both ${\cal V}_{L}$ and ${\cal V}_{R}$. A less diffusive option is given by an HLL solver (Harter-Lax-van Leer, see e.g. Toro, 2013) introduced by González et al. (2007). In this case, fluxes are computed as

[TABLE]

where $\lambda_{L}$ and $\lambda_{R}$ are, respectively, the minimum and maximum characteristic signal speeds, taking into account both ${\cal V}_{L}$ and ${\cal V}_{R}$ states. Finally, a third option is given by an HLLC solver that estimates the HD (MHD) fluxes as described in Mignone & Bodo (2005) (see also Mignone & Bodo, 2006), and the radiation fluxes as described in Section 3.3.

From Eqs. (20)-(25) we can see that, if interaction terms are disregarded, the equations of Rad-RMHD can be divided into two independent systems, one corresponding to the equations of relativistic MHD and the other to those of radiation transport. Hence, we can expect the maximum and minimum signal speeds of both systems to be, in the frozen limit111 In the theory of stiff relaxation systems, the frozen limit refers to the small time step regime, when the effect of source terms on the characteristic waves is still negligible., different. In view of this, we compute the fluxes independently for each subsystem of equations obtaining the speeds shown in Appendix A. In this manner, as it is pointed out in Sadowski et al. (2013), we avoid the excessive numerical diffusion that appears when the same signal speeds are used to update both radiation and MHD fields. This has been verified in our tests.

Once the fluxes are obtained, we can compute the operator ${\cal R}$ which, in the direction $d$, reads

[TABLE]

where $A^{d}_{\pm}$ are the cell’s right ($+$) and left ($-$) interface areas and $\Delta V^{d}$ is the cell volume in that direction (see Mignone et al., 2007a). Here ${\cal S}^{d}_{e}({\cal U})$ accounts for geometrical terms that arise when the divergence is written in different coordinate systems. The full operator ${\cal R}$ is in the end computed as $\sum_{d}\mathcal{R}_{d}$.

Once the update of the conserved variables is completed, the time step is changed using the maximum signal speed computed in the previous step, according to the Courant-Friedrichs-Lewy condition (Courant et al., 1928):

[TABLE]

where $\Delta l_{\min}^{d}$ and $\lambda_{\max}^{d}$ are, respectively, the minimum cell width and maximum signal speed along the direction $d$, and $C_{a}$, the Courant factor, is a user-defined parameter.

Finally, when magnetic fields are included, the divergence-free condition can be enforced using either the constrained transport method (Balsara & Spicer, 1999; Londrillo & Del Zanna, 2004) or hyperbolic divergence cleaning (Dedner et al., 2002; Mignone et al., 2010a, b). Both methods are available in the code.

3.3 HLLC solver for radiation transport

We now present a novel Riemann solver for the solution of the homogeneous radiative transfer equation. To this purpose, we consider the subsystem formed by Eqs. (24)-(25) by neglecting interaction terms and restrict our attention to a single direction, chosen to be the $x$ axis, without loss of generality. In Cartesian coordinates, the resulting equations take the form

[TABLE]

where ${\cal U}_{r}=(E,\,\mathbf{F})^{\intercal}$ while $\Phi=(F_{x},\,P_{xx},\,P_{yx},\,P_{zx})^{\intercal}$ and we have omitted the subscripts $r$ for clarity purposes (we shall maintain that convention throughout this section). From the analysis carried out in Appendix A, we know that the Jacobian $\mathbf{J}^{x}$ of this system has three different eigenvalues $\{\lambda_{1},\lambda_{2},\lambda_{3}\}$, satisfying $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$. Since the system is hyperbolic (see e.g. Toro, 2013), the breaking of an initial discontinuity will involve the development of (at most) as many waves as the number of different eigenvalues. On this basis, we have implemented a three-wave Riemann solver.

Following Berthon et al. (2010), we define the following fields:

[TABLE]

where $\xi$ is given by Eq. (33). With these definitions, the fluxes in Eq. (43) can be written as

[TABLE]

and $F_{x}$ can be shown to satisfy $F_{x}=\left(E+\Pi\right)\beta_{x}$. These expressions are similar to those of relativistic hydrodynamics (RHD henceforth), where $\beta_{x}$, $\Pi$ and $\mathbf{F}$ play, respectively, the role of $v_{x}$, $p_{g}$ and $\mathbf{m}$ while $E$ is tantamount to the total energy. With the difference that there is no field corresponding to density, the equations are exactly the same as those corresponding to energy-momentum conservation of a fluid, with a different closure relation.

With this in mind, we follow analogous steps to those in Mignone & Bodo (2005) in order to construct a HLLC solver for the system defined by Equations (43). In this case, instead of the intermediate constant state considered in the HLL solver, we include an additional middle wave (the analog of a “contact” mode) of speed $\lambda^{*}$ that separates two intermediate states $\mathcal{U}^{*}_{L}$ and $\mathcal{U}^{*}_{R}$, where

[TABLE]

In this way, the full approximate solution verifies

[TABLE]

The corresponding fluxes are

[TABLE]

States and fluxes are related by the Rankine-Hugoniot jump conditions across the outermost waves $\lambda_{S}$ ($S=L,R$),

[TABLE]

A similar condition must also hold across the middle wave so that, when Equation (49) is applied to all three waves, one has at disposal a system of 12 equations for the 17 unknowns (${\cal U}^{*}_{r,L}$, ${\cal U}^{*}_{r,R}$, $\Phi^{*}_{L}$, $\Phi^{*}_{R}$, and $\lambda^{*}$) and therefore further assumptions must be made. From the results of the tests performed with the HLL solver, we have verified that $\beta_{x}$ and $\Pi$ are conserved along the intermediate contact mode for all the obtained solutions. Noting that $\lambda_{2}(E,\mathbf{F})=\beta_{x}(E,\mathbf{F})$, it can be seen that, for a discontinuity of speed $\beta_{x}$ along which $\beta_{x}$ and $\Pi$ are continuous, the jump conditions (49) are satisfied, as pointed out in Berthon et al. (2010) and proven in Hanawa & Audit (2014). Thus, we impose the constraints $\lambda^{*}=\beta^{*}_{x,L}=\beta^{*}_{x,R}$ and $\Pi^{*}_{L}=\Pi^{*}_{R}$. These conditions are analogous to those satisfied by the contact discontinuity in RHD, across which $p_{g}$ and $v_{x}$ are conserved, and where the latter coincides with the propagation speed. Following Mignone & Bodo (2005), we assume that $\Phi^{*}$ can be written in terms of the five variables $(E^{*},\Pi^{*},\beta_{x}^{*},F^{*}_{y},F^{*}_{z})$ in the following way:

[TABLE]

where for consistency we have defined $F^{*}_{x}\equiv(E^{*}+\Pi^{*})\beta^{*}_{x}$. Under these constraints, the jump conditions across the middle wave are automatically satisfied, and Eq. (49) is reduced to the following system of 8 equations in 8 unknowns:

[TABLE]

which holds for both subscripts L and R (we shall maintain this convention in what follows). The first two equations in Eq. (51) can be turned into the following quadratic expression, from which $\lambda^{*}$ can be obtained:

[TABLE]

with

[TABLE]

Once $\lambda^{*}$ is known, we can compute $\Pi^{*}$ as

[TABLE]

and the remaining fields from Eq. (51). Similarly to the RHD counterpart, among the two roots of Equation (52) we must choose the only one that guarantees $\lambda^{*}\in[-1,1]$, which in our case corresponds to that with the minus sign. As shown in Appendix B, this definition of $\lambda^{*}$ satisfies Eq. (46). We have also checked by means of extensive numerical testing that the intermediate states $\mathcal{U}^{*}_{L}$ and $\mathcal{U}^{*}_{R}$ constructed in this way satisfy Equation (34), which guarantees the positivity of our HLLC scheme. However, unlike the RHD case, the coefficients $\{A_{L},B_{L},A_{R},B_{R}\}$ defined in Equations (53) and (54) can simultaneously be equal to zero, meaning that $\lambda^{*}$ can no longer be determined from Equation (52). This happens under the conditions $||\mathbf{F}||=E$ for both L and R, and $F_{xL}/||\mathbf{F}_{L}||\leq F_{xR}/||\mathbf{F}_{R}||$, in which case the jump conditions lead to the formation of vacuum-like intermediate states. We overcome this issue by switching the solver to the standard HLL whenever these conditions are met.

As for the HLL solver, signal velocities must be limited when describing radiation transfer in highly opaque materials in order to reduce numerical diffusion (see Appendix A). Whenever this occurs, we also switch to the standard HLL solver, and limit $\lambda_{L}$ and $\lambda_{R}$ according to Equation (A6). Hence, we can only expect the HLLC solver to improve the accuracy of the obtained solutions in optically thin regions of space, whereas the results should be the same for both HLL and HLLC everywhere else. Finally, although the use of the HLLC solver can reduce the numerical diffusion when compared to the HLL solver, this can cause spurious oscillations around shocks that would be damped with a more diffusive method. As for the HLLC solver for relativistic HD and MHD included in PLUTO, this problem can be reduced by implementing an additional flattening in the vicinity of strong shocks (see e.g. Mignone & Bodo, 2005).

3.4 Implicit step

We now describe the algorithm employed for the implicit integration of the radiation-matter interaction terms. A typical implicit step of an IMEX scheme (see Eqs. 36 and 37) takes the form

[TABLE]

where $s$ is a constant and primed terms denote some intermediate state value. Equation (56) shows that the mass density, computed as $\rho\gamma$, as well as the total energy and momentum densities, defined as $E_{tot}=\mathcal{E}+E_{r}$ and $\mathbf{m}_{tot}=\mathbf{m}+\mathbf{F}_{r}$, must be conserved during this partial update owing to the particular form of the source terms. This yields the following implicit relations between $\mathcal{V}$ and ${\cal U}_{r}$:

[TABLE]

We can then solve Eq. (56) in terms of the following reduced system:

[TABLE]

with $\mathcal{G}\equiv(G^{0},\mathbf{G})^{\intercal}$, where $G^{\mu}$ is given in Eq. (19). In Eq. (60), radiation fields can be regarded as functions of the MHD fields and vice-versa by means of Eq. (59), and therefore the system can be solved in terms of either one of these.

In order to solve Equation (60), we have implemented and compared three different multidimensional root finder algorithms, which we now describe.

Fixed-point method. This method (originally proposed by Takahashi et al., 2013) is based on iterations of ${\cal U}_{r}$ and follows essentially the same approach outlined by Palenzuela et al. (2009) in the context of resistive MHD. In this scheme all of the MHD primitive variables, as well as $D^{ij}$, are written at a previous iteration with respect to ${\cal U}_{r}$. In that manner, $\mathcal{G}$ can be written at a given iteration $m$ as

[TABLE]

where $\mathcal{M}$ is a matrix and $b$ a column vector, both depending on $\mathcal{V}$ and $D^{ij}$, and the numbers between parentheses indicate the iteration in which the fields are evaluated. Inserting this in Equation (60), the updated conserved fields can be computed as

[TABLE]

after which primitive fields can be updated using Eq. (59). 2. 2.

Newton’s method for radiation fields, implemented in Sadowski et al. (2013) and McKinney et al. (2014). This scheme consists in finding the roots of the nonlinear multidimensional function

[TABLE]

updating the radiation variables on each iteration as

[TABLE]

where we have defined the Jacobian matrix $\mathcal{J}$ as $\mathcal{J}_{ij}=\partial\mathcal{Q}_{i}/\partial{\cal U}_{r}^{j}$. The elements of $\mathcal{J}$ are computed numerically, taking small variations of the iterated fields. As in the fixed-point method, matter fields are computed from ${\cal U}_{r}$ for each step by means of an inversion of Eq. (59). 3. 3.

Newton’s method for matter fields, implemented in McKinney et al. (2014). This procedure is identical to the previous one, with the difference that in this case the iterated fields are the fluid’s pressure and the spatial components of its four-velocity, which we denote as $\mathcal{W}=(p_{g},\mathbf{u})^{\intercal}$. These are updated as

[TABLE]

where now $\mathcal{J}_{ij}=\partial\mathcal{Q}_{i}/\partial\mathcal{W}^{j}$ and $\mathcal{Q}$ is regarded as a function of $\mathcal{W}$. This scheme is much faster than the previous one, since the computation of ${\cal U}_{r}$ from $\mathcal{W}$ by means of Eq. (59) is now straightforward, and no longer requires a cumbersome inversion of conserved to primitive fields.

For each of these methods, iterations are carried out until convergence is reached by means of some error function. In the first of them, such function is chosen as the norm of the relative differences between successive values of $\mathcal{V}$, whereas in the last two of them it is defined as the norm of $\mathcal{Q}^{(m+1)}$. If $\mathcal{E}\ll E_{r}$, the errors of the matter fields can be large even when radiation fields converge, since Eq. (59) implies that $\mathcal{E}$ and $E_{r}$ have the same absolute error, as well as $\mathbf{m}$ and $\mathbf{F}_{r}$. Therefore, having small relative differences of $E_{r}$ does not guarantee the same for $\mathcal{E}$, which can lead to non-negligible inaccuracies if the second method is used. Equivalently, the same problem can occur whenever $\mathcal{E}\gg E_{r}$ if method 3 is chosen (see also McKinney et al., 2014). To overcome this issue, we have included in the code the option of adding to the convergence function the norm of the relative differences of $\mathcal{E}$ when using the second method, and of $E_{r}$ when using the third one. We have seen in the performed tests that the fixed-point method converges rather fast, meaning that the number of iterations that it requires frequently coincides with that obtained with the last two methods. This scheme has sufficed to perform all the tests carried out in this work, being often the fastest one when compared to the other two, having been overcome only by method 3 in a few cases.

4 Numerical Benchmarks

We show in this section a series of numerical benchmarks to verify the code performance, as well as the correctness of the implementation under different physical regimes and choices of coordinates. Unless otherwise stated, we employ the HLLC solver introduced in Section 3.3, the domain is discretized using a fixed uniform grid and outflow boundary conditions are imposed for all the fields. Magnetic fields are neglected in all the considered problems, except in Section 4.7. Furthermore, all the tests have been run with both the IMEX-SSP2(2,2,2) and IMEX1 methods, obtaining equivalent results unless indicated otherwise.

4.1 Riemann Problem for optically-thin radiation transport

We first validate the implemented radiation transport schemes when any interaction with matter is neglected. To this end, we have run several one-dimensional Riemann problems setting all the interaction terms to zero and focusing only on the evolution of the radiation fields. The initial setup of these consists of two regions of uniform $E_{r}$ and $\mathbf{F}_{r}$, separated by a discontinuity at $x=0$. The full domain is defined as the interval $[-20,20]$. We show here two of such tests, exploring the case $||\mathbf{F}_{r}||<E_{r}$ (test 1) and the free-streaming limit, $||\mathbf{F}_{r}||\simeq E_{r}$ (test 2).

In the first test, initial states are assigned at $t=0$ as

[TABLE]

The solution, plotted in Fig 1 at $t=20$ with a resolution of $2^{14}$ zones (solid black line), shows a three-wave pattern as it is expected from the eigenstructure of the radiation transport equations (see Section 3.3 and Appendix B). The left and right outermost waves are, respectively, a left-facing shock and a right-going expansion wave, while the middle wave is the analog of a contact wave. The fields $\Pi$ and $\beta_{x}$, defined in Section 3.3, are constant across the contact mode. On the same Figure, we show the solution obtained with the HLL and HLLC solvers at the resolution of 256 zones using a $1^{\rm st}$ order reconstruction scheme (see Mignone & Bodo, 2005). As expected, the employment of the HLLC solver yields a sharper resolution of the middle wave.

For the second test, the initial condition is defined as

[TABLE]

Results obtained with the $1^{\rm st}$-order scheme and the HLL and HLLC solvers are plotted in Fig. 2 together with the reference solution (solid black line) at $t=20$. As for the previous case, a three-wave pattern emerges, formed by two left- and right-going shocks and a middle contact wave. It can be also seen that $\Pi$ and $\beta_{x}$ are again continuous across the contact wave. Differences between HLLC and HLL are less pronounced than the previous case, with the HLL (HLLC) overestimating the left-going shock position by 50% (30%).

For both tests, we have conducted a resolution study covering the range $[2^{6},2^{10}]$ using $1^{\rm st}$- as well as $2^{\rm nd}$-order reconstructions making use of the second-order harmonic mean limiter by van Leer (1974). In Figure 3, we plot the error in L1-norm of $E_{r}$ (computed with respect to the reference solution) as functions of the resolution. The Courant number is $C_{a}=0.4$ for both cases. Overall, the HLLC yields smaller errors when compared to HLL, as expected. This discrepancy is more evident in the $1^{\rm st}-$order case and it is mitigated in the case of a $2^{\rm nd}$ order interpolant (a similar behavior is also found in Mignone & Bodo, 2005).

4.2 Free-streaming beam

A useful test to investigate the code’s accuracy for multidimensional transport is the propagation of a radiation beam oblique to the grid (see e.g. Richling et al., 2001; González et al., 2007). This problem is also useful to quantify the numerical diffusion that may appear when fluxes are not aligned with the axes. We again neglect the gas-radiation interaction terms, and follow solely the evolution of the radiation fields.

The initial setup consists of a square Cartesian grid of side $L=5$ cm, where the radiation energy density is set to $E_{r,0}=10^{4}$ erg cm*-3*. At the $x=0$ boundary, a radiation beam is injected by fixing $E_{r}=10^{8}\,E_{r,0}$ and $\mathbf{F}_{r}=(1/\sqrt{2},1/\sqrt{2})\,E_{r}$ for $y\in[0.30,0.44]$ cm. Thus, the injected beam satisfies the equality $||\mathbf{F}_{r}||=E_{r}$, which corresponds to the free-streaming limit. Outflow conditions are imposed on the remaining boundaries.

Again we compare the performance of the HLL and HLLC solvers, using the fourth-order linear slopes by Miller & Colella (2001) and resolutions of $150\times 150$ and $300\times 300$ zones. The Courant number is $C_{a}=0.4$. The energy density distribution obtained with the HLLC solver at the largest resolution is shown in Fig. 4 at $t=5\times 10^{-10}$ s. In every case, a beam forms and reaches the upper boundary between $x=4$ cm and $x=5$ cm, after crossing a distance equivalent to roughly $\sim 64$ times its initial width. Since no interaction with matter is considered, photons should be transported in straight lines. As already mentioned, the free-streaming limit corresponds to a delta-like specific intensity parallel to $\mathbf{F}_{r}$. Hence, photons are injected in only one direction, and the beam’s structure should be maintained as it crosses the computational domain. However, in the simulations, the beam broadens due to numerical diffusion before reaching the upper boundary. For this particular test, due to its strong discontinuities, we have seen that this effect is enhanced by the flattening applied during the reconstruction step in order to satisfy Equation (34), which is necessary for stability reasons.

In order to quantify this effect and its dependecy on the numerical resolution, we have computed several time-averaged $E_{r}(y)$ profiles along vertical cuts at different $x$ values. As an indicator of the beam’s width, we have computed for each $x$ the standard deviation of these profiles as

[TABLE]

with

[TABLE]

where the weighting function $f(y)$ is defined as

[TABLE]

being $\overline{E}_{r}$ the time-averaged value of $E_{r}$. We have then divided the resulting values of $\sigma_{y}$ by $\sigma_{y0}\equiv\sigma_{y}(x=0)$, in order to show the relative growth of the dispersion. The resulting values of $\sigma_{y}/\sigma_{y0}$ are shown in Fig. 4, where it can be seen that the beam’s dispersion grows with $x$. The difference between $\sigma_{y}/\sigma_{y0}$ and its ideal value ($\sigma_{y}/\sigma_{y0}\equiv 1$) gets reduced by a factor between 2 and 2.5 when the highest resolution is used. In the same figure, it can be seen that the dispersion is only slightly reduced when the HLLC solver is used instead of HLL. A similar plot of $\sigma_{y}/\sigma_{y0}$ is obtained with the second-order limiter by van Leer (1974), where the values of the relative dispersion increase roughly between $30\%$ and $40\%$, showing as in Section 4.1 that the accuracy of these methods not only depends on the chosen Riemann solver but it is also extremely sensitive to the chosen reconstruction scheme.

4.3 Radiation-matter coupling

In order to verify the correct integration of the interaction terms, we have run a test proposed by Turner & Stone (2001), in which matter and radiation approach thermal equilibrium in a homogeneous system. This is achieved by solving the Rad-RHD equations in a single-cell grid, thus removing any spatial dependence. In this configuration, due to the form of Equations (20)-(25), all the fields but the energy densities of both radiation and matter remain constant for $t>0$. Using conservation of total energy, the resulting equation for the evolution of the gas energy density (in cgs units) is

[TABLE]

This can be simplified if the chosen initial conditions are such that $E_{r}$ is constant throughout the system’s evolution. In that case, Equation (71) can be solved analytically, leading to an implicit relation between $\mathcal{E}$ and $t$ that can be inverted using standard methods.

We have run this test for two different initial conditions, using in both cases $\rho=10^{-7}$ g cm*-3*, $E_{r}=10^{12}$ erg cm*-3*, opacities $\kappa=0.4$ cm2 g*-1* and $\sigma=0$, and a mean molecular weight $\mu=0.6$. A constant-gamma EoS has been assumed, with $\Gamma=5/3$. We have chosen the initial gas energy density to be either $\mathcal{E}=10^{10}$ erg cm*-3* or $\mathcal{E}=10^{2}$ erg cm*-3*, which are, respectively, above and below the final equilibrium value, of around $7\times 10^{7}$ erg cm*-3*.

The gas energy density is plotted as a function of time for both conditions in Fig. 5. Simulations are started from $t=10^{-10}$ s, with an initial time step $\Delta t=10^{-10}$ s. An additional run between $t=10^{-16}$ s and $10^{-10}$ s is done for each initial condition with an initial $\Delta t=10^{-16}$ s, in order to show the evolution in the initial stage. In every case, the gas radiation energy goes through an initial constant phase that lasts until $t\sim 10^{-14}$ s, after which it varies towards the equilibrium value. Equilibrium occurs when the condition $E_{r}=4\pi B(T)$ is reached (see Eq. (71)), i.e., when the power emitted by the gas equals its energy absorption rate. This happens around $t\approx 10^{-7}$ s for both initial conditions. As shown in Fig. 5, the numerical solutions match the analytical ones in the considered time range.

4.4 Shock waves

We now study the code’s ability to reproduce general shock-like solutions without neglecting the interaction terms. To this purpose, we have reproduced a series of tests proposed by Farris et al. (2008). As in Section 4.1, we place a single initial discontinuity at the center of the one-dimensional domain defined by the interval $[-20,20]$. At $t=0$, both matter and radiation fields are constant on each side of the domain, and satisfy the condition for LTE between matter and radiation, that is, $\tilde{E}_{r}=4\pi B(T)$. Additionally, the fluxes on each side obey $\tilde{F}_{r}^{x}=0.01\times\tilde{E}_{r}$. A constant-gamma EoS is assumed, scattering opacity is neglected, and a Courant factor $C_{a}=0.25$ is used.

Initial conditions are chosen in such a way that the system evolves until it reaches a final stationary state. Neglecting time derivatives, Equations (20)-(25) lead to

[TABLE]

A time-independent solution demands that quantities under derivative in Equations (72)–(74) remain constant, and this condition must also be respected by the initial states. In addition, Equations (75) to (76) show that the final $F_{r}^{x}$ and $P_{r}^{xx}$ must be continuous, although their derivatives can be discontinuous. This does not necessarily imply that the final $E_{r}$ profile must also be continuous, since any value of $P_{r}^{xx}(E_{r},F^{x}_{r})$ can correspond to up to two different $E_{r}$ values for fixed $F^{x}_{r}$. However, in the particular case where $F^{x}_{r}<P_{r}^{xx}$, it can be shown using Eqs. (31)-(33) that the inversion of $P_{r}^{xx}(E_{r},F^{x}_{r})$ in terms of $E_{r}$ leads to unique solutions, and thus $E_{r}$ must be continuous. In the same way, we have verified that this condition is equivalent to $F^{x}_{r}/E_{r}<3/7$.

We have performed four tests for different physical regimes. All the initial values are chosen to coincide with those in Farris et al. (2008). In that work, as in several others where the same tests are performed (see e.g. Zanotti et al., 2011; Fragile et al., 2012; Sadowski et al., 2013), the Eddington approximation, given by $\tilde{P}^{xx}_{r}=\tilde{E}_{r}/3$, is used instead of the M1 closure. Therefore, our results are not comparable with these unless the final state satisfies $\tilde{P}^{xx}_{r}\simeq\tilde{E}_{r}/3$ in the whole domain. We now outline the main features of each test, whose parameters are summarized in Table 1:

Nonrelativistic strong shock. A gas-pressure-dominated shock moves at a nonrelativistic speed in a cold gas ($p_{g}\ll\rho$), with a maximum $u^{x}$ of 0.015. The final profiles of $\rho$, $p_{g}$, $u^{x}$, $\tilde{E}_{r}$, and $\tilde{F}^{x}_{r}$ are shown in Fig. 6. As in the non-radiative case, the first three show an abrupt change at $x=0$, while radiation fields seem continuous. 2. 2.

Mildly relativistic strong shock. The conditions are similar to the previous test, with the difference that a mildly relativistic velocity ($u^{x}\leq 0.25$) is chosen. The final profiles (see Fig. 7) look similar to those in Fig. 6, with the difference that $\tilde{E}_{r}$ exhibits a small discontinuity close to $x=0$. 3. 3.

Highly relativistic wave. Initial conditions are those of a highly relativistic gas-pressure-dominated wave ($u^{x}\leq 10$, $\rho\ll\tilde{P}^{xx}_{r}<p_{g}$). In this case, as it can be seen in Fig. 8, all the profiles are continuous. 4. 4.

Radiation-pressure-dominated wave. In this case we study a situation where the radiation pressure is much higher than the gas pressure, in a shock that propagates at a mildly relativistic velocity ($u^{x}\leq 0.69$). As in the previous case, there are no discontinuities in the final profiles (see Fig. 9).

In order to test the convergence of the numerical solutions, we have performed each simulation twice, dividing the domain in $800$ and in $3200$ zones. In every case, as shown in Figs. 6-9, both solutions coincide. However, our results do not coincide with those obtained in the references mentioned above. The most noticeable case is the test shown in Fig. 7, where the ratio $\tilde{P}^{xx}_{r}/\tilde{E}_{r}$ reaches a maximum value of $0.74$ close to the shock, instead of the value of $1/3$ that would be obtained within the Eddington approximation. The result is a much smoother $\tilde{E}_{r}$ profile than the one shown in, for instance, Farris et al. (2008). Yet, our results show a good agreement with those in Takahashi et al. (2013), where the tests are also performed assuming the M1 closure.

We point out that, in the nonrelativistic strong shock case, characteristic fluid speeds are $\sim 35$ times smaller than those corresponding to radiation transport. Still, computations do not show significant increase of numerical diffusion owing to such scale disparity. The same conclusion holds if computations are done in the downstream reference frame (not shown here).

4.5 Radiation pulse

Following Sadowski et al. (2013), we have tested the evolution of a radiation pulse in the optically thin and optically thick limits. These two regimes allowed us to assess, respectively, the code performance when choosing different coordinate systems and its accuracy in the diffusion limit, as summarized below.

4.5.1 Optically thin case

We considered an initial spherically symmetric radiation energy distribution, contained around the center of a 3D box of side $L=100$. Radiation energy is initially set as $E_{r}=4\pi B(T_{r})$, with

[TABLE]

where $r$ is the spherical radius, while

$T_{0}=10^{6}$ and $w=5$. Similarly, gas pressure is set in such a way that $T(\rho,p_{g})=T_{0}$, which means that the system is initially in thermal equilibrium far from the pulse. We also set $\rho=1$, $v^{x}=0$ and $F_{r}^{x}=0$ in the whole domain, $\Gamma=5/3$, $C_{a}=0.4$, $\kappa=0$, and a small scattering opacity $\sigma=10^{-6}$. In this way, the total optical depth from side to side of the box is $\tau=\rho\,\sigma L=10^{-4}\ll 1$, i.e., the box is transparent to radiation.

We have computed the departure from these conditions using 1D spherical and 3D Cartesian coordinates. In the Cartesian case, we have employed a uniform grid resolution of $200\times 200\times 200$ zones. On the other hand, in spherical geometry, our domain is the region $r\in[0,L/2]$ using a uniformly spaced grid of $100$ zones, in order to have a comparable resolution with the 3D simulations. In this last case, reflective boundary conditions have been set at $r=0$.

As shown in Fig. 10, the pulse expands and forms a nearly isotropic blast wave, which slightly deviates from the spherical shape in the Cartesian case due to grid noise. The evolution of the radiation energy profiles in both simulations is shown in the two upper panels of Figure 11. Since no absorption in the material is considered, the total radiation energy is conserved, and thus the maximum energy density of the formed expanding wave decreases as $1/r^{2}$. As it can be seen in Fig. 11, this dependence is effectively verified once the blast wave is formed. The same kind of analysis is possible if radiation is contained entirely on the plane $z=0$. In this case, the maximum energy density decreases as $1/R$, with $R=\sqrt{x^{2}+y^{2}}$. We have verified this behavior in 1D cylindrical and 2D Cartesian coordinates, employing uniform grids of $100$ zones in the first case and $200\times 200$ in the second (see the two lower panels in Fig. 11). In every case, the same simulations performed with different coordinates systems show a good agreement.

4.5.2 Optically thick case

We now consider the case where the scattering opacity is nine orders of magnitude larger than in the previous simulations, i.e., $\sigma=10^{3}$, and all the other parameters remain unchanged. In that situation, the optical thickness from side to side of the box is $\tau=10^{5}\gg 1$, which means that the box is largely opaque to radiation. Here we solve the evolution equations on a Cartesian one-dimensional grid with uniform spacing. Using a resolution of 101 zones, the optical thickness of a single cell is $\tau\sim 10^{3}$. For this reason, signal speeds are always limited accordingly to Eq. (A6).

Under these conditions, the system evolves in such a way that $|\partial_{t}F_{r}^{x}|\ll|\partial_{x}P_{r}^{xx}|$ and $|F_{r}^{x}|\ll E_{r}$, and therefore $P_{r}^{xx}\simeq E_{r}/3$, as pointed out in Section 2.4. Neglecting the term $\partial_{t}F_{r}^{x}$ in Eq. (25) and assuming $P_{r}^{xx}=E_{r}/3$, the radiation flux can be written as $F_{r}^{x}=-\partial_{x}E_{r}/3\rho\chi$. Hence, assuming the density to remain constant, the radiation energy density should evolve accordingly to the following diffusion equation:

[TABLE]

With the chosen initial conditions, this equation can be solved analitically, e.g., by means of a Fourier transform in the spatial domain. The exact and numerical solution are shown in Fig. 12. Our results show a good agreement between the analytical and numerical solutions. Furthermore, we have verified that, if radiation-matter interaction is not taken into account for the signal speed calculation, i.e., if the limiting given by Eq. (A6) is not applied, the pulse gets damped much faster than what it should be expected from Eq. (78), due to the numerical diffusion that occurs when signal speeds are overestimated.

We have observed that this test leads to inaccurate values of $F_{r}^{x}$ if IMEX-SSP2(2,2,2) is used, although the values of $E_{r}$ remain close to the analytical ones. This problem lies in the fact that both the gradient of the flux of $F^{x}_{r}$ and its source term largely exceed $F^{x}_{r}$ and are not compensated in the last explicit step of the method (see Eq. (36)). When these conditions are met, we have observed that IMEX-SSP2(2,2,2) can lead to inaccuracies and instabilities due to failure in preserving energy positivity (see Section 4.6). On the contrary, IMEX1 shows better performances in those cases, as flux and source terms are more accurately balanced during the implicit steps (see Eq. (37)).

The limiting scheme in Eq. (A6) depends on the optical depth of individual cells, which is inversely proportional to the resolution. Therefore, when AMR is used, there can be situations where this limiting is applied in the coarser levels, but not in the finer ones. Furthermore, when using HLLC, the solver is replaced by HLL for every zone where Eq. (A6) is enforced. To study the code’s performance under these conditions, we have run this test on a static AMR grid using $128$ zones at the coarsest level with 6 levels of refinement with a jump ratio of 2, yielding an equivalent resolution of $8192$ zones. We choose $\sigma=50$ so that levels 0 to 4 are solved with the HLL solver limiting the maximum signal speeds accordingly to Eq. (A6), while levels 5 and 6 are solved using the HLLC solver. The solution thus obtained converges to the analytic solution of Eq. (78) in all the refinement levels (see Fig. 13).

However, we have observed the formation of spurious overshoots at the boundaries between refinement levels. These artifacts are drastically reduced if the order of the reconstruction scheme is increased; for instance, if the weighted essentially non-oscillatory (WENO) method by Jiang & Shu (1996) or the piecewise parabolic method (PPM) by Colella & Woodward (1984) are used, as shown in Fig. 13. We argue that such features, which are not uncommon in AMR codes (Choi et al., 2004; Chilton & Colella, 2010), can be attributed to the refluxing process needed to ensure correct conservation of momentum and total energy (see Mignone et al., 2012). In addition, the presence of source terms requires additional care when solving the Riemann problem betwen fine-coarse grids due to temporal interpolation (Berger & LeVeque, 1998). We do not account here for such modifications and defer these potential issues to future investigations.

4.6 Shadows

One of the main features of the M1 closure is its ability to reproduce the behavior of physical systems in which the angular distribution of the radiation specific intensity has strong spatial variations. One such example is a system where a free-streaming radiation field encounters a highly opaque region of space, casting a shadow behind it. To test the code performance when solving such problems, we have performed a test in which a shadow is formed behind a high-density elliptic cylinder, following Hayes & Norman (2003) and using the same parameters as in González et al. (2007).

Computations are carried out in the two-dimensional domain given by $\left\{(x,y)\in[-0.5,0.5]\,\mbox{cm}\times[0,0.6]\,\mbox{cm}\right\}$. Reflective boundary conditions are imposed at $y=0$. A constant density $\rho_{0}=1$ g cm*-3* is fixed in the whole space, except in the elliptic region, where $\rho=$ $\rho_{1}=10^{3}$ g cm*-3*. In order to have a smooth transition between $\rho_{0}$ and $\rho_{1}$, the initial density field is defined as

[TABLE]

where

[TABLE]

with $(x_{0},y_{0})=(0.10,0.06)$ cm. In such a way, the region with $\rho=\rho_{1}$ is approximately contained in an ellipse of semiaxes $(x_{0},y_{0})$. Initially, matter is set in thermal equilibrium with radiation at a temperature $T_{0}=290$ K, and fluxes and velocities are initialized to zero. The absorption opacity in the material is computed according to Kramers’ law, i.e., $\kappa=\kappa_{0}\left(\frac{\rho}{\rho_{0}}\right)\left(\frac{T}{T_{0}}\right)^{-3.5}$, with $\kappa_{0}=0.1$ g*-1*cm2, while scattering is neglected. Therefore, the cylinder’s optical thickness along its largest diameter is approximately $\tau\approx 2\times 10^{4}$, which means that its width exceeds $\tau\gg 1$ times the photons’ mean free path in that region. On the contrary, above $y>y_{0}$, the optical thickness is $\tau=0.1$, so that the exterior of the cylinder is transparent to radiation while its interior is opaque.

Radiation is injected from the left boundary at a temperature $\left(c\,E_{r}/4\,\sigma_{SB}\right)^{1/4}=1740$ K $>T_{0}$, with a flux $\mathbf{F}_{r}=c\,E_{r}\,\hat{\bm{e}}_{x}$. Hence, the radiation field is in initially in the free-streaming limit, and should be transported at the speed of light in the transparent regions.

We have initially computed the system’s evolution in a fixed uniform grid of resolution $280\times 80$, using a fourth-order reconstruction scheme with a Courant factor $C_{a}=0.4$, and with $\Gamma=5/3$. Simulations show a radiation front that crosses the space at light speed from left to right, producing a shadow behind the cylinder. After interacting with it, the shadow settles into a final stable state that is ideally maintained until the matter distribution is modified due to its interaction with radiation. The radiation energy density distribution is shown in the upper panel of Fig. 14 at $t=10\,t_{c}$, where $t_{c}=1\,\mbox{cm}/c=3.336\times 10^{-11}$ s is the light-crossing time, namely, the time it takes light to cross the domain horizontally in the transparent region. Behind the cylinder, radiation energy is roughly equal to its initial equilibrium value of $(4\,\sigma_{SB}/c)\,T_{0}^{4}$. This value is slightly affected by small waves that are produced in the upper regions of the cylinder, where the matter distribution stops being opaque to radiation along horizontal lines. Above the cylinder, the radiation field remains equal to the injected one. The transition between the shadowed and transparent regions is abrupt, as it can be seen in Fig. 14. The shape of the $E_{r}$ profile along vertical cuts is roughly maintained as radiation is transported away from the central object.

When IMEX-SSP2(2,2,2) is used, we have noticed that $E_{r}$ goes frequently below [math] on the left edge of the cylinder where the radiation field impacts it. Still, the obtained solutions are stable and convergent as long as $E_{r}$ is floored to a small value whenever this occurs. As in Section 4.5.2, the radiation flux is much smaller in those zones than both its flux and the source terms, and the problem does not occur if IMEX1 is used.

We have used this same problem to test the code’s performance when AMR is used in a multidimensional setup. In this case, we have run the same simulation, using an initially coarse grid of resolution $80\times 16$ set to adapt to changes in $E_{r}$ and $\rho$ (see Mignone et al., 2012). We have used $5$ refinement levels, in every case with a grid jump ratio of $2$, which gives an equivalent resolution of $2560\times 512$. The resulting energy profiles are plotted in the lower panels of Figure 14, for $t=0.2\,t_{c}$, $0.6\,t_{c}$, and $10\,t_{c}$, and agree with those computed using a fixed grid. In each panel we have superimposed the refinement level.

4.7 Magnetized cylindrical blast wave

We now examine a case in which matter is affected by both radiation and large-scale EM fields. We consider the case of a cylindrical blast wave, represented in a two-dimensional Cartesian grid as in the planar configurations described in Section 4.5.1. In the context of MHD, this kind of problem has been used formerly to check the robustness of the employed methods when handling relativistic magnetized shocks, as well as their ability to deal with different kinds of degeneracies (see e.g. Komissarov, 1999; Mignone & Bodo, 2006). In our case, we draw on this configuration as an example system that can switch from radiation-dominated to magnetically dominated regimes, depending on the material’s opacity. To this end, we set up a cylindrical explosion from an area where the magnetic pressure is of the same order of the gas pressure, and both are smaller than the radiation pressure. Under this condition, matter dynamics is magnetically dominated when the opacities are low, and radiation-dominated in the opposite case. The latter case also serves to investigate the high-absorption regime in which both the diffusion approximation and LTE are valid.

We consider a square domain defined as $(x,y)\in[-6,6]\times[-6,6]$, initially threaded by a uniform magnetic field, $\mathbf{B}=B_{0}\,\hat{\bm{e}}_{x}$ with $B_{0}=0.1$. Gas pressure and density are initially set as follows:

[TABLE]

where $p_{0}=3.49\times 10^{-5}$, $\rho_{0}=10^{-4}$ are the ambient values while $p_{1}=1.31\times 10^{-2}$, $\rho_{1}=10^{-2}$ identify the over-pressurized region. Here $R=\sqrt{x^{2}+y^{2}}$ is the cylindrical radius while $\delta\equiv\delta(R/R_{0})$ is a taper function that decreases linearly for $R_{0}<R\leq 1$ (we use $R_{0}=0.8$). The ideal equation of state with $\Gamma=4/3$ is used throughout the computations.

A radiation field is introduced initially in equilibrium with the gas. Since $\mathbf{v}=\mathbf{0}$ in the whole domain, the condition of LTE is initially satisfied if $E_{r}=4\pi B(T)$ and $\mathbf{F}_{r}=\mathbf{0}$. These conditions are chosen in such a way that, close to the center of the domain, $p_{g}\sim\mathbf{B}^{2}/2<E_{r}/3$, where $\mathbf{B}^{2}/2$ gives the initial magnetic contribution to the total pressure (see Eq. (26)). To guarantee the condition $\nabla\cdot\mathbf{B}=0$, necessary for the solutions’ stability, we have implemented in every case the constrained transport method.

Figure 15 shows a set of 2D color maps representing the fields’ evolution at $t=4$, using a resolution of $360\times 360$ zones. The two upper rows correspond to computations using $\sigma=0$ and $\kappa=1$ (top) or $1000$ (middle). For $\kappa=1$, the initial optical depth along the central sphere is $\tau\approx\rho_{1}\kappa\Delta x=0.02\ll 1$, and therefore the material’s expansion should not be noticeably affected by the radiation field. Indeed, in this case, the radiation energy profile expands spherically as in Section 4.5.1. The dynamic is magnetically dominated and matter is accelerated up to $\gamma\sim 1.7$ along the magnetic field lines along the $x$ axis, which is why the hydrodynamic variables are characterized by an elongated horizontal shape.

The second row of Fig. 15 shows analog results obtained with $\kappa=1000$, where $\tau\approx 20\gg 1$. In this case, the interaction of the radiation field with the gas during its expansion produces a much more isotropic acceleration. This acceleration is still maximal along the $x$ direction, due to the joint contributions of the magnetic field and the radiation pressure. This is why the Lorentz factor is larger in this case, reaching $\gamma\sim 2.7$. Gas density and pressure reach their maxima along an oblated ring, instead of the previous elongated distributions obtained with $\kappa=1$. As shown in the same figures, the magnetic field lines are pushed away from the center as matter is radially accelerated, producing a region of high magnetic energy density around the area where $\gamma$ is the highest, and a void of lower magnetic energy inside. Also differently from the previous case, the radiation energy distribution is no longer spherically symmetric due to its interaction with the matter distribution.

For high values of $\rho\kappa$, it is expected that the radiation reaches LTE with matter, as Eqs. (18), (24) and (25) lead to $\tilde{E}_{r}\rightarrow 4\pi B(T)$ and $\tilde{F}^{i}_{r}\rightarrow 0$ for smooth field distributions that do not vary abruptly in time. In this limit, Eqs. (20)-(25) can be reduced to those of relativistic MHD, redefining the total gas pressure as

[TABLE]

with $p_{g}^{\prime}=p_{g}+\tilde{p}_{r}$, and the enthalpy density as

[TABLE]

where $\tilde{P}_{r}^{ij}=\tilde{p}_{r}\,\delta^{ij}$, which follows from the M1 closure in this limit. Taking a constant-$\Gamma$ EoS with $\Gamma=4/3$ in every case, the equations of state of both systems of equations coincide in the large-opacity limit, and therefore the results obtained with both of them are comparable.

The third row of Fig. 15 shows the results of an ideal relativistic MHD simulation performed in such a way, using the same initial conditions as before. To compute the gas pressure represented herein, it was assumed that $p_{g}^{\prime}\simeq\tilde{p}_{r}=4\pi B(T)/3$, from where it is possible to extract $T$ and then $p_{g}$. Following the same idea, an effective $E_{r}$ was computed boosting its comoving value, assumed to be equal to $4\pi B(T)$, and taking $\tilde{F}^{i}_{r}=0$. The resulting plots thus obtained are in fact similar to those computed with $\kappa=1000$, with slight differences that can be explained taking into account that $\kappa$ has a finite value, and that close to the shocks the fields do not satisfy the conditions of being smooth and varying slowly with time. Consequently, the value of $\tilde{E}_{r}$ can be different than $4\pi B(T)$ in the regions of space that are close to discontinuities, which means that the hypothesis of LTE, assumed by ideal MHD, is not satisfied in the whole domain. This is verified in Figure 16, where it is shown that, for $\kappa=1000$, the ratio $\tilde{E}_{r}/4\pi B(T)$ differs from $1$ only in the regions that are close to shocks, shown in Fig. 15.

4.8 Parallel Performance

Parallel scalability of our algorithm has been investigated in strong scaling through two- and three-dimensional computations. For simplicity, we have chosen the (unmagnetized) blast wave problem of Section 4.7 with $\kappa=10$ leaving the remaining parameters unchanged. For reference, computations have been carried out with and without the radiation field on a fixed grid of $2304^{2}$ (in 2D) and $288^{3}$ (in 3D) zones, a constant time step and the solver given by Eq. (39). The number of processors - Intel Xeon Phi7250 (KnightLandings) at 1.4 GHz - has been varied from $N_{\rm CPU}=8$ to $N_{\rm CPU}=1024$.

The corresponding speed-up factors are plotted in Fig. 17 as a function of NCPU (solid lines with symbols) together with the ideal scaling-law curve $\propto N_{\rm CPU}$. We compute the speedup factors as $S=T_{\rm ref}/T_{\mathrm{N}_{\mathrm{CPU}}}$ where $T_{\rm ref}$ is a normalization constant while $T_{\mathrm{N}_{\mathrm{CPU}}}$ is the total running time for a simulation using $N_{\mathrm{CPU}}$ processors.

Overall, favorable scaling properties are observed in two and three dimensions with efficiencies that remain above $90\%$ up to 256 cores and drops to $\sim 70\%$ when $N_{\rm CPU}=1024$. Slighlty better results are achieved when radiation is included, owing to the additional computational overhead introduced by the implicit part of the algorithm which uses exclusively local data without requiring additional communication between threads.

Note that, for convenience, we have normalized the curves to the corresponding running time without the radiation field. This demonstrates that, by including radiation, the code is (approximately) four times more expensive than its purely hydro counterpart, regardless of the dimensionality.

5 Summary

We have presented a relativistic radiation transfer code, designed to function within the PLUTO code. Our implementation can be used together with the relativistic HD and MHD modules of PLUTO to solve the equations of radiation transfer under the gray approximation. Integration is achieved through one of two possible IMEX schemes, in which source terms due to radiation-matter interaction are integrated implicitly and flux divergences, as well as every other source term, are integrated explicitly. The transition between optically thick and thin regimes is controlled by imposing the M1 closure to the radiation fields, which allows to handle both the diffusion and free-streaming limits. Opacity coefficients can be arbitrarily defined, depending on problem at hand, as functions of the primitive variables.

In our implementation, a novel HLLC Riemann solver for radiation transport has been introduced. The new solver is designed to improve the accuracy of the solutions with respect to it predecessors (such as HLL) in optically thin regions of space. The module has been designed to function with either Cartesian, cylindrical or spherical coordinates in multiple spatial dimensions and it is suitable for either serial or parallel computations. Extension to adaptive grids, based on the standard implementation of the CHOMBO library within the code, has also been presented.

We have performed a series of numerical benchmarks to assess the module performance under different configurations, including handling of radiation transport, absorption, and emission in systems with different characteristics. Our results demonstrate excellent stability properties under the chosen parameters, in both the free-streaming and diffusion limits. In the latter case, numerical diffusion is successfully controlled by limiting the signal speeds of the radiation transport equations whenever the material is opaque across single grid cells. Overall, the transition between both regimes has been properly captured by the code in all the considered cases. For optically thin transport, our HLLC solver produces more accurate solutions when compared to HLL. Regarding the implemented IMEX schemes, we have seen a similar performance of both IMEX-SSP2(2,2,2) and IMEX1 except in tests where the order of magnitude of the radiation flux is much smaller than both its source terms and the divergence of its own flux, in which IMEX1 seems to have better stability and positivity-preserving properties. When AMR is used, the obtained solutions exhibit a similar overall behavior to those computed using a fixed grid. Good agreement is also shown with standard tests whenever the comparison is possible. Furthermore, parallel performance tests show favorable scaling properties which are comparable to those of the RHD module of PLUTO.

The code presented in this work will be made publicly available as part of future versions of PLUTO, which can currently be downloaded from http://plutocode.ph.unito.it/.

The authors would like to thank the referee for his/her constructive comments which helped to improve the quality of this manuscript. Also, we acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support.

Appendix A Signal speeds

We describe in this section the computation of the characteristic signal speeds used in the explicit step (see Section 3.2). In the particular form of Equations (20)-(25), the MHD fluxes are independent of the radiation variables, and vice-versa. Hence, the Jacobian matrices of the system are block-diagonal, one block corresponding to MHD and the other to radiation transport. Consequently, their corresponding sets of eigenvalues can be obtained by computing the eigenvalues of each of these blocks individually.

For the MHD block, maximum and minimum signal speeds are computed as detailed in Mignone & Bodo (2006) and Mignone et al. (2007b). The treatment needed for the radiation wave speeds is rather simpler, as a short calculation shows that, for every direction $d$, the radiation block depends only on the angle $\theta$ between $\mathbf{F}_{r}$ and $\hat{\bm{e}}_{d}$, and on $f=||\mathbf{F}_{r}||/E_{r}$. This simplifies the calculation, which can be performed analitically as shown in Audit et al. (2002) and Skinner et al. (2013). The full set of eigenvalues of the radiation block, which we denote as $\{\lambda_{r1},\lambda_{r2},\lambda_{r3}\}$, can be computed as

[TABLE]

[TABLE]

[TABLE]

where $\xi(f)$ is defined in Eq. (33), while

[TABLE]

When $f=0$, $\lambda_{r2}$ is replaced by 0, i.e., its limit when $f\rightarrow 0$. It can be seen from Equations (A1)-(A4) that the following inequalities hold for every value of $f$ and $\theta$:

[TABLE]

In the free-streaming limit ($f=1$), all these eigenvalues coincide and are equal to $\cos\theta$, which gives $\lambda_{rj}=\pm 1$ in the parallel direction to $\mathbf{F}_{r}$, and $\lambda_{rj}=0$ in the perpendicular ones for $j=1,2,3$. On the other hand, in the diffusion limit ($f=0$), we have $(\lambda_{r1},\lambda_{r2},\lambda_{r3})=(-1/\sqrt{3},0,1/\sqrt{3})$ in every direction.

The above analysis can be applied to homogeneous hyperbolic systems. Although the equations of Rad-RMHD do not belong to this category, this is not a problem when radiation transport dominates over radiation-matter interaction. On the contrary, in the diffusion limit, the moduli of the maximum and minimum speeds, both equal to $1/\sqrt{3}$, may be too big and lead to an excessive numerical diffusion. In those cases, the interaction terms need to be taken into account to estimate the wave speeds. With this purpose, following Sadowski et al. (2013), we include in the code the option of locally limiting the maximum and minimum speeds by means of the following transformations:

[TABLE]

where $\tau=\rho\,\gamma\,\chi\,\Delta x$ is the optical depth along one cell, being $\Delta x$ its width in the current direction. Hence, this limiting is only applied whenever cells are optically thick. The reduced speeds in Eq. (A6) are based on a diffusion equation like Eq. (78), where the diffusion coefficient is $1/3\rho\chi$.

Appendix B Semi-analytical proof of $\lambda_{L}\leq\lambda^{*}\leq\lambda_{R}$

In order to check the validity of Equation (46), we have verified the following relations:

[TABLE]

As in Section 3.3, we omit the subindex $r$, as it is understood that only radiation fields are here considered.

We begin by proving the positivity of $A_{R}$. From its definition in Equation (53), we have:

[TABLE]

where $\lambda_{3,L/R}=\lambda_{3}(f_{L/R},\theta_{L/R})$. Since $E>0$, we can conclude from Eq. (B3) that $A_{R}\geq 0$ is automatically satisfied if

[TABLE]

From Eq. (A3), this condition can be rewritten as

[TABLE]

where $\Delta=\sqrt{4-3f^{2}}$. Taking squares at both sides and rearranging terms, this condition reads

[TABLE]

where $X(f)=\frac{2}{3}(4-3f^{2}-\Delta)$ and $Y(f)=(1-f^{2})(4-3f^{2}-2\Delta)$. Since only the second of these terms can be smaller than [math], it is enough to prove that (B6) holds for $\cos^{2}\theta=1$, since that yields the minimum value that the left-hand side can take when $Y<0$. Hence, it is enough to prove

[TABLE]

which holds since the last term between parentheses is always greater than or equal to [math]. This finishes the proof of Eq. (B4). Using the same equations, we can see that $\lambda_{3}(f,\theta)-f_{x}=0$ is only satisfied if $f=1$. An analog treatment can be used for $A_{L}$, from which we arrive to the following inequalities:

[TABLE]

We now proceed to verify Equations (B1) and (B2), firstly considering the case $f_{L},f_{R}<1$, in which $A_{L/R}\neq 0$. Under this condition, the ratio $B_{S}/A_{S}$ depends only on $(f_{L},f_{R},\theta_{L},\theta_{R})$ as

[TABLE]

with $S=L,R$. In order to verify Eq. (B1), we must prove $\lambda_{R}\geq B_{R}/A_{R}$ and $\lambda_{R}\geq B_{L}/A_{L}$. Since $\lambda_{R}=\max(\lambda_{3,L},\lambda_{3,R})$, we can write the first of these conditions considering the cases $\lambda_{R}=\lambda_{3,R}$ and $\lambda_{R}=\lambda_{3,L}$, as

[TABLE]

The first of these can be verified from the graph of $\lambda_{3}(f,\theta)-\alpha(\lambda_{3}(f,\theta),f,\theta)$, where it can be seen that this function is always positive and tends to [math] for $f\rightarrow 1$. Similarly, we have checked the second one numerically by plotting $\lambda_{3,L}-\alpha(\lambda_{3,L},f_{R},\theta_{R})$ under the condition $\lambda_{3,R}<\lambda_{3,L}$, taking multiple values of $\lambda_{3,L}$ covering the range $[-1,1]$. The condition $\lambda_{R}\geq B_{L}/A_{L}$ can be proven in a similar fashion, by considering the cases $\lambda_{L}=\lambda_{1,L}$ and $\lambda_{L}=\lambda_{1,R}$. Since $\lambda_{R}\geq\lambda_{3,L}$, it is enough to prove the following conditions:

[TABLE]

which can be verified in the same manner, finishing the proof of Eq. (B1) for the case $f_{L},f_{R}<1$. The same procedure can be used to prove the validity of Eq. (B2).

Unlike the RHD case, the maximum and minimum eigenvalues do not satisfy $\lambda_{L}<0$ and $\lambda_{R}>0$. However, studying the parabolae defined at both sides of Eq. (52), it can be shown that $\lambda^{*}$ is always contained between $B_{R}/A_{R}$ and $B_{L}/A_{L}$, regardless of the order of these two values and of the signs of $\lambda_{L}$ and $\lambda_{R}$. Hence,

[TABLE]

Together with relations (B1) and (B2), this proves Eq. (46) for $f_{L},f_{R}<1$. These results are also valid in the cases $f_{L}=1$ and $f_{R}=1$ whenever the $A$ functions differ from 0. Let us now assume $f_{a}=1$ and $f_{b}\neq 1$. From Eqs. (53) and (54), we have $A_{a}\cos\theta_{a}=B_{a}$ and consequently $A_{a}=0$ implies that $B_{a}=0$. If $A_{a}=0$ and $A_{b}\neq 0$, from (52) we can extract that $\lambda^{*}=B_{b}/A_{b}$. Finally, from the above analysis, we know that $\lambda_{L}\leq B_{b}/A_{b}\leq\lambda_{R}$, from which we conclude that (46) holds even in this case. The only remaining case is that in which $f_{L}=f_{R}=1$ and $A_{L}=A_{R}=0$, already considered in Section 3.3, where the HLLC solver is replaced by the usual HLL solver.