© The Authors 2023

1 Introduction

Shocks are inherently out-of-equilibrium time-dependent phenomena that permeate space. They appear over a wide range of scales, ranging from, for example, accretion onto stars or proto-planetary disks, winds and jets driven by accreting (proto)stars, planetary nebulae, supernova remnants, starburst galaxies, jets from active galactic nuclei (AGN), and to galaxy-galaxy collisions (physical sizes ranging from subastronomical unit to kiloparsec scales; e.g., Bally 2016; Wright et al. 1993; Mouri 1994; Goldader et al. 1997; Appleton et al. 2006). Common to all these phenomena is that the input kinetic energy flux dissipated by the shock accelerates, heats, and compresses the medium.

When the medium cools down, radiation is emitted, which we observe. To understand the physical origin of emission (e.g., preshock density, shock velocity) and the energetic processing taking place in shocks, it is thus necessary to reverse engineer the observed light. Doing so requires models.

One of the often-used tracers of shocks is molecular hydrogen, H₂ (e.g., Hollenbach & McKee 1989; Kaufman & Neufeld 1996; Rosenthal et al. 2000). This is the most abundant molecule in the interstellar medium by some four orders of magnitude over CO and H₂O. The molecule is the lightest, and so it has the most widely spaced rotational levels (J = 1 has E_up/k_B = 170 K and J = 2 has E_up/k_B = 510 K). As such, it is predominantly excited in warm (T ≳ 10² K) and hot (T ≳ 10³ K) molecular gas. This molecule has no permanent dipole moment, and only forbidden electric quadrupole transitions occur, although at low probability. The main reason H₂ emission is still bright is because of its high abundance.

H₂ emission is readily observed from the ground, particularly in higher-excited rovibrational transitions at near-infrared wavelengths (e.g., Froebrich et al. 2015). The brightest of these is typically the υ = 1–0 S(1) line at 2.12 µm. A few pure rotational lines are also accessible from the ground, and the line profiles may even be velocity resolved on telescopes such as the Very Large Telescope (VLT, Santangelo et al. 2014). However, it is necessary to go above the atmosphere to observe the lower-excited pure rotational transitions of H₂. Space-based telescopes such as the Infrared Space Observatory (ISO) and the Spitzer Space Telescope (Spitzer) both observed these transitions toward numerous shocked regions (e.g., Rosenthal et al. 2000; Neufeld et al. 2006; Valentijn & van der Werf 1999; Lutz et al. 2003; Verma et al. 2005), as did the Stratospheric Observatory For Infrared Astronomy (SOFIA, Reach et al. 2019; Neufeld et al. 2019). Now the James Webb Space Telescope (JWST) is doing the same (e.g., García-Bernete et al. 2022; Berné et al. 2022; Yang et al. 2022; Appleton et al. 2023; Álvarez-Márquez et al. 2023). Particularly, the MIRI instrument is observing the rotational H₂ transitions with a gain in sensitivity and spatial resolution of two orders of magnitude compared with Spitzer, and an increase in spectral resolution of a factor five (e.g., Figs. 7 and 8 of Rigby et al. 2023). Similar improvements are reached with the NIRSpec instrument compared with the VLT-SINFONI integral-field unit, allowing deep observations of the rovibrational lines of H₂. The wavelength coverage of NIRSpec, NIRCam, and MIRI are illustrated in Fig. 1, which shows a simulated H₂ spectrum with the instrument wavelength coverages displayed.

Planning and interpreting the abovementioned observations is often done by use of models. With models, it is possible to constrain, for example, the shock velocity and preshock density, which together give the input kinetic energy flux, 1/2 , where ρ is the mass density and υ_s is the shock velocity. In molecular shocks, a comparison reveals that up to 50% of the input energy is radiated away in H₂ emission (Kaufman & Neufeld 1996), depending on shock conditions, making H₂ the dominant coolant in these shocks. Spitzer particularly opened up for characterization of the pure rotational H₂ lines. Observations and subsequent modeling revealed that most H₂ emission could be reproduced by shock models (e.g., in protostellar outflows; Maret et al. 2009; Dionatos et al. 2010). However, when additional constraints, such as the H/H₂ ratio and the cooling length are included for pro-tostellar outflows, a single shock model no longer reproduces observations (Nisini et al. 2010). Instead, as argued, the observational beam likely catches different shocks, or more complex shock geometries than 1D, which is to be expected; this is not just the case for protostellar outflows, but also observations of shocks in the diffuse gas of starburst and colliding galaxies (Kristensen et al. 2008; Gustafsson et al. 2010; Lesaffre et al. 2013; Tram et al. 2018; Lehmann et al. 2022). Irrespective of the specific science case, the first step in comparing observations to models is to have the models available.

The Paris–Durham shock code (e.g., Godard et al. 2019, and references therein) has been developed and maintained for more than 35 yr (Flower et al. 1985). The code can either find jump (J-type shocks) or continuous (C-type shocks) solutions depending on the input physical parameters. Recent developments include the treatment of an external UV radiation field (Godard et al. 2019), and self-irradiation in high-velocity shocks (υ_s ≳ 30 km s⁻¹; Lehmann et al. 2022). Here we present the results of running a large grid of simulations of (externally irradiated) shocks with the goal of exploring how the input energy flux (kinetic and radiative) is reprocessed and ultimately results in H₂ emission. These model predictions can be used directly to interpret, for example, JWST observations of shock emission.

The paper is organized as follows. Section 2 describes the shock model and the model grid, with a particular emphasis on H₂ excitation and emission. The section also describes which physical quantities were extracted from the models, and the methodology applied. Section 3 describes the results and provides a discussion of these results. Finally, the main points are summarized in Sect. 4.

Fig. 1 Synthetic H2 spectrum produced with a shock model with velocity 30 km s−1, preshock density 104 cm−3, a transverse magnetic field strength of 10 µG, and no external UV radiation. Wavelength ranges of the NIRSpec and MIRI spectrographs, as well as the wide-, medium-, and narrow-band filters for NIRCam and the MIRI filters on JWST are indicated as black and gray horizontal bars. The colors are for lines with different vibrational upper levels. The resolving power is assumed to be uniform across the wavelength range at λ/Δλ = 2500.

2 Model and grid description

The current version of the multifluid shock code is extensively described in Godard et al. (2019) and references therein, and only the main relevant points will be described here. These points particularly relate to H₂ emission and other observable diagnostics, but also how the initial shock conditions are calculated. The code is publicly available¹, and the entire grid presented in this paper is also available on the ISM platform². In Appendix A, we provide an introduction to this platform and demonstrate how it can be used.

2.1 Initial conditions

The main focus of this paper is on H₂, and so the chemistry considered in this paper and, more importantly, in the models run, is a gas-phase-only chemistry. That is, grain adsorption and desorption processes are not included. The only exceptions are the formation of H₂ on grains, and grain erosion for the release of elemental Si, Fe, etc. into the gas phase. Photochemistry is included in all steps of the calculation; readers can refer to the text below for more details.

