© The Authors 2024

1 Introduction

Supernova explosions that occur in the vicinity of molecular clouds drive shocks through the cold, dense molecular gas, providing mechanical and energetic feedback over tens of parsecs and thus playing a key role in the evolution of galaxies. About ~30 kyr ago, one such explosion produced the Galactic supernova remnant (SNR) IC443, an appropriate candidate to study this type of SNR-cloud interaction during the late radiative phase (e.g., Reach et al. 2019). In IC443, slow shocks (υ_s ≪ 10² km s⁻¹) are still propagating into nearby molecular regions, where they modify the physical and chemical state of the gas, thereby altering the local star formation rate. Estimating the total mass of shocked clumps of gas in IC443 gained a particular interest with the detection of teraelectronvolt emission along the molecular shell of the remnant (MAGIC, Albert et al. 2007; VERITAS, Acciari et al. 2009; Fermi, Abdo et al. 2010), which is the product of interactions between cosmic rays (CRs) accelerated in the early phase of the SNR and the local ISM (Tavani et al. 2010; Tang & Chevalier 2015; Xu 2021). Molecular hydrogen (H₂) is the most abundant molecule in the universe, and it is one of the main coolants in structures processed by slow (5–25 km s⁻¹) molecular shocks (Godard et al. 2019; Kristensen et al. 2023, and references therein). Therefore, it constitutes an effective tracer of the molecular gas phase during the radiative phase of an SNR as well as a means to evaluate the energetic impact of the SNR on its environment. In addition, measuring the H₂ column density (hereafter $N_{{\text{H}}_2}$) allows one to proportionately constrain chemical abundances in a molecular environment. Consequently, the value of a measurement of $N_{{\text{H}}_2}$ is at least four-fold, providing constraints on the molecular density, mass and temperature for future studies of energetics, star formation, γ-ray emission, and chemistry in IC443. However, since the two lowest rotational transitions of H₂ have upper-level energies E/k ≃ 510 K and 1015 K above the ground level (Dabrowski 1984), the emission from rovibrational lines of H₂ mainly traces warm shocked gas (T > 100 K). For this reason, cold molecular hydrogen is essentially invisible, and it remains a challenge to weigh the total mass of a molecular clump from the emission of H₂ infrared lines only.

In the following, we show that this limitation can be overcome by constraining the temperature distribution of the molecular gas based on ¹²CO and ¹³CO (J = 1−0, J = 2−1, J = 3−2) observations of the cold phase in the same region (Dell’Ova et al. 2020). We first briefly present the SNR IC443 and previous studies of this object (Sect. 2). The observational details of the Spitɀer/IRS spectral-line maps are presented in Sect. 3. In order to measure the column density of H₂ and the temperature of the molecular gas, we used a two-temperature approach, which is described in Sect. 4. Then, we adopted a second method, outlined in Sect. 5, in which we reproduced the emission of H₂ lines with a power-law distribution of gas temperatures. Finally, the conclusions are summarized in Sect. 6.

2 The IC443 supernova remnant

The SNR IC443 (G189.1+3.0) is located at a distance of 1.8 kpc (Ambrocio-Cruz et al. 2017; Yu et al. 2019). Notably, IC443 is one of the most studied mixed-morphology evolved SNR for its interaction with nearby interstellar gas and molecular clouds (see Fig. 1, Lee et al. 2008, 2012; Dell’Ova et al. 2020; Ritchey et al. 2020 and references therein). Recent studies have provided an estimate of approximately 30 kyr for the age of IC443 (Ambrocio-Cruz et al. 2017). Observations of radio, optical, infrared, and X-ray continuum have indicated a shell-like structure where the SNR interacts with an inhomogeneous interstellar environment (Braun & Strom 1986; Castelletti et al. 2011; Alarie & Drissen 2019; Reach et al. 2019; Troja et al. 2006). In the southern part of the remnant, the SNR blast is encountering a molecular cloud, in which shock-excited broad molecular lines are detected (Denoyer 1979a; Burton et al. 1988; Dickman et al. 1992; van Dishoeck et al. 1993). Along the southwestern H2 ridge of the SNR, several shocked clumps have been identified, including the shocked clump G (Dickman et al. 1992). The molecular emission, shock chemistry, and kinematics in IC443G were investigated by Ziurys et al. (1989), Turner et al. (1992), van Dishoeck et al. (1993), Tauber et al. (1994), Snell et al. (2005) and Cosentino et al. (2022). The shock in IC443G is mostly perpendicular to the line of sight, and no single shock type can account for the observations. Several species were detected, including SiO, CS, SO, H₂CO, H₂O, CN, NH₃, and C I. Infrared spectral observations of H₂ in IC443G were performed from the ground by Richter et al. (1995), Kokusho et al. (2020), and Deng et al. (2023); in flight by SOFIA (Reach et al. 2019); and from space by the Infrared Space Observatory (Cesarsky et al. 1999), Spitɀer (Neufeld et al. 2007; Noriega-Crespo et al. 2009), and AKARI (Shinn et al. 2011). However, to this day there is no systematic pixel-per-pixel analysis of the emission from pure rotational lines of H₂ in IC443G.

3 Spitzer/IRS spectral-line maps

3.1 Observations by Spitzer/IRS

The IC443G clump was mapped with Spitɀer/IRS by D.A. Neufeld et al. in 2007 and later by Noriega-Crespo et al. (2009). Neufeld et al. (2007) mapped eight pure rotational transitions of H₂ in a ~2′ × 2′ field of observations (the spectroscopic parameters of the transitions are shown in Table 1). In their initial proposal, Neufeld et al. (2004) aimed to perform spectral-line mapping toward four SNRs (IC443, 3C391, W28 and W44) in order to constrain the molecular shock parameters in these regions. The analysis of the IC443C clump was published by Neufeld et al. (2007) and Yuan & Neufeld (2011), but the IC443G maps obtained with the SL module present several horizontal and vertical stripe-like artifacts and patterns (shown in Fig. A.1). The location of the ~2′ × 2′ field of H₂ observations is shown in Fig. 2, with respect to the IRAM-30m field of ¹²CO observations presented by Dell’Ova et al. (2020). The Spitɀer/IRS field of observation was centered on the shocked clump, with no contribution from the quiescent cloud to the northwest. To measure the line intensities, the lines were fit with a Gaussian plus a first-order baseline. The S(0) to S(7) spectralline maps were projected on 100 × 100 grids (pixel size: 1″.6) that were spatially calibrated (i.e., for any (i, j) a pixel (x₀, y₀) of the Sⁱ map corresponds to the same line of sight as the pixel (x₀, y₀) in the S^j map).

3.2 Morphology

The maps of the spectral emission from lines S(0) to S(7) are represented in Fig. 3. The 2′ × 2′ region was centered on a shocked molecular clump interacting with the SNR blast, labeled ‘IC443G’ (Dickman et al. 1992). The spatial correlation between the different H₂ lines is significant: the bright, shocked structure along the northwest-southeast axis is shown in each panel of the figure, albeit with slightly varying structures at the smallest available angular scales. The shocked clump ‘G’ is encompassed by the southern H₂ ridge described by Burton et al. (1988) along the shell A. The higher spatial resolution offered by the spectralline maps S(4) to S(7) shows that the apparent extended emission in the S(0) map resolves into several knots of shocked gas to the west of the main structure.

