© The Authors 2024

1 Introduction

The local interstellar medium (ISM) is full of chemical mysteries. Since the first detection of molecules in space (e.g. Adams 1941; Douglas & Herzberg 1941), observations of atoms and molecules have revealed chemical, excitation, and kinematic properties that regularly challenge idealised descriptions of the diffuse neutral gas. Among these challenges are the observations of CH⁺ (e.g. Gredel et al. 2002; Pan et al. 2004; Sheffer et al. 2008), of SH+ (Godard et al. 2012), and of the pure rotational levels of H₂ (e.g. Wakker 2006; Shull et al. 2021), whose abundances are several orders of magnitude above those predicted by a chemistry solely driven by ultraviolet (UV) photons, the bimodality of the line profiles of CH (Lambert et al. 1990; Crane et al. 1995), which suggests that this species is formed through two separated chemical pathways, or the tight correlation between the abundances of HCO⁺ and OH or H₂O (Lucas & Liszt 1996; Gerin et al. 2019), which cannot be explained by static clouds at chemical equilibrium. As has been shown in several studies, all of these observations are precious as they may be diagnostics of the dynamical evolution (condensation and evaporation) of the multiphase ISM (Valdivia et al. 2017; Godard et al. 2023), the dissipation processes of mechanical energy in the cold neutral medium (Godard et al. 2014; Lesaffre et al. 2020), or the out-of-equilibrium chemical evolution of molecular clouds (Panessa et al. 2023). More generally, these features detected in the local ISM are of fundamental importance to identifying the energy sources that power atomic and molecular emission lines which, in turn, can be used to interpret extragalactic observations and draw the mechanical and radiative energy balance of the ISM in entire galaxies (Lehmann et al. 2022; Villa-Vélez et al. 2024).

Among these chemical conundrums are the intriguing excitation conditions of neutral carbon observed in the local diffuse interstellar gas. Using the high-resolution Space Telescope Imaging Spectrograph aboard the Hubble Space Telescope, Jenkins & Tripp (2001) and Jenkins & Tripp (2011, hereafter, JT11) conducted several observations of the UV multiplets of CI in absorption against the continuum of nearby OB stars. These absorption lines, which originate from the three fine structure levels of neutral carbon, provide the column-averaged ratios of the level populations, and hence clues regarding the density and temperature of the foreground absorbing material. Through a detailed analysis of the excitation of neutral carbon, JT11 showed that these observations can only be explained if the lines of sight are composed of two distinct environments: a medium at low thermal pressure (P ~ 3.8 × 10³ K cm⁻³) that encompasses ~95% of the mass of the gas and corresponds to the bulk of the ISM (Wolfire et al. 2003), and a medium at a thermal pressure one to three orders of magnitude above. The coexistence along every line of sight of environments with large contrasts in thermal pressure raises a deep interpretative challenge regarding the dynamical state of the diffuse interstellar matter, the origin of the high-pressure component, its survival timescale, and the absence of gas at intermediate thermal pressure.

In their analysis, JT11 originally proposed that the gas at high thermal pressure could originate from the dissipation of turbulent energy in the cold neutral medium (CNM; Joulain et al. 1998; Godard et al. 2009), from the propagation of shock waves in the warm neutral medium (WNM; Bergin et al. 2004), or, alternatively, from the interactions of the background stars with their surrounding environments, which include the expansion of HII regions and the pressure induced by the recoil of H atoms following the photodissociation of H₂ at the edge of molecular clouds (Field et al. 2009). In a follow-up study, Jenkins & Tripp (2021) finally ruled out the influence of the background stars by observing CI in absorption towards extragalactic sources and finding similar fractions of diffuse gas at high thermal pressure along these lines of sight. This lack of influence of the background stars is a key result. It suggests that the excitation of neutral carbon is linked to a more general feature of the dynamical evolution of the diffuse matter.

The fact that shocks associated with supernova remnants (SNRs) produce gas at high thermal pressure is well known. Analysis of UV absorption lines towards stars in the IC 443 and the Vela SNRs (e.g. Ritchey et al. 2020; Ritchey 2023 and reference therein) have long revealed the presence of shocked gas at thermal pressures as high as ~10⁶ K cm⁻³. A step forwards was accomplished by Bergin et al. (2004), who showed that shocks driven in the diffuse atomic gas by spiral arm density waves, supernovae blast waves, or the infall of material onto the galactic disc could explain the excitation state of neutral carbon observed by JT11. Indeed, this work promotes the idea that shocks influence the dynamics of the entire ISM and are detectable towards lines of sight with no obvious associated SNR. If this scenario is confirmed, the observations of CI would provide invaluable information regarding the mechanisms by which mechanical energy is injected into the ISM, which are still uncertain and currently believed to be a mixture of stellar feedback and galactic dynamics including orbital motions and mass accretion onto galactic discs (Brucy et al. 2020). However, and to our knowledge, the scenario proposed by Bergin et al. (2004) has never been studied quantitatively.

In this paper, we present a detailed and quantitative analysis of the production and excitation of neutral carbon by shocks propagating at intermediate (30 < V_S ⩽ 100 km s⁻¹) and high velocities (V_S > 100 km s⁻¹) in the WNM. These shocks, which incidentally induce phase transition between the WNM and the CNM, are modelled using the most recent version of the Paris-Durham shock code designed to follow the dynamical, thermal, and chemical evolutions of shocks with velocities up to 500 km s⁻¹ (Godard et al. 2024, hereafter Paper I). Following JT11, the observations of the absorption lines of CI are decomposed into a low-pressure component and a high-pressure component. The column densities of CI at high pressure are then compared with the predictions of the model to derive the distribution of dissipation rates required to explain the observations.

The observational sample and the decomposition procedure of JT11 are presented in Sec. 2. The model and its predictions regarding the abundance, the excitation, and the line profiles of CI are presented in Sec. 3. The resulting distribution of dissipation rates and its dependence on the model parameters are shown in Sec. 4. The origin of these shocks that likely, but not exclusively, involves SNRs expanding in the diffuse phases of the Galaxy is discussed in Sec. 5, in which we highlight the unavoidable nature of SNR surfaces in the observational sample of JT11. Some limitations of this work are exposed and discussed in Sec. 6. The main conclusions are summarised in Sec. 7.

Fig. 1 Aitoff projection, in Galactic longitude and latitude coordinates, of the background sources of the observational sample of JT11. The colour-coding indicates the distance of background sources in parsecs.

2 Observational sample

2.1 Positions and distances

The observational dataset studied here contains all the lines of sight observed by JT11 and available on the Strasbourg astronomical Data Center¹. The sample is composed of 89 lines of sight in which the UV multiplets of CI have been detected in absorption against the continuum of nearby OB stars. The position of the sources in Galactic coordinates and their distances derived from measurements of their parallaxes are shown in Fig. 1. The background sources cover all Galactic longitudes with two clusters between −120° and −30° and between 90° and 120°, which are associated with the Scutum-Centaurus, the Sagittarius, and the Perseus spiral arms (see Sec. 5). Similarly to the sample studied by Bellomi et al. (2020), most of the background sources are located at Galactic latitudes below 15°, following the distribution of stars in the solar neighbourhood.

The background stars are located at distances between 120 pc and 6.1 kpc. Since the amount of gas exponentially decreases as a function of the distance from the midplane of the Galaxy (Dickey & Lockman 1990), the length, l_los, of the neutral material intercepted by any line of sight is necessarily smaller than the distance of the background source, especially when the source is located at high Galactic latitude. Following the simple prescription adopted by Bellomi et al. (2020) and Godard et al. (2023), we assume that the gas extends above the midplane over a height of 100 pc and estimate this length as $$l_{\mathrm{l}\mathrm{o}\mathrm{s}}=min(d,\frac{100}{\sin (|b|)})\mathrm{p}\mathrm{c},$$(1)

where d is the distance of the background source and b is its Galactic latitude. With this prescription, the length of the intercepted material, l_los, is found to follow a flat distribution in logscale between about 200 pc and 3 kpc. This length will be used in Sects. 4 and 5 to compute the distribution of dissipation rate in shocks propagating in the WNM. The impact of the prescription chosen to compute l_los will be discussed then.

2.2 Excitation conditions deduced from CI lines

The UV spectra observed by JT11 are composed of absorption lines originating from the three fine structure levels of the ground electronic state of neutral carbon, the ³P₀ fundamental level, and the ³P₁ and ³P₂ excited levels lying, respectively, at 23.62 K and 62.46 K above the fundamental level. Each spectrum contains several UV multiplets that can be used to solve the confusion induced by the mixture of line opacities. Using this information and a method devised by Jenkins & Tripp (2001), JT11 derived the column density spectra of each level, per velocity interval of 0.5 km s⁻¹, along any line of sight. These final products of their analysis (see Tables 3 and 4 in their paper) are the data used in this work. Following partially their notation, we define $N\left({\mathrm{C}}_{{}^3{\mathrm{P}}_0}\right),N\left({\mathrm{C}}_{{}^3{\mathrm{P}}_1}\right)$, and $N\left({\mathrm{C}}_{{}^3{\mathrm{P}}_2}\right)$ as the column densities of the three fine structure levels of carbon integrated over the entire line of sight, and $$f_1=\frac{N\left({\mathrm{C}}_{{}^3{\mathrm{P}}_1}\right)}{N(\mathrm{C})}$$(2)

and $$f_2=\frac{N\left({\mathrm{C}}_{{}^3{\mathrm{P}}_2}\right)}{N(\mathrm{C})}$$(3)

as the ratios of the column densities of the two excited levels to the total column density of carbon, $$N(\mathrm{C})=N\left({\mathrm{C}}_{{}^3{\mathrm{P}}_0}\right)+N\left({\mathrm{C}}_{{}^3{\mathrm{P}}_1}\right)+N\left({\mathrm{C}}_{{}^3{\mathrm{P}}_2}\right).$$(4)