Our assumption is that the initial conditions are in equilibrium, that is, thermal and chemical equilibrium with or without an incident radiation field. Running a shock model therefore requires multiple steps, all done using the Paris-Durham code (see Godard et al. 2019, for details). This code simulates steady-state gas equilibrium, photon-dominated regions (PDRs), or shocks. These steps are illustrated in Fig. 2. First, a chemical steady-state calculation is run with the given density and radiation field. For irradiated shocks, the next step is to take the final equilibrium conditions from the chemical steady-state calculation and use these as input for a PDR calculation, where a tracer particle is advected at a small velocity (≤0.01 km s⁻¹) from an A_V of 10⁻⁹ to 10⁻¹. The advection speed is chosen such that the time it takes to cross the PDR front is long enough that equilibrium is reached; this timescale is 10⁵–10⁹ yr for high to low densities. The choice of a final A_V of 0.1 is motivated by two considerations. First, the primary focus of this paper is H₂ and the A_V thus needs to be high enough that the preshock gas is substantially molecular (molecular fraction ≥0.1) for the majority of the G₀ values here, specifically the part of the grid where G₀/n_H < 1. Second, the A_V should be low enough that H₂ is not fully self-shielded. These two conditions are met at an A_V of 0.1. The final conditions, in terms of steady-state abundances, temperature, and H₂ level populations, are then used as the input physical conditions of the shock calculation. The shock is run in the final step.

The initial elemental abundances are provided in Table 1. Of particular importance is the abundance of polycyclic aromatic hydrocarbons (PAHs). In the model, a representative PAH molecule is included, C₅₄H₁₈ and its singly charged ions. Table 1 reports the amount of H and C locked up in this PAH for a PAH abundance of X(PAH) = 10⁻⁶. The grain temperature is kept fixed at 15 K.

We cover a 6D parameter space with preshock density (n_H = 2 n(H₂) + n(H)), shock velocity (υ_s), strength of the transverse magnetic field³ (b), external Uv radiation (G₀ in units of the field from Mathis et al. 1983), H₂ cosmic-ray ionization rate (ζ_H2), and the fractional abundance of the PAHs (X(PAH)). The parameter space is presented in Table 2. Depending on the initial conditions, the code either finds a Jump (J-type) solution or a Continuous (C-type) solution (see below, Sect. 3.1 for more details). Throughout this paper, we use two shock models to illustrate differences when changing b from 0.1 to 1.0; these are referred to as model A and B (Table 3). For the given set of input parameters, model A gives rise to a J-type shock, and model B a C-type shock.

2.2 Molecular hydrogen

Collisional excitation and de-excitation of H₂ is calculated for collisions with H, H₂, and He. The collisional rate coefficients for H₂-H₂ collisions are adopted from Flower & Roueff (1998a) and for H₂–He collisions from Flower et al. (1998). In the case of H₂–H collisions, for the first 49 levels of H₂ the rates are from Flower (1997) and Flower & Roueff (1998b), where the rates have been calculated using a full quantum mechanical approach. For the remaining levels, the rates from Martin & Mandy (1995) are used. They were calculated using a quasi-classical approach. The reactive reaction rates of H₂ with H are from Le Bourlot et al. (1999).

The number of levels has been set to 150 here, and the highest level is υ = 8, J = 3 (E/k_B = 39 000 K). The model assumes that there are no levels between the user-set value and the dissociation level. This may be important when calculating the dissociation rate of H₂, since molecules that are already excited have internal energies that are closer to the dissociation limit, and thus require less energy to dissociate. For the models run here, we find that there is no significant difference in H₂ emission by increasing the number of levels.

Depending on the initial conditions, H₂ may dissociate in the shock through collisions. As the post-shock gas cools, H₂ reforms on the grains (Appendix A of Flower & Pineau des Forêts 2013) and it is necessary to account for the bond energy released (4.5 eV ~5.1 × 10⁴ K). We assume that approximately one third of the energy goes to internal energy of the molecule. This internal energy distribution follows a Boltzmann distribution with a temperature corresponding to ~17 000 K. The remaining energy is equally split between kinetic energy of the newly formed H₂ molecule, and heating of the grain.

The H₂ level populations are used for calculating the local H₂ line emissivities. This is done under the assumption of optically thin emission, which typically applies to H₂ emission because of its lack of a permanent dipole moment. Of these lines, 1000 are output explicitly and stored as emissivity profiles in this grid. About 900 of these H₂ lines are covered by the JWST instruments MIRI and NIRSpec. These two instruments together cover the wavelength range of 0.6–28 µm, that is the υ = 0–0 S(0) ground-state line at 28.3 µm (Fig. 1) is not covered.

Fig. 2 Illustration of the three steps required for running an externally irradiated shock model. The shock model shown here has a preshock density of 104 cm−3, shock velocity of 20 km s−1, and it is irradiated by a UV radiation field with G0 of 10. The strength of the transverse magnetic field is 100 µG. First a chemical steady-state model is run, and the thermal, chemical, and excitation output is used as input for a PDR model. The output of the PDR model is then used as input for the shock model. The top row shows the temperature evolution across the model run, the middle row the abundances of H and H2, while the bottom row shows normalized populations of the first five rotational levels of H2. The “bump” in the temperature profile at t ~ 102 yr in the chemical steady-state model comes from reformation of a small fraction of H2 on the grain, and the release of its binding energy.

2.3 Grid

The total set of grid parameters is presented in Table 2; covering this range of parameter space resulted in ~ 14 000 simulations in total. Each simulation produces a number of outputs that are all stored in human-readable ASCII files and an HDF5 file for easy extraction⁴. These include physical properties of the shock (e.g., temperature, density, velocity) as a function of distance and time through the shock, and chemical properties (e.g., local densities, charge state, column densities), excitation of H₂ (level populations and local emissivities). In this case, the time is calculated as the neutral flow time, t_n = ∫ d_z/υ_n. In total, more than 2600 quantities are stored as profiles through each shock, and 1400 quantities are stored as integrated values.

The model integrates the gas state far downstream in order to ensure that a steady-state solution is contained within the simulation. Therefore, special care needs to be taken when extracting integrated quantities such as column densities or line intensities. We here adopt a similar criterion for the size of the shock as in Godard et al. (2019) based on radiative energy dissipation. We here set that limit as the point where 99.9% of the total radiation has been emitted (see Appendix B). Specifically, this means that the size, z_s is defined as: (1)

where is the sum of the kinetic, magnetic, and thermal energy fluxes.

For ease of use, we provide a number of tables containing already-extracted results at the Centre de Données astronomiques de Strasbourg (CDS). Example tables are provided in Appendix B in Tables B.1–B.7. These tables include:

B.1 Physical parameters such as peak temperature, density, width, and age of the shock;

B.2 Column densities of selected species, particularly H, H₂, O, OH, H₂,C⁺,C, and CO;

B.3 Data required for creating H₂ excitation diagrams, i.e., ln(N/g) and E for each of the 150 levels;

B.4 H₂ integrated intensities of the 1000 lines extracted, along with their wavelength;

B.5 Width of the H₂ emitting zone for the υ = 0−0 S(1), 1−0 S(1), 0−0 S(9), 1−0 O(5), and 2−1 S(1) lines;

B.6 H₂ o/p ratios determined both locally and integrated through the shock;

B.7 Integrated line intensities of 29 transitions arising from C⁺, Si⁺, H, C, Si, O, S⁺, N⁺, N, and S.

On occasion, the model does not converge for numerical reasons; this happens in ~5% of cases. This convergence-failure occurs often in C^*-type shocks, when the flow crosses the first sonic point (see Appendix C in Godard et al. 2019). In these cases, the model output is ignored but the input parameters are still recorded in the tables.

2.4 Model limitations

The model has a number of inherent assumptions, which are discussed in the following. This includes the shock geometry, magnetic field orientation, self-irradiation, stationary shocks, and grain chemistry.

