© The Authors 2024

1 Introduction

The dynamical evolution of the diffuse matter toward the formation of dense structures is a fascinating problem of physics. Because of thermal instability (see Fig. 1), the diffuse neutral medium is known to be separated into two stable thermal states, a warm phase at ~ 10 000 K, the warm neutral medium (WNM), and a cold phase at ~80 K, the cold neutral medium (CNM) (e.g., Field 1965; Heiles & Troland 2003; Wolfire et al. 2003; Murray et al. 2018; Marchal et al. 2019). Although these two phases are found to be at thermal pressure equilibrium (Jenkins & Tripp 2011), they do not peacefully coexist. By inducing pressure variations and shear motions at all scales, magnetized turbulence drives transfers of mass and energy between the WNM and the CNM (e.g., Audit & Hennebelle 2005; Seifried et al. 2011). Under the combined actions of turbulence and thermal instability, the diffuse gas therefore constantly flows from one phase to the other, eventually leading to the formation of structures massive enough to trigger gravitational collapse.

One remarkable aspect is the role of interstellar shocks at each stage of this dynamical evolution. The driving of interstellar turbulence in galaxies occurs through a wealth of supersonic events, which include outflows from young stellar objects, infall of material onto the galactic disk, supernova explosions, outflows from active galactic nuclei, and galaxy collisions. All these events inevitably form shocks, which not only dissipate a fraction of the initial kinetic energy but also generate a turbulent cascade in their wake; for example, through the production of baroclynic vorticity (Elmegreen & Scalo 2004; Vazquez-Semadeni et al. 1996). On the theoretical side, this process is particularly well illustrated by simulations of supernova-driven turbulence (e.g., Pan et al. 2016). On the observational side, one of the most iconic examples of this process is the Stephan’s Quintet galaxy collision, where colliding galaxies at relative velocities of ~1000 km s⁻¹ generate a turbulent flow whose mechanical energy is traced by the rovibrational lines of H₂ (Guillard et al. 2009; Appleton et al. 2023). Interstellar shocks are therefore the fingerprints of the mechanisms by which mechanical energy is injected into the ISM.

By altering the physical state of matter, shocks propagating in the WNM are known to bring the gas in an unstable cooling state (red region in Fig. 6 of Falle et al. 2020). Hydrodynamical simulations show that if the postshock thermal pressure is higher than the maximum pressure of the WNM and the cooling time is shorter than the dynamical time of the perturbation, shocks inevitably induce phase transition between the WNM and the CNM (e.g., Hennebelle & Pérault 1999; Kupilas et al. 2021). This process is mitigated by the presence of a magnetic field. Magnetohydrodynamic simulations of colliding flows (Hennebelle & Pérault 2000; Inoue & Inutsuka 2008) and of interstellar shocks (Falle et al. 2020) show that the magnetic field efficiently prevents the formation of very dense gas and may force the final state to lie on the unstable part of the thermal equilibrium (green curve in Fig. 1). Nevertheless, shocks still induce phase transition in the sense that this unstable state subsequently splits into warm and cold environments (van Loo et al. 2010). Shocks may therefore play a substantial role in the mass transfer between the WNM and the CNM.

Finally, as the turbulence develops, the kinetic energy injected at large scale is transferred toward small scales. This turbulent cascade, which operates in a multiphase environment, involves a complex interplay between compressive and solenoidal modes (e.g., Porter et al. 2002; Vestuto et al. 2003; Pan et al. 2016) that leads to the formation of low-velocity shocks, which propagate in the WNM and the CNM and dissipate a substantial fraction of the kinetic energy (Stone et al. 1998; Lehmann et al. 2016; Lesaffre et al. 2020; Richard et al. 2022). Shocks could therefore be the fingerprints of the dissipation of the mechanical energy stored in interstellar turbulence.

While extensive studies have been performed to identify and interpret atomic and molecular tracers of low-velocity shocks (V_S ⩽ 30 km s⁻¹) propagating in the CNM (e.g., Draine & Katz 1986; Monteiro et al. 1988; Flower & Pineau des Forêts 1998; Lesaffre et al. 2013), little has been done so far regarding the observational tracers of shocks propagating in the WNM at low (V_s ⩽ 30 km s⁻¹), intermediate (30 < V_s < 100 km s⁻¹), and high velocities (V_s ⩾ 100 km s⁻¹) that induce phase transition between the WNM and the CNM. This oversight may be attributed to the numerical challenges associated with the problem. Indeed, one major difficulty is that intermediate- and high-velocity shocks produce high-energy photons that modify the state of the preshock gas and interact with the shock itself (Raymond 1979; Hollenbach & McKee 1989). Numerical codes capable of describing all the intricacies of self-irradiated shocks are rare.

The impact of ionizing radiation has long been included in numerical simulations to study the thermal instability induced by interstellar shocks (e.g., Sutherland et al. 2003) and more recently in simulations of supernova remnants (SNRs) (Sarkar et al. 2021a,b) and Galactic winds (Sarkar et al. 2022). However, as detailed simulations are numerically expensive, they are usually limited to specific cases. It follows that systematic studies involving the complete exploration of a large parameter domain can only be achieved with 1D steady-state models. To the best of our knowledge, the only existing public model that follows the physics of self-irradiated shocks is the MAPPINGS V code. Developed over the past 40 yr (Dopita 1978; Dopita & Sutherland 1995, 1996; Allen et al. 2008; Dopita et al. 2013; Sutherland & Dopita 2017), MAPPINGS is a self-consistent model designed to compute the thermochemical structure and emission of photoionized regions and fast atomic shocks. For instance, MAPPINGS has been successfully applied to study narrow line regions and active galactic nuclei (e.g., Best et al. 2000; Zhu et al. 2023), Herbig-Haro objects (e.g., Dopita & Sutherland 2017), and SNRs (e.g., Kopsacheili et al. 2020; Priestley et al. 2022).

Given the existence of a highly sophisticated and regularly updated model such as MAPPINGS, one might wonder why we would develop another code. In addition to the fact that it is always advisable to have several models at hand to describe the same physical problem, the Paris-Durham shock code has a few advantages for the study of shocks propagating in the WNM. First, the Paris-Durham shock code is a multifluid model capable of computing various types of magnetohydrodynamic stationary solutions with distinctive physical properties, including C-type, C*-type, and CJ-type shocks (Chernoff 1987; Roberge & Draine 1990; Godard et al. 2019). Second, the code was recently improved to take into account the impact of an external UV radiation field (Lesaffre et al. 2013; Godard et al. 2019), and therefore contains all the ingredients required to treat the phase transition process between the WNM and the CNM, which results from the thermal instability that naturally occurs in diffuse irradiated environments. Last, this code is designed to study the formation and excitation of molecules in shocks. Although this aspect is not addressed here, developing the model for higher velocities allows to identify molecular tracers of highvelocity shocks, which are naturally produced during the phase transition and the formation of CNM structures.

In the first paper of this series, we therefore introduce a new version of the Paris-Durham shock code capable of treating the physics of self-irradiated shocks. As a first application, we present a study of shocks propagating at velocities between 10 and 500 km s⁻¹ in the diffuse WNM irradiated by the standard interstellar radiation field. This work is a stepping stone toward two follow-up papers dedicated to specific observational tracers of these shocks (e.g., Godard et al. 2024).

The main physical ingredients of the code and the recent updates performed to model self-irradiated shocks are presented in Sec. 2 and Appendix A. The dynamical, thermal, and chemical structures of a prototypical shock and the induced phase transition between the WNM and the CNM are described in Sec. 3. The reprocessing of the input mechanical energy flux into line and continuum radiation and the corresponding heating and ionizing rates of the interstellar gas are shown and discussed in Sec. 4.

Fig. 1 Thermal equilibrium state of the diffuse interstellar gas obtained with the Paris-Durham code and displayed in a particle density versus thermal pressure diagram. The color code highlights the three branches of the equilibrium state: the dynamically stable WNM (red) and CNM (blue) phases at ~104 K and ~80 K, respectively, and the dynamically unstable branch (green) at intermediate temperatures. The light blue and red zones show the regions of the diagram where the cooling rate is larger or smaller, respectively, than the heating rate. Light gray lines are isothermal contours from 1 to 106 K (from bottom right to top left).

2 Physical ingredients and numerical method

The theoretical model presented in this work is based on the Paris-Durham shock code, a public numerical tool¹ built to follow the dynamical, thermal, and chemical structures of a multi-fluid plane-parallel interstellar shock at steady-state. Initially designed to study low-velocity shocks (V_S ⩽ 30 km s⁻¹) propagating in molecular environments, the code was recently improved to treat shocks irradiated by an external UV radiation field (Godard et al. 2019) and intermediate-velocity shocks at velocities up to 60 km s⁻¹ self-irradiated by UV photons produced through the collisional excitation of atomic hydrogen (Lehmann et al. 2020, 2022). The model used here is built upon these recent improvements with additional updates on the chemistry and the radiative transfer, presented below, that allow the computation of shocks propagating at velocities larger than 60 km s⁻¹ in ionized or partially ionized atomic regions such as the WNM.

The physics of photodissociating and photoionizing shocks is well understood (Dopita & Sutherland 2003) and has been the subject of numerous studies over the past 60 yr (e.g., Shull & McKee 1979; Raymond 1979; Hollenbach & McKee 1989; Dopita & Sutherland 1995, 1996; Allen et al. 2008; Sutherland & Dopita 2017). In a nutshell, high-velocity shocks produce UV, EUV, and soft X-ray photons that not only propagate downstream where they interact with the postshock gas but also propagate upstream where they interact with the preshock medium and induce the formation of a radiative precursor. Depending on the shock speed, the flux of ionizing photons escaping upstream may become large enough to fully ionize the preshock, effectively leading to the formation of an HII region. As it modifies the entrance conditions of the gas, the radiation field generated by the shock itself therefore alters the shock own thermodynamical structure. This brief summary shows that modeling high-velocity shocks requires to take into account all the chemical and cooling processes that occur in HII regions and to implement a numerical scheme that self-consistently computes the physical state of the radiative precursor and its feedback on the shock structure. Following Sutherland & Dopita (2017) and Lehmann et al. (2020), we apply here a simple iterative scheme. The thermochemical structure of a shock is first computed using the Paris-Durham shock code. The propagation of the photons generated by the shock is then solved using a postprocessing radiative-transfer algorithm. Both steps are repeated until convergence.

2.1 Parameters and initial conditions

We consider a plane-parallel shock propagating at a velocity V_S in a magnetized medium with a constant proton density $n_{\mathrm{H}}^0$ and a magnetic field B₀ set in the direction perpendicular to the direction of propagation². In addition to the radiation produced by the shock itself, we assume that the shock is externally irradiated by an isotropic UV and optical photon flux. This flux is set to the standard interstellar radiation field (ISRF, Mathis et al. 1983) and scaled with a parameter G₀. Given the typical distances between OB star associations in the solar neighborhood (~100 pc, e.g., Zari et al. 2018) and the typical density of the WNM (~0.5 cm⁻³), we neglect the absorption of the external radiation field by the surrounding medium in which the shock propagates (Bellomi et al. 2020). Neglecting also the absorption by the shock itself (see Sec. 5.3), the external UV and optical photon flux is therefore assumed constant throughout the entire model. The preshock and shock gas are finally assumed to be pervaded by cosmic ray particles with an equivalent H₂ ionization rate ${\zeta }_{{\mathrm{H}}_2}$. The main parameters of the model, the standard values adopted, and the range of values explored in this work are given in Table 1.

Within this geometry and under the steady-state approximation, the magnetohydrodynamics equations and the chemical network adopted in the model lead to a set of coupled first order differential equations that describe the dynamical, thermal, and chemical evolutions of a fluid particle during its trajectory from the preshock medium to the postshock gas (Flower 2010). This system is solved using the DVODE forward integration algorithm (Hindmarsh 1983) starting from a given position in the preshock with known initial conditions. The initial position of the fluid particle (i.e., its distance from the shock front) is chosen so that the dissociating and ionizing radiation field generated by the shock is entirely absorbed by the radiative precursor. To limit the number of free parameters, we finally assume that the gas at this position is at thermochemical equilibrium. While it entirely fixes the initial conditions of the problem, this choice is highly debatable. As an illustration, we display in Fig. 2 the equilibrium value of the electronic fraction, n(e⁻)/n_H, and the timescale required to reach thermochemical equilibrium in the diffuse interstellar medium as functions of the density of the gas. This figure shows that the WNM (0.1 cm⁻³ ⩽ $n_{\mathrm{H}}^0$ ⩽ 2 cm⁻³) requires ~10 Myr to reach thermochemical equilibrium, a value larger than the dynamical evolution timescale of the diffuse matter (turnover timescale and evaporation and condensation timescales of the turbulent ISM) or the typical timescale separating the arrival of two successive SNRs at a given point in the Galaxy (Draine 2011).

Fig. 2 Timescale required to reach the ionization, thermal, and chemical equilibrium states in the diffuse medium as a function of the proton density nH (red curve), and corresponding equilibrium values of the ionization fraction, xe = n(e−)/nH (blue curve). The equilibrium timescale is estimated as the time at which the kinetic temperature, the electronic fraction, and the chemical abundances of all species are within 10% of their equilibrium values, starting from random initial conditions.

2.2 Dust grains and chemical network

The chemical composition of a fluid particle across an interstellar shock is computed taking into account gas, dust, and PAHs. The dust-to-gas mass ratio is set to 0.01 and the dust size distribution is assumed to follow the classical MRN distribution (Mathis et al. 1977) with a power law index of −3.5. The fractional abundance of PAHs, n_PAH/n_H, is set to 10⁻⁶, which corresponds to standard conditions for unshielded diffuse environments (Draine & Li 2007). Oppositely to our previous work (Godard et al. 2019) we neglect here the formation and destruction of grain mantles. For simplification, we also neglect shattering, vaporization, and coagulation of grain cores (Jones et al. 1996). The grain size distribution is therefore solely affected by erosion through thermal sputtering (Sect. 2.2 and Appendix A of Godard et al. 2019). This is a strong limitation of the model as shattering and vaporization of large grains are expected to occur in high-velocity shocks (see for instance Fig. 7 of Guillet et al. 2009) and therefore modify the surface available for the recombination of ionized species.

