A Reduced-order NLTE Kinetic Model for Radiating Plasmas of Outer Envelopes of Stellar Atmospheres

The material gas under investigation is made up of three chemical components: free electrons ${e}^{-}$, hydrogen atoms ${\rm{H}}$, and protons ${{\rm{H}}}^{+}$. These components constitute the set ${ \mathcal C }=\{{e}^{-},{\rm{H}},{{\rm{H}}}^{+}\}$. The heavy-particle components (i.e., hydrogen atoms and protons) are stored in the subset ${{ \mathcal C }}_{{\rm{h}}}=\{{\rm{H}},{{\rm{H}}}^{+}\}$. The hydrogen bound states (sorted by increasing energy) are denoted by the global index i and stored in set ${ \mathcal I }$. Based on an StS approach, the bound states ${\rm{H}}(i)$ are treated as separate pseudo-species (Bates et al. 1962a, 1962b; Capitelli et al. 2000). The related statistical weights and formation energies are indicated, respectively, by the symbols ${g}_{i}$ and ${E}_{i}$. The corresponding quantities for the protons are ${g}_{+}\,=\,1$ and ${E}_{+}=13.59\,\mathrm{eV}$. The StS species are stored in the species set ${ \mathcal S }=\{{e}^{-},{\rm{H}}(i),{{\rm{H}}}^{+};i\in { \mathcal I }\}$.

In order to make the problem tractable, simplifying assumptions are introduced. For the adopted pressure and temperature conditions, it is possible to use Boltzmann statistics. Each species is modeled as a thermally perfect gas by disregarding collective plasma effects such as pressure ionization and multiple charge–charge interactions (Rogers 1974; Rogers et al. 1996; Hubeny & Mihalas 2014). Further, the velocity distribution function of each species is taken to be Maxwellian at its own translational temperature. For heavy particles, this temperature is taken to be a common heavy-particle temperature ${T}_{{\rm{h}}}$. This assumption is justified in light of (i) the efficient energy transfer in collisions between particles with similar masses and (ii) the large cross-section for resonant charge transfer in hydrogen–proton collisions (i.e., ${\rm{H}}+{{\rm{H}}}^{+}={{\rm{H}}}^{+}+{\rm{H}}$) (Mihalas & Mihalas 1984, Chapter 8). On the other hand, collisions between heavy particles and free electrons are energetically inefficient, due to the large mass disparity. This motivates the introduction of a separate free-electron temperature ${T}_{{\rm{e}}}$. The adopted StS model is therefore a two-temperature plasma model where thermal non-equilibrium effects between heavy particles and free electrons are taken into account through the macroscopic parameters ${T}_{{\rm{h}}}$ and ${T}_{{\rm{e}}}$, respectively (Appleton & Bray 1964; Murty 1971). In this work, it is further assumed that the plasma is locally neutral. Electromagnetic fields and transport phenomena are also neglected.

In view of the assumptions introduced in the Section 2.1, the velocity distribution functions of free electrons, hydrogen bound states, and protons are