Geometry

The model treats a plane-parallel shock front, thus ignoring geometry. The lack of geometry is especially important in J-type shocks, where the gas may be compressed by four orders of magnitude or more. In nature, such a compression would quickly lead to an expansion of the high-pressure post-shock gas into the surrounding low-pressure medium, however, that is not possible in a 1D simulation. As a result, the post-shock density could be overestimated. For the case of H₂ emission, this is less important: most of the H₂ emission is generated in the warm parts of the shock where T > 100 K, prior to where significant post-shock expansion would occur.

Magnetic field orientation

The magnetic field orientation is assumed to be perpendicular to the direction of motion. This may not always be the case in molecular clouds, in fact, there is no a priori reason to assume the shock wave and field orientation are well aligned. If the field is not perpendicular to the direction of motion, the compression will lead to a change in field geometry, as described and discussed in Lehmann & Wardle (2016). These effects are not included here.

Self-irradiation

The model is best suited for molecular shocks. In shocks where H₂ is dissociated and atomic H is excited, the shocks become self-irradiated. While this self-irradiation can be solved iteratively (Lehmann et al. 2020, 2022), it is not included in the present version of the grid. This limits J-type shocks to υ_s ≲ 30 km s⁻¹.

Fig. 3 Energetic reprocessing for four shocks with b = 0.1. The pie charts show the percentage of energy lost relative to the input kinetic energy flux. The kinetic energy flux is primarily converted to heat, which goes to exciting the atoms and molecules that then radiate the energy away. This radiation is either from H2 (rotational and vibrational emission), other molecules (primarily CO, OH and H2O), or atoms (primarily H, O, and S). Some kinetic energy goes into compressing the magnetic field (“mag”), dissociating H2 collisionally (“H2 chem”), or atoms/molecules thermalizing with grains (“grain”). The percentages are shown in each pie slice, and the input shock parameters inside the pies. The input parameters all result in the shocks being J-type shocks, and model A is marked.

Stationary shocks

All the shocks in this paper are stationary shocks. This implies there needs to be enough time for the stationary structure to fully develop. While the code can mimic non-stationary shocks, an additional free parameter, the age of the shock, is needed, and it is deemed beyond the scope of this work to explore the effects of that parameter (e.g., Lesaffre et al. 2004a,b; Gusdorf et al. 2008).

Grain chemistry

Grain-grain interactions are omitted in this grid. For conditions where the velocity is below ~25 km s⁻¹ and the density is below ~10⁵ cm⁻³, this assumption is likely valid (Guillet et al. 2009, 2011). At larger velocities or densities, grains may interact, leading to grain evaporation and fragmentation which changes the size distribution of grains. Finally, in this grid we do not include ice mantles on the grains.

3 Results and discussion

The shock has an initial kinetic energy flux of 1/2 , where ρ = 1.4 n_H m_H is the mass density; most of this energy is radiated away in the shock. Figure 3 shows how the energy is lost in shocks with b = 0.1, velocities of 20 and 30 km s⁻¹, and densities of 10⁴ and 10⁶ cm⁻³. The pie charts are sorted by initial kinetic energy flux going from left to right, and top to bottom. The H₂ fraction decreases with increasing velocity and density because of dissociation. H₂ then reforms on the grains in the postshock gas introducing a heating term which counteracts the cooling of H₂. This is visible in the pie charts as the fraction of H₂ emission decreases monotonically with input kinetic energy flux, from 75% to 0.5%.

Figure 4 is similar to Fig. 3, but for a stronger magnetic field (b = 1.0), i.e., the input kinetic energy fluxes are the same as above. Increasing b to 1 has the consequence that the two 20-km s⁻¹ shocks become C-type shocks; the 30-km s⁻¹ shocks remain J-type shocks. The J-type shocks are dissociative, and the H₂ cooling fraction thus decreases significantly, as also illustrated in Fig. 3.

The distribution of energy flux into emission lines has been described previously (e.g., Kaufman & Neufeld 1996; Flower & Pineau des Forêts 2010, 2015; Lehmann et al. 2020), and a comparison in H₂ cooling fractions of the total input kinetic energy flux reveals broad agreement between different models and previous versions of the Paris-Durham model. These pie charts provide a global view of the energetic reprocessing in these shocks. In the following, the role of the different input parameters on the energetic reprocessing will be discussed in more detail, with a specific emphasis on H₂ emission.

Fig. 4 As Fig. 3 but for b = 1.0. For this change in b, shock models represented in the two left-most pie-charts are C-type, while the two right pie charts are J-type shocks; model B is marked.

3.1 Magnetic field

The strength of the transverse magnetic field, B, sets the ionmagnetosonic speed, c_ims, together with the ion mass density, ρ_i: (2)

where c_s is the sound speed. For υ_s < c_ims, the ionized and neutral fluids are decoupled and a magnetic precursor is present (Mullan 1971; Draine 1980); the code treats these multiple fluids self-consistently. For υ_s > c_ims, the ionized and neutral fluids are coupled, and there is no magnetic precursor (Fig. 5). We refer to Sect. 2.1 of Lehmann et al. (2022) for a more in-depth description of the differences between J- and C-type shocks. Figure 5 shows where the different shock types are as a function of b and υ_s for a density of 10⁴ cm⁻³, Fig. 6 shows the shock type for a part of the grid presented in this paper. For low values of b (≲0.3), the resulting shocks are J-type shocks, while for b ≳ 1.0 the resulting shocks are predominantly C-type shocks.

The effects of the magnetic precursor is that the input kinetic energy flux is deposited over a much larger spatial range (Fig. 5), resulting in lower peak temperatures when compared to shocks with the same input kinetic energy flux but no magnetic precursor. This naturally affects the excitation of H₂, as illustrated in Fig. 7 in the form of the fraction of total integrated intensity to initial kinetic energy flux. The H₂ excitation is illustrated for the two reference shocks (Table 3), both with the same input kinetic energy flux. The figure demonstrates that for both shocks, most of the kinetic energy is radiated away in H₂ emission (see Figs. 3 and 4); the difference in total H₂ integrated intensity from the two shocks is ~15%. However, the integrated intensity from model B (b = 1.0) is dominated by pure rotational emission (>99% of H₂ emission), whereas it is spread over the vibrational levels in model A (b = 0.1).

The differences in H₂ excitation and the origin thereof for different values of b are further explored in Fig. 8 for models A and B in the left and right column, respectively. The first row shows the emerging H₂ spectrum from the two shocks. As was already clear from Fig. 7, most of the H₂ emission in model A is spread over the vibrational transitions, whereas emission in model B predominantly is rotational. To make these artificial spectra, a uniform resolving power of R = λ/Δλ = 2500 is assumed, similar to the resolving powers of the NIRSpec and MIRI instruments on JWST, and the line shapes are Gaussian. That is, the integrated intensity calculated in the models is . A uniform resolving power implies that the emission from longer-wavelength transitions is spread over a larger wavelength range, and thus the peak emission is lower. This stark difference in the H₂ spectra can be understood from the physical structure of the shock.

The kinetic energy flux injected into the two shocks is the same, but the temperature structure is very different. For J-type shocks, such as model A, the maximum temperature can be approximated by (Lesaffre et al. 2013): (3)

For model A, the maximum temperature is ~2 × 10⁴ K (Fig. 8, second row). This high temperature ensures that the vibrational H₂ levels are readily populated. For model B (b = 1.0), on the other hand, the magnetic precursor causes the kinetic energy to be deposited over a much larger scale (~10³ AU vs. ~1 AU), and the resulting peak temperature is much lower (~2000 K). In this case, the temperature is so low that only the rotational levels are significantly excited.

