# C-type shock modelling – the effect of new H2–H collisional rate coefficients

A. Nesterenok, Duncan Bossion, Yohann Scribano, Francois Lique

► To cite this version:

A. Nesterenok, Duncan Bossion, Yohann Scribano, Francois Lique. C-type shock modelling – the effect of new H2–H collisional rate coefficients. Monthly Notices of the Royal Astronomical Society, 2019, 489 (4), pp.4520-4529. 10.1093/mnras/stz2441. hal-02172344

HAL Id: hal-02172344

https://hal.science/hal-02172344

Submitted on 17 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# C-type shock modelling – the effect of new H₂–H collisional rate coefficients

A.V. Nesterenok¹*, D. Bossion², Y. Scribano² and F. Lique³

¹ Ioffe Institute, 26 Polytechnicheskaya St., 194021 Saint Petersburg, Russia

² LUPM-UMR 5299, CNRS-Université de Montpellier, Place E. Bataillon, F-34095 Montpellier, France

³ LOMC-UMR 6294, CNRS-Université du Havre, 25 rue Philippe Lebon, BP 1123, F-76063 Le Havre, France

Accepted XXX. Received YYY; in original form ZZZ

## ABSTRACT

We consider collisional excitation of H₂ molecules in C-type shocks propagating in dense molecular clouds. New data on collisional rate coefficients for (de-)excitation of H₂ molecule in collisions with H atoms and new H₂ dissociation rates are used. The new H₂–H collisional data are state of the art and are based on the most accurate H₃ potential energy surface. We re-examine the excitation of rotational levels of H₂ molecule, the para-to-ortho-H₂ conversion, and H₂ dissociation by H₂–H collisions. At cosmic ray ionization rates ζ ≳ 10⁻¹⁶ s⁻¹ and at moderate shock speeds, the H/H₂ ratio at the shock front is mainly determined by the cosmic ray ionization rate. The H₂–H collisions play the main role in the para-to-ortho-H₂ conversion and, at ζ ≳ 10⁻¹⁵ s⁻¹, in the excitation of vibrationally excited states of H₂ molecule in the shock. The H₂ ortho-to-para ratio (OPR) of the shocked gas and column densities of rotational levels of vibrationally excited states of H₂ are found to depend strongly on the cosmic ray ionization rate. We discuss the applicability of the presented results to interpretation of observations of H₂ emission in supernova remnants.

**Key words:** shock waves – molecular data – molecular processes – ISM: supernova remnants

## 1 INTRODUCTION

The hydrogen molecule has two possible nuclear spin states due to the proton spin of 1/2 – ortho-H₂ and para-H₂. In the electronic ground state, the rotational levels of ortho-H₂ have odd values of the angular momentum *j* while the levels of para-H₂ have even *j* values. During radiative processes and non-reactive collisions of H₂ with other atoms and molecules, only transitions with even values of angular momentum change Δ*j* are permitted, thus preserving the ortho-to-para ratio (OPR). The ortho-/para-H₂ interconversion in the interstellar gas is possible via reactive collisions (which include exchange of protons) of H₂ with H, H⁺, H₃⁺ and other species, as well as via chemical reactions (H₂ formation on dust grains). The conversion between ortho-H₂ and para-H₂ also takes place on the surface of interstellar dust grains (Bron et al. 2016; Bovino et al. 2017; Furuya et al. 2019).

The H₂ OPR is an important parameter of the cold interstellar medium – the H₂ OPR is a controller of the cold cloud deuteration chemistry (e.g., Pineau des Forets et al. 1991; Flower et al. 2006; Pagani et al. 2009), it affects the molecular excitation (Troscompt et al. 2009), determines the heat capacity of the gas (Vaytet et al. 2014). In the cold molecular gas, the H₂ OPR slowly decays to a small value – 0.001 or even less (Pagani et al. 2011). The ortho-H₂ destruction is compensated by the H₂ forming on dust grains and by

destruction of ions (Flower et al. 2006; Le Bourlot 1991). The OPR of H₂ may not be at thermal equilibrium because the time-scale of the ortho-to-para conversion can be significantly longer than the time-scale of dynamic evolution of the molecular gas (Flower et al. 2006). Reactive collisions of H₂ molecule with hydrogen atoms have a substantial activation energy, ≈ 5000 K (Lique et al. 2014). This channel becomes dominant in a hot molecular gas, for example, behind interstellar shocks.

Shock waves are born due to large pressure disturbances in the interstellar medium caused by star-driven jets and winds, supernova explosions, and collisions between molecular clouds. Interstellar shock waves can be distinguished depending on the intensity of magnetic field, the ionization fraction of the gas, and the shock speed (Flower 2007; Draine 2011). Here we focus on magnetohydrodynamic non-dissociative shock waves – C-type shocks. Molecular hydrogen is an important tracer of non-dissociative shocks as it dominates cooling of the shocked gas (e.g., Kaufman & Neufeld 1996; Flower & Pineau des Forêts 2015; Tram et al. 2018). Excitation of molecular hydrogen and evolution of the H₂ OPR in interstellar shocks were studied theoretically by Timmermann (1998); Wilgenbus et al. (2000), and observed in many Galactic sources (e.g., Neufeld et al. 2006, 2007; Shinn et al. 2011, 2012; Neufeld et al. 2019). The main parameters which determine the speed of para-to-ortho-H₂ conversion in a C-type shock are the fraction of atomic hydrogen and the gas temperature. The rate coefficients for

<footnote>* E-mail:alex-n10@yandex.ru</footnote>
collisions involving H₂ molecules and H atoms are crucial for mod-
elling of H₂ OPR in such shocks.

The interstellar ultraviolet (UV) radiation field cannot penetrate into the interiors of dark clouds, and thus, cosmic ray particles and X-ray radiation (if present) appear the main drivers of the gas phase chemistry (Larsson et al. 2012). The abundance of atomic H in dark clouds is determined by the interplay between H₂ destruction by cosmic rays and its formation on dust grains (Goldsmith & Li 2005; Padovani et al. 2018). At high ionization rates ($\zeta \gtrsim 10^{-15}$ s⁻¹), the atomic hydrogen can have considerable abundance of the order of $10^{-3}-10^{-2}$ or even higher. Once a shock passes through molecular gas at such physical conditions, H₂–H collisions become important not only for the ortho-/para-state exchange of H₂ but also for the ro-vibrational excitation of H₂ molecules (e.g., Neufeld & Yuan 2008). In relatively fast C-type shocks, the abundance of H atoms can be higher due to partial dissociation of H₂ molecules.

Lique (2015) has performed nearly exact quantum time-
independent calculations of the rate coefficients for the collisional
(de-)excitation of H₂ by H atoms. These new data are based on a
highly accurate H₃ global potential energy surface, and the reactive
hydrogen exchange channels are taken into account rigorously. Bos-
sion et al. (2018) have performed quasi-classical trajectory (QCT)
calculations of rate coefficients for the collisional (de-)excitation of
H₂ by H atoms (including the three-body collisional dissociation)
for almost all rotational energy levels of the ground electronic state
of H₂ molecule. In the present paper, the new rate coefficients of
H₂–H collisions are incorporated into the model of C-type shock
propagating in a dense molecular cloud (Nesterenok 2018). We re-
examine (i) the excitation of rotational levels of H₂ molecule in
C-type shocks, (ii) the para-to-ortho-H₂ conversion, and (iii) the H₂
dissociation by H₂–H collisions. The effect of the elevated levels of
cosmic ray ionization rate on these processes is also studied.

## 2 DESCRIPTION OF THE CALCULATIONS

2.1 H₂ collisional rate coefficients