$$\begin{eqnarray}&&{f}_{{\rm{e}}}({{\boldsymbol{c}}}_{{\rm{e}}})={n}_{{\rm{e}}}{\left(\displaystyle \frac{{m}_{{\rm{e}}}}{2\pi {k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right)}^{3/2}\exp \left(-\displaystyle \frac{{m}_{{\rm{e}}}| {{\boldsymbol{c}}}_{{\rm{e}}}-{\boldsymbol{v}}{| }^{2}}{2{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{f}_{i}({{\boldsymbol{c}}}_{i})={n}_{i}{\left(\displaystyle \frac{{m}_{{\rm{H}}}}{2\pi {k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right)}^{3/2}\exp \left(-\displaystyle \frac{{m}_{{\rm{H}}}| {{\boldsymbol{c}}}_{i}-{\boldsymbol{v}}{| }^{2}}{2{k}_{{\rm{B}}}{T}_{{\rm{h}}}}\right),\quad i\in { \mathcal I },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{f}_{+}({{\boldsymbol{c}}}_{+})={n}_{+}{\left(\displaystyle \frac{{m}_{{{\rm{H}}}^{+}}}{2\pi {k}_{{\rm{B}}}{T}_{{\rm{h}}}}\right)}^{3/2}\exp \left(-\displaystyle \frac{{m}_{{{\rm{H}}}^{+}}| {{\boldsymbol{c}}}_{+}-{\boldsymbol{v}}{| }^{2}}{2{k}_{{\rm{B}}}{T}_{{\rm{h}}}}\right),\end{eqnarray} \tag{ 3 }$$

where the quantities n_s, $s\in { \mathcal S }$, and m_c, $c\in { \mathcal C }$, denote, respectively, the species number densities and the mass of the chemical components. The symbols ${k}_{{\rm{B}}}$ stands for Boltzmann's constant, whereas the vectors ${{\boldsymbol{c}}}_{s}$, $s\in {\boldsymbol{s}}$, and ${\boldsymbol{v}}$ represent the molecular velocities and the bulk velocity of the material gas, respectively.

Thermodynamic properties (e.g., pressure, energy density) of the single species or the whole mixture can be computed based on suitable moments of the velocity distribution functions (Equations (1)–(3)) (Ferziger & Kaper 1972, Chapter 1). The gas pressure is obtained by taking the average of the sum of the diagonal terms of the pressure tensor and reads

$$\begin{eqnarray}p & = & \displaystyle \frac{1}{3}{\displaystyle \int }_{-\infty }^{+\infty }{m}_{{\rm{e}}}\,{{\boldsymbol{C}}}_{{\rm{e}}}\cdot {{\boldsymbol{C}}}_{{\rm{e}}}\,{f}_{{\rm{e}}}({{\boldsymbol{c}}}_{{\rm{e}}})\,d{{\boldsymbol{c}}}_{{\rm{e}}}\\ & & +\displaystyle \sum _{i\in { \mathcal I }}\displaystyle \frac{1}{3}{\displaystyle \int }_{-\infty }^{+\infty }{m}_{{\rm{H}}}\,{{\boldsymbol{C}}}_{i}\cdot {{\boldsymbol{C}}}_{i}\,{f}_{i}({{\boldsymbol{c}}}_{i})\,d{{\boldsymbol{c}}}_{i}\\ & & +\displaystyle \frac{1}{3}{\displaystyle \int }_{-\infty }^{+\infty }{m}_{{{\rm{H}}}^{+}}\,{{\boldsymbol{C}}}_{+}\cdot {{\boldsymbol{C}}}_{+}\,{f}_{+}({{\boldsymbol{c}}}_{+})\,d{{\boldsymbol{c}}}_{+},\end{eqnarray} \tag{ 4 }$$

where the peculiar velocities are ${{\boldsymbol{C}}}_{s}={{\boldsymbol{c}}}_{s}-{\boldsymbol{v}}$, $s\in { \mathcal S }$. The evaluation of the integrals in Equation (4) leads to Dalton's law for a thermal non-equilibrium plasma:

$$\begin{eqnarray}&&p={p}_{{\rm{e}}}+{p}_{{\rm{h}}},\end{eqnarray} \tag{ 5 }$$

where the free-electron and heavy-particle partial pressures are, respectively, ${p}_{{\rm{e}}}={n}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}}$ and ${p}_{{\rm{h}}}={n}_{{\rm{h}}}{k}_{{\rm{B}}}{T}_{{\rm{h}}}$. Quantity ${n}_{{\rm{h}}}$ stands for the heavy-particle number density and is obtained by summing the contributions of the hydrogen bound states and protons, ${n}_{{\rm{h}}}={\sum }_{i\in { \mathcal I }}{n}_{i}+{n}_{+}$, where ${n}_{+}\,=\,{n}_{{\rm{e}}}$ in light of the assumed charge neutrality (see Section 2.1).

The heavy-particle and free-electron thermal energy densities are defined as

$$\begin{eqnarray}\rho {e}_{{\rm{h}}} & = & \displaystyle \sum _{i\in { \mathcal I }}{\displaystyle \int }_{-\infty }^{+\infty }\left(\displaystyle \frac{1}{2}{m}_{{\rm{H}}}\,{{\boldsymbol{C}}}_{i}\cdot {{\boldsymbol{C}}}_{i}+{E}_{i}\right){f}_{i}({{\boldsymbol{c}}}_{i})\,d{{\boldsymbol{c}}}_{i}\\ & & +{\displaystyle \int }_{-\infty }^{+\infty }\left(\displaystyle \frac{1}{2}{m}_{{{\rm{H}}}^{+}}\,{{\boldsymbol{C}}}_{+}\cdot {{\boldsymbol{C}}}_{+}+{E}_{+}\right){f}_{+}({{\boldsymbol{c}}}_{+})\,d{{\boldsymbol{c}}}_{+},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\rho {e}_{{\rm{e}}}={\int }_{-\infty }^{+\infty }\displaystyle \frac{1}{2}{m}_{{\rm{e}}}\,{{\boldsymbol{C}}}_{{\rm{e}}}\cdot {{\boldsymbol{C}}}_{{\rm{e}}}\,{f}_{{\rm{e}}}({{\boldsymbol{c}}}_{{\rm{e}}})\,d{{\boldsymbol{c}}}_{{\rm{e}}}.\end{eqnarray} \tag{ 7 }$$

Using Equations (4) and (5), the integrals (6) and (7) can be evaluated in a straightforward manner to give

$$\begin{eqnarray}&&\rho {e}_{{\rm{h}}}=\displaystyle \frac{3}{2}{p}_{{\rm{h}}}+\displaystyle \sum _{i\in { \mathcal I }}\,{n}_{i}{E}_{i}+{n}_{+}{E}_{+},\quad \rho {e}_{{\rm{e}}}=\displaystyle \frac{3}{2}{p}_{{\rm{e}}}.\end{eqnarray} \tag{ 8 }$$

The gas total energy density is obtained by summing the thermal contribution from free electrons and heavy particles, and the kinetic contribution of the gas as a whole, $\rho E=\rho {e}_{{\rm{h}}}+\rho {e}_{{\rm{e}}}+\rho {\boldsymbol{v}}\cdot {\boldsymbol{v}}/2$. The mass density is evaluated based on the species number densities as $\rho ={\sum }_{s\in { \mathcal S }}{n}_{s}\,{m}_{s}$. The total enthalpy density follows by adding the pressure to the total energy density, $\rho H=\rho E+p$.

Under complete thermodynamic equilibrium, all temperatures are equal (i.e., ${T}_{{\rm{h}}}={T}_{{\rm{e}}}=T$) and the population of hydrogen bound states is given by the Boltzmann distribution law:

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{i}}{{n}_{{\rm{H}}}}=\displaystyle \frac{{g}_{i}}{{Z}_{{\rm{H}}}^{{\rm{i}}}(T)}\exp \left(-\displaystyle \frac{{E}_{i}}{{k}_{{\rm{B}}}T}\right),\end{eqnarray} \tag{ 9 }$$

$i\in { \mathcal I }$, where the hydrogen number density and internal partition function are ${n}_{{\rm{H}}}={\sum }_{i\in { \mathcal I }}{n}_{i}$ and ${Z}_{{\rm{H}}}^{{\rm{i}}}(T)\ ={\sum }_{i\in { \mathcal I }}{g}_{i}\exp (-{E}_{i}/{k}_{{\rm{B}}}T)$, respectively. The hydrogen number density, ${n}_{{\rm{H}}}$, is related to that of free electrons and protons by Saha's equation and charge conservation (i.e., ${n}_{{\rm{e}}}={n}_{+}$):

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{{\rm{e}}}\,{n}_{+}}{{n}_{{\rm{H}}}}=\displaystyle \frac{{g}_{{\rm{e}}}\,{Z}_{{\rm{e}}}^{{\rm{t}}}(T)\,{g}_{+}\,{Z}_{{{\rm{H}}}^{+}}^{{\rm{t}}}(T)}{{Z}_{{\rm{H}}}^{{\rm{t}}}(T)\,{Z}_{{\rm{H}}}^{{\rm{i}}}(T)}\exp \left(-\displaystyle \frac{{E}_{+}}{{k}_{{\rm{B}}}T}\right),\end{eqnarray} \tag{ 10 }$$

where the translational partition functions per unit volume are ${Z}_{c}^{{\rm{t}}}{(T)=(2\pi {m}_{c}{k}_{{\rm{B}}}T)}^{3/2}/{h}_{{\rm{P}}}^{3}$, $c\in { \mathcal C }$, where ${h}_{{\rm{P}}}$ is Planck's constant. The quantity ${g}_{{\rm{e}}}=2$ is the electron degeneracy accounting for its spin. Under a general NLTE situation, Equations (9) and (10) no longer hold. However, as shown in the results, Section 5, in the internal relaxation region behind the shock, the population of high-lying hydrogen bound states can be approximately described by a Boltzmann distribution (9) at its own internal temperature before equilibrium is reached. This motivates the introduction of reduced-order models, treated in Section 3.1, which represent the main subject of the present work.

Before ending the present discussion on thermodynamics, it should be recalled that, when studying the interaction between matter and radiation, one should in general account for the momentum and energy content of photons (Oxenius 1986, Chapter 7). The former plays an important role, for instance, in stellar interiors where up to 50% of the pressure may be due to radiation (Rogers et al. 1996). However, as shown in the work by Fadeyev & Gillet (1998, 2000), the radiant pressure and energy density are negligible compared to the corresponding gas values for the conditions adopted in this work.

The present section describes the NLTE kinetic model for atomic hydrogen plasmas used in this work. This StS model provides the basis for the reduced-order models developed in Section 3.1 and accounts for both collisional and radiative processes. The following subsections provide a description of their modeling. It is worth mentioning that, in the literature, other StS models for hydrogen plasmas exist such as the one developed by D'Ammando et al. (2010) and Colonna et al. (2012).

The calculation of the mass production and energy transfer terms due to collisional and radiative processes is obtained using the velocity distribution functions (1)–(3) in the moments of the Boltzmann equation for gas particles (for the details, the reader may consult Chapter 4 of the book by Oxenius 1986).

2.3.1. Collisional Processes

Collisional processes account for inelastic/elastic transitions, due to collisions between heavy particles and free electrons and comprise

• (i)  
  excitation (e) and ionization (i) by electron impact:

  $$\begin{eqnarray}{\rm{H}}(i)+{e}^{-} & \leftrightarrows & {\rm{H}}(j)+{e}^{-},\quad i\lt j,\\ {\rm{H}}(i)+{e}^{-} & \leftrightarrows & {{\rm{H}}}^{+}+{e}^{-}+{e}^{-},\quad i,j\in { \mathcal I },\end{eqnarray} \tag{ 11 }$$
• (ii)  
  elastic (el) energy exchange between heavy particles and free electrons.

Excitation and ionization by heavy-particle impact have been neglected since, as shown in Section 5, the number of free electrons produced by photo-ionization in the precursor guarantees that the collisional kinetics is dominated by free electrons in the internal relaxation region.

The species production terms due to excitation and ionization are obtained by taking the zeroth-order moment (i.e., mass) of the collision operators for excitation and ionization in the Boltzmann equation for free electrons, hydrogen bound states, and protons (Oxenius 1986, Chapter 2):

$$\begin{eqnarray}&&{\omega }_{{\rm{e}}}^{\mathrm{col}}=\displaystyle \sum _{i\in { \mathcal I }}{n}_{{\rm{e}}}\,[{n}_{i}\,{k}_{i}^{{\rm{I}}}({T}_{{\rm{e}}})-{n}_{{\rm{e}}}\,{n}_{+}\,{k}_{i}^{{\rm{R}}}({T}_{{\rm{e}}})],\quad {\omega }_{+}^{\mathrm{col}}={\omega }_{{\rm{e}}}^{\mathrm{col}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}{\omega }_{i}^{\mathrm{col}} & = & -\,{n}_{{\rm{e}}}\,[{n}_{i}\,{k}_{i}^{{\rm{I}}}({T}_{{\rm{e}}})-{n}_{{\rm{e}}}\,{n}_{+}{k}_{i}^{{\rm{R}}}({T}_{{\rm{e}}})]\\ & & -\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j\in { \mathcal I }}{j\ne i}}{n}_{{\rm{e}}}\,[{n}_{i}\,{k}_{{ij}}^{{\rm{E}}}({T}_{{\rm{e}}})-{n}_{j}\,{k}_{{ji}}^{{\rm{E}}}({T}_{{\rm{e}}})],\end{eqnarray} \tag{ 13 }$$

$i\in { \mathcal I }$, where the acronym "col" (i.e., collisional) has been introduced to distinguish from the mass production terms due to radiative processes treated in Section 2.3.2. The endothermic rate coefficients for collisional excitation and ionization are evaluated based on the assumption of stationary heavy particles:

$$\begin{eqnarray}{k}_{{ij}}^{{\rm{E}}}({T}_{{\rm{e}}}) & = & \sqrt{\displaystyle \frac{8{k}_{{\rm{B}}}{T}_{{\rm{e}}}}{\pi {m}_{{\rm{e}}}}}\displaystyle \frac{1}{{({k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{2}}\\ & & \times {\displaystyle \int }_{{E}_{{ij}}}^{+\infty }{\sigma }_{{ij}}^{{\rm{E}}}(\varepsilon )\,\varepsilon \,\exp \left(-\displaystyle \frac{\varepsilon }{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right)d\varepsilon ,\quad i\lt j,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{k}_{i}^{{\rm{I}}}({T}_{{\rm{e}}})=\sqrt{\displaystyle \frac{8{k}_{{\rm{B}}}{T}_{{\rm{e}}}}{\pi {m}_{{\rm{e}}}}}\displaystyle \frac{1}{{({k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{2}}{\int }_{{E}_{+i}}^{+\infty }{\sigma }_{i}^{{\rm{I}}}(\varepsilon )\,\varepsilon \,\exp \left(-\displaystyle \frac{\varepsilon }{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right)d\varepsilon ,\end{eqnarray} \tag{ 15 }$$

$i,j\in { \mathcal I }$, where the quantities ${\sigma }_{{ij}}^{{\rm{E}}}$ and ${\sigma }_{i}^{{\rm{I}}}$ denote, respectively, the total cross-section for excitation and ionization, with the corresponding threshold energies being defined as ${E}_{{ij}}={E}_{j}-{E}_{i}$ and ${E}_{+i}={E}_{+}-{E}_{i}$. The symbol ε stands for the free-electron energy. The exothermic rate coefficients for de-excitation and three-body recombination are related to those for excitation and ionization, respectively, via micro-reversibility (Oxenius 1986, Equations (2.3.10) and (2.3.16)):

$$\begin{eqnarray}\displaystyle \frac{{k}_{{ji}}^{{\rm{E}}}({T}_{{\rm{e}}})}{{k}_{{ij}}^{{\rm{E}}}({T}_{{\rm{e}}})} & = & \displaystyle \frac{{g}_{i}}{{g}_{j}}\exp \left(\displaystyle \frac{{E}_{{ij}}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right),\quad i\lt j,\\ \quad \displaystyle \frac{{k}_{i}^{{\rm{R}}}({T}_{{\rm{e}}})}{{k}_{i}^{{\rm{I}}}({T}_{{\rm{e}}})} & = & \displaystyle \frac{{g}_{i}\,{Z}_{{\rm{H}}}^{{\rm{t}}}({T}_{{\rm{e}}})}{{g}_{{\rm{e}}}\,{Z}_{{\rm{e}}}^{{\rm{t}}}({T}_{{\rm{e}}})\,{g}_{+}\,{Z}_{{{\rm{H}}}^{+}}^{{\rm{t}}}({T}_{{\rm{e}}})}\exp \left(\displaystyle \frac{{E}_{+i}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right),\end{eqnarray} \tag{ 16 }$$

$i,j\in { \mathcal I }$. It is worth mentioning that the free-electron temperature, ${T}_{{\rm{e}}}$, is the only temperature appearing in the rate coefficients (14)–(16), due to the hypothesis of stationary heavy particles. In the present work, the cross-section for electron impact excitation and ionization have been modeled according to Drawin's semi-classical formula (Drawin 1963), with the absorption oscillator strengths taken from the NIST atomic database (Kramida 2010; Kramida et al. 2014).

The volumetric time rate of change of free-electron energy due to excitation, ionization, and elastic collisions with heavy particles can be written as ${{\rm{\Omega }}}_{{\rm{e}}}^{\mathrm{col}}={{\rm{\Omega }}}_{{\rm{e}}}^{{\rm{EL}}}+{{\rm{\Omega }}}_{{\rm{e}}}^{{\rm{E}}}+{{\rm{\Omega }}}_{{\rm{e}}}^{{\rm{I}}}$, where the individual contributions are obtained by taking the second-order moment (i.e., energy) of the related collision operators in the Boltzmann equation for free electrons:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{e}}}^{{\rm{EL}}}=3\,{n}_{{\rm{e}}}\,{k}_{{\rm{B}}}\,({T}_{{\rm{e}}}-{T}_{{\rm{h}}})\,{\nu }_{{\rm{e}} \mbox{-} {\rm{h}}}^{{\rm{EL}}},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}{{\rm{\Omega }}}_{{\rm{e}}}^{{\rm{E}}} & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{i,j\in { \mathcal I }}{j\gt i}}{n}_{{\rm{e}}}\,{E}_{{ji}}\,[{n}_{i}\,{k}_{{ij}}^{{\rm{E}}}({T}_{{\rm{e}}})-{n}_{j}\,{k}_{{ji}}^{{\rm{E}}}({T}_{{\rm{e}}})],\\ \quad {{\rm{\Omega }}}_{{\rm{e}}}^{{\rm{I}}} & = & \displaystyle \sum _{i\in { \mathcal I }}\,{n}_{{\rm{e}}}\,{E}_{+i}\,[{n}_{i}\,{k}_{i}^{{\rm{I}}}({T}_{{\rm{e}}})-{n}_{{\rm{e}}}\,{n}_{+}\,{k}_{i}^{{\rm{R}}}({T}_{{\rm{e}}})],\end{eqnarray} \tag{ 18 }$$

where the effective collision frequency for elastic energy transfer in electron-heavy interactions reads (Petschek & Byron 1957; Desloge 1962; Mitchner & Kruger 1973):

$$\begin{eqnarray}&&{\nu }_{{\rm{e}} \mbox{-} {\rm{h}}}^{{\rm{EL}}}=\displaystyle \frac{8}{3}\displaystyle \sum _{c\in {{ \mathcal C }}_{{\rm{h}}}}\left(\displaystyle \frac{{m}_{{\rm{e}}}}{{m}_{c}}\right){n}_{c}\,{{\rm{\Omega }}}_{{\rm{e}}c}^{(1,1)}.\end{eqnarray} \tag{ 19 }$$

The quantities ${{\rm{\Omega }}}_{{\rm{e}}c}^{(1,1)}$ denote the electron-heavy collision integrals in the first-order Laguerre–Sonine polynomial expansion (Devoto 1966, 1967; Ferziger & Kaper 1972). The collision integrals ${{\rm{\Omega }}}_{{\rm{e}}c}^{(1,1)}$ are often referred to as diffusion cross-sections as they appear in the electron-heavy binary diffusion coefficients for ionized gas. For interactions between free electrons and hydrogen, the collision integral ${{\rm{\Omega }}}_{{\rm{e}}{\rm{H}}}^{(1,1)}$ has been evaluated based on the work by Bruno et al. (2010). For interactions between free electrons and protons, a screened Coulomb potential has been used to evaluate ${{\rm{\Omega }}}_{{\rm{e}}{{\rm{H}}}^{+}}^{(1,1)}$ (Devoto 1966, 1967; Bruno et al. 2010).

The functional form of the excitation and ionization energy transfer terms in Equation (18) reflect the physical fact that, under the assumption of stationary heavy particles, the only active source/sink of free-electron energy in inelastic and ionizing collisions is the change of heavy-particle formation energies.

2.3.2. Radiative Processes

The modeling of radiation in plasmas requires accounting for transitions characterized by absorption, emission, and scattering of light. Since the main focus of this work is the complexity reduction of the NLTE kinetics, scattering is not taken into account.

The radiative processes leading to emission and absorption of light can be subdivided into three groups: bound–bound (bb), bound–free/free–bound (bf/fb), and free–free (ff) (Oxenius 1986, Chapter 1; Hubeny & Mihalas 2014, Chapter 5). The radiative processes considered in this work are

• (i)  
  spontaneous emission/absorption and induced emission (bb):

  $$\begin{eqnarray}&&\quad {\rm{H}}(j)\leftrightarrows {\rm{H}}(i)+{h}_{{\rm{P}}}{\nu }_{{ij}},\quad {\rm{H}}(j)+{h}_{{\rm{P}}}{\nu }_{{ij}}\to {\rm{H}}(i)+2{h}_{{\rm{P}}}{\nu }_{{ij}},\\ &&\qquad \quad i\lt j,i,j\in { \mathcal I },\end{eqnarray} \tag{ 20 }$$
• (ii)  
  photo-ionization and spontaneous/induced radiative recombination (bf/fb):

  $$\begin{eqnarray}&&{\rm{H}}(i)+{h}_{{\rm{P}}}\nu \leftrightarrows {{\rm{H}}}^{+}+{e}^{-},\quad i\in { \mathcal I },\end{eqnarray} \tag{ 21 }$$
• (iii)  
  spontaneous/induced and inverse bremsstrahlung for protons (ff).

bb transitions—The monochromatic emission and absorption coefficients due to bb radiation are (Oxenius 1986, Equations (3.2.9) and (3.2.15)):

$$\begin{eqnarray}{\varepsilon }_{\lambda }^{{\rm{BB}}} & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{i,j\in { \mathcal I }}{j\gt i}}\displaystyle \frac{{h}_{{\rm{P}}}c}{4\pi {\lambda }_{{ij}}}\,{n}_{j}\,{A}_{{ji}}\,{\phi }_{\lambda }^{{ji}},\\ \quad {\kappa }_{\lambda }^{{\rm{BB}}} & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{i,j\in { \mathcal I }}{j\gt i}}\displaystyle \frac{{h}_{{\rm{P}}}{\lambda }_{{ij}}}{4\pi }\,[{n}_{i}\,{B}_{{ij}}\,{\psi }_{\lambda }^{{ij}}-{n}_{j}\,{B}_{{ji}}\,{\phi }_{\lambda }^{{ji}}],\end{eqnarray} \tag{ 22 }$$

where the wavelength associated with the transition $j\to i$ is ${\lambda }_{{ij}}={h}_{{\rm{P}}}c/{E}_{{ji}}$, with c being the speed of light and the energy difference ${E}_{{ji}}$ having the same definition as in Equations (14) and (16) (i.e., ${E}_{{ji}}={E}_{j}-{E}_{i}$). The symbols ${A}_{{ji}}$, ${B}_{{ji}}$, and ${B}_{{ij}}$ (with $i\lt j$) denote, respectively, the Einstein coefficients for spontaneous emission, induced emission, and absorption. The Einstein coefficients are not independent of each other and satisfy the Einstein–Milne relations (Oxenius 1986, Equations (1.6.31) and (1.6.32)):

$$\begin{eqnarray}&&{B}_{{ji}}=\displaystyle \frac{{\lambda }_{{ij}}^{3}}{2{h}_{{\rm{P}}}c}{A}_{{ji}},\quad {g}_{j}{B}_{{ji}}={g}_{i}{B}_{{ij}},\quad i\lt j,\end{eqnarray} \tag{ 23 }$$

$i,j\in { \mathcal I }$. The functions ${\phi }_{\lambda }^{{ji}}$ and ${\psi }_{\lambda }^{{ij}}$ in Equation (22) are, respectively, the monochromatic line emission and absorption profiles. In this work, it has been assumed that the emission and absorption profiles coincide (i.e., complete redistribution) and are described using a Voigt function by accounting for natural, Doppler, and collisional broadening (Hubeny & Mihalas 2014, Chapter 8). Doppler broadening has been evaluated at the heavy-particle temperature ${T}_{{\rm{h}}}$ as explained in Hubeny & Mihalas (2014), Equation (8.18). Collisional broadening accounts for collisions between charged particles (i.e., Stark broadening) and has been computed at the free-electron temperature ${T}_{{\rm{e}}}$ based on the work of Cowley (1971). The effects of resonance and pressure broadening and Doppler shifts due to macroscopic gas motion have not been taken into account. In order to speed up the calculations, the numerical evaluation of the Voigt function (i.e., convolution between Gaussian and Lorentzian profiles) has been accomplished using the curve fit by Whiting (1968).

Before proceeding further, it is worth spending some words on the definition of the emission and absorption coefficients given in Equation (22). In the aforementioned expression, induced/stimulated emission is treated as a negative absorption. However, according to the principles of electrodynamics, the true monochromatic emission coefficient should account for both the spontaneous and the induced/stimulated contributions (Oxenius 1986, Chapter 3). Hence, the expressions given in Equation (22) should be interpreted as the phenomenological emission and absorption coefficients. The treatment of induced/stimulated emission as a negative absorption is adopted only with the purpose of isolating the plasma optical properties, which do not explicitly depend on the local radiation field (i.e., monochromatic intensity). The convenience of this approach becomes clear when writing down the RTE (see Section 2.4). In the present work, the same convention is also used in the definition of the plasma optical properties due to bf/fb and ff radiation.

bb transitions lead to a change in the occupation numbers of hydrogen bound states, for which the related mass production terms are (Oxenius 1986, Equation (2.4.13)):

$$\begin{eqnarray}{\omega }_{i}^{{\rm{BB}}} & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j\in { \mathcal I }}{j\gt i}}\,[{n}_{j}\,{A}_{{ji}}^{{\rm{m}}}-({n}_{i}\,{B}_{{ij}}^{{\rm{m}}}-{n}_{j}\,{B}_{{ji}}^{{\rm{m}}})\,{{\rm{\Phi }}}_{{ij}}]\\ & & -\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{j\in { \mathcal I }}{j\lt i}}\,[{n}_{i}\,{A}_{{ij}}^{{\rm{m}}}-({n}_{j}\,{B}_{{ji}}^{{\rm{m}}}-{n}_{i}\,{B}_{{ij}}^{{\rm{m}}})\,{{\rm{\Phi }}}_{{ij}}],\end{eqnarray} \tag{ 24 }$$

$i\in { \mathcal I }$, where the mass production (superscript m) Einstein coefficients, defined by ${A}_{{ji}}^{{\rm{m}}}={A}_{{ji}}$ (with $i\lt j$) and ${B}_{{ij}}^{{\rm{m}}}={B}_{{ji}}\,{\lambda }_{{ij}}^{2}/c$, have been introduced for convenience. The quantity ${{\rm{\Phi }}}_{{ij}}$ is defined by the line-shape integral:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{ij}}={\int }_{0}^{+\infty }{J}_{\lambda }\,{\phi }_{\lambda }^{{ij}}\,d\lambda ,\quad i\ne j,\end{eqnarray} \tag{ 25 }$$

$i,j\in { \mathcal I }$, where the symbol ${J}_{\lambda }$ denotes the average monochromatic intensity. The former is obtained by the integration of the (directionally dependent) monochromatic intensity ${I}_{\lambda }$ over all directions, ${J}_{\lambda }=1/4\pi \oint {I}_{\lambda }d{\rm{\Omega }}$, with Ω being the solid angle (Hubeny & Mihalas 2014, Chapter 11). The changes in the occupation numbers of hydrogen bound states are accompanied by a transfer of energy between the material gas and the radiation field. The volumetric time rate of loss of matter energy due to bb transitions is obtained by integrating the quantity ${\varepsilon }_{\lambda }^{{\rm{BB}}}-{I}_{\lambda }\,{\kappa }_{\lambda }^{{\rm{BB}}}$ over all wavelengths and directions (Oxenius 1986, Chapter 3). Performing the required integrations, one obtains

$$\begin{eqnarray}{{\rm{\Omega }}}^{{\rm{BB}}} & = & {\displaystyle \int }_{0}^{+\infty }\oint ({\varepsilon }_{\lambda }^{{\rm{BB}}}-{I}_{\lambda }\,{\kappa }_{\lambda }^{{\rm{BB}}})\,d\lambda \,d{\rm{\Omega }}\\ & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{i,j\in { \mathcal I }}{j\gt i}}\,[{n}_{j}\,{A}_{{ji}}^{{\rm{e}}}-({n}_{i}\,{B}_{{ij}}^{{\rm{e}}}-{n}_{j}\,{B}_{{ji}}^{{\rm{e}}})\,{{\rm{\Phi }}}_{{ji}}],\end{eqnarray} \tag{ 26 }$$