Figure 2 displays the f₁ and f₂ ratios computed along all lines of sight of the observational sample (coloured points), the ratios computed per velocity interval (small black points), and the total ratios obtained over the entire sample (green point). For comparison, Fig. 2 also shows the theoretical f₁ and f₂ ratios obtained in homogeneous clouds at different proton densities, n_H, and kinetic temperatures, T, calculated at equilibrium taking into account excitation and de-excitation of neutral carbon by radiative decay, collisions with electrons, H, H₂, and He, radiative pumping of the electronic states of C, and chemical pumping of the fine structure levels of C during the recombination of C⁺ (see Sec. 3 for more details). This convoluted figure is similar to Fig. 2 of JT11. Although the exact predictions of the theoretical model depends on the abundances of the collisional partners and on the strength of the impinging UV radiation field², the results obtained here are generic and reveal three main features.

First, Fig. 2 shows that almost none of the observed excitation conditions can be explained with a homogeneous cloud at a single thermal pressure. While it is not explicitly shown in the figure, this result holds regardless of the abundances of the collisional partners and of the strength of the UV radiation field (see JT11). Second, lines of sight with larger column densities of neutral carbon display f₁ and f₂ ratios that lie closer to the total ratios derived over the entire sample (green point). This feature is indicative that all lines of sight intercept similar environments with excitation properties that tend towards the mean values when the column density increases by virtue of the central limit theorem. Last, the theoretical predictions of the f₁ and f₂ ratios obtained at different kinetic temperatures and thermal pressures (from 10² to 10⁷ K cm⁻³) appear to form an envelope around the observational points. All these features imply that any line of sight is necessarily composed of gas with markedly different physical conditions and that the observed column densities result from the sum of these different environments.

Fig. 2 Excitation properties of the neutral carbon along the lines of sight observed by JT11. The filled coloured circles indicate the f1 and f2 column density ratios observed along each line of sight. Each circle is colour-coded according to the total column density of neutral carbon, N(C). The black dots show the f1 and f2 column density ratios obtained per velocity interval of 0.5 km s−1. The filled green circle shows the total f1 and f2 ratios computed over the entire observational sample. The black curves display the theoretical ratios obtained for homogeneous gas with kinetic temperatures of 10, 30, 100, and 1000 K and thermal pressures varying, along each curve, between 102 and 107 K cm−3 (from the bottom left to the top right). The empty white circles indicate integer values of the logarithm of the thermal pressure (see Fig. 3).

Fig. 3 Schematic view of the decomposition method proposed by JT11. The f1 and f2 column density ratios observed along a given line of sight (top panel) are supposed to result from the linear combination of a fraction, glow, of neutral carbon at low pressure, and a fraction, ghigh, of neutral carbon at high pressure. The thermal pressure of the high-pressure components is supposed to be known (106 K cm−3 in this example). glow and ghigh are derived from the geometric position of the observed point along the line that connects the low-pressure and high-pressure components (dashed black line in the top panel). This decomposition procedure is applied to all observations (bottom panel) to derive the distribution of thermal pressure of the low-pressure component.

2.3 Decomposition method of Jenkins & Tripp

Through a detailed analysis, JT11 demonstrated (see Fig. 3 of their paper) that the excitation conditions displayed in Fig. 2 cannot result from a lognormal distribution of thermal pressure expected in a turbulent environment (e.g. Kritsuk et al. 2007) or from a power-law tail of thermal pressure originating from a polytropic gas with an exponent smaller than one (Passot & Vázquez-Semadeni 1998) or self-gravitating or collapsing clouds (Federrath 2013; Girichidis et al. 2014). The reason is that the curvature of the line defined by the f₁ and f₂ combinations as a function of the gas thermal pressure (see Fig. 2) is insufficient to raise the f₂ ratio to the observed values if a lognormal or a power-law distributions of pressure are assumed. They concluded that the most straightforward interpretation of the data is to consider that any line of sight results from a bimodal distribution where most of the gas lies at thermal pressures between ~10³ and 10⁴ K cm⁻³ and a small fraction of the mass at thermal pressures one to three orders of magnitude above.

Expanding on this conclusion, JT11 applied the decomposition procedure schematised in Fig. 3. They assumed that the f₁ and f₂ ratios of the high-pressure gas are known (and correspond, for instance, to the values derived for T = 1000 K and P = 10⁶ K cm⁻³, see Fig. 3). If the kinetic temperature of the low-pressure environment is known (the value used in the theoretical model is fixed in Fig. 3 for pedagogical purposes), any observed line of sight (yellow point on the top panel of Fig. 3) can be seen as as a linear combination of a fraction, g_low, of neutral carbon at low pressure and a fraction, g_high, of neutral carbon at high pressure. JT11 estimated the kinetic temperature of the gas at low pressure from measurements of the rotational temperature of the two first levels of molecular hydrogen (see Table 2 in their paper). When this temperature was not measured, they assumed a canonic value of 80 K close to the mean kinetic temperature of the CNM. Applying this methodology to the entire sample (bottom panel of Fig. 3), JT11 derived the g_low and g_high fractions along every line of sight and the distribution of thermal pressure of the low-pressure environment. In doing so, they not only found that ~95% of the mass of the gas lies at low pressure with a lognormal distribution centred at P ~ 3.8 × 10³ K cm⁻³, but they also demonstrated that these findings weakly depends on the exact f₁ and f₂ ratios assumed for the high-pressure component.

Despite these successes, the observations of the excited levels of neutral carbon still raise several conundrums. Although the exact f₁ and f₂ ratios of the high-pressure component weakly affects the decomposition, it would be interesting to know what the actual physical conditions of this gas are. More importantly, these observations raise the question of the origin of the high-pressure component: why the ISM should contain very small amounts of gas (~5%) at anomalously large pressures and how it is possible that such a component, seen on every line of sight, co-exists with CNM environments at much lower thermal pressure.

2.4 Line profiles

The observed spectra also contain an important piece of information: the velocity structure of the high-pressure component. Using the highest-resolution configuration of the STIS instrument aboard the Herschel Space Telescope, JT11 derived the column density spectra of each level with a resolution of 0.5 km s⁻¹. Since the high-pressure component is characterised by large values of the f₂ ratio, these spectra can be used to put constraints on the kinematic structure of the associated gas. Analysis of the excitation conditions obtained in each channel of the observed spectra shows that, in most cases, the velocity components of the gas at high pressure are blended with those of the gas at low pressure. There are, however, a few exceptions.

Figure 4 displays the column density spectrum of neutral carbon obtained towards HD 210839. While the velocity components above −25 km s⁻¹ are found to have f₂ ratios below 0.1 and are thus mostly associated with the low-pressure component, the two velocity components below −25 km s⁻¹ are characterised by f₂ ratios larger than 0.3 and are therefore entirely composed of gas at high thermal pressure. Expanding on this example, we find 14 lines of sight where the high-pressure gas appears as a separated velocity component. Those are identified as isolated velocity components where the f₂ column density ratio is larger than 0.13 over the entire component; that is, for each velocity channel. The kinematic information are extracted by performing multi-Gaussian fits of these 14 spectra using the curve fit optimisation algorithm provided by the SciPy python library. The names of these 14 lines of sight, the centroid velocities, the full width at half maximum (FWHM), and the column density of neutral carbon associated with the velocity components at high pressure are given in Table 1. This table shows that the gas at high pressure displays Gaussian velocity profiles with a mean FWHM of 3.5 km s⁻¹ and a dispersion of 1.1 km s⁻¹. Interestingly the column density of carbon contained in each velocity profile covers a wide range of values from 5.5 × 10¹² to 2.5 × 10¹⁴ cm⁻².

In this paper, we explore the possibility that the high-pressure component originates from shocks propagating in the WNM. The question is whether interstellar shocks propagating in neutral and diffuse gas can produce neutral carbon in sufficient quantity, with excitation and kinematic properties that match those derived from the observations.

Fig. 4 Column density spectrum of neutral carbon observed towards HD 210839 (from Table 4 of JT11). The column densities are derived over velocity intervals of 0.5 km s−1. Each point along the spectrum is colour-coded according to the value of the f2 column density ratio to highlight velocity components at different thermal pressures.

3 Models of interstellar shocks

3.1 Physical ingredients

Interstellar shocks propagating in the WNM are modelled using the Paris-Durham shock code, a public multi-fluid model³ built to follow the dynamical, thermal, and chemical structures of shock waves, at steady-state, in a plane-parallel geometry (Godard et al. 2019). Designed to study low-velocity shocks (V_S ⩽ 30 km s⁻¹) propagating in irradiated environments, the code was recently improved by Paper I to treat shocks propagating at higher velocities, up to 500 km s⁻¹, which are known to generate UV, extreme UV, and X-ray photons that interact with the pre-shock medium and induce the formation of a radiative precursor (e.g. Raymond 1979; Hollenbach & McKee 1989; Sutherland & Dopita 2017). This latest version of the code includes the chemical evolution of multi-ionised species, the cooling induced by these species⁴, and an exact radiative transfer algorithm for line emission, built on the coupled escape probability formalism (Elitzur & Asensio Ramos 2006), which accurately follows the propagation and interactions of the photons produced by the shock with the shocked and the pre-shocked gas.

The Paris-Durham shock code models a plane-parallel shock propagating at a velocity, V_S, in a magnetised 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. The preshock medium and the shocked gas are assumed to be irradiated by an isotropic UV radiation field set to the standard interstellar radiation field of Mathis et al. (1983) and scaled with a parameter, G₀, and to be pervaded by cosmic ray particles with an equivalent H₂ ionisation rate, ${\zeta }_{{\mathrm{H}}_2}$. Within this geometry and in the frame of reference of the shock front, the code solves the out-of-equilibrium dynamical, thermal, and chemical evolution of a fluid particle during its trajectory from the pre-shock to the post-shock.

The chemical evolution of the abundance of neutral carbon results from a chemical network involving 209 species and 3340 chemical reactions (see Paper I) that include, in particular, photoionisation processes (Heays et al. 2017), radiative and dielectronic recombinations (Badnell et al. 2003; Badnell 2006), and recombinations onto grains and polycyclic aromatic hydrocarbons (PAHs; Draine & Sutin 1987). The populations of the three fine structure levels of C are computed, at each point of the trajectory, taking into account excitation and de-excitation processes by radiative decay, inelastic collisions, radiative pumping, and chemical pumping. Inelastic collisions are treated including the collisional rates of neutral carbon with H (Launay & Roueff 1977), H₂ in its ortho and para states (Schroder et al. 1991), He (Staemmler & Flower 1991), H⁺ (Roueff & Le Bourlot 1990), and e⁻ (Johnson et al. 1987). Excitation rates of the fine structure levels by radiative pumping of the electronic states of carbon followed by fluorescence are computed using the optical pumping rates of JT11 (see Table 6 in their paper). The probability of exciting the neutral carbon during its chemical formation is finally calculated by assuming that the three fine structure levels are equally likely to be populated during the recombination of C⁺, the main formation pathway of C in the models explored here.