The gas-phase composition of the fluid particle results from the setup of an extensive chemical network. Throughout this work, the gas-phase elemental abundances and the depletion in grain cores are set to the standard values derived in the solar neighborhood (see Table 1 of Kristensen et al. 2023). Initially designed to study the formation of molecules in dense environments, the chemical network implemented in the online version of the Paris-Durham shock code includes nine different elements (H, He, C, N, O, S, Mg, Si, and Fe) and 141 species that interact with one another through an ensemble of ~2800 chemical reactions whose rates were updated in 2022. These include cosmic ray ionizations, photodissociation and photoionization processes, radiative and dissociative recombinations, collisional ionizations and dissociations, and ion-neutral and neutral-neutral reactions.

For the purposes of the present paper and the follow-up papers, this chemical network has been further extended to take into account the evolution of the abundances of multi-ionized species. The ionization stages of all the elements listed above with a formation enthalpy below 10⁸ K have been included, along with photoionization processes, electron impact collisional ionizations, radiative and dielectronic recombinations, and recombinations onto dust and PAHs. These modifications have increased the number of species treated in the code to 209. In practice, the photoionization rate of any ion i at a position z along the shock trajectory (see Fig. 1 of Godard et al. 2019) is computed from the integration of its frequency dependent ionization cross section, σ_i(ν), over the specific intensity of the local radiation field, I_ν, as $$k_{\gamma }^i(z)=\int \int {\sigma }_i(\nu )\frac{1}{h\nu }{I}_{\nu }(z)\mathrm{d}\nu \mathrm{d}\mathrm{\Omega },$$(1)

where $I_{\nu }=I_{\nu }^{\mathrm{e}\mathrm{x}\mathrm{t}}+I_{\nu }^{\mathrm{i}\mathrm{n}\mathrm{t}}$ includes the contribution of both the external ISRF (see Sec. 2.1) and the radiation field generated by the shock itself (see Sec. 2.4 and Appendix A). As illustrated in Appendix B, the photoionization cross sections are taken from Verner & Yakovlev (1995). The product of any photoionization process and the probability to produce Auger’s electrons are taken from Kaastra & Mewe (1993). The collisional ionization rates of all atoms and atomic ions via electron impact are computed using the database of Lennon et al. (1988) (Tables 2 and 3 and Eqs. (6) and (7) in their paper). The radiative recombination rates are calculated from the data and prescriptions given by Badnell (2006) (Table 1 and Eqs. (1) and (2) in their paper), Altun et al. (2007) (Table 2 in their paper), and Abdel-Naby et al. (2012) (Table 4 in their paper). The dielectronic recombination rates are derived from the database³ provided by Badnell et al. (2003). In the few cases where this database is incomplete, the dielectronic recombination rates are alternatively taken from Landini & Monsignori Fossi (1990) (Table 1 in their paper). The recombination rates on dust and PAHs are finally computed using the prescription of Draine & Sutin (1987).

2.3 Cooling processes

The cooling processes treated in the public version of the Paris-Durham shock code are described in Flower & Pineau des Forêts (2003), Lesaffre et al. (2013), and Godard et al. (2019). In particular, the code includes the excitation of the fine structure and metastable levels of all neutral and singly ionized atoms, the excitation of the rovibrational levels of several molecules (e.g., H₂, ¹²CO, and ¹³CO), and the energy lost through the ionization and dissociation of atoms and molecules by collisions with H, He, H⁺, H₂, and e⁻. To accurately follow the thermal evolution of a plasma shocked at high temperature (T > 10⁵ K), three additional processes have been implemented in the shock model: the continuum Bremsstrahlung emission, the loss of energy due to the ionization of ionized species by electron impact, and the collisional excitation of multi-ionized species.

The total Bremsstrahlung cooling rate per unit volume (in erg cm⁻³ s⁻¹) is modeled with Eq. (6.17) of Dopita & Sutherland (2003) $${\mathrm{\Lambda }}_{\mathrm{f}\mathrm{f}}=\frac{16}{3\sqrt{3}}{\left(\frac{(2\pi {)}^{3}k{T}_{\mathrm{e}}}{h^2{m}_{\mathrm{e}}^3}\right)}^{1/2}{\left(\frac{q_{\mathrm{e}}^2}{c}\right)}^3{n}_{\mathrm{e}}{\displaystyle \sum _i}{n}_i{Z}_i^2⟨g_{\mathrm{f}\mathrm{f}}\left({\gamma }^2\right)⟩,$$(2)

where k and h are the Boltzmann and the Planck constants, c is the speed of light, q_e, m_e, T_e, and n_e are the charge, mass, temperature, and density of the electrons, and n_i and Z_i are the density and the ionization stage of each ion i. The total freefree Gaunt factors, ⟨g_ff(γ²)⟩, are calculated using a logarithmic interpolation of Table 4 of van Hoof et al. (2014) who computed these factors over a wide range of values of the scaled quantity ${\gamma }^2=Z_i^2{E}_{\mathrm{R}\mathrm{y}}/k{T}_{\mathrm{e}}$, where E_Ry is the energy associated with one Rydberg (2.17987 × 10⁻¹¹ erg).

The cooling induced by electron impact ionization is treated in the same way than the thermal energy lost during any endothermic chemical reaction (Flower et al. 1985). The cooling rate following the ionization of an ion i writes $${\mathrm{\Lambda }}_{\mathrm{i}\mathrm{o}\mathrm{n}}^i=n_{\mathrm{e}}{n}_i{I}_i{k}_{\mathrm{i}\mathrm{o}\mathrm{n}}^i\left(T_{\mathrm{e}}\right),$$(3)

where I_i is the ionization potential of the ion and $k_{\text{ion }}^i(T_{\mathrm{e}})$ is the temperature dependent rate coefficient for collisional ionization computed by Lennon et al. (1988) and used in the chemical network (see Sec. 2.2).

Finally, the cooling induced by the excitation of the electronic levels of any atom in any ionization stage by collision with electrons is computed using CHIANTI⁴ (Dere et al. 1997, 2019; Del Zanna et al. 2021). Originally developed to study solar and stellar atmospheres, the CHIANTI database contains an up-to-date set of atomic data (energy levels, radiative transition probabilities, and collisional excitation rates) for a large number of ions. For the species treated in this work, the database includes 22 948 energy levels and 763 970 radiative transitions. As described in Sec. 2.4, the populations of all the levels and the emissivities of all the transitions are computed self-consistently in the postprocessing radiative transfer. Such a detailed treatment is, however, numerically prohibitive for the Paris-Durham shock code. To simplify, the time-dependent cooling is computed with precalculated functions. Assuming that each species is in its ground state, the cooling induced by any ion writes $${\mathrm{\Lambda }}_{\mathrm{e}\mathrm{x}\mathrm{c}}^i=n_{\mathrm{e}}{n}_i{\displaystyle \sum _j}\mathrm{\Delta }{E}_j{C}_j\left(T_{\mathrm{e}}\right),$$(4)

where ΔE_j is the energy difference between the level j and the ground state, C_j(T_e) are the temperature dependent excitation rate coefficient between the ground state and the level j provided by the CHIANTI database, and the sum is performed over all excited electronic levels. The sum in Eq. (4) is fitted, with an accuracy of 10% between 10 and 10⁸ K, with a sum of modified Arrhenius components using the curve fit optimization algorithm provided by the SciPy python library, and is used as a cooling function in the Paris-Durham shock code. The validity of this approach is checked, at each iteration, by comparing the cooling computed in the Paris-Durham shock code to the total emission of each species computed by the postprocessing radiative transfer tool.

2.4 Postprocessing of the radiative transfer

The temperature reached in high-velocity shocks, the productions of multi-ionized species and free electrons, and the collisional excitation of atoms and ions, lead to the emission of continuum and line radiations. By heating and ionizing the shocked and preshocked materials, the UV, EUV, and X-ray photons exert a negative feedback that acts against their own production. An accurate description of the propagation of these photons is therefore essential to compute the thermodynamical structure of the shock. The difficulty of the task lies in the treatment of photon trapping, that is the resonant absorption and reemission of line photons, which can only escape through frequency scattering in the line wings.

In this paper, the propagation of photons is computed in postprocessing with the following methodology. (1) All the emission lines contained in the CHIANTI database for any atom in any ionization stage are calculated with an exact radiative-transfer algorithm outlined in Appendix A. Following the principles of coupled escape probability (Elitzur & Asensio Ramos 2006), absorption and induced emission processes are handled as nonlocal couplings within the evolution equations of the level populations. The resulting system of equations is then solved, at equilibrium, taking into account all available levels, collisional rates, and radiative transitions. (2) The emission of continuum radiation is calculated including Bremsstrahlung, recombination, and two-photon processes. For the sake of simplicity and for the purposes of the present paper, dust emission is not included. (3) The line fluxes and the flux density of continuum radiation are then computed over the shock profile and the preshock medium taking into account the absorption due to photoionization and the extinction induced by dust and PAHs (see Appendix B).

The emissivity of free-free radiation (in erg cm⁻³ s⁻¹ Hz⁻¹ sr⁻¹) at frequency v is taken from Eq. (6.14) of Dopita & Sutherland (2003), $$\begin{array}{rl}{ϵ}_{\mathrm{f}\mathrm{f}}(\nu )=& \frac{16}{3\sqrt{3}}{\left(\frac{\pi }{2k{m}_{\mathrm{e}}^3{T}_{\mathrm{e}}}\right)}^{1/2}{\left(\frac{q_{\mathrm{e}}^2}{c}\right)}^3{n}_{\mathrm{e}}\exp (-\frac{h\nu }{k{T}_{\mathrm{e}}})\\ & \times \sum _i{n}_i{Z}_i^2⟨g_{\mathrm{f}\mathrm{f}}({\gamma }^2,\frac{h\nu }{k{T}_e})⟩,\end{array}$$(5)

where ⟨g_ff(γ², hν/kT_e)⟩ are the frequency dependent free-free Gaunt factors (taken from Table 3 of van Hoof et al. 2014) and the sum is performed over all ions i. The emissivity of the free-bound radiation induced by the recombination of an ion i onto a bound state b of the recombined ion i′ is modeled using Eq. (10.25) of Draine (2011) as $${ϵ}_{\mathrm{f}\mathrm{b}}^{i,\mathrm{b}}(\nu )=\frac{1}{2}\frac{g_{\mathrm{b}}}{g_i}\frac{h^4{\nu }^3}{{\left(2\pi m_{\mathrm{e}}k{T}_{\mathrm{e}}\right)}^{3/2}{c}^2}{n}_{\mathrm{e}}{n}_i{\sigma }_{i^{\mathrm{\prime }}}^{\mathrm{b}}(\nu )\exp (-\frac{h\nu -I_{\mathrm{b}}}{k{T}_{\mathrm{e}}}),$$(6)

where g_i is the degeneracy of the ground state of the ion i, and g_b, I_b, and ${\sigma }_{i^{\mathrm{\prime }}}^{\mathrm{b}}$ are the degeneracy, the ionization potential, and the photoionization cross section of the recombined ion i′ in the bound state b. The photoionization cross section of the ground state of the recombined ion is taken from Verner & Yakovlev (1995). For all the other states, the photoionization cross sections are derived from the prescription of Karzas & Latter (1961), setting the associated bound-free Gaunt factors to one. Finally, and as described in Appendix A, any transition that occurs through the emission of two photons is treated by the radiative-transfer algorithm as an optically thin line with a radiative decay rate, $A_{2\mathrm{p}\mathrm{h}}$, provided by the CHIANTI database. The emissivity of the two-photon emission process is then computed with the prescription of Shull & McKee (1979) as $${ϵ}_{2\mathrm{p}\mathrm{h}}(\nu )=\frac{h{\nu }_{2\mathrm{p}\mathrm{h}}}{4\pi }{A}_{2\mathrm{p}\mathrm{h}}{n}_{2\mathrm{p}\mathrm{h}}\frac{6\nu }{{\nu }_{2\mathrm{p}\mathrm{h}}^2}(1-\frac{\nu }{{\nu }_{2\mathrm{p}\mathrm{h}}}),$$(7)

where ν_2ph is the frequency associated with the transition and n_2ph is the density of the upper level obtained with the radiative-transfer algorithm.

As shown above and in Appendix A, we made the choice to address the propagation of photons in interstellar shocks with a detailed treatment that involves an exact and convoluted radiative-transfer algorithm. Such a detailed treatment is paramount to predict the fluxes and the line profiles of specific resonant transitions like, for instance, the Lyman series of atomic hydrogen (see Fig. A.2). It should be noted, however, that all these intricacies have a weak influence on the couplings between the self-generated radiation field and the thermochemical state of the gas. Indeed, we find that treating all lines as optically thin transitions provides similar results regarding the thermochemical structure of the shock and the extent and physical state of the radiative precursor. This is not due to the velocity broadening or velocity gradients induced by shocks. It comes from the fact that if a specific line is optically thick, either the photons escape in the wings of the line, or the associated upper level cascade through the emission of other lines which, in total, produce a similar contribution to the ionizing and dissociating radiation field.

3 Thermochemical structure of WNM shocks

3.1 Standard model