The third row of Fig. 8 shows excitation diagrams for the two shocks. For model A, all points fall on a single curved line, indicating that the levels are probing a range of excitation temperatures, T_ex. Particularly, the higher-J and rovibrational transitions probe hotter gas than the lower-J transitions, and the slope is thus shallower (slope = −1/T_ex). In this case, the excitation temperatures is similar to the gas temperature where the local emissivity peaks (second row of Fig. 8). The excitation diagram for model B shows more scatter (caused by the low initial o/p ratio, see below), but the excitation temperatures still match the gas kinetic temperature where the levels are excited. In Appendix C.1 we provide figures showing the extracted excitation temperatures sampling the full range of initial density and shock velocity for b = 0.1 and 1.0, and G₀ = 0 and 1.

Another feature of the excitation diagram for model B is that there is a clear difference between the ortho- and para-levels of H₂. Here the ortho-levels (odd J) are displaced downward compared to the corresponding para-levels (even J), and the resulting zigzag pattern indicates that the ortho/para (o/p) ratio is lower than the high-temperature statistical equilibrium value of 3 (Neufeld et al. 2006).

There are no radiative or collisional transitions between ortho- and para-H₂ levels, only exchange reactions with H, H₂, and protonated ions (e.g., , HCO⁺) can change the spin state (Sect. 2.1 of Le Bourlot et al. 1999). The line emission and resulting excitation diagram is integrated through the shock, and thus does not provide information on the local o/p ratio. This is calculated directly from the level populations as n_o/n_p, and it can be compared to the cumulative column density ratio, N_o/N_p. Both these values are shown in the bottom row of Fig. 8. This column density ratio is often dominated by the column densities of H₂ in the two lowest rotational levels, J = 0 and 1, which are not accessible in emission. Therefore, we also show the o/p ratio as calculated from the column densities of the lowest observable rotational levels, in this case from the J = 2−9 levels (S(0) to S(7) transitions). In model A, the temperature is high enough that the H exchange reaction proceeds efficiently (e.g., Wilgenbus et al. 2000). The resulting o/p ratios are thus close to 3, although the inferred rotational o/p is somewhat lower than 3 (~1). For model B, the temperature never gets high enough that the exchange reactions with H become dominant; instead, the ion-neutral proton-transfer reactions dominate, but they are limited by the low abundances of ions. Thus, the o/p ratios remain at ~0.1. In both models, the initial temperature is 10 K and the gas is dense, which leads to a steady-state o/p ratio of 10⁻³ (see Fig. 1 of Flower et al. 2006). Had the initial temperature been higher or the gas not been in steady state, the initial o/p ratio would have been higher, and the o/p ratio through the shock also correspondingly higher. All in all, however, special care must be taken when interpreting o/p ratios inferred from observations (see also Fig. 4 of Wilgenbus et al. 2000).

As mentioned above, the input kinetic energy flux is deposited over a larger spatial range for increasing values of b. Specifically, a “phase transition” occurs when the resulting shock type goes from being J- to C-type, and a magnetic precursor develops. This typically happens at higher values of b or lower velocities (Fig. 6 shows which physical conditions lead to which shock type). Naturally the ionization fraction also plays a role in setting the shock type (Eq. (2)), but the gas is primarily neutral for the conditions examined here, and effectively this fraction does not play a role here. To measure the width and to make it a usable observational constraint, we have extracted the scale over which 80% of the H₂ emissivity is generated for a subset of lines: the υ = 0−0 S(1), 1−0 S(1), 0−0 S(9), 1−0 O(5), and 2−1 S(1) lines. These widths are shown in Fig. 9 together with the integrated intensity of the lines; here we show the widths of the υ = 0−0 S(1) and 1−0 S(1) emitting regions. The shocks with b = 0.1 all have widths less than 10 AU, whereas the b = 1 shocks have widths up to ~ 10⁵ AU or ~ 1 pc. For these shocks, there is an anticorrelation between the width and the integrated intensity: the wider shocks have lower integrated intensities. The J-type shocks occurring for b =1 and υ_s > 25 km s⁻¹ have larger widths than their b = 0.1 counterparts by one order of magnitude. Even though these are J-type shocks, the magnetic field still plays a significant role.

Fig. 5 Illustration of shock parameter space. Top: various shock type regimes as a function of transverse magnetic field strength and shock velocity for G0 = 0 and preshock density of 104 cm−3 (adapted from Fig. 3, Godard et al. 2019). Here a field strength of 10 µG corresponds to b = 0.1, and 103 µG is b = 10. The blue region is for J-type shocks where the velocity exceeds the ion-magnetosonic velocity (cims). Above 30 km s−1 the shocks start producing UV-photons and become self-irradiated, which is not included in the current grid (Lehmann et al. 2022). The red, green and orange areas are for CJ, C*, and C-type shocks, respectively; these shocks fall between the ion- and neutral-magnetosonic velocities (cims and cnms, where the latter is calculated in a similar manner to cims, Eq. (2), but with the neutral mass density). Models A and B are marked. Middle: gas temperature profiles for models A and B. Bottom: velocity profiles for models A and B. The dotted line shows the local sound speed, cs, the dashed is for the ion speed, υi, and the full is for the neutral speed, υn. All velocities are in the reference frame of the shock, and the line colors are for the same b parameters as above. The difference between the left- and right-hand sides of the middle and bottom panels is that the left panel is on a linear-linear scale, while the right is on a log-log scale.

Fig. 6 Resulting shock type as a function of UV radiation field strength, preshock density, shock velocity, and magnetic field strength for a subset of the grid (for brevity, the shock velocities 2, 3, 4, 15, 25, and 50 km s−1 are not shown). The cell color denotes the shock type, where blue is for J-type shocks, orange is for C-type shocks, red is for CJ-type shocks, and green is for C*-type shocks, as is also written in each cell. The figure shows model results for ζH2 = 10−17 s−1, X(PAH) = 10−6. White cells are for models that did not converge numerically.

Fig. 7 Top: H2 integrated intensities integrated over 4π steradian from each vibrational level compared to the input kinetic energy flux for model A and B. The total H2 intensity is similar in the two models, 5.5 × 10−3 and 4.7 × 10−3 erg s−1 cm−2 sr−1 for model A and B, respectively. Bottom: As the top panel, but for densities of 102,104, and 106 cm−3, b = 1.0, and υs = 20 km s−1. The maximum neutral gas temperature is given in brackets.

3.2 Velocity and density

The shock velocity, υ_s sets the maximum temperature in J-type shocks (Eq. (3)). H₂ excitation is sensitive to temperature, and so the velocity effectively sets the excitation. This is seen in the simulated spectra (Fig. 10). At the lowest velocity (5 km s⁻¹), the integrated intensity is low and only a few rotational lines are seen in the spectrum. On the contrary, at velocities ≳20 km s⁻¹, we see rich vibrational H₂ spectra. At the same time the peak specific intensity increases by a factor of ~10, until the velocity reaches 30 km s⁻¹ and the shock becomes dissociative. In this case, H₂ only contributes to the cooling once it has reformed on the grains. Thus, to a first order, the excitation is set primarily by the velocity in J-type shocks, and the density plays a role in setting the total integrated intensity.