We take into account 298 rotational levels of the ground electronic state of H₂ molecule for which the Einstein coefficients are given by Wolniewicz et al. (1998). The level energies of H₂ are taken from Dabrowski (1984).

The data on collisional (de-)excitation of H₂ by H atoms by Wrathmall et al. (2007), which are usually used in astrophysical modelling, are restricted to non-reactive scattering. The treatment of ortho-/para-H₂ interconversion was usually adopted according to Le Bourlot et al. (1999), that in turn is based on the results of QCT calculations by Mandy & Martin (1993); Martin & Mandy (1995) and laboratory experiments (Schulz & Le Roy 1965; Schofield 1967). Until recently, the contribution of the reactive scattering channels to the H₂ excitation by H atoms remained a significant source of uncertainty, see discussion by Wrathmall et al. (2007); Lique et al. (2014).

Lique (2015) reported time-independent quantum mechanical calculations of rate coefficients for the collisional (de-)excitation of H₂ by H atoms. The calculations were based on a high accuracy global potential energy surface (PES) by Mielke et al. (2002). These new data are the first computed for the ortho-/para-H₂ interconversion at high temperature and using a purely quantum approach. The data were obtained for the lowest 54 rotational levels of H₂ with internal energies up to 22 000 K (highest energy level is $v = 3$, $j = 8$) and for kinetic temperatures ranging from 100 to 5000 K.

<table><caption>**Table 1.** Data on collisional rate coefficients for H₂ molecule.</caption>
  <thead>
    <tr>
      <th>Collisional partner</th>
      <th>Reference</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>H<sub>2</sub></td>
      <td>Wan et al. (2018),<br>Flower & Roueff (1998), Flower & Roueff (1999)</td>
    </tr>
    <tr>
      <td>He</td>
      <td>Flower et al. (1998)</td>
    </tr>
    <tr>
      <td>H</td>
      <td>Lique (2015), Bossion et al. (2018)</td>
    </tr>
    <tr>
      <td>thermal e<sup>-</sup></td>
      <td>Gerjuoy & Stein (1955), Ehrhardt et al. (1968),<br>England et al. (1988), Yoon et al. (2008)</td>
    </tr>
    <tr>
      <td>non-thermal e<sup>-</sup></td>
      <td>Tiné et al. (1997)</td>
    </tr>
    <tr>
      <td>H<sup>+</sup></td>
      <td>González-Lezana & Honvault (2017)</td>
    </tr>
  </tbody>
</table>

Bossion et al. (2018) presented QCT calculations of rate coef-
ficients for the ro-vibrational (de-)excitation of H₂ by H atoms for
temperatures up to 15 000 K, and for almost all rotational energy
levels of the ground electronic state of H₂ molecule. Their calcu-
lations are based on the same PES by Mielke et al. (2002). The
ortho-/para-H₂ interconversion and the three-body collisional dis-
sociation are included in their calculations. For this study, Bossion
et al. (2018) calculations have been extended to all the bound states
of H₂ molecule (internal energy up to dissociation limit ≈ 55 100 K)
and up to 20 000 K of collisional energy for the rate constant (the
maximal temperature attained by the neutral gas in the high speed
shocks in our simulations). For this high temperature regime we
extended the QCT cross sections up to 100 000 K of collisional en-
ergy in order to ensure convergence on the rate constants. We used
an energy-step of 2000 K up to that limit. Those QCT calculations
were performed considering only the ground electronic state of H₃.
Still the non-adiabatic effects are not expected to strongly impact
the results (the first excited electronic state lies for internal energies
over 58 000 K), and integrating the cross sections over the whole
energy domain induces an averaging effect. Moreover Bouakline
et al. (2010) studied the geometric-phase effect on the state-to-state
cross section of the H₃ system and found little contribution of the
non-adiabatic effects on the dynamics. It is to be noted that we
considered the quasi-bound states as pertaining to the three-body
dissociation channel. This assertion remains valid as long as the av-
erage lifetime of the quasi-bound molecules is lower than the typical
collision time, this is true for low to moderate density media.

For H₂–H₂ collisions, the data by Wan et al. (2018) were used
in our calculations. Wan et al. (2018) carried out quantum mechan-
ical close-coupling calculations of collisional (de-)excitation rates
based on the H₂–H₂ PES developed by Patkowski et al. (2008).
Initial rotational levels *j* = 2 – 31 of the ground vibrational state
and kinetic temperatures up to 10 000 K were considered in their
calculations. For the transitions involving other energy levels, the
data by Flower & Roueff (1998, 1999) are used in our simulations.
We take into account the ro-vibrational excitation of H₂ molecule
by fast electrons produced by cosmic ray ionization of the molecular
gas (Tiné et al. 1997). González-Lezana & Honvault (2017) studied
ortho-/para-H₂ conversion in collisions with H⁺ and provided colli-
sional rate coefficients for H₂ transitions between lowest rotational
levels belonging to *v* = 0 – 3 vibrational states. We use their data in
our simulations. Table 1 presents data on collisional rate coefficients
for H₂ molecule used in our calculations.

The dissociation of H₂ molecule by H₂ impact is taken into
account (Martin et al. 1998). Martin et al. (1998) considered H₂
dissociation from the ground vibration–rotation state. The rate of
dissociation of vibrationally excited H₂ molecule by H₂ impact is
estimated using dissociation coefficients given by Ceballos et al.
(2002).
## 2.2 C-type shock model

The simulations of C-type shock consist of two steps: (i) modelling of the chemical and thermal evolution of a static molecular cloud and (ii) the shock simulations. A complete gas-grain chemical network is taken into account – the gas-phase chemistry, the adsorption of gas species on dust grains and several desorption mechanisms, the chemical reactions on the surface of dust grains, the ion neutralization on dust grains. The population densities of energy levels of ions C I, C II and O I and molecules H₂, CO, H₂O are computed in parallel with the dynamical and chemical rate equations. The detailed information on the calculations can be found in the paper by Nesterenok (2018). In the appendix A we give brief overview of the model updates.

Table 2 summarizes the range of parameters of the shock wave explored in simulations. We assume that molecular cloud is shielded from the interstellar radiation field. The ionization rate $\zeta$ of the molecular gas is assumed to be due to cosmic rays, but we notice that X-rays affect molecular gas similarly (Mackey et al. 2019). The cosmic ray ionization rate<fnref n="1" /> is allowed to take values in the range $10^{-17} - 10^{-14}$ s⁻<fnref n="1" />. The highest value of $\zeta$ considered is comparable to the upper estimate of this parameter in molecular clouds in the vicinity of supernova remnants (Vaupré et al. 2014; Shingledecker et al. 2016). Increasing cosmic ray ionization rate even more leads to the decline of the molecular hydrogen abundance and the gas becomes almost entirely atomic at $10^{-12}$ s⁻<fnref n="1" /> (Bayet et al. 2011).

In molecular clouds not far away from sources of ionizing radiation, the ortho-to-para-H₂ conversion in collisions of H₂ with H⁺ and H⁺ is expected to be effective – the time-scale of this process is of the order of 0.1 Myr at the cosmic ray ionization rate $10^{-16}$ s⁻<fnref n="1" /> and the gas density $2 \times 10^4$ cm⁻³ (Pagani et al. 2013). As a result, the H₂ OPR in molecular clouds at such physical conditions is expected to be close to its steady-state value. At gas temperature $T_g = 30$ K, the steady-state H₂ OPR is equal to 0.03 (Flower et al. 2006). In our simulations, the H₂ OPR in the preshock gas has been taken equal to 0.1. The chemistries of ortho- and para-H₂ are not distinguished in the model. The ortho-/para-H₂ interconversion is considered in reactive collisions H₂–H and H₂–H⁺, and via the H₂ formation on dust grains.