The parameters of the standard model are given in Table 1. As a prototypical case, we consider an interstellar shock at a velocity of 200 km s⁻¹ propagating in a medium with an initial density $n_{\mathrm{H}}^0$ = 1 cm⁻³ and illuminated by the standard ISRF (G₀ = 1). The total H₂ cosmic ray ionization rate ${\zeta }_{{\mathrm{H}}_2}$ is set to 3 × 10⁻¹⁶ s⁻¹, as deduced from the observations of ${\mathrm{H}}_3^{+}$, OH⁺, H₂O⁺, and ArH⁺ in the local diffuse ISM (McCall et al. 2003; Indriolo et al. 2007, 2015; Neufeld & Wolfire 2017). Finally, the transverse magnetic field strength B₀ is set to 1 μG, in agreement with the magnetic field plateau obtained from Zeeman measurements for gas at low total column density (Crutcher et al. 2010) and with the Faraday rotation measurements along diffuse lines of sight (Frick et al. 2001). All the parameters are set to the typical conditions of the WNM observed in the local ISM, except for the shock velocity. Indeed, shocks at V_S = 200 km s⁻¹ are not particularly frequent in the diffuse local ISM. This choice is mostly made here to test the validity of the model and the interplay of the physical processes implemented in the code. With these parameters, the standard model is a single fluid J-type shock with no decoupling between the ionized and the neutral flows.

3.2 Thermal profiles and electronic fraction

The outcome of the iterative procedure is illustrated in Fig. 3, which displays the profiles of the kinetic temperature (top panels), the fractional abundances of H and H⁺ (middle panels), and the fractional abundances of He, He⁺, and He⁺⁺ (bottom panels) obtained with the standard model in the preshock (left panels) and postshock medium (right panels). The left panels of this figure are similar to Fig. 2 of Sutherland & Dopita (2017), which shows the ionization and thermal structures of the preshock for a shock at 150 km s⁻¹ with no external UV and optical radiation field. As found by Sutherland & Dopita (2017), we find that about 20 iterations are required to reach convergence on the radiative transfer and therefore on the thermochemical profiles of the preshock gas.

The magnetohydrodynamics and thermochemical equations are solved in the Paris-Durham shock code using a forward integration technique. Figure 3 therefore shows the out-of-equilibrium evolution of a fluid particle traveling from the preshock gas to the postshock medium. The initial conditions of the gas are typical of the WNM, with a kinetic temperature of ~8000 K and an ionization fraction of ~0.01 (see Fig. 2). As it encounters the radiative precursor, the gas becomes almost fully ionized and is heated to ~13 400 K. The flux of ionizing photons generated by a shock at 200 km s⁻¹ is so strong that the gas is found to be at thermochemical equilibrium along its entire trajectory through the preshock. This is in line with the results of Sutherland & Dopita (2017), who show that the preshock gas can be considered at thermochemical equilibrium for shock velocities V_S ⩾ 120–150 km s⁻¹.

As described in previous studies, the radiative precursor has mainly two feedbacks on the shock structure. First, the increase of the ionization state of the gas in the radiative precursor implies that the mean mass of particles drops from ~1.27 to ~0.6 a.m.u. As dictated by the jump conditions of adiabatic magnetohydrodynamic shocks (Roberge & Draine 1990), the postshock maximal temperature is thus reduced by a factor of ~2.1 (see the initial and final trajectories in Fig. 3). Second, the gradient of temperature induced in the preshocked gas exerts a force that slows down and compresses the fluid. The gas therefore enters the shock with a reduced effective velocity. This effect is mitigated here because the initial ram pressure is about three orders of magnitude larger than the maximum thermal pressure in the preshock.

With all these effects self-consistently computed, a shock propagating at 200 km s⁻¹ in the WNM is found to heat the gas to a postshock temperature of ~6 × 10⁵ K. The evolution of the postshock medium is then driven by the cooling mechanisms and the recombination processes. At these temperatures, we find that the cooling time is ~4 × 10⁴ yr, which corresponds to a cooling length of 2 pc (see Fig. 3). Both these results are in fair agreement (within a factor of two) with Eq. (36.33) of Draine (2011) who estimated the cooling time of high-velocity shocks with a simplified cooling function for hot plasmas (T > 10⁵ K).

3.3 Chemical profiles

Figure 3 (right panel) shows that the large temperatures produced by the standard model imply that both hydrogen and helium are fully ionized in the postshock medium. Interestingly, in both cases, an ionization equilibrium is reached between collisional ionization and radiative and dielectronic recombinations that lasts until the gas cools down to temperatures below 10⁵ K. As a comparison, and for future analyses of observational tracers, Fig. 4 displays the abundances of the different ionization stages of carbon, nitrogen, and oxygen obtained at the final iteration of the standard model. In contrast to hydrogen and helium, none of the elements shown in Fig. 4 (right panel) is fully ionized in the standard shock model and none of them reaches ionization equilibrium. This particular result highlights the need to follow the time dependent evolution of the chemistry even in stationary shock structures.

The analysis of the preshock (Fig. 4, left panel) shows that the soft X-ray radiation field generated by the standard model is strong enough to doubly ionize carbon, nitrogen, and oxygen in the radiative precursor. This feature is particularly important for the interpretation of observations of atoms obtained in the diffuse local ISM.

Fig. 3 Convergence of the standard model on the shock thermochemical profiles. Profiles of the kinetic temperature (top panels), the fractional abundances of H and H+ (middle panels), and the fractional abundances of He, He+, and He++ (bottom panels) as functions of the distance in the preshock (left panels) and in the shock (right panels). The black curves show the intermediary profiles computed during the iterative procedure and the colored curves the profiles obtained in the final converged iteration.

Fig. 4 Production of multi-ionized species in the standard model. Abundances profiles of carbon-bearing species (top panels), nitrogen-bearing species (middle panels), and oxygen-bearing species (bottom panels) as functions of the distance in the preshock (left panels) and in the shock (right panels). The high ionization stages of carbon, nitrogen, and oxygen (C6+, N6+, N7+, O7+, and O8+) are not included because their relative abundances are all below the lowest value shown in this figure.

Fig. 5 Trajectory of a shocked fluid particle from the ambient medium to the postshock region obtained in the standard model. The trajectory is displayed in a particle density versus thermal pressure diagram (top panel) and in a proton density versus thermal pressure diagram (bottom panel) to highlight the phase transition from the WNM to the CNM induced by a shock at 200 km s−1. A fluid particle initially evolves through the radiative precursor (pink curve). As it crosses a shock front, it undergoes three successive regimes: an adiabatic heating (red curve) followed by a quasi-isobaric cooling (green curve), and finally a quasi-isochoric cooling (blue curve) as the magnetic pressure becomes dominant in the postshock gas. The black curve indicates the thermal equilibrium state of the diffuse gas obtained for G0 = 1 (see Fig. 1). Light gray lines are isothermal contours from 1 to 107 K (from bottom right to top left). We note that these isocontours are evenly spaced in the top panel but not in the bottom panel. Tics on the trajectories indicate the time, in log scale, starting from the adiabatic phase, ranging from 10−1 yr to 106 yr.

3.4 Phase transition

The impact of an interstellar shock on the thermodynamical state of the diffuse matter is shown in the top panel of Fig. 5, which displays the trajectory of the standard model in a particle density versus thermal pressure diagram. To bring out the physical properties of shocks, leaving aside the variation of the number of particles induced by the ionization in the preshock and the recombination in the postshock, we also display this trajectory in a proton density versus thermal pressure diagram in the bottom panel of Fig. 5. While this figure is drawn for a particular set of parameters, it depicts a very general modus operandi of stationary interstellar shocks. A fluid particle crossing a shock front basically undergoes three successive regimes.

First, the gas is heated adiabatically (red curve in Fig. 5). For strong magnetohydrodynamic shocks, that is shocks where the sonic and the Alfvénic Mach numbers are much higher than one, the Rankine-Hugoniot relations imply that the mass density ρ, velocity, V, and temperature, T, of the gas after the adiabatic jump are (Roberge & Draine 1990) $$\rho =\frac{\gamma +1}{\gamma -1}{\rho }_0\text{, }$$(8) $$V=\frac{\gamma -1}{\gamma +1}{V}_{\mathrm{S}}\text{, }$$(9)

and $$T=\frac{2(\gamma -1)}{(\gamma +1{)}^2}\frac{\mu V_{\mathrm{S}}^2}{k},$$(10)

where ρ₀ is the mass density before the jump, μ is the mean mass of the particles, and γ is the adiabatic index. In such conditions, the thermal pressure after the jump writes $$P=\frac{2}{(\gamma +1)}{\rho }_0{V}_{\mathrm{S}}^2\text{, }$$(11)

hence $$P=5.1\times {10}^6\mathrm{K}{\mathrm{c}\mathrm{m}}^{-3}\left(\frac{n_{\mathrm{H}}^0}{1{\mathrm{c}\mathrm{m}}^{-3}}\right){\left(\frac{V_{\mathrm{S}}}{200\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right)}^2.$$(12)

After the jump, and if the magnetic pressure is small compared to the termal pressure, the gas cools down isobarically (green curve in Fig. 5). This isobaric equation of state comes from the fact that the shock is at steady-state and that the postshock thermal pressure must therefore be at equilibrium with the preshock ram pressure. As the gas cools down, the postshock density increases. Because the magnetic field is frozen in the ionized gas, the magnetic field strength increases accordingly.

The final point of the trajectory is obtained when the postshock gas reaches thermal equilibrium. If, on the one hand, the magnetic pressure remains negligible over the postshock, the final density of the gas is given by the first intersection between the postshock thermal pressure (Eq. (12)) and the thermal equilibrium curve (black curve in Fig. 5). If, on the other hand, the postshock magnetic pressure is no longer negligible compared to the postshock thermal pressure, the gas progressively shifts from an isobaric to an isochoric evolution (blue curve in Fig. 5). In this case, the final density is given by $$n_{\mathrm{H}}^{\mathrm{f}}=146{\mathrm{c}\mathrm{m}}^{-3}{\left(\frac{n_{\mathrm{H}}^0}{1{\mathrm{c}\mathrm{m}}^{-3}}\right)}^{3/2}\left(\frac{V_{\mathrm{S}}}{200\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right){\left(\frac{B_0}{1\mathrm{\mu }\mathrm{G}}\right)}^{-1}.$$(13)

As already shown though theoretical works and numerical simulations (e.g., Hennebelle & Pérault 1999; Falle et al. 2020), all these properties imply that interstellar shocks can induce a phase transition between the WNM and the CNM. This is illustrated in Fig. 5, which shows that the standard model naturally converts the initial warm neutral material into cold neutral material. From the above considerations, it is clear that phase transition occurs only if the thermal pressure is larger than the maximum pressure of the WNM and if the maximum density allowed by magnetic compression is larger than the minimum density of the CNM. For the gas-phase elemental abundances adopted in this work (see Sec. 2.2) and for 0.1 ⩽ G₀ ⩽ 10, we estimate that the maximum pressure of the WNM scales as $$P_{\mathrm{W}\mathrm{N}\mathrm{M}}^{max}\sim 8\times {10}^3\mathrm{K}{\mathrm{c}\mathrm{m}}^{-3}{G}_0{}^{0.5},$$(14)

and the minimum density of the CNM as $$n_{\mathrm{H},\mathrm{C}\mathrm{N}\mathrm{M}}^{min}\sim 8{\mathrm{c}\mathrm{m}}^{-3}{G}_0{}^{0.5}.$$(15)

These scaling relations combined with Eqs. (12) and (13) imply that interstellar shocks induce phase transition if $$V_{\mathrm{S}}⩾8\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}{\left(\frac{n_{\mathrm{H}}^0}{1{\mathrm{c}\mathrm{m}}^{-3}}\right)}^{-1/2}{G}_0{}^{0.25},$$(16)

and if $$V_{\mathrm{S}}⩾11\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}{\left(\frac{n_{\mathrm{H}}^0}{1{\mathrm{c}\mathrm{m}}^{-3}}\right)}^{-3/2}\left(\frac{B_0}{1\text{ μ}\mathrm{G}}\right)G_0{}^{0.5}.$$(17)

The validity of these conditions is discussed in Sec. 4.

3.5 Emitted spectrum

As described in the previous sections, the radiation field generated by the shock has a fundamental impact on the thermochemical structure of the gas. In addition, since most of the initial kinetic energy is reprocessed into line and continuum emission (Allen et al. 2008), this radiation field is a unique signature of the shock itself. To illustrate the conversion process, we display in Fig. 6 the energy flux, $\nu F_{\nu }^{\text{int }}$, emitted toward the preshock, at the origin of the distances (defined as z = 0, see Fig. 3), by the standard model. Figure 6 summarizes the strengths and the relative contributions of the different emission processes, including the Bremsstrahlung, bound-free, and two-photon continuum emissions, and the emissions of resonant and non-resonant lines.

Figure 6 displays the fluxes of all the 763 970 lines contained in the CHIANTI database (see Sec. 2.3), regardless of their strength. For the standard model, we find that only ~2200 lines have an energy flux above 3% of the continuum level. These dominant lines carry more than 90% of the initial flux of mechanical energy. About 56% of the line flux is emitted in the EUV and X-ray domain (λ ⩽ 911 Å), 37% in the UV range (911 < λ ⩽ 2400 Å), and 7% in the optical and the infrared (λ > 2400 Å). All the elements display clear line features that extend from a few tens of Angströms to a few hundreds of micrometers. Yet, the energy budget shows that most of the photon emission is carried out by electronic lines of He⁺, O³⁺, O⁴⁺, O⁵⁺, and H, and to a lesser extent by the electronic lines of C⁺, C²⁺, C³⁺, N²⁺, N³⁺, N⁴⁺, S⁺, and S²⁺. Seemingly prominent in Fig. 6, both the Bremsstrahlung and two-photon continuum emissions only carry a few percent of the initial flux of mechanical energy.

4 Grid of models of WNM shocks