Fig. 1 Overview of the SNR IC443. Top-left panel: schematic drawing of the SNR morphology. The two-shell model proposed by Troja et al. (2006) is shown as a colored background image: the gray torus represents the molecular cloud discovered by Cornett et al. (1977) and is highlighted by white curved lines based on the higher resolution observations by Lee et al. (2012). The blue hemisphere labeled “subshell A” represents the shock front in the eastern region. In the northeast, where it has been confined by the encounter with the neutral H I cloud of Denoyer (1978), the ionic shock front is traced by optical, infrared and very soft X-ray emission (Alarie & Drissen 2019, dashed-dotted black lines). In the southern region, the thick black curved lines represent the areas of interaction with the molecular cloud and the shocked H2 ridges of Burton et al. (1988). The labels and positions of the different molecular clumps identified by Denoyer (1979b) and Dickman et al. (1992) are indicated with capital letters (B, C, D, E, F, G, H). The red hemisphere labeled “subshell B” represents the shock front in the western region, where it is expanding in a homogeneous and less dense medium. The position of the PWN discovered by Olbert et al. (2001) is indicated by a white star. The extension of the radio continuum halo is shown by dashed black lines (Lee et al. 2008, later confirmed by Castelletti et al. 2011). Locations and extensions of the four gamma-ray peaks at different energy ranges are indicated: EGRET centroid (∘), MAGIC centroid (□), VERITAS centroid (◊), and Fermi LAT centroid (∆), respectively measured by Esposito et al. (1996), Albert et al. (2007), Acciari et al. (2009) and Abdo et al. (2010). The respective localization errors are shown as crosses of size 0.5 σ. We note that the high-energy emission is extended, and hence the relation between Fermi LAT, VERITAS, and MAGIC observations may be better characterized by the extent of the emission rather than the centroids. The sizes of typical instrumental beams are indicated in the bottom-right corner of the figure. The other panels represent color composite images of IC443, with H I 21 cm line emission (Lee et al. 2008) in yellow and 24 µm (Spitɀer/MIPS) emission in white in all three panels. The third color channel represents a distinct signal in each panel: top-right panel: 12CO J = 1–0 (Lee et al. 2012, green); bottom-left panel: DSS optical observations (York et al. 2000, blue); bottom-right panel: XMM-Newton observations at 1.4–5 keV (Greco et al. 2018, purple).

3.3 Spitzer/IRS SH versus SL map comparison

The H₂ υ = 0–0 S(2) line was mapped both by the SH and SL modules. Hence, we could estimate the uncertainty introduced by the noisy vertical stripe (see Fig. A.1). The relative difference between the total flux of the S(2)[SH] and S(2)[SL] maps is ~5.5% (with I^SL > I^SH). We performed a pixel-per-pixel comparison of the S(2) spectral-line maps to estimate the uncertainty associated with the data (see Fig. 4). The left panel displays a small number of obvious outliers in the SL map (for I^SL = [100−125] × 10⁻⁵ erg cm⁻² s⁻¹ sr⁻¹ and ƒ^SH = [15–35] × 10⁻⁵ erg cm⁻² s⁻¹ sr⁻¹) that are related to the noisy vertical stripe. These outliers were removed by the preprocessing transformations that we applied to the spectral-line maps (see Appendix A, and right panel of Fig. 4). In order to obtain a first estimate of the systematic and random uncertainties on the integrated intensity, we measured the statistical variables that describe the relation between the SL map and the SH map. First, we modeled the sample of data points ƒ^SL = f (I^SH) by a linear function x ↦ ax + b. Then, we computed the standard deviation between the model I^mod and the SL data points I^SL. The results are presented in Table 2. Our pre-processing transformations enhanced the signal-to-noise ratio and reduced the systematic uncertainties (with respect to the raw SL maps).

Fig. 2 Finderchart of the Spitɀer/IRS field of observations. Left: 12CO J = 2−1 10′ × 10′ map of the peak temperature obtained with the IRAM-30m telescope at an angular resolution of 11″.2 (Dell’Ova et al. 2020). Right: H2 υ = 0–0 S(1) map of the integrated line intensity obtained with Spitɀer/IRS.

Fig. 3 H2 υ = 0–0 S(0) to S(7) pre-processed maps obtained with Spitɀer/IRS toward the region IC443G. The S(1) to S(7) maps were convolved to the spatial resolution of the S(0) map (8″.3, see Sect. 4.2). The S(0) line was mapped by the LH module, the S(1) and S(2) lines were mapped by the SH modules, and the rest were mapped by the SL module. The S(2) line was mapped by both the SH (top row) and SL (second row) modules. The coordinates of the field of observations are indicated in Fig. 2. Raw maps are shown in Fig. A.1. The array of white circles indicates the positions in which we extracted the signals of the S(0) to S(7) spectral-line maps to produce the population diagram mosaics presented in Fig. 5 (model 1: standard LTE approach) and Fig. 9 (model 2: H2 thermal admixture model). The size of the circles corresponds to the angular resolution of the S(0) observations, that is, 8″.3 in diameter.

Fig. 4 Pixel-per-pixel comparison of the integrated intensity S(2) maps measured by the Spitɀer/IRS SH and SL modules before (left panel) and after pre-processing (rightpanel). The solid lines represents the 1:1 relation expected if the SH and SL data arrays were identical, and the dashed lines represent the ±25% relative errors.

3.4 Spitɀer/IRS versus ISO comparison

The H₂ S(2) to S(7) rotational lines were also observed with ISO (Cesarsky et al. 1999). In order to check the reliability of the data, we compared the Spitɀer/IRS and ISO measurements toward the positions A, B, and C (Table 1 of Cesarsky et al. 1999). We found an average relative error of 26% between the IRS and ISO spectra. To take this discrepancy in our study into account, we defined the uncertainties on any IRS measurements I_ν∆ν as follows: $$\sigma =\sqrt{{\left(0.26\times I_v\text{Δ}v\right)}^2+{\sigma }_{\text{MAD}}^2},$$(1)

where σ_MAD is the mean absolute deviation measured in the background of each spectral map. Therefore, the total uncertainty is representative of both the intrinsic noise of the images (σ_MAD) and the absolute flux calibration errors suggested by the discrepancy between the ISO and IRS measurements (26%).

4 Two-temperature model

We aim to produce maps of the column density and excitation temperature of the molecular phase toward the shocked clump IC443G. Our main goal is to provide a measurement of the H₂ column density $\left(N_{{\text{H}}_2}\right)$ that accounts for the mass of cold gas. To this end, we first performed a two-temperature population diagram analysis of H₂ lines in each pixel of the Spitɀer/IRS spectral-line maps.

4.1 Assumptions

(i) Line opacity. We assumed that the emission of the pure rotational lines S(0) to S(7) is in the optically thin regime. Given the low Einstein coefficients A_ul of the quadrupolar transitions of H₂ (see Table 1), the emission of H₂ rovibrational transitions remains optically thin even in proto-planetary disks, where the H₂ column density can be larger than 10²³ cm⁻² (e.g., Bitner et al. 2008).

(ii) Thermalization. We assumed that the gas is at local thermodynamic equilibrium (LTE). This condition is verified if the local medium is sufficiently dense for a collisional equilibrium to settle, that is, if the following condition on the local density n is satisfied: $$n\gg n_{\text{crit }}$$(2) $$n_{\text{crit }}=\frac{A_{\text{ul}}}{{{{\displaystyle \sum }}^{\text{}}}^{\text{}}{C}_{\text{ul}}},$$(3)

where A_ul and C_ul are respectively the Einstein coefficient of spontaneous emission and the collisional de-excitation rate. Wrathmall et al. (2007) presented the results of calculations of the collisional rate coefficients of H₂ for a kinetic temperature T = 1000 K, T = 2500 K, and T = 4500 K. To test our assumption, we combined their calculations of C_ul with the A_ul coefficients shown in Table 1 to estimate the critical density n_crit of the H₂ pure rotational lines. Depending on the kinetic temperature, the critical density of the S(7) line is n_crit ~ 1 × 10⁴−4 × 10⁵ cm⁻³ (C₉₇ (1000 K) = 4.6 × 10⁻¹³ cm³ s⁻¹; C₉₇(4500 K) = 1.9 × 10⁻¹¹ cm³ s⁻¹). The critical densities of the other lines are lower. Cesarsky et al. (1999) estimated a pre-shock density ~10⁴ cm⁻³ in the clump G. Thus, the low-J H₂ lines are likely thermalized, although the high-J H₂ lines might not be (see also van Dishoeck et al. 1993; Shinn et al. 2011). Assuming a density compression ratio of 10, all the lines would be thermalized in the post-shock.

(iii) Extinction. The scattering and absorption of incident photons by interstellar dust is expected to cause the extinction of H₂ lines (e.g., Draine 2011). The S(0) to S(7) lines are characterized by λ_ul = 5.5–28.2 µm (Table 1), a wavelength range in which the extinction is potentially non-negligible, depending on the magnitude of the visual extinction (A_υ). For a given wavelength λ, the observed intensity ƒ is related to the emitted signal I⁰ by the following equation: $$I(\lambda )=I^0(\lambda )\times {10}^{\left(-A_{\lambda }/2.5\right)},$$(4)