Following Godard et al. (2019) we define the transition between the shock and the post-shock medium as the point at which most of the molecular radiation (95 per cent) induced by the shock has been emitted. For the sake of simplicity, we assume that the line of sight is perpendicular to the shock front, and the population densities of H₂ rotational levels are integrated along the gas flow. Once the column densities $N_j$ of rotational levels are obtained, the H₂ OPR for a given ortho-H₂ level j is calculated as described by Wilgenbus et al. (2000):

$$ \mathrm{OPR}_j = \mathrm{OPR}_{\mathrm{LTE},j} (T_{\mathrm{exc},j}) \frac{N_j}{N_{\mathrm{LTE},j} (T_{\mathrm{exc},j})}, \quad (1) $$

<table><caption>Table 2. Parameters of the shock model.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Preshock gas density, n<sub>H,tot</sub></td><td>2 × 10<sup>3</sup> – 2 × 10<sup>5</sup> cm<sup>-3</sup></td></tr><tr><td>Shock speed, u<sub>s</sub></td><td>5 – 100 km s<sup>-1</sup></td></tr><tr><td>Cosmic ray ionization rate, ζ</td><td>10<sup>-17</sup> – 10<sup>-14</sup> s<sup>-1</sup></td></tr><tr><td>Initial ortho-to-para-H<sub>2</sub> ratio</td><td>0.1</td></tr><tr><td>Magnetic field strength, β</td><td>1</td></tr><tr><td>Visual extinction, A<sub>V</sub></td><td>10</td></tr><tr><td>Micro-turbulence speed, v<sub>turb</sub></td><td>0.3 km s<sup>-1</sup></td></tr><tr><td>Grain surface area density</td><td>9.7 × 10<sup>-22</sup> cm<sup>2</sup> per H</td></tr></tbody></table>

$\beta$ is the dimensionless value of the transverse (to the shock velocity) magnetic field, $\beta = B[\mu G]/(n_{\mathrm{H,tot}}[\mathrm{cm}^{-3}])^{1/2}$.

where the excitation temperature between levels j + 1 and j − 1, $T_{\mathrm{exc},j}$, is given by

$$ T_{\mathrm{exc},j} = \frac{E_{j+1} - E_{j-1}}{k_B \ln (N_{j-1}g_{j+1}/N_{j+1}g_{j-1})}, \quad (2) $$

and $N_{\mathrm{LTE},j}$ is the column density obtained by the interpolation between the points j − 1 and j + 1 assuming local thermodynamic equilibrium at temperature $T_{\mathrm{exc},j}$. For a range of rotational levels, for example, for levels $2 \le j \le 6$, the H₂ OPR is calculated using the equation (1) for ortho-H₂ levels j = 3 and 5, and then averaged. The H₂ excitation diagram – the column densities in rotational levels (divided by the statistical weight $g_{v,j}$ of the level) as a function of the level energy – may exhibit curvature, indicating that the emitting gas has a range of kinetic temperatures. The error introduced by this curvature is discussed by Wilgenbus et al. (2000).

# 3 RESULTS

## 3.1 H₂ excitation diagram

Fig. 1 shows the H₂ excitation diagram for the C-type shock models using the different data on collisional rate coefficients for H₂ molecule, other things being equal. The C-type shock models were calculated with the shock speed $u_s = 25$ km s⁻<fnref n="1" />, the preshock gas density $n_{\mathrm{H,tot}} = 2 \times 10^4$ cm⁻³ and the cosmic ray ionization rate $\zeta = 10^{-15}$ s⁻<fnref n="1" />. The abundance of atomic hydrogen in the preshock gas is $1.7 \times 10^{-3}$, and it is approximately unchanged through the shock region – the H atom production in chemical reactions is not significant. At H/H₂ ratio in question, H₂–H collisions are dominant in the excitation of high-lying rotational levels of the ground vibrational state of H₂ molecule. Both H₂–H and H₂–H₂ collisions are important in the ro-vibrational excitation of H₂.

The populations of low-lying rotational levels of the ground vibrational state of H₂ are close to the local thermodynamic equilibrium, and column densities do not depend on the collisional data. But it is no longer valid for the high-lying levels – it is seen from Fig. 1 that populations of the high-lying rotational levels j ≥ 10 of the ground vibrational state strongly depend on the collisional data for both H₂–H and H₂–H₂ collisions. At high gas temperature $T_g \approx 1000$ K, the data by Lique (2015) predict higher rates, by a factor of few, of collisional (de-)excitation of energy levels of the ground vibrational state with high j compared with the data by Wrathmall et al. (2007) (with the inclusion of reactive channels as in Le Bourlot et al. (1999)). For energy levels belonging

<footnote n="1">¹ The rate of electron production in the unit volume by ionization of molecular gas is given by $\zeta n_{\mathrm{H}_2}$, where the He ionization is taken into account, He/H ratio is equal to 0.09. Our definition of the cosmic ray ionization rate is related to the definition of the parameter used in the UMIST Database for Astrochemistry (McElroy et al. 2013).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The H₂ excitation diagram for the C-type shock model using the different data on H₂–H and H₂–H₂ collisions. Parameters of the model are: $n_{H,\text{tot}} = 2 \times 10^4 \text{ cm}^{-3}$, $u_s = 25 \text{ km s}^{-1}$, $\zeta = 10^{-15} \text{ s}^{-1}$. *Circles*: the simulations employ the data for H₂–H collisions by Lique (2015); Bossion et al. (2018), data for H₂–H₂ collisions by Wan et al. (2018) for transitions involving H₂ energy levels within the ground vibrational state ($j \le 31$), and by Flower & Roueff (1998, 1999) for other transitions. *Crosses*: the H₂–H collisional rate coefficients by Martin & Mandy (1995); Le Bourlot et al. (1999); Wrathmall et al. (2007) are used. *Triangles*: the rate coefficients by Flower & Roueff (1998, 1999) are used for all H₂–H₂ collisional transitions. In the latter case, only column densities for energy levels of the ground vibrational state of H₂ are shown. The vibrational state of the energy level is marked by the color of the data point. For the ground vibrational state, the column densities for energy levels with angular momentum $j \ge 2$ are shown.</figcaption></figure>

to the vibrationally excited states the situation is reversed – collisional (de-)excitation rates by Wrathmall et al. (2007) are higher. These differences can be explained by the effect of the reactive hydrogen exchange channels that are not included in the calculations by Wrathmall et al. (2007), see discussion by Lique (2015). As a result, the rate coefficients by Lique (2015); Bossion et al. (2018) yield higher column densities for rotational levels belonging to the ground vibrational state and lower column densities for the ro-vibrationally excited levels than the data by Wrathmall et al. (2007). The higher the vibration quantum number of the energy level, the larger the difference between simulations. Column densities of ro-vibrationally excited levels have smaller scatter from a line on the excitation diagram for new rate coefficients.

The data by Wan et al. (2018) on H₂–H₂ collisions lead to lower (de-)excitation rates of H₂ molecule than the data by Flower & Roueff (1998, 1999) by a factor of about 1.5–2 for low rotational levels and by a factor of about 10 for levels with $j \ge 10$ at $T_g \approx 1000$ K. As a result, the rate of H₂ cooling is less effective in the model using the data by Wan et al. (2018), and the neutral gas temperature at the shock peak is higher by about 10 per cent in this case (the emission in H₂ rotational transitions of the ground vibrational state is the dominant cooling mechanism). The H₂–H₂ collisions do not contribute to the excitation of high-lying rotational

levels. The difference between the results of simulations that are based on the data by Wan et al. (2018) and by Flower & Roueff (1998, 1999) is explained by the difference in the temperature attained by the shocked gas, see Fig. 1.