In the previous section, the Paris-Durham shock code was applied with the standard set of parameters in order to validate the recent updates and highlight the thermochemical state of the gas and the main processes at play in shocks propagating in the WNM. To test the robustness of the code and broaden its applications, we present here an exploration performed over a grid of models that encompass the typical physical conditions of the WNM expected in the solar neighborhood. The grid of models contains 475 runs with an initial preshock density varying between 0.1 and 2 cm⁻³, a transverse magnetic field varying between 0.1 and 10 μG, and a shock velocity varying between 10 and 500 km s⁻¹.

Fig. 6 Energy flux, computed as $\nu F_{\nu }^{\text{int }}$, emitted toward the preshock at z = 0 in the standard model as a function of the wavelength of the emitted photon. The black curve corresponds to the total continuum emission (including Bremsstrahlung, bound-free, and two-photon processes) and the colored spikes to line emissions. Spikes are colored according to the element responsible for the line emission regardless of its ionization state. The top panel displays the full wavelength range while the bottom panel is a zoom onto the range between 10 and 2000 Å. The black vertical line indicates the ionization threshold of hydrogen.

4.1 Critical velocities and expected shock types

The nature of an interstellar shock is the result of a subtle interplay between the strength of the magnetic field, the ionization degree, and the cooling rate (e.g., Chernoff 1987; Roberge & Draine 1990; Godard et al. 2019). The prototypical model presented in the previous section is a single fluid interstellar shock in which the ionized and the neutral fluids are fully coupled with each other. This kind of shock is referred to as a J-type shock because of the discontinuity in the temperature, velocity, and density profiles, which abruptly change at the scale of the particle mean free path (see Fig. 3 and Sec. 5.1). However, different physical conditions of the preshocked gas and different shock velocities may lead to a decoupling between the ionized and the neutral fluids and to the formation of a magnetic precursor. This process gives rise to other kinds of shock structures called CJ-type and C-type shocks depending on whether the gas abruptly becomes subsonic at some point along the trajectory or not. The types of shocks obtained over the grid of models explored here are shown in Fig. 7.

The division between the domains of existence of different types of shocks can be understood by comparing the shock velocity with the characteristic wave speeds of the interstellar gas. The magnetosonic speeds of a medium with a magnetic field strength B are $$c_{\mathrm{n}\mathrm{m}\mathrm{s}}={(c_n^2+B^2/4\pi {\rho }_n)}^{1/2},$$(18)

and $$c_{\text{ims }}={(c_i^2+B^2/4\pi {\rho }_i)}^{1/2}\text{, }$$(19)

where c_n and c_i are the sound speeds in the neutral and the ionized fluids, and ρ_n and ρ_i are the mass densities of the neutrals and the ions. In the limit of strong coupling between the ions and the neutrals, the minimum speed required for a shock to propagate is given by the neutrals magnetosonic speed (Eq. (18)). If the shock velocity is larger than the ions magnetosonic speed (Eq. (19)), the shock is a J-type shock. Below this limit, the shock is either a CJ-type or a C-type shock. Because the cooling rate of the WNM is inefficient to reduce the temperature gradients induced by the compression and the ion-neutral drift, the limit between these two kinds of shocks can be estimated as the point where the maximum postshock thermal pressure (given by Eq. (12)) is equal to the maximum magnetic pressure allowed by an adiabatic compression. This limit is shown as a dashed line in Fig. 7 and is found to be in excellent agreement with the results of the Paris-Durham shock code, which shows that shocks below (respectively, above) this value are C-type (respectively, CJ-type).

Fig. 7 Critical velocities of interstellar shocks for $n_{\mathrm{H}}^0$ = 0.1 cm−3 (top panel) and 1 cm−3 (bottom panel). The solid and dash-dotted lines show the magnetosonic speed computed in the ionized fluid (Eq. (19)) and the neutral fluid (Eq. (18)), respectively. The colored areas highlight the range of velocities over which stationary shocks are C-type (red), CJ-type (blue), and J-type (green) shocks. The transition between C-type and CJ-type shocks (dashed line) is estimated as the point where the maximum postshock thermal pressure is equal to the maximum magnetic pressure allowed by an adiabatic compression. The gray points indicate the values of the shock velocity and transverse magnetic field strength explored in the grid of models. The red dot highlights the position of the standard model.

4.2 Phase trajectories

The trajectories followed by all the models explored in the grid are shown in Fig. 8 in a proton density versus thermal pressure diagram. In many ways, this simple figure is the main outcome of this work. It provides, at a glance, an almost complete summary of all the models explored. It also demonstrates the robustness of the code for velocities between 10 and 500 km s⁻¹ and the coherence of its predictions compared to the theoretical and analytical descriptions given in the previous sections.

For standard values of the magnetic field expected in the WNM (B₀ ≲ 1 μG), most of the models are found to be J-type shocks, except at low velocities (V_S ≲ 20 km s⁻¹) where the code converges toward stationary CJ-type solutions. This result is in agreement with the expectations given in Fig. 7 based on the sole consideration of characteristic wave speeds. In this parameter domain, the predicted trajectories are found to follow the description given in Sec. 3.4. In particular, interstellar shocks display successive adiabatic, isobaric, and isochoric evolutions with a maximum thermal pressure and a final compression that match the expected values given by Eqs. (12) and (13). It follows that interstellar shocks are found to induce phase transition between the WNM and the CNM as long as the shock velocity verifies the conditions given in Eqs. (16) and (17). These criteria are fulfilled by a vast majority of the shocks, except at low preshock density and low velocity where the postshock thermal pressure is below the maximum pressure of the WNM, and at low preshock density and large values of B₀ where the maximum density of the final state is below the minimum density of the CNM.

For larger values of the strength of the magnetic field (B₀ ≳ 3 μG), a substantial fraction of the models are found to be C-type shocks, which is in agreement with the expectation given in Fig. 7. The strength of the preshock magnetic field has mostly two impacts on the shock structure. When the preshock magnetic pressure is comparable to the ram pressure, the magnetic field is strong enough to reduce or even remove the rise of the gas temperature. The initial growth of thermal pressure thus occurs with a lower slope than that predicted by the adiabatic hydrodynamic jump conditions. When the Alfvénic Mach number is ≲3, the magnetic pressure in the postshock becomes strong enough to entirely suppress the isobaric cooling stage. It follows that the gas rapidly reaches an isochoric cooling evolution toward thermal equilibrium. Although many shocks in this domain of parameter occur through ion-neutral decoupling, or with reduced postshock thermal pressure, the criteria established to obtain phase transition are still valid. It is so because, for large values of the magnetic field, the sole criterion that dictates a possible transition between the WNM and the CNM is Eq. (17), which is simply derived from the conversion of the initial ram pressure into magnetic pressure.

Interestingly, Fig. 8 shows that the final states of several shocks explored in the grid lie on the unstable branch of the thermal equilibrium state (black curve in Fig. 8). Because the shocks are described as stationary plane-parallel structures, the model predicts that the postshock gas remains on this unstable branch indefinitely. At first sight, these solutions could appear to be unphysical (Falle et al. 2020). However, when such shocks occur in nature, the gas effectively follows the trajectories displayed in Fig. 8. The unstable nature of the final state just implies that the postshock gas then naturally splits up into a warm and a cold neutral components with mass fractions that depend on the thermal pressure of the surrounding gas. This final evolution doesn’t taint the predictions of this work regarding the thermochemical evolution and the emission properties of the shock.

Conversely, Fig. 8 shows that all the models at $n_{\mathrm{H}}^0$ = 2 cm⁻³ start from an unstable thermal state. This is due to the fact that the preshock conditions are calculated assuming thermochemical equilibrium regardless of possible thermodynamical instabilities. These conditions might seem unrealistic if we consider that the diffuse medium is composed solely of WNM and CNM. This is not the case. Indeed, observations of the 21 cm line of atomic hydrogen show that the unstable phase, sometimes called the Lukewarm Neutral Medium (LNM, Marchal et al. 2019; Bellomi et al. 2020), contains a substantial fraction (~20%) of the mass of the diffuse ISM (e.g., Murray et al. 2018; Marchal et al. 2019). Although its volume filling factor is small compared to that of the WNM, the LNM is still an important component of the diffuse matter where shocks could develop. The unstable initial conditions obtained at $n_{\mathrm{H}}^0$ = 2 cm⁻³ can therefore be seen as an idealized modeling of shocks propagating in the LNM.

Fig. 8 Trajectories followed by interstellar shocks in a proton density versus thermal pressure diagram (see Fig. 5). Each panel corresponds to different values of the initial density $n_{\mathrm{H}}^0$ (top labels) and transverse magnetic field B0 (right labels). The different curves correspond to the following values of the shock velocity ordered from top to bottom trajectories: 500, 450, 400, 350, 300, 250, 200, 150, 120, 100, 90, 80, 70, 60, 50, 40, 30, 20, and 10 km s−1. Colored trajectories correspond to J-type shocks, gray trajectories to CJ-type shocks, and black trajectories to C-type shocks. The black curve common to all panels indicates the thermal equilibrium state of the diffuse gas as a function of the proton density obtained for G0 = 1 (see Fig. 1). Light gray lines are isothermal contours from 1 to 107 K (from bottom right to top left). We note that the trajectories of low-velocity shocks are not always visible as they only slightly perturb the gas from its original state.

Fig. 9 Fractional abundances of H and H+ (top panel), He, He+, and He++ (middle panel), and C+, C2+, C3+, C4+, C5+, and C6+ (bottom panel) in the gas entering the shock as functions of the shock velocity. Except for the velocity, all other parameters are set to their standard values (see Table 1).

4.3 Preshock ionization state

The thermochemical profiles of intermediate- and high-velocity shocks and the potential decoupling between the ionized and the neutral flows are driven, in part, by the impact of the radiation field generated by the shock on the thermochemical state of the preshock gas. To illustrate this aspect, we display in Fig. 9 the chemical composition of the radiative precursor (at z = 0) in hydrogen-, helium-, and carbon-bearing compounds as a function of the velocity of the shock for $n_{\mathrm{H}}^0$ = 1 cm⁻³ and B₀ = 1 μG. This figure shows that the ionization state of the preshock gas displays several transitions depending on the shock velocity, and therefore on the strength and hardness of the emitted radiation field.

At V_S ~ 100 km s⁻¹, the radiative precursor is found to abruptly switch to a medium dominated by H⁺, He⁺ and C²⁺. The reason for this sharp transition is that a shock at 100 km s⁻¹ not only collisionally ionizes neutral helium but is also strong enough to collisionally excite the 300 Å lines of He⁺. These lines, which dominate the EUV spectrum (see Fig. 10), correspond to photon energies of ~40 eV much larger than the ionization potential of H(~13.6 eV), He(~24.6 eV), and C⁺ (~24.4 eV), yet still smaller than the ionization potential of C²⁺(~47.9 eV). This behavior is in line with the Figs. 7, 8, and 9 of Allen et al. (2008) and in excellent agreement with the predictions of Hollenbach & McKee (1989) who show that the high-velocity shocks fully ionize the radiative precursor for shock velocities V_S ≳ 120 km s⁻¹ (Fig. 11 in their paper).

As the shock velocity increases, the radiative precursor displays progressive transitions toward higher ionization stages. For instance, we find that the preshock gas is mainly composed of He²⁺ above ~200 km s⁻¹, and mainly composed of C³⁺ and C⁴⁺ above ~300 km s⁻¹. This smooth transitions toward higher ionization stages comes from the fact that shocks at larger and larger velocities progressively excite a forest of lines of atomic ions at larger and larger energies (see Fig. 10), which contribute together to the ionization of the different elements.

4.4 Reprocessing of kinetic energy

As an illustration of the reprocessing of the initial mechanical energy, we display in Fig. 10 the spectra emitted by six shocks with velocities of 10, 20, 50, 100, 200, and 400 km s⁻¹, and in Fig. 11 the fraction of the radiation (including lines and continuum) emitted over different domains of wavelength. For all the models displayed here, we find that most of the input flux of mechanical energy is reprocessed into line emission. This result extends almost over the entire grid of models, except for low Alfvénic Mach number where most of the energy goes in the compression of the initial magnetic field. Shocks at low velocity (V_S ≲ 30 km s⁻¹) mostly emit photons in the infrared and submillimeter domains. Shocks at intermediate velocities (30 < V_S < 100 km s⁻¹) emit most of their radiation in the UV and EUV range, while shocks at high velocity (V_S ⩾ 100 km s⁻¹) preferentially emit in the EUV and X-ray domains.

At low and intermediate velocities, the amount of ionizing radiation (>13.6 eV) is found to be almost entirely carried by a few lines of He and He⁺ (red spikes in Fig. 10). This feature is responsible for the abrupt jump in the chemical state of the preshock medium obtained at V_S ~ 100 km s⁻¹ (see Fig. 9). Conversely, the ionizing radiation emitted by shocks at high velocity results from a combination of electronic lines of helium, carbon, oxygen, nitrogen, and sulfur in various stages of ionization. The gas-phase elemental abundances adopted here imply that the cooling is dominated by the different ionization stages of oxygen. However, Fig. 10 suggests that this would not necessarily be the case for a different set of elemental abundances.

The continuum flux is found to be largely dominated by two-photon emission of hydrogen and helium as long as V_S ≲ 200 km s⁻¹. Above this value, most of the continuum photons originate from Bremsstrahlung emission. Even though the continuum radiation becomes more prominent at larger velocities, it remains subdominant compared to line cooling over the entire grid of models. Given the dependences of the cooling rates induced by collisional excitation and the free-free radiation on the gas temperature, we estimate that continuum radiation should dominate over line emission for shock velocities V_s ≳ 1000 km s⁻¹. This domain of parameter, which dangerously starts to fall in the realm of relativistic shocks, is far outside the scope of the present study.