3.2 Standard set-up and grid of models

The range of parameters explored in this work and the values adopted in the standard model are given in Table 2. As a prototypical case, we consider a shock propagating at a velocity, V_S = 80 km s⁻¹, in a medium with physical conditions close to those of the WNM observed in the local ISM. The pre-shock gas is assumed to have a proton density of $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$, a transverse magnetic field strength of B₀ = 1 μG (Frick et al. 2001; Crutcher et al. 2010), and to be irradiated by the standard UV radiation field of Mathis et al. (1983) (G₀ = 1). The total cosmic ray ionisation rate of H₂ is set to 3 × 10⁻¹⁶ s⁻¹ (e.g. Indriolo et al. 2015; Neufeld & Wolfire 2017). To explore the impact of various types of shocks on the production and excitation of neutral carbon, we also consider a grid of 325 models, with an initial pre-shock density varying between 0.1 and 2 cm⁻³, a transverse magnetic field varying between 0.1 and 10 μG, and a shock velocity ranging between 10 and 200 km s⁻¹. The grid of models explored here is similar to the grid presented in Paper I, except for the shock velocity, which is limited to 200 km s⁻¹. The reason is that the dynamical timescales of SNRs with terminal velocities above 200 km s⁻¹ are shorter than the shock cooling times (see Sec. 5 for more details). This implies that shocks above 200 km s⁻¹ driven by SNRs should not be considered at steady-state.

The thermodynamical structure of the standard model is summarised in Fig. 5, which displays the trajectory of a shocked fluid particle in a proton density–thermal pressure diagram. As is described in Paper I, a fluid particle crossing a shock front undergoes three successive regimes. For strong magnetohydrodynamic shocks – that is, shocks where the sonic and the Alfvénic Mach numbers are much higher than one – the fluid particle is first heated adiabatically (red curve in Fig. 5) to a thermal pressure that scales as $$P=8.1\times {10}^5\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}}}{80\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right)}^2.$$(5)

To ensure an equilibrium between the post-shock thermal pressure and the pre-shock ram pressure, the gas then cools down isobarically (green curve in Fig. 5). Because the magnetic field is frozen in the ionised fluid, the strength of the magnetic field increases. When the magnetic pressure becomes comparable to the post-shock thermal pressure, the gas progressively shifts from an isobaric to an isochoric cooling (blue curve in Fig. 5). The final density of the gas, $n_{\mathrm{H}}^{\mathrm{f}}$, is therefore given by the equilibrium between the post-shock magnetic pressure and the pre-shock ram pressure: $$n_{\mathrm{H}}^{\mathrm{f}}=61{\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}}}{80\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right){\left(\frac{B_0}{1\mu \mathrm{G}}\right)}^{-1}.$$(6)

This simple description shows that interstellar shocks with velocities larger than 30 km s⁻¹ bring the gas to thermal pressure larger than 10⁵ K cm⁻³, and that the gas at high pressure reaches densities comparable to the typical densities of the CNM. Both features are of paramount importance for the production and the excitation of neutral carbon.

It should be noted that this description only holds for strong J-type shocks and cannot be blindly applied to C-type and CJ-type shocks, or even weak J-type shocks where the built up of magnetic pressure during the adiabatic jump suppresses the iso-baric cooling stage (see Fig. 8 of Paper I). Keeping this warning in mind, Eq. (5) is found to accurately describe a large majority of the models explored here, while Eq. (6) holds for almost all models.

Fig. 5 Trajectory of a shocked fluid particle from the ambient medium to the post-shock region obtained in the standard model. The trajectory is displayed in a proton density–thermal pressure diagram to highlight the phase transition from the WNM to the CNM induced by a shock at 80 km s−1. A fluid particle initially evolves through the radiative precursor (pink). 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 post-shock gas. The black curve indicates the thermal equilibrium state of the diffuse gas obtained for G0 = 1 (see Fig. 1 of Paper I). Light grey lines are isothermal contours from 1 to 106 K (from bottom right to top left). Note that these isocontours are not evenly spaced due to the change in the number of particles in a gas at high temperature (see Fig. 5 of Paper I).

3.3 Shock size and timescale

An important issue for the comparison with observations is the calculation of the extent of the shock. By construction, one-dimensional (1D) stationary models postulate that the gas is confined and cannot escape in the directions parallel to the shock front. As a result, all the shocked gas cools down homogeneously towards thermal equilibrium and remains on this state indefinitely. Such a prescription raises two problems. First, it prevents the development of thermal and dynamical instabilities, which are known to occur in the wake of shocks propagating in multiphase environments (e.g. Falle et al. 2020; Raymond et al. 2020; Kupilas et al. 2021; Markwick et al. 2021). Second, it implies that all integrated quantities, such as column densities or line intensities, depend on the extent of the post-shock gas (e.g. Gusdorf et al. 2008; Leurini et al. 2014).

Since the idealised trajectory computed by the model becomes, at some point, unrealistic, it is necessary to apply a criterion to cut the contribution of the infinite post-shock material. Several criteria have been proposed in the past, based on the abundances and temperature profiles or on energetic considerations (e.g. Wardle 1999; Melnick & Kaufman 2015; Godard et al. 2019). In this work, we propose a new criterion based on the physical evolution of the gas. As was described in the previous section, shocks propagating in the WNM produce a post-shock medium that eventually becomes dominated by magnetic pressure. When it happens, the gas effectively hits a magnetic wall and loses the freedom to be compressed in the direction of propagation of the shock. It is therefore likely that flows along the directions parallel to the shock front can no longer be neglected and that the subsequent isochoric evolution predicted by the model (see Fig. 5) is somehow artificial. In this work, we thus assume that the model is reliable and can be integrated until the density reaches a fraction, χ = 0.95, of the final density, $n_{\mathrm{H}}^{\mathrm{f}}$. Such a criterion is found to prevent the integration of an unrealistic amount of post-shock material and to provide reliable results on the column densities and the width of the line profiles, which vary by less than a factor of four and less than a factor of two, respectively, for χ, varying between 0.8 and 0.99.