In C-type shocks, the combination of density and velocity is what affects the excitation and the integrated intensity (Fig. 7, bottom panel). This is illustrated in the top row of Fig. 11, which shows the total H₂ integrated intensity emitted as well as the brightest line. Here, the brightest line serves as a proxy for the excitation in the sense that the higher excited the brightest line is, the higher the excitation is. For the C-type shocks (orange dots), there is a clear intensity and excitation gradient which depends on both density and velocity. The brightest lines are rotational over the bulk of parameter space (from 0−0 S(0) to S(6)), and they are typically para-H₂ transitions (even J). For the case of J-type shocks (blue dots), the intensity gradient is dominated by the density, as discussed above. However, the brightest lines quickly become vibrational; the υ = 1−0 Q(1) line (2.41 µm) is predicted to be particularly bright, as is the υ = 1−0 S(3) line (1.96 µm). Thus, identifying the brightest line in the H₂ spectrum provides constraints on where in parameter space the shock is located. Appendix D provides an overview of the dominant cooling lines across the grid.

The H₂ fraction in the gas is highest at the lower densities and lower velocities where H₂ does not dissociate. However, for a given velocity, the total H₂ integrated intensity increases mono-tonically with density, as shown in Fig. 11. This is in spite of the fraction of input kinetic energy flux radiated by H₂ is mono-tonically decreasing. Thus, for the shocks with the brightest H₂ emission, other molecules and atoms are needed to trace the bulk deposition of kinetic energy. Examples include emission from CO and H₂O at lower velocities, and O, S, and H at higher velocities.

Fig. 8 Top row: resulting H2 spectrum in model A and B (Table 3). The wavelength ranges of NIRSpec and MIRI on the JWST are labeled. The resolving power is uniform across the spectrum, implying that the linewidths increase with wavelength (not visible on this plot). Second row: H2 emissivities for three lines. The temperature profile is shown in black. Third row: H2 excitation diagrams color-coded according to vibrational level, and with four extracted excitation temperatures shown. The excitation temperatures are obtained from a linear fit to all the indicated levels. Fourth row: o/p ratio measured from local densities (no/np), cumulative column densities (No/Np), and from the column densities of the υup = 0, J = 2–9 levels (No/Np(rot)).

Fig. 9 Width of the H2 emitting zone for b = 0.1 (left) and 1.0 (right) for the transitions υ = 0−0 S(1) (17.05 µm; top) and 1−0 S(1) (2.12 µm; bottom). The width is measured as the extent of the region where 80% of the emission is radiated away. In each cell is written the log of the integrated intensity of the line in units of erg cm−2 s−1 sr−1. The colored dots indicate the resulting shock type, where blue is for J-type shocks and orange for C-type shocks.

3.3 UV radiation field

In an externally UV-irradiated shock, the UV photons lead to increased gas ionization and thus higher density of the charged fluid. This increase causes a tighter coupling between the neutral and charged fluids, which in turn leads to the kinetic energy typically being deposited over shorter scales compared to in the absence of external UV radiation. Thus, the temperature typically increases and the shocks become narrower (see Fig. 6 in Godard et al. 2019). The increased temperature naturally causes higher excitation of H₂, as is illustrated in the H₂ spectra in Fig. 12. Here, the shock in model B, showing pure rotational excitation of H₂, is exposed to increasing strengths of an external UV-field, from G₀ = 0 to 10³. The increase in temperature (from 1700 K to 2800 K) leads to an increase in excitation, and the vibrational levels start to become populated.

The second effect of the UV field is to deposit additional energy into the shock (Fig. 12 in Godard et al. 2019). Either this energy deposition is indirect in the form of ionization followed by recombination and release of binding energy, or the energy deposition is direct, where UV photons excite H₂ electronically, from which the molecules can de-excite radiatively. It is clear that for the highest values of G₀, the additional energetic input is significant. This is illustrated in Fig. 13. Here, the energy radiated away by H₂ as a function of vibrational level is shown for model B, similar to Fig. 7. In this case, model B is exposed to stronger UV fields, and the higher vibrational levels are excited, as also seen in Fig. 12. The total fraction of energy lost in H₂ emission increases almost monotonically from 0.63 to 1.07 of the input kinetic energy flux. Thus, at least 7% of the excitation is caused by the UV field, and likely more as there are other channels of energy loss (Fig. 4). For a quantitative description of the role of UV pumping on the H₂ level populations, we refer to Fig. 8 of Godard et al. (2019).

Even for relatively weak UV field strengths (e.g., G₀ = 1), the UV photons may play a significant role. Figure 14 is similar to Fig. 11 in that the top panels show the total amount of H₂ emission and the strongest H₂ line. For the weak shocks (low density, low velocity), one major difference is seen when the UV field is turned on: in the absence of external UV radiation, the brightest lines are all para-H₂ lines (even J) because there is no significant para- to ortho-H₂ conversion. For the weak UV field, the strongest lines are predominantly ortholines (odd J), which is consistent with observations of the diffuse gas in colliding galaxies (Ingalls et al. 2011; Guillard et al. 2012; Appleton et al. 2017). This suggests that interstellar shocks in general are not fully shielded, but exposed to some UV radiation.

Fig. 10 H2 spectra for velocities of 5, 10, 20, and 30 km s−1 and a density of 104 cm−3, b = 0.1, G0 = 0; the shock with υS = 20 km s−1 is model A, and the 30 km-s−1 shock is shown in Fig. 1. The colors are for different vibrational levels as in Fig. 8. The complete coverage of NIRSpec and MIRI are shown; we refer to Fig. 1 for the NIRCam and MIRI filter coverage.

3.4 H₂ excitation for JWST observers

JWST represents an increase in sensitivity, spatial and spectral resolution by more than an order of magnitude over previous infrared space-based telescopes (Rigby et al. 2023). We here outline some of the ways in which the models may be used to plan and interpret the JWST observations of shocked regions, keeping in mind the model limitations listed in Sect. 2.4.

H₂ spectroscopy

The spectroscopic capabilities of NIR-Spec and MIRI make them perfectly suited for observing H₂ line emission. The excitation of H₂ is the result of a complex interplay between various input parameters, as discussed above, with some degeneracies, especially between the density and shock velocity. This is for example illustrated in Fig. 13 of Kristensen et al. (2007), where observations of H₂ emission from the explosive Orion-KL protostellar outflow are analyzed. With high enough spectral resolution, independent constraints can be made on the shock velocity, thus directly breaking the degeneracy (Santangelo et al. 2014).

It will likely not be possible to strongly constrain shock conditions from H₂ observations alone, unless the observers only consider subgrids of physical parameters relevant to their studies. An example could be that if shocks in diffuse clouds are studied, only the lowest densities in the grid would be relevant. Furthermore, in a large number of cases, G₀ can be independently constrained, for example, by studying ionized gas lines, UV continuum observations, or PAH features at infrared wavelengths. Observers should also be aware that, in shock-dominated environments, the total H₂ line emission in a given beam is likely the product of a distribution of shocks arising from a multiphase medium with different conditions. Such an example of shock probability distributions convolved with the use of grids of shock models have been used to interpret H₂ observations in the intra-group shocked diffuse gas in colliding galaxies (e.g., Guillard et al. 2009; Lesaffre et al. 2013).

Fig. 11 Top row: total H2 integrated intensity in two subgrids of models with b = 0.1 (left) and 1.0 (right), G0 = 0. The colored dots indicate the resulting shock type, where blue is for J-type shocks and orange for C-type shocks (the single green point is a C* shock). The brightest line is written in each cell. Bottom row: as the top row, but for the total H2 integrated intensity radiated away compared to the input kinetic energy flux. The actual percentage is provided in each cell.