where the energy transfer (e) Einstein coefficients are defined as ${A}_{{ji}}^{{\rm{e}}}={E}_{{ji}}\,{A}_{{ji}}$ (with $i\lt j$) and ${B}_{{ji}}^{{\rm{e}}}={h}_{{\rm{P}}}{\lambda }_{{ij}}{B}_{{ji}}$.

BF/FB transitions—The monochromatic emission and absorption coefficients due to bf and fb radiation are (Oxenius 1986, Equations (3.2.28) and (3.2.32)):

$$\begin{eqnarray}{\varepsilon }_{\lambda }^{{\rm{FB}}} & = & \displaystyle \frac{{h}_{{\rm{P}}}^{4}\,{c}^{2}\,{n}_{{\rm{e}}}\,{n}_{+}}{{\lambda }^{5}{(2\pi {m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}}\\ & & \times \displaystyle \sum _{i\in { \mathcal I }}\,{\sigma }_{i}^{{\rm{PI}}}(\lambda )\displaystyle \frac{{g}_{i}}{{g}_{+}}\exp \left(\displaystyle \frac{{E}_{+i}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right),\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}{\kappa }_{\lambda }^{{\rm{BF}}} & = & \displaystyle \sum _{i\in { \mathcal I }}\,{\sigma }_{i}^{{\rm{PI}}}(\lambda )\left[{n}_{i}-\displaystyle \frac{1}{2}\displaystyle \frac{{h}_{{\rm{P}}}^{3}\,{n}_{{\rm{e}}}\,{n}_{+}}{{(2\pi {m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}}\displaystyle \frac{{g}_{i}}{{g}_{+}}\right.\\ & & \times \left.\exp \left(\displaystyle \frac{{E}_{+i}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)\right],\end{eqnarray} \tag{ 28 }$$

where the ionization threshold is defined, as in Equations (15) and (16), via the relation ${E}_{+i}={E}_{+}-{E}_{i}$. The quantity ${\sigma }^{{\rm{PI}}}$ denotes the total photo-ionization cross-section. In the present work, the former has been evaluated based on Kramer's formula (Zel'dovich & Raizer 1967, Equation (5.34); Hubeny & Mihalas 2014, Equation (7.84)).

The mass production terms for free electrons, protons, and hydrogen bound states due to photo-ionization and radiative recombination can be expressed as (Oxenius 1986, Equations (2.4.22) and (2.4.27)):

$$\begin{eqnarray}&&{\omega }_{{\rm{e}}}^{{\rm{BF}}/{\rm{FB}}}=\displaystyle \sum _{i\in { \mathcal I }}[{n}_{i}\,{k}_{i}^{{\rm{PI}}}({J}_{\lambda })-{n}_{{\rm{e}}}\,{n}_{+}\,{k}_{i}^{{\rm{RR}}}({T}_{{\rm{e}}},{J}_{\lambda })],\\ &&\quad {\omega }_{+}^{{\rm{BF}}/{\rm{FB}}}={\omega }_{{\rm{e}}}^{{\rm{BF}}/{\rm{FB}}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\omega }_{i}^{{\rm{BF}}/{\rm{FB}}}=-[{n}_{i}\,{k}_{i}^{{\rm{PI}}}({J}_{\lambda })-{n}_{{\rm{e}}}\,{n}_{+}\,{k}_{i}^{{\rm{RR}}}({T}_{{\rm{e}}},{J}_{\lambda })],\end{eqnarray} \tag{ 30 }$$

$i\in { \mathcal I }$, where the rate coefficients for photo-ionization and radiative recombination (accounting for both spontaneous (s) and induced (i) contributions) are

$$\begin{eqnarray}{k}_{i}^{{\rm{PI}}}({J}_{\lambda }) & = & \displaystyle \frac{4\pi }{{h}_{{\rm{P}}}c}{\displaystyle \int }_{0}^{{\lambda }_{+i}}{\sigma }_{i}^{{\rm{PI}}}(\lambda )\,{J}_{\lambda }\,\lambda \,d\lambda ,\\ {k}_{i}^{{\rm{RR}}}({T}_{{\rm{e}}},{J}_{\lambda }) & = & {k}_{i}^{{\rm{RR}} \mbox{-} {\rm{s}}}({T}_{{\rm{e}}})+{k}_{i}^{{\rm{RR}} \mbox{-} {\rm{i}}}({T}_{{\rm{e}}},{J}_{\lambda }),\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}{k}_{i}^{{\rm{RR}} \mbox{-} {\rm{s}}}({T}_{{\rm{e}}}) & = & \sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{{g}_{i}}{{g}_{+}}\displaystyle \frac{{h}_{{\rm{P}}}^{3}c}{{({m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}}\exp \left(\displaystyle \frac{{E}_{+i}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right)\\ & & \times {\displaystyle \int }_{0}^{{\lambda }_{+i}}\displaystyle \frac{{\sigma }_{i}^{{\rm{PI}}}(\lambda )}{{\lambda }^{4}}\,\exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)d\lambda ,\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}{k}_{i}^{{\rm{RR}} \mbox{-} {\rm{i}}}({T}_{{\rm{e}}},{J}_{\lambda }) & = & \sqrt{\displaystyle \frac{1}{2\pi }}\displaystyle \frac{{g}_{i}}{{g}_{+}}\displaystyle \frac{{h}_{{\rm{P}}}^{2}}{{({m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}c}\exp \left(\displaystyle \frac{{E}_{+i}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right)\\ & & \times {\displaystyle \int }_{0}^{{\lambda }_{+i}}{\sigma }_{i}^{{\rm{PI}}}(\lambda )\,\lambda \,{J}_{\lambda }\,\exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)d\lambda ,\end{eqnarray} \tag{ 33 }$$

$i\in { \mathcal I }$, with the threshold wavelength for photo-ionization/radiative recombination being ${\lambda }_{+i}={h}_{{\rm{P}}}c/{E}_{+i}$.

By analogy with the procedure outlined above for bb radiation, the volumetric time rate of loss of matter energy due to bf/fb radiation is obtained through the integration of the quantity ${\varepsilon }_{\lambda }^{{\rm{FB}}}-{I}_{\lambda }\,{\kappa }_{\lambda }^{{\rm{BF}}}$ over all wavelengths and directions:

$$\begin{eqnarray}{{\rm{\Omega }}}^{{\rm{BF}}/{\rm{FB}}} & = & \displaystyle \sum _{i\in { \mathcal I }}\left[\sqrt{\displaystyle \frac{1}{2\pi }}\displaystyle \frac{{g}_{i}}{{g}_{+}}\displaystyle \frac{{h}_{{\rm{P}}}^{3}\,{n}_{{\rm{e}}}\,{n}_{+}}{{({m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}}\right.\\ & & \times \exp \left(\displaystyle \frac{{E}_{+i}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right){\displaystyle \int }_{0}^{{\lambda }_{+i}}{\sigma }_{i}^{{\rm{PI}}}(\lambda )\,{J}_{\lambda }\\ & & \times \exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)d\lambda \\ & & -\left.4\pi {n}_{i}{\displaystyle \int }_{0}^{{\lambda }_{+i}}{\sigma }_{i}^{{\rm{PI}}}(\lambda )\,{J}_{\lambda }\,d\lambda \right]\\ & & +\displaystyle \sum _{i\in { \mathcal I }}\sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{{g}_{i}}{{g}_{+}}\displaystyle \frac{{h}_{{\rm{P}}}^{4}{c}^{2}\,{n}_{{\rm{e}}}\,{n}_{+}}{{({m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}}\\ & & \times \exp \left(\displaystyle \frac{{E}_{+i}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right){\displaystyle \int }_{0}^{{\lambda }_{+i}}\displaystyle \frac{{\sigma }_{i}^{{\rm{PI}}}(\lambda )}{{\lambda }^{5}}\\ & & \times \exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)d\lambda .\end{eqnarray} \tag{ 34 }$$

It should be noted that the volumetric energy loss term (34) refers to the whole gas, which includes free electrons and heavy particles. The adoption of an additional energy equation for free electrons (see Section 2.4) requires the evaluation of the corresponding quantity ${{\rm{\Omega }}}^{{\rm{BF}}/{\rm{FB}}}$ for the free-electron gas alone. This is accomplished in a straightforward manner by taking the second-order moment of the collision operator for photo-ionization/radiative recombination of the Boltzmann equation for free electrons. Under the assumption of Maxwellian free electrons at temperature ${T}_{{\rm{e}}}$ and stationary heavy particles, the final result is

$$\begin{eqnarray}{{\rm{\Omega }}}_{{\rm{e}}}^{{\rm{BF}}/{\rm{FB}}} & = & {{\rm{\Omega }}}^{{\rm{BF}}/{\rm{FB}}}\\ & & +\displaystyle \sum _{i\in { \mathcal I }}{E}_{+i}\,[{n}_{i}\,{k}_{i}^{{\rm{PI}}}({J}_{\lambda })-{n}_{{\rm{e}}}\,{n}_{+}\,{k}_{i}^{{\rm{RR}}}({T}_{{\rm{e}}},{J}_{\lambda })].\end{eqnarray} \tag{ 35 }$$

ff transitions—The monochromatic emission and absorption coefficients due to ff radiation produced by encounters between free electrons and protons are (Oxenius 1986, Equations (3.2.41) and (3.2.42)):

$$\begin{eqnarray}&&{\varepsilon }_{\lambda }^{{\rm{FF}}}=\displaystyle \frac{8}{3}\sqrt{\displaystyle \frac{2\pi }{3{m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}}}}\displaystyle \frac{{n}_{{\rm{e}}}\,{n}_{+}\,{q}_{{\rm{e}}}^{6}}{{\lambda }^{2}\,{(4\pi {\epsilon }_{0})}^{3}{m}_{{\rm{e}}}{c}^{2}}\exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right),\\ &&{\kappa }_{\lambda }^{{\rm{FF}}}=\displaystyle \frac{4}{3}\sqrt{\displaystyle \frac{2\pi }{3{m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}}}}\displaystyle \frac{{\lambda }^{3}\,{n}_{{\rm{e}}}\,{n}_{+}\,{q}_{{\rm{e}}}^{6}}{{(4\pi {\epsilon }_{0})}^{3}{m}_{{\rm{e}}}{h}_{{\rm{P}}}{c}^{4}}\left[1-\exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)\right],\end{eqnarray} \tag{ 36 }$$

where the quantities ${\epsilon }_{0}$ and ${q}_{{\rm{e}}}$ denote, respectively, the vacuum permittivity and the electron charge.

The volumetric time rate of loss of matter energy due to ff radiation is obtained, as done before, by integrating quantity ${\varepsilon }_{\lambda }^{{\rm{FF}}}-{I}_{\lambda }\,{\kappa }_{\lambda }^{{\rm{FF}}}$ over all wavelengths and directions:

$$\begin{eqnarray}{{\rm{\Omega }}}^{{\rm{FF}}} & = & \displaystyle \frac{32}{3}\pi \sqrt{\displaystyle \frac{2\pi }{3{m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}}}}\displaystyle \frac{{n}_{{\rm{e}}}\,{n}_{+}\,{q}_{{\rm{e}}}^{6}}{{(4\pi {\epsilon }_{0})}^{3}{m}_{{\rm{e}}}{h}_{{\rm{P}}}}\\ & & \times \left\{{k}_{{\rm{B}}}{T}_{{\rm{e}}}-\displaystyle \frac{1}{2c}{\displaystyle \int }_{0}^{+\infty }{\lambda }^{3}\left[1-\exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)\right]\right.\\ & & \times {J}_{\lambda }\,d\lambda \}.\end{eqnarray} \tag{ 37 }$$

By collecting the results through Equations (22)–(37), it is possible to write down the mass production terms due to radiative transitions for free electrons, protons, and hydrogen bound states as ${\omega }_{{\rm{e}}}^{\mathrm{rad}}={\omega }_{{\rm{e}}}^{{\rm{BF}}/{\rm{FB}}}$, ${\omega }_{+}^{\mathrm{rad}}={\omega }_{+}^{{\rm{BF}}/{\rm{FB}}}$, and ${\omega }_{i}^{\mathrm{rad}}={\omega }_{i}^{{\rm{BB}}}+{\omega }_{i}^{{\rm{BF}}/{\rm{FB}}}$, respectively. The corresponding energy loss rates for the whole gas and free electrons alone are ${{\rm{\Omega }}}^{\mathrm{rad}}={{\rm{\Omega }}}^{{\rm{BB}}}+{{\rm{\Omega }}}^{{\rm{BF}}/{\rm{FB}}}+{{\rm{\Omega }}}^{{\rm{FF}}}$ and ${{\rm{\Omega }}}_{{\rm{e}}}^{\mathrm{rad}}={{\rm{\Omega }}}_{{\rm{e}}}^{{\rm{BF}}/{\rm{FB}}}+{{\rm{\Omega }}}^{{\rm{FF}}}$, respectively.

The steady flow across a normal shock wave of a two-temperature radiating plasma is governed by the species continuity equations, the global momentum and energy equations, and the free-electron energy equation. In the absence of transport phenomena, the former set of equations reads (Zel'dovich & Raizer 1967, Equation (7.30); Murty 1971, Equation (1); see also Appendix C):

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial x}\left(\begin{array}{c}{n}_{s}u\\ p+\rho {u}^{2}\\ \rho {Hu}\\ \rho {e}_{{\rm{e}}}u\end{array}\right)=\left(\begin{array}{c}{\omega }_{s}\\ 0\\ -\,{{\rm{\Omega }}}^{\mathrm{rad}}\\ -{p}_{{\rm{e}}}\displaystyle \frac{\partial u}{\partial x}-({{\rm{\Omega }}}_{{\rm{e}}}^{\mathrm{col}}+{{\rm{\Omega }}}_{{\rm{e}}}^{\mathrm{rad}})\end{array}\right),\end{eqnarray} \tag{ 38 }$$

$s\in { \mathcal S }$, where the mass production terms are given by the sum of the collisional and radiative contributions, ${\omega }_{s}={\omega }_{s}^{\mathrm{col}}+{\omega }_{s}^{\mathrm{rad}}$. The quantity u stands for the flow velocity measured in the shock wave reference frame. Gravity terms have not been taken into account as the length scales of the radiative shock waves investigated in this work (e.g., 10¹–10² m) are much smaller than those typical of stellar atmospheres (e.g., 10⁶ m).

It is known that electron heat conduction, neglected in Equation (38), may play an important role in shaping the temperature profile in the precursor (i.e., conduction precursor) (Jukes 1957; Shafranov 1957; Imshennik 1962; Jaffrin & Probstein 1964; Fadeyev et al. 2002). However, for the standing shocks studied in this work, electron heat conduction is influenced little by the dynamics of excited electronic states,⁸ and thus is expected to play a minor role compared to radiation and chemistry on the accuracy of a reduced-order NLTE model.

It is worth noticing that, by assuming zero velocity and adding gravity terms in Equation (38), one may retrieve the structural equations for the material gas within a one-dimensional plane-parallel (and static) stellar atmosphere. As a matter of fact, under the above assumption, the species continuity equations reduce to ${\omega }_{s}=0$, which is precisely the condition of statistical equilibrium, whereas the global momentum and energy equations reduce to the requirements of hydrostatic and radiative equilibrium, respectively (i.e., $\partial p/\partial x=-\rho g$ and ${{\rm{\Omega }}}^{\mathrm{rad}}=0$).

For a radiating gas, the flow-governing Equations (38) must be coupled with the RTE (Mihalas & Mihalas 1984, Chapters 7 and 8). The former can be thought of as the kinetic equation for a photon gas and describes the evolution in space and time of a radiation field, due to emission, absorption, and scattering of light. In the case of a plane-parallel non-scattering medium under steady-state conditions, the RTE reads

$$\begin{eqnarray}&&\mu \displaystyle \frac{\partial {I}_{\lambda \mu }}{\partial x}={\varepsilon }_{\lambda }^{}-{\kappa }_{\lambda }^{}{I}_{\lambda \mu },\end{eqnarray} \tag{ 39 }$$

where the quantity $\mu \in [-1,1]$ stands for the cosine of the angle between the line of sight and the x-axis. The (total) emission and absorption coefficients in the RTE (39) are obtained by summing the individual contributions due to the bb, bf/fb, and ff transitions as ${\varepsilon }_{\lambda }^{}={\varepsilon }_{\lambda }^{{\rm{BB}}}+{\varepsilon }_{\lambda }^{{\rm{FB}}}+{\varepsilon }_{\lambda }^{{\rm{FF}}}$ and ${\kappa }_{\lambda }^{}={\kappa }_{\lambda }^{{\rm{BB}}}+{\kappa }_{\lambda }^{{\rm{BF}}}+{\kappa }_{\lambda }^{{\rm{FF}}}$, respectively. In radiative transfer problems within plane-parallel media, it is often convenient to transform the RTE (39) to a second-order differential equation as proposed by Feautrier (1964):

$$\begin{eqnarray}&&{\mu }^{2}\displaystyle \frac{{\partial }^{2}{P}_{\lambda \mu }}{\partial {\tau }_{\lambda }^{2}}={P}_{\lambda \mu }-{S}_{\lambda },\end{eqnarray} \tag{ 40 }$$

where the quantity ${P}_{\lambda \mu }$ is defined as ${P}_{\lambda \mu }=({I}_{\lambda \mu }+{I}_{\lambda -\mu })/2$, and the monochromatic source function is given by the ratio of emission and absorption coefficients, ${S}_{\lambda }={\varepsilon }_{\lambda }^{}/{\kappa }_{\lambda }^{}$. The symbol $d{\tau }_{\lambda }$ denotes the infinitesimal monochromatic optical thickness increment, $d{\tau }_{\lambda }={\kappa }_{\lambda }^{}\,{dx}$. For the RTE (40), the range of the angular variable μ is restricted to $[0,1]$ (Hubeny & Mihalas 2014, Chapter 11). When written in terms of the newly introduced unknown ${P}_{\lambda \mu }$, the average monochromatic intensity becomes simply

$$\begin{eqnarray}&&{J}_{\lambda }=\displaystyle \frac{1}{4\pi }\oint {I}_{\lambda }\,d{\rm{\Omega }}={\int }_{0}^{1}{P}_{\lambda \mu }\,d\mu .\end{eqnarray} \tag{ 41 }$$

Following the work by Liu et al. (2015), the logarithm of the normalized population within a given group k is written as a polynomial in the internal energy ${E}_{i}$:

$$\begin{eqnarray}&&\mathrm{ln}\left(\displaystyle \frac{{n}_{i}}{{g}_{i}}\right)={\tilde{\alpha }}_{k}+{\tilde{\beta }}_{k}{E}_{i}+{\tilde{\gamma }}_{k}{E}_{i}^{2}+{\rm{H}}.{\rm{O}}.{\rm{T}},\end{eqnarray} \tag{ 42 }$$

$i\in {{ \mathcal I }}_{k}$, $k\in { \mathcal K }$, where the sets ${{ \mathcal I }}_{k}$ and ${ \mathcal K }$ denote, respectively, the energy levels within group k and the group indices. In the present work, terms that are second and higher order in energy are neglected. Under these circumstances, one needs to determine only the quantities ${\tilde{\alpha }}_{k}$ and ${\tilde{\beta }}_{k}$. These are related to the group populations, ${\tilde{n}}_{k}$, and average energies, ${\tilde{E}}_{k}$, by the following moment constraints on particle number and energy (Liu et al. 2015):

$$\begin{eqnarray}&&{\tilde{n}}_{k}=\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{n}_{i},\quad {\tilde{n}}_{k}{\tilde{E}}_{k}=\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{n}_{i}{E}_{i},\end{eqnarray} \tag{ 43 }$$

$k\in { \mathcal K }$. Substituting Equation (42) into the moment constraints (43) gives the implicit relation linking ${\tilde{\alpha }}_{k}$ and ${\tilde{\beta }}_{k}$ to the group populations and energies:

$$\begin{eqnarray}&&{\tilde{n}}_{k}=\exp {\tilde{\alpha }}_{k}\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{g}_{i}\,\exp ({\tilde{\beta }}_{k}{E}_{i}),\quad {\tilde{E}}_{k}=\displaystyle \frac{\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{g}_{i}\,{E}_{i}\,\exp ({\tilde{\beta }}_{k}{E}_{i})}{\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{g}_{i}\exp ({\tilde{\beta }}_{k}{E}_{i})},\end{eqnarray} \tag{ 44 }$$

$k\in { \mathcal K }$. After introducing, for the sake of convenience, group temperatures ${T}_{k}=-1/{k}_{{\rm{B}}}{\tilde{\beta }}_{k}$ and partition functions ${\tilde{Z}}_{k}^{}({T}_{k})\ ={\sum }_{i\in {{ \mathcal I }}_{k}}{g}_{i}\exp ({\tilde{\beta }}_{k}{E}_{i})$, it is possible to re-write Equation (44) as

$$\begin{eqnarray}{\tilde{\alpha }}_{k} & = & \mathrm{ln}{\tilde{n}}_{k}-\mathrm{ln}[{\tilde{Z}}_{k}({T}_{{\rm{k}}})],\\ {\tilde{E}}_{k} & = & \displaystyle \frac{1}{{\tilde{Z}}_{k}^{}({T}_{k})}\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{g}_{i}\,{E}_{i}\,\exp \left(-\displaystyle \frac{{E}_{i}}{{k}_{{\rm{B}}}{T}_{k}}\right),\end{eqnarray} \tag{ 45 }$$

$k\in { \mathcal K }$. Substituting the first of Equation (45) into Equation (42) leads finally to

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{i}}{{g}_{i}}=\displaystyle \frac{{\tilde{n}}_{k}}{{\tilde{Z}}_{k}^{}({T}_{k})}\exp \left(-\displaystyle \frac{{E}_{i}}{{k}_{{\rm{B}}}{T}_{k}}\right),\end{eqnarray} \tag{ 46 }$$

$i\in {{ \mathcal I }}_{k}$, $k\in { \mathcal K }$. Equation (46) shows that retaining terms up to first order in energy in Equation (43) is equivalent to assuming a Maxwell–Boltzmann distribution within each group. The Maxwell–Boltzmann distribution is the LTE distribution for which the entropy is maximum (Groot & Mazur 2011, Chapter 1). This is the reason that reduced-order models developed based on Equation (42) are called ME models (Liu et al. 2015). The model corresponding to Equation (46) is the Maximum Entropy Linear (MEL) model as only linear terms are retained in the energy polynomial (40). The MEL model reduces to the Maximum Entropy Uniform (MEU) model when taking the limit of infinite group temperatures (i.e., ${\tilde{\beta }}_{k}=0$) or, equivalently, when retaining only the zeroth-order energy term in Equation (42). In this case, Equation (46) reduces to

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{i}}{{g}_{i}}=\displaystyle \frac{{\tilde{n}}_{k}}{{\tilde{Z}}_{k}},\end{eqnarray} \tag{ 47 }$$

$i\in {{ \mathcal I }}_{k}$, $k\in { \mathcal K }$, where the group partition function is now the sum of the statistical weights of the energy levels within the group, ${\tilde{Z}}_{k}={\sum }_{i\in {{ \mathcal I }}_{k}}{g}_{i}$. It is worth mentioning that the MEU distribution (47) does not allow equilibrium (i.e., Maxwell–Boltzmann distribution) to be retrieved.

The governing equations for the ME model are obtained by taking the moments with respect to the energy ${E}_{i}$ of the species continuity equations (Liu et al. 2015). As explained in Section (3.1), in the present work only the zeroth and first-order moments are needed:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial x}\left(\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{n}_{i}u\right)=\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{\omega }_{i},\quad \displaystyle \frac{\partial }{\partial x}\left(\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{n}_{i}{E}_{i}u\right)=\displaystyle \sum _{i\in {{ \mathcal I }}_{k}}{\omega }_{i}{E}_{i},\end{eqnarray} \tag{ 48 }$$

$k\in { \mathcal K }$. Using the moment constraints (43), Equation (48) becomes

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial x}\left(\begin{array}{c}{\tilde{n}}_{k}u\\ {\tilde{n}}_{k}{\tilde{E}}_{k}u\end{array}\right)=\left(\begin{array}{c}{\tilde{\omega }}_{k}\\ {\tilde{{\rm{\Omega }}}}_{k}\end{array}\right),\end{eqnarray} \tag{ 49 }$$