Fig. 2 shows excitation diagrams of H₂ for shock models computed for a range of preshock gas densities, shock speeds and cosmic ray ionization rates. There is a strong sensitivity of H₂ excitation diagram (especially for ro-vibrationally excited levels) to the cosmic ray ionization rate. The width of C-type shock (and, hence, the temperature attained by the gas) depends on the momentum transfer rate between neutral and ionized fluids. At low gas densities ($n_{H,\text{tot}} \lesssim 10^4 \text{ cm}^{-3}$), there is a strong dependence of the shock width on the ionization fraction of the gas, which is determined by cosmic ray ionization. This effect was discussed by Wardle (1999); Gusdorf et al. (2012). Hence, the cosmic ray ionization affects the excitation of H₂ rotational levels in two ways – it increases the abundance of atomic hydrogen, and enhances the temperature of the shocked gas, see Fig. 2. At $n_{H,\text{tot}} \gtrsim 10^5 \text{ cm}^{-3}$, the charged dust grains dominate the momentum transfer and heating of the neutrals (Wardle 1998; Chapman & Wardle 2006). Ionization fraction of the gas has small effect on the structure of C-type shock at such densities. However, high cosmic ray fluxes affect hydrogen atom abundance that has significant effect on the H₂ excitation diagram at low shock speeds, see Fig. 2. At high shock speeds, the abundance of atomic hydrogen is increased by H₂ dissociation, and H₂ level populations approach local thermodynamic equilibrium. The effect of cosmic rays diminishes in this case. The total energy emitted by H₂ molecule at the given shock speed is almost independent on the cosmic ray ionization rate, see also Godard et al. (2019). Note, that Godard et al. (2019) conducted a similar study of the H₂ excitation in C-type shocks depending on different physical conditions, but in the presence of a strong external UV radiation field.

The non-thermal electrons produced by cosmic rays are an important excitation mechanism of H₂ vibrational states (along with the H₂ excitation in formation process) at low shock velocities, $u_s \le 15 \text{ km s}^{-1}$, and at moderate cosmic ray ionization rates, $\zeta \le 10^{-15} \text{ s}^{-1}$. However, column densities of vibrationally excited levels of H₂ are low in this case, see Fig. 2. At higher shock speeds, the gas temperature in the shock is high enough that collisional excitation by gas species dominates the level populations of low vibrational states. The same holds for higher cosmic ray ionization rates. This agrees with the conclusions by Tiné et al. (1997), who pointed out that H₂ excitation by non-thermal electrons is important at gas temperatures $T_g \le 1000$ K.

The effect of H₂–H⁺ collisions on the column densities of some vibrationally excited levels of H₂ is of the order of 10 per cent at low shock speeds and at moderate cosmic ray ionization rates. We note that effect of H₂ collisions with ions may be substantial if all abundant ions (which might undergo proton-exchange reactions with H₂) will be taken into account. For example, the abundance of H₃⁺ can be 1–2 orders of magnitude higher than that of H⁺.

## 3.2 Para-to-ortho-H₂ conversion

For the shock model presented in Fig. 1, the H₂ OPR – deduced from column densities integrated over the shock width and averaged over the energy levels $2 \le j \le 6$ of the ground vibrational state – is equal to 1.8. The same H₂ OPR is for energy levels of the first vibrational state ($0 \le j \le 6$). The data by Wrathmall et al. (2007) yield about 10 per cent higher OPR for energy levels of the ground vibrational state and about 5 per cent lower for the first vibrationally excited state. The abundance of atomic hydrogen in the preshock gas
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The H₂ excitation diagrams for C-type shock models with different speeds, preshock gas densities and cosmic ray ionization rates. The plots in one row correspond to the fixed density of the preshock gas, results for gas densities $n_{H,\text{tot}} = 2 \times 10^3, 2 \times 10^4, 2 \times 10^5 \text{ cm}^{-3}$ are shown. The preshock gas density in cm⁻³ and shock speed in km s⁻¹ are indicated in each plot. The highest shock speed for which the results are shown are 45, 50, 40 km s⁻¹ at $2 \times 10^3, 2 \times 10^4$ and $2 \times 10^5 \text{ cm}^{-3}$, respectively.</figcaption></figure>

is relatively high at the cosmic ray ionization rate in question, and the para-to-ortho-H₂ conversion is effective. The H₂ OPR reaches the high-temperature limit 3 in the hot shocked gas. However, there is a non-negligible contribution to H₂ emission of the initial part of the shock, where the gas temperature is high enough to excite vibrational states, but H₂ OPR has not enough time to reach the high-temperature limit, see also Wilgenbus et al. (2000). The spin

conversion process in H₂–H collisions is few times slower than the spin conserving process at high gas temperature (Lique 2015).

The correlation between the H₂ OPR and the excitation temperature of H₂ energy levels is often considered as a tracer of physical parameters of the shock wave (e.g., Wilgenbus et al. 2000; Neufeld et al. 2006). Fig. 3 shows the mean H₂ OPR for rotational levels belonging to the $v = 0$ and 1 vibrational states plotted against the mean excitation temperature of para-H₂. The data by Lique (2015)
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The H₂ OPR plotted against the mean excitation temperature of para-H₂ for a range of shock speeds, cosmic ray ionization rates and different preshock densities. The H₂ OPR and the mean excitation temperature of para-H₂ are calculated for rotational levels *j* = 2–6 belonging to the ground vibrational state (top panels), for levels *j* = 0–6 belonging to the vibrational state *v* = 1 (bottom). The preshock gas density *n*<sub>H,tot</sub> is indicated at the top of the panel column. Each point on a given curve represents a different shock speed, starting from 10 km s<sup>-1</sup> with an interval 5 km s<sup>-1</sup> (the shock speeds are indicated on the upper-left panel for ζ = 10<sup>-17</sup> s<sup>-1</sup> curve). For *n*<sub>H,tot</sub> = 2 × 10<sup>4</sup> cm<sup>-3</sup>, the simulations based on the data by Wrathmall et al. (2007); Le Bourlot et al. (1999) are shown (dashed line). The stars on the top panels show the H₂ OPR and the excitation temperature for supernova remnants W28, 3C 391, W44, IC 443C (in the order of increasing H₂ OPR) that are calculated based on the observational data by Neufeld et al. (2007). The plus markers show the data on supernova remnants G348.5-0.0, 3C 396, G349.7+0.2, Kes 69, Kes 17, G346.6-0.2 that are calculated based on the observational data by Hewitt et al. (2009). The filled cross markers on the bottom panels show the data on supernova remnant IC 443 (regions B, C, G). In this case, the H₂ OPR and the excitation temperature are estimated based on the observed column densities of energy levels *v* = 1, *j* = 1–3 (Shinn et al. 2011).</figcaption></figure>

predict lower para-/ortho-H₂ interconversion rates than the rates es-
timated based on Wrathmall et al. (2007); Le Bourlot et al. (1999).
In particular, for the transitions *j* = 2 → 1 and *j* = 4 → 3 within
the ground vibrational state, the data by Lique (2015) predict about
2–3 times lower de-excitation rates than the rate of para-to-ortho-
H₂ conversion given by Schulz & Le Roy (1965); Schofield (1967)
at 500–1000 K. As a result, data by Lique (2015) predict lower H₂
OPR for levels of the ground vibrational state, see Fig. 3. New rate
coefficients by Lique (2015); Bossion et al. (2018) predict higher
excitation temperature and higher integrated H₂ OPR (at moderate
and high shock speeds) for the vibrationally excited energy levels
than the analogous simulations with the data by Wrathmall et al.
(2007); Le Bourlot et al. (1999), see Fig. 3.