Within this prescription, and for strong J-type shocks, the crossing time is equivalent to the isobaric cooling time. Indeed, for most of the models explored in this work, the crossing time is found to weakly depend on the strength of the magnetic field as long as B₀ ⩽ 3 μG and to scale (for 10 ⩽ V_S ⩽ 200 km s⁻¹) as $$t_{\mathrm{c}}={10}^4\mathrm{y}\mathrm{r}{\left(\frac{n_{\mathrm{H}}^0}{1{\mathrm{c}\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{V_{\mathrm{S}}}{100\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right)}^{-3}\text{ for }{V}_{\mathrm{S}}⩽100\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$$(7)

and $$t_{\mathrm{c}}={10}^4\mathrm{y}\mathrm{r}{\left(\frac{n_{\mathrm{H}}^0}{1{\mathrm{c}\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{V_{\mathrm{S}}}{100\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right)}^2\text{ for }{V}_{\mathrm{S}}>100\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1},$$(8)

as is expected from the dependence of the cooling rate on the gas kinetic temperature. This crossing time corresponds to a shock size that ranges between ~10⁻³ pc and ~1 pc over the entire grid of models.

Fig. 6 Chemical and excitation profiles of neutral carbon across the standard model. Top panel: relative abundance of neutral carbon as a function of the distance (bottom axis) or time (top axis) in the shock. The curve is colour-coded according to the value of the thermal pressure along the trajectory. Bottom panel: populations of the three fine structure levels of C as a function of the distance (bottom axis) or time (top axis) in the shock. The vertical grey line in both panels indicates the distance and time at which the gas reaches a fraction, χ = 0.95, of the final density, $n_{\mathrm{H}}^{\mathrm{f}}$.

3.4 Chemistry and excitation of CI

The evolutions of the abundance of neutral carbon and of the populations of its three fine structure levels obtained in the standard model are shown in Fig. 6. As the gas is converted from the WNM to the CNM, the compression induced by the shock wave leads to an increase in the recombination of C⁺, and hence in the relative abundance of C. Unexpectedly, shocks induce an overproduction of neutral carbon, compared to the relative abundances found in the WNM and the CNM, which occurs before the thermal pressure decreases. This fundamental feature is explained by comparing the chemical timescales to the shock cooling time. Because interstellar shocks partly or even fully ionise neutral hydrogen, the recombination of C⁺ is dominated by radiative and dielectronic recombinations⁵ with a characteristic timescale (Badnell et al. 2003; Badnell 2006) $$t_{\mathrm{r}}\sim 4\times {10}^2\mathrm{y}\mathrm{r}{\left(\frac{n_{\mathrm{H}}^0}{1{\mathrm{c}\mathrm{m}}^{-3}}\right)}^{3/2}\left(\frac{V_{\mathrm{S}}}{80\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}\right){\left(\frac{B_0}{1\text{ μ}\mathrm{G}}\right)}^{-1}.$$(9)

The destruction of C, in the part of the trajectory that contributes the most to the column density of neutral carbon, is dominated by photoionisation with a characteristic timescale (Heays et al. 2017): $$t_{\gamma }\sim {10}^2\mathrm{y}\mathrm{r}{\left(\frac{G_0}{1}\right)}^{-1}.$$(10)

These two scalings (Eqs. (9) and (10)) imply that the time required for neutral carbon to reach chemical equilibrium is always at least 100 times smaller than the isobaric cooling time (see Eqs. (7) and (8)). It follows that interstellar shock waves generate large column densities of neutral carbon with excitation conditions typical of a gas at high thermal pressure.

To illustrate this result, we display in Figs. 7 and 8 the total column densities of neutral carbon in shocks for different velocities and transverse magnetic field and the f₁ and f₂ column density ratios obtained in shocks at different velocities. In both figures, the column densities are calculated from the integration of the local densities, in the direction perpendicular to the shock front, up to the point where the proton density reaches a fraction χ = 0.95 of its final value $n_{\mathrm{H}}^{\mathrm{f}}$ (see, for instance, the grey vertical line in Fig. 6). Figure 7 shows that interstellar shocks propagating in the WNM produce column densities of neutral carbon that range between a few 10¹² and 10¹⁴ cm⁻² for B₀ ⩽ 3 μG. This range of values obtained in single shocks is in remarkable agreement with the column densities of neutral carbon derived along the 14 lines of sight where the gas at high thermal pressure is identified as a separated velocity component (see Table 1).

As is shown in Fig. 8, the excitation conditions of C in interstellar shocks trace the post-shock thermal pressure and are therefore highly dependent on the shock velocity (see Eq. (5)). For the models displayed here, the f₁ and f₂ column density ratios are found to cover a large range of values with (f₁, f₂) combinations of (0.18, 0.04) and (0.38, 0.5) for V_S = 20 and 150 km s⁻¹, respectively. When compared with Fig. 3, the combinations of f₁ and f₂ ratios shown in Fig. 8 are found to be either those required to explain individual lines of sight or those required to interpret the lines of sight as combinations of environments at low and high thermal pressure. All these results indicate that shocks propagating in the WNM could be a natural solution to the presence of neutral carbon at high pressure in the local diffuse ISM, provided that they produce kinematic signatures similar to the observations.

Fig. 7 Coloured table of the column densities of neutral carbon (in cm−2) computed across shocks for different velocities and different strengths of the pre-shock magnetic field. The column densities are integrated in the direction perpendicular to the shock front over gas with a proton density, $n_{\mathrm{H}}⩽0.95n_{\mathrm{H}}^{\mathrm{f}}$. All other parameters are set to their standard values (see Table 2). Numbers in parenthesis are powers of ten.

Fig. 8 f1 and f2 column density ratios computed across shocks at different velocities. The column densities are integrated in the direction perpendicular to the shock front over gas with a proton density, $n_{\mathrm{H}}⩽0.95n_{\mathrm{H}}^{\mathrm{f}}$. The colour of each point indicates the velocity of the shock. All other parameters are set to their standard values (see Table 2).

3.5 Line profiles

The kinematic signature of shocks are explored by calculating synthetic spectra of the column density of neutral carbon in its second excited level across shocks in the direction parallel to the direction of propagation. The velocity dependent column density spectrum is $$N\left({\mathrm{C}}_{{}^3{\mathrm{P}}_2}\right)(v)={\int }_0^{L}n\left({\mathrm{C}}_{{}^3{\mathrm{P}}_2}\right)\varphi (v)\mathrm{d}z,$$(11)

where $n\left({\mathrm{C}}_{{}^3{\mathrm{P}}_2}\right)$ is the local density of neutral carbon in the second excited level, ϕ(υ) is the local line profile, z is the direction of propagation of the shock, and the integral is performed up to the distance, L, at which the local proton density, $n_{\mathrm{H}}=0.95n_{\mathrm{H}}^{\mathrm{f}}$. The local line profile is assumed to be Gaussian and is computed as $$\varphi (v)=\frac{1}{\sqrt{2\pi }{\sigma }_v}\exp (-\frac{1}{2}{\left[\frac{v-v_0}{{\sigma }_v}\right]}^2),$$(12)

where υ₀ is the local velocity of the gas along the direction of propagation of the shock and σ_υ is the local 1D velocity dispersion, $${\sigma }_v={(\frac{kT}{m}+{\sigma }_{\text{tur }}^2)}^{1/2},$$(13)

set by the local kinetic temperature, T, the mass of neutral carbon, m, and the micro-turbulent velocity dispersion, σ_tur. For simplicity, and because there is no information regarding the amplitude of the micro-turbulence generated in the post-shock gas, σ_tur is set to 0.

The calculations applied to all models show that the synthetic spectra of interstellar shocks display Gaussian line profiles. This is due to the fact that the column density of neutral carbon is mostly produced in the post-shock medium, in which the velocity of the gas is roughly constant. The line profile integrated over the shock structure is therefore dominated by the thermal broadening set by the evolution of the kinetic temperature during the isobaric cooling of the post-shock gas, with a weak, yet non-negligible, contribution of the velocity gradient. To quantify further this feature, the synthetic spectra are fitted with a Gaussian component using the curve fit optimisation algorithm provided by the SciPy python library. The result of this procedure is shown in Fig. 9, which displays the FWHM of the column density spectra of shocks for different velocities and transverse magnetic field.

As intuition dictates, the dependence of the predicted FWHM on the shock velocity, the transverse magnetic field, and the pre-shock density, is found to roughly follow the dependence expected from the square root of the ratio of the post-shock thermal pressure (Eq. (5)) and the final density (Eq. (6)), except for the dependence on B₀, which is somehow shallower. Keeping in mind that the predictions could be affected by the turbulent motions generated in the wake of interstellar shocks, the FWHM are found to range from 1 to 6 km s⁻¹, for a transverse magnetic field of B₀ ⩽ 3 μG, in remarkable agreement with the range of FWHM derived from the observations in velocity components solely associated with the high-pressure gas (see Table 1). This result, combined with those presented in the previous section, implies that shocks propagating in the WNM meet all the criteria required to explain the presence of neutral carbon at high pressure. It also clearly rule out shocks propagating in highly magnetised WNM (B₀ > 3 μG). The question now is whether this possible scenario is plausible or not.

Fig. 9 Coloured table of the FWHM (in km s−1) of the absorption lines originating from the 3P2 level of neutral carbon computed across shocks for different velocities and different strengths of the pre-shock magnetic field. The FWHMs are obtained by fitting Gaussian profiles on the velocity profiles predicted by the shock model in the direction perpendicular to the shock front taking only into account gas with a proton density $n_{\mathrm{H}}⩽0.95n_{\mathrm{H}}^{\mathrm{f}}$. All other parameters are set to their standard values (see Table 2).

4 Dissipation rate

4.1 Identification of individual shocks

The fact that the line profiles predicted by the model are Gaussian with line widths similar to those observed for individual velocity components at high pressure suggests that each of these velocity components is associated with a single shock. The number of velocity components clearly identified as gas at high pressure also suggests that each line of sight observed by JT11 intercepts only a few shocks propagating in the WNM and not a collection of tens or hundreds of shocks. It would therefore be tempting to use the predictions given in the previous section to analyse the lines of sight observed by JT11 and identify the properties of individual shocks.

For instance, the column density spectrum of HD 210839 (see Fig. 4) displays two velocity components at high pressure, centred at V = −34 and −31 km s⁻¹, with FWHM of 2.6 and 2.9 km s⁻¹, and total column densities of neutral carbon of 7.1 × 10¹³ and 5 × 10¹³ cm⁻² (see Table 1). The results given in Figs. 7 and 9 indicate that the observations could be explained by two shocks propagating along the line of sight at velocities between 80 and 100 km s⁻¹ in a medium with a pre-shock density of 1 cm⁻³ and a transverse magnetic field between 0.1 and 0.3 μG. There are, however, two problems with this approach. The first problem is that such a solution is not unique. Although not shown in Figs. 7 and 9, equally acceptable solutions could be found with different sets of shock parameters (V_S, $n_{\mathrm{H}}^0$, and B₀). Moreover, this method depends on the angle of inclination of the shock along the line of sight, which affects the predicted column density, and on the amplitude of turbulent motions in the post-shock flow, which may affect the width of the line profiles, two unknown quantities. The second problem is that this method can only be applied on the 14 lines of sight where velocity components at high thermal pressure are identified and not on the entire sample of JT11.

The various degeneracies cited above implies that the sole observations of the excited levels of neutral carbon are insufficient to identify the exact properties of the shocks and of the medium in which they propagate. The search for additional observational tracers will be carried in forthcoming papers. Short of these additional tracers, we devise, in the following, a method to analyse every line of sight and derive a more meaningful and robust quantity: the distribution of the dissipation rate of mechanical energy required to explain the observations.

4.2 Decomposition method

To illustrate the method, we considered first a simple case and assume that all the shocks propagating in the WNM have the same physical properties: the same pre-shock density, $n_{\mathrm{H}}^0$, shock velocity, V_S, and transverse magnetic field, B₀. If those parameters are fixed, the f₁ and f₂ column density ratios integrated across this shock are known and given by the model (see Fig. 8). If the gas at high thermal pressure corresponds to shocked material, the decomposition procedure of JT11 (see Sec. 2.3 and Fig. 3) can be applied adopting the (f₁, f₂) combination predicted by the model for the gas at high thermal pressure instead of an arbitrary combination (see Sec. 2.3). This decomposition procedure provides the fraction of gas at high thermal pressure and the associated column density of neutral carbon, $N_{\text{high }}^{\mathrm{o}}$(C), along any given line of sight. Therefore, the number of shocks required to explain each observation is $${\mathcal{N}}_{\text{shk }}=\frac{N_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}^{\mathrm{o}}(\mathrm{C})}{N^{\mathrm{m}}(\mathrm{C})},$$(14)

where N^m(C) is the column density of carbon predicted by the model (see Fig. 7). It follows that the rate of dissipation of mechanical energy (in erg cm⁻³ s⁻¹), averaged over the length of the neutral gas, is given by $${\epsilon }_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{s}}=\frac{1}{2}\mu n_{\mathrm{H}}^0{V}_{\mathrm{S}}^3\frac{{\mathcal{N}}_{\mathrm{s}\mathrm{h}\mathrm{k}}}{l_{\mathrm{l}\mathrm{o}\mathrm{s}}\phi },$$(15)

where μ is the mean particle mass of the atomic gas (in g) and ϕ is the fraction of volume occupied by the ionised gas along the line of sight set to the conservative value of 0.5 (see Sect. 5.3 of Bellomi et al. 2020). l_los is the length of diffuse gas (ionised or not) intercepted by the line of sight and estimated with Eq. (1).

The result of this methodology is shown in Fig. 10 which displays the probability distribution function of the logarithm of the dissipation rate obtained after applying the decomposition procedure to the entire observational sample, assuming that all the shocks are identical and correspond to the standard model ($n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$, B₀ = 1 μG, and V_s = 80 km s⁻¹). Under this strong assumption, the dissipation rate is found to cover more than three orders of magnitude, with a median of 3 × 10⁻²⁵ erg cm⁻³ s⁻¹, and a dispersion of about 0.4 dex (corresponding to a factor of 2.5). This distribution weakly depends on the prescription adopted in Sec. 2.1 to compute the length of the intercepted diffuse material (see Eq. (1)). This is because most of the lines of sight are located at low Galactic latitude where l_los is close or equal to the distance of the background source (see Fig. 1).

Fig. 10 Probability distribution function of the dissipation rate required to reproduce the column density of carbon at high pressure along the N lines of sight of JT11, assuming that all the shocks intercepted by these lines of sight are identical and defined by the parameters of the standard model: $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$, B0 = 1 μG, and VS = 80 km s−1. The probability distribution function is normalised so that its integral is equal to one.

4.3 Impact of the parameters

Although interesting, the distribution of dissipation rates obtained in Fig. 10 has a limited value because it is derived assuming that all the shocks propagating in the WNM are identical. Such an assumption is obviously not realistic. Even if the density of the WNM in the local ISM is known to follow a log-normal distribution around a mean value of about 0.7 cm⁻³ (Marchal & Miville-Deschênes 2021), shocks are expected to occur in a variety of warm neutral environments with different densities and magnetisation. Moreover, the probable driving sources of interstellar shocks (see Sec. 5) imply that the observed lines of sight necessarily intercept shocks with different velocities. These various distributions of physical conditions need to be taken into account to derive a meaningful distribution of dissipation rates.

The impacts of the model parameters on the decomposition procedure described in the previous section are shown in Fig. 11 which displays the distributions of the logarithm of the dissipation rate obtained for various sets of pre-shock density, transverse magnetic field, and shock velocity. As is shown in Fig. 8, the (f₁, f₂) combination highly depends on the shock velocity. In particular, shocks at low velocity lead to f₁ and f₂ ratios which are not sufficient to perform the decomposition of all the lines of sight. For $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$ and B₀ = 1 μG, we find that shocks with velocities larger than 30–40 km s⁻¹ are required to allow the decomposition of most of the observational sample (see Figs. 2 and 8). Each distribution displayed in Fig. 11 is thus computed and normalised taking only into account the lines of sight which can be decomposed for the associated set of shock parameters (V_S, $n_{\mathrm{H}}^0$, and B₀).

Figure 11 reveals two fundamental and surprising results. Firstly, the logarithm of the dissipation rate almost systematically follows a Gaussian distribution (plus a tail at low energy) with a dispersion of ~0.4, regardless of the model parameter. Secondly, most of the distributions shown in Fig. 11 appear to be superimposed. This last result is quantified further in Fig. 12, which displays the median value of these distributions for $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$. As long as B₀ ⩽ 3 μG, the median dissipation rate is found to weakly depend on the model parameters and to vary by only a factor of ~10 for input fluxes of kinetic and magnetic energies varying by a factor of ~10 000. For shock velocities between ~40 and ~150 km s⁻¹ (turquoise and blue regions in Fig. 12), the median dissipation rate obtained for $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$ varies by a factor of ~6, which corresponds to only twice the dispersion of the distribution.

This weak dependence of the dissipation rate on the model parameters has no simple analytical form and basically results from two antagonistic effects. On the one side, larger velocities and pre-shock densities lead to higher column densities of neutral carbon, N^m(C), with higher f₂ ratios, and hence lower $N_{\text{high }}^{\mathrm{o}}$(C). The number of shocks required to explain a given line of sight therefore drastically decreases when the velocity or the density increases (see Eq. (14)). On the other side, shocks with larger velocities or pre-shock densities have larger individual dissipation rates. These antagonistic effects compensate each other almost perfectly, leading to a stable distribution of dissipation rates over a wide range of parameters.

The result displayed in Figs. 11 and 12 is one of the main outcome of this work. It implies that the dissipation rate required to explain the observations of JT11 and the presence of carbon at high pressure is roughly independent of the exact nature of the shocks along the lines of sight and the medium in which they propagate. It also implies that the distribution displayed in Fig. 10 is more meaningful than originally thought because a similar distribution would be obtained assuming a convoluted distribution of shock velocities, pre-shock densities, and transverse magnetic field. It finally means that the validity of this scenario can be estimated by considering the possible driving sources of interstellar shocks in the WNM and compare the associated distributions of dissipation rates with Fig. 10.

Fig. 11 Same as Fig. 10 but for a pre-shock density varying between 0.1 and 2 cm−3 (from left to right), a strength of pre-shock magnetic field varying between 0.1 and 10 μG (from top to bottom), and six values of shock velocity: 20 (red), 40 (orange), 60 (green), 80 (turquoise), 100 (light blue), and 120 (dark blue) km s−1. For coherence, the colour code is the same than that used in Fig. 8. In many cases, only the last distribution (corresponding to VS = 120 km s−1) is visible, because those obtained at lower velocities are identical and hidden behind it. The bottom left panels are empty because the distributions fall outside the range of dissipation rates displayed.

5 A scenario of shocks driven by supernovae

5.1 Origin of high velocity shocks

As was shown in the previous sections, the presence of neutral carbon at high pressure could be explained by shocks propagating at velocities V_S ≳ 30 km s⁻¹ in a WNM environment with a typical density, $n_{\mathrm{H}}^0⩾0.3{\mathrm{c}\mathrm{m}}^{-3}$, and a magnetic field strength, B₀ ⩽ 3 μG. In theory, shocks with such velocities could be driven by a variety of phenomena in the local ISM including stellar winds, bipolar outflows, supernovae explosions, or the infall of matter onto the Galactic disc. The relative contributions of these events depend, however, on their individual rates, the deceleration time of the shock they produce, and the volume of the ISM they affect. Both analytical studies and numerical simulations show that supernovae explosions should dominate the injection rate of mechanical energy in the Milky Way (e.g. Norman & Ferrara 1996; Brucy et al. 2020). We therefore discuss here the scenario where shocks propagating in the WNM are driven by the expansion of SNRs in the diffuse phase of the local ISM.

5.2 Expansion of supernova remnants

To estimate the physical properties of shocks driven by SNRs, we consider the classical picture of a supernovae explosion with an ejecta mass, M_ej, and an initial kinetic energy, E, expanding in an homogeneous medium with a proton density $n_{\mathrm{H}}^0$. In an idealised spherical expansion, the resulting SNR is known to evolve through four different stages, designated as the free expansion stage (FE), the Sedov-Taylor stage (ST), the pressure-driven stage (PD), and the momentum conserving stage (MC) (e.g. Cioffi et al. 1988; Truelove & McKee 1999; Kim & Ostriker 2015). During each of these stages, the radius, R_b, and the velocity, V_b, of the terminal blast wave can be approximated by the following expressions (e.g. Draine 2011; Vink 2020) $$\begin{array}{rl}& R_B=\{\begin{array}{cl}3.1\mathrm{p}\mathrm{c}{\left(\frac{M_{\mathrm{e}\mathrm{j}}}{M_{\odot }}\right)}^{1/3}{\left(\frac{n_{\mathrm{H}}^0}{n_0}\right)}^{-1/3}\left(\frac{t}{t_{\mathrm{S}\mathrm{T}}}\right)& \text{ if }0<t⩽t_{\mathrm{S}\mathrm{T}}\\ 3.1\mathrm{p}\mathrm{c}{\left(\frac{M_{\mathrm{e}\mathrm{j}}}{M_{\odot }}\right)}^{1/3}{\left(\frac{n_{\mathrm{H}}^0}{n_0}\right)}^{-1/3}{\left(\frac{t}{t_{\mathrm{S}\mathrm{T}}}\right)}^{2/5}& \text{ if }{t}_{\mathrm{S}\mathrm{T}}<t⩽t_{\mathrm{P}\mathrm{D}}\\ 23.8\mathrm{p}\mathrm{c}{\left(\frac{E}{E_{51}}\right)}^{0.29}{\left(\frac{n_{\mathrm{H}}^0}{n_0}\right)}^{-0.42}{\left(\frac{t}{t_{\mathrm{P}\mathrm{D}}}\right)}^{2/7}& \text{ if }{t}_{\mathrm{P}\mathrm{D}}<t⩽t_{\mathrm{M}\mathrm{C}}\\ 46.0\mathrm{p}\mathrm{c}{\left(\frac{E}{E_{51}}\right)}^{0.29}{\left(\frac{n_{\mathrm{H}}^0}{n_0}\right)}^{-0.42}{\left(\frac{t}{t_{\mathrm{M}\mathrm{C}}}\right)}^{1/4}& \text{ if }{t}_{\mathrm{M}\mathrm{C}}<t⩽t_{\mathrm{f}\mathrm{a}\mathrm{d}\mathrm{e}}\end{array}\\ & \end{array}$$(16)

and $$V_B=\{\begin{array}{cl}\mathrm{10\; 000}\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}{\left(\frac{E}{E_{51}}\right)}^{1/2}{\left(\frac{M_{\mathrm{e}\mathrm{j}}}{M_{\odot }}\right)}^{-1/2}& \text{ if }0<t⩽t_{\mathrm{S}\mathrm{T}}\\ 3988\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}{\left(\frac{E}{E_{51}}\right)}^{1/2}{\left(\frac{M_{\mathrm{e}\mathrm{j}}}{M_{\odot }}\right)}^{-1/2}{\left(\frac{t}{t_{\mathrm{S}\mathrm{T}}}\right)}^{-3/5}& \text{ if }{t}_{\mathrm{S}\mathrm{T}}<t⩽t_{\mathrm{P}\mathrm{D}}\\ 135\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}{\left(\frac{E}{E_{51}}\right)}^{0.07}{\left(\frac{n_{\mathrm{H}}^0}{n_0}\right)}^{0.13}{\left(\frac{t}{t_{\mathrm{P}\mathrm{D}}}\right)}^{-5/7}& \text{ if }{t}_{\mathrm{P}\mathrm{D}}<t⩽t_{\mathrm{M}\mathrm{C}}\\ 26\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}{\left(\frac{E}{E_{51}}\right)}^{0.07}{\left(\frac{n_{\mathrm{H}}^0}{n_0}\right)}^{0.13}{\left(\frac{t}{t_{\mathrm{M}\mathrm{C}}}\right)}^{-3/4}& \text{ if }{t}_{\mathrm{M}\mathrm{C}}<t⩽t_{\mathrm{f}\mathrm{a}\mathrm{d}\mathrm{e}},\end{array}$$(17)

where E₅₁ = 10⁵¹ erg and n₀ = 1 cm⁻³ are normalisation factors, t_ST, t_PD, and t_MC are the onset times of the ST, PD, and MC stages, respectively, and t_fade is the fade away time of the supernova. The onset time of the Sedov-Taylor stage is set to ensure the continuity of the radius between the FE and the ST stages⁶, $$t_{\mathrm{S}\mathrm{T}}=0.303\mathrm{k}\mathrm{y}\mathrm{r}{\left(\frac{E}{E_{51}}\right)}^{-1/2}{\left(\frac{M_{\mathrm{e}\mathrm{j}}}{M_{\odot }}\right)}^{5/6}{\left(\frac{n_{\mathrm{H}}^0}{n_0}\right)}^{-1/3}.$$(18)

Following Draine (2011), the onset time of the pressure-driven stage is set to the moment when the integrated radiative losses of the hot gas become comparable to its thermal energy (ΔE_th/E_th = −1/3), $$t_{\mathrm{P}\mathrm{D}}=49.3\mathrm{k}\mathrm{y}\mathrm{r}{\left(\frac{E}{E_{51}}\right)}^{0.22}{\left(\frac{n_{\mathrm{H}}^0}{n_0}\right)}^{-0.55}.$$(19)

The onset of the momentum conserving stage is hard to estimate (Kim & Ostriker 2015). To simplify we set it here to 10 times the onset time of the PD stage. The fade away time is finally taken as the moment where the terminal velocity falls below the 1D velocity dispersion of the WNM, set to 10 km s⁻¹, $$t_{\text{fade }}=1.76\mathrm{M}\mathrm{y}\mathrm{r}{\left(\frac{E}{E_{51}}\right)}^{0.31}{\left(\frac{n_{\mathrm{H}}^0}{n_0}\right)}^{-0.38}.$$(20)

Figure 13 displays the evolutions of the blast radius and velocity considering the prototypical case of a supernova explosion with M_ej = 1.4 M_⊙, E = 10⁵¹ erg, and $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$ (Draine 2011). The discontinuities of the terminal velocity come from the idealised descriptions adopted above which neglect the deceleration of the blast wave by the surrounding environment at the end of the FE stage and the cooling of the hot shell at the end of the ST stage. Figure 13 shows that a typical SNR produces a terminal shock with a velocity above ~30 km s⁻¹ over long periods that last from a few tens to a few hundreds kyr and is therefore a likely source of intermediate and high velocity shocks in the WNM. Interestingly, the terminal velocity at the end of the Sedov-Taylor stage is ~188 km s⁻¹ with small dependences on the model parameters. Since the blast deceleration timescale is smaller than the cooling timescale of the hot gas in the ST stage, this result implies that shocks driven by SNRs at velocities larger than ~200 km s⁻¹ are not radiative (Vink 2012) and cannot be treated at steady-state.

Fig. 12 Coloured table of the median value of the dissipation rate (see Fig. 11) as a function of the shock velocity and the strength of the magnetic field for $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$. All the other parameters are set to their standard values (see Table 2). Models with a median dissipation rate within a factor of two of that obtained with the standard model are shown in turquoise. Models with a median dissipation rate between two and four times larger (resp. lower) than that obtained with the standard model are shown in green (resp. blue). All other models are shown in red.

Fig. 13 Basic analytical descriptions of the evolutions of the radius (blue curve) and the terminal velocity (red curve) of a 1D SNR (e.g. Draine 2011; Vink 2020), characterised by an ejecta mass of Mej = 1.4 M⊙ and an initial kinetic energy energy of E = 1051 erg, and expanding in a homogeneous medium with a density of 1 cm−3. The idealised successive phases known as the free expansion stage, the Sedov-Taylor stage, the pressure-driven stage, and the momentum conserving (or snowplow) stage are highlighted in red, orange, green, and blue, respectively.

5.3 Dissipation rate of a simplistic distribution of SNR

Estimating the dissipation rate of mechanical energy induced by SNRs along the lines of sight observed by JT11 requires a description of the spatial distribution of supernovae in the Galaxy. To perform a back of the envelope calculation of the distribution of the dissipation rate, we adopt the following simplified picture. All supernovae are assumed to have the same ejecta mass, M_ej = 1.4 M_⊙, and initial kinetic energy, E = 10⁵¹ erg, and to expand in a medium defined by a single proton density, $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$ (see Fig. 13). The rate of supernova explosions in the Galaxy, k_SN, is set to 46 kyr⁻¹ (Adams et al. 2013). In a permanent regime, the number of SNRs in the Galaxy therefore is $${\mathcal{N}}_{\mathrm{S}\mathrm{N}\mathrm{R}}=k_{\mathrm{S}\mathrm{N}}{t}_{\text{fade }},$$(21)

where t_fade is the SNR fade away time (see Eq. (20)). Based on the model of Adams et al. (2013), we assume that 70% of the supernovae correspond to core collapses (SNII, SNIb, and SNIc) and are distributed along the spiral arms, and that 30% of the supernovae are thermonuclear (SNIa) and are homogeneously distributed in the Galactic disc. The Galaxy is modelled as a thin disc with a radius of 12 kpc and a vertical Gaussian structure with a dispersion of 35 pc set from the observed scale heights of very young open clusters (Hao et al. 2021) and of O-B₅ stars (Maíz-Apellániz 2001). Following Vallée (2022), the spiral arms are described as four logarithmic spirals equally spaced in azimuth by 90°. The pitch angle is set to 13.4°. The starting point of the Sagittarius spiral arm is set at an angle of 50° and at a galactocentric distance of 2 kpc. The distribution of core-collapse supernovae in each spiral arm is finally assumed to follow a Gaussian distribution with a dispersion of 500 pc (Pohl & Esposito 1998).

This simple model is used to generate 3D random distributions of SNRs in the Milky Way with random ages homogeneously drawn between 0 and t_fade. An example of such a distribution is shown in Fig. 14 (panels a and b) which displays the positions, the sizes, and the shock velocities of all blast waves with terminal velocities between 30 and 200 km s⁻¹ in projection onto the Galactic disc. As is shown in Fig. 14 (panel c), any line of sight observed by JT11 (red lines in Fig. 14) intercepts either N = 0, 1, or 2 surfaces of a given SNR, depending on the 3D position and the radius of the SNR and the 3D position of the background star. The dissipation rate induced by shocks with velocities between 30 and 200 km s⁻¹ along this line of sight is therefore given by $${\epsilon }_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{s}}={\displaystyle \sum _{j=1}^{{\mathcal{N}}_{\mathrm{S}\mathrm{N}\mathrm{R}}}}\frac{1}{2}\mu n_{\mathrm{H}}^0{V}_{\mathrm{B},j}^3\frac{N_j}{l_{\mathrm{los}}\phi }\frac{1}{\cos \left(i_j\right)},$$(22)

