# Supernova Shocks in Molecular Clouds: Shocks Driven into Dense Cores in IC 443 and 3C 391

William Reach, Ngoc Le Tram, Curtis Dewitt, Pierre Lesaffre, Benjamin

Godard, Antoine Gusdorf

► To cite this version:

William Reach, Ngoc Le Tram, Curtis Dewitt, Pierre Lesaffre, Benjamin Godard, et al. Supernova Shocks in Molecular Clouds: Shocks Driven into Dense Cores in IC 443 and 3C 391. The Astrophysical Journal, 2024, 977 (2), pp.149. 10.3847/1538-4357/ad8d59. hal-04779786

HAL Id: hal-04779786

https://hal.science/hal-04779786v1

Submitted on 8 Jan 2025

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
<figure><img src="image_2.png" /></figure>

# Supernova Shocks in Molecular Clouds: Shocks Driven into Dense Cores in IC 443 and 3C 391

William T. Reach<fnref n="1" />, Le Ngoc Tram<fnref n="2" />,<fnref n="3" />, Curtis DeWitt<fnref n="1" />, Pierre Lesaffre<fnref n="4" />, Benjamin Godard<fnref n="4" />, and Antoine Gusdorf<fnref n="4" />,<fnref n="5" />

Received 2024 September 26; revised 2024 October 29; accepted 2024 October 29; published 2024 December 10

## Abstract

Supernova shocks into dense molecular cores in IC 443 (clumps B, C, and G) and 3C 391 were observed using the Stratospheric Observatory for Infrared Astronomy and complemented by archival data from the Herschel Space Observatory. The pure rotational transitions 0–0 S(1) and S(5) of H$_{2}$, and the ground-state 1$_{10}$–1$_{01}$ transition of H$_{2}$O, are all broadened, arising from molecules that survive the passage of the shock front. Theoretical models from the Paris–Durham shock code were analyzed to generate velocity profiles that approximately match the observations. The observations can be fit with two shock conditions, which approximate the range of densities in the preshock molecular cloud. The width and brightness of the S(5) lines require shocks into gas with a density of order 2000 cm$^{-3}$, into which the IC 443 blast wave drives shocks with speed 60 km s$^{-1}$. The brightness and narrower width of the S(1) lines requires different shocks, into gas with density of order 10$^{5}$ cm$^{-3}$, with shock speeds of 10 km s$^{-1}$. The H$_{2}$O velocity distribution is also consistent with these shocks. The existence of shocks into dense gas shows that the bright shocked clumps in IC 443 were prestellar cores. It is unlikely that they will form stars soon after the passage of the shock front, given the input of kinetic and thermal energy from the shocks.

*Unified Astronomy Thesaurus concepts: Supernova remnants (1667); Shocks (2086); Molecular clouds (1072); Infrared spectroscopy (2285)*

## 1. Introduction

The two supernova remnants (SNRs), IC 443 and 3C 391, present some of the most accessible sites of supernova-shocked molecular gas. The southern portion of IC 443 is impacting a molecular cloud along an extensive and clumpy ridge of shocked molecular gas first seen in millimeter-wave CO emission (L. K. Denoyer 1979; Y. L. Huang et al. 1986) and near-infrared rovibrational 1–0 S(1) emission (M. G. Burton et al. 1988). Submillimeter lines reveal a wealth of chemistry as well as dynamic changes across the shock ridge (E. F. van Dishoeck et al. 1993). A pure rotational line of H$_{2}$, 0–0 S(4), expected to be a strong coolant for many types of molecular shock, has been detected in ground-based observations (M. J. Richter et al. 1995), but lines from the lowest energy levels are inaccessible at ground level. The Infrared Space Observatory detected S(0) through S(8) lines from the molecular ridge in IC 443, and the results indicated young (1000–2000 yr) shocks into moderate-density (10$^{4}$ cm$^{-3}$) gas (D. Cesarsky et al. 1999). Spitzer detected and imaged H$_{2}$ emission lines (D. A. Neufeld et al. 2007; Y. Yuan & D. A. Neufeld 2011; P. Dell’Ova et al. 2024). The Submillimeter-Wave Astronomical Satellite detected H$_{2}$O emission from the shocked gas, with a remarkably wide velocity distribution (R. L. Snell et al. 2005). The SNR 3C 391 impinges upon a giant molecular cloud (D. J. Wilner et al. 1998), and shocked clumps are evident in infrared

(W. T. Reach & J. Rho 1996, 1998) and millimeter waves (W. T. Reach & J. Rho 1999; D. A. Neufeld et al. 2014).

Figure 1 summarizes the context for IC 443 and 3C 391. The stellar remnant for IC 443 has been found (C. M. Olbert et al. 2001), but none has yet been found for 3C 391. The explosion sites were likely closer to the surfaces of their molecular clouds than the geometric centers of the X-ray emission or the radio shell that include much faster shocks into lower-density gas. The ages and association with molecular clouds are similar for the two SNRs. A significant difference is their distance: ~2 kpc for IC 443, and ~7 kpc for 3C 391 (Y.-H. Lee et al. 2020). Both SNRs are very radio bright (D. A. Green 2019); IC 443 is brighter as seen from Earth, but 3C 391 is actually more radio luminous.

The IC 443 and 3C 391 supernova blast waves are impinging upon dense (10$^{3}$–10$^{5}$ cm$^{-3}$) molecular gas, providing a detailed view of shock fronts where the parent molecular clouds of the progenitors are being disrupted. The warm, shocked gas and compressed magnetic field are also an apparent origin for cosmic rays, which make IC 443 a prominent γ-ray source (A. A. Abdo et al. 2010). The γ-rays from IC 443 have been interpreted as interaction between cosmic rays and a molecular cloud in front of the SNR (D. F. Torres et al. 2008), where fast (~150 km s$^{-1}$) shock fronts act as a “hadronic accelerator” (M. Tavani et al. 2010).

In this paper, we study the emission from molecular hydrogen and water using new observations of the velocity-resolved spectral line profiles of transitions among low energy levels. We concentrate on two H$_{2}$ lines that arise from very different energy levels. The first part of this project (W. T. Reach et al. 2019, hereafter referred to as Paper I) was to measure the 0–0 S(5) line of H$_{2}$, a ground vibrational

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Space Science Institute, 4765 Walnut St, Suite B, Boulder, CO 80301, USA; wreach@spacesciece.org</footnote>

<footnote n="2">²Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121, Bonn, Germany</footnote>

<footnote n="3">³Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="4">⁴Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 75005, Paris, France</footnote>

<footnote n="5">⁵Observatoire de Paris, PSL University, Sorbonne Université, LERMA, 75014 Paris, France</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Diagram showing the context for SNRs IC 443 and 3C 391, produced by supernovae near or just inside the edges of molecular clouds. The blast wave propagates through lower-density gas making X-rays (blue), then into the molecular clumps (purple) and dense cores (red, labeled for those of IC 443).</figcaption></figure>

state transition of ortho-H₂ between the rotational *J* = 7 and *J* = 5 levels. The S(5) line is the brightest line emerging from a wide range of molecular shocks with the pressure of the blast wave of a 10⁴ yr old SNR. (The age and pressure of the blast wave are derived in Paper I.) The new observational data in this paper are velocity-resolved spectra of the S(1) line, which we combine with data from Paper I and archival observations of a ground-state H₂O line to better constrain the nature of the shock fronts that SNRs drive into molecular clouds.

## 2. Observations

### 2.1. SOFIA

Observations were made using the Stratospheric Observatory for Infrared Astronomy (SOFIA), which was an infrared-optimized telescope on a 747 aircraft, operated by the United States and German space agencies (E. T. Young et al. 2012) through 2022. The infrared spectra were taken with the Echelon-Cross-Echelle Spectrograph (EXES; M. J. Richter et al. 2018), with a resolving power of 50,000 (6 km s⁻¹). The spectral configuration was centered on the 587 cm⁻¹ (17.03 µm) pure rotational transition of ortho-molecular-hydrogen, i.e., the 0-0 S(1) line. The 3′′2 slit was centered on a known location of bright H₂ emission for each of four supernova-remnant molecular shocks.

The observations of IC 443 focus on the well-known shocked clumps B, C, and G (Y. L. Huang et al. 1986). Each of these shocked regions is large (up to an arcmin), so the precise observing locations were derived from more recent imaging observations. For IC 443 C, we observed the peak in the H₂ S(5) emission line map from Spitzer (D. A. Neufeld et al. 2007). For the other three sources, we observed the same positions as were previously observed in the S(5) line with SOFIA previously (Paper I). Specifically, we observed a 2MASS Ks peak (also used in Paper I) due to vibrationally

<table><caption>**Table 1**
SOFIA Observing Log<sup>a</sup></caption><thead><tr><th>Date</th><th>Flight</th><th>Altitude<br>(feet)</th><th>Elevation<br>(°)</th><th>Time<br>(minutes)</th></tr></thead><tbody><tr><td colspan="5">IC 443 C (06:17:42.4 +22:21:27)</td></tr><tr><td>2022 Apr 27</td><td>860</td><td>38,020</td><td>37</td><td>22</td></tr><tr><td>2022 Apr 28</td><td>861</td><td>39,000</td><td>30</td><td>12</td></tr><tr><td>2022 Apr 29</td><td>862</td><td>39,010</td><td>25</td><td>10</td></tr><tr><td colspan="5">IC 443 G (06:16:43.6 +22:32:37)</td></tr><tr><td>2022 Apr 27</td><td>860</td><td>39,000</td><td>26</td><td>6</td></tr><tr><td>2022 Apr 28</td><td>861</td><td>39,010</td><td>36</td><td>13</td></tr><tr><td colspan="5">IC 443 B (06:17:16.1 +22:25:15):</td></tr><tr><td>2022 Apr 29</td><td>862</td><td>39,010</td><td>35</td><td>138</td></tr><tr><td colspan="5">3C 391 (18:49:22.6 -00:57:17)</td></tr><tr><td>2022 Apr 27</td><td>860</td><td>43,020</td><td>53</td><td>20</td></tr></tbody></table>

