Non-LTE radiation hydrodynamics in PLUTO

TL;DR

This paper introduces an extended non-LTE radiation hydrodynamics module for the PLUTO code, enabling more accurate modeling of astrophysical phenomena where LTE assumptions are invalid, validated through shock tests and error analysis.

Contribution

The authors developed and implemented a non-LTE radiation hydrodynamics module in PLUTO, extending previous LTE-based models to handle non-LTE regimes with validated accuracy.

Findings

Good agreement with semi-analytical solutions, maximum error of 7%

Non-LTE approach is crucial for accurate accretion shock modeling

The module is validated through radiative shock tests

Abstract

Modeling the dynamics of most astrophysical structures requires an adequate description of the radiation-matter interaction. Several numerical (magneto)hydrodynamics codes were upgraded with a radiation module to fulfill this request. However, those among them that use either the flux-limited diffusion (FLD) or the M1 radiation moment approaches are restricted to the local thermodynamic equilibrium (LTE). This assumption may be not valid in some astrophysical cases. We present an upgraded version of the LTE radiation-hydrodynamics module implemented in the PLUTO code, originally developed by Kolb et al. (2013), which we have extended to handle non-LTE regimes. Starting from the general frequency-integrated comoving-frame equations of radiation hydrodynamics (RHD), we have justified all the assumptions made to obtain the non-LTE equations actually implemented in the module, under the FLD approximation. An operator-split method is employed, with two substeps: the hydrodynamic part is solved with an explicit method by the solvers already available in PLUTO, the non-LTE radiation diffusion and energy exchange part is solved with an implicit method. The module is implemented in the PLUTO environment. It uses databases of radiative quantities that can be provided independently by the user: the radiative power loss, the Planck and Rosseland mean opacities. Our implementation has been validated through different tests, in particular radiative shock tests. The agreement with the semi-analytical solutions (when available) is good, with a maximum error of 7%. Moreover, we have proved that non-LTE approach is of paramount importance to properly model accretion shock structures. Our radiation FLD module represents a step toward the general non-LTE RHD modeling. The module is available, under request, for the community.

Tables

Table 1. Table 1: Relative sizes of terms in the comoving-frame radiation energy equation ( 17d ) ( erg ​ cm − 3 ​ s − 1 ) erg superscript cm 3 superscript s 1 \left({\rm erg}\;{\rm cm^{-3}}\;{\rm s^{-1}}\right)

Radiation energy		$\frac{\partial E_0}{\partial t}$	$\nabla \cdot {𝑭}_0$	$\text{varv}\cdot \frac{\partial }{\partial t}\left(\frac{{𝑭}_0}{c^2}\right)$	$\nabla \cdot \left(E_0\text{varv}\right)$	$P_0:\nabla \text{varv}$	$2{𝑭}_0\cdot \frac{𝒂}{c^2}$	$-c{G}_0^0$	Normalization
(tag)		(a)	(b)	(c)	(d)	(e)	(f)	(g)
free streaming	$t_r\sim \frac{\ell }{c}$	1	1	$\beta$	$\beta$	$\beta$	${\beta }^2$	$L_0\frac{t_r}{E_0}\left(\underset{LTE}{\to }\tau \right)$	$\times \left(\frac{E_0}{t_r}\right)$
	rank	1	1	2	2	2	3	1
static diffusion	$t_d\sim \frac{{\ell }^2}{c{\lambda }_p}$	1	1	$\frac{\beta }{\tau }$	$\beta \tau$	$\beta \tau$	${\beta }^2$	1	$\times \left(\frac{E_0}{t_d}\right)$
$t_d\ll t_f$	rank	1	1	3	2	2	4	1
dynamic diffusion	$t_f\sim \frac{\ell }{\text{varv}}$	1	$\frac{1}{\beta \tau }$	$\frac{\beta }{\tau }$	1	1	$\frac{\beta }{\tau }$	1	$\times \left(\frac{E_0}{t_f}\right)$
$t_f\ll t_d$	rank	1	2	3	1	1	3	1

Table 2. Table 2: Relative sizes of terms in the comoving-frame radiation momentum equation ( 17e ) ( dyn ​ cm − 3 ) dyn superscript cm 3 \left({\rm dyn}\;{\rm cm^{-3}}\right)

Radiation momentum		$\frac{\partial }{\partial t}\left(\frac{{𝑭}_0}{c^2}\right)$	$\nabla \cdot P_0$	$\left(\text{varv}\cdot \nabla \right)\left(\frac{{𝑭}_0}{c^2}\right)$	$\left(\frac{{𝑭}_0}{c^2}\cdot \nabla \right)\text{varv}$	$\left(\nabla \cdot \text{varv}\right)\left(\frac{{𝑭}_0}{c^2}\right)$	$\frac{\partial P_0}{\partial t}\frac{\text{varv}}{c^2}$	$\frac{1}{c^2}\left(P_0𝒂\right)$	$\frac{1}{c^2}\left(E_0𝒂\right)$	$-{𝑮}_0$	Normalization
(tag)		(a)	(b)	(c)	(d)	(e)	(f)	(g)	(h)	(i)
free streaming	$t_r\sim \frac{\ell }{c}$	1	1	$\beta$	$\beta$	$\beta$	$\beta$	${\beta }^2$	${\beta }^2$	$\tau$	$\times \left(\frac{F_0}{c\ell }\right)$
	rank	1	1	2	2	2	2	3	3	¿ 1
static diffusion	$t_d\sim \frac{{\ell }^2}{c{\lambda }_p}$	$\frac{1}{{\tau }^2}$	1	$\frac{\beta }{\tau }$	$\frac{\beta }{\tau }$	$\frac{\beta }{\tau }$	$\frac{\beta }{\tau }$	${\beta }^2$	${\beta }^2$	1	$\times \left(\frac{F_0}{c{\lambda }_p}\right)$
$t_d\ll t_f$	rank	2	1	3	3	3	3	4	4	1
dynamic diffusion	$t_f\sim \frac{\ell }{\text{varv}}$	$\frac{\beta }{\tau }$	1	$\frac{\beta }{\tau }$	$\frac{\beta }{\tau }$	$\frac{\beta }{\tau }$	${\beta }^2$	${\beta }^2$	${\beta }^2$	1	$\times \left(\frac{F_0}{c{\lambda }_p}\right)$
$t_f\ll t_d$	rank	3	1	3	3	3	2	2	2	1

Table 3. Table 3: Comparisons of the characteristic temperatures of a subcritical shock, T 2 subscript 𝑇 2 T_{2} , T − subscript 𝑇 T_{-}  , T + subscript 𝑇 T_{+}   (Eq.  44 a, b, c), obtained from analytical estimates, from our non-LTE code, and from the LTE initial code by Kolb et al. ( 2013 ) . The last column shows the relative deviation between our numerical solution and the analytical estimate.

	Analytical	Our numerical	Kolb	Deviation
	estimate	solution	et al. (2013)
$T_2$	$\approx$ 865 K	817 K	816.6 K	$\approx 5\%$
$T_{-}$	$\approx$ 315 K	332 K	331.9 K	$\approx 5\%$
$T_{+}$	$\approx$ 1075 K	1151 K	1147.1 K	$\approx 7\%$

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \text{varv}\right)=0,$$

$$\frac{\partial }{\partial t}\left(\rho \text{varv}\right)+\nabla \cdot \left(\rho \text{varv}\otimes \text{varv}+p\mathbb{I}\right)=\rho \mathbf{g}+{\mathbf{G}}_0,$$

$$\frac{\partial }{\partial t}\left(\mathcal{E}+\frac{1}{2}\rho {\text{varv}}^2\right)+\nabla \cdot \left[\left(\mathcal{E}+\frac{1}{2}\rho {\text{varv}}^2+p\right)\text{varv}\right]=\rho \mathbf{g}\cdot \text{varv}$$

$$+{\mathbf{G}}_0\cdot \text{varv}+c{G}_0^0.$$