$k\in { \mathcal K }$, where the group mass production and energy transfer terms are defined as ${\tilde{\omega }}_{k}={\sum }_{i\in {{ \mathcal I }}_{k}}{\omega }_{i}$ and ${\tilde{{\rm{\Omega }}}}_{k}={\sum }_{i\in {{ \mathcal I }}_{k}}{\omega }_{i}{E}_{i}$, respectively. In the case of the MEU model, the group number densities ${\tilde{n}}_{k}$ are the only unknowns. Thus, only the first of Equation (49) is needed. The complete set of flow-governing equations for the MEU/MEL models is obtained based on those for the StS model (38) by replacing the species continuity equations for the hydrogen bound states with the moment equations (49). The related expressions for thermodynamic properties, mass/energy production terms, and emission/absorption coefficients are given in Appendix B.

The MEU and MEL models have already been successfully applied to the study of collisional excitation, dissociation, and ionization in atomic and molecular gases (Panesi & Lani 2013; Munafò et al. 2014, 2015; Liu et al. 2015). In that work, the determination of the states contained in a given group can be based, for instance, on an even subdivision of the internal energy ladder. This is justified by the fact that the rate coefficients for inelastic collisional processes (e.g., electron impact ionization) are larger for states with similar energy. The inclusion of radiative transitions (in particular bb radiation) completely changes the picture for states that are close in energy but are strongly coupled through radiative transitions. Such states would be inaccurately modeled by being placed within the same energy group.