One last remarkable aspect of Fig. 10 is the behavior of the infrared and submillimeter parts of the spectra as a function of the shock velocity. Indeed, the fluxes of the dominant infrared and submillimeter lines are modified from 10 to 400 km s⁻¹, but only slightly compared to the rise of the flux of mechanical energy (6.4 × 10⁴ over this range of velocities). This feature highlights the fact that these lines are emitted when the gas cools from ~10⁴ to ~70 K. This evolution is common to all the shocks displayed here (see Fig. 8). The resulting infrared and submillimeter spectrum therefore weakly depends on the shock velocity.

Fig. 10 Energy flux, computed as $\nu F_{\nu }^{\text{int }}$, emitted toward the preshock at z = 0 as a function of the wavelength of the emitted photon. The black curve corresponds to the total continuum emission (including Bremsstrahlung, bound-free, and two-photon processes) and the colored spikes to line emissions. Spikes are colored according to the element responsible for the line emission regardless of its ionization state. Predictions are shown for shocks at different velocities, ranging from 10 km s−1 (top left panel) to 400 km s−1 (bottom right panel). All other parameters are set to their standard values (see Table 1). The dashed gray vertical line indicates the ionization threshold of hydrogen.

Fig. 11 Fraction of the radiation emitted by shocks in the X-ray and EUV domains (red curve), in the ultraviolet range (green curve), and in the optical and the infrared domains (blue curve) as functions of the shock velocity for $n_{\mathrm{H}}^0$ = 1 cm−3 and B0 = 1 μG.

4.5 Induced photo-ionization rates

Figure 12 displays the strengths of the UV radiation field and the ionization rate of hydrogen resulting from the EUV and X-ray photons generated by the shocks for different velocities and preshock densities, and for B₀ = 1 μG. As already suggested by Figs. 7 and 8, this result is found to be independent of the strength of the transverse magnetic field as long as B₀ ⩽ 1 μG, which correspond to the standard values expected in the WNM.

As intermediate-velocity shocks mostly emit UV photons, the strength of the UV radiation field generated by shocks is found to be proportional to the input flux of mechanical energy and to scale as $V_{\mathrm{S}}^3$ as long as V_S < 100 km s⁻¹. Above this value, the self-generated UV radiation field still increases as a function of the shock velocity but with a shallower slope ($\propto V_{\mathrm{S}}^2$), as most of the energy is radiated away in the X-ray domain. We find that the equivalent G₀ is below the fiducial value associated with the external UV radiation field (G₀ = 1) as long as V_S < 200 km s⁻¹ ${\left(0.5{\mathrm{c}\mathrm{m}}^{-3}/n_{\mathrm{H}}^0\right)}^{1/3}$. This indicates that only high-velocity shocks propagating in dense environments may have a substantial impact on the UV irradiation of the surrounding interstellar gas.

Intermediate-velocity shocks are found to generate a photoionization rate of hydrogen that scales as $V_{\mathrm{S}}^6$ for V_S < 100 km s⁻¹. This rate suddenly jumps at V_S ~ 100 km s⁻¹ because of the excitation of the electronic levels of He⁺ and reaches a plateau for higher velocities, ${\zeta }_{\mathrm{H}}={10}^{-10}{\mathrm{s}}^{-1}\left(n_{\mathrm{H}}^0/1{\mathrm{c}\mathrm{m}}^{-3}\right)$. This plateau comes from the fact that high-velocity shocks radiate the input mechanical energy by hardening the X-ray radiation field rather than increasing the output flux of EUV photons (see Fig. 10). Since the photoionization cross-section of hydrogen quickly drops as the photon energy increases (σ_H(ν) ∝ ν⁻³, see Fig. B.1), the hardening of the X-ray radiation field has no impact on the ionization rate. All these behaviors only apply to atomic hydrogen. Ionization rates of other species in different ionization stages display different scalings as functions of the velocity of the shock.

The spatial extensions of the ionization rates of the main constituents of the WNM predicted by the standard model are shown in Fig. 13. While the induced ionization of atomic hydrogen roughly extends over the size of the radiative precursor (~10 pc for the standard model), the ionization of other elements extends over much larger distances. This is due to the fact that these species are effectively ionized by photons at larger energies (see Fig. B.1), which propagate over higher distances in the preshock (see Fig. B.2).

Fig. 12 Strength of the self-generated UV radiation field (top panel) and ionization rate of hydrogen induced by the self-generated EUV and X-ray radiation field (bottom panel) in the radiative precursor, at z = 0. Both quantities are shown as functions of the shock velocity for different values of the preshock density and for B0 = 1 μG. The strength of the UV radiation field is calculated from the integral of $I_{\nu }^{\text{int }}$, over solid angle and frequency between 911 and 2400 Å, normalized to the value obtained by Mathis (Mathis et al. 1983). The strength of the external UV radiation field used in this work (G0 = 1) is shown as a gray line on the top panel. The cosmic ray ionization rate of atomic hydrogen used in this work $({\zeta }_{\mathrm{H}}={\zeta }_{{\mathrm{H}}_2}/2.3=1.3\times {10}^{-16}{\mathrm{s}}^{-1})$ is shown as a gray line on the bottom panel.

Fig. 13 Ionization rates of H, He, C+, O, N, S+, Si+, and Fe+ induced by the standard shock in the preshock gas as functions of the preshock distance. These species are selected because they are the main carriers of the elements included in the model in the WNM.

5 Discussion

5.1 Viscosity and jump conditions

The Paris-Durham shock code has the particularity of treating the viscous stresses in the evolution equation of the gas velocity for J-type shocks (Lesaffre et al. 2013). In practice, this is done by setting a constant viscous length $\lambda =1/\sigma n_{\mathrm{H}}^0$, where σ = 3 × 10⁻¹⁵ cm² is the cross section of elastic collisions between H₂ and H₂ computed by Monchick & Schaefer (1980). This prescription is at the root of the viscous jump shown, for instance, in Fig. 3 which has a typical size of 10⁻⁴ pc for $n_{\mathrm{H}}^0$ = 1 cm⁻³. While the resulting evolutions of the thermodynamical state of the gas are found to be in perfect agreement with the jump conditions expected for adiabatic magnetohydrodynamic shocks, it is important to note, however, that this viscous jump is somehow artificial for the models explored in this paper.

First, the viscous length should be set from the cross section of elastic collision between hydrogen atoms and not molecular hydrogen. In addition, this viscous length is relevant only if the radiative precursor is not fully ionized. Indeed, when the radiative precursor is fully ionized, the adiabatic jump of interstellar shocks is known to be mediated by long range interactions and plasmoid instabilities leading to collisionless shock waves.

Since the adiabatic jump conditions are always fulfilled and have no impact on the subsequent cooling trajectory of the gas, regardless of the microphysics mechanism responsible for the jump, we deliberately decided to keep the prescription of the Paris-Durham shock code as it is. The only impact is that the length of the jump is somehow larger than it should be.

5.2 Approximation of equilibrium for the level populations

The postprocessing radiative-transfer algorithm (see Appendix A) and the prescription adopted for the cooling induced by collisional excitation of ions (see Sec. 2.3) are built on the hypothesis that the populations of the levels of atoms and atomic ions are at equilibrium. To test this approximation, we developed two separated versions of the Paris-Durham shock code, one which solves the time dependent evolution of the populations of the first ten electronic levels of atomic hydrogen, and the other one the time dependent evolution of the three fine structure levels of atomic carbon. In both cases, the evolution of the populations of all the levels considered were found to be identical to those computed at equilibrium. It is so because the typical deexcitation timescales of electronic levels and of fine structure levels are orders of magnitude smaller than the cooling timescales. This statement also holds for levels that relax through the emission of a resonant photon which, because of photon trapping, have a local effective deexcitation coefficient several orders of magnitude below the spontaneous decay rate.

5.3 Absorption of the external radiation field

We consider in this work that the shocked gas is irradiated by two radiation fields: the radiation field emitted by the shock itself, which extends at all wavelengths, and an external UV radiation field, which extends between 911 Å and 2400 Å. As described in Appendix A.5, the internal radiation field is calculated taking into account all absorption processes including dust and PAHs. As stated in Sec. 2.1, the external radiation field is considered constant throughout the entire shock trajectory.

It is important to note that the Paris-Durham shock code is capable of treating the absorption of the external radiation field by dust and PAHs and the self-shielding of H₂ (Godard et al. 2019). We choose, however, to keep this field constant here to simplify the problem and to focus on the sole propagation of the internal radiation field. This is equivalent to assume that the shock is surrounded by external UV radiation sources, which is probably the case in the diffuse ISM. The validity of this approximation depends on the typical visual extinction of the shocks. The 1D visual extinction computed across shocks over the entire grid of models explored in this work is found to range between a few 10⁻³ and 0.3. These values indicate that dust absorption can indeed be neglected and has a negligible or weak impact on the emission of atomic species and on the immediate follow-up studies.

However, the range of visual extinction also suggest that self-shielding processes might be important. Indeed, as shown by Hollenbach & McKee (1989), the self-shielding of H₂ leads to the reformation of molecular hydrogen in the postshock when the temperature falls below 10⁴ K. As H₂ reforms, it efficiently cools the gas and triggers the formation of a wealth of molecules. The spectrum emitted by the shock not only contains strong rovibrational lines of H₂ but also other molecular lines. This was recently illustrated by Lehmann et al. (2020) who showed that non-irradiated shocks at intermediate velocities induce substantive molecular emission. The formation of molecules and the molecular tracers of intermediate- and high-velocity shocks propagating in the WNM, including the prediction of the full spectrum of molecular hydrogen, could be addressed in future works.

5.4 X-ray emission induced by charge exchange

While charge exchanges between singly ionized atoms and neutral species are included in the chemical network, we do not include the processes of charge exchange between strongly multi-ionized atoms and neutral species. Indeed, the recombination of strongly multi-ionized species is expected to be dominated by capture of free electrons rather than charge exchange with neutral species in stationary shocks. Moreover, the X-ray emission induced by this process has been shown to be negligible in shocks compared to the X-ray emission induced by collisional excitation by electron impact (Wise & Sarazin 1989). X-ray emission induced by charge exchange processes requires astrophysical environments where a neutral phase coexist with a highly ionized phase (e.g., Lallement 2004). Such mixing layers are out of the scope of this paper.

5.5 Limitations of 1D models

The model presented here is limited by its one-dimensional nature which constrains the final state of the postshock gas. Indeed, the postshock gas cannot escape in the x and y directions and is forced to be confined. It follows that the postshock trajectory predicted by 1D models necessarily becomes, at some point, unrealistic. This is true for shocks propagating in a single phase medium but also for shocks propagating in multiphase environments where the postshock gas is known to be subject to thermal instabilities (e.g., Falle et al. 2020). As shown by Lesaffre et al. (2020, see Fig. 21 in their paper), 1D models still provide robust predictions regarding the thermochemical state and the emission properties of shocks, providing that the postshock gas is not integrated indefinitely but up to a final point which, depending on the problem, can be chosen from the temperature or density profiles (e.g., Wardle 1999) or energetic considerations (e.g., Godard et al. 2019).

Another limitation of the model is that it cannot treat shocks propagating in inhomogeneous media which are known to trigger various dynamical instabilities (Falle et al. 2020; Raymond et al. 2020; Kupilas et al. 2021; Markwick et al. 2021).

6 Summary and conclusions

In this paper, we present a new version of the Paris-Durham shock code, which extends its domain of application to self-irradiated shocks at intermediate (30 < V_S < 100 km s⁻¹) and high (V_S ⩾ 100 km s⁻¹) velocity. The code is updated to take into account the production of multi-ionized species and the cooling induced by free-free emission and by the collisional excitation of the electronic levels of any atom in any ionization stage. This last part is done through the inclusion of all the levels, radiative transitions, and collisional rates contained in the CHIANTI database. The interactions of the photons produced by the shock with the radiative precursor and the shock itself are computed with a new, exact, and fast radiative-transfer algorithm for line emission based on the formalism of coupled escape probabilities. The chemical and thermodynamical structures of the radiative precursor and the shock are calculated self-consistently at steady-state with an iterative procedure that progressively accounts for the feedback of the self-generated radiation field on the thermal and chemical states of the gas.

The resulting model is explored over a large grid of parameters including low-, intermediate-, and high-velocity shocks propagating in poorly and highly magnetized (0.1 ⩽ B₀ ⩽ 10 μG) diffuse environments ($n_{\mathrm{H}}^0$ ⩽ 2 cm⁻³). As far as we know, this is the first time that multifluid shocks and shocks irradiated by themselves and by an external UV radiation field are studied over this range of velocities and densities. The exploration reveals a diversity of stationary solutions including C-type, CJ-type, and J-type shocks whose existence depends on the shock velocity compared to the speeds of the magnetosonic waves in the radiative precursor.

As the gas is shocked, it undergoes a well-defined trajectory in a proton density versus thermal pressure diagram. Depending on the physical conditions of the preshock gas and on the shock velocity, shocks can induce phase transition, that is the conversion of the initial WNM material into CNM material. For standard conditions of the WNM ($n_{\mathrm{H}}^0$ ⩾ 0.5 cm⁻³, B₀ ⩽ 1 μG), the criteria required to obtain phase transition (Eqs. (16) and (17)) are met by most of the models, indicating that interstellar shocks may play a fundamental role in the exchange of mass and energy between the warm and cold phases of the diffuse interstellar gas.