Shock width

The NIRCam instrument on JWST is well-suited for observing H₂ emission. The instrument contains three categories of filters, narrow-, medium-, and wide-band filters. Their wavelength coverages are illustrated in Fig. 1. Of the narrowband filters, three center on H₂ lines: F212N (υ = 1−0 S(1)), F323N (υ = 1−0 O(5)), and F470N (υ = 0−0 S(9)). The spatial resolution ranges from 0″.07 to 0″.16, corresponding to linear scales of 14 and 32 AU at a distance of 200 pc, a typical distance to nearby star-forming regions. As illustrated in Fig. 9, the width of shocks with b = 1.0 is typically resolvable if the shock is observed close to edge on, except at the highest densities (≳10⁷ cm⁻³ for C-type shocks, and ≳10⁶ cm⁻³ for J-type shocks). Shocks with b = 0.1 are not resolvable at a distance of 200 pc. Having a measured shock width puts additional constraints on the shock models: the width is sensitive to the strength of the transverse magnetic field and thus serves as an independent constraint of this parameter (Figs. 8 and 9 in Kristensen et al. 2008, for observations of a spatially resolved bow shock in Orion-KL). Besides NIRCam, the MIRI IFU offers the possibility of producing spectral line maps of H₂ emission at 160 AU (0″.5) spatial resolution at a distance of 200 pc of the 0−0 S(1) line at 17 µm. Emission from this line traces colder gas, and so is typically more extended than the higher-excited lines shown in Fig. 9. This resolution is therefore still enough to resolve shock-dominated line emission from dissipative regions in nearby star-forming clouds (Richard et al. 2022).

Fig. 12 H2 spectra for model B but with G0 of 0 (Model B), 1, and 103. The colors are for different vibrational levels as in Fig. 8. The complete coverage of NIRSpec and MIRI are shown; we refer to Fig. 1 for the NIRCam and MIRI filter coverage.

H₂ photometry

As shown in Fig. 1, the NIRCAM and MIRI imaging filters includes multiple ro-vibrational and rotational H₂ lines, so the use of a those filters may prove to be efficient as far as exposure time and mapping area are concerned. Such observations may be used for constraining shock conditions. As an example, Figs. 11 and 14 show the brightest lines for a given set of initial conditions. Thus, if an observed region is dominated by shocked H₂ emission, then it might be possible to broadly constrain the range of parameter space where the emission is generated. That is, with the model results in hand, the user can construct “H₂ photometry” which can be compared to observations, assuming H₂ emission dominates the spectrum and the contribution from, e.g., PAH emission is negligible, or assuming that a combination of filters can be used to remove the contribution of the continuum emission. A similar approach has been shown to work efficiently for the wideband MIRI filters for observations of the colliding galaxies in Stephan’s Quintet (Appleton et al. 2023).

Fig. 13 Emission distribution for model B with increasing G0 (similar to Fig. 7). The total fraction of H2 emission to input kinetic energy flux is provided in the legend.

H₂ summary

Table 4 summarizes what sets the H₂ integrated intensity and the excitation. This table is by no means exhaustive, but may be used as an overview guide of H₂ emission in shocks. To constrain the excitation properly, it is necessary to cover as large a wavelength range as possible, and to cover both rotational and rovibrational lines. The former are predominantly excited in C-type shocks, and the latter in J-type shocks. Once a solution has been found that approximately reproduces observations, we recommend the user to fine-tune the grid further for more precise solutions. This can be done either by interpolating the grid values; in this case care must be taken when going from one shock type to another. Alternatively the user can download the model and run their own shock models, in which case we recommend benchmarking their results against the models presented here in a first step. Finally, we recommend that the total integrated intensity of the H₂ lines is compared to the total available mechanical energy output from a given source, to ensure that the best-fit shock model is physical (see Lehmann et al. 2022, for the methodology).

Atomic lines

Apart from H₂ emission, the model calculates line emission from several other atomic and ionic species. As an example, JWST-MIRI will observe the [S I] line at 25 µm (e.g., toward the nearby protostellar outflow from IRAS15398, Yang et al. 2022), and the integrated line intensity of this line is calculated and tabulated from the grid. The same applies to lines from other species, e.g., O, and C. Naturally, these lines light up in different parts of parameter space compared to H₂, and thus provide complementary information.

Other emission lines

The abundances of some 140 other species have been calculated through the shock. Examples of particular relevance to JWST and shocks include Fe⁺, OH and H₂O, because these species have a number of transitions visible in the NIRSpec and MIRI wavelength ranges and these species are some of the dominant coolants (e.g., Figs. 3 and 4). The abundance, temperature, and density profiles are calculated through the shock, which means that the profiles can be post-processed to calculate integrated line intensities using for example a large velocity gradient (LVG) radiative transfer code (e.g., Gusdorf et al. 2011), which has not been done for this grid of models. Just as for the atomic lines, these will provide complementary observational constraints.

Fig. 14 As Fig. 11, but for G0 = 1.

4 Summary

Here we present the results of an extensive grid of planeparallel steady-state shock models. The grid was constructed by varying six parameters: the preshock density, shock velocity, strength of the transverse magnetic field, strength of the UV field impinging on the shock, the cosmic-ray-ionization rate, and the PAH abundance. This is the first time such an extensive grid of shock models has been run and made publicly available.

The purpose of running this grid of models was to examine under which shock conditions H₂ is efficiently excited, and how shock conditions affect the H₂ excitation and integrated line intensities. H₂ is already being extensively observed with JWST, and the coming years will see a flood of H₂ observations. At the moment it is therefore critical for planning and interpreting JWST observations.

We find that the strength of the transverse magnetic field, as quantified by the magnetic scaling factor, b, plays a key role in the excitation of H₂. At low values of b (≲0.3, J-type shocks), H₂ excitation is dominated by vibrationally excited lines; whereas, at higher values (b ≳ 1, C-type shocks), rotational lines dominate the spectrum for shocks without an external radiation field. Shocks with b ≥ 1 can potentially be spatially resolved with JWST for nearby objects, which serves as an additional constraint.

H₂ is typically the dominant coolant at lower densities (≲10⁴ cm⁻³); at higher densities, other molecules such as CO, OH, and H₂O take over at velocities ≲20 km s⁻¹ and atoms, for example, H, O, and S, dominate at higher velocities. Together, the velocity and density set the input kinetic energy flux. When this increases, the excitation and integrated intensity of H₂ increases similarly.

An external UV field mainly serves to increase the excitation, particularly for shocks where the input radiation energy is comparable to or greater than the input kinetic energy flux. Together, these results provide an overview of the energetic reprocessing of input energy and the resulting H₂ line emission observable by JWST.

Acknowledgements

We would like to thank F. Boulanger and S. Cabrit for simulating discussions, particularly at the beginning of this project, as well as J. A. Villa Vélez. The research leading to these results has received funding 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 simulations used in this work has been run on the computing cluster Totoro of the ERC MIST, administered by MesoPSL. 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. The research of LEK is supported by a research grant (19127) from VILLUM FONDEN. PG would like to thank the Sorbonne University, the Institut Universitaire de France, the Centre National d’Etudes Spatiales (CNES), the “Programme National de Cosmologie and Galaxies” (PNCG). This work has made use of the Paris-Durham public shock code V1.1, distributed by the CNRS-INSU National Service “ISM Platform” at the Paris Observatory Data Center⁵.