For convenience, the shock, which is treated as a discontinuous surface, is placed at x = 0. The left and right boundaries are placed at $x=-{x}_{{\rm{L}}}$ and $x={x}_{{\rm{R}}}$, respectively. The lengths ${x}_{{\rm{L}}}$ and ${x}_{{\rm{R}}}$ are set to values of the order of 10²–10³ m in order to include the whole extent of the precursor and radiative relaxation regions, respectively. In order to properly resolve the smaller scales of the internal relaxation region (≃10⁻²–10¹ m), an exponential stretching is applied to reduce the grid size around the shock location.

The wavelength domain is discretized as suggested by Fadeyev & Gillet (1998, 2000). When accounting only for continuum radiation, the procedure goes as follows. After prescribing the minimum and maximum wavelengths (${\lambda }_{\min }$ and ${\lambda }_{\max }$, respectively), the interval $[{\lambda }_{\min },{\lambda }_{\max }]$ is divided into sub-intervals determined by the photo-ionization thresholds ${\lambda }_{+i}$. Each of these sub-intervals is then discretized using Gauss–Legendre quadrature points. This procedure is slightly modified to account for line radiation by adding additional sub-intervals for each atomic line. This is motivated by the rapid variation of emission and absorption coefficients (and source functions as well) over the lines. For this reason, the wavelength domain close to an atomic line l is discretized by adopting a Gauss–Legendre or uniform grid over the interval $[{\lambda }_{l}-\delta {\lambda }_{l},{\lambda }_{l}+\delta {\lambda }_{l}]$, where the width $\delta {\lambda }_{l}$ is set to $5$–$10\,\mathring{\rm A}$. The remaining sub-intervals are discretized as done for continuum radiation.

The angular variable μ is also discretized by using Gauss–Legendre quadrature points.

The flow-governing equations (38) are solved using a space-marching approach in both the pre-shock and the post-shock regions. This requires the specification of initial conditions. For the pre-shock region, LTE conditions are assumed during the first iteration. However, the gas in the far precursor is not in LTE due to non-equilibrium excitation caused by absorption of resonant radiation in atomic line wings (Murty 1968; Foley & Clarke 1973). To take this into account, the occupation numbers of free electrons, protons, and hydrogen bound states are computed (starting from the second iteration) by solving the statistical equilibrium equations (Hubeny & Mihalas 2014, Chapter 9). These equations are obtained by setting to zero the convective term in the species continuity equations (i.e., ${\omega }_{s}=0$) and are solved iteratively by a Newton–Raphson procedure. Once the shock location is reached, the integration of Equation (38) is stopped and the Rankine–Hugoniot jump relations (Zel'dovich & Raizer 1967, Chapter 1) are applied to determine the kinematic and thermo-chemical state of the gas just behind the shock. This provides the initial solution for the post-shock region. The jump relations are solved under the assumption of frozen kinetics and by neglecting the effects of radiant energy fluxes (Marshak 1958). Free electrons are assumed isothermal within the shock. Alternatively, one could consider a slightly more accurate method by treating the compression of free electrons as isentropic (Zel'dovich & Raizer 1967, Chapter 7; Fadeyev & Gillet 1998). Preliminary calculations indicated, however, that this second approach does not lead to appreciable improvements compared to the first one (isothermal compression). In the present work, the flow-governing equations (38) are numerically integrated by using a fifth-order Backward Differentiation Formula method (Gear 1971, Chapter 7) implemented in the lsode library for stiff initial value problems (Radhakrishnan & Hindmarsh 1993). For convenience, the numerical integration is performed by using the mass fractions, velocity, and temperatures as solution variables (Panesi et al. 2009, 2011; Munafò et al. 2014, 2015). This is motivated by the particularly simple form assumed by the flow-governing equations (38) when rewritten in terms of these variables.

Radiative transfer is treated by using Feautrier's method (in the improved version proposed by Ribicky & Hummer 1991) because of its straightforward implementation for one-dimensional plane-parallel problems. The Feautrier form of the RTE (40) is discretized using a conventional second-order finite difference method under the assumption of no incoming radiation from both boundaries (Auer 1967; Hubeny & Mihalas 2014). For a non-scattering medium where Doppler shifts due to bulk motions are neglected, the above procedure leads to a set of uncoupled tridiagonal systems of equations (one for each discrete angle-wavelength point). These tridiagonal systems are solved using the elimination scheme proposed by Ribicky & Hummer (1991) for the sake of better numerical conditioning.