We find that most of the flux of mechanical energy is reprocessed into line emission over the entire grid of models and that the resulting spectrum hardens as the shock velocity increases. While the spectra of shocks propagating at low velocity (V_S ⩽ 30 km s⁻¹) is dominated by the emission of infrared photons, high-velocity shocks (V_S ⩾ 100 km s⁻¹) mostly emit EUV and X-ray photons through the deexcitation of the electronic levels of He⁺, and of different ionization stages of oxygen, carbon, nitrogen, and sulfur. These photons efficiently heat and modify the chemical state of the radiative precursor. In particular, the strength of the ionizing radiation field is found to be sufficient to fully ionize H, He, and C⁺ for V_S ≳ 100 km s⁻¹, fully ionize He⁺ for V_s ≳ 200 km s⁻¹, and fully ionize C²⁺ for V_S ≳ 300 km s⁻¹.

The thermal pressure reached by the shocked gas and the subsequent cooling trajectories lead to specific thermal and chemical properties and specific emission spectra whose signatures may already be present in ancillary observations of the diffuse interstellar matter. A detailed modeling of these phenomena therefore provides new clues for the interpretation of observations and the understanding of several fundamental mechanisms, including the injection and dissipation of mechanical energy in the ISM and the processes by which the gas is converted from the WNM to the CNM. Follow-up papers are currently underway to interpret ancillary data of atoms and ions observed in emission and absorption in the solar neighborhood. The production of molecules and molecular lines – which naturally arises during the phase transition process and the cooling of the postshock gas toward the CNM – will be addressed in future works.

Acknowledgements

We are very grateful to the referee for their thorough reading and their comments. The research leading to these results has received fundings from the European Research Council, under the European Community’s Seventh framework Programme, through the Advanced Grant MIST (FP7/2017–2022, No. 742719). The grid of models used in this work has been run on the computing cluster Totoro of the ERC MIST, administered by MesoPSL. The models were developed using the atomic data currently available on the CHIANTI database. CHIANTI is a collaborative project involving George Mason University, the University of Michigan (USA), University of Cambridge (UK) and NASA Goddard Space Flight Center (USA). We would also like to acknowledge the support from the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES.

Appendix A Radiative-transfer algorithm

Fig. A.1 Schematic view of the subdivision of a plane-parallel shock into zones. The red curve displays a typical temperature profile induced by an interstellar shock wave in arbitrary units as a function of the distance in logscale. A shock profile is divided into N zones. Each zone is labelled with an index i and is supposed to extend between the boundaries zi−1 and zi which are not homogeneously separated (as the scheme would lead to believe). The classical definitions of the angle θ and its cosine μ are used to account for the angular dependance of the radiation field in a plane-parallel geometry.

The iterative procedure used to compute the thermochemical structure of self-irradiated shocks relies on the proper calculation of the propagation of photons emitted by the shock and their interactions with the shocked and preshocked material. To perform this step we developed an exact radiative transfer code based on the formalism of coupled escape probability (CEP) for line emission (Elitzur & Asensio Ramos 2006).

Appendix A.1 Choice of the method

One of the main issue of radiative transfer is the treatment of the chicken-or-egg dilemma which arises from the fact that the radiation field propagating in a given direction in a cloud depends on the state of matter everywhere which itself depends on the radiation field propagating in every direction (Dullemond 2013). An exact and proven solution to this problem is the Accelerated Lambda Iteration (ALI) formalism in which the radiation field is computed through the iterative application of an operator Λ onto the source function (e.g., Olson et al. 1986). Because ALI algorithms and their developments for multilevel systems (MALI, Rybicki & Hummer 1991) can be numerically expansive, Elitzur & Asensio Ramos (2006) proposed a new and exact numerical method for solving the radiative transfer for line emission. In this method, designated as coupled escape probability, the effect of the radiation field is encoded into nonlocal couplings between the level populations at different positions in the cloud. At steady-state, the level populations are thus fully described by a set of non linear algebraic equations. Elitzur & Asensio Ramos (2006) showed that solving this system provides an exact solution to the radiative transfer without having to compute the radiation field itself, and for a computational cost substantially smaller than those of MALI algorithms.

In their original paper, Elitzur & Asensio Ramos (2006) developed the CEP for a plane-parallel geometry assuming a fixed line profile and applied it to standard problems with constant physical conditions. To solve the radiative transfer in plane-parallel interstellar shocks we propose here a modified version of the algorithm that takes into account spatial variations of the line profile and a new conditioning of the nonlocal coupling terms. To further reduce the computation cost, we finally apply additional simplifications that hold for diffuse material.

Appendix A.2 Formal solution for an emission line

Let’s consider the radiative transfer within a line u − l of a chemical species, where u designates the upper energy state and l the lower energy state of the line. In a plane-parallel geometry, the specific intensity at frequency ν along a ray, defined by its angle θ = acos(μ) with respect to the normal direction (see Fig. A.1), obeys the radiative transfer equation $$\mu \frac{d{I}_{\nu }(\mu ,z)}{dz}=j_{\nu }(\mu ,z)-{\kappa }_{\nu }(\mu ,z)I_{\nu }(\mu ,z),$$(A.1)

The emissivity j_ν and absorption coeffcient κ_ν write $$j_{\nu }(\mu ,z)=\frac{h{\nu }_{ul}}{4\pi }{A}_{ul}{n}_u(z){\varphi }_{\nu }(\mu ,z)$$(A.2)

and $${\kappa }_{\nu }(\mu ,z)=\frac{h{\nu }_{ul}}{4\pi }[B_{lu}{n}_l(z)-B_{ul}{n}_u(z)]{\varphi }_{\nu }(\mu ,z),$$(A.3)

where A_ul, B_ul, and B_lu are the Einstein coefficients of the line, n_u(z) and n_l(z) are the local densities of the upper and lower energy states of the species under consideration, and $${\varphi }_{\nu }(\mu ,z)=\frac{1}{\sqrt{2\pi }{\sigma }_{\nu }(z)}\exp [-\frac{1}{2}{\left(\frac{\nu -{\nu }_0(\mu ,z)}{{\sigma }_{\nu }(z)}\right)}^2]$$(A.4)

is the line profile. For the sake of simplicity, we consider here a Gaussian line profile centered at the local rest frequency $${\nu }_0(\mu ,z)={\nu }_{ul}(1+\frac{{\nu }_z(z)\mu }{c})$$(A.5)

and with a thermal broadening (in frequency space) $${\sigma }_{\nu }(z)=\sqrt{\frac{kT(z)}{m}}\frac{{\nu }_{ul}}{c}$$(A.6)

where υ_z(z) and T(z) are the local velocity and temperature of the fluid and m is the mass of the species under consideration. Assuming complete frequency redistribution, the source function of the line writes $$S(z)=\frac{\epsilon (z)}{\gamma (z)}=\frac{A_{ul}{n}_u(z)}{B_{lu}{n}_l(z)-B_{ul}{n}_u(z)},$$(A.7)

where $$\epsilon (z)=\frac{h{\nu }_{ul}}{4\pi }{A}_{ul}{n}_u(z)$$(A.8)

and $$\gamma (z)=\frac{h{\nu }_{ul}}{4\pi }[B_{lu}{n}_l(z)-B_{ul}{n}_u(z)]$$(A.9)

are the integrals over frequency of the emissivity and the absorption coefficient, respectively.

With these notations, and assuming that the cloud is not irradiated by an external radiation source, the radiative transfer equation (A.1) has the following formal solution $$\begin{array}{rl}& I_{\nu }(\mu ,z)={\int }_0^z{\mathrm{e}}^{-{\int }_y^z\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx}\frac{\epsilon (y){\varphi }_{\nu }(\mu ,y)}{\mu }dy\text{ for }\mu >0,\\ & I_{\nu }(\mu ,z)={\int }_{z_N}^z{\mathrm{e}}^{-{\int }_y^z\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }}dx\frac{\epsilon (y){\varphi }_{\nu }(\mu ,y)}{\mu }dy\text{ for }\mu <0\end{array}$$(A.10)

where z_N is the final distance of the cloud (see Fig. A.1). Now, let’s define the net radiative bracket of the line (Athay & Skumanich 1971) $$p(z)=1-\frac{\overline{J}(z)}{S(z)},$$(A.11)

where $$\overline{J}(z)=\frac{1}{2}\int d\nu {\int }_{-1}^1{\varphi }_{\nu }(\mu ,z)I_{\nu }(\mu ,z)d\mu$$(A.12)

is the profile integrated averaged intensity. The formal solution of the radiative-transfer equation (Eq. A.10) implies that $$\begin{array}{r}p(z)=1-\frac{1}{2S(z)}\int d\nu {\int }_0^1{\varphi }_{\nu }(\mu ,z)d\mu \\ \times {\int }_0^{z_N}{\mathrm{e}}^{-\left|{\int }_y^z\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx\right|}\frac{\epsilon (y){\varphi }_{\nu }(\mu ,y)}{\mu }dy.\end{array}$$(A.13)

This equation is equivalent to Eq. 7 of Elitzur & Asensio Ramos (2006) except that the integration is performed over z rather than the opacity and that the line profile is supposed to have an explicit dependance on μ and z. This expression of the net radiative bracket is at the root of the CEP approach. On the one hand, p(z) solely depends on the densities n_u(z) and n_l(z) across the cloud. On the other hand, the equations of evolution n_u(z) and n_l(z) solely depends on p(z) (see Sect. A.4). The net radiative bracket therefore contains all the nonlocal terms that couple the level populations at different positions in the cloud. In the following, we therefore develop the expression of p(z) over a discretized grid of positions (Sect. A.3) and then explicit the system of equations we need to solve for a multilevel species (Sect. A.4).

Appendix A.3 Net radiative bracket over a discretized grid

As described in Fig. A.1 and following Elitzur & Asensio Ramos (2006), the plane-parallel cloud (an interstellar shock in our case) is divided into N zones. The boundaries of the zones are designated by z_i. The i^th zone is defined as the zone occupying the range of distance z_i−1 < z ⩽ z_i. The distance between two boundaries z_i and z_j is designated as $$z^{i,j}=|z_i-z_j|\text{. }$$(A.14)

All the physical and chemical conditions within a zone i are assumed to be constant and are calculated as averaged of the actual profile between z_i−1 and z_i.

With these assumptions, the net radiative bracket for a specific line u − l in the zone i writes $$p^i=\frac{1}{z^{i-1,i}}{\int }_{z_{i-1}}^{z_i}p(z)dz.$$(A.15)

The combination of this equation with Eq. A.13 gives $$\begin{array}{r}{p}^i=1-\frac{1}{2S^i{z}^{i-1,i}}\sum _{j=1}^N{\epsilon }^j{\int }_{z_{i-1}}^{z_i}d{y}^{\mathrm{\prime }}{\int }_{z_{j-1}}^{z_j}dy\int d\nu \\ \times {\int }_0^{1}d\mu {\varphi }_{\nu }^i(\mu ){\mathrm{e}}^{-\left|{\int }_y^{y^{\mathrm{\prime }}}\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx\right|}\frac{{\varphi }_{\nu }^j(\mu )}{\mu },\end{array}$$(A.16)

where all the superscripts refer to physical quantities averaged over the corresponding zone. This last equation can be separated into $$p^i={\beta }^{i,i}+{\displaystyle \sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N}{\beta }^{i,j}$$(A.17)

where β^i,i corresponds to the contribution of the zone i to pⁱ and β^i,j corresponds to the contribution of the zone j to pⁱ.

Let’s first consider β^i,i. From Eqs. A.16 & A.17, we have $$\begin{array}{rl}{\beta }^{i,i}& =1-\frac{{\epsilon }^i}{2S{}^i{z}^{i-1,i}}{\int }_{z_{i-1}}^{z_i}d{y}^{\mathrm{\prime }}{\int }_{z_{i-1}}^{z_i}dy\int d\nu \\ & \times {\int }_0^{1}d\mu \frac{{({\varphi }_{\nu }^i(\mu ))}^2}{\mu }{\mathrm{e}}^{-\left|{\int }_y^{y^{\mathrm{\prime }}}\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx\right|}.\end{array}$$(A.18)

Performing the integration over y (treating successively the cases y ⩽ y′ and y > y′) and then over y′ leads to the following expression $${\beta }^{i,i}=\frac{1}{z^{i-1,i}{\gamma }^i}\int d\nu [\frac{1}{2}-{\int }_0^1\mu d\mu {\mathrm{e}}^{-\frac{{\gamma }^i{\varphi }_{\nu }^i(\mu )z^{i-1,i}}{\mu }}].$$(A.19)

Assuming that the line profile is independent of μ, this expression can be simplified further into $${\beta }^{i,i}=\frac{1}{z^{i-1,i}{\gamma }^i}\int d\nu [\frac{1}{2}-E_3\left({\gamma }^i{\varphi }_{\nu }^i{z}^{i-1,i}\right)],$$(A.20)

where E₃ is the third exponential integral. This final equation for β^i,i is completely equivalent to Eq. 22 of Elitzur & Asensio Ramos (2006).

Let’s now consider the nonlocal coupling term β^i,j. From Eqs. A.16 & A.17, we have $$\begin{array}{rl}& {\beta }^{i,j}=-\frac{{\epsilon }^j}{2S^i{z}^{i-1,i}}{\int }_{z_{i-1}}^{z_i}d{y}^{\mathrm{\prime }}{\int }_{z_{j-1}}^{z_j}dy\int d\nu \\ & \times {\int }_0^{1}d\mu \frac{{\varphi }_{\nu }^i(\mu ){\varphi }_{\nu }^j(\mu )}{\mu }{\mathrm{e}}^{-\left|{\int }_y^{y^{\mathrm{\prime }}}\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx\right|}.\\ & \end{array}$$(A.21)

Performing the successive integrations over y and y′ for j < i leads to $${\beta }^{i,j}=\frac{S^j}{2S^i{\gamma }^i{z}^{i-1,i}}\int dv{\int }_0^{1}d\mu \mu [{\alpha }^{i,j}-{\alpha }^{i-1,j}-{\alpha }^{i,j-1}+{\alpha }^{i-1,j-1}]$$(A.22)

