Issue 
A&A
Volume 631, November 2019



Article Number  A41  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201935991  
Published online  23 October 2019 
NonLTE radiation hydrodynamics in PLUTO^{⋆}
^{1}
Dipartimento di Fisica e Chimica, Université degli Studi di Palermo, via Archirafi 36, Palermo, Italy
email: salvatore.colombo@inaf.it
^{2}
LERMA, Observatoire de Paris, Sorbonne Université, Université de CergyPontoise, CNRS, Paris, France
^{3}
INAF – Osservatorio Astronomico di Palermo, Palermo, Italy
^{4}
Universidad de Las Palmas de Gran Canaria, Gran Canaria, Spain
^{5}
Université Paris Diderot, Sorbonne Paris Cité, AIM, UMR 7158, CEA, 91191 GifsurYvette, France
Received:
30
May
2019
Accepted:
5
July
2019
Context. Modeling the dynamics of most astrophysical structures requires an adequate description of the interaction of radiation and matter. Several numerical (magneto) hydrodynamics codes were upgraded with a radiation module to fulfill this request. However, those that used either the fluxlimited diffusion (FLD) or the M1 radiation moment approaches are restricted to local thermodynamic equilibrium (LTE). This assumption may not be valid in some astrophysical cases.
Aims. We present an upgraded version of the LTE radiationhydrodynamics (RHD) module implemented in the PLUTO code, which we have extended to handle nonLTE regimes.
Methods. Starting from the general frequencyintegrated comovingframe equations of RHD, we have justified all the assumptions that were made to obtain the nonLTE equations that are implemented in the module under the FLD approximation. An operatorsplit method with two substeps was employed: the hydrodynamics part was solved with an explicit method by the solvers that are currently available in PLUTO, and the nonLTE radiation diffusion and energy exchange part was solved with an implicit method. The module was implemented in the PLUTO environment. It uses databases of radiative quantities that can be provided independently by the user: the radiative power loss, and the Planck and Rosseland mean opacities. In our case, these quantities were determined from a collisionalradiative steadystate model, and they are 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 semianalytical solutions (when available) is good, with a maximum error of 7%. Moreover, we have proved that a nonLTE approach is of paramount importance to properly model accretion shock structures.
Conclusion. Our radiation FLD module represents a step toward a general nonLTE RHD modeling.
Key words: radiation: dynamics / hydrodynamics / opacity
The module is available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/631/A41 and upon request to the first author.
© ESO 2019
1. Introduction
The description of fluid motion in many astrophysical systems requires considering 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 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 disks (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), and 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 or even the magnetohydrodynamics equations (RMHD) requires solving at each step the radiative transfer equation (RTE), and in nonlocal thermodynamic equilibrium (nonLTE) regimes, the kinetic equilibrium equations, in order to infer the radiation fourforce density vector that describes the momentum and energy coupling between radiation and matter. However, this is a challenge that still is far beyond the current capabilities of computers. The main reason is that the RTE is an integrodifferential equation, whose unknown, the specific intensity, depends on seven variables in a 3D 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 to generate synthetic spectra without coupling to hydrodynamic evolution.
A workaround is to directly solve the equations that involve the radiation moments. These are obtained from angular moments of the RTE. This approach entails the creation of a hierarchy of moments, and necessitates the use of closure relations. 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 frequencyintegrated, radiation moment equations. The Eddington tensor is defined as the ratio of the radiation pressure to the radiation energy. The method consists of solving the moment equations and assumes that the Eddington tensor is known. When the structure of the medium is determined, the Eddington tensor is updated by solving the RTE with the shortcharacteristics 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 socalled 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 frequencyintegrated 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 fluxlimited diffusion (FLD) approximation (Alme & Wilson 1973; Levermore & Pomraning 1981). Although 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, most robust, and most efficient method in terms of computational cost, and it 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 comovingframe (also known as the fluid frame; i.e., the frame moving with the macroscopic velocity of the fluid) radiation flux is determined, through the Fick law, as a quantity that is antiparallel to the gradient of the radiation energy density, and that recovers through an ad hoc function (the flux limiter), the correct relations between flux and energy in the two asymptotic regimes, the freestreaming and the (static or dynamic) diffusion regime. The FLD method has been implemented in many multidimensional RHD or RMHD codes in astrophysical contexts, in the frequencyintegrated 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 these 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 postshock region of an accretion column (e.g., Ardila 2007). It is therefore important to become independent of the restrained LTE assumption, and to develop more general RHD algorithms that are capable of handling a nonLTE regime, while the LTE regime could just be a particular case to which the nonLTE regime would tend in parts of a flow that would have appropriate density and temperature conditions. We here present an extended version of the radiation module for PLUTO that was originally developed by Kolb et al. (2013) in the LTE approximation, using the FLD method with the frequencyintegrated comovingframe radiation quantities. We have expanded the capabilities of this module, so that it can now allow for nonLTE regimes. We have adopted an implicit scheme to couple the gas and radiation energy exchange after linearization of a nonLTE 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 the radiation hydrodynamics. We justify all the assumptions we use, and we present the approximated equations that are solved by the code. In Sect. 3 we explain the theoretical model we used to generate the opacity and the radiative powerloss 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 we obtained by letting a given flow evolve in a forced LTE regime and in a nonLTE regime. Finally, we draw our conclusions in Sect. 7.
2. Equations
The originality of our radiation module coupled to the 3D magnetohydrodynamics (MHD) code PLUTO lies in the fact that it considers the nonLTE radiation regime: we have removed the LTE assumption that is commonly adopted in radiation hydrodynamics implementations. Several RHD modules and codes exist and are explained in the literature, with a wide variety of approximations and numerical models. It therefore appeared to be 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 made in order to define the simplified equations that are solved by radiation modules. Our purpose is to start from the most general RHD equations and to successively introduce step by step all the simplifying assumptions that we have considered, until we obtained the RHD equations that are solved by our module and were coupled to the MHD equations solved by PLUTO: these are synthesized in the next section (Sect. 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. General comovingframe RHD equations
We start with the general comovingframe frequencyintegrated RHD equations for a nonrelativistic flow, a nonviscous and nonconducting fluid that is subject to external gravity, but has no nuclear reactions. For all applications of our module throughout, we always restrain ourselves to such flows. The equations are to O(v/c) (e.g., Mihalas & Mihalas 1984)
Equations (1a)–(1c) describe the conservation of mass (continuity equation), matter momentum, and total energy of matter, respectively. Equations (1d) and (1e) describe the total, or frequencyintegrated, radiation energy conservation and radiation momentum conservation, respectively.
We use the following notations: ρ is the mass density, v is the fluid (macroscopic) velocity, a ≡ ∂v/∂t is the local fluid acceleration, p is the matter pressure, g is the external gravity, and ℰ 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”. and G_{0} are the time component and the space components of the radiation fourforce density vector , respectively, and are defined by (Mihalas & Mihalas 1984)
where I_{0}(r,t;n_{0},ν_{0}) is the comovingframe specific intensity defined at position r, time t, for the comovingframe direction unit vector n_{0} and frequency ν_{0}, and where χ_{0}(r,t;n_{0},ν_{0}) and η_{0}(r,t;n_{0},ν_{0}) are the extinction coefficient (cm^{−1}) and emission coefficient or emissivity (erg cm^{−3} s^{−1} sr^{−1} Hz^{−1}), respectively,
where , , , are the comovingframe thermal absorption coefficient, scattering coefficient, thermal emission coefficient, and scattering emission coefficient, respectively.
In addition, E_{0}, F_{0}, and ℙ_{0} are the total comovingframe energy density, the radiation flux, and the radiation pressure:
The radiation moment approach couples the three matterrelated conservation Eqs. (1a)–(1c), an equation of state (EOS) for matter, p = p(ρ, T), an equation for the internal energy of matter, or caloric equation of state, ℰ = ℰ(ρ, T), and the two radiation moment Eqs. (1d) and (1e). This system needs to be closed by relations that link the three radiation moments (cf. Sect. 2.2.4).
2.2. NonLTE approximate RHD equations
In this section we start from the RHD equations presented in Sect. 2 and we go through all the assumptions made to obtain the approximated RHD equations. Then, we analyze the relative sizes of the terms present in the equations under different physical regimes. Finally, we present the final RHD set of equations actually solved by the module where the negligible terms are omitted.
2.2.1. The approximate radiation fourforce density
Our objective is to determine the comovingframe radiation fourforce density vector (cf. Eq. (2a) and (2b)), without having to calculate the specific intensity. We focus on . We assume coherent and isotropic scattering in the comoving frame. We can show (Mihalas & Auer 2001) that
Then, the scattering terms in Eq. (2a) cancel out. We also assume isotropic thermal absorption . Finally, we obtain
where L_{0} is the total radiation emission per unit volume and time (or radiative power loss),
and where κ_{0E} is the energyweighted (or absorption) mean opacity, which is defined in the comoving frame as follows:
where we use the simplified notation .
We focus on G_{0}. We assume the additional assumption of isotropic thermal emission . Then, the emission terms vanish in Eq. (2b). We find
where is the fluxweighted (or radiation momentum) mean opacity, defined in the comoving frame as (Mihalas & Auer 2001)
2.2.2. Approximate mean opacities
The mean opacities κ_{0E} Eq. (8) and Eq. (10a) are not known in advance because they depend on the unknown radiation energy E_{ν0} and radiation flux F_{ν0}. Our objective is to provide approximate expressions for these quantities.
When we consider an optically thick medium in the equilibrium diffusion regime, we know that (Mihalas & Mihalas 1984)
where B_{ν0} = B(ν_{0},T_{0}) is the Planck function at material temperature T_{0}. As a consequence, using Eq. (8),
where κ_{0P} is the Planck mean opacity, defined as
From Eq. (11b), we find after frequency integration that
where is the frequencyintegrated Planck function (a_{R}: radiation density constant), and where χ_{0R} is the Rosseland mean opacity, defined as
In the equilibrium diffusion regime, relations (6) and (9) can be replaced by
We consider now the optically thin regime. In this case, F_{ν0} → cE_{ν0} (Mihalas & Mihalas 1984; Krumholz et al. 2007). Then, we have . When we consider Thomson scattering by free electrons of density n_{e}, for instance, with a (frequencyindependent) cross section σ_{e}, we have . We show (Mihalas & Mihalas 1984) that in this optically thin regime, and assuming LTE thermal emission (Kirchhoff’s law), the Planck mean κ_{0P} is the appropriate mean to use in order to best approximate . When we neglect scattering, we could under the LTE thermal emission assumption adopt in optically thin regime, . However, this assumption is not physically consistent because in an optically thin medium, we are in a nonLTE regime. Moreover, in our approach, we do not assume LTE a priori. As a result, we choose not to set in the optically thin regime because there is no reason for such a specification to improve the accuracy.
In conclusion, we first approximate κ_{0E} with the Planck mean κ_{0P}; second, we consider the flux spectrum to be the same in each direction, and then we replace the tensor with a scalar χ_{0F} that we approximate with the Rosseland mean χ_{0R}. We use the approximate expressions (16a) and (16b) to calculate the components and G_{0} of the radiation fourforce density vector. In 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 κ_{0P}, the Rosseland mean opacity χ_{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 nonLTE regime (Rodríguez et al. 2018). A summary of the adopted physical assumptions is provided in Sect. 3
We note that our radiation module never assumes LTE a priori, but always considers from the outset the nonLTE equations. If at a given time and position, the properties of the simulated medium are such that the LTE conditions prevail, then the calculations will provide the same results as would be obtained from LTE equations: .
Within the above approximations, we can rewrite the system (1a)–(1e) as follows:
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 comovingframe equations (Mihalas & Mihalas 1984; Stone et al. 1992), and in the mixedframe equations (Mihalas & Mihalas 1984; Krumholz et al. 2007; Skinner & Ostriker 2013). We here apply the same approach to the nonLTE equations.
2.2.3.1. 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 ℓ the characteristic structural length at a given position and in a given direction in the flow, and λ_{p}(ν_{0})≡1/χ_{0}(ν_{0}) the photon mean free path (distance traveled by a photon before it is thermally absorbed or scattered). The optical depth at the distance ℓ in the flow at frequency ν_{0} is τ(ν_{0})≡χ_{0}(ν_{0}) ℓ = ℓ/λ_{p}(ν_{0}). We distinguish the freestreaming regime in an optically thin medium (τ(ν_{0})≪1), where a photon can move freely in the medium at the speed of light, and the diffusion regime in an optically thick medium (τ(ν_{0})≫1), where a photon travels one free path between interactions (thermal absorption or scattering) with matter. In the diffusion case, we distinguish two subcases. The first is the static diffusion regime for media with no or low enough velocity. Photons that are trapped in the material diffuse through a random walk process. The second subcase is the dynamic diffusion regime for media with high enough velocity. Photons that 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 fluidflow timescale t_{f} ≡ ℓ/v_{f}, the typical time for a fluid particle to cross a distance ℓ at characteristic velocity v_{f}, the diffusion timescale t_{d} ≡ ℓ^{2}/cλ_{p}, the typical time for a photon to cross a distance ℓ through a random walk composed of steps of length λ_{p} (mean free path) covered at the speed of light (Mihalas & Mihalas 1984). The corresponding characteristic photon diffusion velocity is v_{d} ≡ ℓ/t_{d}. We can then quantitatively express the physical distinction, discussed above, between the static diffusion regime, v_{d} ≫ v_{f} (equivalent to t_{d} ≪ t_{f}), and the dynamic diffusion regime, v_{f} ≫ v_{d} (equivalent to t_{f} ≪ t_{d}). We introduce a characteristic β ≡ v_{f}/c ≪ 1 (nonrelativistic flow). In this way, we can synthesize the quantitative classification of the three asymptotic regimes as follows (Krumholz et al. 2007):
where τ(ν_{O}) is the monochromatic optical depth as defined above. We extend this classification to the frequencyintegrated optical depth.
Relative sizes of the terms. We analyze the relative sizes of the terms of the comovingframe radiation moment Eqs. (17d) and (17e) by evaluating the order of magnitude of the evolution of each term on the characteristic structural length ℓ. We introduced above the characteristic times for a photon to cross ℓ: the diffusion timescale t_{d} in a static diffusion regime, and the fluidflow timescale t_{f} in a dynamic diffusion regime. In a freestreaming regime, the characteristic time for a photon to cross ℓ is the radiationflow timescale t_{r} ≡ ℓ/c.
The three total radiation moments, E_{0}, F_{0}, and ℙ_{0} can be linked in the three asymptotic regimes described above. We show that (Mihalas & Mihalas 1984) in the free streaming regime,
We also show that (Mihalas & Mihalas 1984) in either the static or dynamic firstorder equilibrium diffusion regime,
In addition, within the secondorder equilibrium diffusion approximation (Mihalas & Mihalas 1984; Castor 2004), the comovingframe radiation energy E_{0} has the following orders of magnitude in the static and the dynamic diffusion regime:
using the results above, we can now estimate the orders of magnitude of each term of Eqs. (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) and Skinner & Ostriker (2013). As in Krumholz et al. (2007) and Skinner & Ostriker (2013), we adopt the characteristic length ℓ 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 freestreaming regime, t_{d} in the static diffusion regime, and t_{f} in the dynamic diffusion regime. The spatial derivatives are characterized by 1/ℓ. The time derivatives are characterized by 1/(characteristic timescale).
Tables 1 and 2 display the relative scaling of each term in Eqs. (17d) and (17e), respectively. In each asymptotic regime, the normalization parameter is indicated in the last column. We also sorted in each regime the terms according to their relative orders of magnitude, from the largest (rank 1) to the smallest magnitude (ranks 3 or 4). The tables show that the terms with acceleration ( (f) in Table 1, and (g) and (h) in Table 2) are from one to three orders of magnitude smaller than the leading terms. They can therefore be always safely neglect, as has been 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. They can also always safely be neglected. In Table 1 the relative order of magnitude of the radiation energy source – sink term, , cannot be estimated a priori in nonLTE (it is τ in LTE). Its rank is therefore 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. We can therefore neglect them to obtain a solution valid to O(1).
2.2.3.2. 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 nonLTE radiation moment equations are
We note that if we wish to solve the radiation moment equations to the next order, that is, by keeping the terms that are O(β) relative to the leading terms, we need to 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. Fluxlimited diffusion approximation
To close the system, we apply the FLD approximation (Alme & Wilson 1973). The FLD approach, which has been widely used in many RHD codes (cf. Sect. 1 for a review) consists of replacing the radiation momentum equation with the Fick 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:
where λ (not to be confused with the photon mean free path λ_{p}(ν_{0}) defined in Sect. 2.2.3) is the socalled 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 λ so that
where the dimensionless quantity R is defined as follows:
In this way, we recover the asymptotic relations between E_{0} and 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 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):
Finally, a closure relation is required to relate the radiation pressure to the radiation energy. The most commonly used relation (e.g., by Turner & Stone 2001; Krumholz et al. 2007; Zhang et al. 2011, 2013) is provided by Levermore (1984),
where f is the Eddington factor, a dimensionless quantity, related to λ and R as follows:
where the limits are inferred from Eqs. (24) and (25). Then, we verify that Eq. (27) recovers the following asymptotic relations: Eq. (19b) in an optically thin medium, and Eq. (20b) in an optically thick medium.
2.3. NonLTE RMHD equations solved by PLUTO
We implemented in a module that is coupled to PLUTO the radiation terms that account for a nonLTE RHD description of a flow by expanding the LTE RHD equations (Kolb et al. 2013) to the more general nonLTE regime. Because PLUTO is a MHD code, we have thus enhanced its capabilities, so that it is now a 3D nonLTE RMHD code.
We use the RHD equations within the approximations detailed throughout Sect. 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 Eqs. (17a)–(17c), (22a), and (23). The system is written as
where B is the magnetic field, λ is the flux limiter (defined in Sect. 2.2.4), k_{B} is the Boltzmann constant, γ is the ratio of specific heats, μ is the mean molecular weight, m_{H} is the standard mass of one atom of hydrogen, F_{c} is the conductive flux, is the viscous tensor, η is the magnetic resistivity tensor, and J ≡ ∇ × B is the current density (for more information on these quantities, see Mignone et al. 2007, 2012; Orlando et al. 2008).
The opacities κ_{0P}, χ_{0R}, and the radiative losses L_{0} can be provided separately. In our case, we use databases that are precalculated as functions of density and temperature in a nonLTE regime (cf. Sect. 2.2.2).
We note that an additional approximation is made in the radiation energy Eq. (29d) solved by PLUTO. The terms ∇ ⋅ (E_{0} v) and ℙ_{0} : ∇ are missing, compared to Eq. (22a). They were already missing in the LTE RHD version by Kolb et al. (2013, cf. their Eq. (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 freestreaming regime. At this stage, our module therefore 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. NonLTE opacities and radiative power loss: theoretical model
Our module needs the following input data at a given set of density and temperature (ρ, T) and for a given composition of the plasma: the radiation emission L_{0}, the Planck mean opacity κ_{0P}, and the Rosseland mean opacity χ_{0R}.
Databases have been generated by Rodríguez et al. (2018) in a nonLTE regime, but also in the LTE regime. Our module uses three tables (L_{0}, k_{0P}, and k_{0R}) as functions of (ρ, T) for a plasma with solarlike abundances, where k_{0P} and k_{0R} are (cm^{2} g^{−1})
In this section we briefly summarize the features of the theoretical model, whose details are explained in Rodríguez et al. (2018). We 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. We calculated the atomic quantities of the different chemical elements of the multicomponent plasma, such as the relativistic energy levels, wave functions, oscillator strengths and photoionization cross sections, with 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 twice the ionization energy of the ground configuration of the ion.
The atomic level populations were obtained assuming that the plasma is in steady state. This approach is valid when the characteristic time of the most relevant atomic process in the plasma is considerably shorter than the time associated with changes in the plasma density and temperature, that is, 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 that are required to calculate radiative property databases for radiationhydrodynamics simulations. In the steadystate approximation, the population density of the atomic level i of the ion with charge state ζ, denoted N_{ζi}, is obtained by solving the set of rate equations that are implemented in a collisionalradiative steadystate (CRSS) model, given by
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), and take into account all the atomic processes that contribute to populating and depopulating the atomic configuration ζ_{i}. The atomic processes included in the CRSS model were collisional ionization (Lotz 1968) and threebody 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 to be 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 that are available. The formulation developed by Stewart & Pyatt (1966) was applied. In our CRSS model, the ions were considered to be at rest. On the other hand, a MaxwellBoltzmann distribution for the free electrons was assumed when we calculated the rates of the atomic processes. 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 × 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 FermiDirac distribution is not necessary. The CRSS model we described is implemented in the MIXKIP code (Espinosa et al. 2017).
As said before, we simulated a multicomponent plasma. 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 solarlike abundances provided by Asplund et al. (2009) were used. For a given electron density and temperature, the CRSS model was 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 because 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).
After we obtained the atomic data and level populations, the monochromatic absorption coefficient and emissivity of each chemical element m, κ_{m}(ν), and j_{m}(ν), where ν is the photon frequency, were calculated using the RAPCAL code (Rodríguez et al. 2008, 2010). Both coefficients include the boundbound, boundfree, and freefree contributions. The radiative transition rates were calculated in the electric dipole approximation using the FAC code. The oscillator strengths included a correction to take configuration interaction effects into account that are due to the mix between relativistic configurations that belong to the same nonrelativistic 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, electronimpact (Dimitrijevic & Konjevic 1987), and UTA broadenings (Bauche et al. 1987). The lineshape function was applied with the Voigt profile that incorporated all these broadenings. For the freefree contributions, the Kramers semiclassical expression for the inverse bremsstrahlung cross section was used (Rose 1992). In order to determine the opacity, k(ν), 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
where ρ 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} versus free electron density (cm^{−3}) and temperature (K) in logarithmic scale in a nonLTE regime, and in an LTE regime.
Fig. 1. Total radiation emission (erg cm^{−3} s^{−1}) (top panels), Planck opacities (cm^{2} g^{−1}) (mid panels), and Rosseland opacities (cm^{2} g^{−1}) (bottom panels) in a nonLTE regime (left) and in an LTE regime (right). See Sect. 3 for more details. 
Both MIXKIP and RAPCAL codes have been successfully tested with experimental results and numerical simulations for plasmas of single elements that were included in the mixture we analyzed here (Espinosa et al. 2017; Rodríguez et al. 2008, 2010), both in LTE and nonLTE, 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 nonLTE 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, 2013). We used an operatorsplit method. As detailed in this section, advancing the variables during a time step is made in two substeps, an explicit (Sect. 4.2), and an implicit step (Sect. 4.3). The latter involves the radiation emission L_{0}, which is an analytical function of temperature in LTE (L_{0} = k_{0P} ρ c a_{R} T^{4}: cf. Sect. 2.2.2), but has no analytical expression in nonLTE, and is stored in databases (cf. Sect. 3). Our original contribution, with respect to the usual LTE implementations, consists of handling this nonanalytical radiation emission term in the treatment of the implicit scheme, as explained below in Sect. 4.3.
4.1. Reformulation of the equations
For the sake of simplicity, we removed from the full system solved by PLUTO, Eqs. (29a)–(29b), the terms depending on the magnetic field, and rewrote it in a simplified version for an inviscid fluid without heatconduction, but subject to external gravity (even though the following discussion can be applied to the full RMHD system):
where we have omitted the subscript 0 in the radiation quantities for clarity.
The above system contains 11 unknowns: 6 principal variables (ρ,v,p,and E) that are determined by solving the 6 Eqs. (33a)–(33d), and 5 variables (T,ℰ,F) that can be inferred from the 6 principal ones from the FLD relation (33e), the EOS (33f), and the caloric EOS (33g).
At each time step Δ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 (ρ^{n},v^{n},p^{n},and E^{n}) and (T^{n},ℰ^{n},and F^{n}) the values of the variables at time t^{n}, and (ρ^{n + 1},v^{n + 1},p^{n + 1},and E^{n + 1}) and (T^{n + 1},ℰ^{n + 1},and F^{n + 1}) 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 Eqs. (33a)–(33c), but without the radiation source (also known as sourcesink) terms k_{P} ρ c E − L in the gas energy Eq. (33c). The system to be solved is
where we have introduced Eqs. (33e) and (33g) in (33b) and (33c). A Godunovtype algorithm is applied, with several possibilities of Riemann solvers (see Mignone et al. 2007 for details).
We use an asterisk superscript to denote the quantities that are obtained at the completion of this substep. Some of these quantities have intermediate results that will be updated in the next substep (Sect. 4.3). Starting from (ρ^{n},v^{n},and p^{n}), PLUTO determines (ρ^{*},v^{*},and p^{*}), the terms involving the radiation energy density E^{n} that is in the source part of the solved system, and is therefore not updated. The corresponding intermediate temperature is obtained from the EOS (33f):
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 Eq. (33d) and the gas energy Eq. (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 nonLTE is
where the gas energy density ℰ is directly related to the gas temperature through the caloric equation of state (33g).
Because the timescales for radiation are much shorter than timescales for hydrodynamics, time steps that would be inferred from the CFL condition if applied to radiation would lead to impracticable computations. Consequently, we use an implicit scheme for solving the system of Eqs. (36a–36b), with the time step Δt^{n} provided by the explicit hydrodynamics substep 1.
Starting from (E^{n},ρ^{*},T^{*}), we obtain at the completion of this substep, (E^{n + 1},T^{n + 1}). The time discretization of Eqs. (36a) and (36b) is made according to the following implicit scheme:
where while we relate ℰ to ρ and T with (33g), we assume that the mass density ρ is not modified by this implicit substep (following Commerçon et al. 2011b; Zhang et al. 2011, 2013, and Kolb et al. 2013), and where
The radiation diffusion coefficient K, the flux limiter λ, and the dimensionless quantity R are defined in Sect. 2.2.4.
The space discretization of Eq. (37a) is obtained from a finitevolume method (now adapted to nonLTE regimes):
where
represents the discretized radiation diffusion term, a linear function of E^{n + 1} given by Eq. (14) of Kolb et al. (2013), and where the indices i, j, and 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 firstorder Taylor expansion, starting from the state (ρ^{*}, T^{*}) obtained from substep 1:
where we have struck out the term in mass density because we have assumed that the mass density ρ is not modified by this substep 2. Moreover, the dependence of the radiation emission L on T is much higher than that on ρ (see Rodríguez et al. 2018 for more details).
Using Eq. (40a) inside Eqs. (37b) and (39a), we obtain after some algebra (we removed the subscripts i, j, and k from each term for simplicity)
Equation (41a) represents a linear system whose unknown is the radiation energy density , 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; in addition, the boundary conditions for the radiation energy density can be periodic, symmetric, reflective, or with fixed values.
This linearization procedure of a radiative loss function coupled to an implicit scheme was previously implemented by Cunningham et al. (2011, cf. their Eqs. (7) and (8), to be compared to our Eqs. (36b) and (36a), respectively). 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).
After we determine E^{n + 1}, it is straightforward to calculate the temperature T^{n + 1} by simply applying formula (41b). Then, using the approximation ρ^{n + 1} = ρ^{*}, we update the Rosseland opacity , and then the flux limiter λ(R^{n + 1}), therefore the radiation contribution to the righthand side of Eq. (34b) and (34c). We also immediately obtain the pressure p^{n + 1} from the EOS (33f),
The velocity v is not involved in the equations solved by this substep 2. Then, v^{n + 1} = v^{*}. Finally, from the above quantities, we can infer the gas energy ℰ^{n + 1} by applying Eq. (33g) and the radiation flux F^{n + 1} by applying Eq. (33e).
5. Tests
To validate our implementation of the nonLTE equations in the radiation module in PLUTO, we simulated some test cases, and compared our solutions either to analytical or semianalytical 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 these tests, we always used our nonLTE discretization scheme (cf. Sect. 4.3), but imposed the radiation emission (or radiative losses) L to be equal to the LTE emission (L = k_{P} ρ c a_{R} T^{4}). These tests are described below in Sects. 5.1 and 5.2.
The tests are 1D problems. Even though PLUTO works in 1D, 2D, or 3D, the radiation module was developed only in 3D. To model the test cases, we therefore used quasi1D domains, which are cuboids with a length much longer than the width or height. The tests were made only in Cartesian grid. Because 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 blockJacobi (bjacobi) preconditioner.
5.1. Radiation matter coupling
Our objective was to test the correctness of the implementation of Eq. (36b), which couples the evolution of the matter (or gas) internal energy ℰ with the radiation source terms k_{P} ρ c E − L. This equation is solved with an implicit method, in substep 2 of the operatorsplit scheme (cf. Sect. 4.3). In particular, this test enabled us to verify our linearization procedure of the radiation emission L (Eq. (40a)), which 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 wish to freeze the evolution of any other quantity but the temperature (and, therefore, the gas internal energy). To do so, we first built our test model as a static and uniform fluid that is initially out of radiative equilibrium, and we suppressed external gravity. Second, we assumed the radiation energy E to be high enough compared to the gas internal energy for E to be considered as a constant quantity during the energy exchange between matter and radiation throughout the evolution process. In this respect, the radiation flux, which is related to the gradient of radiation energy through the FLD relation of Eq. (23), is negligible, and the continuity, matter momentum, and radiation energy equations (Eqs. (33a), (33b), and (33d)) become irrelevant. The only remaining relevant equation to be solved is the matter energy equation, Eq. (33c), but in the form of Eq. (36b).
Even though we used our nonLTE discretization scheme with the radiation emission (or radiative losses) quantity L in the equations (Sect. 4.3), we set L to be equal to the LTE emission (L = k_{P} ρ c a_{R} T^{4}). In this way, Eq. (36b) can be recast as the following ordinary differential equation:
ℰ_{f} is the gas internal energy in the final radiative equilibrium state.
5.1.1. Setup
For this test, we set the density ρ = 10^{−7} g cm^{−3}, the radiation energy density E = 10^{10} erg cm^{−3}, the Planck^{1} mean opacity ρ k_{P} = 4 × 10^{−8cm} ^{−1}, the mean molecular weight μ = 0.6, and the ratio of specific heats γ = 5/3. The domain consisted of a cuboid whose length far exceeds the other two dimensions. In particular, the cuboid had a width and height of 3 cm, and a length of 100 cm. The grid consisted of 3 × 3 × 100 points. All the boundaries were periodic. For this test the hydrodynamic solver was turned off.
The simulation started at t = 0 s with an initial time step of δt = 10^{−20} s. After each step, the time step increased by 0.1%. The test was performed using three different initial gas internal energy densities ℰ_{0} = (6.4 × 10^{3}; 6.4 × 10^{7}; 6.4 × 10^{8}) erg cm^{−3}.
5.1.2. Results
Figure 2 shows the comparison between the solution found with our model (red dots) and with the semianalytical reference model (black line) obtained by solving Eq. (43a) with a fourthorder RungeKutta scheme (alternatively, a full analytical solution is provided by Swesty & Myra 2009), for the three different initial energy radiation densities ℰ_{0}. The agreement between the two solutions is excellent for all the three cases.
Fig. 2. Coupling test with three different initial gas internal energy densities. The black line represents the reference solution, and the orange dots are the results obtained with the radiation module. 
5.2. Radiative shocks
We tested our implementation of the full set of nonLTE radiation hydrodynamics equations in PLUTO (Eqs. (33a)–(33g)) (without external gravity) by assessing the ability of the code to reproduce the structure of a radiative shock. We simulated a simple shock case in a quasi1D 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 simulated an initially uniform fluid with a density ρ = 7.78 × 10^{−10} g cm^{−3} and temperature T = 10 K that moved with a velocity v along the zaxis. We set the ratio of specific heats γ = 7/5, and the mean molecular weight μ = 1, in analogy to Ensman (1994). We imposed a constant opacity k_{R} × ρ = k_{P} × ρ = 3.1 × 10^{−10} cm^{−1}, and the initial radiation energy density was set by the equation E = a_{R}T^{4}. As in the previous test case (Sect. 5.1), we imposed the radiative losses L = k_{P} ρ c a_{R} T^{4} (LTE radiation emission) and used our nonLTE discretization scheme with L in the equations.
The computational domain had a width and height of 3.418 × 10^{7} cm, and a length (zaxis) of 7 × 10^{10} cm. Following Kolb et al. (2013), we chose a grid composed of 4 × 4 × 2048 cells, the Minerbo flux limiter (Eq. (26b) in Sect. 2.2.4), a LaxFriedrichs scheme, 0.4 as CFL (CourantFriedrichsLewy) value, and a relative tolerance ϵ = 10^{−5} for the matrix solver. The lateral boundaries were periodic. In the direction of the fluid flow (zaxis), we used a reflective boundary at the bottom of the domain (z = 0) to generate the shock, and a zerogradient condition at the top (z = z_{max} = 7 × 10^{10} cm).
Moving from z_{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 (high velocity) (Mihalas & Mihalas 1984). We simulated both cases: we set v = 6 × 10^{5} cm s^{−1} to obtain a subcritical shock, and v = 20 × 10^{5} 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 introduced the quantity s = z − v t. Figure 3 shows the gas temperature (red) and the radiation temperature (blue), T_{rad} = (E/a_{R})^{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.
Fig. 3. Gas temperature (red) and radiation temperature (blue) vs. s = z − v t for the sub (top panel) and supercritical (bottom panel) shocks. The subcritical shock is shown at time t = 3.8 × 10^{4} s and the supercritical shock at t = 7.5 × 10^{3} s. 
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 postshock temperature, reached after radiative cooling. We have
where R_{G} = k_{B}/μm_{H} is the perfect gas constant, σ_{SB} is the StefanBoltzmann 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 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.
Comparisons of the characteristic temperatures of a subcritical shock, T_{2}, T_{−}, T_{+} (Eqs. (44a)–(44c)) obtained from analytical estimates from our nonLTE code and from the initial LTE code by Kolb et al. (2013).
6. LTE versus nonLTE radiative shocks
The purpose of this section is to show the crucial importance of considering the appropriate regime (LTE or nonLTE) 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 the two regimes, as exemplified by Fig. 1. These deviations can have a great effect on the momentum and energy exchanges between matter and radiation and thus on the structure of the flow.
We modified the shock test from Ensman (1994), described in Sect. 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 used the nonLTE discretization scheme with the radiation emission (or radiative losses) quantity L in the equations (Sect. 4.3). In one case, referred to as the “nonLTE case”, we used the nonLTE radiative database (calculated by Rodríguez et al. 2018: cf. Sect. 3), and let the system evolve following the flow conditions, for which at a given time and position, either the nonLTE or the LTE regime prevails and is selfconsistently taken into account by the database. In the other case, referred to as the “LTE case”, we used the LTE database (still calculated by Rodríguez et al. 2018), and therefore force an LTE regime, regardless of the physical conditions.
6.1. Setup
We simulated a uniform fluid with an initial density n = 10^{12} cm^{−3} and temperature T = 2 × 10^{4} K that moves with a velocity v = 5 × 10^{7} cm^{−1} along the zaxis. We set γ = 7/5 and μ ≈ 1.29, that is, we assumed solar abundances. Unlike in the preceding test cases, we did not impose a constant value for k_{P}, k_{R}, and L, but used the radiative databases. The initial radiation energy E was chosen in order to start with a fluid in radiative equilibrium.
The computational domain describes a box of length 10^{8} cm in the nonLTE case, and 10^{7} cm in the LTE case, which are of equal width and height, 3.418 × 10^{7} cm in both cases. The grid is composed of 3 × 3 × 1024 cells.
We adopted the Minerbo flux limiter, and boundary conditions identical to those in the shock tests (Sect. 5.2). For the CFL condition, we used 0.01. For the matrix solver we set a relative tolerance ϵ = 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 and forms a postshock region. Figure 4 shows the profiles of temperature, radiative gains (G = k_{P} ρ c E), which represent the energy gained by the fluid after absorbing radiation, and radiative losses (L) after 6 s of evolution.
Fig. 4. Gas temperature (black dashed line), radiative gains G (red line), and radiative losses L (blue line) vs. z for the LTE (top panel) and nonLTE (bottom panel) shock cases. Both cases are shown after 6 s of evolution. In each plot the gray solid line represents the initial temperature in the domain. 
The two cases present several differences while we follow their evolutions. This is related to the different regimes that are taken into account. In the LTE case (Fig. 4, top panel), the shock reaches a maximum temperature of ∼5 × 10^{4} K (dashed line curve). At this peak, the radiative losses (blue curve: L ∼ 5 × 10^{6} erg cm^{−3} s^{−1}) are higher than the radiative gains (red curve: G ∼ 5 × 10^{5} erg cm^{−3} s^{−1}). As a consequence, the material rapidly cools down until radiative equilibrium (G = L), and the postshock region remains relatively cold at ∼2 × 10^{4} K, close to the initial temperature.
Even though the radiative losses are extremely high compared to the nonLTE 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 that the radiation energy emitted per unit time is not enough to heat up the unshocked plasma; second, according to Fig. 1 (middle right panel), the Planck opacity, k_{P}, in LTE regime is lower by several orders of magnitudes than that in the nonLTE case. Therefore, matter absorbs far less radiation in LTE than in nonLTE (we recall that the gain of radiation energy by matter is G = k_{P} ρ c E).
In the nonLTE regime (Fig. 4, bottom panel), the radiative losses in the shocked region are lower by around three orders of magnitude than those in the LTE regime: L_{LTE} ∼ 5 × 10^{6} erg cm^{−3} s^{−1} (see above), L_{Non − LTE} ∼ 10^{3}−10^{4} erg cm^{−3} s^{−1}. In this case, the shock heats the material up T ∼ 5 × 10^{6} K, and generates a hot postshock region. Because the radiative losses are lower than in the LTE case, the shockheated material needs more time to cool down, and forms a hot slab. After 6 s, the radiative losses trigger the thermal instabilities at the base of the postshock: 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 ∼ 10^{5} K (dashed 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 considers the more general nonLTE 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 nonLTE regimes (including, selfconsistently, the particular case of LTE regime, depending on the physical plasma conditions). We used an operatorsplit method. The system was solved in two substeps, an explicit step for the hyperbolic subsystem, and an implicit step 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 nonLTE conditions.
Starting from the general frequencyintegrated comovingframe radiation hydrodynamics equations, we have reviewed all the assumptions and approximations that have led to the equations that are coded in PLUTO. In particular, we used the fluxlimited diffusion approximation. Moreover, our implementation is valid for plasma in conditions ranging from the freestreaming regime to the 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 nonLTE can be a further improvement of our module. These possible developments will be the subject of future works.
The module needs the following radiative quantities as input data versus density and temperature: the radiation emission, the Planck mean opacity, and the Rosseland mean opacity. Our module currently uses nonLTE databases generated with a CRSS model for a plasma with solarlike abundances. The user may provide any other set of databases that would be more appropriate to the problem to be investigated, however. In our CRSS approach, the absorption processes, such as photoabsorption and photoionization, are ignored in the determination of the atomic level populations. Considering these 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 nonLTE equations, but in LTE conditions, and compared our results with semianalytical 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 the importance of considering the appropriate regime in LTE and in nonLTE regimes by comparing the structure of a radiative shock. This is required 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 to the scientific community upon request at the website of the observatory of Palermo^{2}.
Acknowledgments
The authors are grateful to 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 INAFCINECA (2017)” the CINECA Award HP10B1GLGV and the HPC facility (SCAN) of the INAF – Osservatorio Astronomico di Palermo, for the availability of highperformance computing resources and support. This work was supported by the Programme National de Physique Stellaire (PNPS) of CNRS/INSU cofunded 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 ANR11IDEX000402.
References
 Alme, M. L., & Wilson, J. R. 1973, ApJ, 186, 1015 [NASA ADS] [CrossRef] [Google Scholar]
 Ardila, D. R. 2007, in StarDisk Interaction in Young Stars, eds. J. Bouvier, & I. Appenzeller, IAU Symp., 243, 103 [NASA ADS] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Aubert, D., & Teyssier, R. 2008, MNRAS, 387, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Bauche, J., BaucheArnoult, C., & Klapisch, M. 1987, Adv. At. Mol. Phys., 23, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I. 2004, Radiation Hydrodynamics (West Nyack: Cambridge University Press), 368 [Google Scholar]
 Colombo, S., Ibgui, L., Orlando, S., et al. 2019, A&A, 629, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Hennebelle, P., & Henning, T. 2011a, ApJ, 742, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011b, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Costa, G., Orlando, S., Peres, G., Argiroffi, C., & Bonito, R. 2017, A&A, 597, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cunningham, A. J., Frank, A., & Blackman, E. G. 2006, ApJ, 646, 1059 [NASA ADS] [CrossRef] [Google Scholar]
 Cunningham, A. J., Klein, R. I., Krumholz, M. R., & McKee, C. F. 2011, ApJ, 740, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, S. W., Stone, J. M., & Jiang, Y.F. 2012, ApJS, 199, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, S. W., Jiang, Y.F., Stone, J. M., & Murray, N. 2014, ApJ, 796, 107 [NASA ADS] [CrossRef] [Google Scholar]
 de Sá, L., Chièze, J. P., Stehlé, C., et al. 2019, A&A, 630, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dimitrijevic, M. S., & Konjevic, N. 1987, A&A, 172, 345 [NASA ADS] [Google Scholar]
 Dubroca, B., & Feugeas, J. 1999, Acad. Sci. Paris Comptes Rendus Ser. Sci. Math., 329, 915 [Google Scholar]
 Ensman, L. 1994, ApJ, 424, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Espinosa, G., Rodríguez, R., Gil, J. M., et al. 2017, Phys. Rev. E, 95, 033201 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Gehmeyr, M., & Mihalas, D. 1994, Phys. D Nonlinear Phenom., 77, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Gnedin, N. Y., & Abel, T. 2001, New Astron., 6, 437 [NASA ADS] [CrossRef] [Google Scholar]
 González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 González, M., Vaytet, N., Commerçon, B., & Masson, J. 2015, A&A, 578, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Griem, H. 1997, Principles of Plasma Spectroscopy (Cambridge University Press: Cambridge) [CrossRef] [Google Scholar]
 Gu, M. F. 2008, Can. J. Phys., 86, 675 [NASA ADS] [CrossRef] [Google Scholar]
 Hayes, J. C., & Norman, M. L. 2003, ApJS, 147, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Hirose, S., Krolik, J. H., & Blaes, O. 2009, ApJ, 691, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., & Burrows, A. 2007, ApJ, 659, 1458 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton University Press) [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jiang, Y.F., Stone, J. M., & Davis, S. W. 2012, ApJS, 199, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, Y.F., Stone, J. M., & Davis, S. W. 2014, ApJ, 796, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Klapisch, M., & Busquet, M. 2013, New J. Phys., 15, 015012 [NASA ADS] [CrossRef] [Google Scholar]
 Klassen, M., Kuiper, R., Pudritz, R. E., et al. 2014, ApJ, 797, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
 Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kolb, S. M., Stute, M., Kley, W., & Mignone, A. 2013, A&A, 559, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kramers, H. A. 1923, Mag. J. Sci., 46, 836 [CrossRef] [Google Scholar]
 Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007, ApJ, 667, 626 [NASA ADS] [CrossRef] [Google Scholar]
 Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Krumholz, M. R., Klein, R. I., & McKee, C. F. 2011, ApJ, 740, 74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D. 1984, J. Quant. Spec. Rad. Transf., 31, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Lotz, W. 1968, Z. Phys., 216, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Melon Fuksman, J. D., & Mignone, A. 2019, ApJS, 242, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7 [Google Scholar]
 Mihalas, D., & Auer, L. H. 2001, J. Quant. Spec. Rad. Transf., 71, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., & Klein, R. I. 1982, J. Comput. Phys., 46, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (New York: Oxford University Press) [Google Scholar]
 Minerbo, G. N. 1978, J. Quant. Spec. Rad. Transf., 20, 541 [Google Scholar]
 Orlando, S., Bocchino, F., Reale, F., Peres, G., & Pagano, P. 2008, ApJ, 678, 274 [NASA ADS] [CrossRef] [Google Scholar]
 Pomraning, G. C. 1973, The Equations of Radiation Hydrodynamics (Oxford: Pergamon Press) [Google Scholar]
 Raskutti, S., Ostriker, E. C., & Skinner, M. A. 2016, ApJ, 829, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Rodríguez, R., Florido, R., Gil, J., et al. 2008, Laser Part. Beams, 26, 433 [CrossRef] [Google Scholar]
 Rodríguez, R., Florido, R., Gil, J., et al. 2010, Comput. Phys, 8, 185 [Google Scholar]
 Rodríguez, R., Espinosa, G., & Gil, J. M. 2018, Phys. Rev. E, 98, 033213 [NASA ADS] [CrossRef] [Google Scholar]
 Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188 [NASA ADS] [CrossRef] [Google Scholar]
 Rose, S. J. 1992, J. Phys. B At. Mol. Opt. Phys., 25, 1667 [Google Scholar]
 Rutten, R. J. 1995, Radiative Transfer in Stellar Atmospheres (Utrecht: Sterrekundig Instituut Utrecht) [Google Scholar]
 Sacco, G. G., Argiroffi, C., Orlando, S., et al. 2008, A&A, 491, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Skinner, M. A., & Ostriker, E. C. 2013, ApJS, 206, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Skinner, M. A., & Ostriker, E. C. 2015, ApJ, 809, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Skinner, M. A., Dolence, J. C., Burrows, A., Radice, D., & Vartanyan, D. 2019, ApJS, 241, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Stewart, J. C., & Pyatt, Jr., K. D. 1966, ApJ, 144, 1203 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819 [NASA ADS] [CrossRef] [Google Scholar]
 Swesty, F. D., & Myra, E. S. 2009, ApJS, 181, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]
 van der Holst, B., Tóth, G., Sokolov, I. V., et al. 2011, ApJS, 194, 23 [NASA ADS] [CrossRef] [Google Scholar]
 van Regemorter, H. 1962, ApJ, 136, 906 [NASA ADS] [CrossRef] [Google Scholar]
 Vaytet, N., Audit, E., Chabrier, G., Commerçon, B., & Masson, J. 2012, A&A, 543, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vaytet, N., Chabrier, G., Audit, E., et al. 2013a, A&A, 557, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vaytet, N., González, M., Audit, E., & Chabrier, G. 2013b, J. Quant. Spec. Rad. Transf., 125, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Zel’dovich, Y. B., & Raizer, Y. P. 1967, Physics of Shock Waves and Hightemperature Hydrodynamic Phenomena (Mineola: Dover Publications Inc.) [Google Scholar]
 Zhang, W., Howell, L., Almgren, A., Burrows, A., & Bell, J. 2011, ApJS, 196, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, W., Howell, L., Almgren, A., et al. 2013, ApJS, 204, 7 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Comparisons of the characteristic temperatures of a subcritical shock, T_{2}, T_{−}, T_{+} (Eqs. (44a)–(44c)) obtained from analytical estimates from our nonLTE code and from the initial LTE code by Kolb et al. (2013).
All Figures
Fig. 1. Total radiation emission (erg cm^{−3} s^{−1}) (top panels), Planck opacities (cm^{2} g^{−1}) (mid panels), and Rosseland opacities (cm^{2} g^{−1}) (bottom panels) in a nonLTE regime (left) and in an LTE regime (right). See Sect. 3 for more details. 

In the text 
Fig. 2. Coupling test with three different initial gas internal energy densities. The black line represents the reference solution, and the orange dots are the results obtained with the radiation module. 

In the text 
Fig. 3. Gas temperature (red) and radiation temperature (blue) vs. s = z − v t for the sub (top panel) and supercritical (bottom panel) shocks. The subcritical shock is shown at time t = 3.8 × 10^{4} s and the supercritical shock at t = 7.5 × 10^{3} s. 

In the text 
Fig. 4. Gas temperature (black dashed line), radiative gains G (red line), and radiative losses L (blue line) vs. z for the LTE (top panel) and nonLTE (bottom panel) shock cases. Both cases are shown after 6 s of evolution. In each plot the gray solid line represents the initial temperature in the domain. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.