where A_λ is the extinction at wavelength λ. This equation is correct when the extinction is caused by a screen external to the source, which is the case here. The quantity A_λ/A_υ can be found in Rieke & Lebofsky (1985) for λ = 4.8, 7.0, 9.7, 14.3, 19, 30 µm. We performed a bicubic interpolation of this sample of data points to obtain the quantity A_λ/A_υ for the exact wavelength corresponding to each H2 pure rotational line. Following Shinn et al. (2011), we adopted a visual extinction A_υ = 10.8, inferred from A_2.12 = 1.3, obtained by Richter et al. (1995), employing the ‘Milky Way, R_υ = 3.1’ extinction curve. Our correction of the extinction alters the flux by a few percentage points only.

(iv). We assumed a standard H₂ ortho-to-para ratio (o/p = 3). We discuss the impact of choosing a different ortho-to-para ratio below.

4.2 Pixel-per-pixel method

Using Gaussian kernels, we first convolved the S(1) to S(7) maps (λ = 5.5−17 µm) to the spatial resolution of the S(0) map (8″.3 at λ = 28.2 µm). We then applied the following steps in each pixel (x, y) of the maps: (1) Correct the integrated intensity of each pixel for dust extinction (Eq. (4)). (2) Convert the integrated intensity I into an upper-level population N_up (Eq. (B.3)). (3) Build the population diagram, that is, log(N_up/𝑔_up) with respect to E_up, where 𝑔_up includes the nuclear spin degeneracy (2I + 1)(2J + 1), where I = 0/1 for J even/odd (cf. Table 1), and infer the column density and temperature from the linear fit of the data points.

4.3 Results: Population diagrams

We performed a systematic population diagram analysis in each pixel of the Spitɀer/IRS spectral-line maps. We performed two separate linear fits of the H₂ energy levels respectively traced by the S(0), S(1), S(2) and S(3), S(4), S(5), S(6), S(7) lines to take into account temperature stratification (see the following developments). The results of our analysis are presented in Fig. 5. We extracted nine population diagrams in different locations of the IC443G clump (these locations are indicated in Fig. 3, first panel). We selected regions with low signal-to-noise ratios (first column of the mosaic) and regions with high signal-to-noise ratios (third column of the mosaic). The diagrams are shown in Fig. 5. In order to check the quality of the best-fits, we performed a chi-square test on the data points, defined by the reduced chi-square (or chi-square per degree of freedom): $${\chi }^2=\frac{1}{N-p}{\displaystyle \underset{J_{\text{up }}=2}{\overset{9}{{{\sum }^{\text{}}}^{\text{}}}}}\frac{{\left(N_{\text{up}}-f\left(E_{\text{up}}\right)\right)}^2}{\sigma {\left(N_{\text{up}}\right)}^2},$$(5)

where N_up is the column density; f is the adopted linear model (which is distinct for the energy levels traced by the S(0)–S(2) and S(3)–S(7) lines, following our decomposition into two components); σ(N_up) is the uncertainty on N_up derived from Eq. (1); and N = 8 is the number of data points. Lastly, p = 4 is the number of free parameters (i.e., the excitation temperature and total column density for the warm and cold components). In the extracted population diagrams, χ² varies between 0.95 and 2.67 (see Fig. 5). The value of χ² increases for the high signal-to-noise beams (3, 6, 9) since the small error bars make the data harder to reproduce with our model. The linear fit of the J_up = 5–9 upper levels (gray curve in Fig. 5) systematically fails to reproduce the column density of the population level J_up = 2. Conversely, the linear fit of the remaining data points (J_up = 2–4, red curve) does reproduce the lower energy levels, although it is completely off for higher energy levels. Classically, at least two excitation temperatures are required to reproduce the H2-level populations with a ‘cold’ component (dominating levels J_up = 2–4) and a ‘warm’ component (dominating levels J_up = 3–9). We found, however, that the levels J = 3 and J = 4 are not well reproduced by our two-temperature model, which we interpret as evidence that a continuous distribution of temperature would better represent the data than a bimodal distribution. Naturally, the total column density is given by the sum of the column densities inferred from each component (‘warm’ and ‘cold’).

4.4 Discussion

We have obtained an estimate of the total H₂ column density inferred from the emission of the S(0) to S(7) lines (see Table 3, first and second rows). The fit of the S(0), S(1), S(2) (in red in Fig. 5) traces the ‘cold component’ of H₂ (T_ex = 370–470 K, $N_{{\text{H}}_2}~{10}^{20.6}$ cm⁻²) and the bulk of the molecular mass. Similarly, the fit of the S(3), S(4), S(5), S(6), S(7) (in gray in Fig. 5) traces the ‘warm component’ of H₂ (T_ex = 740–840 K, $N_{{\text{H}}_2}~{10}^{19.7}$ cm⁻²). As expected, the transitions of H₂ cannot be reproduced by a single excitation temperature; hence, our results indicate that the medium is characterized by thermal stratification along the line of sight. Both the column density and excitation temperature increase toward the shocked clump: low signal-to-noise beams (1, 4, 7) are characterized by a column density $N_{{\text{H}}_2}~{10}^{20.5}$ cm⁻² and temperature T_ex ~ 370 K (for the cold component), while the high signal-to-noise beams (3, 6, 9) are characterized by a column density $N_{{\text{H}}_2}~{10}^{21}$ cm⁻² and temperature T_ex ~ 470 K (for the cold component).

The energies E_up of the upper levels (J = 2–9) are 510–7200 K (Table 1). It is unlikely that the result of our analysis yields the total H₂ column density since the rotational energy levels J = 0 and J = 1 are expected to hold the greater part of the ‘cold mass’ of molecular hydrogen and perhaps another component should be considered for the first two levels. In other words, it is likely that these two levels are not thermalized with the others. It is well known that the ‘standard’ population diagram analysis of H₂ pure rotational lines is a poor method to measure the total column density. Roussel et al. (2007) presented results of the analysis of warm H₂ lines in the Spitzer SINGS (Spitzer Infrared Nearby Galaxy Survey) galaxy sample. Under a conservative assumption about the distribution of temperatures, they found that the column densities of warm H2 derived from Spitzer/IRS observations amount to between 1% and 30% of the total mass of H2 (their estimate of the mass of cold H2 was derived from velocity-integrated ¹²CO J = 1–0 observations). Therefore, it is possible that our estimate of the H₂ column density from the sum of the 'gray' and 'red' fits (Fig. 5) lies between 1% and 30% of the true column density.

Fig. 5 Population diagrams extracted from the positions shown in Fig. 3. The corresponding positions are indicated in the top-left corner of each panel. The black data points represent the Nup measurements and their uncertainties. The gray curve represents the linear model obtained for the S(3) to S(7) transitions. The red curve represents the linear model obtained for the S(0) to S(2) transitions. Gray circles represent the sum of the two models for each transition. The errors associated with the best linear models are represented by a filled area around the curves. The corresponding excitation temperatures are shown in the bottom-left corner of the diagrams ('cold' component in red and 'warm' component in gray). The total column density is also shown (in black). The results of the χ2 test (see text and Eq. (5)) are presented in the top-right corner of the diagrams. The order of the diagrams (from left to right and top to bottom) is the same as the order of the white circles in Fig. 3.

5 Thermal admixture model

5.1 Assumptions

Following the method presented by Neufeld & Yuan (2008, referred to as the H₂ thermal admixture model), we assumed that the molecular hydrogen probed along a line of sight has a column density that can be described by the following power-law distribution with respect to the kinetic temperature of H₂ (denoted T): $$\text{d}N=\text{Λ}{T}^{-\text{Γ}}\text{d}T,$$(6)

where Λ and Γ are two constants. If Γ is known, Λ can be determined by integrating Eq. (6) between T_min and T_max: $$\text{Λ}=\frac{N_{\text{tot}}\left(T\ge T_{\min }\right)(\text{Γ}-1)}{T_{\text{min}}^{1-\text{Γ}}-T_{\text{max}}^{1-\text{Γ}}}.$$(7)