**Note.**

excited H₂ in IC 443 G; the Herschel H₂O and [O I] peak in IC 443 B (from archival observation 1342227773); and a SOFIA H₂ S(5) peak in 3C 391 from Paper I.

Table 1 summarizes the observing conditions, with the aircraft altitude and telescope elevation specified because telluric absorption depends significantly on those two factors. The telescope was nodded between the target position and a nearby position with no trace of H-2 in the 2MASS Ks band or Spitzer line images. But the sky subtraction is never perfect, leading to contamination of the nod-differenced spectra by emission lines from the Earth's atmosphere. We devised a method to remove the problematic residual emission by adding and subtracting small fractions of the sky-radiance spectrum (the sky flux before nod subtraction). In all H₂ S(1) observations there is only ever one order in the echellogram containing the 587 cm⁻¹ H₂ emission line, even considering the various Doppler shifts to the different targets. At most, two atmospheric lines, from N₂O and CO₂ contaminate the celestial H₂ emission. The sky radiance was scaled by a factor (±0.5%) that matches the amplitude of the features observed in the nod difference, then it was subtracted. Figure 2 illustrates the before and after results in the 2D echellogram for IC 443 B.

The S(5) line observations (from Paper I) were not as significantly affected by atmospheric effects. Figure 3 shows the echelle spectral image for two orders that include part of the velocity range where IC 443 B shocked gas is detected. No telluric features are visible on this scale. There is some evidence for a change in the velocity distribution across the 10" slit, with the upper part of the slit somewhat more concentrated to lower velocities, and the lower part of the slit having a broader velocity distribution.

We used 8" portions of the S(1) and S(5) observed slits to extract the spectra used for analysis in this paper. During the observation, the center of the slit for the S(1) observation was adjusted to the peak of S(5) emission, which is why the H₂ S(1) emission is centered in the slit in Figure 2 but the S(5) emission is near the top of the slit in Figure 3.

Figure 4 shows the S(1) and S(5) profiles. The S(5) spectrum is significantly different from S(1), in particular for IC 443 B. To determine whether there are strong localized variations, subslit spectra were extracted from the top and bottom half of

<footnote><sup>a</sup> Target names and J2000 coordinates, with observing conditions for each flight on which they were observed</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Slit spectrum of the H₂ S(1) line from IC 443 B showing subtraction of residual sky emission. The upper panel shows one echelle order from standard nod subtraction. Sky lines from CO₂ and N₂O remain, with the N₂O lines partially overlapping the region of celestial H₂ S(1) emission. The lower panel shows the same order, after subtracting an additional small fraction of the sky spectrum. The H₂ line is labeled below the lower panel. While the H₂ line is visible in both panels, the suppression of sky features significantly improves the accuracy of measuring the line profile. The spectrum used for analysis was extracted from an 8" portion of the slit near its center.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Spectral image for the IC 443 B observations of the H₂ S(5) line. Two orders of the echelle spectrum are shown, because they both cover part of the very broad celestial H₂ line. The minimum and maximum wavenumber and velocity of each order are labeled on their lower and upper corners, respectively. An 8" portion of the slit, toward the top, was extracted to match the region extracted for the S(1) line.</figcaption></figure>

the S(5) slit for IC 443 B. There is in fact a difference, but only in the amount of lower-velocity emission, between −17 and 0 km s⁻¹. Thus it appears that the high-velocity wing (−20 to −60 km s⁻¹) is widespread across the 10" slit, and the lower-velocity emission is more spatially concentrated. The S(1) emission concentrated to lower velocities and the center of the slit.

## 2.2. Archival

Observations from the ESA Herschel Space Observatory (G. L. Pilbratt et al. 2010) taken with the HIFI (T. de Graauw et al. 2010) instrument were obtained from the ESA science archive. The specific observation identifiers were 11342205310 (IC 443 C), 11342205311 (IC 443 B), and 11342205312 (IC 443 G), and 1342218432 (3C 391). Each observation was an on-the-fly spectral map, with the receiver tuned to the H₂O 1₁₀–1₀₁ spectral line at 556.936 GHz in the lower sideband of the dual-sideband system. The spectral resolution was 0.6 km s⁻¹, and the diffraction resolution was 38". The pixel in the HIFI spectral cube containing the SOFIA-observed position was extracted, and the horizontal and vertical polarizations were averaged.

Figure 4 shows the H₂O spectra. The velocity distribution is generally similar to that of the H₂ S(5) line, except for IC 443 G, where there is a deep absorption dip at −4 km s⁻¹ (discussed in Section 4.6 below).

Table 2 summarizes the observed properties of the spectral lines. The line centroid and FWHM were measured using the

first and second velocity moments of the velocity spectra, respectively. The peak intensity is in the observational units of Jy arcsec⁻², where 1 Jy is 10⁻²³ erg s cm⁻² Hz⁻¹.

## 3. Modeling the Supernova Shocks into Molecular Cores

Theoretical models of shocks into dense clumps were calculated with the Paris-Durham shock code (P. Lesaffre et al. 2013; D. R. Flower & G. Pineau des Forêts 2015; B. Godard et al. 2019). The models are similar to those in the extensive grid (L. E. Kristensen et al. 2023) that is available from the interstellar medium database<fnref n="6" /> of the Paris Observatory. In summary, the models begin with interstellar conditions calculated for at equilibrium for a specific ($n_0$), external radiation field ($\chi$), cosmic-ray ionization rate $\zeta_{CR}$, and magnetic field B. The output from the static calculation is used to initialize a photodissociation region, which is slowly advected to the column density $2 \times 10^{20} \text{cm}^{-2}$, to calculate the preshock abundances of all species. The models considered in this paper all have radiation field equal to the average interstellar value ($\chi = 1$) based on J. S. Mathis et al. (1983), with an intervening extinction $A_V = 1$. The magnetic field strength is set for each preshock density according to flux-freezing scaling, with a normalization such that $B = \beta n_0^{-0.5} \mu G$ (with $n_0$ in cm⁻³).

<footnote n="6">⁶ http://ism.obspm.fr/</footnote>
<figure><img src="image_7.png" /></figure>

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

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

<figure><img src="image_10.png" /><figcaption>**Figure 4.** SOFIA spectra of H₂ S(1) and S(5) and H₂O 1₀⁻¹₀₁ lines from IC 443 clumps B, C, G, and from 3C 391 at the positions listed in Table 1. The H₂ lines are in units of Jy arcsec⁻² and the H₂O line is in units of brightness temperature (K). The S(5) and H₂O spectra are vertically offset by 1 and 2 for clarity. All brightnesses for 3C 391 were doubled to be on the same scale.</figcaption></figure>

<table><caption>**Table 2**