where $${\alpha }^{i,j}={\mathrm{e}}^{-{\int }_{z_j}^{z_i}\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx}.$$(A.23)

To simplify this expression further, we make the approximation that the differences in Eq. A.22 can be replaced by derivatives, that is $$\begin{array}{rl}{\alpha }^{i,j}-{\alpha }^{i,j-1}& \approx z^{j-1,j}\frac{d}{dz}\left[{\mathrm{e}}^{-{\int }_z^{z_i}\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx}\right]\left(\frac{z_{j-1}+z_j}{2}\right)\\ & \approx z^{j-1,j}\frac{{\gamma }^i{\varphi }_{\nu }^j(\mu )}{\mu }{\mathrm{e}}^{-{\int }_{\frac{z_j-1+z_j}{2}}^{z_i}\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx}\end{array}$$(A.24)

and $$\begin{array}{rl}{\alpha }^{i-1,j}-{\alpha }^{i-1,j-1}& \approx z^{j-1,j}\frac{d}{dz}\left[{\mathrm{e}}^{-{\int }_z^{z_{i-1}}\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx}\right]\left(\frac{z_{j-1}+z_j}{2}\right)\\ & \approx z^{j-1,j}\frac{{\gamma }^i{\varphi }_{\nu }^j(\mu )}{\mu }{\mathrm{e}}^{-{\int }_{\frac{z_{j-1}+z_j}{2}}^{z_{i-1}}\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx}.\end{array}$$(A.25)

Performing the same approximation on the difference between Eqs. A.24 and A.25 and inserting the result in Eq. A.22, we obtain $$\begin{array}{rl}{\beta }^{i,j}=& -\frac{S^j{\gamma }^j{z}^{j-1,j}}{2S^i}\int d\nu {\int }_0^{1}d\mu \frac{{\varphi }_{\nu }^i(\mu ){\varphi }_{\nu }^j(\mu )}{\mu }\\ & \times {\mathrm{e}}^{-{\int }_{\frac{z_{j-1}+z_j}{2}}^{\frac{z_{i-1}+z_i}{2}}\frac{\gamma (x){\varphi }_{\nu }(\mu ,x)}{\mu }dx}.\end{array}$$(A.26)

We note that this equation is valid for j < i. The same reasoning applied to the case j > i virtually leads to the same expression except that the boundaries in the final integral are obviously inverted. This last equation has a deep physical meaning. It basically states that the amplitude of the coupling induced by zone j on zone i is proportional to the source function and the opacity⁵ of zone j, to the overlap of the line profiles of zones i and j, and to the probability of nonabsorption of the photons between zone i and j which is given by the exponential term. Assuming that the line profile is independent of μ, Eq. A.26 can be written in the more concise form, valid for any couple i and j, $${\beta }^{i,j}=-\frac{S^j{\gamma }^j{z}^{j-1,j}}{2S^i}\int d\nu {\varphi }_{\nu }^i{\varphi }_{\nu }^j{E}_1\left(\right|{\int }_{\frac{z_{j-1}+z_j}{2}}^{\frac{z_{i-1}+z_i}{2}}\gamma (x){\varphi }_{\nu }(x)dx\left|\right)$$(A.27)

where E₁ is the first exponential integral. It should be noted that this equation is very different from the expression of β^i,j given by Elitzur & Asensio Ramos (2006). This comes from the facts that we take into account the spatial variation of the line profile and that we replace the differences in Eq. A.22 by derivatives. Beside the physical insight described above, such an expression has the numerical advantage that β^i,j and the induced nonlocal coupling term are necessarily negative. It therefore provides a well behaved conditioning of the numerical problem.

Appendix A.4 Multilevel system and fluxes

Let’s consider a chemical species with L energy levels. In any zone i, the density of the species, nⁱ, naturally verifies $$n^i={\displaystyle \sum _{k=1}^L}{n}_k^i=n^i{\displaystyle \sum _{k=1}^L}{\mathrm{x}}_k^i$$(A.28)

where $n_k^i$ and ${\mathrm{x}}_k^i$ are the density and fractional population of the species in the level k and zone i. Taking into account excitations and de-excitations by collisions with electrons, spontaneous and stimulated emissions, and absorption processes, the fractional population of level k in zone i obeys the following evolution equation $$\begin{array}{rl}\frac{d{\mathrm{x}}_k^i}{dt}& =\sum _{u=k+1}^L{\mathrm{x}}_u^i{A}_{uk}{p}_{uk}^i+C_{uk}^i{n}_e^i({\mathrm{x}}_u^i-\frac{g_u}{g_k}{\mathrm{x}}_k^i{\mathrm{e}}^{-E_{uk}/k{T}^i})\\ & -\sum _{l=1}^{k-1}{\mathrm{x}}_k^i{A}_{kl}{p}_{kl}^i+C_{kl}^i{n}_e^i({\mathrm{x}}_k^i-\frac{g_k}{g_l}{\mathrm{x}}_l^i{\mathrm{e}}^{-E_{kl}/k{T}^i}),\end{array}$$(A.29)

where $p_{ul}^i$ is the net radiative bracket of the transition $u-l,n_e^i$ is the density of electrons, g_k is the degeneracy of level k, and $C_{ul}^i$ is the de-excitation rate (in cm³ s⁻¹) from level u to level l by collision with electrons. At steady-state, the above equations form a large set of non linear algebraic equations. Because all the nonlocal couplings are contained in $p_{ul}^i$, solving this system provides an exact solution to the radiative transfer. The computation of the intensity and the flux of any line is then straightforward. From the formal solution of the radiative transfer (Eq. A.10), the specific intensity of a line u − l at boundary z_j writes $$I_{\nu ,ul}(\mu ,z_j)={\displaystyle \sum _{i=1}^j}{S}_{ul}^i[{\alpha }_{ul}^{j,i}-{\alpha }_{ul}^{j,i-1}]\text{ for }\mu >0$$(A.30)

and $$I_{\nu ,ul}(\mu ,z_j)=-{\displaystyle \sum _{i=j+1}^{z_N}}{S}_{ul}^i[{\alpha }_{ul}^{j,i}-{\alpha }_{ul}^{j,i-1}]\text{ for }\mu <0.$$(A.31)

Again, replacing differences by derivatives (as we did from Eq. A.22 to Eq. A.26), we obtain $$\begin{array}{r}{I}_{\nu ,ul}(\mu ,z_j)=\sum _{i=1}^j{S}_{ul}^i{z}^{i-1,i}\frac{{\gamma }_{ul}^i{\varphi }_{\nu ,ul}^i(\mu )}{\mu }\\ \times {\mathrm{e}}^{-{\int }_{\frac{z_{i-1}+z_i}{2}}^{z_j}\frac{{\gamma }_{ul}(x){\varphi }_{\nu ,ul}(\mu ,x)}{\mu }dx}\text{ for }\mu >0.\end{array}$$(A.32)

and $$\begin{array}{r}{I}_{\nu ,ul}(\mu ,z_j)=-\sum _{i=j+1}^N{S}_{ul}^i{z}^{i-1,i}\frac{{\gamma }_{ul}^i{\varphi }_{\nu ,ul}^i(\mu )}{\mu }\\ \times {\mathrm{e}}^{-{\int }_{\frac{z_{i-1}+z_i}{2}}^{z_j}\frac{{\gamma }_{ul}(x){\varphi }_{\nu ,ul}(\mu ,x)}{\mu }dx}\text{ for }\mu <0\text{. }\end{array}$$(A.33)

If the line profile is independent of μ, the flux densities in the forward and backward directions at boundary z_j, defined as $$F_{\nu ,ul}^{+}\left(z_j\right)=2\pi {\int }_0^1{I}_{\nu ,ul}(\mu ,z_j)\mu d\mu$$(A.34)

and $$F_{\nu ,ul}^{-}\left(z_j\right)=2\pi {\int }_{-1}^0{I}_{\nu ,ul}(\mu ,z_j)\mu d\mu ,$$(A.35)

simply write $$F_{\nu ,ul}^{+}\left(z_j\right)=2\pi {\displaystyle \sum _{i=1}^j}{S}_{ul}^i{z}^{i-1,i}{\gamma }_{ul}^i{\varphi }_{\nu ,ul}^i{E}_2({\int }_{\frac{z_{i-1}+z_i}{2}}^{z_j}{\gamma }_{ul}(x){\varphi }_{\nu ,ul}(x)dx)$$(A.36)

and $$F_{\nu ,ul}^{-}\left(z_j\right)=2\pi {\displaystyle \sum _{i=j+1}^N}{S}_{ul}^i{z}^{i-1,i}{\gamma }_{ul}^i{\varphi }_{\nu ,ul}^i{E}_2({\int }_{z_j}^{\frac{z_{i-1}+z_i}{2}}{\gamma }_{ul}(x){\varphi }_{\nu ,ul}(x)dx),$$(A.37)

where E₂ is the second exponential integral. These previous equations for the intensities and the flux densities are different from Eqs. 19 and 20 of Elitzur & Asensio Ramos (2006). This is, again, because we approximate differences by derivatives and because we include a position dependent line profile. We stress that for an optically thin line, the second exponential integral is close to 1. In this case, it is remarkable to note that, while the specific intensity diverges for small values of μ, the forward and backward flux densities are well behaved and give exactly the value we would obtain by integrating the emissivity along the μ = 1 direction and multiplying the result by 2π.

A final remark should be made regarding the balance of energy. All the equations above have been derived considering that the sole source of energy is thermal. The level populations calculated over the entire cloud therefore necessarily fullfil the energy conservation equation $$\begin{array}{r}\sum _{i=1}^N\sum _{u=1}^L\sum _{l=1}^{u-1}{E}_{ul}{z}^{i-1,i}{C}_{ul}^i{n}_e^i(n_l^i\frac{g_u}{g_l}{\mathrm{e}}^{-E_{ul}/k{T}^i}-n_u^i)=\\ \sum _{u=1}^L\sum _{l=1}^{u-1}{\int }_{\nu }(F_{\nu ,ul}^{-}(0)+F_{\nu ,ul}^{+}\left(z_N\right))d\nu \end{array}$$(A.38)

where the left hand side corresponds to the total cooling flux induced by collisions, and the right hand side corresponds to the total line flux escaping the cloud. This equation provides a fundamental sanity check of the results obtained from the numerical scheme used to solved the multilevel system.

Appendix A.5 Approximations for diffuse interstellar shocks

The general methodology presented in the previous sections is used to postprocess the radiative transfer in interstellar shocks within all the lines of any atom in any ionization state. As described above, the method relies on the subdivision of the plane-parallel structure into zones. In the case of interstellar shocks, the number of zones N and their boundaries z_i are chosen so that the density of the species under consideration, nⁱ, the increment of total column density, $n_{\mathrm{H}}^i{z}^{i-1,i}$, the density of electrons, $n_e^i$, the temperature of the gas, Tⁱ, and its velocity, $v_z^i$, vary by less than 10% from one zone to the next.

In practice the nonlinear system given by Eqs. A.17, A.20, A.27, A.28, and A.29 could be solved, at steady-state, using a Newton-Raphson algorithm. However, an additional simplification can be made by noticing that in the diffuse interstellar medium, only the ground electronic state of atoms and atomic ions are substantially populated, regardless of the gas temperature. It implies that the fractional population of the ground electronic state can be set to one and that only the dipole-allowed electronic lines connected to this level should be considered as resonant and treated with a detailed calculation of their net radiative bracket. All the other transitions, including the two-photon emission processes, can be considered as optically thin with a net radiative bracket $p_{ul}^i$ = 1 for any zone i. With this approximation, the original system of N × L nonlinear equations is transformed into L systems of N linear equations that can be solved sequentially, starting from the highest level, L, and going downward.

Fig. A.2 Flux densities of the Lyα (top panel) and Lyβ (bottom panel) lines of atomic hydrogen obtained by applying the exact radiative-transfer algorithm of Appendix A on an interstellar shock. The flux densities are computed in the backward direction, that is in the direction of the preshock material, at the position z = 0. The model displayed corresponds to a shock at 50 km s−1 propagating in a medium with a density of 1 cm−3 and a magnetic field strength of 1 μG, and irradiated by an external UV radiation field with a scaling factor G0 = 1 (in Mathis units).

The method described in the previous sections has been developed neglecting the impacts of the external radiation field and the shock continuum emission (through Bremsstrahlung and recombination radiation) on the radiative pumping of excited states. As discussed by Elitzur & Asensio Ramos (2006) (in their Appendix A), these effects are not particularly difficult to include and simply imply to add pumping rates from the ground electronic states in Eq. A.29. In this paper, we deliberately ignore this aspect in order to check that the method ensures the conservation of energy (Eq. A.38). The impact of external and continuum radiations onto the level excitations will be included in future works if the method proves to be effective in a wide variety of environments.

It may appear surprising that all the equations are written without including the absorption induced by photoionization and dust and PAH extinctions. Solving this problem self-consistently is, in fact, straightforward and requires to include the opacity of each process computed in the μ = 1 direction, over the zone i in the exponential integral of Eq. A.20, and between zones i and j in the exponential integrals of Eqs. A.27, A.36, and A.37. Despite the simplicity of the task, we choose here to adopt the following approach: the level populations and the line fluxes are first calculated without these processes over the entire shock profile ; the line fluxes are then computed, in the shock and the preshock media, including gas and dust absorptions. This choice is done, again, to check that the method effectively ensures the conservation of energy (Eq. A.38). The self-consistent treatment will be adopted in future works once the method has proven successful over a larger set of parameters than that explored here. We stress, however, that this self-consistent treatment should not impact the results of this paper. It is so because the absorption lengths of gas and dust are not homogeneous but display extremely sharp transitions in the shock and the preshock materials (see Appendix B, Fig. B.2).

Appendix A.6 Example of line profiles