$$\{\begin{array}{rl}& \frac{\partial E_0}{\partial t}+\nabla \cdot {\mathbf{F}}_0+\text{varv}\cdot \frac{\partial }{\partial t}\left(\frac{{\mathbf{F}}_0}{c^2}\right)+\nabla \cdot \left(E_0\text{varv}\right)+{\mathbb{P}}_0:\nabla \text{varv}\\ & +2\left(\frac{{\mathbf{F}}_0}{c^2}\right)\cdot \mathbf{a}=-c{G}_0^0,\end{array}$$

$$\{\begin{array}{rl}& \frac{\partial }{\partial t}\left(\frac{{\mathbf{F}}_0}{c^2}\right)+\nabla \cdot {\mathbb{P}}_0+\left(\text{varv}\cdot \nabla \right)\left(\frac{{\mathbf{F}}_0}{c^2}\right)+\left(\frac{{\mathbf{F}}_0}{c^2}\cdot \nabla \right)\text{varv}\\ & +\left(\nabla \cdot \text{varv}\right)\left(\frac{{\mathbf{F}}_0}{c^2}\right)+\frac{\partial {\mathbb{P}}_0}{\partial t}\frac{\text{varv}}{c^2}+\frac{1}{c^2}\left({\mathbb{P}}_0\mathbf{a}+E_0\mathbf{a}\right)=-{\mathbf{G}}_0,\end{array}$$

$$G_0^0\left(\mathbf{r},t\right)\equiv \frac{1}{c}{\int }_0^{\infty }d{\nu }_0{\oint }_{4\pi }\left({\chi }_0{I}_0-{\eta }_0\right)d{\Omega }_0,$$

$${\mathbf{G}}_0\left(\mathbf{r},t\right)\equiv \frac{1}{c}{\int }_0^{\infty }d{\nu }_0{\oint }_{4\pi }\left({\chi }_0{I}_0-{\eta }_0\right){\mathbf{n}}_{0}d{\Omega }_0,$$

$${\chi }_0\left(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0\right)={\kappa }_0^{\mathrm{th}}\left(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0\right)+{\kappa }_0^{\mathrm{sca}}\left(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0\right),$$

$${\eta }_0\left(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0\right)={\eta }_0^{\mathrm{th}}\left(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0\right)+{\eta }_0^{\mathrm{sca}}\left(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0\right),$$

$$E_0\left(\mathbf{r},t\right)\equiv \frac{1}{c}{\int }_0^{\infty }{\oint }_{4\pi }{I}_0\left(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0\right)d{\Omega }_{0}d{\nu }_0,$$

$${\mathbf{F}}_0\left(\mathbf{r},t\right)\equiv {\int }_0^{\infty }{\oint }_{4\pi }{I}_0\left(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0\right){\mathbf{n}}_{0}d{\Omega }_{0}d{\nu }_0,$$

$${\mathbb{P}}_0\left(\mathbf{r},t\right)\equiv \frac{1}{c}{\int }_0^{\infty }{\oint }_{4\pi }{I}_0\left(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0\right){\mathbf{n}}_0\otimes {\mathbf{n}}_{0}d{\Omega }_{0}d{\nu }_0.$$

$${\eta }_0^{\mathrm{sca}}(\mathbf{r},t;{\nu }_0)={\kappa }_0^{\mathrm{sca}}(\mathbf{r},t;{\nu }_0){\oint }_{4\pi }{I}_0(\mathbf{r},t;{\mathbf{n}}_0^{\prime },{\nu }_0)\frac{d{\Omega }_0^{\prime }}{4\pi }$$

$$c{G}_0^0=c{\kappa }_{0E}{E}_0-L_0$$

$$L_0\left(\mathbf{r},t\right)\equiv {\int }_0^{\infty }d{\nu }_0{\oint }_{4\pi }{\eta }_0^{\mathrm{th}}(\mathbf{r},t;{\mathbf{n}}_0,{\nu }_0)d{\Omega }_0$$

$${\kappa }_{0E}\equiv \frac{{\int }_0^{\infty }{\kappa }_{{\nu }_0}^{\mathrm{th}}{E}_{{\nu }_0}d{\nu }_0}{{\int }_0^{\infty }{E}_{{\nu }_0}d{\nu }_0}$$

$${\mathbf{G}}_0=\frac{1}{c}{\text{bbchi}}_{0F}{\mathbf{F}}_0$$

$${\text{bbchi}}_{0F}{\mathbf{F}}_0={\int }_0^{\infty }{\chi }_{{\nu }_0}{\mathbf{F}}_0({\nu }_0)d{\nu }_0$$

$${\text{bbchi}}_{0F}=\left(\begin{array}{ccc}{\chi }_{0F}^{(1)}& 0& 0\\ 0& {\chi }_{0F}^{(2)}& 0\\ 0& 0& {\chi }_{0F}^{(3)}\end{array}\right)$$

$$E_{{\nu }_0}$$

$${\mathbf{F}}_{{\nu }_0}$$

$${\kappa }_{0E}={\kappa }_{0P}$$

$${\kappa }_{0P}\equiv \frac{{\int }_0^{\infty }{\kappa }_{{\nu }_0}^{\mathrm{th}}{B}_{{\nu }_0}d{\nu }_0}{{\int }_0^{\infty }{B}_{{\nu }_0}d{\nu }_0}$$

$${\mathbf{F}}_0=-\frac{4\pi }{3{\chi }_{0R}}\nabla B_0$$

$$\frac{1}{{\chi }_{0R}}\equiv \frac{{\int }_0^{\infty }\frac{1}{{\chi }_{{\nu }_0}}\frac{\partial B_{{\nu }_0}}{\partial T}d{\nu }_0}{{\int }_0^{\infty }\frac{\partial B_{{\nu }_0}}{\partial T}d{\nu }_0}$$

$$c{G}_0^0$$

$${\mathbf{G}}_0$$

$$\{\begin{array}{rl}& \frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \text{varv}\right)=0,\end{array}$$

$$\{\begin{array}{rl}& \frac{\partial }{\partial t}\left(\rho \text{varv}\right)+\nabla \cdot \left(\rho \text{varv}\otimes \text{varv}+p\mathbb{I}\right)=\rho \mathbf{g}+\frac{1}{c}{\chi }_{0R}{\mathbf{F}}_0,\end{array}$$

$$\{\begin{array}{rl}& \frac{\partial }{\partial t}\left(\mathcal{E}+\frac{1}{2}\rho {\text{varv}}^2\right)+\nabla \cdot \left[\left(\mathcal{E}+\frac{1}{2}\rho {\text{varv}}^2+p\right)\text{varv}\right]=\rho \mathbf{g}\cdot \text{varv}\\ & +\frac{1}{c}{\chi }_{0R}{\mathbf{F}}_0\cdot \text{varv}+c{\kappa }_{0P}{E}_0-L_0,\end{array}$$

$$\{\begin{array}{rl}& \frac{\partial E_0}{\partial t}+\nabla \cdot {\mathbf{F}}_0+\text{varv}\cdot \frac{\partial }{\partial t}\left(\frac{{\mathbf{F}}_0}{c^2}\right)+\nabla \cdot \left(E_0\text{varv}\right)+{\mathbb{P}}_0:\nabla \text{varv}\\ & +2\left(\frac{{\mathbf{F}}_0}{c^2}\right)\cdot \mathbf{a}=L_0-c{\kappa }_{0P}{E}_0,\end{array}$$

$$\{\begin{array}{rl}& \frac{\partial }{\partial t}\left(\frac{{\mathbf{F}}_0}{c^2}\right)+\nabla \cdot {\mathbb{P}}_0+\left(\text{varv}\cdot \nabla \right)\left(\frac{{\mathbf{F}}_0}{c^2}\right)+\left(\frac{{\mathbf{F}}_0}{c^2}\cdot \nabla \right)\text{varv}\\ & +\left(\nabla \cdot \text{varv}\right)\left(\frac{{\mathbf{F}}_0}{c^2}\right)+\frac{\partial {\mathbb{P}}_0}{\partial t}\frac{\text{varv}}{c^2}+\frac{1}{c^2}\left({\mathbb{P}}_0\mathbf{a}+E_0\mathbf{a}\right)\\ & =-\frac{1}{c}{\chi }_{0R}{\mathbf{F}}_0.\end{array}$$

$$\tau \ll 1$$

$$\tau \gg 1$$

$$\beta$$

$${\mathbf{F}}_0$$

$${\mathbb{P}}_0$$

$${\mathbf{F}}_0$$

$${\mathbb{P}}_0$$

Full text

*Aims. *We present an upgraded version of the LTE radiation-hydrodynamics module implemented in the PLUTO code, originally developed by Kolb et al. (2013), which we have extended to handle non-LTE regimes.

*Methods. *Starting from the general frequency-integrated comoving-frame equations of radiation hydrodynamics (RHD), we have justified all the assumptions made to obtain the non-LTE equations actually implemented in the module, under the FLD approximation. An operator-split method is employed, with two substeps: the hydrodynamic part is solved with an explicit method by the solvers already available in PLUTO, the non-LTE radiation diffusion and energy exchange part is solved with an implicit method. The module is implemented in the PLUTO environment. It uses databases of radiative quantities that can be provided independently by the user: the radiative power loss, the Planck and Rosseland mean opacities. In our case, these quantities were determined from a collisional-radiative steady-state model (CRSS), and tabulated as functions of temperature and density.

*Results. *Our implementation has been validated through different tests, in particular radiative shock tests. The agreement with the semi-analytical solutions (when available) is good, with a maximum error of $7\%$. Moreover, we have proved that non-LTE approach is of paramount importance to properly model accretion shock structures.

*Conclusions. *Our radiation FLD module represents a step toward the general non-LTE RHD modeling. The module is available, under request, for the community.

1 Introduction

The description of fluid motion in many astrophysical systems requires to consider the effects of radiation through its momentum and energy exchanges with matter, which are described by the radiation hydrodynamics (RHD) equations. Examples of relevant systems include, to name just a few, star formation structures (e.g., Krumholz et al. 2009; Commerçon et al. 2011a; Vaytet et al. 2012, 2013a; Davis et al. 2014; Skinner & Ostriker 2015; Raskutti et al. 2016), protoplanetary discs (e.g., Flock et al. 2013, 2017), accretion flows around young stellar objects such as T Tauri stars (Costa et al. 2017; Colombo et al. 2019; de Sá et al. 2019), accretion around black holes (e.g., Hirose et al. 2009; Jiang et al. 2014). In all these cases, radiation is coupled with matter from the dynamical and energetic points of view.

A direct coupling of the radiation effects to the hydrodynamic equations (RHD) or even the magnetohydrodynamic equations (RMHD) requires to solve the radiative transfer equation (RTE) and, in non-local thermodynamic equilibrium (non-LTE or NLTE) regimes, the kinetic equilibrium equations, at each time step, in order to infer the radiation four-force density vector that describes the momentum and energy coupling between radiation and matter. However, this is a challenge that is, still to date, far beyond the current capabilities of computers. The main reason is that the RTE is an integro-differential equation, whose unknown, the specific intensity, depends on seven variables in a three-dimensional description (e.g., Pomraning 1973; Mihalas & Mihalas 1984; Castor 2004; Hubeny & Mihalas 2014). Fully solving the 3D radiative transfer equation itself requires a dedicated approach (e.g., Ibgui et al. 2013) that can be used for the generation of synthetic spectra, without coupling to hydrodynamic evolution.

A work around is to solve directly the equations that involve the radiation moments. These are obtained from angular moments of the RTE. Such an approach entails the creation of a hierarchy of moments, and necessitates the use of closure relations. In fact, the equations describing the conservation of mass, matter momentum, and matter energy, are solved together with the equations describing the conservation of radiation momentum and radiation energy.

The most accurate of the radiation moment techniques is the variable Eddington tensor (VET) method (Stone et al., 1992; Gehmeyr & Mihalas, 1994; Gnedin & Abel, 2001; Hayes & Norman, 2003; Hubeny & Burrows, 2007; Jiang et al., 2012). The method has been applied to the total, or frequency-integrated, radiation moment equations. The Eddington tensor is defined as the ratio of the radiation pressure to the radiation energy. The method consists in solving the moment equations, while assuming that the Eddington tensor is known. Then, once the structure of the medium is determined, the Eddington tensor is updated by solving the RTE with the short-characteristics method (Davis et al., 2012). Despite its good precision, this method is very costly from the computational point of view.

Another technique, less accurate than the VET method, but much cheaper from the computational point of view, is the so-called M1 approximation (Levermore, 1984; Dubroca & Feugeas, 1999). The difference with the VET method, for the total radiation moment equations, is that the Eddington tensor is provided by an analytical closure relation; we do not solve the RTE. The M1 method has been implemented in several multidimensional RHD or RMHD codes, not only in the frequency-integrated approach (HERACLES: González et al. 2007, ATON: Aubert & Teyssier 2008, RAMSES: Rosdahl et al. 2013, ATHENA: Skinner & Ostriker 2013, PLUTO: Melon Fuksman & Mignone 2019), but also in the multigroup approach (HERACLES: Vaytet et al. 2013b, FORNAX: Skinner et al. 2019).

Finally, we mention the third, radiation moment based, technique: the Flux-limited Diffusion approximation (FLD) (Alme & Wilson, 1973; Levermore & Pomraning, 1981). Though less accurate than M1, especially in the optically thin regions, it is the most widely used method in the R(M)HD codes, because it is the simplest one, the most robust, the most efficient in terms of computational cost, and provides very good results in flows where optically thick regions are of paramount importance. The difference with M1 is that the radiation flux equation is not solved; instead, the comoving-frame (also known as the fluid frame; i.e., the frame moving with the macroscopic velocity of the fluid) radiation flux is determined, through a Fick’s law, as a quantity that is antiparallel to the gradient of the radiation energy density, and that recovers, thanks to an ad-hoc function (the flux-limiter), the correct relations between flux and energy in the two asymptotic regimes: free-streaming and (static or dynamic) diffusion. The FLD method has been implemented in many multidimensional RHD or RMHD codes in astrophysical context, in the frequency-integrated approach (ZEUS: Turner & Stone 2001, ORION: Krumholz et al. 2007, V2D: Swesty & Myra 2009, NIRVANA: Kley et al. 2009, RAMSES: Commerçon et al. 2011b, CASTRO: Zhang et al. 2011, PLUTO: Kolb et al. 2013; Flock et al. 2013), and in the multigroup approach (CRASH: van der Holst et al. 2011, CASTRO: Zhang et al. 2013, FLASH: Klassen et al. 2014, RAMSES: González et al. 2015).

All the above R(M)HD codes that use the M1 or the FLD approximations are restricted to the Local Thermodynamic Equilibrium (LTE) approximation. While the LTE regime is a good approximation in the optically thick parts of a flow where densities are high enough, it can no longer be advocated in optically thin parts, such as in the post-shock region of an accretion column (e.g., Ardila 2007). It is therefore important to get free from the restrained LTE assumption, and to develop more general RHD algorithms, that would be capable of handling a non-LTE regime, while the LTE regime could just be a particular case to which the non-LTE regime would tend in parts of a flow that would have appropriate density and temperature conditions. In this work, we present an extended version of the radiation module for PLUTO, originally developed by Kolb et al. (2013) in the LTE approximation, using the FLD method with the frequency-integrated comoving-frame radiation quantities. We have expanded the capabilities of this module, so that it can now allow for non-LTE regimes. We have adopted an implicit scheme to couple the gas and radiation energy exchange, after linearization of a non-LTE radiative emission function, as implemented for a line emission function in ORION (Cunningham et al., 2011; Krumholz et al., 2011, 2012).

The paper is structured as follows: In Sect. 2, we introduce the full system of equations to describe radiation hydrodynamic. We justify all the assumptions used in this work, and we present the approximated equations solved by the code. In Sect. 3, we explain the theoretical model used to generate the opacity and the radiative power loss databases used by our module. In Sect. 4, we explain the numerical implementation used in the code. In Sect. 5, we validate the modifications made in the code through some test cases. In Sect. 6, we compare the two radiative shock structures obtained by letting evolve a given flow under, respectively, forced LTE regime, and non-LTE regime. Finally, in Sect. 7, we draw our conclusions.

2 The equations

The originality of our radiation module coupled to the three-dimensional (3D) MHD code PLUTO lies in the fact that it considers the non-LTE radiation regime: we have removed the LTE assumption that is commonly adopted in radiation hydrodynamics implementations. A bunch of RHD modules and codes exist and are explained in the literature, with a wide variety of approximations and numerical models. It has, therefore, appeared useful to us to start from the basic general physical RHD equations, formulated in the comoving frame, then to clarify and list the approximations that can be done in order to end up with the simplified equations that are actually solved by radiation modules. Our purpose is to start from the most general RHD equations and to introduce, successively, and step by step, all the simplifying assumptions that we have considered, until we obtain the RHD equations that are actually solved by our module and coupled to the MHD equations solved by PLUTO: these are synthesized in the next section (§ 2.3). We also discuss each of these assumptions and their consequence on the RHD modeling of a flow. Kolb et al. (2013) have chosen to express the radiation quantities in the comoving frame. We have followed their choice of reference frame.

2.1 The general comoving-frame RHD equations

We start with the general comoving-frame frequency-integrated RHD equations, for a nonrelativistic flow, a nonviscous and nonconducting fluid, but subject to external gravity, and without nuclear reactions. In the whole paper, and for all applications of our module, we will always restrain to such flows. The equations are, to $O\left(\varv/c\right)$ (e.g., Mihalas & Mihalas 1984):

[TABLE]

Equations (1a), (1b), and (1c) respectively describe the conservation of mass (continuity equation), matter momentum, and matter total energy. Equations (1d) and (1e) respectively describe the total, or frequency integrated, radiation energy conservation and radiation momentum conservation.

We use the following notations: $\rho$ is the mass density, $\boldsymbol{\varv}$ is the fluid (macroscopic) velocity, $a\equiv\partial\boldsymbol{\varv}/\partial t$ is the local fluid acceleration, $p$ is the matter pressure, $\boldsymbol{g}$ is the external gravity, $\mathscr{E}$ is the matter internal energy density. In our approach, all radiation quantities are considered in the comoving frame. Following Mihalas & Mihalas (1984), we denote them with a subscript “0”. $G_{0}^{0}$ and $\boldsymbol{G}_{0}$ are, respectively, the time component and the space components of the radiation four-force density vector $G^{\alpha}_{0}\equiv\left(G_{0}^{0},\boldsymbol{G}_{0}\right)$, and are defined by (Mihalas & Mihalas, 1984):

[TABLE]

where $I_{0}\left(\boldsymbol{r},t;\boldsymbol{n}_{0},\nu_{0}\right)$ is the comoving-frame specific intensity defined at position $\boldsymbol{r}$, time $t$, for the comoving-frame direction unit vector $\boldsymbol{n}_{0}$ and frequency $\nu_{0}$, and where $\chi_{0}\left(\boldsymbol{r},t;\boldsymbol{n}_{0},\nu_{0}\right)$ and $\eta_{0}\left(\boldsymbol{r},t;\boldsymbol{n}_{0},\nu_{0}\right)$ are, respectively, the extinction coefficient (${\rm cm^{-1}}$) and emission coefficient or emissivity (${\rm erg}\;{\rm cm^{-3}}\;{\rm s^{-1}}\;{\rm sr^{-1}}\;{\rm Hz^{-1}})$:

[TABLE]

where $\kappa^{\rm th}_{0}\left(\boldsymbol{r},t;\boldsymbol{n}_{0},\nu_{0}\right)$, $\kappa^{\rm sca}_{0}\left(\boldsymbol{r},t;\boldsymbol{n}_{0},\nu_{0}\right)$, $\eta^{\rm th}_{0}\left(\boldsymbol{r},t;\boldsymbol{n}_{0},\nu_{0}\right)$, $\eta^{\rm sca}_{0}\left(\boldsymbol{r},t;\boldsymbol{n}_{0},\nu_{0}\right)$ are, respectively, the comoving-frame thermal absorption coefficient, scattering coefficient, thermal emission coefficient, and scattering emission coefficient.

Also, $E_{0}$, $\boldsymbol{F}_{0}$, and $\mathbb{P}_{0}$ are the total, comoving-frame, energy density, radiation flux, and radiation pressure:

[TABLE]

The radiation moment approach couples the three matter-related conservation equations (1a), (1b), (1c), an equation of state (EOS) for matter, $p=p(\rho,T)$, an equation for the matter internal energy, or caloric equation of state, $\mathscr{E}=\mathscr{E}(\rho,T)$, and the two radiation moment equations (1d) and (1e). This system needs to be closed by relations that link the three radiation moments (cf. § 2.2.4).

2.2 The non-LTE approximate RHD equations

2.2.1 The approximate radiation four-force density

Our objective is to determine the comoving-frame radiation four-force density vector $\left(G_{0}^{0},\boldsymbol{G}_{0}\right)$ (cf. Eq. (2a) and (2b)), without having to calculate the specific intensity.

We focus on $G_{0}^{0}$. We assume coherent and isotropic scattering in the comoving frame. Then, one can show (Mihalas & Auer, 2001) that

[TABLE]

Then, the scattering terms in Eq. (2a) cancel out. We also assume isotropic thermal absorption $\kappa^{\rm th}_{0}(\boldsymbol{r},t;\nu_{0})$. Finally, we obtain:

[TABLE]

where $L_{0}$ is the total radiation emission per unit volume and time (or radiative power loss):

[TABLE]

and where $\kappa_{0E}$ is the energy-weighted (or absorption) mean opacity, which is defined in the comoving frame as follows:

[TABLE]

where we use the simplified notation $\kappa^{\rm th}_{\nu_{0}}=\kappa^{\rm th}_{0}(\boldsymbol{r},t;\nu_{0})$.

We focus on $\boldsymbol{G}_{0}$. We assume the additional assumption of isotropic thermal emission $\eta^{\rm th}_{0}(\boldsymbol{r},t;\nu_{0})$. Then, the emission terms vanish in Eq. (2b). We end up with:

[TABLE]

where $\bbchi_{0F}$ is the flux-weighted (or radiation momentum) mean opacity, defined in the comoving frame, as follows (Mihalas & Auer, 2001):

[TABLE]

where we use the simplified notation $\chi_{\nu_{0}}=\chi_{0}\left(\boldsymbol{r},t;\nu_{0}\right)$.

2.2.2 The approximate mean opacities

The mean opacities $\kappa_{0E}$ (8) and $\bbchi_{0F}$ (10a) are not known in advance, since they depend on the unknown radiation energy $E_{\nu_{0}}$ and radiation flux $\boldsymbol{F}_{\nu_{0}}$. Our objective is to provide approximate expressions for these quantities.

If we consider an optically thick medium in the equilibrium diffusion regime, we know that (Mihalas & Mihalas, 1984):

[TABLE]

where $B_{\nu_{0}}=B\left(\nu_{0},T_{0}\right)$ is the Planck function at material temperature $T_{0}$. As a consequence, using Eq. (8):

[TABLE]

where $\kappa_{0P}$ is the Planck mean opacity, defined as:

[TABLE]

From Eq. (11b), we find, after frequency integration, that:

[TABLE]

where $B_{0}=ca_{R}T_{0}^{4}/4\pi$ is the frequency-integrated Planck function ($a_{R}$: radiation density constant), and where $\chi_{0R}$ is the Rosseland mean opacity, defined as:

[TABLE]

In the equilibrium diffusion regime, relations (6) and (9) can be replaced with:

[TABLE]

Let us consider now the optically thin regime. In that case, $F_{\nu_{0}}\to cE_{\nu_{0}}$ (Mihalas & Mihalas, 1984; Krumholz et al., 2007). Then, we have $\chi^{(i),\,i=1,3}_{0F}\to(\kappa_{0E}\text{ + a scattering term})$. If, for example, we consider Thomson scattering by free electrons of density $n_{e}$, with (frequency independent) cross section $\sigma_{e}$, we have $\chi^{(i),\,i=1,3}_{0F}\to\kappa_{0E}+\sigma_{e}\,n_{e}$. Now, one can show (Mihalas & Mihalas, 1984) that, in this optically thin regime, and assuming LTE thermal emission (Kirchhoff’s law), the Planck mean $\kappa_{0P}$ is the appropriate mean to use in order to best approximate $G_{0}^{0}$. Then, if we neglect scattering, we could, under the LTE thermal emission assumption, adopt, in optically thin regime, $\chi^{(i),\,i=1,3}_{0F}=\kappa_{0E}=\kappa_{0P}$. However, this assumption is not physically consistent, because, in an optically thin medium, we are in a non-LTE regime. And, in our approach, we do not assume LTE a priori. As a result, we choose not to set $\chi^{(i)}_{0F}=\kappa_{0P}$ in the optically thin regime, because there is no reason for such a specification to improve the accuracy.

In conclusion, first, we approximate $\kappa_{0E}$ with the Planck mean $\kappa_{0P}$; second, we consider the flux spectrum to be the same along each direction, and then replace the tensor $\bbchi_{0F}$ with a scalar $\chi_{0F}$ that we approximate with the Rosseland mean $\chi_{0R}$. We use the approximate expressions (16a) and (16b) to calculate the components $G_{0}^{0}$ and $\boldsymbol{G}_{0}$ of the radiation four-force density vector. This way, we follow the idea of Krumholz et al. (2007) of giving accuracy priority to optically thick parts of the flow. The Planck mean opacity $\kappa_{0P}$, the Rosseland mean opacity $\chi_{0R}$, and the total radiation emission, or radiative losses, $L_{0}$, have been tabulated as functions of density and temperature, for a solar composition of the plasma, assuming a non-LTE regime (Rodríguez et al., 2018). A summary of the adopted physical assumptions is provided in Sect. 3

Note that our radiation module never assumes a priori LTE, but always considers, from the outset, the non-LTE equations. Now, if, at a given time and position, the properties of the simulated medium are such as the LTE conditions prevail, then the calculations will provide the same results as one would obtain from LTE equations: $c\,G_{0}^{0}\xrightarrow[\text{LTE}]{}\,c\,\kappa_{0P}\,E_{0}-\kappa_{0P}\,4\,\pi\,B_{0}$.

Within the above approximations, we can rewrite the system (1a) - (1e) as follows:

[TABLE]

2.2.3 Scaling of terms in the radiation moment equations

The orders of magnitude of the terms involved in the LTE total radiation moment equations have been evaluated in the comoving-frame equations (Mihalas & Mihalas, 1984; Stone et al., 1992), and in the mixed-frame equations (Mihalas & Mihalas, 1984; Krumholz et al., 2007; Skinner & Ostriker, 2013). We apply here the same approach to the non-LTE equations.

2.2.3.1 The three asymptotic physical regimes

We summarize the characteristics of the three asymptotic physical regimes, as classified by Mihalas & Mihalas (1984) and used by Krumholz et al. (2007) and Skinner & Ostriker (2013). We denote by $\ell$ the characteristic structural length at a given position and along a given direction in the flow, and $\lambda_{p}(\nu_{0})\equiv 1/\chi_{0}(\nu_{0})$ the photon mean free path (distance traveled by a photon before being thermally absorbed or scattered). The optical depth along the distance $\ell$ in the flow, at frequency $\nu_{0}$, is $\tau(\nu_{0})\equiv\chi_{0}(\nu_{0})\,\ell=\ell/\lambda_{p}(\nu_{0})$. We distinguish the free-streaming regime in an optically thin medium ($\tau(\nu_{0})\ll 1$), where a photon can move freely in the medium at the speed of light, and the diffusion regime in an optically thick medium ($\tau(\nu_{0})\gg 1$), where a photon travels one free path between interactions (thermal absorption or scattering) with matter. In the diffusion case, we differentiate two subcases. The first one is the static diffusion regime for media with no or small enough velocity. Photons, which are trapped in the material, diffuse through a random walk process. The second one is the dynamic diffusion regime for media with large enough velocity. Photons, which are also trapped in the material, can be advected by the material motion faster than they can diffuse. The distinction between these two subcases in an optically thick medium can be quantified by appropriate parameters (Mihalas & Klein, 1982; Mihalas & Mihalas, 1984; Krumholz et al., 2007; Skinner & Ostriker, 2013). We introduce the fluid-flow timescale $t_{f}\equiv\ell/\varv_{f}$, the typical time for a fluid particle to cross distance $\ell$ at characteristic velocity $\varv_{f}$, the diffusion timescale $t_{d}\equiv\ell^{2}/c\lambda_{p}$, the typical time for a photon to cross distance $\ell$ through a random walk composed of steps of length $\lambda_{p}$ (mean free path) covered at the speed of light (Mihalas & Mihalas, 1984). The corresponding characteristic photon diffusion velocity is $\varv_{d}\equiv\ell/t_{d}$. We can then express quantitatively the physical distinction, discussed above, between the static diffusion regime, $\varv_{d}\gg\varv_{f}$ (equivalent to $t_{d}\ll t_{f}$), and the dynamic diffusion regime, $\varv_{f}\gg\varv_{d}$ (equivalent to $t_{f}\ll t_{d}$). We introduce a characteristic $\beta\equiv\varv_{f}/c\ll 1$ (nonrelativistic flow). This way, we can synthesize the quantitative classification of the three asymptotic regimes as follows (Krumholz et al., 2007):

[TABLE]

where $\tau(\nu_{0})$ is the monochromatic optical depth as defined above. We extend this classification to the frequency-integrated optical depth.

2.2.3.2 Relative sizes of the terms

We analyse the relative sizes of the terms of the comoving-frame radiation moment equations (17d) and (17e), by evaluating the order of magnitude of the evolution of each term on the characteristic structural length $\ell$. We have already introduced, above, the characteristic times for a photon to cross $\ell$: the diffusion timescale $t_{d}$ in a static diffusion regime, and the fluid-flow timescale $t_{f}$ in a dynamic diffusion regime. In a free-streaming regime, the characteristic time for a photon to cross $\ell$ is the radiation-flow timescale $t_{r}\equiv\ell/c$.

One can link the three total radiation moments, $E_{0}$, $\boldsymbol{F}_{0}$, $\mathbb{P}_{0}$ in the three asymptotic regimes described above. One can show that (Mihalas & Mihalas, 1984) in the free streaming regime:

[TABLE]

One can also show that (Mihalas & Mihalas, 1984) in the, either static of dynamic, first order equilibrium diffusion regime:

[TABLE]

In addition, within the second-order equilibrium diffusion approximation (Mihalas & Mihalas, 1984; Castor, 2004), the comoving-frame radiation energy $E_{0}$ has the following orders of magnitude, in, respectively, the static diffusion regime, and the dynamic diffusion regime:

[TABLE]

Using the results above, we can now estimate the orders of magnitude of each term of Eq. (17d) and (17e). We globally follow the usual approaches described in Mihalas & Klein (1982); Mihalas & Mihalas (1984); Stone et al. (1992); Krumholz et al. (2007); Skinner & Ostriker (2013). As did Krumholz et al. (2007) and Skinner & Ostriker (2013), we adopt the characteristic length $\ell$ in the three asymptotic regimes, for the fluid and radiation quantities. We adopt the characteristic timescale $t_{f}$ for the fluid quantities. We adopt the following characteristic timescales for the radiation quantities: $t_{r}$ in the free-streaming regime, $t_{d}$ in the static diffusion regime, and $t_{f}$ in the dynamic diffusion regime. The spatial derivatives are characterized by $1/\ell$. The time derivatives are characterized by $1/\text{(characteristic timescale)}$.

Table 1 and Table 2 display the relative scaling of each term in, respectively, Eq. (17d) and Eq. (17e). In each asymptotic regime, the normalization parameter is indicated in the last column. We also have sorted, in each regime, the terms according to their relative orders of magnitude, from the largest ones (rank 1), to the smallest ones (rank 3 or rank 4). The tables show that the terms with acceleration ( (f) in Table 1, (g) and (h) in Table 2) are from one to three orders of magnitude smaller than the leading terms. Then, one can safely always neglect them, as already stated by Mihalas & Mihalas (1984). Moreover, the terms that account for the time derivative of the radiation momentum ( (c) in Table 1), and of the radiation pressure ( (f) in Table 2), are from one to two orders of magnitude smaller than the leading terms. Then, one can also safely always neglect them. In Table 1, the relative order of magnitude of the radiation energy source - sink term, $-\,c\,G_{0}^{0}$, cannot be estimated a priori in non-LTE (it is $\tau$ in LTE). That is why, its rank is considered to be one by default. Finally, in Table 2, the terms (c), (d), and (e) are from one to two orders of magnitude smaller than the leading terms. So, we can neglect them, to get a solution valid to $O(1)$.

2.2.3.3 The radiation moment equations to $O(1)$

The discussion above on the orders of magnitude suggests that, if we keep only the terms that are leading in at least one of the three asymptotic regimes, then the total non-LTE radiation moment equations are:

[TABLE]

Note that, if we want to solve the radiation moment equations to the next order, i.e., by keeping the terms that are $O(\beta)$ relative to the leading terms, then, we should add term $(c)$ from Table 1 in the radiation energy equation, and terms $(c)$, $(d)$, $(e)$, and $(f)$ from Table 2 in the radiation momentum equation.

2.2.4 The flux-limited diffusion approximation

To close the system, we apply the flux-limited (FLD) approximation (Alme & Wilson, 1973). The FLD approach, which has been widely used in many RHD codes (cf. § 1 for a review) consists in replacing the radiation momentum equation with a Fick’s law of diffusion that links the comoving total radiation flux to the comoving total radiation energy through a radiation diffusion coefficient $K$, as written below:

[TABLE]

where $\lambda$ (not to be confused with the photon mean free path $\lambda_{p}(\nu_{0})$ defined in § 2.2.3.1) is the so-called flux-limiter, a dimensionless quantity that should be defined so that the relation between the radiation flux and the radiation energy is correct in the optically thin and optically thick asymptotic regimes. In other words, we must choose $\lambda$ so that:

[TABLE]

where the dimensionless quantity $R$ is defined as follows:

[TABLE]

This way, we recover the asymptotic relations between $E_{0}$ and $\boldsymbol{F}_{0}$, in an optically thin medium, Eq. (19a), and in an optically thick medium, Eq. (20a). The flux-limiter is then defined as a function of $R$. Different functions are proposed in the literature. The following three ones, respectively suggested by Levermore & Pomraning (1981) (Eq. 26a), Minerbo (1978) (Eq. 26b), and Kley (1989) (Eq. 26c), were implemented in PLUTO by Kolb et al. (2013):

[TABLE]

Finally, one needs a closure relation to relate the radiation pressure to the radiation energy. The most commonly used one (e.g., by Turner & Stone 2001, Krumholz et al. 2007, Zhang et al. 2011, Zhang et al. 2013) is that provided by Levermore (1984):

[TABLE]

where $f$ is the Eddington factor, a dimensionless quantity, related to $\lambda$ and $R$ as follows:

[TABLE]

where the limits are inferred from (24) and (25). Then, we verify that (27) recovers the following asymptotic relations: (19b) in an optically thin medium, and (20b) in an optically thick medium.

2.3 The non-LTE RMHD equations solved by PLUTO

We have implemented in a module, which is coupled to PLUTO, the radiation terms that account for a non-LTE radiation hydrodynamics (RHD) description of a flow, by expanding the LTE RHD equations (Kolb et al., 2013) to the more general non-LTE regime. Since PLUTO is a MHD code, we have actually enhanced its capabilities, so that it is now a 3D non-LTE radiation magnetohydrodynamics (RMHD) code.

We use the RHD equations within the approximations detailed throughout § 2.2. The full set of 3D nonrelativistic RMHD equations solved by PLUTO is based on the full MHD equations (Mignone et al., 2007, 2012), and on Equations (17a), (17b), (17c), (22a), and (23). The system is written as

[TABLE]

where $B$ is the magnetic field, $\lambda$ is the flux-limiter (defined in § 2.2.4), $k_{\rm B}$ is the Planck constant, $\gamma$ is the ratio of specific heats, $\mu$ is the mean molecular weight, $m_{\rm H}$ is the standard mass of one atom of hydrogen, $\boldsymbol{F}_{c}$ is the conductive flux, $\bbtau_{\text{visc}}$ is the viscous tensor, $\eta$ is the magnetic resistivity tensor, and $\boldsymbol{J}\equiv\nabla\times\boldsymbol{B}$ is the current density (for more information on these quantities, see Mignone et al. 2007, 2012; Orlando et al. 2008).

The opacities $\kappa_{0P}$, $\chi_{0R}$, and the radiative losses $L_{0}$ can be provided separately. In our case, we use databases pre-calculated as functions of density and temperature in a non-LTE regime (cf. § 2.2.2).

Note that an additional approximation is made in the radiation energy equation (29d) solved by PLUTO. The terms $\nabla\cdot\left(E_{0}\,\boldsymbol{\varv}\right)$ and $\mathbb{P}_{0}:\nabla$ are missing, compared to (22a). They were already missing in the LTE RHD version by Kolb et al. (2013) (cf. their equation 4). Table 1 reveals that these two terms are leading terms only in the dynamic diffusion regime, but are of second order in the static diffusion regime and in the free streaming regime. So, at this stage, our module cannot be used in the dynamic diffusion regime. We plan to implement these two terms in a future upgraded version of our RHD module.

3 Non-LTE opacities and radiative power loss: theoretical model

Our module needs, at a given set of density and temperature $(\rho,T)$, and for a given composition of the plasma, the following input data: the radiation emission $L_{0}$, the Planck mean opacity $\kappa_{0P}$, and the Rosseland mean opacity $\chi_{0R}$.

Databases have been generated by Rodríguez et al. (2018) in a non-LTE regime, but also in LTE regime. Our module uses three tables, ($L_{0},k_{0P},k_{0R}$) as functions of ($\rho$, $T$), for a plasma with solar-like abundances, where $k_{0P}$ and $k_{0R}$ are (${\rm cm^{2}}\;{\rm g^{-1}}$):

[TABLE]

In this section we briefly summarize the features of the theoretical model, whose details are explained in Rodríguez et al. (2018). Note that our module can read any other set of data ($L_{0},k_{0P},k_{0R}$) that would be provided by the user.

Plasma radiative properties depend on the plasma level populations and atomic properties. In the present work, the atomic quantities of the different chemical elements of the multi-component plasma, such as the relativistic energy levels, wave functions, oscillator strengths and photoionization cross sections, were calculated using the FAC code (Gu, 2008), in which a fully relativistic approach based on the Dirac equation is solved. The atomic calculations were carried out in the relativistic detailed configuration account (Bauche et al., 1987). The atomic configurations selected for each ion in the plasma mixture were those with energies within two times the ionization energy of the ground configuration of the ion.

The atomic level populations were obtained assuming the plasma in steady-state. This approach is valid when the characteristic time of the must relevant atomic process in the plasma is considerably shorter than the time associated to changes in the plasma density and temperature, i.e., the characteristic time of the plasma evolution. When this criterion is fulfilled, the atomic processes are fast enough to distribute the atomic level populations in the plasma before the density and temperature of the plasma change. This approach is commonly used to obtain the atomic level populations in the plasma needed in the calculation of radiative properties databases for radiation-hydrodynamics simulations. In the steady-state approximation, the population density of the atomic level $i$ of the ion with charge state $\zeta$, denoted as $N_{\zeta i}$ is obtained by solving the set of rate equations implemented in a collisional-radiative steady-state (CRSS) model, given by

[TABLE]

Two complementary equations have to be satisfied together with Eq. (31). First, the requirement that the sum of all the partial densities equals the total ion density and, second, the charge neutrality condition in the plasma. In Eq. (31), $R^{+}_{\zeta^{\prime}j\rightarrow\zeta i}$ and $R^{-}_{\zeta^{\prime}j\rightarrow\zeta i}$ take into account all the atomic processes which contribute to populate and depopulate the atomic configuration $\zeta_{i}$, respectively. The atomic processes included in the CRSS model were collisional ionization (Lotz, 1968) and three-body recombination, spontaneous decay (Gu, 2008), collisional excitation (van Regemorter, 1962) and deexcitation, radiative recombination (Kramers, 1923), autoionization and electron capture (Griem, 1997). The rates of the inverse processes were obtained through the detailed balance principle. In the simulations carried out in this work, the plasma was assumed as optically thin. The effect of the plasma environment on the population of the atomic levels was modeled through the depression of the ionization potential or continuum lowering, which can reduce the number of bound states available. The formulation developed by Stewart & Pyatt (1966) was applied. In the present CRSS model, the ions were considered to be at rest. On the other hand, in the calculation of the rates of the atomic processes, a Maxwell-Boltzmann distribution for the free electrons was assumed. For electron densities between $10^{11}$ and $10^{14}$ cm*-3* and electron temperatures lower than $200$ eV, the electron mean free paths range between $3.33\times 10^{5}$ and $25.8$ cm (Rodríguez et al., 2018). This property provides an estimation of the average volume needed for the free electrons to thermalize. For the range of plasma conditions analyzed in this work, the Fermi-Dirac distribution is not necessary. The CRSS model described is implemented in MIXKIP code (Espinosa et al., 2017).

As said before, in this work a multi-component plasma was simulated. The chemical elements considered in the mixture were H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca and Fe, and the solar-like abundances provided by Asplund et al. (2009) were used. For a given electron density and temperature, the CRSS model is solved for each element without considering any atomic processes that connect different chemical elements. However, they were coupled through the common pool of free electrons since the plasma level populations of each chemical element have to be consistent with the common electron density. This fact ensures the plasma neutrality (Klapisch & Busquet, 2013).

Once the atomic data and level populations were obtained, the monochromatic absorption coefficient and emissivity of each chemical element $m$, $\kappa_{m}(\nu)$ and $j_{m}(\nu)$, respectively, with $\nu$ the photon frequency, were calculated using the RAPCAL code (Rodríguez et al., 2008, 2010). Both coefficients include the bound-bound, bound-free and free-free contributions. The radiative transition rates were calculated in the electric dipole approximation using the FAC code. The oscillator strengths included a correction to take into account configuration interaction effects due to the mix between relativistic configurations that belong to the same non-relativistic one. The photoionization cross sections were calculated in the distorted wave approach. Complete redistribution hypothesis was assumed for the line profile which included natural, Doppler, electron-impact (Dimitrijevic & Konjevic, 1987) and UTA broadenings (Bauche et al., 1987). The line-shape function was applied with the Voigt profile that incorporated all these broadenings. For the free-free contributions, the Kramers semi-classical expression for the inverse bremsstrahlung cross section was used (Rose, 1992). In order to determine the opacity, $k(\nu)$, from the absorption coefficient, the scattering of photons was also taken into account, and, in RAPCAL, this was approximated using the Thomson scattering cross section (Rutten, 1995). The monochromatic opacities and emissivities of the mixture were obtained from the weighted contributions of the different chemical elements

[TABLE]

where $\rho$ is the mass density of the mixture. From the monochromatic opacities of the mixture, the Planck $k_{0P}$ and Rosseland $k_{0R}$ mean opacities were calculated, as well as the radiative power loss, or radiation emission, $L_{0}$, from the monochromatic emissivity. Figure 1 portrays maps of $L_{0}$, $k_{0P}$, and $k_{0R}$, vs. free electron density $({\rm cm^{-3}})$ and temperature $({\rm K})$ in logarithmic scale, in a non-LTE regime, and in LTE regime.

Both MIXKIP and RAPCAL codes have been successfully tested with experimental results and numerical simulations for plasmas of single elements included in the mixture analyzed in this work (Espinosa et al., 2017; Rodríguez et al., 2008, 2010), both in LTE and non-LTE, and for the plasma mixture in LTE simulations (Rodríguez et al., 2018). A more detailed explanation of both codes can be found in Rodríguez et al. (2018).

4 Implementation of the radiation terms

We explain, here, how the radiative terms of the physical equations are implemented in the code. We have followed and upgraded to the non-LTE case the techniques described for the LTE regime in Kolb et al. (2013) (see also Commerçon et al. 2011b, van der Holst et al. 2011, Zhang et al. 2011, and Zhang et al. 2013). We use an operator-split method. As detailed in this section, advancing the variables during a time step is made in two substeps, an explicit one (§ 4.2), and an implicit one (§ 4.3). The latter involves the radiation emission $L_{0}$, which is an analytical function of temperature in LTE ($L_{0}=k_{0P}\,\rho\,c\,a_{R}\,T^{4}$: cf. § 2.2.2), but has no analytical expression in non-LTE, and is stored in databases (cf. § 3). Our original contribution, with respect to the usual LTE implementations, consists in handling this non-analytical radiation emission term in the treatment of the implicit scheme, as explained below in § 4.3.

4.1 Reformulation of the equations

For the sake of simplicity, we remove from the full system solved by PLUTO, (29a) - (29h), the terms depending on the magnetic field, and rewrite it in a simplified version, for an inviscid fluid without heat-conduction, but subject to external gravity (even though the following discussion can be applied to the full RMHD system):

[TABLE]

where we have omitted, for economy of notation, the subscript [math] in the radiation quantities.

The above system contains 11 unknowns: six principal variables $\left(\rho,\boldsymbol{\varv},p,E\right)$ that are determined by solving the six equations (33a), (33b), (33c), and (33d), and five variables $\left(T,\mathscr{E},\boldsymbol{F}\right)$ that can be inferred from the six principal ones from the FLD relation (33e), the EOS (33f), and the caloric EOS (33g).

At each time step $\Delta t^{n}=t^{n+1}-t^{n}$, determined by the CFL condition (Mignone et al., 2007), we advance the gas and radiation variables from $t^{n}$ to $t^{n+1}$, by solving the above system in two consecutive substeps, as described below. We denote as $\left(\rho^{n},\boldsymbol{\varv}^{n},p^{n},E^{n}\right)$ and $\left(T^{n},\mathscr{E}^{n},\boldsymbol{F}^{n}\right)$ the values of the variables at time $t^{n}$, and $\left(\rho^{n+1},\boldsymbol{\varv}^{n+1},p^{n+1},E^{n+1}\right)$ and $\left(T^{n+1},\mathscr{E}^{n+1},\boldsymbol{F}^{n+1}\right)$ their values at time $t^{n+1}$, after the two substeps are completed.

4.2 Substep 1: explicit hydrodynamics

In this first substep, PLUTO solves the hyperbolic subsystem made of Equations (33a), (33b), (33c), but without the radiation source (also known as source-sink) terms $k_{P}\,\rho\,c\,E-L$ in the gas energy equation (33c). The system to be solved is:

[TABLE]

where we have introduced (33e), (33g) in (33b), (33c). A Godunov type algorithm is applied, with several possibilities of Riemann solvers (see Mignone et al. 2007 for details).

We use a $\ast$ superscript to denote the quantities obtained at the completion of this substep. Some of these quantities have intermediate results that will be updated in the next substep (§ 4.3). Starting from $\left(\rho^{n},\boldsymbol{\varv}^{n},p^{n}\right)$, PLUTO determines $\left(\rho^{\ast},\boldsymbol{\varv}^{\ast},p^{\ast}\right)$, the terms involving the radiation energy density $E^{n}$ being in the source part of the solved system, and, therefore, not updated. The corresponding intermediate temperature is obtained from the EOS (33f):

[TABLE]

4.3 Substep 2: implicit radiation diffusion and source terms

In this second substep, we determine, at time $t^{n+1}$, the radiation energy density $E^{n+1}$, and the gas temperature $T^{n+1}$, by solving the radiation energy equation (33d), and the gas energy equation (33c) without any velocity term (we couple the gas internal energy density evolution rate with the radiation source terms). The system to be solved in non-LTE is:

[TABLE]

where the gas energy density $\mathscr{E}$ is directly related to the gas temperature through the caloric equation of state (33g).

Since timescales for radiation are much smaller than timescales for hydrodynamics, time steps that would be inferred from the CFL condition would lead, if applied to radiation, to impracticable computations. Consequently, we use an implicit scheme for solving the system (36a) - (36b), with the time step $\Delta t^{n}$ provided by the explicit hydrodynamics substep 1.

Starting from $\left(E^{n},\rho^{\ast},T^{\ast}\right)$, we obtain, at the completion of this substep, $\left(E^{n+1},T^{n+1}\right)$. The time discretization of (36a) and (36b) is made according to the following implicit scheme:

[TABLE]

where, while relating $\mathscr{E}$ to $\rho$ and $T$ with (33g), we assume that the mass density $\rho$ is not modified by this implicit substep (following Commerçon et al. 2011b, Zhang et al. 2011, Zhang et al. 2013, and Kolb et al. 2013), and where

[TABLE]

The radiation diffusion coefficient $K$, the flux-limiter $\lambda$, and the dimensionless quantity $R$ have been defined in § 2.2.4.

The space discretization of (37a) is obtained from a finite volume method (now adapted to non-LTE regimes):

[TABLE]

represents the discretized radiation diffusion term, a linear function of $E^{n+1}$ given by Equation (14) of Kolb et al. (2013), and where the indices $i,j,k$ identify the positions of the cell centers in the 3D computational grid.

The radiation emission at time $t^{n+1}$, $L^{n+1}$ (Eq. 38b), is determined from a $1^{\rm st}$ order Taylor expansion, starting from the state $(\rho^{\ast},T^{\ast})$ obtained from substep 1:

[TABLE]

where we have struck out the term in mass density. This is because we have assumed that the mass density $\rho$ is not modified by this substep 2. Moreover, the dependence of the radiation emission $L$ on $T$ is much higher than that on $\rho$ (see Rodríguez et al. 2018 for more details).

Using (40a) inside (37b) and (39a), we obtain, after some algebra (we removed the subscript $i,j,k$ from each term for simplicity):

[TABLE]

Eq. (41a) represents a linear system whose unknown is the radiation energy density $E^{n+1}_{ijk}$, at time $t^{n+1}$, and at all positions $(i,j,k)$ in the computational domain. The system is solved using the PETSc solver, which is already implemented in the original version of the module; also, the boundary conditions for the radiation energy density can be periodic, symmetric, reflective, or with fixed values.

Note that such a linearization procedure of a radiative loss function, coupled to an implicit scheme, was previously implemented by Cunningham et al. (2011) (cf. their equations 7 and 8, to be respectively compared to Eq. (36b) and (36a) above). In their application, the formation of star clusters, they use a tabulated function that represents line emission (Cunningham et al., 2006) superimposed on the LTE opacity due to dust (Semenov et al., 2003). The detailed system of equations, implemented in the ORION code (Krumholz et al., 2007), is presented in Krumholz et al. (2011, 2012).

Once $E^{n+1}$ is determined, the calculation of the temperature $T^{n+1}$ is straightforward by simply applying formula (41b). Then, using the approximation $\rho^{n+1}=\rho^{\ast}$, we update the Rosseland opacity $k_{R}^{n+1}$, and then the flux-limiter $\lambda(R^{n+1})$, therefore the radiation contribution to the right-hand side of Eq. (34b) and (34c). We also immediately obtain the pressure $p^{n+1}$ from the EOS (33f):

[TABLE]

The velocity $\boldsymbol{\varv}$ is not involved in the equations solved by this substep 2. Then, $\boldsymbol{\varv}^{n+1}=\boldsymbol{\varv}^{\ast}$. Finally, from the above quantities, we can infer the gas energy $\mathscr{E}^{n+1}$ by applying (33g), and the radiation flux $\boldsymbol{F}^{n+1}$ by applying (33e).

5 Tests

To validate our implementation of the non LTE equations in the radiation module in PLUTO, we simulated some test cases, and compared our solutions either to analytical, or semi-analytical results, when they exist, or, when appropriate, to the LTE version of the code. All the test cases present in the literature assume the LTE regime. To have a direct comparison with those tests, we always used our non-LTE discretization scheme (cf. § 4.3), but imposed a radiation emission (or radiative losses) $L$ to be equal to the LTE emission ($L=k_{P}\,\rho\,c\,a_{R}\,T^{4}$). Such tests are described below in § 5.1 and § 5.2.

The tests are one-dimensional problems. Even though PLUTO works in one, two or three dimensions, the radiation module was developed only in 3D. For this reason, to model the test cases we used quasi-one-dimensional domains, which are cuboids with a length much longer than the width or the height. The tests were made only in Cartesian grid. Since, we did not change the geometrical terms in Eq. (41a), we did not need to check the results using different grids.

To compare the results with Kolb et al. (2013), we used the solver based on the PETSc library for all test cases, with the GMRES iteration scheme and a block-Jacobi (bjacobi) preconditioner.

5.1 Radiation matter coupling

Our objective is to test the correctness of the implementation of Equation (36b), which couples the evolution of the matter (or gas) internal energy $\mathscr{E}$ with the radiation source terms $k_{P}\,\rho\,c\,E-L$. This equation is solved with an implicit method, in substep 2 of the operator-split scheme (cf. § 4.3). In particular, this test will enable us to check our linearization procedure of the radiation emission $L$ (Eq. 40a), wich results in the expression of the temperature $T^{n+1}$ versus $T^{n}$ (Eq. 41b). Such a test was first proposed by Turner & Stone (2001).

We want to freeze the evolution of any other quantity but the temperature (and, therefore, the gas internal energy). To do so, first, we build our test model as a static and uniform fluid that is initially out of radiative equilibrium, and we suppress external gravity. Second, we assume the radiation energy $E$ to be large enough compared to the gas internal energy, so that, during the energy exchange between matter and radiation throughout the evolution process, $E$ can be considered as a constant quantity. In this respect, the radiation flux, which is related to the gradient of radiation energy through the FLD relation (23), is negligible, and the continuity, matter momentum, and radiation energy equations (33a, 33b, 33d) become irrelevant. The only remaining relevant equation to be solved is the matter energy equation (33c), but in the form (36b).

Even though we use our non-LTE discretization scheme with the radiation emission (or radiative losses) quantity $L$ in the equations (§ 4.3), we set $L$ to be equal to the LTE emission ($L=k_{P}\,\rho\,c\,a_{R}\,T^{4}$). This way, Eq. (36b) can be recast as the following ordinary differential equation:

[TABLE]

$\mathscr{E}_{\rm f}$ is the gas internal energy in the final radiative equilibrium state.

5.1.1 Setup

For this test, we set the density \rho=10^{-7}\leavevmode\nobreak\g cm*-3*, the radiation energy density E=10^{10}\leavevmode\nobreak\erg cm*-3*, the Planck111The choice of the Rosseland mean opacity is irrelevant since it is associated to the radiation flux that is negligible in our test; however, for computational convenience, we also adopted $\rho\,k_{R}=4\times 10^{-8}\leavevmode\nobreak\ \rm{cm}^{-1}$. mean opacity $\rho\,k_{P}=4\times 10^{-8}\leavevmode\nobreak\ \rm{cm}^{-1}$, the mean molecular weight $\mu=0.6$, and the ratio of specific heats $\gamma=5/3$. The domain consists of a cuboid that has a length much larger than the other two dimensions. In particular, the cuboid has a width and a height of 3 cm, and a length of 100 cm. The grid consists of $3\times 3\times 100$ points. All the boundaries are periodic. For this test the hydrodynamic solver is turned off.

The simulation starts at t=0\leavevmode\nobreak\s with an initial time step of $\delta t=10^{-20}$ s. After each step, the time step increases by $0.1\%$. The test is performed using 3 different initial gas internal energy densities \mathscr{E}_{0}=(6.4\times 10^{3};\,6.4\times 10^{7};\,6.4\times 10^{8})\leavevmode\nobreak\erg cm*-3*.

5.1.2 Results

Fig. 2 shows the comparison between the solution found with our model (red dots) and the semi-analytical reference one (black line) obtained by solving (43a) with a fourth order Runge-Kutta scheme (alternatively, a full analytical solution is provided by Swesty & Myra 2009), for the three different initial energy radiation densities $\mathscr{E}_{0}$. The agreement between the two solutions is excellent for all the three cases.

5.2 Radiative shocks

We test our implementation of the full set of non-LTE radiation hydrodynamics equations in PLUTO (33a-33g) (without external gravity), by assessing the ability of the code to reproduce the structure of a radiative shock. We simulate a simple shock case in a quasi-one-dimensional domain; this test follows from Ensman (1994). It is possible to compare some characteristic quantities derived from the test, with analytical estimates (Mihalas & Mihalas, 1984). We compare our results with the original version of the code.

5.2.1 Setup

For this test, we simulate an initially uniform fluid with density \rho=7.78\times 10^{-10}\leavevmode\nobreak\g cm*-3*, and temperature T=10\leavevmode\nobreak\K, moving with velocity $\varv$ along z-axis. We set the ratio of specific heats $\gamma=7/5$, and the mean molecular weight $\mu=1$, in analogy to Ensman (1994). We impose a constant opacity k_{R}\times\rho=k_{P}\times\rho=3.1\times 10^{-10}\leavevmode\nobreak\cm*-1*, and the initial radiation energy density is set by the equation $E=a_{R}T^{4}$. As in the previous test case (§ 5.1), we impose the radiative losses $L=k_{P}\,\rho\,c\,a_{R}\,T^{4}$ (LTE radiation emission), while using our non-LTE discretization scheme with $L$ in the equations.

The computational domain has a width and a height of 3.418\times 10^{7}\leavevmode\nobreak\cm, and a length (z-axis) of 7\times 10^{10}\leavevmode\nobreak\cm. Along with Kolb et al. (2013), we have chosen a grid composed of $4\times 4\times 2048$ cells, the Minerbo flux-limiter (Eq. 26b in § 2.2.4), a Lax-Friedrichs scheme, 0.4 as CFL (Courant-Friedrichs-Lewy) value, and a relative tolerance $\epsilon=10^{-5}$ for the matrix solver. The lateral boundaries are periodic. In the direction of the fluid flow (z-axis), we use a reflective boundary at the bottom of the domain ($z=0$), to generate the shock, and a zero-gradient condition at the top (z=z_{\rm max}=7\times 10^{10}\leavevmode\nobreak\cm).

Moving from $z_{\rm max}$ to $z=0$, the fluid impacts onto the reflective boundary, which generates a shock that propagates back into the fluid. The fluid velocity is calculated with respect to the computational domain, which is the physical frame of reference. The shock can be subcritical (low velocity) or supercritical (large velocity) (Mihalas & Mihalas, 1984). We simulate both cases: we set \varv=6\times 10^{5}\leavevmode\nobreak\cm s*-1* to obtain a subcritical shock, and \varv=20\times 10^{5}\leavevmode\nobreak\cm s*-1* to obtain a supercritical shock.

5.2.2 Results

In order to compare our results with those obtained by Ensman (1994) or Kolb et al. (2013), we introduce the quantity $s=z-\varv\,t$. Fig. 3 shows the gas temperature (red) and the radiation temperature (blue), $T_{\rm rad}=\left(E/a_{R}\right)^{1/4}$, versus $s$ for the subcritical shock (top panel) and the supercritical shock (bottom panel). In the supercritical shock the temperature after and before the shock are equal, as expected.

In the case of a subcritical shock, some characteristic gas temperatures can be estimated analytically (Mihalas & Mihalas, 1984): $T_{-}$, the temperature immediately ahead of the shock front; $T_{+}$, the temperature immediately behind the shock front (Zel’dovich spike: Zel’dovich & Raizer 1967); $T_{2}$, the final equilibrium post-shock temperature, reached after radiative cooling. We have:

[TABLE]

where $R_{G}=k_{B}/\mu m_{\rm H}$ is the perfect gas constant, $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant, and $u$ is the velocity of the shock relative to the upstream fluid.

Table 3 shows the comparison between our numerical solution, the analytical estimate, and the solution reported by Kolb et al. (2013). Our numerical solution agrees very well with the original version of the code. The relative deviation with respect to the analytical estimate is no larger than $7\%$. Moreover, the position and the shape of the shocks are very well reproduced. We can conclude that our modifications in the radiation module have maintained the accuracy of the code.

6 LTE vs non-LTE radiative shocks

The purpose of this section is to show the crucial importance of considering the appropriate regime (LTE vs non-LTE) for given physical conditions, in order to correctly model the structure and dynamics of a radiating fluid. This is because opacities and radiation emissions can differ by several orders of magnitudes between both regimes, as exemplified by Fig. 1. Such deviations can have a major impact on the momentum and energy exchanges between matter and radiation, and, therefore, on the structure of the flow.

We modify the shock test from Ensman (1994), described in § 5.2, to have physical conditions quite similar to those in accretion shocks in young stars (Sacco et al. 2008; Colombo et al. 2019). As in the two test cases, we use the non-LTE discretization scheme with the radiation emission (or radiative losses) quantity $L$ in the equations (§ 4.3). In one case, referred to as the “non-LTE case”, we use the non-LTE radiative database (calculated by Rodríguez et al. 2018: cf. § 3), and let the system evolve following the flow conditions, for which, at a given time and position, either non-LTE or LTE regime prevails and is self-consistently taken into account by the database. In the other case, referred to as the “LTE case”, we use the LTE database (still calculated by Rodríguez et al. 2018), and, therefore, force the LTE regime, no matter the physical conditions.

6.1 Setup

We simulate a uniform fluid with initial density n=10^{12}\leavevmode\nobreak\cm*-3*, and temperature T=2\times 10^{4}\leavevmode\nobreak\K, moving with velocity $\varv=5\times 10^{7}\leavevmode\nobreak\ \rm{cm}^{-1}$ along z-axis. We set $\gamma=7/5$ and $\mu\approx 1.29$, i.e., assume solar abundances. Note that here, unlike the preceding test cases, we do not impose a constant value for $k_{P}$, $k_{R}$ and $L$, but we use the radiative databases. The initial radiation energy $E$ is chosen in order to start with a fluid in radiative equilibrium.

The computational domain describes a box of length 10^{8}\leavevmode\nobreak\cm in the non-LTE case, and 10^{7}\leavevmode\nobreak\cm in the LTE case, and of equal width and height, 3.418\times 10^{7}\leavevmode\nobreak\cm in both cases. The grid is composed of 3x3x1024 cells.

We adopt the Minerbo flux-limiter, and boundary conditions identical to those in the shock tests (§ 5.2). For the CFL condition, we use 0.01. For the matrix solver we set a relative tolerance $\epsilon=10^{-4}$.

6.2 Results

The fluid impacts onto the reflective boundary, which generates a shock that propagates back into the fluid. The shock heats up the material forming a post-shock region. Fig. 4 shows the profiles of temperature, radiative gains ($G=k_{P}\,\rho\,c\,E$), which represent the energy gained by the fluid after absorbing radiation, and radiative losses ($L$) after 6 s of evolution.

The two cases present several differences while we follow their evolutions. This is related to the different regimes taken into account. In the LTE case (Fig.4, top panel), the shock reaches a maximum temperature of \sim 5\times 10^{4}\leavevmode\nobreak\K (dashed line curve). At this peak, the radiative losses (blue curve: $L\sim 5\times 10^{6}\leavevmode\nobreak\ {\rm erg}\;{\rm cm^{-3}}\;{\rm s^{-1}}$) are higher than the radiative gains (red curve: $G\sim 5\times 10^{5}\leavevmode\nobreak\ {\rm erg}\;{\rm cm^{-3}}\;{\rm s^{-1}}$). As a consequence, the material rapidly cools down until radiative equilibrium ($G=L$), and the post-shock region remains relatively cold at \sim 2\times 10^{4}\leavevmode\nobreak\K, close to the initial temperature.

Even though the radiative losses are extremely high compared to the non-LTE case (see below, the last paragraph of this section), there is no precursor region. We can invoke two reasons for this: first, the region that emits is quite small, so the radiation energy emitted per unit time is not enough to heat up the unshocked plasma; second, according to Fig.1 (mid-right panel), the Planck opacity, $k_{P}$, in LTE regime is smaller by several order of magnitudes compared to that in the non-LTE case. Therefore, matter absorbs far less radiation in LTE than in non-LTE (a reminder: the gain of radiation energy by matter is $G=k_{P}\,\rho\,c\,E$).

In the non-LTE regime (Fig.4, bottom panel), the radiative losses are, in the shocked region, smaller by around three orders of magnitude compared to those in the LTE regime: $L_{\,\rm LTE}\sim 5\times 10^{6}\leavevmode\nobreak\ {\rm erg}\;{\rm cm^{-3}}\;{\rm s^{-1}}$ (see above), $L_{\,\rm NLTE}\sim 10^{3}-10^{4}\leavevmode\nobreak\ {\rm erg}\;{\rm cm^{-3}}\;{\rm s^{-1}}$. In this case, the shock heats the material up T\sim 5\times 10^{6}\leavevmode\nobreak\K, and generates a hot post-shock region. Since the radiative losses are smaller than in the LTE case, the shock-heated material needs more time to cool down, and forms a hot slab. After 6\leavevmode\nobreak\s, the radiative losses trigger the thermal instabilities at the base of the post-shock: the material rapidly loses thermal energy through radiation emission. This drop in temperature produces enough radiation energy to heat up the unshocked plasma, thereby generating a precursor region with a temperature of T\sim 10^{5}\leavevmode\nobreak\K (dashed line curve on the right in Fig.4).

7 Conclusions

Including the effects of radiation in HD and MHD models is a mandatory task to fully describe many astrophysical systems. Several codes fulfill this request, but none of them consider the more general non-LTE regime.

Here, we have presented our extended version of the LTE radiation module developed by Kolb et al. (2013) and implemented in the PLUTO code. The upgraded module is now able to handle non-LTE regimes (including, self-consistently, the particular case of LTE regime, depending on the physical plasma conditions). We use an operator-split method. The system is solved in two substeps, an explicit one for the hyperbolic subsystem, and an implicit one for the subsystem that involves radiation diffusion and radiation source terms. It is this second subsystem that we have upgraded so that it can now handle non-LTE conditions.

Starting from the general frequency-integrated comoving-frame radiation hydrodynamics equations, we have reviewed all the assumptions and approximations that have led to the equations that are actually coded in PLUTO. In particular, we use the flux-limited diffusion approximation. Also, our implementation is valid for plasma in conditions ranging from free streaming regime to static diffusion regime. It cannot describe the dynamic diffusion regime: to do so, we have to include the two advection terms in the radiation energy equation. Moreover, the multigroup implementation in non-LTE can be a further improvement of our module. These possible developments will be the subject of future works.

The module needs, as input data versus density and temperature, the following radiative quantities: the radiation emission, the Planck mean opacity, and the Rosseland mean opacity. Our module currently uses non-LTE databases generated with a CRSS model for a plasma with solar-like abundances. But the user may provide any other set of databases that would be more appropriate to the problem to be investigated. In our CRSS approach, the absorption processes, such as photoabsorption and photoionization, are ignored in the determination of the atomic level populations. Considering such processes is a much more difficult problem that leads to stiff rate equations and stiff coupling with the thermodynamic state.

We have tested the new implementation of the non-LTE equations, but in LTE conditions, and have compared our results with semi-analytical solutions, and with the results given by the previous version of the code. These tests have established the validity of our implementation.

We also have demonstrated, through the comparison of the structure of a radiative shock, in LTE and in non-LTE regimes, the importance of considering the appropriate regime in order to correctly describe the dynamics of a radiating fluid.

We have already successfully applied this new upgraded version of PLUTO to demonstrate the existence of a radiative precursor in the accreting stream onto the surface of a classical T Tauri Star (Colombo et al., 2019).

The radiation module has been implemented in version 4.0 of PLUTO. The code is available, under request, on the observatory of Palermo website222http://cerere.astropa.unipa.it/progetti_ricerca/HPC/resources.htm, to the scientific community.

Acknowledgements.

The authors are grateful to Dr. Mario Flock for useful discussions on the implementation of radiation equations in the MHD PLUTO code. PLUTO is developed at the Turin Astronomical Observatory in collaboration with the Department of Physics of Turin University. We acknowledge the ”Accordo Quadro INAF-CINECA (2017)” the CINECA Award HP10B1GLGV and the HPC facility (SCAN) of the INAF – Osservatorio Astronomico di Palermo, for the availability of high performance computing resources and support. This work was supported by the Programme National de Physique Stellaire (PNPS) of CNRS/INSU co-funded by CEA and CNES. This work has been done within the LABEX Plas@par project, and received financial state aid managed by the Agence Nationale de la Recherche (ANR), as part of the programme ”Investissements d’avenir” under the reference ANR-11-IDEX-0004-02.