Fig. 3 shows that the H₂ OPR of the shocked gas strongly
depends on the cosmic ray ionization rate at moderate shock speeds
– there is a substantial increase of the H₂ OPR as the cosmic ray
ionization rate increases from 10⁻¹⁷ to 10⁻¹⁵ s⁻¹. On the other
hand, the excitation temperature of H₂ levels depends both on the

shock speed and the cosmic ray ionization rate. The H₂ OPR is
less than 3 for the ground and vibrationally excited states for most
of the shock models. At high cosmic ray ionization rates and at
high shock speeds, the integrated H₂ OPR may be slightly higher
than the local thermodynamic equilibrium value of 3, see Fig. 3.
At such physical conditions, the local H₂ OPR may attain values
up to 3.5 – 4 in the initial part of the shock. The reason of H₂
OPR values higher than 3 is that level populations are not in local
thermodynamic equilibrium. The critical densities (at which the
probabilities of collisional and radiative de-excitation are equal) are
n<sub>H₂</sub> ~ 10⁴ – 10⁵ cm⁻³ for high-lying rotational levels (v, j) > (0, 7)
(at T<sub>g</sub> ≈ 5000 K). At low gas density, most of the molecules populate
low-lying rotational levels even at high gas temperature. There is a
high ratio of statistical weights for energy levels with low angular
momentum in the limit of high temperature, for example, for levels
j = 5 and 4 this ratio is equal to 3.7. On the other hand, if the H atom
concentration is relatively high (at high cosmic ray ionization rate
or due to dissociation in the high speed shock), the characteristic
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The rate of H₂ dissociation in the hot shocked gas at the cosmic ray ionization rate of 10⁻¹⁶ s⁻¹, the preshock gas density of 2 × 10⁴ cm⁻³ and the shock speed of 65 km s⁻¹, the shock speed is close to critical. The contribution of each of the vibrational states *v* = 0, 2, 4, 6, 8 to the H₂–H dissociation channel is shown. The dissociation rate in H₂–H collisions calculated using the formula given by Le Bourlot et al. (2002) is also shown.</figcaption></figure>

time of para-to-ortho-H₂ conversion is lower than the time needed for the gas temperature reach its peak in the shock, and the H₂ OPR attains its local equilibrium value. As the gas cools and becomes dense downstream, the H₂ OPR approaches local thermodynamic equilibrium value of 3. We note that para-to-ortho-H₂ transitions are 3–5 times faster than ortho-to-para-H₂ transitions (Lique et al. 2012). The effect of high H₂ OPR (higher than 3) in low density gas is also discussed by Le Bourlot et al. (1999).

The H₂ OPR of ro-vibrationally excited levels may be higher or lower than that of the ground vibrational state depending on shock parameters. The effect of H₂–H⁺ collisions on the H₂ OPR and the excitation temperature is negligibly small.

### 3.3 Dissociation of H₂ molecules in H₂–H collisions

At high shock speeds, the main route of H₂ destruction in the hot gas is H₂ dissociation by collisions with H atoms. According to our calculations, the dissociation of H₂ in H₂–H₂ collisions contributes about 1–5 per cent to the total H₂ dissociation rate, dissociation from both the ground and vibrationally excited states being important. The H₂ destruction in collisions with ions and electrons is negligibly small at high shock speeds.

At some 'critical' speed $u_{crit}$ of the shock wave, there is almost all dissociation of molecular hydrogen. In this case the C-type solution ceases to exist, the shock wave becomes a J-type (Le Bourlot et al. 2002). Fig. 4 shows the dissociation rate of H₂ in H₂–H collisions at shock speed just below the critical speed. At high shock speeds, there is a fast rise of neutral gas temperature up to 15 000 – 20 000 K at the beginning of the shock as populations of H₂ and, hence, the H₂ cooling do not respond instantaneously to the changes of physical conditions. As vibrationally excited states

of H₂ become populated, the H₂ cooling is effective and the gas temperature falls. The first sharp peak in the H₂ dissociation curve in Fig. 4 is explained by the fast rise of the neutral gas temperature. While the second peak – by the increase of the number densities of atomic hydrogen and H₂ molecule on vibrationally excited states, that is lagged behind the temperature peak. The H₂ dissociation rate in H₂–H collisions depends on the number density of H₂ molecule on the vibrationally excited states (see Fig. 4). Thus, the inclusion of collisional excitation of high-lying rotational levels is important in calculations of H₂ dissociation rates. Le Bourlot et al. (2002) provided an approximate formula for the H₂ dissociation rate in collisions with H atoms which is based on the QCT calculations by Dove & Mandy (1986). The data by Bossion et al. (2018) predict 3–5 times higher dissociation rate than an estimate based on Le Bourlot et al. (2002), other things being equal.

At preshock gas density $n_{H,tot} = 2 \times 10^3$ cm⁻³ and low cosmic ray ionization rates, our calculations predict $u_{crit} \approx 100$ km s⁻¹ that is about 25 per cent higher than critical speed at this density published by Le Bourlot et al. (2002). At $n_{H,tot} = 2 \times 10^4$ cm⁻³ and $2 \times 10^5$ cm⁻³, our simulations provide $u_{crit} \approx 65$ km s⁻¹ and $u_{crit} \approx 40$ km s⁻¹, respectively, that are close to critical speeds found by Le Bourlot et al. (2002). Elevated levels of ionization rate lead to higher degree of ionization of the gas and higher H atom abundance in the preshock gas – critical velocities are lower in this case, see also discussion by Le Bourlot et al. (2002); Melnick & Kaufman (2015); Godard et al. (2019). At $\zeta = 10^{-14}$ s⁻¹, the maximal shock speed, consistent with the existence of a C-type shock solution, is about 30, 45, 40 km s⁻¹ at $n_{H,tot} = 2 \times 10^3, 2 \times 10^4, 2 \times 10^5$ cm⁻³, respectively. The dependence of critical speed on ionization rate is stronger at low gas density. Here we consider only solutions where the neutral flow stays supersonic over its entire trajectory. This condition may be broken at low gas densities and at high intensities of ionizing radiation (Godard et al. 2019).

## 4 DISCUSSION

### 4.1 Physical conditions in the preshock gas

The rate of para-to-ortho-H₂ conversion in the shock is determined by the abundance of atomic hydrogen. Usually, the gas-phase production of water from atomic oxygen is considered as the main source of H atoms in warm gas at moderate shock speeds (e.g., Neufeld et al. 2006). But most of the oxygen is locked up in H₂O on icy mantles of dust grains in our model, and this channel is not so effective in H atom formation as considered in previous works. The oxygen atoms presented in the gas-phase are adsorbed on dust grains and become hydrogenated on the grain surface during the chemical evolution of the dark cloud. The gas-phase chemical reactions – production of water, ammonia, methane – have noticeable effect on the abundance of atomic hydrogen in the shock only at low cosmic ray ionization rates, $\zeta \sim 10^{-17}$ s⁻¹. The H atom abundance in the shocked gas is low in this case, H/H₂ ~ 10⁻⁵ – 10⁻⁴. This explains low H₂ OPR at low cosmic ray ionization rates and at low and moderate shock speeds, see Fig. 3. At $\zeta \gtrsim 10^{-16}$ s⁻¹, ionization rate is the main parameter controlling the abundance of atomic hydrogen and, hence, it determines the rate of para-to-ortho-H₂ conversion. It explains strong dependence of the H₂ OPR on the cosmic ray ionization rate, see Fig. 3. In addition to the elevated fluxes of low-energy cosmic rays, molecular gas near supernova remnants also may be subject to the X-ray emission from the supernova interior
and to the UV radiation produced by nearby fast shocks, see also discussion by Yuan & Neufeld (2011).