As an example of application of the postprocessing radiative-transfer algorithm, we display in Fig. A.2 the flux densities of the Lyα and Lyβ lines of atomic hydrogen in the backward direction (i.e., escaping toward the preshocked gas) obtained for a shock propagating at 50 km s⁻¹ in a medium with a density of 1 cm⁻³ and a magnetic field strength of 1 μG. As expected in self-irradiated shocks, both the Lyα and Lyβ lines display complex and asymmetric profiles. Because of the large opacity of atomic hydrogen across the shock, Lyα and Lyβ photons can only escape in the wings of the line. Moreover, the flux density of the Lyβ line is found to be about 20 times smaller than that of the Lyα line. This is the signature of the conversion of Lyβ photons into Hα and Lyα photons which occur in environments with large opacities. All these features are in agreement with the line profiles we recently computed using a MALI algorithm for shocks at similar velocities, yet at significantly larger densities (see Fig. 6 of Lehmann et al. 2020).

Appendix B Photon-matter interactions

The propagation of photons in the shock and the preshock media is controlled by their interactions with the interstellar matter and the cross sections of the associated processes. The main interactions of UV, EUV, and X-ray photons in atomic environments are the absorption by PAHs and dust and the photoionization of atoms and atomic ions. These processes not only drive the physical and chemical states of the shocked and the preshocked material but also determine the propagation scale of the radiation. Following our previous work (Godard et al. 2019), the extinction by PAHs and dust is modeled using the absorption coefficients of spherical graphite grains derived by Draine & Lee (1984) and Laor & Draine (1993). All PAHs are assumed to have a constant size of 6.4 Å. The size of dust particles is computed from the mean square radius of grains taking into account the erosion of grain cores (see Appendix A of Godard et al. 2019). The photoionization of gas phase species is modeled using the prescription of Verner & Yakovlev (1995) who performed analytical fits of the partial ionization cross sections of atoms and atomic ions for all subshells nl of the associated ground electronic state. The total ionization cross section of any species is the sum of these partial ionization cross sections. As an illustration, we display in Fig. B.1 the total photoionization cross sections of the several elements, in their neutral form and in their first ionization stage, as functions of the photon energy. In addition to the expected asymptotic behavior (in ν^−3.5 at large ν; e.g., Dopita & Sutherland 2003), this figure highlights the sudden jumps of the cross sections at specific energy thresholds. These jumps are due to the openings of new ionization channels corresponding to the removal of one electron (and possibly several Auger’s electrons) from the inner subshells of the species.

Fig. B.1 Photoionization cross sections of neutral atoms (top panel) and their first ionization stage (bottom panel) as functions of the photon energy, derived from the analytical formula of Verner & Yakovlev (1995) (Eq. 1 and Table 1 in their paper). The discontinuities are due to jumps of the cross sections when the photon energy becomes larger than the different thresholds required to extract an electron from internal subshells of the species.

The impacts of PAHs, dust, and gas phase particles on the absorption of the radiation field is shown in Fig. B.2 which displays the local absorption length of photons in the preshocked and shocked materials as a function of the photon wavelength, computed with the standard shock model. In high-velocity shocks, photons with wavelength below 911 Å are subject to two walls of absorption, one in the preshock, and one in the postshock. This feature is the signature of the chemical composition of the gas. Between these walls, hydrogen and helium are mostly ionized and EUV and X-ray photons are mostly absorbed by ionized metals (in the preshock), and PAHs and dust (in the shock). Beyond these walls, the gas is mostly neutral and high energy photons are quickly absorbed by both H and He. In contrast, photons with wavelength above 911 Å are subjected to a single wall of absorption in the postshock. This comes from the fact that UV photons are preferentially absorbed by dust and PAHs. Since the density is roughly constant in the preshock medium, so is the absorption length. The wall seen in the postshock medium corresponds to the position where the gas cools from a few 10⁵ K to a few tens of K (see Fig. 3) and is induced by the isobaric increase of the density up to the point where the gas is supported by the magnetic pressure.

Fig. B.2 Local absorption length of photons obtained in the final converged iteration of the standard shock model. The absorption length, given in pc with the color scale, is shown as a function of the photon energy and the distance in the preshock gas (left panel) and in the shock (right panel). Above 911 Å, photons are preferentially absorbed by dust and PAHs. Below this limit, photons are preferentially absorbed by atoms, in particular H, He, and He+. The walls seen by the photons below 911 Å correspond to the regions in the preshock and the postshock gas where neutral hydrogren is reformed.

References

All Tables

All Figures

Fig. 1 Thermal equilibrium state of the diffuse interstellar gas obtained with the Paris-Durham code and displayed in a particle density versus thermal pressure diagram. The color code highlights the three branches of the equilibrium state: the dynamically stable WNM (red) and CNM (blue) phases at ~104 K and ~80 K, respectively, and the dynamically unstable branch (green) at intermediate temperatures. The light blue and red zones show the regions of the diagram where the cooling rate is larger or smaller, respectively, than the heating rate. Light gray lines are isothermal contours from 1 to 106 K (from bottom right to top left).

In the text

Fig. 2 Timescale required to reach the ionization, thermal, and chemical equilibrium states in the diffuse medium as a function of the proton density nH (red curve), and corresponding equilibrium values of the ionization fraction, xe = n(e−)/nH (blue curve). The equilibrium timescale is estimated as the time at which the kinetic temperature, the electronic fraction, and the chemical abundances of all species are within 10% of their equilibrium values, starting from random initial conditions.

In the text

Fig. 3 Convergence of the standard model on the shock thermochemical profiles. Profiles of the kinetic temperature (top panels), the fractional abundances of H and H+ (middle panels), and the fractional abundances of He, He+, and He++ (bottom panels) as functions of the distance in the preshock (left panels) and in the shock (right panels). The black curves show the intermediary profiles computed during the iterative procedure and the colored curves the profiles obtained in the final converged iteration.

In the text

Fig. 4 Production of multi-ionized species in the standard model. Abundances profiles of carbon-bearing species (top panels), nitrogen-bearing species (middle panels), and oxygen-bearing species (bottom panels) as functions of the distance in the preshock (left panels) and in the shock (right panels). The high ionization stages of carbon, nitrogen, and oxygen (C6+, N6+, N7+, O7+, and O8+) are not included because their relative abundances are all below the lowest value shown in this figure.

In the text

Fig. 5 Trajectory of a shocked fluid particle from the ambient medium to the postshock region obtained in the standard model. The trajectory is displayed in a particle density versus thermal pressure diagram (top panel) and in a proton density versus thermal pressure diagram (bottom panel) to highlight the phase transition from the WNM to the CNM induced by a shock at 200 km s−1. A fluid particle initially evolves through the radiative precursor (pink curve). As it crosses a shock front, it undergoes three successive regimes: an adiabatic heating (red curve) followed by a quasi-isobaric cooling (green curve), and finally a quasi-isochoric cooling (blue curve) as the magnetic pressure becomes dominant in the postshock gas. The black curve indicates the thermal equilibrium state of the diffuse gas obtained for G0 = 1 (see Fig. 1). Light gray lines are isothermal contours from 1 to 107 K (from bottom right to top left). We note that these isocontours are evenly spaced in the top panel but not in the bottom panel. Tics on the trajectories indicate the time, in log scale, starting from the adiabatic phase, ranging from 10−1 yr to 106 yr.

In the text

Fig. 6 Energy flux, computed as $\nu F_{\nu }^{\text{int }}$, emitted toward the preshock at z = 0 in the standard model as a function of the wavelength of the emitted photon. The black curve corresponds to the total continuum emission (including Bremsstrahlung, bound-free, and two-photon processes) and the colored spikes to line emissions. Spikes are colored according to the element responsible for the line emission regardless of its ionization state. The top panel displays the full wavelength range while the bottom panel is a zoom onto the range between 10 and 2000 Å. The black vertical line indicates the ionization threshold of hydrogen.

In the text

Fig. 7 Critical velocities of interstellar shocks for $n_{\mathrm{H}}^0$ = 0.1 cm−3 (top panel) and 1 cm−3 (bottom panel). The solid and dash-dotted lines show the magnetosonic speed computed in the ionized fluid (Eq. (19)) and the neutral fluid (Eq. (18)), respectively. The colored areas highlight the range of velocities over which stationary shocks are C-type (red), CJ-type (blue), and J-type (green) shocks. The transition between C-type and CJ-type shocks (dashed line) is estimated as the point where the maximum postshock thermal pressure is equal to the maximum magnetic pressure allowed by an adiabatic compression. The gray points indicate the values of the shock velocity and transverse magnetic field strength explored in the grid of models. The red dot highlights the position of the standard model.

In the text

Fig. 8 Trajectories followed by interstellar shocks in a proton density versus thermal pressure diagram (see Fig. 5). Each panel corresponds to different values of the initial density $n_{\mathrm{H}}^0$ (top labels) and transverse magnetic field B0 (right labels). The different curves correspond to the following values of the shock velocity ordered from top to bottom trajectories: 500, 450, 400, 350, 300, 250, 200, 150, 120, 100, 90, 80, 70, 60, 50, 40, 30, 20, and 10 km s−1. Colored trajectories correspond to J-type shocks, gray trajectories to CJ-type shocks, and black trajectories to C-type shocks. The black curve common to all panels indicates the thermal equilibrium state of the diffuse gas as a function of the proton density obtained for G0 = 1 (see Fig. 1). Light gray lines are isothermal contours from 1 to 107 K (from bottom right to top left). We note that the trajectories of low-velocity shocks are not always visible as they only slightly perturb the gas from its original state.

In the text

	Fig. 9 Fractional abundances of H and H+ (top panel), He, He+, and He++ (middle panel), and C+, C2+, C3+, C4+, C5+, and C6+ (bottom panel) in the gas entering the shock as functions of the shock velocity. Except for the velocity, all other parameters are set to their standard values (see Table 1).
In the text

Fig. 10 Energy flux, computed as $\nu F_{\nu }^{\text{int }}$, emitted toward the preshock at z = 0 as a function of the wavelength of the emitted photon. The black curve corresponds to the total continuum emission (including Bremsstrahlung, bound-free, and two-photon processes) and the colored spikes to line emissions. Spikes are colored according to the element responsible for the line emission regardless of its ionization state. Predictions are shown for shocks at different velocities, ranging from 10 km s−1 (top left panel) to 400 km s−1 (bottom right panel). All other parameters are set to their standard values (see Table 1). The dashed gray vertical line indicates the ionization threshold of hydrogen.

In the text

	Fig. 11 Fraction of the radiation emitted by shocks in the X-ray and EUV domains (red curve), in the ultraviolet range (green curve), and in the optical and the infrared domains (blue curve) as functions of the shock velocity for $n_{\mathrm{H}}^0$ = 1 cm−3 and B0 = 1 μG.
In the text

Fig. 12 Strength of the self-generated UV radiation field (top panel) and ionization rate of hydrogen induced by the self-generated EUV and X-ray radiation field (bottom panel) in the radiative precursor, at z = 0. Both quantities are shown as functions of the shock velocity for different values of the preshock density and for B0 = 1 μG. The strength of the UV radiation field is calculated from the integral of $I_{\nu }^{\text{int }}$, over solid angle and frequency between 911 and 2400 Å, normalized to the value obtained by Mathis (Mathis et al. 1983). The strength of the external UV radiation field used in this work (G0 = 1) is shown as a gray line on the top panel. The cosmic ray ionization rate of atomic hydrogen used in this work $({\zeta }_{\mathrm{H}}={\zeta }_{{\mathrm{H}}_2}/2.3=1.3\times {10}^{-16}{\mathrm{s}}^{-1})$ is shown as a gray line on the bottom panel.

In the text

	Fig. 13 Ionization rates of H, He, C+, O, N, S+, Si+, and Fe+ induced by the standard shock in the preshock gas as functions of the preshock distance. These species are selected because they are the main carriers of the elements included in the model in the WNM.
In the text

Fig. A.1 Schematic view of the subdivision of a plane-parallel shock into zones. The red curve displays a typical temperature profile induced by an interstellar shock wave in arbitrary units as a function of the distance in logscale. A shock profile is divided into N zones. Each zone is labelled with an index i and is supposed to extend between the boundaries zi−1 and zi which are not homogeneously separated (as the scheme would lead to believe). The classical definitions of the angle θ and its cosine μ are used to account for the angular dependance of the radiation field in a plane-parallel geometry.

In the text

Fig. A.2 Flux densities of the Lyα (top panel) and Lyβ (bottom panel) lines of atomic hydrogen obtained by applying the exact radiative-transfer algorithm of Appendix A on an interstellar shock. The flux densities are computed in the backward direction, that is in the direction of the preshock material, at the position z = 0. The model displayed corresponds to a shock at 50 km s−1 propagating in a medium with a density of 1 cm−3 and a magnetic field strength of 1 μG, and irradiated by an external UV radiation field with a scaling factor G0 = 1 (in Mathis units).

In the text

Fig. B.1 Photoionization cross sections of neutral atoms (top panel) and their first ionization stage (bottom panel) as functions of the photon energy, derived from the analytical formula of Verner & Yakovlev (1995) (Eq. 1 and Table 1 in their paper). The discontinuities are due to jumps of the cross sections when the photon energy becomes larger than the different thresholds required to extract an electron from internal subshells of the species.

In the text

Fig. B.2 Local absorption length of photons obtained in the final converged iteration of the standard shock model. The absorption length, given in pc with the color scale, is shown as a function of the photon energy and the distance in the preshock gas (left panel) and in the shock (right panel). Above 911 Å, photons are preferentially absorbed by dust and PAHs. Below this limit, photons are preferentially absorbed by atoms, in particular H, He, and He+. The walls seen by the photons below 911 Å correspond to the regions in the preshock and the postshock gas where neutral hydrogren is reformed.

In the text