where i is the inclination angle of an SNR surface compared to the plane perpendicular to the line of sight (see Fig. 14, panel c) and where the sum is performed over all SNRs with terminal velocities between 30 and 200 km s⁻¹.

Figure 15 (top panel) displays the distribution of dissipation rates expected from SNR surfaces along the lines of sight observed by JT11 averaged over 500 realisations of the model described above. Interestingly, we find that shocks driven by SNRs are an unavoidable scenario. Indeed, as is suggested by the projection view displayed in Fig. 14 and as is shown in Fig. 15, the lines of sight observed by JT11 are bound to intercept SNR surfaces, and thus interstellar environments at high thermal pressure. The dissipation rate induced by SNRs is found to follow a log-normal distribution with a median value of ~ 3 × 10⁻²⁵ erg cm⁻³ s⁻¹ and a dispersion of ~0.6 dex (corresponding to a factor of ~4). It is remarkable to see that this distribution, derived from a simplistic statistical model of SNRs, covers the same orders of magnitude than the distribution of dissipation rates required to explain the observations of the UV lines of CI (see Sec. 4). These results imply that shocks propagating in the WNM are not only a possible but also a probable scenario to account for the presence of gas at high thermal pressure detected in the local ISM. It also suggests that these shocks propagating at intermediate and high velocities are probably mainly driven by SNRs.