Properties of the Observed Spectral Lines</caption><thead><tr><th rowspan="2">Source<br>Model</th><th colspan="3">H<sub>2</sub> S(1)</th><th colspan="3">H<sub>2</sub> S(5)</th><th colspan="3">H<sub>2</sub>O 1<sub>10</sub>-1<sub>01</sub></th></tr><tr><th>⟨v⟩<br>(km s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>I<sub>peak</sub><br>(Jy arcsec<sup>-2</sup>)</th><th>⟨v⟩<br>(km s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>I<sub>peak</sub><br>(Jy arcsec<sup>-2</sup>)</th><th>⟨v⟩<br>(km s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>I<sub>peak</sub><br>(Jy arcsec<sup>-2</sup>)</th></tr></thead><tbody><tr><td>IC 443 B</td><td>-19.6</td><td>28.6</td><td>0.91</td><td>-36.6</td><td>48.0</td><td>0.6</td><td>-34.3</td><td>36.7</td><td>0.06</td></tr><tr><td>IC 443 C</td><td>-20.7</td><td>33.5</td><td>0.60</td><td>-19.9</td><td>37.1</td><td>1.4</td><td>-27.2</td><td>37.0</td><td>0.04</td></tr><tr><td>IC 443 G</td><td>-9.9</td><td>23.2</td><td>1.8</td><td>-11.2</td><td>24.0</td><td>2.0</td><td>-11.4</td><td>25.0</td><td>0.12</td></tr><tr><td>3C 391</td><td>111.4</td><td>24.1</td><td>0.39</td><td>110.0</td><td>21.8</td><td>0.4</td><td>111.2</td><td>28.1</td><td>0.03</td></tr></tbody></table>

The cosmic-ray ionization rate was set to $\zeta_{CR} = 5 \times 10^{-16} \text{ s}^{-1}$, which is 17 times higher than the models presented in the initial study of irradiated shocks (B. Godard et al. 2019), inspired by the high cosmic-ray ionization rate inferred from the H$_3^+$ abundance (N. Indriolo et al. 2010) and the detection of TeV γ-rays toward IC 443 (M. Tavani et al. 2010; V. G. Sinit-syna & V. Y. Sinitsyna 2023). Models with $\zeta$ a factor of 5 higher and lower than this were examined for differences, but they were not significant for H$_2$ lines, which are the primary subject of this paper.

The physics of shocks into dense, cold gas is largely determined by the ratio of the shock speed to the Alfvén speed in the charged fluid taking into account charged dust grains,

$$V_{A,c} = \frac{B}{\sqrt{4\pi\rho_c}} \approx 21\beta \text{ km s}^{-1}, \quad (1)$$

for the estimated abundance and charge of small grains in molecular clouds (as in V. Guillet et al. 2007). Shocks with speed $V_s > V_{A,c}$ have an abrupt “jump” in physical conditions
at the shock front, followed by gradual cooling; these are commonly referred to as J-type shocks. Molecules can be destroyed in fast shocks, and by the time they reform, the gas velocities are low. The slower, C-type shocks occur at $V_s < V_{A,c}$ because the neutral species decouple from the ions, which remain strongly coupled to the magnetic field (B. T. Draine & C. F. McKee 1993). This allows neutrals to flow past the initial shock front and avoid the abrupt heating that ions experience, allowing some molecules (notably H₂) to avoid dissociation.

The velocities at which molecular emission arises depend on the shock type and speed. In this discussion and throughout this paper, velocities will be measured in the rest frame of the observer ("local standard of rest"; J. Binney & M. Merrifield 1998). This is important because it is shared by all gas struck by the blast wave: we assume that all gas starts at rest and is impacted by the blast wave at the same time. Other frames, such as the commonly used one comoving with the shock front, depend upon the shock velocity, which is not known a priori and differs for gas of different density, making it difficult to compare to observational data. Let us call the preshock velocity $V_0$ (with respect to the observer), with temperature $T_0$, and the shock velocity $V_s$ (with respect to the gas being shocked). The peak temperature $T_s = 53V_s^2$ K for a J-type shock, with $V_s$ in $km s^{-1}$ (P. Lesaffre et al. 2013). For a J-type shock, a very thin layer of hot molecules ($T \sim T_s$) can exist at the shock front, but the bulk of the molecules only return after molecules reform in the cooled ($T \sim T_0$) and decelerated gas; see A. Lehmann et al. (2020) for a modern description of J-type shocks in dense gas. The observed emission from a dissociative J-type shock would comprise highly excited molecules at velocity $V_s + V_0$ and cold (reformed) molecules at $\sim V_0$, with narrow (thermal) velocity distributions. Emission from the colder, slower layer would be from lower energy levels. For a C-type shock, the molecules persist throughout the postshock gas, and molecular emission can arise from velocities ranging from $V_s + V_0$ to $V_0$. The shape of the velocity distribution depends upon the cooling and molecule reformation rate in the postshock gas. For very young shocks, a CJ-type shock regime exists (P. Lesaffre et al. 2004), for which the region very close to the shock front is treated as a J-type shock and the larger-scale evolution as a C-shock. Molecular emission from a broad distribution of velocities is a unique signature of C-type shocks.

The magnetic field strength, parameterized by $\beta$, is crucial for determining the shock structure and the distribution of molecules behind the shock front. Observations indicate $\beta \sim 1$ in interstellar clouds; however, there is significant uncertainty and true spread in $\beta$, because only the radial component is measured, and the conversion to the full magnetic field strength is challenging (R. M. Crutcher 2012). For $\beta = 1$, shocks faster than $21 \text{ km s}^{-1}$ are J type. Using the ram pressure of a supernova blast wave at $10^4$ yr after explosion, this corresponds to shocks into gas with density $n_0 < 10^4 \text{ cm}^{-3}$. Zeeman observations of OH masers from supernova shocks indicate very high field strengths $\sim 1$ mG (C. L. Brogan et al. 2000) in the compressed, postshock gas. For a shock into gas with preshock density $10^4 \text{ cm}^{-3}$ and compression factor 10 (which matches the density required in the masing region), we infer $\beta \sim 3$. We ran elevated magnetic field models with a modest increase to $\beta = 5$. The supernova blast-wave ram pressure is expected to drive $80 \text{ km s}^{-1}$ shocks into $n_0 = 10^3 \text{ cm}^{-3}$ gas, which would be C type if $\beta > 4$, so in order to model C-type

shocks at lower densities, we modeled shocks into $n_0 \le 10^4 \text{ cm}^{-3}$ with both $\beta = 1$ and $\beta = 5$. Another consideration is that only the component of the magnetic field vector parallel to the shock affects the 1D models, so the effective $\beta$ is scaled by the cosine of the angle between the magnetic field and shock velocity vectors. The shocks into high-density ($n > 10^4$) gas with $\beta = 1$ are C type, so to model the case with a decreased or misaligned magnetic field, we generated models for $\beta = 0.1$ and $\beta = 1$.

The phase space of critical model parameters—preshock density ($n_0$), shock velocity ($V_s$), and magnetic field ($\beta$), and the preshock density, shock velocity, and magnetic $\beta$ for the set of fiducial models—is listed in Table 3. This phase space is for a fixed ram pressure equal to that of $10^4$ yr old SNRs appropriate for IC 443 and 3C 391, so the combination of $n_0$ and $V_s$ is tied together (see Paper I). The boundary between C-type and J-type shocks (discussed in the previous paragraph) is indicated in Figure 5. Table 3 lists the suite of models in order of increasing $n_0$ from (a) through (g). Models with a nominal magnetic field ($\beta = 1$) have no subscript, and the models with elevated and decreased $\beta$ have subscripts p or m, respectively. The locations of all these models are shown in Figure 5. An important distinction among the J-type shocks is whether H₂ molecules are destroyed. H₂ molecules can survive shocks slower than $25 \text{ km s}^{-1}$ (A. G. G. M. Tielens 2005). For our model phase space, the high density ($n_0 > 10^4 \text{ cm}^{-3}$), low $\beta$ (oblique or reduced magnetic field) are nondissociative, J-type shocks.

The distinctions between the different shock types are evident in the H₂ distribution and emissivity. Figure 6 shows the emissivity of the H₂ S(1) and S(5) lines for the models listed in Table 3, and Figure 7 shows the predicted spectra (discussed in more detail below). The atomic H density, which is effectively emissivity of the 21 cm atomic H hyperfine line, is included in Figure 6 to show where H₂ molecules are dissociated. The distinction between the J- and C-type shocks is clear when comparing the left- and right-hand sides of the figure. The dissociative J-type shocks (models a, b, c, and d) have two regions of H₂ emission, with atomic H filling the gap in between. A hot emitting region close ($3 \times 10^{-7}$ pc) behind the shock front is where the molecules are being destroyed; it is spatially so thin that it does not contribute more than 1% of the integrated emission for the models considered here. (But some species, like CH+, only exist in this hot region because they are quickly chemically converted in the colder region where molecules reform.) The nondissociative J-type shocks (models $e_m$, $f_m$, and $g_m$) look more like C-type shocks, because molecules survive. However, the lack of strong magnetic field means neutrals and ions stay together, which affects the shock structure. Comparing to the stronger magnetic field C-type shocks into the same densities (models e versus $e_m$, f versus $f_m$, and g versus $g_m$), the J-type shocks are much more compact. The H₂ emitting regions of C-type shocks are of order $10^{-3}$ pc wide, while the nondissociative J-type equivalents are only of order $10^{-6}$ pc wide.

Table 3 summarizes the distance of the emitting regions behind the shock fronts, and their temperatures. The average emitting distance was calculated as follows:

$$ \ell_{\text{emiss}} \equiv \frac{\int \ell \mathcal{E}(\ell) d\ell}{\int \mathcal{E}(\ell) d\ell}, \qquad (2) $$
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Phase space of shock parameters. The three critical regimes (nondissociative C-type, dissociative J-type, and nondissociative J-type) are labeled (and shaded purple, green, and cyan, respectively). The input parameters for the shocks in Table 3 are indicated by their corresponding letter designations. Model (e) falls near the C- to J-type boundary and was found to be C type. This phase space is in accord with similar ones from (B. Godard et al. 2019; L. E. Kristensen et al. 2023) diagrams.</figcaption></figure>

<table><caption>**Table 3**
Distance and Temperature of Emitting Region<sup>a</sup></caption><thead><tr><th rowspan="2"></th><th colspan="3">Model</th><th colspan="2">H<sub>2</sub> S(1)</th><th colspan="2">H<sub>2</sub> S(5)</th><th colspan="2">H<sub>2</sub>O 1<sub>10</sub>-1<sub>01</sub></th></tr><tr><th>density (cm<sup>-3</sup>)</th><th>shock speed (km s<sup>-1</sup>)</th><th>β</th><th>⟨T⟩ (K)</th><th>ℓ<sub>emiss</sub> (mpc)</th><th>⟨T⟩ (K)</th><th>ℓ<sub>emiss</sub> (mpc)</th><th>⟨T⟩ (K)</th><th>ℓ<sub>emiss</sub> (mpc)</th></tr></thead><tbody><tr><td colspan="10" style="text-align: center; font-weight: bold;">C-type</td></tr><tr><td>(a<sub>p</sub>)</td><td>10<sup>3</sup></td><td>80</td><td>5</td><td>1400</td><td>24</td><td>4000</td><td>9</td><td>1100</td><td>26</td></tr><tr><td>(b<sub>p</sub>)</td><td>10<sup>3.3</sup></td><td>60</td><td>5</td><td>1300</td><td>23</td><td>2600</td><td>12</td><td>1100</td><td>26</td></tr><tr><td>(c<sub>p</sub>)</td><td>10<sup>3.8</sup></td><td>37</td><td>5</td><td>920</td><td>19</td><td>1300</td><td>13</td><td>820</td><td>31</td></tr><tr><td>(d<sub>p</sub>)</td><td>10<sup>4</sup></td><td>30</td><td>5</td><td>770</td><td>19</td><td>940</td><td>17</td><td>700</td><td>29</td></tr><tr><td>(e)</td><td>10<sup>4.3</sup></td><td>22</td><td>1</td><td>870</td><td>3.4</td><td>1490</td><td>2.3</td><td>410</td><td>5.5</td></tr><tr><td>(f)</td><td>10<sup>4.8</sup></td><td>15</td><td>1</td><td>700</td><td>1.1</td><td>1020</td><td>0.8</td><td>500</td><td>1.4</td></tr><tr><td>(g)</td><td>10<sup>5</sup></td><td>10</td><td>1</td><td>520</td><td>0.5</td><td>350</td><td>0.9</td><td>330</td><td>0.7</td></tr><tr><td colspan="10" style="text-align: center; font-weight: bold;">J-type</td></tr><tr><td>(a)</td><td>10<sup>3</sup></td><td>80</td><td>1</td><td>420</td><td>3.4</td><td>860</td><td>1.5</td><td>410</td><td>11</td></tr><tr><td>(b)</td><td>10<sup>3.3</sup></td><td>60</td><td>1</td><td>150</td><td>2.0</td><td>210</td><td>2.2</td><td>94</td><td>11</td></tr><tr><td>(c)</td><td>10<sup>3.8</sup></td><td>37</td><td>1</td><td>190</td><td>1.5</td><td>290</td><td>1.7</td><td>97</td><td>11</td></tr><tr><td>(d)</td><td>10<sup>4</sup></td><td>30</td><td>1</td><td>150</td><td>1.3</td><td>200</td><td>1.7</td><td>120</td><td>7</td></tr><tr><td>(e<sub>m</sub>)</td><td>10<sup>4.3</sup></td><td>22</td><td>0.1</td><td>480</td><td>0.006</td><td>1940</td><td>0.008</td><td>220</td><td>0.05</td></tr><tr><td>(f<sub>m</sub>)</td><td>10<sup>4.8</sup></td><td>15</td><td>0.1</td><td>1100</td><td>0.001</td><td>1040</td><td>0.002</td><td>290</td><td>0.03</td></tr><tr><td>(g<sub>m</sub>)</td><td>10<sup>5</sup></td><td>10</td><td>0.1</td><td>1720</td><td>0.001</td><td>1710</td><td>0.002</td><td>490</td><td>0.04</td></tr></tbody></table>

**Note.**

where $E(ℓ)$ is the emissivity (erg cm$^{-3}$) of a spectral line versus distance $ℓ$ behind the shock. The C-type shocks become progressively more compact for regions of higher density, and the H$_2$ emissivity for the S(5) line is closer to the shock than that of the S(1) line. This is intuitively obvious because the

H$_2$ molecules survive the entire shock, emitting as they cool, and the S(5) line is from a higher energy level tracing more highly excited gas. Note that if one were to use the ratio of observed brightness of two spectral lines to determine the gas temperature, the result could be misleading because the

<footnote><sup>a</sup> For each spectral line $ℓ_{emiss}$ and $T_{emiss}$ are the emissivity-weighted distance and temperature of the emitting region behind the shock front.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Modeled emissivity of the H₂ S(1) and S(5) lines as a function of distance behind the model shocks listed in Table 3. The emissivities of the H₂ lines are normalized to the peak value for each modeled transition. To show where molecules are dissociated, the fractional abundance of atomic H is shown in each panel, without normalization (so a value of 1 means complete dissociation of H₂). Each panel is labeled with its shock designation from Table 3, abbreviated preshock density (cm⁻³), shock velocity (km s⁻¹), and magnetic parameter β. The left column shows the C-type shock models, which become progressively more compact as the preshock density increases. The right column shows the J-type shock models. The dissociative J-type models (a, b, c, d) have a thin region just behind the shock front (2 × 10⁻⁷ pc), followed by a region where molecules are destroyed, then followed by a region of reformed molecules (around 10⁻³ pc). Note the two columns have different horizontal scales. The vertical dashed line shows the shock duration of 5000 yr, and vertical dotted lines show how far the shock developed in (from right to left) 3000, 1000, and 100 yr.</figcaption></figure>

emission in the two lines do not arise from the same location. The H₂ S(1) line arises from gas cooler than where the S(5) line arises, and the H₂O 1₀⁻¹₀₁ line arises from gas that is cooler still, in keeping with their respective upper-level energies above ground state.

When predicting the emission from a short-lived phenomenon like an SNR, the duration must be considered. The SNRs are only ~10⁴ yr old, which is not much more than the cooling time of the gas. Using Figure 30.1 of B. T. Draine (2011), the cooling time for gas with density 10³ cm⁻³ is of order 10³ yr. If

the shock front took a fraction *f* of the total time 10⁴ yr since the explosion to reach the shocked clump, then the shock duration is (1-*f*) × 10⁴ yr. In both IC 443 and 3C 391, the shocked molecular clumps are close to the edges of their SNR, and they span a range of projected distances from any explosion site, suggesting 1 > *f* ≳ 0.5, which means shock duration from 0 to 5000 yr. Inspecting the temperature profiles of the shocked gas, the gas cools within 3000 yr for all models. For the shocks into gas with *n*₀ > 10⁴ cm⁻³, the shocks are essentially stationary because cooling times are short. Vertical lines in Figure 6
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Modeled spectra of H₂ emission in the S(1) and S(5) pure rotational lines, and the H₂O 1<sub>10</sub>-1<sub>01</sub> line, for C-type MHD shocks. Each panel is labeled with its shock designation from Table 3, abbreviated preshock density (cm<sup>-3</sup>), shock velocity (km s<sup>-1</sup>), and magnetic parameter β. All spectra are normalized by their peak values, which are reported in Table 4. Velocities are in the rest frame of the observer. The dashed line at 0 km s<sup>-1</sup> indicates the preshocked gas velocity, and the modeled shock is propagating toward the observer (along the line of sight). Unlike Figure 6, the J-type models are not shown, because they are all simple Gaussians centered near zero.</figcaption></figure>

illustrate how far the shock front develops in a clump if it arrived 100, 1000, 3000, and 5000 yr ago. Figure 8 shows the emitted spectra for different ages. The gas in very young shocks has not decelerated, so that the velocity distributions skew further from zero. For the youngest shock (100 yr), the emission arises from just behind the shock front, at velocity $V_s$. None of the observations discussed in this paper are likely to be such young shocks, because they all have emission spanning $V_0$ to $V_s$. But the accelerated clump IC 443:B2 (Paper 1) may be such a case, where the blast wave has only recently reached the clump, perhaps because it is farther from the explosion site. In the remainder of this paper, we use the models for a shock duration of 3000 yr, which is plausible in terms of the SNR age and widespread locations of the shocked clumps.

Synthetic emission spectra were calculated for each model using the emissivity versus distance from Figure 6 and a velocity distribution at each distance that is centered on the neutral gas velocity and has a local width that is a combination of sound speed at that distance and a fixed turbulent dispersion (3 km s⁻¹). Figure 7 shows the predicted spectra from the three different shock models in each of the H₂ S(1) and S(5) lines, for the case where the shock propagates along the line of sight and toward the observer. The rest velocity of the preshock gas is at 0 km s⁻¹. If the shock is propagating perpendicular to the line of sight, the emission lines all center on 0 km s⁻¹, and there is no further diagnostic information from the spectral shape. For intermediate inclination angles, the centroid velocity of the shocked emission shifts approximately as the cosine of the inclination angle (if inclination 0 is defined as propagating parallel to the line of sight).

Synthetic H₂O 1<sub>10</sub>-1<sub>01</sub> spectra were calculated using the velocity, density, temperature, and H₂O abundance at each distance behind the shocks. The excitation of the upper level was calculated with the collision rates from S. Green et al. (1993) and radiative rates from H. M. Pickett et al. (2003). Collisional rates were multiplied by 1.35 to include collisions with He. At each distance behind the shock, the collisional rates were evaluated by interpolating between tabulated temperatures for the downward rates. The upward transition rates, γ<sub>lu</sub>, were precisely calculated from the downward rates, γ<sub>lu</sub>, based on detailed collisional balance. This step prevents numerical error that arise from using tabulated upward and downward rates. An important caveat to the predicted H₂O model line profiles is that we assumed the emission is optically thin. A treatment of optically thick emission is beyond our present scope, though we note it has been attempted for CO (N. Fuda et al. 2023). The resulting synthetic H₂O 1<sub>10</sub>-1<sub>01</sub> spectra are very similar to those of the H₂ S(1) line, so we did not include them in the figures. The H₂O emitting-region thickness, line flux, and peak intensity for each model is included in Tables 3 and 4.

# 4. Comparing Observations and Models

## 4.1. *Preshock Conditions*

First we must establish the rest velocity of the preshock gas and the direction of the shock along the line of sight. The ambient, preshock molecular gas into which the IC 443 blast wave is propagating has been mapped in low-J CO emission, and it arises mostly at radial velocities between -6 and
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Same as Figure 7, but for only the three fastest C-type shocks $a_p$, $b_p$, and $c_p$, at three different ages of 3000 yr (left), 1000 yr (center), and 100 yr (right). The velocity distribution skews away from zero if the shocks are younger, because the accelerated gas closer to the shock becomes a larger fraction of the emission.</figcaption></figure>

<table><caption>**Table 4**
Velocity Distributions from Shock Models<sup>a</sup></caption><thead><tr><th rowspan="2">Model</th><th colspan="4">H<sub>2</sub> S(1)</th><th colspan="4">H<sub>2</sub> S(5)</th><th colspan="4">H<sub>2</sub>O 1<sub>10</sub>-1<sub>01</sub></th></tr><tr><th>⟨v⟩<sup>b</sup></th><th>FWHM<sup>c</sup></th><th>I<sub>peak</sub><sup>d</sup></th><th>I<sub>tot</sub><sup>e</sup></th><th>⟨v⟩<sup>b</sup></th><th>FWHM<sup>c</sup></th><th>I<sub>peak</sub><sup>d</sup></th><th>I<sub>tot</sub><sup>e</sup></th><th>⟨v⟩<sup>b</sup></th><th>FWHM<sup>c</sup></th><th>I<sub>peak</sub><sup>d</sup></th><th>I<sub>tot</sub><sup>e</sup></th></tr></thead><tbody><tr><td colspan="13" style="text-align:center; font-weight:bold;">C-type</td></tr><tr><td>(a<sub>p</sub>)</td><td>-21.5</td><td>34</td><td>0.43</td><td>0.23</td><td>-49.5</td><td>40</td><td>1.24</td><td>3.9</td><td>-17.6</td><td>26</td><td>2.34</td><td>0.034</td></tr><tr><td>(b<sub>p</sub>)</td><td>-22.9</td><td>33</td><td>0.64</td><td>0.34</td><td>-41.6</td><td>27</td><td>2.63</td><td>5</td><td>-18.3</td><td>26</td><td>5.05</td><td>0.075</td></tr><tr><td>(c<sub>p</sub>)</td><td>-21.2</td><td>21</td><td>1.20</td><td>0.65</td><td>-27.3</td><td>13</td><td>2.18</td><td>1.9</td><td>-17.7</td><td>19</td><td>15.58</td><td>0.23</td></tr><tr><td>(d<sub>p</sub>)</td><td>-19.8</td><td>17</td><td>1.52</td><td>0.78</td><td>-22.8</td><td>11</td><td>1.01</td><td>0.73</td><td>-16.6</td><td>15</td><td>26.39</td><td>0.38</td></tr><tr><td>(e)</td><td>-10.2</td><td>14</td><td>0.71</td><td>0.26</td><td>-15.4</td><td>11</td><td>5.50</td><td>4</td><td>-6.3</td><td>13</td><td>71.48</td><td>0.68</td></tr><tr><td>(f)</td><td>-8.2</td><td>11</td><td>1.00</td><td>0.31</td><td>-10.7</td><td>8</td><td>3.42</td><td>1.9</td><td>-5.0</td><td>10</td><td>215.62</td><td>1.8</td></tr><tr><td>(g)</td><td>-5.9</td><td>8</td><td>0.84</td><td>0.19</td><td>-7.1</td><td>8</td><td>0.20</td><td>0.11</td><td>-3.9</td><td>9</td><td>144.20</td><td>1.1</td></tr><tr><td colspan="13" style="text-align:center; font-weight:bold;">J-type</td></tr><tr><td>(a)</td><td>-1.3</td><td>8</td><td>0.01</td><td>0.0012</td><td>-1.3</td><td>8</td><td>0.06</td><td>0.033</td><td>-1.3</td><td>8</td><td>0.01</td><td>6.4e-05</td></tr><tr><td>(b)</td><td>-1.3</td><td>8</td><td>0.01</td><td>0.0015</td><td>-1.3</td><td>8</td><td>0.11</td><td>0.057</td><td>-1.3</td><td>8</td><td>0.00</td><td>1.5e-05</td></tr><tr><td>(c)</td><td>-1.3</td><td>7</td><td>0.01</td><td>0.0022</td><td>-1.4</td><td>8</td><td>0.24</td><td>0.12</td><td>-1.3</td><td>8</td><td>0.01</td><td>6.9e-05</td></tr><tr><td>(d)</td><td>-1.3</td><td>7</td><td>0.02</td><td>0.0033</td><td>-1.3</td><td>8</td><td>0.32</td><td>0.16</td><td>-1.3</td><td>8</td><td>0.02</td><td>0.00015</td></tr><tr><td>(e<sub>m</sub>)</td><td>-0.2</td><td>7</td><td>0.04</td><td>0.0081</td><td>-0.5</td><td>7</td><td>0.37</td><td>0.17</td><td>-0.2</td><td>7</td><td>359.92</td><td>2.2</td></tr><tr><td>(f<sub>m</sub>)</td><td>-0.3</td><td>7</td><td>0.06</td><td>0.011</td><td>-0.3</td><td>7</td><td>0.66</td><td>0.3</td><td>-0.2</td><td>7</td><td>256.97</td><td>1.5</td></tr><tr><td>(g<sub>m</sub>)</td><td>-0.7</td><td>7</td><td>0.02</td><td>0.0042</td><td>-0.7</td><td>7</td><td>0.76</td><td>0.35</td><td>-0.4</td><td>7</td><td>95.10</td><td>0.57</td></tr></tbody></table>

**Notes.**

<sup>a</sup> For shocks directed along the line of sight. Models are labeled as in Table 3 and Figures 5, 6, and 7.

<sup>b</sup> Mean velocity of the emission line profile, in km s<sup>-1</sup>.

<sup>c</sup> FWHM intensity of the emission profile, in km s<sup>-1</sup>.

<sup>d</sup> Peak of the spectral intensity in observational units of Jy arcsec<sup>-2</sup>.

<sup>e</sup> Integrated line surface brightness, in units of 10<sup>-3</sup> erg s<sup>-1</sup> cm<sup>-2</sup> sr<sup>-1</sup>.

−1 km s<sup>−1</sup> with respect to the Galactic neighborhood local standard of rest (J.-J. Lee et al. 2012). For 3C 391, the parent molecular cloud is at 90–110 km s<sup>−1</sup> from the CO (1–0) mapping (D. J. Wilner et al. 1998), and the ambient gas absorbs the shocked emission in the CO(2–1) spectra and produces narrow CS(2–1) emission at 105 km s<sup>−1</sup> (W. T. Reach & J. Rho 1999), which we adopt as the rest velocity. The velocity range of the shocked emission relative to the ambient gas, and the relative velocity spreads of the S(1) and S(5) lines, allow us to determine the orientation of the shocks, which is required before comparing to the models. If the shock were propagating perpendicular to the line of

sight, the velocity gradient behind the shock would be masked and all lines would be centered on the preshock velocity. This appears to be approximately the case for IC 443 G and 3C 391, which agrees with the presence of significant self-absorption in the molecular lines from the lowest energy levels or ground state.

If the shock is propagating along the line of sight, then the shocked gas will be displaced in velocity as predicted by the synthetic spectra in Figure 7 with centroids in Table 4. The H<sub>2</sub> emission from IC 443 B and IC 443 C is significantly offset from the preshock velocity, indicating those shocks are along the line of sight. The shocked gas (Figure 4) is at more negative
velocities than the ambient gas for both SNRs, meaning the shocks are propagating toward us along the line of sight.

## 4.2. Velocity Range

For IC 443 B, the H₂ and H₂O emission spans the entire range from 0 to −80 km s⁻¹. As explained above, we can infer from this that $V_0 \sim 0$ km s⁻¹, $V_s \sim 80$ km s⁻¹, and the shock is moving toward the observer. Using only this shock speed, and the IC 443 blast-wave pressure, the preshock density is of order $10^3$ cm⁻³. Therefore, the initial guess is that models (a) and (aₚ) may apply to IC 443 B. In Paper I, we fit models to the S(5) line profiles and found the best-fitting preshock density $2 \times 10^3$ cm⁻³ and shock speed 60 km s⁻¹, which is similar. The C-type shock with elevated magnetic field, model (aₚ), provides a reasonably good description of the observed range of velocities. On the other hand, the corresponding J-type model, (a), is centered near the velocity of the preshock gas and has a width of only 8 km s⁻¹, which is not at all what was observed. No combination of fully dissociating J-type models can explain the observed velocity distribution. Model (aₚ) also correctly predicts that the H₂ S(1) line will be concentrated at velocities lower than the S(5) line.

For IC 443 C, the observed S(1) and S(5) lines are relatively similar in shape, with a two-component structure. There is a bright component at −8 km s⁻¹ with FWHM 15 km s⁻¹, and a wider component centered at −26 km s⁻¹ with FWHM 28 km s⁻¹. The wide velocity spread indicates the shocks are propagating at least partially along the line of sight (again, toward the observer). The observed FWHM is smaller than for IC 443 B, so a lower-speed shock into higher-density gas applies to IC 443 C. In Paper I, the S(5) line profile was well fitted by a combination of a 60 km s⁻¹ shock into density $2 \times 10^3$ cm⁻³ and a 31 km s⁻¹ shock into a higher density of $10^4$ cm⁻³, similar to models (bₚ) and (cₚ). The S(1) line profile shape is more similar to the higher-density shock model. Referring to Figure 7(b), the models predict S(1) and S(5) lines with similar widths and centers, as observed.

For IC 443 G, the observed S(1) and S(5) lines are again relatively similar and have a two-component structure, but the narrow component is much brighter than the wider component. The narrow components of the S(1) and S(5) lines have FWHM of 10 and 14 km s⁻¹, which are similar to the narrow components of IC 443 C. The deep H₂O absorption provides a clue that the preshock velocity is −4 km s⁻¹, and the observed emission centroids of the S(1) and S(5) lines are −6 and −8 km s⁻¹, respectively. The small shift between the shocked emission and the preshocked gas indicates that the shock is slightly inclined from the line of sight. In Paper I, a viewing angle of 65° was inferred. The small shift in velocity between the S(1) and S(5) lines is in the sense predicted by the shock models for this viewing angle.

## 4.3. Line Brightness and S(1)/S(5) Ratio

The absolute brightness of the observed shocks can be a diagnostic of the shock, as the H₂ lines are a primary coolant (D. R. Flower & G. Pineau des Forêts 2015). The brightness scales approximately as the secant of the inclination with respect to the line of sight, limited by the size of the shocked clump relative to the emitting-region thickness. The increased brightness for shocks perpendicular to the line of sight is partially canceled out by the finite angular resolution of the

observations. From Table 3, the emitting region behind a shock into $10^3$ cm⁻³ gas, at the distance of IC 443, would subtend 2′′1 in the S(1) line and 3′′8 in the S(5) line, which is comparable to the SOFIA angular resolution of the observations (3′′7 at the wavelengths of both the S(1) and S(5) lines; P. Temi et al. 2018); thus the beam filling factor should be of order unity for shocks into $10^3$ cm⁻³ gas, whatever the inclination relative to the line of sight. Shocks perpendicular to the line of sight are also limited by the preshock clump size and by the filling factor of the shock in the beam. The size of dense cores and filaments in molecular clouds is of order $\ell_{\text{clump}} = 0.1$ pc (P. André et al. 2014). The angular resolution of the SOFIA observations is 3′′7 at the wavelengths of both the S(1) and S(5) lines (P. Temi et al. 2018), which is $\ell_{\text{beam}} = 0.036$ pc at the distance of IC 443. The line brightness can be increased for shocks inclined with respect to the line of sight by as much as $\ell_{\text{clump}}/\ell_{\text{beam}} = 2.8$. The observed brightnesses are all lower than the predicted ones, by factors of ~10, taking all these factors into account. That the brightnesses are this close to the predictions is a testament to their validity, as the modeled physics includes MHD, thermodynamics, and chemistry of a complicated and likely nonhomogeneous medium.

The ratio of brightness of the S(1) to S(5) lines is a diagnostic of the shocks independent of both the beam filling factor and the inclination angle. For IC 443 B, the ratio of S(1)/S(5) is 1.5 at the peak of S(1), but it decreases to <0.2 at velocities <−50 km s⁻¹, where S(1) is not detected. For IC 443 C, the ratio S(1)/S(5) is 0.3 in the narrow component (and 0.5 in the tentative wide component). For IC 443 G, the ratio S(1)/S(5) is 0.9. These are in the observational units of Jy arcsec⁻² for peaks of the line profiles. For the models in Table 4, the S(1)/S(5) ratio in the same units is 0.3 for model aₚ and 0.7 for model bₚ, and those models are the best single fits to the line ratios for the broad components of the observed spectral lines. The narrow components have significantly higher S(1), and they arise from cooler gas such as from shocks into higher-density gas or gas that has a long time to cool.

A broader view of the H₂ emission is needed, so we include Spitzer (J. R. Houck et al. 2004; M. W. Werner et al. 2004) observations of H₂ lines S(0) through S(7) from prior work on these targets (D. A. Neufeld et al. 2007; Y. Yuan & D. A. Neufeld 2011; P. Dell'Ova et al. 2024). The observed flux is directly proportional to the column density in the upper level. Figure 9 shows the column densities versus energy of the upper level, where column densities were divided by their statistical weight assuming an ortho-to-para H₂ abundance ratio of 3.

Model (aₚ), which is required to explain the width of the S(5) line, can predict the column densities of the upper levels of the S(4) through S(7) lines. But that model, when normalized to match the observed brightness of the S(5) line, underpredicts the observed flux of the S(0) line by a factor of 8. Shocks into denser gas are much more prolific in the lower-J lines, so we added a second model (g): a 10 km s⁻¹ shock into gas with density $10^5$ cm⁻³. Model (g) was normalized to match 90% of the observed flux in the S(0) line. The combination of these two shock models is able to match the H₂ excitation diagram over the full set of S(0) through S(7) lines. The relative contributions model (aₚ):model (g) for the SOFIA-observed S(1) lines is predicted to be 80%:20% for the S(1) line and 8%:92% for the S(5) line. This explains why there is a broad component, fainter
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 9.** Column density divided by statistical weight for the upper levels of the pure rotational S(0) through S(6) lines toward IC 443 C. The horizontal axis is the energy (expressed as a temperature) of the upper level. The observed values are in black. The predictions from shocks into moderate-density gas (combination of models ($b_p$) and ($c_p$) from Table 3), normalized to match the observed brightness of the S(5) line, are in blue. The predictions from a shock into higher-density gas (model (g) from Table 4), normalized to match the observed brightness of the S(0) line, is in orange. The left panel is for a shock into gas with statistical equilibrium populations of H₂ (OPR = 3). The right panel is for a shock into gas that has remained cold and undisturbed (OPR = 0.006).</figcaption></figure>

in S(1) than in S(5), and a narrow component, relatively brighter in S(1) in the line profiles.

The H₂ ortho-to-para ratio (OPR) for the dense, preshock gas can influence the OPR in the postshock gas, potentially affecting our conclusion that a shock into dense gas is required. For shocks faster than about 20 km s⁻¹, the OPR is completely reset to 3 in the shock front, so the initial OPR only slightly changes the shock structure. This is because the peak temperature behind a faster shock is high enough, for a long enough time, that the OPR of the shocked gas is reset to the high-temperature equilibrium value, and the gas retains that value for longer than the 10⁴ yr age of the SNR. D. A. Neufeld et al. (2024) showed that the OPR for the S(4) through S(7) lines is near the high-temperature limit, in agreement with this prediction of the models regardless of their initial OPR. But for a slow, 10 km s⁻¹ shock, the OPR can remain far from equilibrium, if the preshock gas was out of equilibrium. Therefore, we compared model (g) (which had OPR = 3) with alternatives: a very low initial OPR = 0.006 (which arises from 10⁴ cm⁻³ gas cooling for 60 Myr) and a moderate OPR = 1. For the very low initial OPR model, the shock increases the OPR to only 0.15, and the emitted spectrum from the shocked gas has the distinctive, zig-zag signature of odd-versus-even levels shown in the right-hand panel of Figure 9. For the moderate initial OPR = 1 model, the shock increases the OPR slightly, to 1.2, and the gas remains near OPR = 1 throughout the age of the SNR. The observations show only a slight zig-zag excitation diagram from S(0) to S(3), indicating the preshock gas had 1 < OPR < 3. The observed brightness of the ortho-H₂ S(1) and S(3) lines, compared to the para-H₂ S(0) and S(2) lines, constrains the OPR of the shocked gas to be in the range 2 ± 0.5. This result is similar to what was found from fitting other SNR H₂ excitation diagrams (J. W. Hewitt et al. 2009). In conclusion, the amount of dense, shocked gas can be confidently inferred from the model fits, even using the ortho-H₂ S(1) line, because the OPR is of order unity.

### 4.4. Model Fit to Velocity Distribution

The velocity distribution is the best indicator of the shock speed and postshock gas kinematics. To minimize projection

effects, and to maximize radial velocity spread, we consider in some detail the shocked clump with the broadest and most asymmetric velocity distribution (which is most likely to be a shock along the line of sight): IC 443 B. Figure 10 shows how a linear combination of models from Table 3 compares to the observed spectra. Because we only used a sparse model grid, the combined models have “bumps” corresponding to the peak velocities of each model used in the input. Overall, the models can explain the velocity distribution well in both the S(1) and S(5) lines. The model combination was built as follows. First, the lowest-density (but enhanced magnetic) model, ($a_p$), was required in order to explain the highest-velocity observed gas. Then, a high-density model, (g), was required to explain the brightness of the S(1) line. Two intermediate models were then added to span the broad, continuous velocity distribution. The linear combination that matches the observed spectra had filling factors for models ($a_p$), ($b_p$), (d), and (g) of 0.1, 0.1, 0.3, and 0.7, respectively. Thus the dense clump was a significant presence in the beam, and the lower-density gas fills a small fraction of the beam.

### 4.5. Dissociation of Molecular Hydrogen behind the Shock

The models that match the H₂ line profiles and excitation diagrams are C-type shocks, where the H₂ survives the initial heating and emits as it cools. But such shocks only occur when the magnetic field is both strong and oriented perpendicular to the shock velocity. For the same physical conditions, but a magnetic field oriented more parallel to the shock velocity, dissociative J-type shocks are expected in the low-density gas, as delineated in Figure 5. We now test whether there is significant dissociation of H₂ by measuring the amount and velocity distribution of atomic H behind the shock front. Accelerated atomic gas has been noted before toward IC 443 (R. Braun & R. G. Strom 1986; J.-J. Lee et al. 2008). For regions like the northeastern shell, this shocked atomic gas is expected and associated with shocks into lower-density ($n < 10 \text{ cm}^{-3}$) gas. The atomic gas in those regions arises from the recombined and cooling gas behind the very fast ($V_s > 100 \text{ km s}^{-1}$) J-type shocks, which are prominent at optical wavelengths in H recombination and heavier-element
IC443B

<figure><img src="image_17.png" /><figcaption>**Figure 10.** Observed H<sub>2</sub> S(5) and S(1) spectra of IC 443 B (top) and IC 443 C (bottom) with scaled model spectra. The models are labeled with their name from Table 3. The sum of models (black curve) matches the overall velocity distribution. The bulk of the modeled S(5) flux is from shocks into gas with moderate densities, while the S(1) line is mostly from shocks into high-density gas, although all models contribute to both lines.</figcaption></figure>

IC443C

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

<figure><img src="image_19.png" /><figcaption>**Figure 11.** 21 cm line spectra toward the IC 443 shocked clumps B, C, and G. Each spectrum is from a position centered on the clump, with nearby background subtracted. At low velocities, indicated by the vertical gray dashed lines, the spectrum is strongly contaminated by absorption of the SNR synchrotron emission by cold foreground gas and should be ignored.</figcaption></figure>

fine-structure lines. On the other hand, toward the southern molecular ridge that includes the shocked clumps studied in the present project, the preshock gas has much higher density resulting in slower shocks.

To measure the amount of atomic hydrogen that may be associated with shocks into the molecular gas, we utilized the data from the Galactic Arecibo L-Band Feed Array H1 21 cm line survey (J. E. G. Peek et al. 2011). The beam size of the former, great telescope at Arecibo was 3', which samples the scales of the shocked clumps but does not resolve dense cores. IC 443 was included in the Canadian Galactic Plane Survey, with slightly higher angular resolution (1' × 3') showing similar results (H. Bakış et al. 2024), and at even higher resolution (40") with the Very Large Array (J.-J. Lee et al. 2008). We use the single-dish spectra here for their high sensitivity.

Figure 11 shows the 21 cm line spectra for the three clumps in IC 443. The 21 cm line from each shocked clump was evaluated by locating the “on” spectrum closest to the clump

center, then subtracting the average of two “off” spectra (for each clump) located 1.5 beamwidths away from the “on.” For IC 443 B and C, 21 cm emission is present over a wide range of velocities. For IC 443 G, there is no such high-velocity atomic gas. The 21 cm line profiles (both width and centroid) are very similar to the H<sub>2</sub> and H<sub>2</sub>O spectra, which strongly suggests they arise from the same shocks. The atomic H velocity distributions are most similar to model (a) from Table 3. Thus the atomic gas appears to be associated with the shock into the moderate-density gas. This makes physical sense, because for a 10 km s<sup>-1</sup> shock into the high-density gas, the temperature never rises above 5000 K. For a 60 km s<sup>-1</sup> shock into lower-density gas, on the other hand, the peak temperature is above 15,000 K and destruction of H<sub>2</sub> molecules is possible.

The observed H I column densities are well constrained, and they are only at the velocity of the shocked gas. Interpolating over the region where the 21 cm line could not be accurately measured (between the dashed lines in Figure 11), the column density of atomic H in the shocked gas is 4 × 10<sup>20</sup> cm<sup>-2</sup> for IC 443 B and C, and <1 × 10<sup>20</sup> cm<sup>-2</sup> for IC 443 G. In comparison, models (a<sub>p</sub>), (b<sub>p</sub>), and (c<sub>p</sub>) predict an H I column density of 2 × 10<sup>20</sup>, due to partial dissociation of H<sub>2</sub> in the shocked gas as shown by the green curves in Figure 6. In contrast, the amount of H I is smaller by a factor of at least 10 for the other C-type shocks. The J-type shocks, on the other hand, are completely dissociated for models (a) through (d), and they generate H I column densities at least 10 times higher than observed. The low-magnetic-field, high-density, J-type shocks (e<sub>m</sub>), (f<sub>m</sub>), and (g<sub>m</sub>) are nondissociative but have much lower H I column density than observed.

The H I column densities are thus associated with the same shocks into moderate-density gas, both in terms of the atomic column density and the velocity distribution, as the ones that explain the H<sub>2</sub> and H<sub>2</sub>O emission. For the range of velocities corresponding to molecules partially destroyed in model (a<sub>p</sub>), the spatial distribution of the 21 cm emission is shown in Figure 7 of J.-J. Lee et al. (2008) in blue and for models (b<sub>p</sub>) and (c<sub>p</sub>) in green. The shape of the molecular ridge is recognizable in that 21 cm image, with blue peaks at IC 443 B and IC 443 C where the shocks are largely along the line of sight. (IC 443 G is not noticeable in color, due to the shock
being perpendicular to the line of sight.) The 21 cm emission shows a similar velocity and spatial distribution to the H₂.

To assess the relative amount of atomic and molecular gas directly from the observations of the shocked clumps, we estimate the fraction of the hydrogen nuclei that are in atomic gas,

$$f_{\text{atomic}} = \frac{N(\text{H I})}{2N(\text{H}_2) + N(\text{H I})} \quad (3)$$

The atomic H column density is directly traced by the observed line profile, because the gas is in the optically thin limit where the integrated brightness temperature is proportional to column density. The H₂ column density was calculated from the observed brightnesses of the S(0) through S(7) lines, adding in the column density in the ground state (calculated from the excitation temperature of the lowest levels for each of ortho- and para-H₂). The column densities in levels J=2 through J=9 are directly observed, and the excitation diagram is extrapolated to measure the ground ortho- (J=1) and para- (J=0) abundances. For the shocks into dense gas, most H₂ molecules are in the J=2 and J=3 levels; only 28% of the H₂ is in the lower J=0 and J=1 levels. For the shocks into moderate-density gas, most H₂ molecules are in the J=3 and J=4 levels; only 8% are in the J=0 and J=1 levels. The dissociation fractions estimated this way are about 30% if we consider all of the H₂ associated with the same shock as the H I. If we consider only the H₂ from shocks into moderate-density gas, the dissociation fraction is about 60%. The shocks into dense gas could have as little as zero dissociation of H₂. These observational results are in agreement with the shock models used above, which predict the 21 cm line to arise from partially dissociating C-type shocks into moderate-density gas.

## 4.6. Water Abundance

In the preshock gas, the H₂O abundance can be directly estimated from the Herschel observations. For IC 443 G, the H₂O I₁₀⁻I₀₁ spectrum shows a deep, narrow dip at -4 km s⁻¹. Because this transition is from the ground state of ortho-H₂O, the absorption is due to cold H₂O along the line of sight, presumably in the preshock core. The emission profile from the warm, shocked gas is so much broader than the absorption dip that we can approximate the emission as an effective continuum, fitting a straight line from the blue to red side of the absorption dip. In this way, the dip is measured to be approximately 93% ± 2% deep. The optical depth required to cause the dip is 2.6 ± 0.3, which corresponds to a column density of ortho-H₂O of 1 × 10¹⁴ cm⁻². From a multilevel analysis of the H₂ lines, the H₂ column density of IC 443 G is approximately 10²² cm⁻² (P. Dell'Ova et al. 2024). Thus the water abundance in the preshock core (including para-H₂O) is of order 1 × 10⁻⁸. This abundance is in the range found for translucent molecular clouds (E. F. van Dishoeck et al. 2013). In our models, the preshock H₂O abundance was calculated using a steady-state chemical network, which led to initial abundances somewhat lower than the observed value.

The models that match the H₂ lines predict H₂O lines that are much brighter than observed, even after scaling by the filling factors found from the linear combinations of models to the observed H₂ lines. The velocity distribution of the H₂O is also very broad, being more similar to the H₂ S(5) line despite that

line being from a much higher energy level than that of the H₂O line. The model H₂O emission arises from a warm layer where the H₂O abundance determined from the chemical network is tremendously enhanced to x(H₂O) = 3 × 10⁻⁴. For comparison, in the preshock gas the model abundance is 4 × 10⁻¹⁰, and in the cooled postshock gas it is 1 × 10⁻⁹. A likely explanation is that the H₂O I₁₀⁻I₀₁ line is optically thick, so that our calculation (which assumed they were optically thin) overpredicted the emergent intensity. N. Fuda et al. (2023) demonstrated how to calculate the emission of an optically thick CO line from a shock, but that calculation has not been attempted for the H₂O ground-state line. To get the models to match the observations, we would require an escape probability lower than 10⁻², which is plausible but has not been verified. Because the water chemistry is complicated and relates the gas and solid phases of the interstellar medium (ISM), we do not attempt to interpret the results further, and instead only use the H₂O velocity profile as an aid to interpreting the shock dynamics.

## 5. Conclusions

The new observations and modeling in this paper reveal the velocity distribution, molecular dissociation, preshock gas conditions, and shock age of supernova blast waves in dense molecular clouds. The properties of shocks were inferred from high-spectral-resolution observations of the H₂ pure rotational lines. Using H₂ lines, as opposed to those from molecules of trace elements, is crucial because H₂ comprises the bulk of the gas mass, the low-lying lines of H₂ are the primary coolant of the shocked gas, and the H₂ abundance does not depend upon the complicated chemistry of the shocked gas. Using the pure rotational lines is important because the excitation behind the dense shocks excites H₂ primarily into the lower few rotational levels of the molecule, though the first two vibrational levels are excited for shocks into the lowest densities considered in our models. The velocity distribution of the molecular gas was directly observed using the H₂ S(1) and S(5) lines. This type of observation was not previously possible because the H₂ lines are mostly at wavelengths where the Earth's atmosphere is opaque. SOFIA flew to high enough altitude to get above 99.9% of the telluric absorption to allow the line profiles to be measured. Prior measurements of the H₂ velocity distribution used the rovibrational 1-0 S(1) line (M. Rosado et al. 2007), which traces a high energy level but has the advantage of high- angular-resolution 2D imaging, revealing substructures (both spatial and spectral) down to the angular scale of our new observations. There is limited potential to follow-up our H₂ observations, apart from challenging ground-based observa- tions of the H₂ 0-0 S(2) line, which may be possible with a mid-infrared echelle spectrograph like TEXES (J. H. Lacy et al. 2002); those observations would trace the shock into high- density gas and would constrain the ortho-to-para conversion of that gas. We found an OPR for the shock into dense gas less than the equilibrium value of 3, which implies incomplete conversion of the cold, preshock gas from its initial value (if the molecular cloud had indeed been cold and undisturbed prior to the supernova). The shock into higher-density gas can be readily observed in the near-infrared with the 1-0 S(1) line; our models showed that the shock traced by the near-infrared lines will have an OPR of 3 regardless of the preshock value.

Theoretical predictions from the Paris-Durham shock code reproduced many features of the observations, including the
broad distribution of velocities for IC 443 B and IC 443 C that require C-type (nondissociative) shocks largely along the line of sight. For IC 443 G and 3C 391, the shocks are largely perpendicular to the line of sight, making the radial velocities overlap with that of the preshock gas. Indeed, for those two shocked clumps, their broad emission lines have superposed deep, narrow absorption features, in H₂O (1₀⁻₁₀₁) for IC 443 G and HCO⁺(1–0) for 3C 391, due to cold, dense, unshocked gas from the part of the dense core between the shock front and the observer.

To explain the velocity distribution of the H₂ S(5) line, and the brightness of the S(4) through S(7) lines, preshock densities are in the range 1000–6000 cm⁻³, into which the IC 443 blast wave drives shocks of speed 80–37 km s⁻¹, respectively. To explain the velocity distribution of the H₂ S(1) line, and the brightness of S(0) through S(3) lines, shocks into gas with density 10,000–100,000 cm⁻³ gas are required. Shocks into dense gas have been indirectly inferred before from millimeter-wave observations of the similar SNR W 44 (M. Seta et al. 2004; S. Anderl et al. 2014). The need for a range of shock conditions arising from the same ram pressure is also seen on much larger scales from a galaxy intruding into a group of galaxies (P. N. Appleton et al. 2006), and is a natural consequence of the ISM having multiscale structure (as opposed to clouds of constant density).

We used 21 cm line observations to measure the velocity distribution and abundance of atomic H at the shock fronts, indicating the shocks dissociate less than half of the initial H₂. This is one of the first times that partial dissociation of H₂ has been observed in the ISM (see for similar SNR W 51C B. Koo & D. Moon 1997). Complete dissociation (and partial ionization) occurs in shocks into lower-density (<10² cm⁻³) gas; however, the lack of optical emission lines indicates that such shocks are not significant in the molecular-dominated part of SNRs like IC 443 or 3C 391 (Y. Deng et al. 2023). The Arecibo H1 observations show that the atomic H has the same velocity distribution and peaks at the same location as the shocked molecular gas. Future, high-angular-resolution observations of the 21 cm line with observatories like the Square Kilometer Array (R. Braun et al. 2015) could study the dissociation behind shocks into molecular clouds in more detail.

The properties of the shocked clumps can be readily compared to those of star-forming dense cores. The size scales of the shocked clumps B, C, and G in IC 443 are ~0.4 pc, with the mass of clump G being of order 100 M☉ (R. L. Dickman et al. 1992; G. Cosentino et al. 2022; P. Dell'Ova et al. 2024). The clumps are comparable to the widths of filaments in star-forming regions, which have a typical diameter of 0.2 pc (P. André et al. 2014). The shocked clumps contain significant structure. The spatial resolution of our observations is 0.03 pc, which is comparable to (but somewhat larger than) the size of an individual protostar (Bonnor-Ebert condensation) (P. C. Myers et al. 2000) with mass of order 1 M☉. These sizes are at the lower end of prestellar cores, for which lifetimes are estimated to be 1–10 Myr (N. E. Jessop & D. Ward-Thompson 2000). A likely scenario, then, is that the shocked clumps were pre-existing dense cores of the parent molecular cloud of the progenitor star. The blast waves reached the dense cores about 7000 yr after the supernova explosion, and they have been propagating into them for about 3000 yr.

The pre-existing dense cores were nearly massive enough to be self-gravitating, if they had velocity dispersion less than

0.4 km s⁻¹ (though observed line profiles of the preshock molecular gas are typically larger than this). In the postshock gas, where velocity gradients of several km s⁻¹ are present, the range of shock directions and amplification of turbulence are likely to keep them unbound or drive them from being gravitationally bound to being unbound. The compression in the postshock gas could be a trigger for star formation. For the C-type shock models used in this paper (Table 3, top) that can explain the line profiles, the shock compression is predicted to be a factor of 5–10. In contrast, for shocks into the same gas density but with low magnetic field (J-type shocks; Table 3, bottom), the compression is a factor of 30 to 170. The strong magnetic field, which is required to produce the observed, broad line profiles, keeps the compression low and reduces the shock front's potential to trigger star formation. The shocks into IC 443 and 3C 391 molecular clumps are into dense cores that were not imminently forming stars, and after the passage of the shock they may be even less likely to do so.

Future theoretical work on interaction of shock fronts with 2D or 3D structures with configurations consistent with the best high-resolution observations would advance our understanding of the relationship between supernova shocks and pre-existing molecular cores. Observational studies at higher angular resolution and including multiple energy levels, for example in the infrared with JWST (L. E. Kristensen et al. 2023) or submillimeter with ALMA or NOEMA (N. Dello Russo et al. 2008), will trace the shock fronts and provide a progressively more detailed view of how violent events at the end of stellar evolution affect the potential of their parent molecular clouds to form stars.

## Acknowledgments

Based in part on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA was jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. Financial support for the observational and theoretical work in this article was provided by NASA through award #07_0007 issued by USRA to the Space Science Institute. Herschel was an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.

**Facilities:** SOFIA, Herschel, Spitzer.

**Software:** sofia_redux (M. Clarke et al. 2015), astropy (Astropy Collaboration et al. 2022), matplotlib (J. D. Hunter 2007).

## ORCID iDs

William T. Reach https://orcid.org/0000-0001-8362-4094

Le Ngoc Tram https://orcid.org/0000-0002-6488-8227

Curtis DeWitt https://orcid.org/0000-0002-6528-3836

Pierre Lesaffre https://orcid.org/0000-0002-8447-4398

Benjamin Godard https://orcid.org/0000-0003-0060-8887

Antoine Gusdorf https://orcid.org/0000-0002-0354-1684

## References

* Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 712, 459

* Anderl, S., Gusdorf, A., & Gusten, R. 2014, A&A, 569, A81

* André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 27

* Appleton, P. N., Xu, K. C., Reach, W., et al. 2006, ApJL, 639, L51
Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Bakış, H., Paylı, G., Aktekin, E., Sano, H., & Sezer, A. 2024, MNRAS, 532, 2570

Binney, J., & Merrifield, M. 1998, Galactic Astronomy (Princeton, NJ: Princeton Univ. Press)

Braun, R., Bourke, T., Green, J. A., Keane, E., & Wagg, J. 2015, in Advancing Astrophysics with the Square Kilometre Array, Proceedings of Science, ed. R. Braun et al. (Trieste: SISSA), 174

Braun, R., & Strom, R. G. 1986, A&A, 164, 193

Brogan, C. L., Frail, D. A., Goss, W. M., & Troland, T. H. 2000, ApJ, 537, 875

Burton, M. G., Geballe, T. R., Brand, P. W. J. L., & Webster, A. S. 1988, MNRAS, 231, 617

Cesarsky, D., Cox, P., Pineau des Forêts, G., et al. 1999, A&A, 348, 945

Clarke, M., Vacca, W. D., & Shuping, R. Y. 2015, in ASP Conf. Ser. 495, Astronomical Data Analysis Software and Systems XXIV, ed. A. R. Taylor & E. Rosolowsky (San Francisco, CA: ASP), 355

Cosentino, G., Jiménez-Serra, I., Tan, J. C., et al. 2022, MNRAS, 511, 953

Crutcher, R. M. 2012, ARAA, 50, 29

de Graauw, T., Helmich, F. P., Phillips, T. G., et al. 2010, A&A, 518, L6

Dello Russo, N., Vervack, R. J., Jr., Weaver, H. A., et al. 2008, ApJ, 680, 793

Dell'Ova, P., Gusdorf, A., Gerin, M., et al. 2024, A&A, 685, A69

Deng, Y., Zhang, Z.-Y., Zhou, P., et al. 2023, MNRAS, 518, 2320

Denoyer, L. K. 1979, ApJ, 232, L165

Dickman, R. L., Snell, R. L., Ziurys, L. M., & Huang, Y.-L. 1992, ApJ, 400, 203

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton, NJ: Princeton Univ. Press)

Draine, B. T., & McKee, C. F. 1993, ARAA, 31, 373

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

Fuda, N., Tram, L. N., & Reach, W. T. 2023, ApJ, 944, 17

Godard, B., Pineau des Forêts, G., Lesaffre, P., et al. 2019, A&A, 622, A100

Green, D. A. 2019, JAPA, 40, 36

Green, S., Maluendes, S., & McLean, A. D. 1993, ApJS, 85, 181

Guillet, V., Pineau Des Forêts, G., & Jones, A. P. 2007, A&A, 476, 263

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

Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18

Huang, Y. L., Dickman, R. L., & Snell, R. L. 1986, ApJL, 302, L63

Hunter, J. D. 2007, CSE, 9, 90

Indriolo, N., Blake, G. A., Goto, M., et al. 2010, ApJ, 724, 1357

Jessop, N. E., & Ward-Thompson, D. 2000, MNRAS, 311, 63

Koo, B.-C., & Moon, D.-S. 1997, ApJ, 475, 194

Kristensen, L. E., Godard, B., Guillard, P., Gusdorf, A., & Pineau Des Forêts, G. 2023, A&A, 675, A86

Lacy, J. H., Richter, M. J., Greathouse, T. K., Jaffe, D. T., & Zhu, Q. 2002, PASP, 114, 153

Lee, J.-J., Koo, B.-C., Snell, R. L., et al. 2012, ApJ, 749, 34

Lee, J.-J., Koo, B.-C., Yun, M. S., et al. 2008, AJ, 135, 796

Lee, Y.-H., Koo, B.-C., & Lee, J.-J. 2020, AJ, 160, 263

Lehmann, A., Godard, B., Pineau Des Forêts, G., & Falgarone, E. 2020, A&A, 643, A101

Lesaffre, P., Chièze, J. P., Cabrit, S., & Pineau Des Forêts, G. 2004, A&A, 427, 157

Lesaffre, P., Pineau Des Forêts, G., Godard, B., et al. 2013, A&A, 550, A106

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 500, 259

Myers, P. C., Evans, N. J., II., & Ohashi, N. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & N. Ohashi (Tucson, AZ: Univ. Arizona Press), 217

Neufeld, D. A., Dewitt, C., Lesaffre, P., et al. 2024, ApJ, 976, 3

Neufeld, D. A., Gusdorf, A., Gusten, R., et al. 2014, ApJ, 781, 102

Neufeld, D. A., Hollenbach, D. J., Kaufman, M. J., et al. 2007, ApJ, 664, 890

Olbert, C. M., Clearfield, C. R., Williams, N. E., Keohane, J. W., & Frail, D. A. 2001, ApJ, 554, L205

Peek, J. E. G., Heiles, C., Douglas, K. A., et al. 2011, ApJS, 194, 20

Pickett, H. M., Poynter, R. L., Cohen, E. A., Delitzky, M. L., & Pearson, J. C. 2003, Submillimeter, Millimeter, and Microwave Spectral Line Catalog: Revision 4 (JPL Publication 80-23) (Pasadena, CA: JPL), 80

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Reach, W. T., & Rho, J. 1996, A&A, 315, L277

Reach, W. T., & Rho, J. 1998, ApJL, 507, L93

Reach, W. T., & Rho, J. 1999, ApJ, 511, 836

Reach, W. T., Tram, L. N., Richter, M., Gusdorf, A., & DeWitt, C. 2019, ApJ, 884, 81

Richter, M. J., DeWitt, C. N., McKelvey, M., et al. 2018, JAI, 7, 1840013

Richter, M. J., Graham, J. R., Wright, G. S., Kelly, D. M., & Lacy, J. H. 1995, ApJL, 449, L83

Rosado, M., Arias, L., & Ambrocio-Cruz, P. 2007, AJ, 133, 89

Seta, M., Hasegawa, T., Sakamoto, S., et al. 2004, AJ, 127, 1098

Sinitsyna, V. G., & Sinitsyna, V. Y. 2023, Univ, 9, 98

Snell, R. L., Hollenbach, D., Howe, J. E., et al. 2005, ApJ, 620, 758

Tavani, M., Giuliani, A., Chen, A. W., et al. 2010, ApJL, 710, L151

Temi, P., Hoffman, D., Ennico, K., & Le, J. 2018, JAI, 7, 1840011-186

Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge: Cambridge Univ. Press)

Torres, D. F., Rodriguez Marrero, A. Y., & de Cea Del Pozo, E. 2008, MNRAS, 387, L59

van Dishoeck, E. F., Herbst, E., & Neufeld, D. A. 2013, ChRv, 113, 9043

van Dishoeck, E. F., Jansen, D. J., & Phillips, T. G. 1993, A&A, 279, 541

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

Wilner, D. J., Reynolds, S. P., & Moffett, D. A. 1998, AJ, 115, 247

Young, E. T., Becklin, E. E., Marcum, P. M., et al. 2012, ApJL, 749, L17

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