EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 631 (November 2019) A&A, 631 (2019) A41 Full HTML
Free Access
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/0004-6361/201935991
Published online 23 October 2019

© 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 non-local thermodynamic equilibrium (non-LTE) regimes, the kinetic equilibrium equations, 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 still is 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 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 frequency-integrated, 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 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 (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 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 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 free-streaming 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 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 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 post-shock 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 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. 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 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 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 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 we obtained by letting a given flow evolve in a forced LTE regime and in a non-LTE 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 non-LTE 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 comoving-frame RHD equations

We start with the general comoving-frame frequency-integrated 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)

(1a)

(1b)

(1c)

(1d)

(1e)

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 frequency-integrated, 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 G0 are the time component and the space components of the radiation four-force density vector , respectively, and are defined by (Mihalas & Mihalas 1984)

(2a)

(2b)

where I0(r,t;n0,ν0) is the comoving-frame specific intensity defined at position r, time t, for the comoving-frame direction unit vector n0 and frequency ν0, and where χ0(r,t;n0,ν0) and η0(r,t;n0,ν0) are the extinction coefficient (cm−1) and emission coefficient or emissivity (erg cm−3 s−1 sr−1 Hz−1), respectively,

(3a)

(3b)

where , , , are the comoving-frame thermal absorption coefficient, scattering coefficient, thermal emission coefficient, and scattering emission coefficient, respectively.

In addition, E0, F0, and ℙ0 are the total comoving-frame energy density, the radiation flux, and the radiation pressure:

(4a)

(4b)

(4c)

The radiation moment approach couples the three matter-related 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. Non-LTE 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 four-force density

Our objective is to determine the comoving-frame radiation four-force 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

(5)

Then, the scattering terms in Eq. (2a) cancel out. We also assume isotropic thermal absorption . Finally, we obtain

(6)

where L0 is the total radiation emission per unit volume and time (or radiative power loss),

(7)

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

(8)

where we use the simplified notation .

We focus on G0. We assume the additional assumption of isotropic thermal emission . Then, the emission terms vanish in Eq. (2b). We find

where is the flux-weighted (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,T0) is the Planck function at material temperature T0. As a consequence, using Eq. (8),

(12)

where κ0P is the Planck mean opacity, defined as

(13)

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

(14)

where is the frequency-integrated Planck function (aR: radiation density constant), and where χ0R is the Rosseland mean opacity, defined as

(15)

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 ne, for instance, with a (frequency-independent) 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 non-LTE 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 G0 of the radiation four-force 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, L0, 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

We note that our radiation module never assumes LTE a priori, but always considers from the outset the non-LTE 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:

(17a)

(17b)

(17c)

(17d)

(17e)

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 here apply the same approach to the non-LTE 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 free-streaming 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 fluid-flow timescale tf ≡ ℓ/vf, the typical time for a fluid particle to cross a distance ℓ at characteristic velocity vf, the diffusion timescale td ≡ ℓ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 vd ≡ ℓ/td. We can then quantitatively express the physical distinction, discussed above, between the static diffusion regime, vd ≫ vf (equivalent to td ≪ tf), and the dynamic diffusion regime, vf ≫ vd (equivalent to tf ≪ td). We introduce a characteristic β ≡ vf/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 frequency-integrated optical depth.

Relative sizes of the terms. We analyze the relative sizes of the terms of the comoving-frame 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 td in a static diffusion regime, and the fluid-flow timescale tf in a dynamic diffusion regime. In a free-streaming regime, the characteristic time for a photon to cross ℓ is the radiation-flow timescale tr ≡ ℓ/c.

The three total radiation moments, E0, F0, 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 first-order equilibrium diffusion regime,

In addition, within the second-order equilibrium diffusion approximation (Mihalas & Mihalas 1984; Castor 2004), the comoving-frame radiation energy E0 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 tf for the fluid quantities. We adopt the following characteristic timescales for the radiation quantities: tr in the free-streaming regime, td in the static diffusion regime, and tf 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 non-LTE (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).

Table 1.

Relative sizes of terms in the comoving-frame radiation energy Eq. (17d) (erg cm−3 s−1)

Table 2.