The fact that the two distributions shown in the top panel of Fig. 15 are remarkably similar should not be overinterpreted. On the one side, it should be kept in mind that the statistical model of SNRs results from multiple simplifications. It follows that the distribution of dissipation rate derived from this model is sensitive to several parameters including the density of the ambient medium, $n_{\mathrm{H}}^0$, the rate of supernova explosions in the Galaxy, and the geometrical parameters that describe the shape and extent of the spiral arms and of the Galactic disc. Moreover, this model builds an homogeneous distribution of core collapse supernovae along the spiral arms while supernovae are expected to explode in clusters. Since supernovae are probably spatially correlated with OB star clusters and since OB stars are the background sources used by JT11, such a simplification impacts the dissipation rate induced by SNR along the lines of sight observed by JT11. On the other side, the decomposition procedure applied in Sec. 4 to derive the dissipation rates from the observation of neutral carbon has several limitations. Indeed, even if the model parameters have a weak impact on the distribution of dissipation rates deduced from CI (see Sec. 4), adopting a distribution of shock velocities, pre-shock densities, and transverse magnetic field could broaden the distribution of dissipation rates up to a factor of 4 (see Fig. 12) and induce a shift of the median value. Moreover, the small curvature of the line defined by the f₁ and f₂ combinations at small thermal pressure and for low velocity shocks (see Figs. 2 and 8) implies that the decomposition procedure is not sensitive enough to accurately measure events at low dissipation rate. This implies that the left part of the turquoise distribution shown in Fig. 15 (for ϵ_diss ≲ 10⁻²⁵ erg cm⁻³ s⁻¹) is less reliable than the right part.