The H/H₂ ratio in the preshock gas is close to the chemical equilibrium in our model. However, the molecular gas may not have a steady-state H/H₂ ratio at the time of passage of a shock wave. First, the molecular cloud may be younger than the characteristic time of H₂ formation from atomic gas that is of the order of 1 Myr for a cloud core with gas density 10⁴ cm⁻³ (Goldsmith et al. 2007). In this case cloud core will have H/H₂ ratio higher than equilibrium. Secondly, there is a finite response time of the dark cloud chemistry to the increase of the cosmic ray ionization rate. The H/H₂ ratio reaches equilibrium value at about (3 − 5) × 10⁴ yr after the change of cosmic ray ionization rate (at n<sub>H,tot</sub> = 2 × 10⁴ cm⁻³ and using the 'standard' H₂ production rate on dust grains). This time-scale does not depend on the ionization rate but is inversely proportional to the rate of H₂ formation. Regarding the second item, one can expect that H/H₂ ratio is close to its steady-state value for molecular clouds in the vicinity of middle-aged supernova remnants (with ages of about few 10⁴ yr).

Other important parameter is the H₂ OPR in the preshock gas. The ortho-H₂ decay time is about 0.5 – 1 Myr at ζ = 10⁻¹⁷ s⁻¹ and gas density of 10⁴ cm⁻³, and is in inverse ratio with the cosmic ray ionization rate and the gas density (Pagani et al. 2013; Bovino et al. 2017; Furuya et al. 2019). In molecular gas around young protostars, where the ionization rate is relatively low, the H₂ OPR in the preshock gas may not reach its equilibrium value at the time the protostellar outflow arrives the cloud. For example, Nisini et al. (2010); Giannini et al. (2011) analysed the *Spitzer* Infrared Spectrograph maps of H₂ pure rotational lines towards four outflows from Class 0 sources. They estimated the initial H₂ OPR in the preshock gas being close to 1. Note, accreting protostars can accelerate cosmic rays, and the ionization rate in the surrounding cloud may be comparable or even higher than due to Galactic cosmic rays (Podio et al. 2014; Padovani et al. 2016; Gaches & Offner 2018). At high cosmic ray ionization rate, ζ ≈ 10⁻¹⁴ s⁻¹, the steady-state H₂ OPR in the preshock gas is close to 1. In this case, the H₂ OPR quickly reaches the high temperature equilibrium value in the shock, and the results are almost independent of the initial H₂ OPR.

In our model, cosmic rays are the primary ionization agents and drivers of ion-molecule chemistry. However, at low visual extinctions, < 5 mag, interstellar UV radiation becomes the dominant source of ionization, and it significantly affects the chemistry of the gas. Melnick & Kaufman (2015), and more recently Godard et al. (2019), considered the propagation of low velocity molecular shocks in environments illuminated by a strong external UV radiation field. High levels of UV radiation, as cosmic rays, affect the structure of the C-type shock wave via the enhancement of the photoionization mechanisms and the increase of the dissociation of H₂ molecule (Godard et al. 2019). However, the chemical effects of high fluxes of cosmic rays are less pronounced. Even at the highest cosmic ray ionization rates detected in dense clouds in the vicinity of supernova remnants, the gas is molecular (H/H₂ ~ 0.01) and photodesorption is not effective enough to evaporate icy mantles of dust grains. However, we show that there is a strong dependence of the H₂ excitation on the cosmic ray ionization rate in C-type shocks.

## 4.2 Comparison with observations

Neufeld et al. (2007) carried out spectroscopic mapping observations using the Infrared Spectrograph (IRS) of the *Spitzer Space Telescope* towards the regions where a supernova remnant interacts with a molecular cloud – W44 region E, W28 F, 3C 391, and

IC 443 C. These observations led to the detection of the S(0)–S(7) pure rotational lines of molecular hydrogen. Emission in these lines likely originates in molecular gas subject to a slow, nondissociative shock (Neufeld et al. 2007). Hewitt et al. (2009) carried out analogous observations towards another six supernova remnants that show evidence of shocked molecular gas. The calculated H₂ OPR and the mean excitation temperature are plotted in Fig. 3 for these data. Mean excitation temperatures of para-H₂ rotational levels 2 ≤ *j* ≤ 6 lie in the range 400–700 K that are consistent with the presence of non-dissociative shocks of speed about 10–20 km s⁻¹. Tram et al. (2018) modelled the bow shock by a statistical distribution of planar shocks, and found that emission by bow shock is generally dominated by low-velocity shocks.

Vaupré et al. (2014) used molecular line observations to constrain cosmic ray ionization rate in molecular clouds in the vicinity of the W28 supernova remnant. Towards positions located close to the supernova remnant, they found cosmic ray ionisation rates of the order of 10⁻¹⁵ s⁻¹. Indriolo et al. (2010) observed interstellar absorption of H₃⁺ along sightlines that pass through diffuse gas near the supernova remnant IC 443 and deduced ionization rates of the same order. The H₂ OPR of 0.5 – 1.5 and excitation temperatures 400 – 700 K for low-lying pure rotational levels are consistent with the cosmic ray ionization rate ζ ≈ 10⁻¹⁵ s⁻¹ and preshock gas density n<sub>H,tot</sub> ≈ 10⁴ cm⁻³, see Fig. 3. The observations of H₂ excitation must be supplemented with the independent determination of ionization rate or gas density in order to reveal a unique set of physical conditions, see also discussion by Godard et al. (2019). We note, that Gusdorf et al. (2012) discussed the effect of cosmic ray ionization rate on the populations of rotational levels of H₂ in their analysis of the W28 F, but the dependence of H₂ OPR on the cosmic ray ionization rate was not considered in their work.

Shinn et al. (2011, 2012) presented near-infrared spectral studies of the shocked H₂ gas in the supernova remnant IC 443 using the Infrared Camera aboard the *AKARI* satellite. Shinn et al. (2012) estimated the H₂ OPR in IC 443 B equal to 2.4 ± 0.4 and 2.1 ± 0.3 for levels (v, j) = (0, 11–13) and (1, 1–3), respectively. For regions C and G in IC 443, the H₂ OPR for ro-vibrationally excited levels is close to 3. Shinn et al. (2012) suggested that obtained non-equilibrium OPR for region B probably originates from the refraction of H₂ on dust grains behind the dissociative J-type shock. According to our calculations, the H₂ OPR for vibrationally excited rotational levels is lower than 3 for most of the models of C-type shocks at ζ ≈ 10⁻¹⁵ s⁻¹, see Fig. 3. The value obtained for IC 443 B is consistent with non-dissociative shock wave in the medium with high cosmic ray ionization rate. The H₂ OPR close to 3 may be result of the emission from dissociative shock wave or the high initial value of H₂ OPR. The difficulty in interpretation of observations is that multiple shocks, both non-dissociative and dissociative, might be present along the line of sight (e.g., Neufeld et al. 2007).

The important question that is not discussed in the current study is the shapes of spectral lines of H₂ molecule and velocity shifts between lines of para-H₂ and ortho-H₂. Recently, Neufeld et al. (2019) carried out spectrally resolved observations of pure rotational lines of H₂ molecule using the EXES instrument on *SOFIA* towards the shock region HH7. They reported the detection of systematic velocity shifts between the low-lying rotational lines of ortho-H₂ and para-H₂ that provides the first definitive evidence for the conversion of para-H₂ to ortho-H₂ behind shock waves. The shock models that use different data on collisional rate coefficients for H₂ molecule predict similar behaviour of OPR for low-lying rotational levels of H₂, see Fig. 3. In this case, the effect of new collisional data on
the magnitude of the velocity shifts between ortho- and para-H₂ rotational lines is expected to be small.

## 5 CONCLUSIONS