Relative sizes of terms in the comoving-frame radiation momentum Eq. (17e) (dyn cm−3)

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 non-LTE radiation moment equations are

(22a)

(22b)

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. Flux-limited 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:

(23)

where λ (not to be confused with the photon mean free path λp(ν0) defined in Sect. 2.2.3) 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 λ so that

(24)

where the dimensionless quantity R is defined as follows:

(25)

In this way, we recover the asymptotic relations between E0 and F0, 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):

(26a)

(26b)

(26c)

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),

(27)

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

(28)

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. Non-LTE RMHD equations solved by PLUTO

We implemented in a module that is coupled to PLUTO the radiation terms that account for a non-LTE RHD description of a flow by expanding the LTE RHD equations (Kolb et al. 2013) to the more general non-LTE regime. Because PLUTO is a MHD code, we have thus enhanced its capabilities, so that it is now a 3D non-LTE 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

(29a)

(29d)

(29e)

(29f)

(29g)

(29h)

where B is the magnetic field, λ is the flux limiter (defined in Sect. 2.2.4), kB is the Boltzmann constant, γ is the ratio of specific heats, μ is the mean molecular weight, mH is the standard mass of one atom of hydrogen, Fc 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 L0 can be provided separately. In our case, we use databases that are precalculated as functions of density and temperature in a non-LTE 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 ∇ ⋅ (E0 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 free-streaming 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. Non-LTE 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 L0, the Planck mean opacity κ0P, and the Rosseland mean opacity χ0R.

Databases have been generated by Rodríguez et al. (2018) in a non-LTE regime, but also in the LTE regime. Our module uses three tables (L0, k0P, and k0R) as functions of (ρ, T) for a plasma with solar-like abundances, where k0P and k0R are (cm2 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 (L0, k0P, k0R) 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 multi-component 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 radiation-hydrodynamics simulations. In the steady-state 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 collisional-radiative steady-state (CRSS) model, given by

(31)

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 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 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 Maxwell-Boltzmann distribution for the free electrons was assumed when we calculated the rates of the atomic processes. For electron densities between 1011 and 1014 cm−3 and electron temperatures lower than 200 eV, the electron mean free paths range between 3.33 × 105 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 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 solar-like 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 jm(ν), where ν is 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 configuration interaction effects into account that are 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(ν), 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

(32a)

(32b)

where ρ is the mass density of the mixture. From the monochromatic opacities of the mixture, the Planck k0P and Rosseland k0R mean opacities were calculated, as well as the radiative power loss, or radiation emission, L0, from the monochromatic emissivity. Figure 1 portrays maps of L0, k0P, and k0R versus free electron density (cm−3) and temperature (K) in logarithmic scale in a non-LTE regime, and in an LTE regime.

thumbnail Fig. 1.

Total radiation emission (erg cm−3 s−1) (top panels), Planck opacities (cm2 g−1) (mid panels), and Rosseland opacities (cm2 g−1) (bottom panels) in a non-LTE 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 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, 2013). We used an operator-split 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 L0, which is an analytical function of temperature in LTE (L0 = k0P ρ c aR T4: cf. Sect. 2.2.2), but has no analytical expression in non-LTE, and is stored in databases (cf. Sect. 3). Our original contribution, with respect to the usual LTE implementations, consists of handling this non-analytical 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 heat-conduction, but subject to external gravity (even though the following discussion can be applied to the full RMHD system):

(33a)

(33b)

(33c)

(33d)

(33e)

(33f)

(33g)

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 Δtn = tn + 1 − tn, determined by the CFL condition (Mignone et al. 2007), we advance the gas and radiation variables from tn to tn + 1 by solving the above system in two consecutive substeps, as described below. We denote as (ρn,vn,pn,and En) and (Tn,ℰn,and Fn) the values of the variables at time tn, and (ρn + 1,vn + 1,pn + 1,and En + 1) and (Tn + 1,ℰn + 1,and Fn + 1) their values at time tn + 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 source-sink) terms kP ρ 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 Godunov-type 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,vn,and pn), PLUTO determines (ρ*,v*,and p*), the terms involving the radiation energy density En 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):

(35)

4.3. Substep 2: implicit radiation diffusion and source terms