Fig. 14 Schematic view of the calculation of the dissipation rates expected from SNRs along the lines of sight observed by JT11. A distribution of SNRs with terminal velocities between 30 and 200 km s−1 is drawn out of a simple description of supernovae explosions within the disc and the spiral arms of the Galaxy (panel a). The SNRs are colour-coded according to their terminal velocities. The lines of sight observed by JT11 (red lines) are analysed to count the number of SNR surfaces they intercept and compute the dissipation rates induced by the associated shocks (panels b and c). This calculation is done taking into account the 3D spatial distribution of SNRs, which is shown in projection in panels a and b. Panel c displays a schematic view of a single SNR located at a position Aj. An observed line of sight connects the position of the sun, S, to the position of the background star, B. The number of surfaces crossed by the line of sight for this single SNR (0, 1, or 2) and their inclination angle, ij, depend on the position of the point, Mj, along the SB line and its distance to Aj compared to the SNR radius.

Fig. 15 Top panel: probability distribution function of the dissipation rate along the lines of sight observed by JT11 averaged over 500 realisations of a simple and unclustered distribution of SNRs in the Galaxy (see Fig. 14) compared to that derived from the production and the excitation of neutral carbon in the standard model of interstellar shocks (see Fig. 10). Bottom panel: probability distribution function of the kinetic energy transfer rate derived from the observations of CO in molecular clouds located in the local ISM (Hennebelle & Falgarone 2012). All distributions shown in the top and bottom panels are normalised so that their integral is equal to one.

5.4 Connection with the turbulent cascade

Supernovae explosions are often considered to be the main driving source of interstellar turbulence in the Milky Way. It is therefore interesting to discuss the results obtained in the previous sections in light of the properties of interstellar turbulence derived in the local ISM. One of the main properties of turbulence is the kinetic energy transfer rate of the turbulent cascade which measures the amount of energy per unit time and volume injected at large scale in turbulent motions and transferred to smaller and smaller scales where it is eventually dissipated (Hennebelle & Falgarone 2012; Miville-Deschênes et al. 2017). The distribution of the kinetic energy transfer rate deduced from the velocity dispersion of molecular clouds in the solar neighbourhood by Hennebelle & Falgarone (2012) is shown on the bottom panel of Fig. 15. In the local ISM, this transfer rate is found to vary over four orders of magnitude, with a median value ~10⁻²⁵ erg cm⁻³ s⁻¹ and a dispersion ~0.8 dex (corresponding to a factor of 6). Unexpectedly, these values are found to be similar to the median and dispersion values of the dissipation rate induced by SNRs.

It is important to realise that the top and bottom panels of Fig. 15 do not display the same quantities: the top panel shows the amount of energy allegedly dissipated in shocks driven by SNRs, while the bottom panel displays the amount of energy which is not dissipated in the original shocks and is allegedly injected in interstellar turbulence. The fact that these two distributions are similar may reveal a fundamental property of the mechanisms by which the mechanical energy is injected in the interstellar matter. Indeed, it suggests that the amount of energy dissipated in SNRs and the amount of energy used to feed interstellar turbulence are comparable. It has been proposed in the past that turbulence might be generated in the wake of shocks through the production of baroclynic vorticity (Vazquez-Semadeni et al. 1996; Elmegreen & Scalo 2004) and various thermodynamical instabilities. Studying the detailed connection between the distributions shown in the top and bottom panels of Fig. 15 requires dedicated 3D numerical simulations of supernova blast waves. This is out of the scope of the present paper.

6 Discussion

6.1 Radial velocities

This paper demonstrates that shocks propagating in the WNM are a viable scenario to explain the presence of gas at high thermal pressure detected by JT11. To promote this idea, we showed that shocks efficiently produce neutral carbon with column densities, excitation conditions, and line profiles in agreement with the observations and that the energy required to power these shocks corresponds to that expected from Galactic SNRs. It is important to stress, however, that this paper does not take full advantage of all the information contained in the observations.

One additional piece of information is the distribution of radial velocities of the gas at high thermal pressure and the relation between these velocities and the (f₁, f₂) ratios. Table 1 shows that the 20 components identified as gas at high thermal pressure display a slight asymmetry of radial velocities. Indeed, eight of these 20 components have positive radial velocities up to 28 km s⁻¹, eleven have negative radial velocities down to −36 km s⁻¹, and one is detected at −83 km s⁻¹. Moreover the positive velocity components display, on average, lower (f₁, f₂) combinations compared to the negative velocity components. This last result was already shown by JT11, who found that environments with the highest thermal pressure, and hence the highest f₂ ratio, are systematically associated with velocity components at least 5 km s⁻¹ below the minimum value permitted by the differential Galactic rotation.

At first sight, this slight asymmetry would seem to contradict the simple scenario depicted in Fig. 14 of shocks propagating in random directions and expected to have roughly equal proportions of negative and positive velocities, and to favour a scenario where part of the high-pressure gas is produced through the interaction of the target stars with their immediate environments. It should be kept in mind, however, that Fig. 14 displays an idealised picture in which SNRs are distributed homogeneously along the spiral arms and the Galactic disc. Since supernovae are spatially correlated with OB star clusters, the OB stars used as background sources by JT11 may preferentially reside inside one or several SNRs rather than beyond. This effect would favour the occurrence of negative velocities corresponding to the expansion of the foreground shell, and explain the fact that the gas at highest thermal pressure is preferentially associated with negative velocity components.

Interestingly, Table 1 shows that the gas at high thermal pressure has moderate radial velocities compared to the shock velocity explored in this work. This is in line with previous observations of CI towards stars behind known SNRs which show that environments at high pressure are not always associated with significantly large displaced velocities (e.g. Jenkins et al. 1981, 1984; Raymond et al. 1991; Wallerstein et al. 1995; Nichols & Slavin 2004). This feature could, however, be coherent with shocks driven by SNRs. Indeed, in the framework of expanding SNRs, lines of sight are more likely to intercept the borders of the expanding bubble, and hence to have significant transverse motions, than the centre of the SNR where the gas has a prominent radial motion.

A complete modelling of the radial velocities and the asymmetry described above would require to take into account the dynamics of the Galaxy, including the rotation curve, the velocity dispersion and the vertical motions, and the clustering of supernovae compared to OB star associations. This is out of the scope of the present paper.

6.2 Limitations and alternative contributions

The distributions depicted in Fig. 15 suggest that shocks driven by SNRs are not only likely but also unavoidable contributors to the high-pressure gas observed in CI. However, this conclusion is moderated by the simplicity of the SNR model presented in Sec. 5, in particular the fact that the ambient medium in which SNRs expand is considered as homogeneous. Indeed, the ISM is highly inhomogeneous and contains CNM and dense clouds that cover a great variety of spatial scales and densities. As is shown in hydrodynamic, magnetohydrodynamic, adiabatic, and non-adiabatic numerical simulations, the presence of these clouds profoundly modifies the dynamics of the flow (e.g. Pittard et al. 2010 and references therein). The interaction of a shock with an interstellar cloud generates a forward and a reverse shocks in the cloud, a bow shock (or bow wave) upstream of the cloud structure, and reflected shocks and vorticity sheets downstream (e.g. Poludnenko et al. 2002; Yirak et al. 2010). The cloud itself is subject to Rayleigh-Taylor, Kelvin-Helmholtz, and thermal instabilities which lead to its disruption and fragmentation in the post-shock turbulent flow. It follows that the expansion of an SNR in a clumpy medium could be highly non spherical (see Figs. 2 and 3 of Korolev et al. 2015) and that the initial blast energy not only drives shocks in the dense clouds but is also partly converted into turbulent motions.

On the one side, the fact that part of the SNR energy is used to feed complex dynamical motions implies that the distribution of dissipation rates deduced from our simple SNR model (see Fig. 15) is an overestimation of the actual dissipation rate of shocks propagating in the WNM. On the other side, the propagation of low-velocity shocks in CNM structures induces an increase in the cloud thermal pressure on timescales smaller than its disruption timescale and could therefore contribute to the presence of high-pressure gas observed in CI. A complete modeling of all these aspects is of fundamental importance to understand the feedback of supernovae on the thermodynamical evolution of the multiphase ISM and will be carried out in future studies.

7 Conclusions

Observations of UV absorption lines of neutral carbon in the local ISM performed by JT11 revealed that a significant fraction of the mass of the ISM lies at thermal pressures one to three orders of magnitude above the thermal pressure of the bulk of the gas. In this paper, we propose that this enigmatic component originates from shocks propagating in the WNM. To explore this idea, we perform the first detailed and quantitative study of the production and excitation of neutral carbon in atomic shocks using the latest version of the Paris-Durham shock code presented in Paper I. The model is analysed over a wide range of parameters that cover the typical conditions of the WNM in the local ISM. We propose that these shocks are mainly driven by SNRs and discuss the validity of this scenario based on the expected distribution of the dissipation rate of mechanical energy. This analysis highlights the following results.

1. Shocks propagating in the WNM naturally generate high-pressure environments. The compression of the gas favours the production of neutral carbon through the increase in the radiative and dielectronic recombination of C⁺. Because the associated chemical timescales are much shorter than the cooling timescale, the production of neutral carbon occurs before the thermal pressure decreases.

2. The shocks explored in this work individually produce a column density of neutral carbon between 10¹² and 10¹⁴ cm⁻², distributed over a Gaussian line profile with a FWHM that ranges between 1 and 6 km s⁻¹. These results, which hold as long as the transverse magnetic field B₀ ⩽ 3 μG, are found to be in excellent agreement with the column densities and line widths observed in velocity components solely associated with gas at high thermal pressure.

3. Although the observed line profiles suggest that individual lines of sight intercept only a few shocks, the current set of observations is found to be insufficient to identify their exact physical properties, such as the shock velocities, the pre-shock densities, the strengths of the transverse magnetic field, and the inclination angles.