We study the effect of new H₂–H collisional rate coefficients published by Lique (2015); Bossion et al. (2018) on the excitation of molecular hydrogen in C-type shocks. The new H₂–H collisional data are state of the art and are based on the most accurate H₃ potential energy surface. The main results of the paper are summarized below:

(i) At high gas temperatures, the data by Lique (2015) predict lower rates, by a factor of 1 – 3, of collisional (de-)excitation of H₂ rotational levels belonging to vibrationally excited states compared with the data by Wrathmall et al. (2007) with the inclusion of reactive channels. As a result, the rate coefficients for H₂–H collisions by Lique (2015); Bossion et al. (2018) yield lower column densities of H₂ rotational levels of vibrationally excited states than the data by Wrathmall et al. (2007).

(ii) The effect of new collisional data is substantial in the determination of H₂ OPR for rotational levels of the first vibrationally excited state. The H₂ OPR integrated over the shock length is lower than high-temperature equilibrium value for most of the shock models at $\zeta \lesssim 10^{-15} \text{ s}^{-1}$. It is valid for rotational levels of both the ground and vibrationally excited states.

(iii) The para-to-ortho-H₂ conversion rate in the shock strongly depends on the cosmic ray ionization rate at moderate shock speeds – there is a substantial increase of the H₂ OPR of the shocked gas as the cosmic ray ionization rate increases from $10^{-17}$ to $10^{-15} \text{ s}^{-1}$. The H₂ OPR and the excitation temperature can be used to constrain physical parameters from observations of molecular hydrogen in interstellar shocks.

(iv) The H₂–H collisions are the main H₂ dissociation channel at high shock speeds. The data by Bossion et al. (2018) predict several times higher dissociation rate in H₂–H collisions than an estimate based on calculations by Dove & Mandy (1986). The H₂ dissociation takes place from high-lying ro-vibrationally excited H₂ levels. The data on collisional excitation of high-lying H₂ levels are as important as dissociation rates.

(v) The H₂–H₂ collisions contribute about 1 – 5 per cent to the total dissociation rate. The dissociation of H₂ in H₂–H₂ collisions from the ground vibrational state and from vibrationally excited states is taken into account.

(vi) The data by Wan et al. (2018) on H₂–H₂ collisions lead to lower (de-)excitation rates of H₂ molecule than the data by Flower & Roueff (1998, 1999) by a factor of about 1.5 – 2 for low rotational levels. As a result, the rate of H₂ cooling is less effective in the model using new data on H₂–H₂ collisions.

(vii) The effect of the excitation of H₂ by H⁺ ions and by non-thermal electrons produced by cosmic rays is non-negligible only at low shock speeds and at moderate cosmic ray ionization rates.

## ACKNOWLEDGEMENTS