In this second substep, we determine at time tn + 1 the radiation energy density En + 1 and the gas temperature Tn + 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 non-LTE 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 Δtn provided by the explicit hydrodynamics substep 1.

Starting from (En,ρ*,T*), we obtain at the completion of this substep, (En + 1,Tn + 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 finite-volume method (now adapted to non-LTE regimes):

(39a)

where

(39b)

represents the discretized radiation diffusion term, a linear function of En + 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 tn + 1, Ln + 1 (Eq. (38b)), is determined from a first-order 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 tn + 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 En + 1, it is straightforward to calculate the temperature Tn + 1 by simply applying formula (41b). Then, using the approximation ρn + 1 = ρ*, we update the Rosseland opacity , and then the flux limiter λ(Rn + 1), therefore the radiation contribution to the right-hand side of Eq. (34b) and (34c). We also immediately obtain the pressure pn + 1 from the EOS (33f),

(42)

The velocity v is not involved in the equations solved by this substep 2. Then, vn + 1 = v*. Finally, from the above quantities, we can infer the gas energy ℰn + 1 by applying Eq. (33g) and the radiation flux Fn + 1 by applying Eq. (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 these tests, we always used our non-LTE discretization scheme (cf. Sect. 4.3), but imposed the radiation emission (or radiative losses) L to be equal to the LTE emission (L = kP ρ c aR T4). 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 quasi-1D 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 block-Jacobi (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 kP ρ c E − L. This equation is solved with an implicit method, in substep 2 of the operator-split 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 Tn + 1 versus Tn (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 non-LTE 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 = kP ρ c aR T4). 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 = 1010 erg cm−3, the Planck1 mean opacity ρ kP = 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 × 103;  6.4 × 107;  6.4 × 108) erg cm−3.

5.1.2. Results

Figure 2 shows the comparison between the solution found with our model (red dots) and with the semi-analytical reference model (black line) obtained by solving Eq. (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 ℰ0. The agreement between the two solutions is excellent for all the three cases.

thumbnail 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 non-LTE 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 quasi-1D 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 z-axis. 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 kR × ρ = kP × ρ = 3.1 × 10−10 cm−1, and the initial radiation energy density was set by the equation E = aRT4. As in the previous test case (Sect. 5.1), we imposed the radiative losses L = kP ρ c aR T4 (LTE radiation emission) and used our non-LTE discretization scheme with L in the equations.

The computational domain had a width and height of 3.418 × 107 cm, and a length (z-axis) of 7 × 1010 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 Lax-Friedrichs scheme, 0.4 as CFL (Courant-Friedrichs-Lewy) value, and a relative tolerance ϵ = 10−5 for the matrix solver. The lateral boundaries were periodic. In the direction of the fluid flow (z-axis), we used a reflective boundary at the bottom of the domain (z = 0) to generate the shock, and a zero-gradient condition at the top (z = zmax = 7 × 1010 cm).

Moving from zmax 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 × 105 cm s−1 to obtain a subcritical shock, and v = 20 × 105 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), Trad = (E/aR)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.

thumbnail 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 × 104 s and the supercritical shock at t = 7.5 × 103 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); T2, the final equilibrium post-shock temperature, reached after radiative cooling. We have

where RG = kB/μmH is the perfect gas constant, σ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 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.

Table 3.

Comparisons of the characteristic temperatures of a subcritical shock, T2, T−, T+ (Eqs. (44a)–(44c)) obtained from analytical estimates from our non-LTE code and from the initial LTE code by Kolb et al. (2013).

6. LTE versus non-LTE radiative shocks

The purpose of this section is to show the crucial importance of considering the appropriate regime (LTE or 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 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 non-LTE discretization scheme with the radiation emission (or radiative losses) quantity L in the equations (Sect. 4.3). In one case, referred to as the “non-LTE case”, we used the non-LTE 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 non-LTE or the LTE regime prevails and is self-consistently 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 = 1012 cm−3 and temperature T = 2 × 104 K that moves with a velocity v = 5 × 107 cm−1 along the z-axis. 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 kP, kR, 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 108 cm in the non-LTE case, and 107 cm in the LTE case, which are of equal width and height, 3.418 × 107 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 post-shock region. Figure 4 shows the profiles of temperature, radiative gains (G = kP ρ c E), which represent the energy gained by the fluid after absorbing radiation, and radiative losses (L) after 6 s of evolution.

thumbnail 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 non-LTE (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 × 104 K (dashed line curve). At this peak, the radiative losses (blue curve: L ∼ 5 × 106 erg cm−3 s−1) are higher than the radiative gains (red curve: G ∼ 5 × 105 erg cm−3 s−1). As a consequence, the material rapidly cools down until radiative equilibrium (G = L), and the post-shock region remains relatively cold at ∼2 × 104 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 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, kP, in LTE regime is lower by several orders of magnitudes than that in the non-LTE case. Therefore, matter absorbs far less radiation in LTE than in non-LTE (we recall that the gain of radiation energy by matter is G = kP ρ c E).

In the non-LTE 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: LLTE ∼ 5 × 106 erg cm−3 s−1 (see above), LNon − LTE ∼ 103−104 erg cm−3 s−1. In this case, the shock heats the material up T ∼ 5 × 106 K, and generates a hot post-shock region. Because the radiative losses are lower than in the LTE case, the shock-heated 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 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 ∼ 105 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 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 used an operator-split 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 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 coded in PLUTO. In particular, we used the flux-limited diffusion approximation. Moreover, our implementation is valid for plasma in conditions ranging from the free-streaming 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 non-LTE 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 non-LTE databases generated with a CRSS model for a plasma with solar-like 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 non-LTE equations, but in LTE conditions, and 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 the importance of considering the appropriate regime in LTE and in non-LTE 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 Palermo2.

1

The choice of the Rosseland mean opacity is irrelevant because it is associated with the radiation flux that is negligible in our test; however, for computational convenience, we also adopted ρkR = 4 × 10−8 cm−1.

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 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 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 ANR-11-IDEX-0004-02.

References

  1. Alme, M. L., & Wilson, J. R. 1973, ApJ, 186, 1015 [NASA ADS] [CrossRef] [Google Scholar]
  2. Ardila, D. R. 2007, in Star-Disk Interaction in Young Stars, eds. J. Bouvier, & I. Appenzeller, IAU Symp., 243, 103 [NASA ADS] [Google Scholar]
  3. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  4. Aubert, D., & Teyssier, R. 2008, MNRAS, 387, 295 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bauche, J., Bauche-Arnoult, C., & Klapisch, M. 1987, Adv. At. Mol. Phys., 23, 131 [NASA ADS] [CrossRef] [Google Scholar]
  6. Castor, J. I. 2004, Radiation Hydrodynamics (West Nyack: Cambridge University Press), 368 [Google Scholar]
  7. Colombo, S., Ibgui, L., Orlando, S., et al. 2019, A&A, 629, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Commerçon, B., Hennebelle, P., & Henning, T. 2011a, ApJ, 742, L9 [NASA ADS] [CrossRef] [Google Scholar]
  9. Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011b, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Costa, G., Orlando, S., Peres, G., Argiroffi, C., & Bonito, R. 2017, A&A, 597, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Cunningham, A. J., Frank, A., & Blackman, E. G. 2006, ApJ, 646, 1059 [NASA ADS] [CrossRef] [Google Scholar]
  12. Cunningham, A. J., Klein, R. I., Krumholz, M. R., & McKee, C. F. 2011, ApJ, 740, 107 [NASA ADS] [CrossRef] [Google Scholar]
  13. Davis, S. W., Stone, J. M., & Jiang, Y.-F. 2012, ApJS, 199, 9 [NASA ADS] [CrossRef] [Google Scholar]
  14. Davis, S. W., Jiang, Y.-F., Stone, J. M., & Murray, N. 2014, ApJ, 796, 107 [NASA ADS] [CrossRef] [Google Scholar]
  15. de Sá, L., Chièze, J. P., Stehlé, C., et al. 2019, A&A, 630, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Dimitrijevic, M. S., & Konjevic, N. 1987, A&A, 172, 345 [NASA ADS] [Google Scholar]
  17. Dubroca, B., & Feugeas, J. 1999, Acad. Sci. Paris Comptes Rendus Ser. Sci. Math., 329, 915 [Google Scholar]
  18. Ensman, L. 1994, ApJ, 424, 275 [NASA ADS] [CrossRef] [Google Scholar]
  19. Espinosa, G., Rodríguez, R., Gil, J. M., et al. 2017, Phys. Rev. E, 95, 033201 [NASA ADS] [CrossRef] [Google Scholar]
  20. Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131 [NASA ADS] [CrossRef] [Google Scholar]
  22. Gehmeyr, M., & Mihalas, D. 1994, Phys. D Nonlinear Phenom., 77, 320 [NASA ADS] [CrossRef] [Google Scholar]
  23. Gnedin, N. Y., & Abel, T. 2001, New Astron., 6, 437 [NASA ADS] [CrossRef] [Google Scholar]
  24. González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. González, M., Vaytet, N., Commerçon, B., & Masson, J. 2015, A&A, 578, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Griem, H. 1997, Principles of Plasma Spectroscopy (Cambridge University Press: Cambridge) [CrossRef] [Google Scholar]
  27. Gu, M. F. 2008, Can. J. Phys., 86, 675 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hayes, J. C., & Norman, M. L. 2003, ApJS, 147, 197 [NASA ADS] [CrossRef] [Google Scholar]
  29. Hirose, S., Krolik, J. H., & Blaes, O. 2009, ApJ, 691, 16 [NASA ADS] [CrossRef] [Google Scholar]
  30. Hubeny, I., & Burrows, A. 2007, ApJ, 659, 1458 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton University Press) [Google Scholar]
  32. Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2012, ApJS, 199, 14 [NASA ADS] [CrossRef] [Google Scholar]
  34. Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014, ApJ, 796, 106 [NASA ADS] [CrossRef] [Google Scholar]
  35. Klapisch, M., & Busquet, M. 2013, New J. Phys., 15, 015012 [NASA ADS] [CrossRef] [Google Scholar]
  36. Klassen, M., Kuiper, R., Pudritz, R. E., et al. 2014, ApJ, 797, 4 [NASA ADS] [CrossRef] [Google Scholar]
  37. Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
  38. Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Kolb, S. M., Stute, M., Kley, W., & Mignone, A. 2013, A&A, 559, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Kramers, H. A. 1923, Mag. J. Sci., 46, 836 [CrossRef] [Google Scholar]
  41. Krumholz, M. R., Klein, R. I., McKee, C. F., & Bolstad, J. 2007, ApJ, 667, 626 [NASA ADS] [CrossRef] [Google Scholar]
  42. 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]
  43. Krumholz, M. R., Klein, R. I., & McKee, C. F. 2011, ApJ, 740, 74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71 [NASA ADS] [CrossRef] [Google Scholar]
  45. Levermore, C. D. 1984, J. Quant. Spec. Rad. Transf., 31, 149 [NASA ADS] [CrossRef] [Google Scholar]
  46. Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
  47. Lotz, W. 1968, Z. Phys., 216, 241 [NASA ADS] [CrossRef] [Google Scholar]
  48. Melon Fuksman, J. D., & Mignone, A. 2019, ApJS, 242, 20 [NASA ADS] [CrossRef] [Google Scholar]
  49. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
  50. Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7 [Google Scholar]
  51. Mihalas, D., & Auer, L. H. 2001, J. Quant. Spec. Rad. Transf., 71, 61 [NASA ADS] [CrossRef] [Google Scholar]
  52. Mihalas, D., & Klein, R. I. 1982, J. Comput. Phys., 46, 97 [NASA ADS] [CrossRef] [Google Scholar]
  53. Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (New York: Oxford University Press) [Google Scholar]
  54. Minerbo, G. N. 1978, J. Quant. Spec. Rad. Transf., 20, 541 [Google Scholar]
  55. Orlando, S., Bocchino, F., Reale, F., Peres, G., & Pagano, P. 2008, ApJ, 678, 274 [NASA ADS] [CrossRef] [Google Scholar]
  56. Pomraning, G. C. 1973, The Equations of Radiation Hydrodynamics (Oxford: Pergamon Press) [Google Scholar]
  57. Raskutti, S., Ostriker, E. C., & Skinner, M. A. 2016, ApJ, 829, 130 [NASA ADS] [CrossRef] [Google Scholar]
  58. Rodríguez, R., Florido, R., Gil, J., et al. 2008, Laser Part. Beams, 26, 433 [CrossRef] [Google Scholar]
  59. Rodríguez, R., Florido, R., Gil, J., et al. 2010, Comput. Phys, 8, 185 [Google Scholar]
  60. Rodríguez, R., Espinosa, G., & Gil, J. M. 2018, Phys. Rev. E, 98, 033213 [NASA ADS] [CrossRef] [Google Scholar]
  61. Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188 [NASA ADS] [CrossRef] [Google Scholar]
  62. Rose, S. J. 1992, J. Phys. B At. Mol. Opt. Phys., 25, 1667 [Google Scholar]
  63. Rutten, R. J. 1995, Radiative Transfer in Stellar Atmospheres (Utrecht: Sterrekundig Instituut Utrecht) [Google Scholar]
  64. Sacco, G. G., Argiroffi, C., Orlando, S., et al. 2008, A&A, 491, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Skinner, M. A., & Ostriker, E. C. 2013, ApJS, 206, 21 [NASA ADS] [CrossRef] [Google Scholar]
  67. Skinner, M. A., & Ostriker, E. C. 2015, ApJ, 809, 187 [NASA ADS] [CrossRef] [Google Scholar]
  68. Skinner, M. A., Dolence, J. C., Burrows, A., Radice, D., & Vartanyan, D. 2019, ApJS, 241, 7 [NASA ADS] [CrossRef] [Google Scholar]
  69. Stewart, J. C., & Pyatt, Jr., K. D. 1966, ApJ, 144, 1203 [NASA ADS] [CrossRef] [Google Scholar]
  70. Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819 [NASA ADS] [CrossRef] [Google Scholar]
  71. Swesty, F. D., & Myra, E. S. 2009, ApJS, 181, 1 [NASA ADS] [CrossRef] [Google Scholar]
  72. Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]
  73. van der Holst, B., Tóth, G., Sokolov, I. V., et al. 2011, ApJS, 194, 23 [NASA ADS] [CrossRef] [Google Scholar]
  74. van Regemorter, H. 1962, ApJ, 136, 906 [NASA ADS] [CrossRef] [Google Scholar]
  75. Vaytet, N., Audit, E., Chabrier, G., Commerçon, B., & Masson, J. 2012, A&A, 543, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Vaytet, N., Chabrier, G., Audit, E., et al. 2013a, A&A, 557, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Vaytet, N., González, M., Audit, E., & Chabrier, G. 2013b, J. Quant. Spec. Rad. Transf., 125, 105 [NASA ADS] [CrossRef] [Google Scholar]
  78. Zel’dovich, Y. B., & Raizer, Y. P. 1967, Physics of Shock Waves and High-temperature Hydrodynamic Phenomena (Mineola: Dover Publications Inc.) [Google Scholar]
  79. Zhang, W., Howell, L., Almgren, A., Burrows, A., & Bell, J. 2011, ApJS, 196, 20 [NASA ADS] [CrossRef] [Google Scholar]
  80. Zhang, W., Howell, L., Almgren, A., et al. 2013, ApJS, 204, 7 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Relative sizes of terms in the comoving-frame radiation energy Eq. (17d) (erg cm−3 s−1)

In the text
Table 2.

Relative sizes of terms in the comoving-frame radiation momentum Eq. (17e) (dyn cm−3)

In the text
Table 3.

Comparisons of the characteristic temperatures of a subcritical shock, T2, T−, T+ (Eqs. (44a)–(44c)) obtained from analytical estimates from our non-LTE code and from the initial LTE code by Kolb et al. (2013).

In the text

All Figures

thumbnail Fig. 1.

Total radiation emission (erg cm−3 s−1) (top panels), Planck opacities (cm2 g−1) (mid panels), and Rosseland opacities (cm2 g−1) (bottom panels) in a non-LTE regime (left) and in an LTE regime (right). See Sect. 3 for more details.
In the text
thumbnail 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
thumbnail 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 × 104 s and the supercritical shock at t = 7.5 × 103 s.
In the text
thumbnail 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 non-LTE (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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.