4. Each line of sight is decomposed into a high-pressure component and a low-pressure component to compute the number of shocks required to reproduce the observations and their dissipation rates. The column densities of neutral carbon observed at high thermal pressure require a dissipation rate of mechanical energy of ~3 × 10⁻²⁵ erg cm⁻³ s⁻¹ with a dispersion of about a factor of three. Surprisingly, this distribution of the dissipation rate weakly depends on the model parameters as long as V_S ⩾ 30 km s⁻¹, $n_{\mathrm{H}}^0$ ⩾ 0.3 cm⁻³, and B₀ ⩽ 3 μG. This result implies that the observations of the UV absorption lines of neutral carbon provide a direct measurement of the dissipation rate induced by shocks, regardless of the properties of these shocks and their distributions along the lines of sight.

5. These shocks are likely driven by supernovae explosions and the expansion of the associated SNRs in the diffuse ionised and neutral phases of the Galaxy. Indeed, a simplistic statistical distribution of SNRs located in the spiral arms and in a thin Galactic disc shows that the lines of sight observed by JT11 inevitably cross SNR surfaces. This model leads to an expected distribution of the dissipation rate with a median value of ~3 × 10⁻²⁵ erg cm⁻³ s⁻¹ and a dispersion of about a factor of four, in remarkable agreement with the distribution derived from the UV lines of CI.

6. Interestingly, the distribution of the dissipation rate induced by SNRs expanding in the WNM is found to be similar to the distribution of the kinetic energy transfer rate measured in molecular clouds located in the local ISM. This result implies that the amount of mechanical energy dissipated by shocks in SNRs and the amount of energy injected in the wake of shocks in interstellar turbulence are somehow comparable.

This work uncovers a potential direct tracer of shocks driven by SNRs in the nearby ISM. It also presents a methodology to estimate the rate of dissipation of kinetic energy in events that inject mechanical energy into the ISM, feed the interstellar turbulence, and induce phase transition between the WNM and the CNM. The fact that shocks propagating in the WNM have a distinctive imprint on the excitation properties of neutral carbon suggests that these shocks may produce additional observational tracers. Those could be used to identify the physical properties of the shocks or their distributions and, more generally, to study the injection of mechanical energy and the phase transition process in Galactic and extragalactic environments. The search for additional observational tracers of shocks propagating in the WNM will be carried out in forthcoming papers.

Acknowledgements

We are very grateful to the referee, Edward Jenkins, for his thorough reading of the manuscript and his valuable comments. We dedicate this work to his memory. 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.

References

All Tables

All Figures

	Fig. 1 Aitoff projection, in Galactic longitude and latitude coordinates, of the background sources of the observational sample of JT11. The colour-coding indicates the distance of background sources in parsecs.
In the text

Fig. 2 Excitation properties of the neutral carbon along the lines of sight observed by JT11. The filled coloured circles indicate the f1 and f2 column density ratios observed along each line of sight. Each circle is colour-coded according to the total column density of neutral carbon, N(C). The black dots show the f1 and f2 column density ratios obtained per velocity interval of 0.5 km s−1. The filled green circle shows the total f1 and f2 ratios computed over the entire observational sample. The black curves display the theoretical ratios obtained for homogeneous gas with kinetic temperatures of 10, 30, 100, and 1000 K and thermal pressures varying, along each curve, between 102 and 107 K cm−3 (from the bottom left to the top right). The empty white circles indicate integer values of the logarithm of the thermal pressure (see Fig. 3).

In the text

Fig. 3 Schematic view of the decomposition method proposed by JT11. The f1 and f2 column density ratios observed along a given line of sight (top panel) are supposed to result from the linear combination of a fraction, glow, of neutral carbon at low pressure, and a fraction, ghigh, of neutral carbon at high pressure. The thermal pressure of the high-pressure components is supposed to be known (106 K cm−3 in this example). glow and ghigh are derived from the geometric position of the observed point along the line that connects the low-pressure and high-pressure components (dashed black line in the top panel). This decomposition procedure is applied to all observations (bottom panel) to derive the distribution of thermal pressure of the low-pressure component.

In the text

	Fig. 4 Column density spectrum of neutral carbon observed towards HD 210839 (from Table 4 of JT11). The column densities are derived over velocity intervals of 0.5 km s−1. Each point along the spectrum is colour-coded according to the value of the f2 column density ratio to highlight velocity components at different thermal pressures.
In the text

Fig. 5 Trajectory of a shocked fluid particle from the ambient medium to the post-shock region obtained in the standard model. The trajectory is displayed in a proton density–thermal pressure diagram to highlight the phase transition from the WNM to the CNM induced by a shock at 80 km s−1. A fluid particle initially evolves through the radiative precursor (pink). 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 post-shock gas. The black curve indicates the thermal equilibrium state of the diffuse gas obtained for G0 = 1 (see Fig. 1 of Paper I). Light grey lines are isothermal contours from 1 to 106 K (from bottom right to top left). Note that these isocontours are not evenly spaced due to the change in the number of particles in a gas at high temperature (see Fig. 5 of Paper I).

In the text

Fig. 6 Chemical and excitation profiles of neutral carbon across the standard model. Top panel: relative abundance of neutral carbon as a function of the distance (bottom axis) or time (top axis) in the shock. The curve is colour-coded according to the value of the thermal pressure along the trajectory. Bottom panel: populations of the three fine structure levels of C as a function of the distance (bottom axis) or time (top axis) in the shock. The vertical grey line in both panels indicates the distance and time at which the gas reaches a fraction, χ = 0.95, of the final density, $n_{\mathrm{H}}^{\mathrm{f}}$.

In the text

Fig. 7 Coloured table of the column densities of neutral carbon (in cm−2) computed across shocks for different velocities and different strengths of the pre-shock magnetic field. The column densities are integrated in the direction perpendicular to the shock front over gas with a proton density, $n_{\mathrm{H}}⩽0.95n_{\mathrm{H}}^{\mathrm{f}}$. All other parameters are set to their standard values (see Table 2). Numbers in parenthesis are powers of ten.

In the text

	Fig. 8 f1 and f2 column density ratios computed across shocks at different velocities. The column densities are integrated in the direction perpendicular to the shock front over gas with a proton density, $n_{\mathrm{H}}⩽0.95n_{\mathrm{H}}^{\mathrm{f}}$. The colour of each point indicates the velocity of the shock. All other parameters are set to their standard values (see Table 2).
In the text

Fig. 9 Coloured table of the FWHM (in km s−1) of the absorption lines originating from the 3P2 level of neutral carbon computed across shocks for different velocities and different strengths of the pre-shock magnetic field. The FWHMs are obtained by fitting Gaussian profiles on the velocity profiles predicted by the shock model in the direction perpendicular to the shock front taking only into account gas with a proton density $n_{\mathrm{H}}⩽0.95n_{\mathrm{H}}^{\mathrm{f}}$. All other parameters are set to their standard values (see Table 2).

In the text

Fig. 10 Probability distribution function of the dissipation rate required to reproduce the column density of carbon at high pressure along the N lines of sight of JT11, assuming that all the shocks intercepted by these lines of sight are identical and defined by the parameters of the standard model: $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$, B0 = 1 μG, and VS = 80 km s−1. The probability distribution function is normalised so that its integral is equal to one.

In the text

Fig. 11 Same as Fig. 10 but for a pre-shock density varying between 0.1 and 2 cm−3 (from left to right), a strength of pre-shock magnetic field varying between 0.1 and 10 μG (from top to bottom), and six values of shock velocity: 20 (red), 40 (orange), 60 (green), 80 (turquoise), 100 (light blue), and 120 (dark blue) km s−1. For coherence, the colour code is the same than that used in Fig. 8. In many cases, only the last distribution (corresponding to VS = 120 km s−1) is visible, because those obtained at lower velocities are identical and hidden behind it. The bottom left panels are empty because the distributions fall outside the range of dissipation rates displayed.

In the text

Fig. 12 Coloured table of the median value of the dissipation rate (see Fig. 11) as a function of the shock velocity and the strength of the magnetic field for $n_{\mathrm{H}}^0=1{\mathrm{c}\mathrm{m}}^{-3}$. All the other parameters are set to their standard values (see Table 2). Models with a median dissipation rate within a factor of two of that obtained with the standard model are shown in turquoise. Models with a median dissipation rate between two and four times larger (resp. lower) than that obtained with the standard model are shown in green (resp. blue). All other models are shown in red.

In the text

Fig. 13 Basic analytical descriptions of the evolutions of the radius (blue curve) and the terminal velocity (red curve) of a 1D SNR (e.g. Draine 2011; Vink 2020), characterised by an ejecta mass of Mej = 1.4 M⊙ and an initial kinetic energy energy of E = 1051 erg, and expanding in a homogeneous medium with a density of 1 cm−3. The idealised successive phases known as the free expansion stage, the Sedov-Taylor stage, the pressure-driven stage, and the momentum conserving (or snowplow) stage are highlighted in red, orange, green, and blue, respectively.

In the text

Fig. 14 Schematic view of the calculation of the dissipation rates expected from SNRs along the lines of sight observed by JT11. A distribution of SNRs with terminal velocities between 30 and 200 km s−1 is drawn out of a simple description of supernovae explosions within the disc and the spiral arms of the Galaxy (panel a). The SNRs are colour-coded according to their terminal velocities. The lines of sight observed by JT11 (red lines) are analysed to count the number of SNR surfaces they intercept and compute the dissipation rates induced by the associated shocks (panels b and c). This calculation is done taking into account the 3D spatial distribution of SNRs, which is shown in projection in panels a and b. Panel c displays a schematic view of a single SNR located at a position Aj. An observed line of sight connects the position of the sun, S, to the position of the background star, B. The number of surfaces crossed by the line of sight for this single SNR (0, 1, or 2) and their inclination angle, ij, depend on the position of the point, Mj, along the SB line and its distance to Aj compared to the SNR radius.

In the text

Fig. 15 Top panel: probability distribution function of the dissipation rate along the lines of sight observed by JT11 averaged over 500 realisations of a simple and unclustered distribution of SNRs in the Galaxy (see Fig. 14) compared to that derived from the production and the excitation of neutral carbon in the standard model of interstellar shocks (see Fig. 10). Bottom panel: probability distribution function of the kinetic energy transfer rate derived from the observations of CO in molecular clouds located in the local ISM (Hennebelle & Falgarone 2012). All distributions shown in the top and bottom panels are normalised so that their integral is equal to one.

In the text