Figure 2 shows the evolution of the temperatures and electron mole fraction across the internal relaxation and the near-precursor regions for a shock propagating at 40 km s⁻¹. The absorption of Lyman continuum radiation photo-ionizes the gas in the free stream and leads to an increase in the electron concentration, which is about 1% at the shock location. The spatial extent of the near-precursor region is two orders of magnitude larger than that of the internal relaxation region. The temperature evolution within the latter shows the typical structure observed in atomic plasmas (Zel'dovich & Raizer 1967, Chapter 7), where an initial sudden increase of the free-electron temperature is followed by a region where global ionization is indirectly controlled by elastic energy transfer in electron-heavy collisions. This may be explained as follows. Right behind the gasdynamic jump, the temperature of free electrons is of the order of 1–2 eV. The ionization of hydrogen in the internal relaxation region occurs via a two-step process: (i) excitation of high-lying electronic states and (ii) subsequent ionization from these states. Direct ionization from the ground state is also possible, although this may occur only in collisions with very energetic electrons (i.e., 10–20 eV). In both situations, free electrons lose on average 10 eV to produce a proton. Thus, the formation of a proton requires an energy equal to the thermal energy of 10 electrons. It is clear that if the free-electron energy were not replenished, their temperature would drop very rapidly. This will in turn hinder ionization as the rate of this process is proportional to the Boltzmann factor $\exp (-{E}_{+i}/{k}_{{\rm{B}}}{T}_{{\rm{e}}})$. The energy lost by free electrons in excitation/ionization collisions is replenished via elastic collisions with hot heavy particles. The inefficient energy transfer in these interactions (due to the large mass disparity) limits the amount of energy recovered by free electrons, and thus indirectly controls macroscopic ionization. Within the internal relaxation region, ionization through electron–atom collisions stimulates the further ionization of the gas, leading to avalanche (or cascade) processes (see Equation (11)). The electron avalanche causes an exponential increase in the free-electron concentration and is accompanied, in the final stages, by a rapid drop of the heavy-particle temperature (see Figure 2(b)). The avalanche terminates when the amount of protons produced is such that three-body recombination balances macroscopic ionization. After thermal equilibrium is reached, the plasma enters the radiative relaxation region where the degree of ionization decreases due to radiative recombination (see Figure 2(a)).

Zoom In Zoom Out Reset image size

Figure 2 also reports the temperature and electron mole fraction when neglecting radiative processes. This is done to show the effects of radiation on the shock structure. Neglecting radiative transitions (in particular photo-ionization) leads to a larger internal relaxation region, due to the smaller concentration of free electrons at the shock location. In the present case, the length of the above zone is almost five times larger. Once ionization is completed, the gas properties (e.g., chemical composition, temperatures) maintain their post-shock LTE values.

To investigate the number of iterations needed to obtain a converged solution, the L₂ norm of the relative error on the heavy-particle translational temperature has been monitored (see Figure 3). The relative error is defined as

$$\begin{eqnarray}&&{E}_{{L}_{2}}^{n}(T)=\sqrt{\displaystyle \frac{1}{N}\displaystyle \sum _{j=1,N}| {T}_{j}^{n}-{T}_{j}^{n-1}{| }^{2}},\end{eqnarray} \tag{ 50 }$$

where the symbol Tⁿ_j denotes the temperature at the discrete node j and global iteration n, and N = 1690 stands for the number of spatial nodes. For both situations analyzed in Figure 3 (with and without bb radiation), a converged solution is essentially obtained after 40 lambda iterations. As expected, the rate of convergence is slower when accounting for line radiation. Similar trends to those shown in Figure 3 have been observed for the other cases investigated in the paper.

Zoom In Zoom Out Reset image size

Increasing the shock speed leads to more energetic photons from the radiative cooling region. This causes a larger degree of ionization in the precursor (see Figure 4(a)), which has, in turn, the effect of reducing the length of the internal relaxation region by more than one order of magnitude as revealed by the evolution of the heavy-particle translational temperature shown in Figure 4(b). To help in highlighting the above reduction in length, a logarithmic scale is used for the x-axis in Figure 4(b). For free-stream velocities larger than 60 km s⁻¹, the hydrogen plasma in the precursor is fully ionized. Under these circumstances, when the gas enters the radiative cooling region, the electron mole fraction remains almost constant at the beginning, and then decreases due to radiative recombination. The extent of the region of (almost) constant degree of ionization increases with the shock speed and is reflected in the plateau observed for the heavy-particle temperature in Figure 4(b). At large shock speeds, the ionization avalanche is followed by a fast cooling zone, which is slowed by the onset of radiative recombination. This is demonstrated in Figure 5 which superimposes the electron mole fraction and temperature evolution at 40 and 70 km s⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The pronounced temperature drop at 70 km s⁻¹ (Figure 5(b)) is due to the absence of neutral absorbers, which makes the plasma completely transparent to continuum radiation (see Equation (34)). The temperature decrease favors radiative recombination as the rate coefficient for this process is larger at low temperatures (see Equations (32) and (33)). Once the amount of ground-state hydrogen atoms is large enough, the radiation emitted in the Lyman continuum is partially re-absorbed, causing the temperature inflection point observed in Figure 5(b). At lower speeds (see Figure 5(a)), the above temperature drop is not observed, as the shock strength is not sufficient to fully ionize the incoming gas. The trends shown in Figures 4 and 5 are similar to those reported in the work by Fadeyev & Gillet (1998, 2000).

The evolution of the free-electron temperature in Figure 5 shows a non-monotonic behavior in the precursor. The initial rise is due to energy deposition by photo-ionization (see Equation (35)). During this stage, elastic collisions between charged particles are unable to bring heavy particles and free electrons into equilibrium due to the low degree of ionization. This trend continues until the further increase of ionization causes elastic energy transfer to dominate over photo-ionization. As a result, the free-electron temperature reaches a maximum and then relaxes toward the heavy-particle temperature. The observed precursor behavior of the free-electron temperature is consistent with the findings of Foley & Clarke (1973), who investigated radiative shock waves in helium and argon plasmas. The only significant difference from the above reference is that, in the present work, free electrons are essentially in thermal equilibrium with heavy particles at the shock location. This is most probably due to the more efficient energy transfer in electron–atom collisions caused by the lower atomic weight of hydrogen (which is, respectively, 4 and 40 times lighter than helium and argon). This fact also has a strong influence on the length of the internal relaxation region (where elastic collisions play a crucial role), as recognized in the early theoretical and experimental work by Belozerov & Measures (1969).

The precursor heating for the whole gas becomes substantial only at large shock speeds, as indicated by the rise of the heavy-particle temperature and pressure in Figures 6(a) and (b), and by the volumetric radiant loss term plotted in Figure 6(c). It is interesting to note that at large speeds (i.e., 70 km s⁻¹), the free-electron temperature reaches a local minimum, and then rises again before reaching the shock. The observed behavior is due to the combined effect of the energy absorbed by radiation, which tends to increase the thermal energy of heavy particles, and elastic electron–proton collisions, which tend to keep the heavy-particle and free-electron temperatures together. The absorption of radiation in the precursor, in addition to altering the temperature and concentration profiles with respect to their free-stream values, also modifies the post-shock conditions due to the decrease and increase, respectively, of the Mach number and the pressure of the incoming gas. This is shown in detail in Table 1, which shows the pre- and post-shock conditions with and without radiation. At large speeds, the absorption of radiation reduces the post-shock temperature by more than 20%. As opposed to pressure and temperature, the mass density is quite insensitive to precursor radiation (see Table 1). This is a consequence of the (near) constancy of the flow velocity (not provided in Table 1), which can be explained by recalling that the absorbed radiation goes mainly into the thermal energy of the material gas. For a standing shock, global mass conservation requires that mass flux $\rho u$ must be constant, so that the density cannot vary appreciably when the flow velocity is changed by a small amount.

Zoom In Zoom Out Reset image size

Table 1.  Pressure, Temperature, and Density at the Pre- and Post-shock Locations Switching Radiation On/Off  (StS Model)

	Without Radiation			With Radiation
	Post-shock			Pre-shock			Post-shock
${u}_{\infty }\,(\mathrm{km}\,{{\rm{s}}}^{-1})$	$p\,(\mathrm{Pa})$	$T\,({\rm{K}})$	$\rho \,(\mathrm{kg}\,{{\rm{m}}}^{-3})$	$p\,(\mathrm{Pa})$	$T\,({\rm{K}})$	$\rho \,(\mathrm{kg}\,{{\rm{m}}}^{-3})$	$p\,(\mathrm{Pa})$	$T\,({\rm{K}})$	$\rho \,(\mathrm{kg}\,{{\rm{m}}}^{-3})$
40	144.086	40683.75	4.29 × 10−7	5.087	5034.509	1.21 × 10−7	142.79	40331.96	4.28 × 10−7
45	182.69	50347.42	4.40 × 10−7	5.22	5082.56	1.21 × 10−7	178.46	49209.53	4.38 × 10−7
50	225.83	6114.46	4.48 × 10−7	5.46	5171.12	1.21 × 10−7	214.91	58231.51	4.45 × 10−7
55	273.52	73076.09	4.54 × 10−7	5.89	5320.73	1.21 × 10−7	249.15	66616.85	4.45 × 10−7
60	325.74	86142.43	4.58 × 10−7	6.56	5548.98	1.22 × 10−7	278.84	73744.47	4.52 × 10−7
65	382.51	100343.87	4.62 × 10−7	7.50	5892.28	1.22 × 10−7	306.054	80208.88	4.54 × 10−7
70	443.82	115680.72	4.65 × 10−7	8.61	6322.32	1.22 × 10−7	333.62	86745.66	4.54 × 10−7

Note. The pre-shock pressure, temperature, and density are not provided in the case without radiation as these quantities are given by their free-stream LTE values (${p}_{\infty }=5\,{Pa}$, ${T}_{\infty }=5000\,K$, ${\rho }_{\infty }=1.21e-7\,\mathrm{kg}\,{{\rm{m}}}^{-3}$).

Download table as:  ASCIITypeset image

Before moving to the comparison with the reduced-order ME models discussed in Section 5.2, it is worth briefly analyzing the effect of including/excluding bb transitions from the calculations. Figure 7 shows a sample of this investigation for the temperatures and electron mole fraction profiles obtained for a shock propagating at 40 km s⁻¹. The inclusion of line radiation has the effect of shrinking the internal relaxation region and speeding up radiative recombination. The observed behavior originates from emission in the Balmer and Paschen lines, which are optically thin in the post-shock region, as shown in Figure 8, which shows the average monochromatic intensity at $x=2\,{\rm{m}}$ and $x=10\,{\rm{m}}$. The effect of an optically thin line (in particular the ${H}_{\alpha }$ line) is to provide an efficient channel to depopulate the higher bound electronic states. This is accompanied by a loss of energy which, in turn, induces the faster temperature drop observed in Figure 7(b). At the same time, in the radiative cooling region, the depopulation of high-lying states due to line emission indirectly favors radiative recombination because it makes the quantity $({n}_{i}\,{k}_{i}^{{\rm{PI}}}-{n}_{{\rm{e}}}\,{n}_{+}\,{k}_{i}^{{\rm{RR}}})$ more negative for these states (see Equation (29)). For large shock speeds (e.g., 60 km s⁻¹), the effects of bbtransitions are of less importance, due to the larger amount of photo-ionization in the precursor. This general trend has been observed in all cases investigated in this work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After analyzing the general features of the radiative shock waves considered in this work, predictions obtained by means of the ME models were systematically compared with the StS results. The comparison was first performed by excluding line radiation (Section 5.3), which was then added back into the final runs (Section 5.4).

Before discussing the results, it is worth saying a few words on how the grouping was done in practice. As anticipated at the end of Section 3.2, one should avoid lumping together states that are coupled by strong radiative transitions. This is particularly true for optically thin lines (e.g., ${\rm{H}}\alpha$, ${\rm{H}}\beta$). For this reason, the first and second energy groups (when two or more are used) were assigned to the n = 1 and n = 2 states, respectively, where n stands for the principal quantum number (Pauling & Wilson 1935, Chapter 5; Hubeny & Mihalas 2014, Chapter 7). No internal temperatures are used for these groups in the case of the MEL model, allowing for a further reduction of the number of unknowns. The higher states ($n\geqslant 3$) are lumped based on Equation (46) (or Equation (47) for the MEU model). In what follows, the notations MEU(k) and MEL(k) are used to indicate the reduced-order model used and the number of groups, k.

Figure 9 shows the temperature evolution for a shock propagating at 40 km s⁻¹ when using the MEL(3) model. The (only) internal temperature, indicated by the dotted–dashed line, refers to the states with $n\geqslant 3$. In the precursor, the high-lying levels are not in thermal equilibrium with free electrons since ${T}_{{\rm{i}}}\ne {T}_{{\rm{e}}}$. This behavior is due to the concurrent action of photo-ionization from high-lying states and de-excitation collisions (Foley & Clarke 1973). The internal temperature replicates the observed non-monotonic trend of that of the free electrons. When the incoming gas reaches the shock, the high-lying states are almost in thermal equilibrium with free electrons.

Zoom In Zoom Out Reset image size

The results of Figure 9 were compared with the StS predictions in Figure 10 for the electron mole fraction and the heavy-particle and free-electron temperatures. The MEU(3) solution is also reported. The MEL(3) model is in excellent agreement with the StS results. The two solutions are essentially indistinguishable. These results indicate that using three energy groups plus one internal temperature is sufficient to achieve an accurate prediction of radiative shock waves in hydrogen plasmas. In practical terms, this means that the number of unknowns is reduced by two orders of magnitude, which makes the MEL model very attractive for Computational Fluid Dynamics (CFDs) applications. On the other hand, the MEU(3) predictions are in strong disagreement with the StS results.

Zoom In Zoom Out Reset image size

To gain more insight from the comparison of Figure 10, the normalized populations of the internal energy states, n_i/g_i, have been extracted at the locations $x=-10,2$ and $20\,{\rm{m}}$ (see Figures 11–13). The first location (Figure 11) refers to the near precursor, the second (Figure 12) to the internal relaxation region, and the third (Figure 13) to the radiative cooling region. In the same figures, the corresponding average monochromatic intensity is also provided. For the MEL and MEU models, the population distributions have been reconstructed based on the group number densities and internal temperatures using Equations (46) and (47), respectively. In the present work, the normalized populations, n_i/g_i, have been preferred over the conventional population ratios (Hubeny & Mihalas 2014, Chapter 9), ${n}_{i}/{n}_{i}^{* }$ (where the symbol ${n}_{i}^{* }$ denotes the Saha equilibrium (10) occupation numbers), for the following reason. The population ratios give information only about the deviation from the local Saha distribution. On the other hand, the inspection of the ratio n_i/g_i on a semi-log plot (as done in Figures 11–13) allows for a deeper understanding of NLTE effects since it shows by how much the energy population deviates from the local Boltzmann distribution. This is immediately realized by taking the log of Equation (46):

$$\begin{eqnarray}&&\mathrm{ln}\left(\displaystyle \frac{{n}_{i}}{{g}_{i}}\right)=\mathrm{ln}\left[\displaystyle \frac{{\tilde{n}}_{k}}{{\tilde{Z}}_{k}^{}({T}_{k})}\right]-\displaystyle \frac{{E}_{i}}{{k}_{{\rm{B}}}{T}_{k}},\end{eqnarray} \tag{ 51 }$$

$k\in {{ \mathcal I }}_{k}$, $k\in { \mathcal K }$. Equation (51) shows that in the plot $\mathrm{ln}({n}_{i}/{g}_{i})$ versus E_i, the equilibrium population follows a straight line whose slope is indirectly proportional to the local temperature. This provides an effective way to quantify the departure from the local equilibrium distribution. This fact is of particular importance for high-lying states, which are the ones where most of the ionization occurs. The results in Figures 11–13 show that the superior description of the MEL model lies in its ability to (almost) replicate the StS behavior of the high-lying states. This fact is of great importance, especially for the ionization within the internal relaxation region. The MEU model gives a poor description of the population dynamics due to the assumption of infinite internal temperatures. As already stated in Section 3.1, this hypothesis does not allow the retrieval of the equilibrium (i.e., Boltzmann distribution) and has a strong impact on the temperature evolution in the radiative relaxation region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this zone, where the degree of ionization decreases due to radiative recombination, the collisional rate among the gas particles is large enough to maintain equilibrium at the local free-electron temperature (see Figure 13(a)). The same holds true for the radiation field only in the Lyman continuum, as shown by the observed departure from the local Planck function, ${B}_{\lambda }({T}_{{\rm{e}}})$, in the Balmer and Paschen continuum (see Figure 13(b)). The use of a uniform grouping prevents obtaining a Boltzmann distribution and leads to a higher population of the high-lying states and an underpredicted population for the n = 2 and n = 3 states (see Figure 13(a)). This produces lower temperatures compared to the MEL(3) and StS solutions. The observed difference is not negligible and is also enhanced by the large statistical weight of the high-lying states.⁹ The poor description of the level dynamics by the MEU(3) model also has an adverse effect on the predicted radiative signature of the plasma. This is demonstrated by the comparison in terms of the average monochromatic intensity given in Figures 11(b), 12(b), and 13(b). At all locations, the difference is higher at the Lyman and Balmer edges. The above findings on the MEU(3) and MEL(3) models have also been obtained when analyzing cases at larger shock speeds. For the sake of brevity, these results have not been added to the manuscript.

Before including the effects of line radiation on the reduced-order models (see Section 5.4), it was decided that the sensitivity of the solution to the number of groups will be investigated. The study included both the MEU and MEL models.

Figure 14 shows the free-electron mole fraction and all temperatures for the MEL(1), MEL(2), and MEL(3) models. When all the hydrogen bound states are lumped into one group (i.e., MEL(1) model), it is not possible to account for the under-population experienced by the high-lying states in the precursor and within the internal relaxation region. As a matter of fact, the MEL(1) model forces the high-lying states to be in thermal equilibrium with free electrons, resulting in a larger ionization rate (see Figures 14(a) and (b)). It is worth mentioning that, despite the use of one group, the MEL(1) model is still a NLTE model since (i) the hydrogen number density is not prescribed via Saha's equation (10), and (ii) the group internal temperature ${T}_{{\rm{i}}}$ is allowed to depart from the translational temperature of heavy particles and free electrons. The inclusion of one additional group (i.e., MEL(2) model) for the states with $n\geqslant 2$ greatly improves the solution and leads essentially to the same result obtained with the MEL(3) model. The MEL(3) model, however, gives a superior description of the population distribution as demonstrated by the comparison between the three models in Figure 15. It is worth noting that the lower accuracy for the population of the n = 2 and n = 3 states for the MEL(2) model plays a negligible role on gas quantities such as chemical composition and temperatures. This is no longer true when line radiation is taken into account (see Section 5.4). One may conclude that using only two energy groups is already enough to obtain an accurate prediction when neglecting line radiation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 16 illustrates the results of the sensitivity study for the MEU model. In analogy with the MEL model, increasing the number of groups improves the solution. However, in the present case, this is achieved with a larger number of groups (11), which makes the MEU model less attractive for CFD applications. The use of 11 groups allows the accurate prediction of both temperature and species concentrations in the precursor and internal relaxation region. However, some discrepancies with the StS results still appear in the radiative cooling region. As demonstrated earlier, these are caused by the intrinsic defect of the MEU model, namely, the impossibility of retrieving a Boltzmann distribution for the bound states. In view of these results, the MEU model has not been considered in the calculations performed, including the effects of line radiation as discussed in Section 5.4.

Zoom In Zoom Out Reset image size

The sensitivity analysis at the end of Section 5.3 has shown that two energy groups are sufficient to obtain an accurate solution when accounting only for continuum radiation. To test whether this finding is retrieved when atomic lines are added, the StS and MEL(2) solutions have been compared in Figure 17 for the electron mole fraction and heavy-particle and free-electron temperatures. In the same image, the MEL(3) solution is also shown. Overall, the MEL(2) solution now shows a sensible departure from the StS prediction, the difference being larger in the internal relaxation region (see Figures 17(a) and (b)). In the precursor, the MEL(2) model slightly overestimates photo-ionization, while in the radiative cooling region the disagreement with the StS solution is barely noticeable. The MEL(3) model is again in excellent agreement with the StS prediction (see Figures 17(c) and (d)). The performance degradation of the MEL(2) model is not surprising and is the result of grouping together the n = 2 and n = 3 states. This annihilates the effect of radiative decay through the optically thin ${\rm{H}}\alpha$ line, which, as shown in Figure 7, has a non-negligible impact on the internal relaxation and radiative cooling regions. On the other hand, in the MEL(3) model, the n = 2 and n = 3 states are grouped separately, thus allowing for a more accurate prediction of bb radiative losses. The above qualitative arguments are confirmed when comparing the StS and MEL(2-3) population distributions as done in Figure 18. The above picture refers to the location $x=2\,{\rm{m}}$ behind the shock (internal relaxation region), and for the sake of completeness, also shows results obtained when neglecting line radiation. It is readily seen that, when the n = 2 and n = 3 states are lumped together, the accuracy of the reconstructed distribution strongly degrades. In particular, the population of the above states (and also those close to the ionization limit) is underpredicted by at least one order of magnitude. It is worth noting that the observed departure between the StS and MEL(3) predictions for the population of the n = 3 and n = 4 states (see Figure 18(b)) has essentially no effect on the solution accuracy. This can be explained by recalling that in general, the Einstein coefficient for the ${\rm{H}}\alpha$ line is one order of magnitude larger than of the ${\rm{H}}\beta$ and ${\rm{H}}\gamma$ lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 19 compares the average monochromatic intensity around the ${\rm{H}}\alpha$, ${\rm{H}}\beta$, and ${\rm{H}}\gamma$ lines predicted by the StS and MEL(2-3) models at the location $x=10\,{\rm{m}}$ (radiative cooling region). Despite the good agreement observed for gas quantities (such as the electron mole fraction shown in Figure 17), the MEL(2) model systematically overestimates the peak value of the intensity of all lines. The MEL(3) model is, on the other hand, in excellent agreement with the StS results.

Zoom In Zoom Out Reset image size

The accuracy of the MEL(3) model has been further confirmed by repeating the calculations treated in this section at larger shock speeds. This is demonstrated in Figure 20, which compares the StS and MEL(3) electron mole fraction profiles for increasing shock speeds. In analogy with what is observed for the 40 km s⁻¹ case, the two solutions overlap.

Zoom In Zoom Out Reset image size

The gas pressure is always given by Dalton's law, $p={n}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}}+{n}_{{\rm{h}}}{k}_{{\rm{B}}}{T}_{{\rm{h}}}$, where the heavy-particle number density is obtained by summing the group and proton contributions, ${n}_{{\rm{h}}}={\sum }_{k\in { \mathcal K }}{\tilde{n}}_{k}+{n}_{+}$. In view of Equation (46), the heavy-particle thermal energy density (8) becomes

$$\begin{eqnarray}&&\rho {e}_{{\rm{h}}}=\displaystyle \frac{3}{2}{p}_{{\rm{h}}}+\displaystyle \sum _{k\in { \mathcal K }}{\tilde{n}}_{k}\langle {E}_{i}{\rangle }_{k}+{n}_{+}{E}_{+}.\end{eqnarray} \tag{ 53 }$$

The group mass production terms due to collisional excitation and ionization are obtained via the substitution of Equation (46) into Equation (13) and then summing the result obtained over all of the levels within each group (i.e., ${\tilde{\omega }}_{k}={\sum }_{i\in {{ \mathcal I }}_{k}}{\omega }_{i}$). After some algebraic manipulation, one obtains

$$\begin{eqnarray}{\omega }_{{\rm{e}}}^{\mathrm{col}} & = & \displaystyle \sum _{k\in { \mathcal K }}{n}_{{\rm{e}}}\,[{\tilde{n}}_{k}\,{\tilde{k}}_{k}^{{\rm{I}}}({T}_{{\rm{e}}},{T}_{k})-{n}_{{\rm{e}}}\,{n}_{+}\,{\tilde{k}}_{k}^{{\rm{R}}}({T}_{{\rm{e}}})],\\ {\omega }_{+}^{\mathrm{col}} & = & {\omega }_{{\rm{e}}}^{\mathrm{col}},\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}{\tilde{\omega }}_{k}^{\mathrm{col}} & = & -\,{n}_{{\rm{e}}}\,[{\tilde{n}}_{k}\,{\tilde{k}}_{k}^{{\rm{I}}}({T}_{{\rm{e}}},{T}_{k})-{n}_{{\rm{e}}}\,{n}_{+}{\tilde{k}}_{k}^{{\rm{R}}}({T}_{{\rm{e}}})]\\ & & -\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{l\in { \mathcal K }}{l\ne k}}{n}_{{\rm{e}}}\,[{\tilde{n}}_{k}\,{\tilde{k}}_{{kl}}^{{\rm{E}}}({T}_{{\rm{e}}},{T}_{k})-{\tilde{n}}_{l}\,{\tilde{k}}_{{lk}}^{{\rm{E}}}({T}_{{\rm{e}}},{T}_{l})],\end{eqnarray} \tag{ 55 }$$