AVN acknowledges the Russian Science Foundation (RSF) grant 16-12-10225. DB, YS and FL are supported by the 'Programme National Physique et Chimie du Milieu Interstellaire' (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. YS would like to thank support from the Computing Center of the 'Institut National de Physique Nucléaire' (IN2P3). We are grateful

to Dr. Wan for providing the collisional rate coefficients for H₂–H₂ collisions.

## REFERENCES

Bayet E., Williams D. A., Hartquist T. W., Viti S., 2011, MNRAS, 414, 1583

Black J. H., van Dishoeck E. F., 1987, ApJ, 322, 412

Bossion D., Scribano Y., Lique F., Parlant G., 2018, MNRAS, 480, 3718

Bouakline F., Althorpe S., Larregaray P., Bonnet L., 2010, Mol. Phys., 108, 969

Bovino S., Grassi T., Schleicher D. R. G., Caselli P., 2017, ApJ, 849, L25

Bron E., Le Petit F., Le Bourlot J., 2016, A&A, 588, A27

Ceballos A., Garcia E., Laganà A., 2002, J. Phys. Chem. Ref. Data, 31, 371

Chabot M., Béroff K., Gratier P., Jallat A., Wakelam V., 2013, ApJ, 771, 90

Chapman J. F., Wardle M., 2006, MNRAS, 371, 513

Dabrowski I., 1984, Can. J. Phys., 62, 1639

Dove J. E., Mandy M. E., 1986, ApJ, 311, L93

Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium. Princeton University Press, Princeton

Ehrhardt H., Langhans L., Linder F., Taylor H. S., 1968, Phys. Rev., 173, 222

England J. P., Elford M. T., Crompton R. W., 1988, Aust. J. Phys., 41, 573

Flower D., 2007, Molecular Collisions in the Interstellar Medium. Cambridge University Press, New York

Flower D. R., Pineau des Forêts G., 2015, A&A, 578, A63

Flower D. R., Roueff E., 1998, J. Phys. B: At. Mol. Opt. Phys., 31, 2935

Flower D. R., Roueff E., 1999, J. Phys. B: At. Mol. Opt. Phys., 32, 3399

Flower D. R., Roueff E., Zeippen C. J., 1998, J. Phys. B: At. Mol. Opt. Phys., 31, 1105

Flower D. R., Pineau Des Forêts G., Walmsley C. M., 2006, A&A, 449, 621

Furuya K., Aikawa Y., Hama T., Watanabe N., 2019, arXiv e-prints, p. arXiv:1908.01966

Gaches B. A. L., Offner S. S. R., 2018, ApJ, 861, 87

Gerjuoy E., Stein S., 1955, Phys. Rev., 97, 1671

Giannini T., Nisini B., Neufeld D., Yuan Y., Antoniucci S., Gusdorf A., 2011, ApJ, 738, 80

Godard B., Pineau des Forêts G., Lesaffre P., Lehmann A., Gusdorf A., Falgarone E., 2019, A&A, 622, A100

Goldsmith P. F., Li D., 2005, ApJ, 622, 938

Goldsmith P. F., Li D., Krčo M., 2007, ApJ, 654, 273

González-Lezana T., Honvault P., 2017, MNRAS, 467, 1294

Gusdorf A., Anderl S., Güsten R., Stutzki J., Hübers H.-W., Hartogh P., Heyminck S., Okada Y., 2012, A&A, 542, L19

Heays A. N., Bosman A. D., van Dishoeck E. F., 2017, A&A, 602, A105

Hewitt J. W., Rho J., Andersen M., Reach W. T., 2009, ApJ, 694, 1266

Hocuk S., Szűcs L., Caselli P., Cazaux S., Spaans M., Esplugues G. B., 2017, A&A, 604, A58

Indriolo N., Blake G. A., Goto M., Usuda T., Oka T., Geballe T. R., Fields B. D., McCall B. J., 2010, ApJ, 724, 1357

Kalvāns J., 2018, MNRAS, 478, 2753

Kaufman M. J., Neufeld D. A., 1996, ApJ, 456, 611

Larsson M., Geppert W. D., Nyman G., 2012, Rep. Prog. Phys., 75, 066901

Le Bourlot J., 1991, A&A, 242, 235

Le Bourlot J., Pineau des Forêts G., Flower D. R., 1999, MNRAS, 305, 802

Le Bourlot J., Pineau des Forêts G., Flower D. R., Cabrit S., 2002, MNRAS, 332, 985

Lique F., 2015, MNRAS, 453, 810

Lique F., Honvault P., Faure A., 2012, J. Chem. Phys., 137, 154303

Lique F., Honvault P., Faure A., 2014, International Reviews in Physical Chemistry, 33, 125

Mackey J., Walch S., Seifried D., Glover S. C. O., Wünsch R., Aharonian F., 2019, MNRAS, 486, 1094

Mandy M. E., Martin P. G., 1993, ApJS, 86, 199

Martin P. G., Mandy M. E., 1995, ApJ, 455, L89

Martin P. G., Keogh W. J., Mandy M. E., 1998, ApJ, 499, 793

McElroy D., Walsh C., Markwick A. J., Cordiner M. A., Smith K., Millar T. J., 2013, A&A, 550, A36
Melnick G. J., Kaufman M. J., 2015, ApJ, 806, 227

Mielke S. L., Garrett B. C., Peterson K. A., 2002, J. Chem. Phys., 116, 4142

Murga M. S., Varakin V. N., Stolyarov A. V., Wiebe D. S., 2019, Astronomy Reports, 63, 633

Nesterenok A. V., 2018, Ap&SS, 363, 151

Neufeld D. A., Yuan Y., 2008, ApJ, 678, 974

Neufeld D. A., et al., 2006, ApJ, 649, 816

Neufeld D. A., Hollenbach D. J., Kaufman M. J., Snell R. L., Melnick G. J., Bergin E. A., Sonnentrucker P., 2007, ApJ, 664, 890

Neufeld D. A., DeWitt C., Lesaffre P., Cabrit S., Gusdorf A., Tram L. N., Richter M., 2019, ApJ, 878, L18

Nisini B., Giannini T., Neufeld D. A., Yuan Y., Antoniucci S., Bergin E. A., Melnick G. J., 2010, ApJ, 724, 69

Öberg K. I., 2016, Chem. Rev., 116, 9631

Padovani M., Marcowith A., Hennebelle P., Ferrière K., 2016, A&A, 590, A8

Padovani M., Galli D., Ivlev A. V., Caselli P., Ferrara A., 2018, A&A, 619, A144

Pagani L., et al., 2009, A&A, 494, 623

Pagani L., Roueff E., Lesaffre P., 2011, ApJ, 739, L35

Pagani L., Lesaffre P., Jorfi M., Honvault P., González-Lezana T., Faure A., 2013, A&A, 551, A38

Patkowski K., Cencek W., Jankowski P., Szalewicz K., Mehl J. B., Garberoglio G., Harvey A. H., 2008, J. Chem. Phys., 129, 094304

Pineau des Forets G., Flower D. R., McCarroll R., 1991, MNRAS, 248, 173

Podio L., Lefloch B., Ceccarelli C., Codella C., Bachiller R., 2014, A&A, 565, A64

Schofield K., 1967, Planet. Space Sci., 15, 643

Schulz W. R., Le Roy D. J., 1965, J. Chem. Phys., 42, 3869

Shingledecker C. N., Bergner J. B., Le Gal R., Öberg K. I., Hincelin U., Herbst E., 2016, ApJ, 830, 151

Shinn J.-H., Koo B.-C., Seon K.-I., Lee H.-G., 2011, ApJ, 732, 124

Shinn J.-H., Lee H.-G., Moon D.-S., 2012, ApJ, 759, 34

Timmermann R., 1998, ApJ, 498, 246

Tiné S., Lepp S., Gredel R., Dalgarno A., 1997, ApJ, 481, 282

Tram L. N., Lesaffre P., Cabrit S., Gusdorf A., Nhung P. T., 2018, MNRAS, 473, 1472

Troscompt N., Faure A., Maret S., Ceccarelli C., Hily-Blant P., Wiesenfeld L., 2009, A&A, 506, 1243

Vaupré S., Hily-Blant P., Ceccarelli C., Dubus G., Gabici S., Montmerle T., 2014, A&A, 568, A50

Vaytet N., Tomida K., Chabrier G., 2014, A&A, 563, A85

Wakelam V., et al., 2017, Molecular Astrophysics, 9, 1

Wan Y., Yang B. H., Stancil P. C., Balakrishnan N., Parekh N. J., Forrey R. C., 2018, ApJ, 862, 132

Wardle M., 1998, MNRAS, 298, 507

Wardle M., 1999, ApJ, 525, L101

Wilgenbus D., Cabrit S., Pineau des Forêts G., Flower D. R., 2000, A&A, 356, 1010

Wolniewicz L., Simbotin I., Dalgarno A., 1998, ApJS, 115, 293

Wrathmall S. A., Gusdorf A., Flower D. R., 2007, MNRAS, 382, 133

Yoon J.-S., Song M.-Y., Han J.-M., Hwang S. H., Chang W.-S., Lee B., Itikawa Y., 2008, J. Phys. Chem. Ref. Data, 37, 913

Yuan Y., Neufeld D. A., 2011, ApJ, 726, 76

see Nesterenok (2018). A set of branching ratios for the reactions involving carbon-chain species provided by Chabot et al. (2013) is used. Photodissociation and photoionisation rates of chemical species in the interstellar radiation field and in the cosmic-ray induced UV flux are updated according to Heays et al. (2017). The dissociation rate of molecular hydrogen by cosmic-rays was taken according to Padovani et al. (2018). The photochemistry experiments suggest that the photodissociation of adsorbed molecules is less efficient than the corresponding process in the gas phase (Öberg 2016; Murga et al. 2019). Here we use solid/gas photodissociation coefficient ratio equal to 0.1 (Kalvāns 2018). Note that the rate of dissociative gas-phase reaction H₂–H calculated using the rate coefficients provided by UDfA is more than an order of magnitude higher than the accurate one (Bession et al. 2018).

The properties and size distribution of dust grains determine the rate of H₂ formation and, as a result, the abundance of atomic hydrogen in the cold molecular gas (along with the cosmic ray ionization). A single-size grain model is considered, the grain radius is taken equal to 0.05 µm and the dust–gas mass ratio is equal to 0.01. The corresponding grain surface area is about 10⁻²¹ cm² per H. The dust temperature in simulations of the chemical evolution of the dark cloud is equal to 10 K, that approximately corresponds to starless cold dense regions at high visual extinctions A_V ~ 10 (Hocuk et al. 2017).

Once the H atom adsorbs on the surface of a dust grain, it quickly reacts with one of the icy species (CO, HCO, NO, HNO and others), and, as a result, H₂ molecules form through hydrogen abstraction reactions. The rate of H₂ formation mostly depends on the adsorption rate of H atoms, which in turn depends on the specific surface area of dust grains and the gas temperature. The effective formation rate of H₂ can be expressed as R = γnH,tot nH, where the rate coefficient γ ≈ (2 − 3) × 10⁻¹⁷ cm³ s⁻¹ in our simulations. The rate of H₂ formation in a dense interstellar gas is poorly known. In diffuse clouds, H₂ forms with a rate of R_H₂ = (3 − 5) × 10⁻¹⁷ cm³ s⁻¹ (Wakelam et al. 2017). Here, we assume that H₂ molecule quickly leaves dust grain after the formation, and that one-third of the binding energy of H₂ is deposited statistically as internal excitation of the newly formed H₂ molecule (Black & van Dishoeck 1987). However, the effect of H₂ formation on the H₂ excitation in C-type shocks is negligibly small.

In our model, the hydrogen atom concentration in the gas is approximately equal to 0.5 cm⁻³ at ζ = 10⁻¹⁷ s⁻¹ independent of the gas density. Goldsmith & Li (2005) studied H I self-absorption in five dark clouds and found central number densities of H atoms in the range 2 − 6 cm⁻³ that corresponds to the cosmic ray ionization rate of ζ ~ 10⁻¹⁶ s⁻¹.

## APPENDIX A: CHEMISTRY

In simulations of chemical evolution of the molecular cloud, species are assumed to be initially in atomic form except for hydrogen, which is assumed to be molecular. A simple zero-dimensional model is considered in these simulations, see Nesterenok (2018). The calculated chemical composition at *t* = 0.5 Myr is chosen for the shock wave modelling.

The gas-phase chemical network used here is based on the UMIST Database for Astrochemistry (UDfA), 2012 edition (McElroy et al. 2013). Several updates were done to the chemical network,