Appendix A The ISM platform

The ISM platform⁶ is a web portal that contains a series of services developed for the diffusion of state-of-the-art astro-chemical models and the preparation and interpretation of observations. Regarding the Paris-Durham shock code, the platform provides access to the numerical code and its previous versions, a full documentation of the physical processes implemented, a tutorial to learn how to run the code locally, and a series of selected references. The platform also provides two analysis tools, IDAT and the Chemistry Analyzer tool, which can be used to study the output of the shock code and identify the processes responsible for the thermochemical evolution of the gas in a simulation. Finally, the platform contains a numerical database (InterStellar Medium DataBase or ISMDB) that provides an easy access to recalculated grid of theoretical models.

On this platform it is possible to “Search models in ISMDB” and from there “Browse models.” This leads to a page where combinations of input shock parameters can be specified, and once the selection has been made, it is possible to “Get model.” The resulting page shows the input parameters as well as some of the resulting quantities (e.g., shock type). The entire model output can be downloaded for further analysis, or the model can be quickly inspected directly through “Online analysis with IDAT.” This tool allows the user to select different quantities and plot them against distance through the shock on one or two different y-axes if so desired. An example could be the velocities through the shock as well as the temperature.

Appendix B Tables with extracted parameters

We here provide example tables of the physical quantities already extracted from the grid (Tables B.1–B.7). These tables are available on CDS in electronic format. These tables include:

B.1 Physical quantities such as peak temperature, density, width, and age of the shock;

B.2 Column densities of relevant species, particularly H, H₂, O, OH, H₂,C+,C, and CO;

B.3 Data required for creating H₂ excitation diagrams, i.e., ln(N/g) and E for each of the 150 levels;

B.4 H₂ integrated intensities of the 1000 lines extracted, along with their wavelength;

B.5 Width of the H₂ emitting zone for the υ = 0−0 S(1), 1−0 S(1), 0−0 S(9), 1−0 O(5), and 2−1 S(1) lines;

B.6 H₂ o/p ratios determined both locally and integrated through the shock;

B.7 Integrated line intensities of 29 transitions arising from C⁺, Si⁺, H, C, Si, O, S⁺, N⁺, N, and S.

An energy cutoff of 99.9% was used to define the point at which integrated quantities (e.g., line intensities, column densities) were integrated to (Sect. 2.3). Tests were performed using cutoffs at 95%, 99%, 99.9%, 99.99%, and 99.999%. The two lower values (95 and 99%) did not capture the H₂-emitting zone, particularly in strong CJ-type shocks where the temperature exceeds 10⁵ K. The difference between 99.9% and 99.99% cutoffs were on the order of a few percent in terms of H₂ integrated line intensities for the υ = 0−0 S(1), 1−0 S(1), and 2−1 S(1) transitions for most shock conditions. Thus, a threshold of 99.9% ensured that most of the H₂ radiative cooling zone was encompassed.

Appendix C Additional figures

Appendix C.1 Excitation temperatures

Excitation temperatures have been extracted and calculated from a subset of the grid. Figures C.1 and C.2 show these temperatures calculated from the υ = 0, J = 3 to 5 levels (S(1) to S(3)) and the υ = 0, J = 6 to 11 levels (S(4) to S(9)) levels, respectively. The excitation temperatures are shown for b = 0.1 and 1, and G₀ = 0 and 1. Figures C.3 and C.4 show excitation temperatures for the υ =1, J = 0–8 and υ = 2, J = 0–8 vibrationally excited levels.

Appendix C.2 Cosmic ray ionization rate

In the model, cosmic rays may ionize H₂ and other species. When these species recombine, primarily H₂, secondary UV photons are emitted. Direct excitation by cosmic rays is not included. In this manner, cosmic rays serve as an additional source of both ionization and thus energy input. The expectation is that they will impact the H₂ emission to a similar degree as an external UV field. Their impact, however, is smaller than that of UV radiation. This is illustrated in Fig. C.5, where the integrated line intensity of three representative lines are shown as a function of the cosmic ray ionization rate, ζ_H2, for Model B. In this case, the PAH abundance is set to 10⁻⁸. For no external UV radiation, the integrated intensity increases by ~ one order of magnitude when ζ_H2 increases by two orders of magnitude. For G₀ = 1, there is practically no change in intensity over the same range of ζ_H2, however, the vibrationally excited lines are significantly brighter than for the shocks without an external radiation field.

Fig. C.1 Excitation temperature determined from the υ = 0, J = 3 to 5 level populations (corresponding to the S(1) to S(3) transitions) for b = 0.1 and 1, and G0 = 0 and 1.

Fig. C.2 Excitation temperature determined from the υ = 0, J = 6 to 11 level populations (corresponding to the S(4) to S(9) transitions) for b = 0.1 and 1, and G0 = 0 and 1.

Fig. C.3 Excitation temperature determined from the υ =1, J = 0 to 8 level populations for b = 0.1 and 1, and G0 = 0 and 1.

Fig. C.4 Excitation temperature determined from the υ = 2, J = 0 to 8 level populations for b = 0.1 and 1, and G0 = 0 and 1.

Fig. C.5 Integrated intensity of select H2 lines as a function of the cosmic-ray ionization rate, ζH2. The panel on the left is for G0 = 0, and on the right for G0 = 1. Both are for a C-type shock with shock velocity of 20 km s−1, density of 104 cm−3, and b = 1.0. The PAH abundance is 10−8.

Appendix D Dominant cooling lines

It is natural, when examining such a large grid, to identify the dominant H₂ cooling lines, that is, the H₂ lines that are most likely to be observed for a given set of input parameters. One way of identifying these lines for the entire grid, is to go through each model and tabulate the lines with integrated intensities that are greater than 25% of the maximum intensity. This arbitrary cutoff is chosen from the perspective that if the strongest line is detected at 20σ, then these lines would also be detectable at the 5σ level. Next, the lines are sorted according to which ones are present in the largest number of models, i.e., which are typically the dominant cooling lines in a global perspective. The lines that are present in at least 25% of models are tabulated in Table D.1.

Twenty-four lines are present in at least 25% of models. The lines are either υ = 0−0 or 1−0 transitions; the higher-excited levels are clearly not sufficiently populated over the majority of the grid. Some of the lines in Table D.1 are observable from the ground, for example, the often bright υ = 1−0 S(1) line at 2.12 µm, but the majority of the lines are not (17/24 lines). All lines are, however, observable with the JWST. Eighteen lines are observable with NIRSpec, while seven are observable with MIRI. At 5.06 µm, the υ = 0−0 S(8) line is observable with both instruments, and could serve as a cross-calibrator between the two instruments.

References

All Tables

All Figures

Fig. 1 Synthetic H2 spectrum produced with a shock model with velocity 30 km s−1, preshock density 104 cm−3, a transverse magnetic field strength of 10 µG, and no external UV radiation. Wavelength ranges of the NIRSpec and MIRI spectrographs, as well as the wide-, medium-, and narrow-band filters for NIRCam and the MIRI filters on JWST are indicated as black and gray horizontal bars. The colors are for lines with different vibrational upper levels. The resolving power is assumed to be uniform across the wavelength range at λ/Δλ = 2500.

In the text