$k\in { \mathcal K }$, where the group endothermic rate coefficients for excitation and ionization, and those for the related exothermic de-excitation and three-body recombination processes, are

$$\begin{eqnarray}&&{\tilde{k}}_{k}^{{\rm{I}}}({T}_{{\rm{e}}},{T}_{k})=\langle {k}_{i}^{{\rm{I}}}({T}_{{\rm{e}}}){\rangle }_{k},\quad \quad \,\\ &&{\tilde{k}}_{k}^{{\rm{R}}}({T}_{{\rm{e}}})={\tilde{k}}_{k}^{{\rm{I}}}({T}_{{\rm{e}}},{T}_{{\rm{e}}})\displaystyle \frac{{Z}_{{\rm{H}}}^{{\rm{t}}}({T}_{{\rm{e}}})\,{\tilde{Z}}_{k}^{}({T}_{{\rm{e}}})}{{g}_{{\rm{e}}}\,{Z}_{{\rm{e}}}^{{\rm{t}}}({T}_{{\rm{e}}})\,{g}_{+}\,{Z}_{{{\rm{H}}}^{+}}^{{\rm{t}}}({T}_{{\rm{e}}})}\exp \left(\displaystyle \frac{{E}_{+}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right),\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}{\tilde{k}}_{{kl}}^{{\rm{E}}}({T}_{{\rm{e}}},{T}_{k}) & = & \langle\langle {k}_{{ij}}^{{\rm{E}}}({T}_{{\rm{e}}}){ \rangle\rangle }_{{kl}},\\ {\tilde{k}}_{{lk}}^{{\rm{E}}}({T}_{{\rm{e}}},{T}_{l}) & = & \langle\langle {k}_{{ji}}^{{\rm{E}}}({T}_{{\rm{e}}}){ \rangle\rangle }_{{lk}},\end{eqnarray} \tag{ 57 }$$

$k,l\in { \mathcal K }$, where $k\lt l$ in Equation (57). It is worth noting that the group rate coefficients for ionization and three-body recombination (56) are related by a detailed balance relation among the groups. The same holds true for excitation/de-excitation only in the case of thermal equilibrium with free electrons (i.e., ${T}_{k}={T}_{{\rm{e}}},\forall \,k\in { \mathcal K }$).

By applying the same procedure as above, one arrives at the following relations for the group energy transfer source terms due to ionization and excitation:

$$\begin{eqnarray}{\tilde{{\rm{\Omega }}}}_{k}^{{\rm{I}}} & = & -\displaystyle \sum _{k\in { \mathcal K }}{n}_{{\rm{e}}}\,[{\tilde{n}}_{k}\,{\tilde{G}}_{k}^{{\rm{I}}}({T}_{{\rm{e}}},{T}_{k})-{n}_{{\rm{e}}}\,{n}_{+}{\tilde{G}}_{k}^{{\rm{R}}}({T}_{{\rm{e}}})],\\ \quad {\tilde{{\rm{\Omega }}}}_{k}^{{\rm{E}}} & = & -\displaystyle \sum _{k,l\in { \mathcal K }}{n}_{{\rm{e}}}\,[{\tilde{n}}_{k}\,{\tilde{G}}_{{kl}}^{{{\rm{E}}}_{f}}({T}_{{\rm{e}}},{T}_{k})-{\tilde{n}}_{l}\,{\tilde{G}}_{{lk}}^{{{\rm{E}}}_{b}}({T}_{{\rm{e}}},{T}_{l})],\end{eqnarray} \tag{ 58 }$$

$k\in { \mathcal K }$, where the related transfer rates due to ionization, recombination, and excitation/de-excitation are

$$\begin{eqnarray}&&{\tilde{G}}_{k}^{{\rm{I}}}({T}_{{\rm{e}}},{T}_{k})=\langle {E}_{i}\,{k}_{i}^{{\rm{I}}}({T}_{{\rm{e}}}){\rangle }_{k},\quad \quad \,\\ &&{\tilde{G}}_{k}^{{\rm{R}}}({T}_{{\rm{e}}})={\tilde{G}}_{k}^{{\rm{I}}}({T}_{{\rm{e}}},{T}_{{\rm{e}}})\,\displaystyle \frac{{Z}_{{\rm{H}}}^{{\rm{t}}}({T}_{{\rm{e}}})\,{\tilde{Z}}_{k}^{}({T}_{{\rm{e}}})}{{g}_{{\rm{e}}}\,{Z}_{{\rm{e}}}^{{\rm{t}}}({T}_{{\rm{e}}})\,{g}_{+}\,{Z}_{{{\rm{H}}}^{+}}^{{\rm{t}}}({T}_{{\rm{e}}})}\exp \left(\displaystyle \frac{{E}_{+}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right),\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}{\tilde{G}}_{{kl}}^{{{\rm{E}}}_{f}}({T}_{{\rm{e}}},{T}_{k}) & = & \langle\langle {E}_{i}\,{k}_{{ij}}^{{\rm{E}}}({T}_{{\rm{e}}}){ \rangle\rangle }_{{kl}},\\ {\tilde{G}}_{{lk}}^{{{\rm{E}}}_{b}}({T}_{{\rm{e}}},{T}_{l}) & = & \langle\langle {E}_{i}\,{k}_{{ji}}^{{\rm{E}}}({T}_{{\rm{e}}}){ \rangle\rangle }_{{lk}},\end{eqnarray} \tag{ 60 }$$

$k,l\in { \mathcal K }$. The use of the Boltzmann grouping relation (46) in Equation (18) allows the ionization and excitation energy transfer terms for the free-electron gas to be rewritten as

$$\begin{eqnarray}&&{\tilde{{\rm{\Omega }}}}_{{\rm{e}}}^{{\rm{I}}}=-\displaystyle \sum _{k\in { \mathcal K }}{\tilde{{\rm{\Omega }}}}_{k}^{{\rm{I}}}-{E}_{+}\,{\omega }_{+}^{{\rm{I}}},\quad {\tilde{{\rm{\Omega }}}}_{{\rm{e}}}^{{\rm{I}}}=-\displaystyle \sum _{k\in { \mathcal K }}{\tilde{{\rm{\Omega }}}}_{k}^{{\rm{E}}},\end{eqnarray} \tag{ 61 }$$

where ${\omega }_{+}^{{\rm{I}}}={\omega }_{+}^{\mathrm{col}}$.

bb transitions—The substitution of Equation (46) into Equation (22) allows the monochromatic emission and absorption coefficients due to bb radiation to be rewritten as

$$\begin{eqnarray}{\varepsilon }_{\lambda }^{{\rm{BB}}} & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{k,l\in { \mathcal K }}{l\geqslant k}}\displaystyle \frac{{h}_{{\rm{P}}}c}{4\pi }\,{\tilde{n}}_{l}\,{\tilde{A}}_{{lk}}(\lambda ,{T}_{l}),\\ {\kappa }_{\lambda }^{{\rm{BB}}} & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{k,l\in { \mathcal K }}{l\geqslant k}}\displaystyle \frac{{h}_{{\rm{P}}}}{4\pi }\,[{\tilde{n}}_{k}\,{\tilde{B}}_{{kl}}(\lambda ,{T}_{k})-{\tilde{n}}_{l}\,{\tilde{B}}_{{lk}}(\lambda ,{T}_{l})],\end{eqnarray} \tag{ 62 }$$

where the group wavelength and temperature-dependent Einstein coefficients are defined as

$$\begin{eqnarray}&&{\tilde{A}}_{{lk}}(\lambda ,{T}_{l})= \langle\langle {A}_{{ji}}\,{\phi }_{\lambda }^{{ji}}{ \rangle\rangle }_{{lk}},\quad {\tilde{B}}_{{lk}}(\lambda ,{T}_{l})= \langle\langle {\lambda }_{{ji}}\,{B}_{{ji}}\,{\phi }_{\lambda }^{{ji}}{ \rangle\rangle }_{{lk}},\\ &&{\tilde{B}}_{{kl}}(\lambda ,{T}_{k})= \langle\langle {\lambda }_{{ij}}\,{B}_{{ij}}\,{\phi }_{\lambda }^{{ji}}{ \rangle\rangle }_{{kl}},k\leqslant l\end{eqnarray} \tag{ 63 }$$

$k,l\in { \mathcal K }$, where the obvious symmetry relations satisfied by the line emission/absorption profile and wavelength (i.e., ${\phi }_{\lambda }^{{ij}}={\phi }_{\lambda }^{{ji}}$ and ${\lambda }_{{ij}}={\lambda }_{{ji}}$, respectively) are used. The application of the above procedure to Equations (24) and (26) gives, respectively, the group mass production rates and the volumetric rate of energy loss of matter energy due to bb radiation:

$$\begin{eqnarray}{\tilde{\omega }}_{k}^{{\rm{BB}}} & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{l\in { \mathcal K }}{l\gt k}}\,\{{\tilde{n}}_{l}\,{\tilde{A}}_{{lk}}^{{\rm{m}}}({T}_{l})-[{\tilde{n}}_{k}\,{\tilde{B}}_{{kl}}^{{\rm{m}}}({T}_{k})-{\tilde{n}}_{l}\,{\tilde{B}}_{{lk}}^{{\rm{m}}}({T}_{l})]\}\\ & & -\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{l\in { \mathcal K }}{l\lt k}}\,\{{\tilde{n}}_{k}\,{\tilde{A}}_{{kl}}^{{\rm{e}}}({T}_{k})-[{\tilde{n}}_{l}\,{\tilde{B}}_{{lk}}^{{\rm{e}}}({T}_{l})-{\tilde{n}}_{k}\,{\tilde{B}}_{{kl}}^{{\rm{e}}}({T}_{k})]\},\\ & & \times k\in { \mathcal K },\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}^{{\rm{BB}}}=\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{k,l\in { \mathcal K }}{l\geqslant k}}\,\{{\tilde{n}}_{l}\,{\tilde{A}}_{{lk}}^{{\rm{e}}}({T}_{l})-[{\tilde{n}}_{k}\,{\tilde{B}}_{{kl}}^{{\rm{e}}}({T}_{k})-{\tilde{n}}_{l}\,{\tilde{B}}_{{lk}}^{{\rm{e}}}({T}_{l})]\},\end{eqnarray} \tag{ 65 }$$

where the group mass production and energy transfer Einstein coefficients are given by the following Boltzmann averages:

$$\begin{eqnarray}&&{\tilde{A}}_{{lk}}^{{\rm{m}}}({T}_{l})= \langle\langle {A}_{{ji}}^{{\rm{m}}}{ \rangle\rangle}_{{lk}},\quad {\tilde{A}}_{{lk}}^{{\rm{e}}}({T}_{l})= \langle\langle {A}_{{ji}}^{{\rm{e}}}{ \rangle\rangle }_{{lk}},\end{eqnarray} \tag{ 66 }$$

$$\begin{eqnarray}&&{\tilde{B}}_{{lk}}^{{\rm{m}}}({T}_{l})= \langle\langle {B}_{{ji}}^{{\rm{m}}}{ \rangle\rangle }_{{lk}},\quad {\tilde{B}}_{{lk}}^{{\rm{e}}}({T}_{l})= \langle\langle {B}_{{ji}}^{{\rm{e}}}{ \rangle\rangle }_{{lk}},\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray}&&{\tilde{B}}_{{kl}}^{{\rm{m}}}({T}_{k})= \langle\langle {B}_{{ij}}^{{\rm{m}}}{ \rangle\rangle }_{{kl}},\quad {\tilde{B}}_{{kl}}^{{\rm{e}}}({T}_{k})= \langle\langle {B}_{{ij}}^{{\rm{e}}}{ \rangle\rangle }_{{kl}},\end{eqnarray} \tag{ 68 }$$

$k,l\in { \mathcal K }$, where $k\leqslant l$.

bf/fb transitions—The substitution of Equation (46) into Equations (27) and (28) allows the monochromatic emission and absorption coefficients due to bf and fb radiation to be rewritten as

$$\begin{eqnarray}{\varepsilon }_{\lambda }^{{\rm{FB}}} & = & \displaystyle \frac{{h}_{{\rm{P}}}^{4}\,{c}^{2}\,{n}_{{\rm{e}}}\,{n}_{+}}{{\lambda }^{5}{(2\pi {m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}{g}_{+}}\\ & & \times \exp \left(\displaystyle \frac{{E}_{+}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)\displaystyle \sum _{k\in { \mathcal K }}{\tilde{Z}}_{k}^{}({T}_{{\rm{e}}})\,{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{{\rm{e}}}),\end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray}{\kappa }_{\lambda }^{{\rm{BF}}} & = & \displaystyle \sum _{k\in { \mathcal K }}{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{k})\left[{\tilde{n}}_{k}-\displaystyle \frac{1}{2}\displaystyle \frac{{h}_{{\rm{P}}}^{3}\,{n}_{{\rm{e}}}\,{n}_{+}}{{(2\pi {m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}}\displaystyle \frac{{\tilde{Z}}_{k}^{}({T}_{{\rm{e}}})}{{g}_{+}}\right.\\ & & \times \left.\displaystyle \frac{{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{{\rm{e}}})}{{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{k})}\exp \left(\displaystyle \frac{{E}_{+}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)\right],\end{eqnarray} \tag{ 70 }$$

where the quantity ${\tilde{\sigma }}_{k}^{{\rm{PI}}}$ stands for the temperature-dependent group photo-ionization cross-section. The former is computed via the Boltzmann average ${\tilde{\sigma }}_{k}^{{\rm{PI}}}=\langle {\sigma }_{i}^{{\rm{PI}}}{\rangle }_{k}$. The notation ${\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{{\rm{e}}})$ indicates that the Boltzmann average has to be evaluated at the free-electron temperature.

The mass production terms are

$$\begin{eqnarray}{\omega }_{{\rm{e}}}^{{\rm{BF}}/{\rm{FB}}} & = & \displaystyle \sum _{k\in { \mathcal K }}[{\tilde{n}}_{k}\,{\tilde{k}}_{k}^{{\rm{PI}}}({J}_{\lambda },{T}_{k})-{n}_{{\rm{e}}}\,{n}_{+}\,{\tilde{k}}_{k}^{{\rm{RR}}}({T}_{{\rm{e}}},{J}_{\lambda })],\\ {\omega }_{+}^{{\rm{BF}}/{\rm{FB}}} & = & {\omega }_{{\rm{e}}}^{{\rm{BF}}/{\rm{FB}}},\end{eqnarray} \tag{ 71 }$$

$$\begin{eqnarray}&&{\tilde{\omega }}_{k}^{{\rm{BF}}/{\rm{FB}}}=-[{\tilde{n}}_{k}\,{\tilde{k}}_{k}^{{\rm{PI}}}({J}_{\lambda },{T}_{k})-{n}_{{\rm{e}}}\,{n}_{+}\,{\tilde{k}}_{k}^{{\rm{RR}}}({T}_{{\rm{e}}},{J}_{\lambda })],\end{eqnarray} \tag{ 72 }$$

$k\in { \mathcal K }$, where the group rate coefficients for photo-ionization and radiative recombination read ${\tilde{k}}_{k}^{{\rm{PI}}}=\langle {k}_{i}^{{\rm{PI}}}{\rangle }_{k}$ and ${\tilde{k}}_{k}^{{\rm{RR}}}={\sum }_{i\in {{ \mathcal I }}_{k}}{k}_{i}^{{\rm{RR}}}$, respectively. By setting the elementary photo-ionization cross-sections to zero for wavelengths above the photo-ionization limit (i.e., ${\sigma }_{i}^{{\rm{PI}}}=0$ for $\lambda \gt {\lambda }_{+i}$), it is possible to extend the wavelength integrals to infinity and exchange the order between summation and integration when computing the rate coefficients ${\tilde{k}}_{k}^{{\rm{PI}}}$ and ${\tilde{k}}_{k}^{{\rm{RR}}}$. Performing these operations, one obtains

$$\begin{eqnarray}{\tilde{k}}_{k}^{{\rm{PI}}}({J}_{\lambda },{T}_{k}) & = & \displaystyle \frac{4\pi }{{h}_{{\rm{P}}}c}{\displaystyle \int }_{0}^{+\infty }{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{k})\,{J}_{\lambda }\,\lambda \,d\lambda ,\\ \quad {\tilde{k}}_{k}^{{\rm{RR}}}({T}_{{\rm{e}}},{J}_{\lambda }) & = & {\tilde{k}}_{k}^{{\rm{RR}} \mbox{-} {\rm{s}}}({T}_{{\rm{e}}})+{\tilde{k}}_{k}^{{\rm{RR}} \mbox{-} {\rm{i}}}({T}_{{\rm{e}}},{J}_{\lambda }),\end{eqnarray} \tag{ 73 }$$

$$\begin{eqnarray}{\tilde{k}}_{k}^{{\rm{RR}} \mbox{-} {\rm{s}}}({T}_{{\rm{e}}}) & = & \sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{{\tilde{Z}}_{k}^{}({T}_{{\rm{e}}})}{{g}_{+}}\displaystyle \frac{{h}_{{\rm{P}}}^{3}c}{{({m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}}\\ & & \times \exp \left(\displaystyle \frac{{E}_{+}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right){\displaystyle \int }_{0}^{+\infty }\displaystyle \frac{{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{{\rm{e}}})}{{\lambda }^{4}}\\ & & \times \exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)d\lambda ,\end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray}{\tilde{k}}_{k}^{{\rm{RR}}-{\rm{k}}}({T}_{{\rm{e}}},{J}_{\lambda }) & = & \sqrt{\displaystyle \frac{1}{2\pi }}\displaystyle \frac{{\tilde{Z}}_{k}^{}({T}_{{\rm{e}}})}{{g}_{+}}\displaystyle \frac{{h}_{{\rm{P}}}^{2}}{{({m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}c}\\ & & \times \exp \left(\displaystyle \frac{{E}_{+}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right){\displaystyle \int }_{0}^{+\infty }{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{{\rm{e}}})\,\lambda \,{J}_{\lambda }\\ & & \times \exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)d\lambda ,\end{eqnarray} \tag{ 75 }$$

$k\in { \mathcal K }$. Equations (73)–(75) show that the group rate coefficients for photo-ionization and radiative recombination are formally identical to those for the StS model (32) and (33) and that the former can be obtained from the latter by replacing the statistical weights with the group partition functions, and the elementary cross-sections with the group cross-sections.

The group energy transfer terms due to photo-ionization and radiative recombination are

$$\begin{eqnarray}&&{\tilde{{\rm{\Omega }}}}_{k}^{{\rm{BF}}/{\rm{FB}}}=-[{\tilde{n}}_{k}\,{\tilde{G}}_{k}^{{\rm{PI}}}({J}_{\lambda },{T}_{k})-{n}_{{\rm{e}}}\,{n}_{+}\,{\tilde{G}}_{k}^{{\rm{RR}}}({J}_{\lambda },{T}_{{\rm{e}}})],\end{eqnarray} \tag{ 76 }$$

$k\in { \mathcal K }$, where energy transfer rates due to photo-ionization and radiative recombination are ${\tilde{G}}_{k}^{{\rm{PI}}}=\langle {E}_{i}\,{k}_{i}^{{\rm{PI}}}{\rangle }_{k}$ and ${\tilde{G}}_{k}^{{\rm{RR}}}\,={\sum }_{i\in {{ \mathcal I }}_{k}}{E}_{i}\,{k}_{i}^{{\rm{RR}}}$, respectively.

By repeating the above procedure, the net volumetric energy loss rates for the material gas (34) and free electrons (35) become

$$\begin{eqnarray}{{\rm{\Omega }}}^{{\rm{BF}}/{\rm{FB}}} & = & \displaystyle \sum _{k\in { \mathcal K }}\left[\sqrt{\displaystyle \frac{1}{2\pi }}\displaystyle \frac{{\tilde{Z}}_{k}^{}({T}_{{\rm{e}}})}{{g}_{+}}\displaystyle \frac{{h}_{{\rm{P}}}^{3}\,{n}_{{\rm{e}}}\,{n}_{+}}{{({m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}}\right.\\ & & \times \exp \left(\displaystyle \frac{{E}_{+}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right){\displaystyle \int }_{0}^{+\infty }{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{{\rm{e}}})\,{J}_{\lambda }\\ & & \times \exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)d\lambda \\ & & -\left.4\pi \,{\tilde{n}}_{k}{\displaystyle \int }_{0}^{+\infty }{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{k})\,{J}_{\lambda }\,d\lambda \right]\\ & & +\displaystyle \sum _{k\in { \mathcal K }}\sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{{\tilde{Z}}_{k}^{}({T}_{{\rm{e}}})}{{g}_{+}}\\ & & \times \displaystyle \frac{{h}_{{\rm{P}}}^{4}{c}^{2}\,{n}_{{\rm{e}}}\,{n}_{+}}{{({m}_{{\rm{e}}}{k}_{{\rm{B}}}{T}_{{\rm{e}}})}^{3/2}}\exp \left(\displaystyle \frac{{E}_{+}}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}}\right){\displaystyle \int }_{0}^{+\infty }\displaystyle \frac{{\tilde{\sigma }}_{k}^{{\rm{PI}}}(\lambda ,{T}_{{\rm{e}}})}{{\lambda }^{5}}\\ & & \times \exp \left(-\displaystyle \frac{{h}_{{\rm{P}}}c}{{k}_{{\rm{B}}}{T}_{{\rm{e}}}\lambda }\right)d\lambda ,\end{eqnarray} \tag{ 77 }$$

$$\begin{eqnarray}{{\rm{\Omega }}}_{{\rm{e}}}^{{\rm{BF}}/{\rm{FB}}} & = & {{\rm{\Omega }}}^{{\rm{BF}}/{\rm{FB}}}+\displaystyle \sum _{k\in { \mathcal K }}[{\tilde{n}}_{k}\,{\tilde{G}}_{k}^{{\rm{PI}}}({J}_{\lambda },{T}_{k})\\ & & -{n}_{{\rm{e}}}\,{n}_{+}\,{\tilde{G}}_{k}^{{\rm{RR}}}({J}_{\lambda },{T}_{{\rm{e}}})]+{E}_{+}\,{\omega }_{+}^{{\rm{BF}}/{\rm{FB}}}.\end{eqnarray} \tag{ 78 }$$

ff transitions—The emission and absorption coefficients, and the net volumetric energy loss rate due to ff transitions are not affected by the grouping as only charge–charge interactions are taken into account in this work. Hence, Equations (36) and (37) are used in the same way as for the StS model.