This power-law description is driven by the observation that the population diagrams shown in Fig. 5 are better described by two components than a single T_ex (the 'cold' and 'warm' components). Equation (6) is a continuous generalization of this multi-temperature analysis. We assumed that the emission of the H₂ lines is produced by a mixture of gas temperatures (both along the line of sight and possibly across the plane of sky probed by the beam of Spitzer/IRS). Our adopted description (Eq. (6)) relies on the physical assumption that the shocked medium is stratified in temperature in a way that can be reproduced by a power-law distribution between a minimum gas temperature (T_min) and a maximum gas temperature (T_max). This assumption is supported by the expected temperature profile in a parabolic C-type bow shock in steady state (Shinn et al. 2011; Smith & Brand 1990). If we fix T_min and T_max, we have the same amount of parameters to fit as in the previous section (a_lin, b_lin in the standard population diagram model and Λ, Γ in this model). Once T_min and T_max are fixed, the two remaining parameters of this model are the steepness of the distribution (the power-law index Γ) and the total H₂ column density N_tot(T ≥ T_min). The value of T_max has a marginal effect on the resulting H₂ column density since most of the mass is held by cold gas. Based on the excitation temperature maps found previously (see Fig. 5), we first fixed T_max = 1500 K. We then verified that higher temperatures do not modify the results (up to T_max = 4000 K, the temperature at which H2 is rapidly dissociated in a shock, Le Bourlot et al. 2002). Determining the correct value of T_min is paramount since it controls the column density of H₂. We can estimate the lower boundary of the temperature distribution based on the results of our previous pixel-per-pixel radiative transfer analysis of ¹²CO and ¹³CO J = 1−0, J = 2−1 and J = 3−2 lines (E_up = 5.5, 16.6, 33.2 K), where we found that the temperature of the cold gas is ${25}_{-10}^{+15} \text{K}$ across the inner regions of the shocked clump (Dell'Ova et al. 2020). Hence, we fixed T_min(H₂) = 25 K accordingly with the average temperature of the molecular gas determined toward the shocked clump, and then we estimated the error on the derived H₂ column density by comparing our results with complementary models obtained with T_min = 15 K and T_min = 40 K.

Fig. 6 Average chi-square (‘χ2’, i.e., the average chi-square across the map, indicated by red diamonds) and average column density ($\overline{‘N_{H_2}’}$, indicated by turquoise squares) with respect to the adopted local density. The best-fit model was found with a density n = 2.5 × 105 cm−3, corresponding to χ2 = 1.13 and $\overline{N_{{\text{H}}_2}}=1.64\times {10}^{22} {\text{cm}}^{-2}$

5.2 Method: LTE and non-LTE admixture models

Based on the power-law description of the temperature distribution, we experimented with two distinct admixture models. The first model assumes LTE. The assumptions and discretization method are described in Appendix B. The second model is based on statistical equilibrium (hereafter SE). We used the code developed by Gusdorf et al. (2008) and Godard et al. (2019). This code computes the excitation and radiative transfer in a collection of gas slabs defined by the user. Our modeled gas consists of H₂, H, He, and electrons only, with n(H) = 10⁻² n(H₂), n(He) = 0.2 n(H₂), and n(e⁻) = 10⁻⁷ n(H₂). In order to account for possible deviations from LTE, we performed a series of calculations with variable H₂ local density, between $n_{{\text{H}}_2}\text{\hspace{0.17em}=\hspace{0.17em}}{10}^4{\text{cm}}^{-3}$ and $n_{{\text{H}}_2}\text{\hspace{0.17em}=\hspace{0.17em}}{10}^7{\text{cm}}^{-3}$ (see Fig. 6, where we show the variation of the average chi-square and average column density within the map with respect to the adopted local density). For each model, the remaining free parameters were the column density and the kinetic temperature.

We produced a discrete parameter grid [Γ, N_tot(T ≥ T_min)] for each model. We then performed a chi-square test to compare the observations of H₂ lines to the integrated intensity models. We measured the following quantity in each cell of the model grid: $${\chi }^2=\frac{1}{N-p}{\displaystyle \sum _{J_{\text{up}}=2}^9\frac{{\left(s_{\text{J}\to \text{J}-2}-m_{\text{J}\to \text{J}-2}\right)}^2}{{\sigma }_{J\to J-2}^2}}$$(8)

where s_J→J–2 is the integrated intensity measurement corrected for dust extinction (see previous section), m_J→J–2 is the modeled integrated intensity (Eq. (B.7)), σ_J→J–2 is the uncertainty (derived via Eq. (1)), and N = 8 is the number of data points. Lastly, p = 2 is the number of free parameters (the power-law index and total column density). The best-fit was then found by locating the minimum;χ² value in the array resulting from the chi-square test.

Fig. 7 Pixel-per-pixel comparison between the column densities derived from our best LTE and non-LTE thermal admixture models. The blue color map is a linear representation of the number of counts in the bins. The non-LTE model shown here was computed assuming a local density of n = 2.5 × 105 cm−3 (we explored a range: n = 104−107 cm−3).

5.3 Results

We performed a chi-square test based on Spitzer/IRS mapping observations of the S(0) to S(7) lines and their comparison with model grids, adopting the parameter boundaries Γ = 0,5-5 (linearly spaced) and N_tot(T ≥ T_min) = 10²⁰–10²⁴ cm^–2 (logarithmically spaced), and we explored a range of local densities $n_{{\text{H}}_2}\text{\hspace{0.17em}=\hspace{0.17em}}{10}^4-{10}^7{\text{cm}}^{-3}$ with our non-LTE approach. A direct comparison between our LTE and non-LTE results is shown in Fig. 7. The results of the non-LTE admixture model are presented in Figs. 6, 8, 9, and 10. All the pixels across the spectral maps have been reproduced within a confidence level of 3σ by the LTE and non-LTE models.

Our best-fit non-LTE model corresponds to a local density n = 2.5 × 10⁵ cm^–3 (see Fig. 6); hence, we fixed the input density to this value, which explains the small discrepancies between our LTE and non-LTE results (see Fig. 7). In Fig. 8, we present the chi-square grids resulting from the application of Eq. (8) toward the same locations as presented in Fig. 3 and which we initially analyzed in the previous section. The chi-square grids shown in the first column of the mosaic display the lowest ;χ² values because they correspond to low signal-to-noise regions of the map; hence, s_J→J–2/σ_J→J–2 is systematically lower and so is the product of Eq. (8). In all cases, the first black contour level (corresponding to $1.5{\chi }_{\min }^2,$, where ${\chi }_{\min }^2$ is the minimal value of ;χ² across the model grid) indicates the values of Γ and N_tot that best reproduce the observations. We found that column densities of the order 10²² cm^–2 and power-law indices Γ = 2.2–4.7 reproduce the observations, although with some degeneracy between the parameters. In IC443C, Shinn et al. (2011) found variations of Γ over the range 3-6, with the majority of sight lines in the range 4-5. In addition, they showed that spatially unresolved parabolic C-type bow shocks in steady state (Smith & Brand 1990) would produce a temperature stratification dN ∝ T^–3.8dT (i.e., Γ = 3.8). One would expect Γ to get closer to this predicted value within the shocked clump, but our measurements show the opposite: Γ decreases toward the shocked clump, suggesting a variation in the shock physical conditions. In fact, it is likely that there is a mixing of physical conditions along the line of sight, including both cold quiescent gas and warm high-velocity shocked gas. In addition, using different temperature boundaries (T_min, T_max) has an impact on the value of Γ, and Shinn et al. (2011) performed their analysis with T_min = 100 K, T_max = 4000 K, so it contributes to the discrepancies with our measurements of Γ in IC443G.

Based on our determinations of the best-fit parameters [Γ, N_tot(T ≥ T_min)] for each location of the map shown in Fig. 3, we built synthetic population diagrams (see Fig. 9) to compare our non-LTE results with the population diagrams presented in the previous section (Fig. 5). Toward these locations, the minimized;χ² varies between 0.85 and 2.13. The systematic deviation of the J_up = 5 level population above our best-fit model indicates that our power-law distribution might be unable to reproduce the break between the dynamic of the lower-J level populations (tracing cold molecular hydrogen) and the rest of the data points.

Additional tests and remarks. (i) We tried to modify the fixed parameter T_max from 1500 K to 2000 K, 3000 K and 4000 K to check how the best-fit is modified. These modifications increase the average ;χ², and they have little effect on the total H₂ column density (the relative difference is lower than 10%). Therefore, we kept our initial T_max. (ii) Similarly, we tried to adopt a lower value of T_min (down to 10 K), which increased the average ;χ² as well. (iii) We adopted a standard H₂ ortho-to-para ratio (o/p = 3) and checked that a two-fold variation in the ortho-to-para ratio has a negligible impact on the best-fit column density (less than 2%).

(i) Total column density map. The column density map obtained from the comparison of the S(0) to S(7) transitions with our model grid is shown in the left panel of Fig. 10. The column density varies between N ~ 10²¹ cm^–2 and N ~ 6 × 10²² cm^–2 between the outer edges and the core of the clump G. The peak column density is approximately 40 times higher than with the two-temperature population diagrams.

(ii) Average temperature map. The average temperature map obtained from the comparison of the S(0) to S(7) transitions with our model grid is shown in the right panel of Fig. 10. The average temperature map is produced by determining the contribution of each layer of molecular hydrogen to the temperature measured along the line of sight and dividing by the total column density: $$T_{\text{average }}={\displaystyle \sum _{T=T_{\min }}^{T_{\text{max }}}\frac{N_{\text{tot }}\left(T=T_{\text{layer }}\right)T_{\text{layer }}}{N_{\text{tot }}\left(T\ge T_{\min }\right)}}.$$(9)

The average temperature varies between T ~ 60 K and T ~ 140 K; hence, we succeeded in tracing colder gas. The spatial distribution of the temperature is uncorrelated to the column density map, and several localized ‘warm spots’ were found both within the clump and on the outer edges. We discuss the nature of these warm spots in Sect. 6.4.

(iii) H₂ temperature tomography. Since we had determined the best-fit parameters (Λ, Γ), Eq. (B.4) allowed us to build the partial column density of a layer of molecular hydrogen at any temperature T_min ≤ T_layer ≤ T_max. Using ΔT = 15 K and ΔT = 150 K temperature bins (respectively in the intervals 25–925 K and 925–1500 K), we built a partial column density map mosaic, shown in Fig. 11. This figure shows the variation in column density with respect to the temperature of distinct ther-malized layers of molecular hydrogen along the line of sight. We observed that the very high temperature molecular hydrogen (fourth row: T = 925–1500 K) is more clumpy than the warmer molecular hydrogen (second row: T = 325–475 K). The coldest molecular hydrogen (first row: T = 25–100 K) is also clumpy and more extended, and it has a higher column density. These maps show that the cold molecular hydrogen holds the bulk of the mass. The morphology of the warm gas (second and third rows: T = 325–925 K) is well correlated to the S(0) to S(2) lines, whereas the clumpy morphology of the hot gas (fourth row: T = 925–1500 K) is more correlated to the S(3) to S(7) lines (see Fig. 3 for comparison).

Fig. 8 Results of the chi-square test applied to our observations and the non-LTE model grid of the S(0) to S(7) H2 lines integrated intensities (see Eq. (8)). Each panel of the mosaic corresponds to the panels shown in Fig. 9 (positions indicated in Fig. 3). The term Γ is the power-law index, and Ntot is the total column density for T > Tmin = 25 K (see Eqs. (6) and (7)). The background color indicates the absolute value of χ2, and the set of black contours indicate the inner areas corresponding to the 1σ, 2σ, and 3σ levels, that is, ${\chi }^2\le \left[1.17{\chi }_{\text{min}}^2,1.67{\chi }_{\text{min}}^2,2.5{\chi }_{\text{min}}^2\right]$, where ${\chi }_{\min }^2$ is the minimum value of χ2 across the model grid $(\text{Δ}{\chi }^2(n\sigma )=[1+n/(N-p)]{\chi }_{\min }^2$ where σ is the measurement error, N is the number of measurements, and p is the number of fit parameters).

6 Discussion

We have estimated the total H2 column density and temperature maps from the analysis of the S(0) to S(7) pure rotational lines using three distinct models: a two-temperature description, an LTE admixture, and a non-LTE admixture (results are shown in Table 3, in which the uncertainty on T_min has been taken into account for the estimate of the column density $N_{{\text{H}}_2}$ and mass $M_{{\text{H}}_2}$). Comparison between our three models shows that with respect to the thermal admixture model, the ‘standard’ LTE analysis of the Spitzer/IRS H₂ lines might underestimate the total column density by ~97% in the region IC443G, which would be consistent with the conclusions of Roussel et al. (2007) in nearby galaxies (they found that between 1% and 30% of the total column density is traced by the allowed rotational transitions). Nonetheless, the true temperature distribution of the gas remains unknown, and the discrepancy between our models is mainly representative of the uncertainty associated with our assumptions.

Fig. 9 Population diagrams extracted from the positions shown in Fig. 3 and corresponding to the chi-square test shown in Fig. 8. The black data points represent the Nup measurements and their uncertainties. The red circles represent the best-fit non-LTE model obtained for the S(0) to S(7) transitions. The corresponding power-law indices and total column densities are shown in the bottom-left corner of the diagrams. The results of the χ2 test (see text) are presented in the top-right corner of the diagrams. The order of the diagrams (from left to right and top to bottom) is the same as the order of the white circles in Fig. 3.

6.1 Local density in IC443G

We aim to infer an estimate of the local density from the total H₂ column density, and we propose an approximate geometrical model for the shocked clump (for the sake of simplicity). We describe this structure by a cylinder of length l = 3′ and diameter D = 1′. At a distance of 1.8 kpc (Ambrocio-Cruz et al. 2017; Yu et al. 2019), this corresponds to 1–1.5 pc and D = 0.5 pc. Assuming that other sources of emission along the line of sight are negligible, the local density is given by the ratio of the measured column density and the length of shocked clump crossed. Hence, we have $n_{{\text{H}}_2}=N_{{\text{H}}_2}/D=6.0\times {10}^{22}/1.5\times {10}^{18}\simeq 4\times {10}^4 {\text{cm}}^{-3}$. The local density would be larger toward unresolved knots of H₂. In fact, if we assume that in any line of sight, all H₂ molecules are confined in a single quasi-spherical knot of diameter ~0.1 pc (~10″), then we would have $n_{{\text{H}}_2}\simeq 2\times {10}^5{\text{cm}}^{-3}$. Assuming T_kin = 1000 K, our first estimate ($n_{{\text{H}}_2}=4\times {10}^5{\text{cm}}^{-3}$) implies that the S(6) and S(7) lines might not be fully thermalized. The second estimate ($n_{{\text{H}}_2}~2\times {10}^5 {\text{cm}}^{-3}$) is consistent with the minimization of the local density in our non-LTE models (see Fig. 6), and it implies a thermalization of the S(0) to S(6) lines at T_kin = 1000 K.

Fig. 10 Resulting images from the LTE and SE analyses. Top panel: column density (left) and excitation temperature (right) derived from the analysis of the S(0) to S(7) H2 pure rotational lines, following the LTE method described in Sect. 5.2. Bottom panel: column density (left) and kinetic temperature (right) derived from the SE analysis described in Sect. 5.2.

6.2 Gas mass in IC443G

We estimated the H₂ molecular gas mass in IC443G by summing up all pixels in our column density map (see Fig. 10): $$M_{{\text{H}}_2}=2m_{\text{u}}{\displaystyle \sum _{i,j}{N}_{{\text{H}}_2}}[i,j]\text{Δ}s,$$(10)

where $N_{{\text{H}}_2}[i,j]$ is the column density of a pixel (i, j), Δs[cm²] ≃ 1.5 × 10¹⁰ ΔΩ[″]d_kpc is the physical (sky) surface of a pixel, and m_u ≃ 1.66 × 10^–27 kg is the atomic mass unit. Assuming d_kpc = 1.8 ± 0.2 (Ambrocio-Cruz et al. 2017; Yu et al. 2019) and taking into account the uncertainty on the cold gas temperature ($T_{\min }={25}_{-10}^{+15} \text{K}$), we obtained a mass $M_{{\text{H}}_2}={220}_{-80}^{+110}{\text{M}}_{\odot }$ (solar masses). This mass measurement does not correspond to the total mass of the shocked clump since the Spitzer/IRS field of observations does not cover the entire structure (see Fig. 2). We note, however, that our measurement is consistent with the shocked clump mass estimate inferred from our ¹²CO observations (performed in a 10′ × 10′ box), $M_{\text{CO}-\text{to}-{\text{H}}_2}=210\pm 50M_{\odot }$ (Dell’Ova et al. 2020).

6.3 H₂-to-¹²CO abundance

Using a Gaussian kernel, we re-convolved and re-sampled our H₂ column density map into the same angular resolution and spatial grid as our IRAM-30m and APEX ¹²CO column density maps (Fig. 2; Dell'Ova et al. 2020). A total of 54 pixels can be directly compared between the two column density maps, all characterized by high signal-to-noise ratios (s/n ≥ 6). The pixel-per-pixel comparison between the H₂ and ¹²CO column density maps is shown in Fig. 12. We found a quasi-linear relation between the column densities, and the H₂-to-¹²CO column density ratio estimated from the linear fit is 3.6 × 10⁴. The inferred ratio is in agreement with the standard value derived from visual extinction studies in molecular clouds, albeit a bit higher ($N_{{\text{H}}_2}/N_{{}^{12}\text{CO}}~{10}^4,$ Dickman 1978; Frerking et al. 1982). Toward the shocked clump, it is possible that ionization might decrease the ¹²CO abundance with respect to H2 via the conversion of ¹²CO into C⁺, C, CH₃OH, or CH₄, resulting in a departure from $N_{{\text{H}}_2}/N_{12}{}_{{}_{\text{CO}}}={10}^4.$ In order to prepare a future study of the chemistry in IC443G, we report the values of $N_{{\text{H}}_2}$ and $N_{{}^{12}\text{CO}}$ toward two positions using variable beam sizes in Tables 4 and 5.

6.4 Young stellar object candidates

In Fig. 13, we compare the spatial distribution of young stellar object candidates (YSOs) with the H₂ temperature map derived from our non-LTE analysis. This sample of YSO candidates was selected using color-color filtering of WISE and 2MASS infrared point sources (Dell’Ova et al. 2020). We used the Python package mistree to produce the minimal spanning trees corresponding to our samples of point sources (Naidoo 2019). The minimal spanning tree is defined as the network of lines, or branches, that connect a set of points together such that the total length of the branches is minimized and there are no closed loops (Cartwright & Whitworth 2004; Gutermuth et al. 2009). We found signs of clustering along the warm shocked clump and that a few point sources are correlated with hot knots of gas. We assumed that a few protostars could contribute to the localized heating shown in our map. In Fig. 13, there is only one warm spot that is not associated with a YSO candidate (to the eastern end of the Spitzer/IRS maps). We checked the XMM-Newton observations between 0.4 keV and 7.2 keV (Troja et al. 2006) for a local source of X-ray radiation in this area, and we found a source with a ~20″ shift with respect to the warm spot (see Fig. B.1). It is still possible that one or several embedded Class 0 protostar(s), which cannot be detected from color-color filtering of WISE and 2MASS infrared point sources, could be the cause of the heating toward this eastern position. It is important to note that associations between these point sources and the SNR IC443 can be coincidental. Additionally, considering the protostellar collapse phase timescale of approximately ~10⁵ yr (e.g., Lefloch & Lazareff 1994), we cannot conclude that the mid-infrared point sources are unequivocally linked to YSOs triggered by the IC443 SNR. In fact, these YSO candidates may belong to the same generation as the progenitor of the supernova, which may have evolved more rapidly due to its increased mass.

Fig. 11 Tomographic representation of $N_{{\text{H}}_2}\left(T\right)$ produced using the best-fit parameters (Λ, Γ). Each panel represents the partial column density map (defined by Eq. (B.4)) corresponding to a temperature bin ∆T (∆T = 75 K for the first three rows and ∆T = 150 K for the last row). The temperature T indicated in the top-right corner of each panel is the lower bound of the temperature interval on which the total column density is integrated (i.e., the upper bound is T + ∆T). Note: the color scale changes between rows.

Fig. 12 Pixel-per-pixel comparison of the 12CO and H2 total column density maps. The maps are restricted to the field of observations mapped by both Spitzer/IRS, IRAM-30m, and APEX telescopes. The 12CO data points were obtained with opacity-corrected population diagrams Dell’Ova et al. (2020). The H2 data points correspond to the results obtained with the thermal admixture SE description (see Sect. 5.2). We show the best linear fit associated with the sample of data points and the corresponding H2-to-12CO molecular abundance.

Fig. 13 Gas temperature map derived from our non-LTE thermal admixture model. The 2MASS and WISE young stellar object candidates are represented by a minimal spanning tree.

7 Summary

We performed an analysis of the Spitzer/IRS H₂ spectral-line maps toward the IC443G region. Our main findings are the following:

(a) The level populations traced by the H₂ S(0) to S(7) pure rotational lines cannot be reproduced by a single excitation temperature. This is as expected. When a mixture of shock velocities is involved, the shocked clump is stratified in temperature. Thus, using a power-law distribution of temperatures better represents the data, although it remains an approximation that cannot account for the complex mixing of physical conditions in IC443G, where the emission is driven by both cold quiescent gas and warm high-velocity shocked gas. In fact, even a single shock model may struggle to precisely replicate the observations, primarily because different shocks can be present within the observational beam, as discussed by Kristensen et al. (2023).

(b) Assuming that the temperature distribution of H₂ spans the range 25–1500 K, our thermal admixture model allows (within a confidence level of 3σ) the level populations to be successfully reproduced, with a power-law index Γ ~ 2.2–4.7, although with a substantial degeneracy between Γ and the total column density $N_{{\text{H}}_2}$. In IC443C, Shinn et al. (2011) showed that parabolic C-type bow shocks in a steady state can produce such power-law temperature stratification, as had been proposed for other shocked regions by Smith & Brand (1990). In this scenario, spatial variation of Γ across our field of observation would correspond to variations in the distributions of velocities at the head of unresolved bow shocks. Our average estimate of Γ is below the value predicted by shock models (Γ ~ 3.8). We note, however, that our degenerate model (see Fig. 8) is consistent with Γ ~ 3.8. More importantly, Shinn et al. (2011) used different temperature boundaries to predict this value of Γ.

(c) The total column density of H₂ inferred from our two-.temperature excitation diagram accounts for approximately 3% of the total column density deduced from our thermal admixture model. Our analysis shows that an estimate of the cold molecular gas temperature is paramount to constraining the total H₂ mass. The ¹²CO rotational transitions, thanks to their low-temperature energy levels (E₁ = 5.53 K, E₂ = 16.6 K, E₃ = 33.2 K), can provide such a temperature estimate. We acknowledge, however, that small inaccuracies in the assumed boundaries and shape of the temperature distribution may result in significant errors for the derived column densities.

(d) The agreement between the observations and our thermal admixture model suggests that the S(0) to S(7) H₂ lines could be thermalized. However, the local density inferred from our column density measurement is in the range 0.2−2 × 10⁵ cm⁻³, depending on the adopted geometry. Our χ² -minimization also suggests a local density of n = 2.5 × 10⁵ cm⁻³. This range of. local densities is not consistent with a full thermalization of the S(6) and S(7) lines.

(e) Taking into account the uncertainty on the minimal temperature of the gas ($T_{\min }={25}_{-10}^{+15}$ K, estimated from ¹²CO observations), we measured a molecular gas mass $M_{{\text{H}}_2}={220}_{-80}^{+110}$ M_⊙ across the Spitzer/IRS field of observations. This measurement represents a lower bound of the mass of protons available for CRs to interact with in IC443G since there is also atomic gas along the line of sight (Lee et al. 2008).

(f) We compared our findings with the analysis of ¹²CO rotational lines toward IC443G (Dell’Ova et al. 2020). Our mea. surements are consistent with a H₂-to-¹²CO abundance ratio of 3.6 × 10⁴, indicating that some of the carbon may exist in the form of C or C⁺ in this region.

The analysis of the Spitzer/IRS H₂ spectral-line maps toward the IC443G region provides an insight into the complex dynamics of shocked molecular gas. The observed temperature distribution, reflecting both quiescent and shock-excited components, emphasizes the need for shock models to capture these dynamics fully. Additionally, our findings stress the importance of accounting for multi-temperature environments when attempting to measure the gas mass. Future work should focus on comprehensive shock modeling and higher resolution observations, to disentangle the emission of protostellar outflows from the large-scale shock structure.

Acknowledgements

This project has received funding from the French Agence Nationale de la Recherche (ANR) through the project COSMHIC (ANR-20-CE31-0009), and the French Programme National de Physique et Chimie du Milieu Interstellaire (PCMI) of CNRS/INSU (with INC/INP/IN2P3).

Appendix A Description of the pre-processing transformations: Systematic noise reduction

Using Python, we applied the following pre-processing steps to. the Spitzer/IRS spectral-line maps:

• 1.

  The S(0) spectral-line map is strongly affected by the continuum emission of two stars at α_[J2000] = 6^h16^m41^s, δ_[J2000] = 22º30′58″ and ${\alpha }_{[\text{J}2000]}=6^{\text{h}}{16}^{\text{m}}{44}^{\text{s}}.5$, δ_[J2000] = 22º31′33″. We used two 2D Gaussian functions (G₁, G₂, re-normalized between zero and one) and the S(1) map as an emission model to remove these components in the pixels (x, y): $$\begin{array}{cc}& S^0(x,y)=S^0(x,y)\left[1-\left(G_1(x,y)+G_2(x,y)\right)\right]\\ & +M\times S^1(x,y)\left[G_1(x,y)+G_2(x,y)\right],\end{array}$$(A.1)

  where S⁰ and S¹ are respectively the S(0) and S(1) maps, and M = median(S⁰/S¹) is a re-scaling factor.

• 2.

  The S(6) spectral-line map has the most severe horizontal stripes (see Fig. A.1). We used the S(5) map as an emission model to isolate the emission from the stripes and generate a noise model N(x, y) from a data array initially filled with. zeros: $$\alpha (x,y)=\left[S^6(x,y)/\max \left(S^6\right)-S^5(x,y)/\max \left(S^5\right)\right]$$(A.2) $$N(x,y)=S^6(x,y)\text{\hspace{0.17em}where\hspace{0.17em}}\alpha (x,y)>{\alpha }_{\text{treshold}},$$(A.3)

  where S⁶ and S⁵ are respectively the S(6) and S(5) maps, and α_treshold is an arbitrary threshold that controls the separation between the signal and the systematic stripes. The pre-processed S(6) map is then given by S⁶(x,y) = S ⁶(x,y) −N (x,y).

• 3.

  All the SL maps (i.e., S(2) to S(7)) were systematically presenting a strong vertical stripe. This stripe was located on the columns 52 and 53 of the data arrays. To correct this artifact, we performed the following interpolation: $$S^i(52,y)=\frac{1}{2}\left(S^i(51,y)+S^i(54,y)\right),i=0,\dots ,7$$(A.4) $$S^i(53,y)=\frac{1}{2}\left(S^i(51,y)+S^i(54,y)\right),i=0,\dots ,7.$$(A.5)

• 4.

  Using Gaussian kernels, we convolved the S(1) to S(7) maps to the resolution of the S(0) map $\left({8^{″}}_{·}3\right)$. This convolution lowers the contribution from the horizontal stripes in all maps..

Fig. A.1 H2 v = 0 − 0 S(0) to S(7) raw maps obtained with Spitzer/IRS toward the region IC443G (Neufeld et al. 2004). The S(0) line was mapped by the LH module, the S(1) and S(2) lines were mapped by the SH modules, and the rest by the SL module. The S(2) line was mapped by both the SH (top row) and SL (second row) modules. The coordinates of the field of observations are indicated in Fig. 2. Pre-processed maps are shown in Fig. 3

The results of these pre-processing steps are presented in Fig. 3. The systematic vertical stripe as well as the horizontal stripes in the S(6) map have both been removed. The average flux variation across the map between the raw and pre-processed images is ~1.02%; hence it should not significantly modify the average results, and it is lower than the relative difference between the raw S(2)[SH] and S(2)[SL] maps (~5.5%).

Appendix B Thermal admixture model (LTE)

We built model grids of the integrated intensity for the S(0) to S(7) lines based on the H₂ thermal admixture description. The H₂ total column density at a temperature T = T_layer (referred to as N_tot(T = T_layer)) can be determined from Eq. 6 and Eq. 7. If we adopt the same LTE assumptions as in Section 4.1 for each layer of molecular hydrogen at a temperature T_layer, then the upper-level populations can be derived from the Boltzmann distribution for N_tot(T = T_layer): $$N_{\text{up}}(T)=\frac{N_{\text{tot}}}{Z(T)}{g}_{\text{up}}{\text{e}}^{-E_{\text{up}}/k_{\text{B}}T},$$(B.1)

where 𝑔_up is the degeneracy, E_up is the energy of the level, k_B is the Boltzmann constant and T is the temperature of the gas, and Z(T) is the partition function defined by: $$Z(T)={\displaystyle \sum _s}{g}_s{\text{e}}^{-E_s/k_{\text{B}}T}.$$(B.2)

The integrated intensity produced by the layer of molecular hydrogen can then be inferred from the following relation: $${{\displaystyle \sum }}^{\text{}}{I}_v\text{Δ}v=\frac{h{v}_{\text{ul}}}{4\pi }{N}_{\text{up}}{A}_{\text{ul}},$$(B.3)

where h is the Planck constant, ν_ul is the frequency of the spectral line, and A_ul is the Einstein coefficient for spontaneous. de-excitation. The sum of the integrated intensity contributions over the column density distribution yields an estimate of the integrated intensity measured along the line of sight.

Appendix B.1 Discretization of the column density distribution N(T)

The discretization of Eq. 6 yields $$\text{Δ}N=\text{Λ}{T}^{-\text{Γ}}\text{Δ}T,$$(B.4)

where ∆N = N_tot(T = T_layer) is the column density of a layer of molecular hydrogen at temperature T and ∆T = (T_max − T_min)/N_bin is a temperature bin (N_bin is the total number of bins). We produced a discretized column density distribution that satisfies the following condition: $$N_{\text{tot}}\left(T\ge T_{\min }\right)={\displaystyle \sum _{T=T_{\min }}^{T_{\max }}}\text{Λ}{T}^{-\text{Γ}}\text{Δ}T.$$(B.5)

Appendix B.2 I(Λ, Γ) model grid

For a choice of parameter (Γ, N_tot(T ≥ T_min)), we produced a model of the S(0) to S(7) lines integrated intensity by performing the following computations:

• 1.

  Assuming LTE, we computed the upper-level populations of each molecular hydrogen layer, applying Eq. B.1: $$N_{\text{up}}\left(T=T_{\text{layer}}\right)=\frac{N_{\text{tot}}\left(T=T_{\text{layer}}\right)}{Z\left(T_{\text{layer}}\right)}{g}_{\text{up}}{\text{e}}^{-E_{\text{up}}/k_{\text{B}}{T}_{\text{layer}}},$$(B.6)

  where Z(T) is the partition function defined by Eq. B.2.

• 2.

  Then, we computed the integrated intensity of the S(0) to S(7) lines produced by each layer of gas, applying Eq. B.3: $${{\displaystyle \sum }}^{\text{}}{I}_{\text{ul}}\text{Δ}v=\frac{h{v}_{\text{ul}}}{4\pi }{N}_{\text{up}}\left(T=T_{\text{layer}}\right)A_{\text{ul}}.$$(B.7)

Fig. B.1 Gas temperature map derived from our non-LTE thermal admixture model. White contours represent hard X-ray emission (2–7.2 keV) detected by XMM-Newton (Troja et al. 2006)).

References

All Tables

All Figures

Fig. 1 Overview of the SNR IC443. Top-left panel: schematic drawing of the SNR morphology. The two-shell model proposed by Troja et al. (2006) is shown as a colored background image: the gray torus represents the molecular cloud discovered by Cornett et al. (1977) and is highlighted by white curved lines based on the higher resolution observations by Lee et al. (2012). The blue hemisphere labeled “subshell A” represents the shock front in the eastern region. In the northeast, where it has been confined by the encounter with the neutral H I cloud of Denoyer (1978), the ionic shock front is traced by optical, infrared and very soft X-ray emission (Alarie & Drissen 2019, dashed-dotted black lines). In the southern region, the thick black curved lines represent the areas of interaction with the molecular cloud and the shocked H2 ridges of Burton et al. (1988). The labels and positions of the different molecular clumps identified by Denoyer (1979b) and Dickman et al. (1992) are indicated with capital letters (B, C, D, E, F, G, H). The red hemisphere labeled “subshell B” represents the shock front in the western region, where it is expanding in a homogeneous and less dense medium. The position of the PWN discovered by Olbert et al. (2001) is indicated by a white star. The extension of the radio continuum halo is shown by dashed black lines (Lee et al. 2008, later confirmed by Castelletti et al. 2011). Locations and extensions of the four gamma-ray peaks at different energy ranges are indicated: EGRET centroid (∘), MAGIC centroid (□), VERITAS centroid (◊), and Fermi LAT centroid (∆), respectively measured by Esposito et al. (1996), Albert et al. (2007), Acciari et al. (2009) and Abdo et al. (2010). The respective localization errors are shown as crosses of size 0.5 σ. We note that the high-energy emission is extended, and hence the relation between Fermi LAT, VERITAS, and MAGIC observations may be better characterized by the extent of the emission rather than the centroids. The sizes of typical instrumental beams are indicated in the bottom-right corner of the figure. The other panels represent color composite images of IC443, with H I 21 cm line emission (Lee et al. 2008) in yellow and 24 µm (Spitɀer/MIPS) emission in white in all three panels. The third color channel represents a distinct signal in each panel: top-right panel: 12CO J = 1–0 (Lee et al. 2012, green); bottom-left panel: DSS optical observations (York et al. 2000, blue); bottom-right panel: XMM-Newton observations at 1.4–5 keV (Greco et al. 2018, purple).

In the text

	Fig. 2 Finderchart of the Spitɀer/IRS field of observations. Left: 12CO J = 2−1 10′ × 10′ map of the peak temperature obtained with the IRAM-30m telescope at an angular resolution of 11″.2 (Dell’Ova et al. 2020). Right: H2 υ = 0–0 S(1) map of the integrated line intensity obtained with Spitɀer/IRS.
In the text

Fig. 3 H2 υ = 0–0 S(0) to S(7) pre-processed maps obtained with Spitɀer/IRS toward the region IC443G. The S(1) to S(7) maps were convolved to the spatial resolution of the S(0) map (8″.3, see Sect. 4.2). The S(0) line was mapped by the LH module, the S(1) and S(2) lines were mapped by the SH modules, and the rest were mapped by the SL module. The S(2) line was mapped by both the SH (top row) and SL (second row) modules. The coordinates of the field of observations are indicated in Fig. 2. Raw maps are shown in Fig. A.1. The array of white circles indicates the positions in which we extracted the signals of the S(0) to S(7) spectral-line maps to produce the population diagram mosaics presented in Fig. 5 (model 1: standard LTE approach) and Fig. 9 (model 2: H2 thermal admixture model). The size of the circles corresponds to the angular resolution of the S(0) observations, that is, 8″.3 in diameter.

In the text

	Fig. 4 Pixel-per-pixel comparison of the integrated intensity S(2) maps measured by the Spitɀer/IRS SH and SL modules before (left panel) and after pre-processing (rightpanel). The solid lines represents the 1:1 relation expected if the SH and SL data arrays were identical, and the dashed lines represent the ±25% relative errors.
In the text

Fig. 5 Population diagrams extracted from the positions shown in Fig. 3. The corresponding positions are indicated in the top-left corner of each panel. The black data points represent the Nup measurements and their uncertainties. The gray curve represents the linear model obtained for the S(3) to S(7) transitions. The red curve represents the linear model obtained for the S(0) to S(2) transitions. Gray circles represent the sum of the two models for each transition. The errors associated with the best linear models are represented by a filled area around the curves. The corresponding excitation temperatures are shown in the bottom-left corner of the diagrams ('cold' component in red and 'warm' component in gray). The total column density is also shown (in black). The results of the χ2 test (see text and Eq. (5)) are presented in the top-right corner of the diagrams. The order of the diagrams (from left to right and top to bottom) is the same as the order of the white circles in Fig. 3.

In the text

	Fig. 6 Average chi-square (‘χ2’, i.e., the average chi-square across the map, indicated by red diamonds) and average column density ($\overline{‘N_{H_2}’}$, indicated by turquoise squares) with respect to the adopted local density. The best-fit model was found with a density n = 2.5 × 105 cm−3, corresponding to χ2 = 1.13 and $\overline{N_{{\text{H}}_2}}=1.64\times {10}^{22} {\text{cm}}^{-2}$
In the text

	Fig. 7 Pixel-per-pixel comparison between the column densities derived from our best LTE and non-LTE thermal admixture models. The blue color map is a linear representation of the number of counts in the bins. The non-LTE model shown here was computed assuming a local density of n = 2.5 × 105 cm−3 (we explored a range: n = 104−107 cm−3).
In the text

Fig. 8 Results of the chi-square test applied to our observations and the non-LTE model grid of the S(0) to S(7) H2 lines integrated intensities (see Eq. (8)). Each panel of the mosaic corresponds to the panels shown in Fig. 9 (positions indicated in Fig. 3). The term Γ is the power-law index, and Ntot is the total column density for T > Tmin = 25 K (see Eqs. (6) and (7)). The background color indicates the absolute value of χ2, and the set of black contours indicate the inner areas corresponding to the 1σ, 2σ, and 3σ levels, that is, ${\chi }^2\le \left[1.17{\chi }_{\text{min}}^2,1.67{\chi }_{\text{min}}^2,2.5{\chi }_{\text{min}}^2\right]$, where ${\chi }_{\min }^2$ is the minimum value of χ2 across the model grid $(\text{Δ}{\chi }^2(n\sigma )=[1+n/(N-p)]{\chi }_{\min }^2$ where σ is the measurement error, N is the number of measurements, and p is the number of fit parameters).

In the text

Fig. 9 Population diagrams extracted from the positions shown in Fig. 3 and corresponding to the chi-square test shown in Fig. 8. The black data points represent the Nup measurements and their uncertainties. The red circles represent the best-fit non-LTE model obtained for the S(0) to S(7) transitions. The corresponding power-law indices and total column densities are shown in the bottom-left corner of the diagrams. The results of the χ2 test (see text) are presented in the top-right corner of the diagrams. The order of the diagrams (from left to right and top to bottom) is the same as the order of the white circles in Fig. 3.

In the text

	Fig. 10 Resulting images from the LTE and SE analyses. Top panel: column density (left) and excitation temperature (right) derived from the analysis of the S(0) to S(7) H2 pure rotational lines, following the LTE method described in Sect. 5.2. Bottom panel: column density (left) and kinetic temperature (right) derived from the SE analysis described in Sect. 5.2.
In the text

Fig. 11 Tomographic representation of $N_{{\text{H}}_2}\left(T\right)$ produced using the best-fit parameters (Λ, Γ). Each panel represents the partial column density map (defined by Eq. (B.4)) corresponding to a temperature bin ∆T (∆T = 75 K for the first three rows and ∆T = 150 K for the last row). The temperature T indicated in the top-right corner of each panel is the lower bound of the temperature interval on which the total column density is integrated (i.e., the upper bound is T + ∆T). Note: the color scale changes between rows.

In the text

Fig. 12 Pixel-per-pixel comparison of the 12CO and H2 total column density maps. The maps are restricted to the field of observations mapped by both Spitzer/IRS, IRAM-30m, and APEX telescopes. The 12CO data points were obtained with opacity-corrected population diagrams Dell’Ova et al. (2020). The H2 data points correspond to the results obtained with the thermal admixture SE description (see Sect. 5.2). We show the best linear fit associated with the sample of data points and the corresponding H2-to-12CO molecular abundance.

In the text

	Fig. 13 Gas temperature map derived from our non-LTE thermal admixture model. The 2MASS and WISE young stellar object candidates are represented by a minimal spanning tree.
In the text

Fig. A.1 H2 v = 0 − 0 S(0) to S(7) raw maps obtained with Spitzer/IRS toward the region IC443G (Neufeld et al. 2004). The S(0) line was mapped by the LH module, the S(1) and S(2) lines were mapped by the SH modules, and the rest by the SL module. The S(2) line was mapped by both the SH (top row) and SL (second row) modules. The coordinates of the field of observations are indicated in Fig. 2. Pre-processed maps are shown in Fig. 3

In the text

	Fig. B.1 Gas temperature map derived from our non-LTE thermal admixture model. White contours represent hard X-ray emission (2–7.2 keV) detected by XMM-Newton (Troja et al. 2006)).
In the text