Fig. 2 Illustration of the three steps required for running an externally irradiated shock model. The shock model shown here has a preshock density of 104 cm−3, shock velocity of 20 km s−1, and it is irradiated by a UV radiation field with G0 of 10. The strength of the transverse magnetic field is 100 µG. First a chemical steady-state model is run, and the thermal, chemical, and excitation output is used as input for a PDR model. The output of the PDR model is then used as input for the shock model. The top row shows the temperature evolution across the model run, the middle row the abundances of H and H2, while the bottom row shows normalized populations of the first five rotational levels of H2. The “bump” in the temperature profile at t ~ 102 yr in the chemical steady-state model comes from reformation of a small fraction of H2 on the grain, and the release of its binding energy.

In the text

Fig. 3 Energetic reprocessing for four shocks with b = 0.1. The pie charts show the percentage of energy lost relative to the input kinetic energy flux. The kinetic energy flux is primarily converted to heat, which goes to exciting the atoms and molecules that then radiate the energy away. This radiation is either from H2 (rotational and vibrational emission), other molecules (primarily CO, OH and H2O), or atoms (primarily H, O, and S). Some kinetic energy goes into compressing the magnetic field (“mag”), dissociating H2 collisionally (“H2 chem”), or atoms/molecules thermalizing with grains (“grain”). The percentages are shown in each pie slice, and the input shock parameters inside the pies. The input parameters all result in the shocks being J-type shocks, and model A is marked.

In the text

	Fig. 4 As Fig. 3 but for b = 1.0. For this change in b, shock models represented in the two left-most pie-charts are C-type, while the two right pie charts are J-type shocks; model B is marked.
In the text

Fig. 5 Illustration of shock parameter space. Top: various shock type regimes as a function of transverse magnetic field strength and shock velocity for G0 = 0 and preshock density of 104 cm−3 (adapted from Fig. 3, Godard et al. 2019). Here a field strength of 10 µG corresponds to b = 0.1, and 103 µG is b = 10. The blue region is for J-type shocks where the velocity exceeds the ion-magnetosonic velocity (cims). Above 30 km s−1 the shocks start producing UV-photons and become self-irradiated, which is not included in the current grid (Lehmann et al. 2022). The red, green and orange areas are for CJ, C*, and C-type shocks, respectively; these shocks fall between the ion- and neutral-magnetosonic velocities (cims and cnms, where the latter is calculated in a similar manner to cims, Eq. (2), but with the neutral mass density). Models A and B are marked. Middle: gas temperature profiles for models A and B. Bottom: velocity profiles for models A and B. The dotted line shows the local sound speed, cs, the dashed is for the ion speed, υi, and the full is for the neutral speed, υn. All velocities are in the reference frame of the shock, and the line colors are for the same b parameters as above. The difference between the left- and right-hand sides of the middle and bottom panels is that the left panel is on a linear-linear scale, while the right is on a log-log scale.

In the text

Fig. 6 Resulting shock type as a function of UV radiation field strength, preshock density, shock velocity, and magnetic field strength for a subset of the grid (for brevity, the shock velocities 2, 3, 4, 15, 25, and 50 km s−1 are not shown). The cell color denotes the shock type, where blue is for J-type shocks, orange is for C-type shocks, red is for CJ-type shocks, and green is for C*-type shocks, as is also written in each cell. The figure shows model results for ζH2 = 10−17 s−1, X(PAH) = 10−6. White cells are for models that did not converge numerically.

In the text

Fig. 7 Top: H2 integrated intensities integrated over 4π steradian from each vibrational level compared to the input kinetic energy flux for model A and B. The total H2 intensity is similar in the two models, 5.5 × 10−3 and 4.7 × 10−3 erg s−1 cm−2 sr−1 for model A and B, respectively. Bottom: As the top panel, but for densities of 102,104, and 106 cm−3, b = 1.0, and υs = 20 km s−1. The maximum neutral gas temperature is given in brackets.

In the text

Fig. 8 Top row: resulting H2 spectrum in model A and B (Table 3). The wavelength ranges of NIRSpec and MIRI on the JWST are labeled. The resolving power is uniform across the spectrum, implying that the linewidths increase with wavelength (not visible on this plot). Second row: H2 emissivities for three lines. The temperature profile is shown in black. Third row: H2 excitation diagrams color-coded according to vibrational level, and with four extracted excitation temperatures shown. The excitation temperatures are obtained from a linear fit to all the indicated levels. Fourth row: o/p ratio measured from local densities (no/np), cumulative column densities (No/Np), and from the column densities of the υup = 0, J = 2–9 levels (No/Np(rot)).

In the text

Fig. 9 Width of the H2 emitting zone for b = 0.1 (left) and 1.0 (right) for the transitions υ = 0−0 S(1) (17.05 µm; top) and 1−0 S(1) (2.12 µm; bottom). The width is measured as the extent of the region where 80% of the emission is radiated away. In each cell is written the log of the integrated intensity of the line in units of erg cm−2 s−1 sr−1. The colored dots indicate the resulting shock type, where blue is for J-type shocks and orange for C-type shocks.

In the text

	Fig. 10 H2 spectra for velocities of 5, 10, 20, and 30 km s−1 and a density of 104 cm−3, b = 0.1, G0 = 0; the shock with υS = 20 km s−1 is model A, and the 30 km-s−1 shock is shown in Fig. 1. The colors are for different vibrational levels as in Fig. 8. The complete coverage of NIRSpec and MIRI are shown; we refer to Fig. 1 for the NIRCam and MIRI filter coverage.
In the text

Fig. 11 Top row: total H2 integrated intensity in two subgrids of models with b = 0.1 (left) and 1.0 (right), G0 = 0. The colored dots indicate the resulting shock type, where blue is for J-type shocks and orange for C-type shocks (the single green point is a C* shock). The brightest line is written in each cell. Bottom row: as the top row, but for the total H2 integrated intensity radiated away compared to the input kinetic energy flux. The actual percentage is provided in each cell.

In the text

	Fig. 12 H2 spectra for model B but with G0 of 0 (Model B), 1, and 103. The colors are for different vibrational levels as in Fig. 8. The complete coverage of NIRSpec and MIRI are shown; we refer to Fig. 1 for the NIRCam and MIRI filter coverage.
In the text

	Fig. 13 Emission distribution for model B with increasing G0 (similar to Fig. 7). The total fraction of H2 emission to input kinetic energy flux is provided in the legend.
In the text

	Fig. 14 As Fig. 11, but for G0 = 1.
In the text

	Fig. C.1 Excitation temperature determined from the υ = 0, J = 3 to 5 level populations (corresponding to the S(1) to S(3) transitions) for b = 0.1 and 1, and G0 = 0 and 1.
In the text

	Fig. C.2 Excitation temperature determined from the υ = 0, J = 6 to 11 level populations (corresponding to the S(4) to S(9) transitions) for b = 0.1 and 1, and G0 = 0 and 1.
In the text

	Fig. C.3 Excitation temperature determined from the υ =1, J = 0 to 8 level populations for b = 0.1 and 1, and G0 = 0 and 1.
In the text

	Fig. C.4 Excitation temperature determined from the υ = 2, J = 0 to 8 level populations for b = 0.1 and 1, and G0 = 0 and 1.
In the text

	Fig. C.5 Integrated intensity of select H2 lines as a function of the cosmic-ray ionization rate, ζH2. The panel on the left is for G0 = 0, and on the right for G0 = 1. Both are for a C-type shock with shock velocity of 20 km s−1, density of 104 cm−3, and b = 1.0. The PAH abundance is 10−8.
In the text