Interstellar anatomy of the TeV gamma-ray peak in the IC443 supernova remnant

Supernovae remnants (SNRs) represent a major feedback source from stars on the interstellar medium of galaxies. During the latest stage of supernovae explosions, shock waves produced by the initial blast modify the chemistry of gas and dust, inject kinetic energy in the surroundings, and may alter star formation characteristics. Simultaneously, gamma-ray emission is generated by the interaction between the ambiant medium and the cosmic rays. We study the stellar and interstellar contents of IC443, an evolved shell type SNR at a distance of 1.9 kpc, with an estimated age of 30 kyr. We aim to measure the mass of the gas within the extended G region, which corresponds to the peak of gamma-ray emission detected by VERITAS and Fermi. We performed 10'x10' mapped observations of 12CO and 13CO J=1-0, J=2-1 and J=3-2 pure rotational lines, as well as C18O J=1-0 and J=2-1 obtained with the IRAM-30m and APEX telescopes. We first compared our data with local thermodynamic equilbrium (LTE) models. We estimated the optical depth of each line from the emission of the isotopologues 13CO and C18O. We used the population diagram and large velocity gradient (LVG) assumption to measure the column density, mass, and kinetic temperature of the gas using 12CO and 13CO lines. We used complementary data (stars, gas, and dust at multiple wavelengths) and infrared point source catalogues to search for protostar candidates. Our results emphasize how the mass associated with the ring-like structure and the cloudlet cannot be overlooked when quantifying the interaction of cosmic rays with the dense local medium. Additionally, the presence of numerous possible protostars in the region might represent a fresh source of CR, which must also be taken into account in the interpretation of gamma-ray observations in this region.


Introduction
The violent end of some stellar objects, a supernova (SN) explosion, is the beginning of an incredible sequence of energy injection in the surrounding interstellar medium (ISM). A SN explosion ejects material with mass ranging from ∼1.4 to ∼20 M and a typical energy of 10 50−52 ergs (see e.g. Draine 2011) and drives fast shock waves, about 10 4 km s −1 ahead of the ejecta (the expelled stellar material) through the interstellar medium. Already at early stages, SNe play a crucial, multifaceted role in the evolution of galaxies. The explosion and subsequent dispersion of matter is itself the most important source of elements heavier than nitrogen in the gas phase (François et al. 2004). The fast shocks inject kinetic energy that gradually decays in Article number, page 1 of 36 A&A proofs: manuscript no. 3473992 turbulence, a key process in the regulation of star formation on galactic scales (Mac Low & Klessen 2004). The fast shock fronts are also a recognized site for the production of the bulk of cosmic rays (CRs, at least at GeV energies; see e.g. Bykov et al. 2018a;Tatischeff & Gabici 2018 and references therein for recent reviews), although alternative origins draw more and more attention (such as superbubbles and Fermi bubbles, see e.g. Grenier et al. 2015;Gabici et al. 2019 and references therein for recent reviews). Finally, the hot (10 6 − 10 8 K) and relatively dense (1 − 10 cm −3 ) conditions in the ejecta could be favourable to the synthesis of cosmic dust up to hundred years after the explosion (see recent reviews by Cherchneff 2014;Sarangi et al. 2018;Micelotta et al. 2018). After the first phases of expansion (free, then adiabatic), the temperature of the shock front drops to 10 6 K, allowing the gas to radiatively cool down. At this stage, the so-called supernova remnant (SNR) resembles a spherical shell of 10 − 20 pc radius, delimitated by regions where shocks interact with the ambient medium.
This evolved SNR stage is also of fundamental importance to many aspects of galactic evolution. The fastest remaining shocks (between ∼30 and a few hundreds of km s −1 ) dissociate and ionize the pre-and post-shock medium, and generate far-ultraviolet photons (e.g. Hollenbach & McKee 1989b). More generally, shocks heat, accelerate, and compress the ambient medium. They inject energy and trigger specific dust and gasphase chemical processes (e.g. van Dishoeck et al. 1993), hence significantly participating to the cycle of matter in galaxies. CRs accelerated in the earlier stages of the explosion and trapped in the shock fronts now interact with the dense medium, producing observable X-ray to γ-ray photons (Gabici et al. 2009a;Celli et al. 2019;Tang 2019 and references therein). Cosmic rays of Galactic origin can also be re-accelerated in the shock regions (like in W44, e.g. Cardillo et al. 2016). Finally, evolved SNRs play a key role on star formation. Like all massive stars, the progenitor of a SN explosion has formed in a cluster, in a dense and inhomogeneous environment, where lower-mass stellar companions have a greater life expectancy (Montmerle 1979a). During its life on the main sequence, the progenitor has driven stellar winds in the surrounding medium, possibly triggering a second generation of star formation in the neighbouring molecular clouds (e.g. Koo et al. 2008a). The compression and cooling caused by SNR shocks might also generate star formation (eg. Herbst & Assousa 1977). In any case, the injection of energy exerted by SNRs in all possible forms (CRs, energetic photons, shocks) likely locally alters the characteristics of all possible star formation events over significant spatial scales and times.
With the present study we aim to start characterizing as precisely as possible and on fields as large as possible the mechanisms of energy injection (shocks, photons, CRs) exerted by an evolved SNR, and its consequences on the local star formation. Such a study must be performed on an evolved object, since the energy injection effects can spread over the full duration of the SNR phase, and are all the more visible when the SNR is old. In particular, we want to provide support for the study of CRs properties (acceleration, composition, diffusion) in evolved SNRs. In these objects, CRs interact with the local medium through four processes that all generate γ-ray photons: pion decay from the collision of hadronic CRs with the dense material, Bremsstrahlung from the interaction of leptonic CRs with the local dense medium, inverse Compton scattering of leptonic CRs with the local radiation field, and synchrotron emission of leptonic CRs gyrating around the local magnetic fields. Our first aim is constraining the properties of the local medium that is the target of these interactions: mass and density of all ob-served components, magnetic field structure, and radiation field structure.
Our second aim is to identify all possible sources of on-going acceleration of 'fresh' CRs additional to the 'old' injection of CRs previously accelerated by the SNR. These sources can be of two kinds: ionized regions where kinetic energy is deposited and where the magnetic field structure and the ionization fraction make the acceleration possible (like [H ii] regions, see Padovani et al. 2019), or protostellar jets and outflows (where these conditions can be naturally combined). Our work can thus provide support for further studies of CR-related questions only if the study of local star formation is performed simultaneously. Indeed, Padovani et al. (2015Padovani et al. ( , 2016 have shown that jets can accelerate low-energy CRs, that can be re-accelerated in the shock fronts of the remnant. Other studies have confirmed that supermassive star clusters neighbouring SNRs can be a source of CRs (Hanabata et al. 2014). Conversely an optimal characterization of CR action on the local formation is mandatory to better understand star formation in older galaxies, up to z ∼ 2 for which the star formation efficiency is higher (Madau & Dickinson 2014). Indeed, such a star formation regime is reminiscent of starburst galaxies, where SNe from a given generation of stars affect the next one.
The threefold and intertwined goals of our study (interstellar medium, star formation and CR) make the study of large fields mandatory. Indeed, CR studies rely on γ-ray spectra obtained with telescopes that only provide a limited spatial resolution, typically a few arcminutes. This extent is the minimum field where we have to characterize the interstellar medium and the star formation as best as we can. This is why we have chosen to study a 10 ×10 field in the relatively evolved IC443 SNR (see Fig. 1). More particularly, with the present paper we investigate the physical conditions and dynamical structure of the molecular gas and its association with protostars in such a field located at the peak of γ-ray emission detected in the IC443 SNR, based on observations of the CO emission lines and its isotopologues. In Sect. 2 we present a summarized review of the source. In Sect. 3 we present our observations and propose a description of the morphology and kinematics of the region, emphasising three distinct molecular components. Sect. 4 focuses on the measure of the gas mass for these components. First we perform a local thermodynamical equilibrium (LTE) analysis of the 12 CO, 13 CO and C 18 O emission lines and we build pixel-per-pixel, channel-perchannel population diagrams corrected for optical depth. Then we propose a second method using a radiative transfer code based on the Large Velocity Gradient (LVG) approximation. In Sect. 5 we study the spectral energy distribution of point sources identified in infrared survey catalogs, as well as the spatial distribution of optical point sources detected with GAIA. Finally we summarize our findings in Sect. 6.  (Claussen et al. 1997a). The white dot marks the position of the pulsar wind nebula (Olbert et al. 2001b). The red box represents one of the fields observed by Spitzer-IRS (Neufeld et al. 2004), corresponding roughly to the G region defined by Huang et al. (1986). The white box represents the 10 ×10 field of our observations, that we name 'the extended G region'. The IRAM-30m instrumental beam diameter corresponding to the 12 CO(1-0) transition and the size of the typical PSF of γ-ray telescopes (5 ) are indicated by white disks in the lower left corner of the figure.

The supernova remnant IC443
IC443 is a mixed-morphology supernova remnant, located at a distance of 1.5-2 kpc (Denoyer 1978, Welsh & Sallmen 2003, with recent measures suggesting a kinematic distance of 1.9 kpc (Ambrocio-Cruz et al. 2017). IC443 is an evolved SNR, yet its exact age is a matter of debate. The literature contains two kinds of value (∼3 and ∼30 kyr), depending on the type of data that is analysed. A compelling finding concerning the age and origin of the IC443 SNR was the discovery of the CXOU J061705.3+222127 Pulsar Wind Nebula (PWN), based on Chandra X-ray observatory and Very Large Array (VLA) images (Olbert et al. 2001a). The motion of the PWN is consistent with an age of 30 kyr for the SN event, and its detection strongly supports a core-collapse formation scenario for the SNR. IC443 displays a shell morphology in radio, with two atomic sub-shells (shells A and B, Braun & Strom 1986). It is one of the most striking example of a SNR interacting with neighbouring molecular clouds. The most up-to-date and complete description of the structure and kinematics of both the atomic and dense molecular environment of the SNR was offered by Lee et al. (2008Lee et al. ( , 2012. Their ∼ 1 • ×1 • map of the J=1-0 transition of 12 CO allowed to characterize both the incomplete molecular shell interacting with shocks (toward the southern part of the SNR) and the molecular cloud that is associated with the remnant. Continuum radio emission in IC443 is partly correlated with the molecular shell and the secondary [H i] shell (see Castelletti et al. 2011 for a detailed description of the low-frequency radio emission in IC443 at 74 and 330 MHz, and Egron et al. 2017, Loru et al. 2019 for high and very-high frequency studies, respectively at 7 GHz and 21.4 GHz).
CO emission was observed by Denoyer (1979a) towards the SNR, revealing three shocked CO clumps along the southern molecular ridge (labelled A, B and C). Follow-up observations of CO J=1-0 over a 50 ×50 field by Huang et al. (1986) allowed to detect new areas of shock-cloud interaction and to identify 5 previously unknown CO clumps, extending the classification started by Denoyer (1979a) and providing the first mention of the G knot. The OH 1720 MHz line is a powerful diagnostic for the classification of SNRs interacting with molecular clouds (Frail et al. 1996). 6 masing spots have been identified in IC443 Table 1: Spectroscopic parameters corresponding to the observed lines. J up is the rotational quantum number, ν i j and A i j are respectively the frequency and the Einsten coefficient of the pure rotational transition ( = 0 where v is the vibrational quantum number). E up and g up are respectively the energy and degeneracy corresponding to the upper level. The values given are taken from the Cologne Database for Molecular Spectroscopy (Müller et al. 2001, Müller et al. 2005, Endres et al. 2016 and Jet Propulsion Laboratory database (Pickett et al. 1998), and are the numeric values used for all our measures.  Claussen et al. (1997b), all located in the G region delineated by Huang et al. (1986). He proposed that OH masers spots could be promising candidates for the sites of CR acceleration. Lockett et al. (1999) improved the modelling of the shock origin for these maser lines, associated with moderate temperatures (50-125 K), local densities (∼ 10 5 cm −3 ) and OH column densities of the order of 10 16 cm −2 , then Wardle (1999) added the effect of the dissociation of molecules by far ultraviolet (FUV) photons in molecular clouds subject to CR and X-ray ionization (see Hoffman et al. 2003, Hewitt et al. 2006, Hewitt et al. 2008, Hewitt et al. 2009 for recent studies). From J=1-0 12 CO and HCO + emission, Dickman et al. (1992) measured a mass of 41.6 M for clump G, and estimated that a total molecular mass of 500-2000 M is interacting with the SNR shocks, which corresponds to 5%-10% of the SN energy considering that the average velocity of the clumps is 25 km s −1 . Zhang et al. (2010) showed that two distinct structures are resolved in the region G, labelling G1 the strongest 13 CO peak and G2 the previously mentionned shocked clump. Lee et al. (2012) measured a mass of 57.7±0.9 M for the clump G2. Oddly, Xu et al. (2011) measured a mass of 2.06×10 3 M for the 'cloud G', which is much higher than the previous estimates. Several molecular shocks were mapped within the SNR using 12 CO lines (e.g. White et al. 1987, Wang & Scoville 1992 for clumps A, B and C). In particular, the kinematics of clump G were characterized in details by van Dishoeck et al. (1993) who presented observations of the rotational transition J=3-2 of CO along the shocked molecular ring at a spatial resolution of 20 -30 . In the last ten years, the large scale molecular contents of IC443 have been scrutinized with increasing precision and completeness, since several authors have mapped the J=1-0 transition of the isotopologues 12 CO, 13 CO and C 18 over large fields (from 40 ×45 to 1.5 • ×1.5 • , Zhang et al. 2010, Lee et al. 2012, Su et al. 2014. Several sub-millimeter observations of IC443 were performed to characterize the shocked molecular gas (e.g. van Dishoeck et al. 1993 for a study of the shock chemistry in the southern ridge, towards the clumps B, C and G). Notably, the ground state of shocked ortho-H 2 O was detected towards the clumps B, C and G (Snell et al. 2005). The incomplete shell-morphology is also observed in the J, H, K bands observed by 2MASS (Two-Micron All-Sky Survey, Rho et al. 2001), as well as in infrared and far infrared observa-tions by Spitzer-MIPS (Pinheiro Gonçalves et al. 2011) and WISE (Wide-field Infrared Survey Explorer, Wright et al. 2010). Excitation by shocks was suggested as the most likely scenario for the emission lines detected in IR, instead of X-ray or FUV mechanisms (e.g. [O i], Burton et al. 1990). Rho et al. (2001) analysed the emission of [O i] with shock models, suggesting a fast J-type shock (∼100 km s −1 ) in the NE atomic shell, and a C shock (v s =∼30 km s −1 ) propagating in the southern molecular ridge. In the effort to study molecular shocks, H 2 pure rotational transitions were mapped towards the clump G by ISOCAM (ISO, Cesarsky et al. 1999) and compared to non-stationary shock models, as well as towards the clumps C and G by Spitzer-IRS (Infrared Spectrograph) (Neufeld et al. 2007). The molecular clumps B, C and G were also observed by AKARI (Shinn et al. 2011) and by the The Stratospheric Observatory for Infrared Astronomy (SOFIA) (Reach et al. 2019). All these studies were carried out in small fields (∼ 1 ×1 ) and allowed to put constraints on the shock velocity (∼ 30 km s −1 ) and pre-shock density (∼ 10 4 cm −3 ), and to outline similarities with protostellar shocks in the southern ridge where we aim to focus on the extended G region.
The optical emission is well correlated with radio and [H i] features, reproducing shells A and B. In particular, the SNR displays bright, filamentary structures towards the northeastern part of the remnant (Fesen & Kirshner 1980, Alarie & Drissen 2019 . IC443 was fully mapped by the Sloan Digital Sky Survey (SDSS, York et al. 2000). There are no bright features towards the extended G region, but optical studies offered constraints on the global characteristics of IC443. Ambrocio-Cruz et al. (2017) estimated an age of ∼30 kyr and an energy of 7.2×10 51 erg injected in the environment by the SNR from the comparison of observations of [Hα] with SNR models (Chevalier 1974).
The shell-like structure of IC443 in radio, centrally filled in X-rays puts the remnant into the category of mixed-morphology SNRs (Petre et al. 1988, Rho & Petre 1998. Observations of the hard X-ray contents of IC443 (up to 100 keV) by Bep-poSAX show hints of shock-cloud interaction (Bocchino & Bykov 2000). XMM-Newton mapped IC443 in the ranges 0.3-0.5 keV and 1.4-5.0 keV with an unprecedented field of view and spatial resolution (Bocchino & Bykov 2003), and showed that the soft X-ray emission is partly absorbed by the nearby molecular cloud (Troja et al. 2006). Troja et al. (2008) reported the detection of a ring-shaped ejecta encircling the PWN, associated with a hot metal rich plasma which abundances are in agreement with a core-collapse scenario.
The CR content and γ-ray emission in IC443 have been scrutinized by multiple observatories. Fermi-LAT (Large Area Telescope) detected an extended γ-ray source in the 200 MeV -50 GeV energy band (Abdo et al. 2010, Ackermann et al. 2013. They showed that the spectrum is well reproduced by the decay of neutral pions and pinpointed the clouds B, C, D, F and G as targets of interaction. Tavani et al. (2010) reported the detection of γ-ray enhancement towards the NE shell in the 100 MeV -GeV range observed by AGILE (Astro-Rivelatore Gamma a Immagini Leggero). In their hadronic model, the cloud 'E' is the suggested target for the interaction with CRs. Interestingly, the location of the γ-ray peak differs between the lower energy emission (Fermi and EGRET detections, Esposito et al. 1996) and the VHE emission (VERITAS and MAGIC detections, Albert et al. 2007). This tendancy is also verified by Tavani et al. (2010) who located the 100 MeV γ-ray peak close to the Fermi/EGRET position and further towards NE. Recently, Humensky & VERITAS Collaboration (2015) presented an updated Fermi map where the position of the peak was shifted and consistent with VERITAS, Table 2: Observed lines and corresponding telescope parameters for the observations of the extended G region of IC443: ν is the frequency of the transition, FWHM corresponds to the full width at half maximum of the instrumental beam, F eff is the forward efficiency and B eff the beam efficiency of the dish. T sys is the average system noise temperature and ∆ the nominal spectral resolution. PWV is the precipitable water vapor, and r.m.s is the standard deviation measured on the baseline resampled with a spectral resolution of 0.5 km s −1 to allow direct comparison between the data from the IRAM-30m and APEX telescopes. exposing the uncertainty on the localization of the peak from the analysis of γ-ray observations. The choice of the region studied in this paper is based on the TeV γ-ray significance map from Humensky & VERITAS Collaboration (2015), that emphasizes the extended G region as a favorable target of high-energy CR interaction with dense molecular gas.
Explicit magnetic field studies towards IC443 are scarce. Wood et al. (1991) performed 6.1 cm polarimetric observations of the northeast rim of IC443. They found that the local magnetic field is rather correlated with the rim structure, but with no clear orientation (i.e. parallel or perpendicular) to it. Toward the clump G, only Hezareh et al. (2013) conducted a polarization study. They observed circular and linear polarization of the CO (1-0) and (2-1) lines with the Institut de Radioastronomie Millimétrique 30m antenna (hereafter IRAM-30m), and linear polarization maps from the dust continuum with the Acatama Pathfinder EXperiment (hereafter APEX). Their study constitute a crucial step towards the characterization of the interaction of CRs with the magnetic field in the extended G region.
Star formation in IC443 was the focus of a few studies, but early investigations lacked sufficient data to make clear findings. Recently, Xu et al. (2011) used color criteria for the point sources of IRAS and the Two Micron All Sky Survey (2MASS) to identify protostellar objects and Young Stellar Object (YSO) candidates and showed their association to several regions of the SNR interacting with neighbouring molecular structures, including clump G. Su et al. (2014) confirmed the shell structure of the distribution of YSO candidates using a selection method based on color-color diagrams inferred from WISE band 1, 2 and 3 and 2MASS band K. Both studies concluded that the formation of this YSO population is likely to have been triggered by the stellar winds of the progenitor.

APEX
A mosaic of the IC443G extended region was carried out with APEX 1 . APEX observations towards IC443G were conducted on September 11, 2018. The heterodyne receivers PI230 and FLASH345 (First Light APEX Submillimeter Heterodyne receiver, Heyminck et al. 2006), operating at 230 GHz and 345 GHz respectively, were used in combination with the FFTS4G and the MPIfR fast Fourier transform spectrometer backend (XFFTS, Klein et al. 2012). This setup allowed to cover a total field of 10 ×10 towards the center of the molecular region G. The observations were performed in position-switching/onthe-fly mode using the APECS software (Muders et al. 2006). Tab. 1 contains the spectroscopic parameters of the observed lines and Tab. 2 contains the corresponding observing setups 2 .
The central position of all observations was α [J2000] =6 h 16 m 37.5 s , δ [J2000] =+22 • 35 00 . The off-position used was α [J2000] =6 h 17 m 35.8 s , δ [J2000] =+22 • 33 08 , in the inner region of the SNR. We checked the focus during the observing session on the stars IK Tau and R Dor. We checked line and continuum pointing every hour locally on V370 Aur, Y Tau, IK Tau and R Dor. The pointing accuracy was better than ∼3 rms, regardless of which receiver we used. The absolute flux density scale was also calibrated on these sources. The absolute flux calibration uncertainty was estimated to be ≈ 15% during our observations. We used the GILDAS 3 package to calibrate and 1 This publication is partly based on data acquired under project M9508A_102 with the Atacama Pathfinder EXperiment (APEX). APEX is a collaboration between the Max-Planck-Institut für Radioastronomie (MPIfR), the European Southern Observatory, and the Onsala Space Observatory. 2 APEX telescope efficiencies are taken from http://www. apex-telescope.org/telescope/efficiency/index.php, and IRAM telescope efficiencies are taken from https://www. iram-institute.org/medias/uploads/eb2013-v8.2.pdf. 3 The Grenoble Image and Line Data Analysis Software is developed and maintained by IRAM to reduce and analyze data obtained with the 30m telescope and Plateau de Bure interferometer. See www.iram.fr/ IRAMFR/GILDAS Article number, page 5 of 36 A&A proofs: manuscript no. 3473992 merge the data of all sub-fields to produce a mosaic, and extract the spectral bands containing the signal corresponding to the (2-1) rotational transitions of 12 CO, 13 CO, C 18 O, and the (3-2) rotational transitions of 12 CO and 13 CO. The reduction included 1 st order baseline subtraction, spatial, and spectral regridding. The final products of this reduction process are spectral cubes centered on the previously cited rotational lines with a resolution of 0.5 km s −1 to increase the signal-to-noise ratio, which is more than enough for the expected linewidths (nominal spectral resolutions are indicated in Tab. 2). The observed area is showed in the white box in Fig. 1.

IRAM-30m
The same mosaic of the IC443G extended region was carried out with the IRAM-30m 4 . IRAM observations towards IC443G were conducted during one week, from February 20, 2019 to February 24, 2019. The heterodyne receiver EMIR (Eight MIxer Receiver, Carter et al. 2012), operating at 115 GHz and 230 GHz simultaneously, was used in combination with the FTS200 and VESPA backends. This setup allowed to cover the same total field of 10 × 10 towards IC443. The observations were performed in position-switching/on-the-fly mode. Tab. 2 contains the corresponding observing set-ups.
The central position of all observations was the same as the one used for APEX observations, α [J2000] =6 h 16 m 37.5 s , δ [J2000] =+22 • 35 00 . The off-position used was α [J2000] =6 h 17 m 54 s , δ [J2000] =+22 • 47 40 , in the northeastern ionized region of the SNR. We checked the continuum pointing and focus every hour during the observing sessions on several bright stars. The pointing accuracy was better than ∼3 rms. The absolute flux calibration uncertainty was estimated to be ≈ 15%. We also used the GILDAS package to calibrate and merge the data of all sub-fields to produce a mosaic, and extract the spectral bands containing the signal corresponding to the (1-0) rotational transitions of 12 CO, 13 CO, C 18 O and (2-1) rotational transition of 12 CO. The reduction process and final products are identical to what was presented in the previous subsection. The observed area is showed in the white box in Fig.  1.

Comparison between the IRAM-30m and APEX telescope data cubes
We estimated the systematic errors between the IRAM-30m and APEX based on the data cubes corresponding to the observation of the rotational transition 12 CO(2 − 1) by both telescopes. For visual comparison of the two data cubes, the channel maps are given in Fig. 3 (IRAM-30m) and Fig. D.3 (APEX). We found no evidence of systematic pointing error between the two data cubes. We also performed a quantitative comparison of the two spectral cubes. First, we resampled the data cubes in order to get the same spatial and spectral resolution. The spatial resolution was set to the nominal resolution of APEX θ = 28.7 , and the spectral resolution was set to the nominal resolution of the IRAM-30M ∆ = 0.5 km s −1 . Then, in each frequency channel and every single pixel of the mosaic we compared the signal detected by the two telescopes where the signal is greater than 3σ. The results of this complete investigation are represented on a T APEX mb vs T 30m mb 2D histogram shown in Fig. 2. We determined vs T 30m mb 2D histogram representing the complete comparison between the data cubes obtained with the IRAM-30m and APEX telescopes observations of the transition 12 CO(2−1) in the extended G region. The two cubes were resampled to the same spatial and spectral resolution. Every spectral channel in every single pixel is compared and displayed in the histogram. The dashed black line represents the 1:1 relation expected between the two spectral cubes if the data was identical. The solid black line represents represents the empirical relation that is measured by determining the best linear fit corresponding to the data dispersion, using a treshold of 1.5 K in order to best describe the high signal-to-noise bins.
the best linear fit x → a · x + b corresponding to the data dispersion between the two telescopes using a treshold of 5σ in order to best describe the high signal-to-noise data points. The parameters given by the χ 2 -minimization are a = 0.88; b = -0.35, indicating a slight overestimate in the measurement of the flux by the IRAM-30m with respect to APEX measurements, at least in the scope of our observations. Our statistical analysis of the two data cubes shows that: 1. Our measurements are characterized by a systematic error of approximately 12%. This could be due to inaccurate correction of telescope efficiencies and/or absolute flux calibration. 2. According to the dispersion around the instrumental linear model, our measurements are also affected by a random error characterized by a standard deviation of 305 mK.

Morphology of the region
Our 10 ×10 , ∼10-30 resolution maps of 12 CO and 13 CO in the extended G region provide a detailed picture of the morphology of molecular clumps. Fig. 4 (left panel) shows the emission in 12 CO(2-1) integrated between = −40.0 km s −1 and = +30 km s −1 mapped with the IRAM-30m. This wide interval of velocity includes all the components of the signal that are detected within the bandpass of our observations towards the extended G region. This map reveals a structured region with two main molecular structures that are spatially separated. The first structure (bottom-center in Fig. 4)   of 578 K km s −1 , a magnitude higher than the peak integrated intensity of the second structure (top-right in Fig. 4) that is around 63 K km s −1 . Fig. 3 shows the channel maps corresponding to the IRAM-30m observations of the 12 CO J=2-1 transition, which gives the finest spatial resolution of all our observations: θ=11.2 , or ∼0.1pc at the adopted distance of 1.9 kpc. The channel maps corresponding to other transitions of 12 CO and its isotopologue 13 CO are also available in the appendix:  (Hewitt et al. 2006). These structures, noticeable between = −7.5 km s −1 and = −1.5 km s −1 might either be corresponding to a slice of turbulent medium driven by the SN shockwave and/or belonging to the ambient gas associated with the NW-SE molecular cloud in which IC443G is embedded (Lee et al. 2012). Other than that, the description of the region probed by our observations can be divided into six distinct structures: 1. Cloudlet. In the upper part of the field we observe a large (∼5'×2', i.e. ∼2.8×1.1 pc) elongated cloudlet detected between = −7.0 km s −1 and = −5.5 km s −1 (indicated by the letter 'A' on Fig. 3), which is also detected in 12 CO(3-2). This structure was labelled G1 by Zhang et al. (2010), as part of the double peaked morphology of the extended G region. The 13 CO J=1-0 counterpart of this structure is much  brighter than the other main structures in the field, and it is also detected in the transitions J=2-1 and J=3-2, as well as in C 18 O J=1-0 and J=2-1. This structure was also presented and characterized by Lee et al. (2012) who proposed the label SC 03, among a total of 12 SCs (of size ∼1') found in IC443. 2. Ring-like structure. A ring-like structure seemingly lying in the center of the field (indicated by the letter 'B' on Fig. 3), appearing between = −5.5 km s −1 and = −4.5 km s −1 and also detected in 12 CO(3-2). It has a semi-major axis of 1.5', or 0.8 pc. This structure might be spurious and is likely to be physically connected to the elongated cloudlet as both are spatially contiguous and their emission lines are spectrally close. It is partially detected as well in our observations of 13 CO J=1-0, J=2-1 and J=3-2, and also has a faint, partial counterpart in C 18 O J=1-0 and J=2-1. To understand the nature of this region we searched for counterparts in Spitzer-MIPS, WISE, DSS, XMM-Newton, as well as in near-infrared and optical point source catalogues (Sect. 5), without success. Due to projection effects, this apparent circular shape could also be explained by an unresolved and clumpy distribution of gas. 3. Shocked clump. In the lower part of the field we identify a very bright clump emitting between = −31.0 km s −1 and = 16 km s −1 . This structure of size ∼2'×0.75' (∼1.1×0.4 pc), which is detected in the 12 CO(3-2) transition as well, belongs to the southwestern ridge of the molecular shell of the SNR and has been described as a shocked molecular structure by several studies (van Dishoeck et al. 1993, Cesarsky et al. 1999, Snell et al. 2005, Shinn et al. 2011, Zhang et al. 2010. The core of the shocked clump (indicated by the letter 'C' on Fig. 3) is also detected in 13 CO in the transitions J=1-0, J=2-1 and J=3-2, as well as in C 18 O J=1-0 and J=2-1. 4. Shocked knot. An additional shocked knot (indicated by the letter 'D' on Fig. 3) is also detected to the west of the previ-ously described structure. This fainter and smaller structure is spatially separated from the main shocked clump. 5. At the same position as the shocked clump and extending southward and westward, we find a faint, elongated clump emitting between = 5.0 km s −1 and = 7.5 km s −1 . This structure (indicated by the symbol '*' on Fig. 3) is spatially coinciding with the shocked clump , yet the peak velocity is not the exact same (also see below developments on kinematics of the region in Sect. 3.3). It has a faint counterpart in 13 CO(1-0). Observations of the ambient molecular cloud by Lee et al. (2012) indicate this structure as part of a faint NE-SW complex of molecular gas in the velocity range +3 km s −1 < LSR < +10 km s −1 . 6. Finally, the 13 CO(1-0) map (Fig. 4,right panel and Fig. D.2) indicates a large clump of gas extending from the bottomcenter to the right end of the field, with a bright knot in the bottom-right corner of the field. However, this structure has no bright, well-defined counterpart in any of the 12 CO transitions maps. It is spatially and kinematically correlated with the faint and diffuse 12 CO J=1-0 and J=2-1 emission seen in the velocity range -5.5 km s −1 < LSR < -2 km s −1 . From the comparison with the 12 CO observations of Lee et al. (2012) and 13 CO observations of Su et al. (2014) towards the SNR, we conclude that this structure is part of the western molecular complex observed in the velocity range -10 km s −1 < LSR < 0 km s −1 .

Kinematics of the region
Using the nominal spectral resolution of 0.5 km s −1 attained with our IRAM-30m 12 CO J=2-1 observations, we identified several velocity components of the molecular gas in IC443G based on the determination of the first and second moment maps (Fig. 5, left and center) of the 12 CO(2-1) data cube. We also produced Fig. 5: Left: first moment map of the IRAM-30m 12 CO(2 − 1) data cube. Center: second moment map of the same data cube. Right, upper panel: zoom into the dashed box on the first moment map of the APEX 12 CO(3 − 2) data cube to enhance the spectral resolution. Right top and bottom pannels display the first and second moment maps of the APEX 12 CO(3-2) data cube into the dashed box shown in the right and middle pannels, respectively. The colorbar of the left figure is centered on the velocity of IC443 in the local standard of rest, LSR = −4.5 km s −1 ). Table 3: Summary of the spectral characteristics of the lines of 12 CO and 13 CO within the boxes corresponding to each structures, defined in Fig. 6. T * mb is the peak temperature of the line, 0 is the centroid and v FWHM the linewidth. T * mb and the area are directly measured from reading the spectral data. v FWHM and 0 are measured by modelling the line either by a single or double gaussian function that best fits the data. When two gaussian functions are needed to model the sum of the signal emitted by the ambient cloud and structure of interest, we discard the contribution of the ambient cloud and only take into account the parameters 0 and FWHM associated with the broad component. The asymetric high-velocity wings of 12 CO lines corresponding to the shocked clump require two gaussians, thus we measure 0 = (c 1 + c 2 )/2 and FWHM = ((2.355 σ 1 ) 2 + (2.355 σ 2 ) 2 ) where c 1 , c 2 , σ 1 and σ 2 are the centroids and standard deviations corresponding to each gaussian function, respectively.
the moments maps of the 12 CO(3-2) data cube towards the ringlike structure (Fig. 5, right) to profit from the spectral resolution of 0.1 km s −1 .
1. The cloudlet has a mean velocity of about LSR = −5.5 km s −1 that is remarkably uniform throughout the structure. We measured LSR = −5.7 ± 0.3 km s −1 from the cen-troid of the 13 CO lines, contrasting with the velocity of IC443 in the local standard of rest by more than 1 km s −1 . It is likely that this discrepancy is due to a distinct velocity component with respect to the rest of the molecular gas in the extended G region. If this is not the case, this velocity shift could correspond to a maximum displacement of the kinematic distance Article number, page 9 of 36 A&A proofs: manuscript no. 3473992 ∆d ≈ 300 pc (following Wenger et al. 2018). Yet, the velocity wings of the 12 CO lines are still within the velocity range of the maser source in IC443G. The second moment map reveals a much lower velocity dispersion within the cloudet than for the shocked clump. It varies between ∼5 km s −1 and ∼7 km s −1 , which is slightly higher than the velocity dispersion across the background field, around ∼ 4 km s −1 . 2. The apparent ring-like structure is further analysed in the two right panels of the Fig. 5 where the first and second moment maps are determined for the 12 CO(3-2) data cube obtained with APEX. The superior spectral resolution of the APEX data cube offers a better precision in the determination of the moment maps, at the cost of a lower spatial resolution.
In the first moment map the mean velocity gradient within the ring is suggesting that the structure is rotating or expanding isotropically, as the mean velocity field varies between LSR = −4.7 km s −1 and LSR = −5.8 km s −1 from the western to the eastern arc of the ring. This apparent velocity gradient could be also due to systematic velocity variations between two or more distinct sub-clumps that are not well resolved by our J=3-2 observations (θ = 19.2 , ∼0.2 pc). The velocity dispersion measured within the ring-like structure varies between 1 km s −1 and 3 km s −1 , with a positive gradient from the eastern part to the western part of the structure where it spatially connects to the cloudlet. 3. The shocked clump has a mean velocity varying between LSR = −6 km s −1 and LSR = −8.5 km s −1 throughout its structure. These mean velocities are subject to caution as the self-absorption and asymmetric wings characterizing the line emissions of 12 CO might bias the value of the centroid. In fact, careful measurement of the centroid of 13 CO lines using a single gaussian function favors a velocity centroid of -4.4±0.2 km s −1 for the shocked clump, which is consistent with the velocity v LSR = −4.5 km s −1 of the maser in IC443G (Hewitt et al. 2006) . The second moment map displays important velocity dispersions, spanning from ∼15 km s −1 and up to ∼36 km s −1 within the shocked gas, increasing towards the center of the clump. 4. The shocked knot has a mean velocity v = −9 km s −1 that is slightly shifted with respect to the shocked clump. The second moment map shows an uniform velocity dispersion of ∼ 25 km s −1 , similar to the dispersion measured within the main shocked molecular structure.
The rest of the field of observation has a quasi-uniform mean velocity of ∼-4.1 km s −1 to ∼-2.8 km s −1 , which is slightly different than the mean velocity of IC443G in the local standard of rest but consistent with the ambient NW-SE molecular cloud in which IC443G is embedded (Lee et al. 2012). The velocity dispersion of this ambient gas is spanning from <1 km s −1 to ∼10 km s −1 in a few areas where the velocity dispersion is locally enhanced, with an average of ∼4 km s −1 . Excluding the contribution of the shocked structures and localized high-velocity dispersion knots, velocity dispersions measured from the 12 CO J=2-1 line in the extended G field span a range of r.m.s velocity σ v = 0.4 to 1.3 km s −1 in the ambient gas, and σ v = 1.2 to 1.8 km s −1 towards the cloudlet. At a temperature of 10 K, the thermal contribution is σ v = 0.32 km s −1 and it is likely that small-scale motions within the complex of molecular gas contribute to the measured dispersion, hence the ambient cloud is mostly quiescent, with turbulent motions smaller than 1 km s −1 . The velocity dispersion measured towards the cloudlet with the 13 CO lines is σ v = 0.8±0.1 km s −1 , which is consistent with typical molecular condensations (Larson 1981). Thus, we dot not find any kinematic signature of interaction of the cloudlet nor the ambient cloud with the SNR shocks in the extended G region, except for the few localized high-velocity dispersion knots.

Spatial separation of spectrally uniform structures
We aim to measure the mass associated with each molecular structures described in the previous section. We defined spatial boundaries enclosing these structures and independently studied the spectral data corresponding to each sub-region of the field. The spatial boxes defined for the cloudlet (A), ring-like structure (B), shocked clump (C) and shocked knot (D) are shown in Fig. 6, and the average spectra obtained in these boxes are presented in the Fig. 7 and Fig. 8 for every line from CO and its isotopologues that are available in our IRAM-30m and APEX data cubes. The choice of the boundaries is based on our morphological classification, but we carefully checked that the brightest spectral features are coherent across the different boxes that we defined (coordinates of these boxes are given in Tab. A.1). We performed that selection manually, as the size of our sample is not large enough to apply statistical methods (e.g. clustering, see Bron et al. 2018). Based on the analysis of the emission of 12 CO, 13 CO and C 18 O lines, our description of these spectral features is the following: 1. Cloudlet. Towards box A (Fig. 8, left-panel), the line profile of 12 CO and 13 CO lines are similarly double-peaked, and best modeled by the sum of two gaussian functions centered on the systemic velocities LSR = −5.7 ± 0.3 km s −1 (associated with the cloudlet) and LSR = −3.3 ± 0.1 km s −1 (associated with the ambient cloud).  Fig. 6 for the following lines: 12 CO(1-0), 13 CO(1-0) and C 18 O(1-0) (in black, IRAM-30m); 12 CO(2-1), 13 CO(2-1) and C 18 O(2-1) (in blue, APEX); 12 CO(3-2) and 13 CO(3-2) (in gray, APEX). Spectral cubes were resampled to allow direct comparison between the different spectra. Spatial resolutions of all transitions were modified to the nominal resolution of C 18 O(2-1), θ = 30.2 . Spectral resolutions were set to 0.5 km s −1 for IRAM-30m data, and 0.25 km s −1 for APEX data. On both panels, the LSR of IC443 is indicated with a vertical dashed line (at −4.5 km s −1 ).
2. Ring-like structure. Towards box B (Fig. 8, right-panel), the 12 CO and 13 CO lines are double-peaked as well. The use of two gaussian functions to model the line profile yields the systemic velocities LSR = −5.6 ± 0.2 km s −1 (associated with the ring-like structure) and LSR = −3.3 ± 0.1 km s −1 (associated with the ambient cloud). The gaussian decomposition is very similar to that of the cloudlet, suggesting that the apparent ring-like structure might be incidental despite its remarkable features in the first moment map (Fig. 5).
3. Shocked clump. Considering the geometry of the SNR and the locally perpendicular direction of propagation of the SNR shockwave (van Dishoeck et al. 1993), the highvelocity emission arises from at least two shock waves, if not a collection of transverse shocks propagating along the molecular shell. In other words, the projection along the line of sight of several distinct shocked knots with distinct systematic velocities could contribute to the broadening of the 12 CO lines. We measure s 27 km s −1 and s 21 km s −1 respectively for the blueshifted and redshifted transverse shocks. Except for the J=1-0 spectrum where the emission of the ambient gas contributes to the average spectra, all spectra of 12 CO lines exhibit a significant absorption feature around the LSR of IC443G, suggesting that there is strong self-absorption of the emission lines. Evidence of line absorption is found in the velocity range -6 km s −1 < v LSR < -2 km s −1 , which is where we detect the spatially extended features associated with the NW-SE complex of molecular gas described by Lee et al. (2012). Hence, it is possible that the foreground cold molecular cloud is at the origin of the absorption of the 12 CO J=1-0 and J=2-1 lines. A faint and thin emission line is detected around = 6.5 km s −1 both in the 12 CO and 13 CO spectra. This signal is associated with the NE-SW complex of molecular gas described in section 3.2. 4. Shocked knot. The shock signature of this line is distinct from the shocked clump. As hinted by the moments map ( Fig. 5), its fainter high-velocity wings are displaced towards negative velocities. A self-absorption feature is also observed in this structure. Between v = −5.5 km s −1 and v = −2 km s −1 a bright and thin feature traces the ambient gas shown on the channel maps on Fig. 3 (second row from bottom; first, second and third panels from left).
The linewidth measured from the 13 CO line profiles when we consider only the spectral component that are physically associated with the cloudlet and ring-like structure (discarding the contribution of the ambient cloud) are respectively 2.0 ± 0.3 km s −1 and 1.6 ± 0.3 km s −1 , measured by carefully defining much more constrained spatial boundaries around the structures. From the average spectra presented in this section, there is no spectral evidence for the propagation of molecular shocks and/or outflows within these two structures, except for the faint wings displayed by the 12 CO lines in the box associated with the ring-like structure. These extended wings arise from the contamination by the high-velocity emission of the shocked sub-structures that are contained in box B (Fig. 5). We measured the peak temperature, velocity centroid, FWHM and area of the 12 CO and 13 CO lines   Fig. 6 for the following lines: 12 CO(1-0), 13 CO(1-0) and C 18 O(1-0) (in black, IRAM-30m); 12 CO(2-1), 13 CO(2-1) and C 18 O(2-1) (in blue, APEX); 12 CO(3-2) and 13 CO(3-2) (in gray, APEX). Spectral cubes were resampled to allow direct comparison between the different spectra. Spatial resolutions of all transitions were modified to the nominal resolution of C 18 O(2-1), θ = 30.2 . Spectral resolutions were set to 0.5 km s −1 for IRAM-30m data, and 0.25 km s −1 for APEX data. On both panels, the LSR of IC443 is indicated with a vertical dashed line (at −4.5 km s −1 ).

2D histograms of CO data
In the next section we aim to build population diagrams in which we correct the effect of optical depth on the column density of upper levels. To measure the optical depth of CO lines, we relied on several strong assumptions, in particular the adopted isotopic ratios and the identity of excitation temperature for 12 CO and 13 CO (see section 4.2.2 for a description of our method, and Roueff et al. 2020 for a complete discussion of the corresponding assumptions). To assess the validity of this approach and estimate the key parameters (isotopic ratios, excitation temperature ratios), we built 2D histograms from the 12 CO J=1-0, J=2-1, J=3-2 data cubes, as well as 13 CO and C 18 O J=1-0 data cubes to compare the line intensity from different isotopologues and rotational transitions of CO with modified LTE models (Bron et al. 2018). These assumptions are also discussed in section 4.2.3. We examined four different 2D data histograms: In order to build the first three data histograms, we convolved all IRAM-30m data cubes to the nominal spatial resolution of C 18 O(1-0), i.e. 23.6 and to the nominal spectral resolution of the FTS backend, i.e. 0.5 km s −1 . To build the fourth histogram, we resampled IRAM-30m 12 CO(1-0), APEX 12 CO(2-1) and APEX 12 CO(3-2) data cubes to the nominal spatial resolution of 12 CO(1-0), i.e. 22.5 and to the nominal spectral resolution given by the FTS backend, i.e. 0.5 km s −1 . We used a treshold of 3σ to select data points where the signal is significantly above the noise level. The resulting 2D data histograms are shown in Fig. 9, where we emphasize the high signal-to-noise areas of the third histogram with black contours at 3σ and 4σ.

Results
1. The first histogram is characterized by a high signal-to-noise ratio and represents a large statistical sample (n = 20531). The 13 CO(1-0) vs 12 CO(1-0) relation presents at least two distinct branches. The lower quasi-horizontal branch, tracing bright 12 CO(1-0) emission associated with faint 13 CO(1-0) line emission (T 13 < 1 K km s −1 ). This branch is spatially correlated with the shock structure and spectrally correlated with the high-velocity wings of the 12 CO line that have no bright 13 CO counterpart because of insufficient integration time. It is spatially correlated mainly with the quiescent molecular gas that is found within the cloudlet and ring-like structure, as well as the ambient cloud.  hardly detected at a 3σ confidence level within our data cube. Still, we find evidence of at least one stastistically significant branch. Due to the small size of the sample it is not possible to identify any spatial or spectral correlation with certainty. 3. The third histogram has a poor statistical sample for the same reason as the second one (n = 833). The [ 13 CO(1-0)] / [C 18 O(1-0)] ratio vs [ 12 CO(1-0)] / [ 13 CO(1-0)] relationship is localized in an area with little dispersion. Hence, for high signal-to-noise measurements the isotopic ratios are almost uniform in the field of observations. Nonetheless, the lower signal-to-noise data bins display a statistically well-defined comet-shaped branch extending from this area. This branch might correspond to distinct physical conditions and/or isotopic ratios for a fraction of the field of observations. The high signal-to-noise ratio area of this branch is spatially correlated with the cloudlet and ring-like structure, whereas the 'tail' of the branch is spatially and spectrally associated with the shocked clump.
4. The fourth histogram is built on a large statistical sample with high signal-to-noise data (n = 9793), as the rotational lines J=1-0, J=2-1 and J=3-2 are well detected in our data cubes. 1.3. The crescent-shaped branch is highly correlated to the high-velocity wings of the 12 CO lines tracing the emission of shocked gas, where it is expected to measure high, enhanced J=2-1 / J=1-0 ratio (e.g. Seta et al. 1998).
We compared the observational histograms with synthetic families of curves generated using modified LTE radiative transfer models of the observed rotational transition for 12 CO, 13 CO, and C 18 O. Assuming that CO is at thermal equilibrium, the intensity of the line integrated over an element of spectral resolution Article number, page 13 of 36 A&A proofs: manuscript no. 3473992 (∆ = 0.5 km s −1 ) is given by: where α = 12, 13, 18 (histograms 1-3) or α =1-0, 2-1, 3-2 (histogram 4) respectively for the 12 CO, 13 CO, C 18 O isotopologues and J=1-0, J=2-1, J=3-2 rotational lines. T α is the integrated intensity of the line, defined as: where T 0 α = hν/k B is the temperature of the transition or energy of the upper level E up , given in Tab. 1 for all studied lines. τ α is the optical depth of the line, T cmb is the cosmic microwave background temperature, and T exc α is the excitation temperature of the line. The optical depth τ α and its dependence on the excitation temperature is described by: where τ T 0 α is the optical depth at a kinetic temperature T 0 . We use the optical depth at T 0 = 20K as a reference. We made the assumption that the excitation temperature of the different isotopologues (T 12 exc , T 13 exc , T 18 exc ) and rotational lines (T 1−0 exc , T 2−1 exc , T 3−2 exc ) are distinct and controled their ratios as four supplementary independent parameters of the modified LTE models (Roueff et al. 2020). The opacities of 12 CO and C 18 O isotopologues were determined from the optical depth of 13 CO using the corresponding isotopic ratio: Lastly, we assumed a hierarchy in the optical depth of the different rotational lines such that τ 1−0 ≥ τ 2−1 ≥ τ 3−2 and control their ratio as two supplementary parameters of the modified LTE models (Roueff et al. 2020). Hence, there are a total of 9 control parameters that we can set: the optical depth τ 13 , the isotopic ratio . We used the kinetic temperature of 13 CO as a variable of the parametric equation I i j (T 13 kin ) to synthesize LTE models that we can plot on top of each 2D data histogram. Each curve in Fig. 9 corresponds to a given set of control parameters, with the kinetic temperature growing linearly along a curve. We produced a family of LTE models with linearly increasing optical depth and a kinetic temperature varying between 10 K and 25 K, corresponding to the typical temperatures of cold molecular clouds. We adopted the isotopic ratio [ 12 CO]/[ 13 CO] and [ 13 CO]/[C 18 O] that offers the best visual match with our data, and similarly fine-tuned the other control parameters to reproduce the branches observed in each 2D data histogram.
All other parameters remaining constants, each control parameter has an effect on the pattern and location of a curve in the histogram space: -The modification of the temperature ratios alters the shape and orientation of the curves in histograms 3 and 4, and determines the boundaries of the 12 CO and C 18 O curves.
The LTE models that visually best match our data are represented in Fig. 9 4) it was necessary to set their temperature ratios as greater than 1, such that we infer an average temperature of up to ∼55 K traced by the J=2-1 lines, and up to ∼120 K for the J=3-2 lines in the high-velocity wings. This area of the histogram 4 also required lower τ 2−1 /τ 1−0 and τ 3−2 /τ 1−0 ratios, that we set to 0.2 for our models, suggesting that optical depth decreases significantly in the high-velocity wings of the lines. The higher boundary of our 12 C/ 13 C ratio estimate is anomalous with respect to the local interstellar medium average of 62±4 (Langer & Penzias 1993), however Wilson & Rood (1994) provide an estimate of the elemental abundances based on their distance to the Galactic center (GC): 12 C/ 13 C = (7.5 ± 1.9)D GC + (7.6 ± 12.9) 16 O/ 18 O = (58.8 ± 11.8)D GC + (37.1 ± 82.6) If we adopt the distance of 1.9 kpc for IC443 (Ambrocio-Cruz et al. 2017) the SNR is located at a distance of ∼10 kpc from the GC. Hence our estimate of the 12 CO/ 13 CO ratio is in agreement with the expected abundance of 80±30, but our measured 13 CO/C 18 O is much higher than the predicted ratio of 8±2. The measured enhancement of the 13 CO/C 18 O ratio might be the product of the fractionation of carbon monoxide by photodissociation, as the shielding of the various isotopologues of CO provides a mechanism for the rarefaction of C 18 O with respect to the expected abundance (e.g. Glassgold et al. 1985). We checked the spatial distribution of data points corresponding to particular areas of the histograms presented in Fig. 9.
-Histogram 3. Interestingly, the data points corresponding to high (> 20) 13 CO/C 18 O line ratio are spatially correlated with the cloudlet only. If no protostars are positively identified towards the cloudlet (section 5), radiative decay from X-ray irradiation seen in the SNR cavity with XMM-Newton might provide a source to account for the selective photodissociation of CO within the cloudlet (Troja et al. 2006, Troja et al. 2008). -Histogram 4. Data points corresponding to high (> 2) [J =  Fig. 6) used to measure specific intensity corresponding to each structure. Three velocity ranges are given for the shocked clump, as we distinguish three components within the molecular lines: i. the left (blueshifted) high-velocity wing, ii. the core, iii. the right (redshifted) high-velocity wing.
[-1.5; +25] C shocked knot [-35; +25] D ambient cloud [-6.5; -1.5] A ∩ B ∩ C ∩ D ratio apparently tracing the shock front. Although, this combination of line ratios is not predicted by RADEX models and could be due to an inaccurate registration between the data sets.

Population diagrams
We determined the physical conditions (column density N CO , kinetic temperature T kin ) of the molecular gas in the extended G region assuming that the emission lines are thermalized. We used pixel-per-pixel population diagrams (Goldsmith & Langer 1999) corrected for optical depth for the 12 CO(3-2), 12 CO(2-1) and 12 CO(1-0) transitions to measure the upper level column density N up (hereafter level population) for J up = 3, J up = 2, J up = 1. Systematically, the 13 CO data is used to correct the effect of optical depth in the population diagrams, for each element of resolution (i.e. for each of line of sight) and for all velocity channels.
To perform this task we choosed to use the APEX data cube over the IRAM-30m data cube to measure the specific intensity corresponding to the 12 CO(2-1) transition. The IRAM-30m has the advantage of having a lower noise (see Tab. 2). However, the 13 CO(2-1) transition was only observed with the APEX telescope due to receiver setup constraints. Hence to avoid intercalibration issues, we used the APEX data for 12 CO(2-1) and 13 CO(2-1). First, for consistency we convolved all maps to the same spatial resolution, using the nominal resolution of 13 CO(2-1), such that we have a beam diameter of 30.1 for the six maps considered. The spectral resolution was also modified and set to 2 km s −1 for each transition in order to increase the signal-tonoise ratio. Then, the following measurements were performed pixel-per-pixel using a Python 5 algorithm: 1. Sigma-clipping was applied to the spectra of all transitions of the two isotopologues 12 CO and 13 CO, using a threshold of 3σ. 2. In every single velocity channel remaining after sigmaclipping, the N J=3 , N J=2 and N J=1 level populations were measured both for 12 CO and 13 CO to determine the optical depth of the three transitions using the adopted value of isotopic ratio. We measured the optical depths τ 3−2 , τ 2−1 and τ 1−0 assuming that 13 CO lines are optically thin and that the excitation temperatures of the two isotopologues are approximately equal (this second assumption is supported by the re-sults presented in section 9). Under this assumption we have: I ν, 12 CO dν = hν i j, 12 CO 4π N up, 12 CO A i j, 12 CO · β i j (6) I ν, 13 CO dν = hν i j, 13 CO 4π N up, 13 CO A i j, 13 CO (7) where β i j = [1 − e −τ i j ]/τ i j is the escape probability of a photon, ν i j the frequency of the transition, N up the level population, A i j the probability of the transition and I ν the specific intensity. Based on the comparison of our data with modified LTE models (Sect. 4.2.1) we adopted different values for the expected isotopic ratio towards each region (indicated in Tab. 6). 3. The estimates of 12 CO level populations were corrected for optical depth using the correction factor C τ i j = τ i j /(1 − e −τ i j ), such that the corrected column density is given by N * up = C τ i j N up . 4. We put the corrected level populations N * up and their corresponding energies E up = hν i j /k B into the population diagram N * up = f (E up ), and used a χ 2 -minimization algorithm 6 to determine the best linear fit y = ax + b. The excitation temperature T ex is deduced from the slope, and the total column density N CO is determined by computing the partition function Q(T ex ) 7 and measuring the offset. 5. We finally derived the mass of the gas from the measured total column density N CO , using the expected [H 2 ]/[ 12 CO] ratio towards dense molecular regions and assuming that the distance of IC443 is 1.9 kpc (Ambrocio-Cruz et al. 2017). We adopted the value of 10 4 for the H 2 -to-12 CO abundance ratio (Gordon & Burton 1976, Frerking et al. 1982), thus we have M H 2 = (2m H )AN CO × 10 4 where A is the area of the box and m H the mass of hydrogen.
Uncertainties. The errorbars on total CO column density, gas mass and kinetic temperature are determined by: (i.) instrumental errors (dominated by the systematic uncertainty on the flux, Sect. 3.1.3); (ii.) the quality of the linear fit applied to the population diagram; (iii.) systematic errors on the adopted distance of IC443 and the 12 CO/ 13 CO isotopic ratio.
Results. Average population diagrams are presented on Fig.  10 and Fig. 11. The measurements obtained are presented in Tab. 6, in which the minimum and maximum boundaries are given for the total CO column density, gas mass and kinetic temperature. We measured the column densities and masses for the regions A (cloudlet), B (ring-like structure), C (shocked clump) and D (shocked knot) as indicated in Fig. 6, and also for the ambient cloud by averaging the signal over the entire field (regions A, B, C and D excluded). We also measured the mass of the foreground clump that is spatially coinciding with the shocked clump using a gaussian model for the average CO lines in the velocity range [2.5; 12] km s −1 . The result of this measurement is presented at the end of Tab. 6. 6 The Python module scipy.optimize.curve_fit uses the Levenberg-Marquardt least square curve fitting algorithm. 7 We computed the partition function taking into account the first 40 upper rotational levels of 12 CO taken from the Cologne Database for Molecular Spectroscopy and Jet Propulsion Laboratory database. The thermal energy of the 40 th level is E 40 = 4513K.
Article number, page 15 of 36 A&A proofs: manuscript no. 3473992 Table 6: Results of the population diagrams analysis for each structure. N CO is the average column density and M the mass in solar unit. The adopted isotopic ratio 12 CO/ 13 CO for each region is indicated in the fourth column. In the mass column, we indicate in parenthesis the 'cold mass' estimated by considering only the upper levels J=1,2. T exc is the average excitation temperature. The quality of the linear fits is indicated by the average value of χ 2 . T w exc is the average excitation temperature of the warm component, obtained by considering the levels J=2,3 separately from the level J=1. C τ = M /M gives an indication of the enhancement of the mass estimate when we correct the optical depth of the transition J=1-0 independently of the detection of 13 CO for the transitions J=2-1 and J=3-2. For comparison, we indicate the estimate of the mass obtained by using the CO-to-H 2 conversion factor X CO = 2 × 10 20 cm −2 (K km s −1 ) −1 applied to the raw J=1-0 data (without optical depth correction).

Discussion
As a first-order verification we compared our measurement of the mass based on the population diagram with a rough estimate of the mass using the CO(J=1-0)-to-H 2 conversion factor X CO (Dame et al. 2001, Bolatto et al. 2013. We adopted the following conversion factor, with a ±30% uncertainty: The H 2 column density was determined by the product of X CO with the area of the raw 12 CO J=1-0 line (without optical depth correction), and the mass was inferred from N H 2 in the same manner as described previously. We report the results of this measurement in the last column of Tab. 6 where they are referred to as M J=1−0 X CO . Within errorbars, almost all our measurements are consistent with this rough estimate of the mass, except for the line wings towards the shocked clump. With respect to the population diagram method, this method systematically overestimates the mass. The enhanced J=2-1/J=1-0 emission of the high-velocity gas accounts for a lower mass in the population diagram, as it traces a warmer gas for which the X CO conversion factor yields an overestimate.
We built a single population diagram for each spatial box, using a single value of N up that is the spatial and spectral average of our sample of measurements. The resulting population diagrams are shown in Fig. 10, as well as the statistical informations on the spread of the sample around these average values of N up . As an additional information on the thermalization of carbon monoxide, in Fig. 11 we show the average population diagrams obtained using the 13 CO lines with the same method, without correction of the optical depth.
1. Thermalization. At a temperature of 10 K, the critical density n i j = A i j /C i j of the 12 CO J=1-0, J=2-1 and J=3-2 lines are respectively n 1−0 = 2.2 × 10 3 cm −3 , n 2−1 = 2.3 × 10 4 cm −3 and n 3−2 = 3.5 × 10 4 cm −38 . Towards the shocked clump, it is likely that these critical densities are attained (Cesarsky et al. 1999), hence the lines should be thermalized. In the population diagrams displayed in Fig. 10, the average data points are generally in satisfactory agreement with the assumption that the emission lines are thermalized, as there is no significant divergence from the Boltzmann distribution for any of the structure studied. The highest value of χ 2 is obtained towards the ambient cloud that also presents an abnormally low kinetic temperature, down to ∼5 K. This might indicate that the J=2-1 and J=3-2 are subthermal, which is expected if some parts of the NW-SE molecular cloud we are probing have a density lower than 10 4 cm −3 . There is a systematic discrepancy between the linear fit and the measured column density of the upper level J=2, which is lower than expected in each structure. This anomaly can be solved if we adopt two distinct kinetic temperature to model the relative distribution of the level populations. In the high-velocity wings of the lines, it is expected and it was hinted by our LTE analysis (Sect. 4.2.1) that the J=2-1 and J=3-2 lines trace a warmer gas than the J=1-0 line that is primarily tracing the quiescent and cold phase. In this case, the distribution in these population diagrams should be modeled by a linear fit of kinetic temperature T 1 for the upper levels J up =1,2 (cold component) and a linear fit of temperature T 2 for levels J up =2,3 (warm component), hence assuming that along a line of sight we observe two layer of gas thermalized at distinct temperatures, with T 1 < T 2 . In Tab. 6 we present the temperature obtained for the warm component, as well as the mass of the cold component. The statistical spread of our sample of measured level populations around the average linear distribution provides a strong motive for the use of pixel-per-pixel, channel-per-channel population diagrams, as it is evidence of a large range of physical conditions that require distinct LTE models.
2. Filling factor. The filling factor used to infer column densities from the main beam temperature measured during our observations is set to 1. From the morphology of the gas mapped in 12 CO(2-1) with the IRAM-30m with a nominal resolution of 11.2 (Fig. 4, Fig. 3) we consider that we are probing extended, clumpy structures with dimensions that are greater than the beam diameter characterizing our are indicated for two adopted values of the 12 CO/ 13 CO isotopic ratio ( 12 CO/ 13 CO = 60 in black, 12 CO/ 13 CO = 100 in red). The black solid line represents the linear fit of the black data points. The boundaries of the filled areas are defined by the linear fits corresponding to the first and third quartiles for each value of the isotopic ratio adopted. The 5 th percentile and 95 th percentile are indicated by the caps, the first quartile and third quartile are indicated by the ends of the box and the average is indicated by the cap inside of the box. Each diagram accounts for a sample of measurements performed in one of the spatial boxes defined in Fig. 6 and integration ranges defined in Tab. 5: a) the shocked clump, b) the cloudlet, c) the ring, d) the ambient cloud and e) the shocked knot.
observations, between and 19.2 and 30.1 for the data cubes used to measure the gas mass.
3. Optical depth. The measurement of the optical depth from the comparison of the fluxes between 12 CO and 13 CO is displayed in the appendix , as well as the measurement based on the 12 CO/C 18 O isotopic ratio (Fig. F.1). The optical depth is non negligible both in the ring-like structure and cloudlet. It is important in the self-absorbed component of the emission lines corresponding to the shocked clump where the cold, ambient cloud intervenes, and much lower in the high-velocity wings. We checked the assumption that 13 CO is optically thin using the C 18 O J=1-0 and J=2-1 lines to estimate its optical depth based on the expected 13 CO/C 18 O isotopic ratio. In our pixel-per-pixel, channel-per-channel population diagrams the correction of the optical depth is partially biased, as the detection of 13 CO is limited by the noise level. At a certainty level of 5σ, we cannot measure a flux lower than 125 mK, 420 mK and 400  mK respectively for the J=1-0, J=2-1 and J=3-2 transitions. As a consequence, towards a significant number of line of sights and velocity channels the optical depth is measured for the J=1-0 transition but not for J=2-1 and J=3-2. In the population diagram this results in the biased displacement of the level J=1 with respect to the levels J=2 and J=3, hence the slope of the linear fit is artifically enhanced and the excitation temperature is lowered (Fig. 10, Tab. 6). 4. Self-absorption of CO lines towards the shocked clump. We performed an independant measure of the temperature and column density of the absorbing envelope of the shocked clump from the absorption features of CO lines. Towards the position (6 h 16 m 42.5 s , +22 • 32 12 ), we renormalized the signal in each frequency channel with respect to the profile of the emission lines and measured the optical depth for all transitions in a beam of radius 22.5 . From the equivalent widths W of the transitions J=1-2 and J=0-1, we measured the temperature of the gas solving the following equation: Where T 0 = hν/k and F(u) = e −3u cosh(u) (see appendix B). We obtained an excitation temperature of 18 ± 3K. We then built a population diagram of the lower population N l levels from the measures of the optical depth, using the following Article number, page 17 of 36 A&A proofs: manuscript no. 3473992 relation: Where f ul is the oscillator force, λ the wavelength of the transition and σ v the velocity dispersion of the absorbing medium, which we assume equal to the velocity dispersion of the lines associated with the shocked clump that we measured in Tab. 3. This population diagram yield a total column density of the absorbing medium N CO = 3.4 × 10 16 cm −2 , corresponding to a total mass of 20 M for the absorbing envelope of the shocked clump.

LVG method
The results obtained in the previous section are based on several strong assumptions, in particular that the excitation temperature of 12 CO and 13 CO are equal. To adress this issue, we investigated the application of a second method that does not require this assumption to determine the mass of the gas towards our field of observations using CO lines. We used the LVG model for an expanding spherical shell (e.g. Sobolev 1960, Surdej 1977 to model the lines observed. This model require a local, large velocity gradient ∇v along the line of sight. We consider that this assumption is verified throughout the field of observation as the linewidth of 12 CO lines suggest typical turbulent broadening (Sect. 3.3). We used the non-LTE radiative transfer RADEX code (van der Tak et al. 2007), a computer program solving the radiative transfer equation based on the escape probability formulation in an isothermal and homogeneous medium characterized by a large velocity gradient. Three geometries are available, we selected the expanding spherical shell for which the escape probability is related to the optical depth τ by the following formula (Mihalas 1978): In the framework of the LVG assumption, the value of the optical depth at line center determined by RADEX is proportional to N/∆ where N is the column density and ∆ the full width at half-maximum of the line profile. Hence when using RADEX the choice of the parameter ∆ is critical to model the expected intensities of atomic and molecular lines. The program solves the statistical equilibrium equations taking into account up to seven collisional partners, allowing to analyze spectral line observations if the molecular collisional data is known. RADEX models require a total of 7 input parameters: the column density N, local density n H 2 , kinetic temperature T kin , line width ∆ and background temperature T b . In return, RADEX estimates the excitation temperature T ex , optical depth τ i j , peak temperature T R , upper and lower level populations N up , N low , and flux T dv of an arbitrary number of transitions defined by the user. We applied the LVG method to average spectra measured over the extent of each box defined in Fig. 6 and velocity ranges defined in Tab. 5. Using a χ 2 -minimization, we compared RADEX outputs with our measurements in order to constrain the physical conditions. We compared the line intensities of the three first rotational transitions of 12 CO and 13 CO with a grid of RADEX models. In the same manner as in the population diagram approach, for consistency we convolved all maps to the same spatial resolution, using the nominal resolution of 13 CO(2-1), such that we have a beam diameter of 30.1 for the six maps considered. The spectral resolution was also set to 2 km s −1 for each transition in order to increase the signal-to-noise ratio. Then we produced several three-dimensional grids of RADEX models with the varying parameters (N, n H 2 , T ) and fixed parameters (T b , ∆ ). These parameters are specified in Tab. 7 for the different grids used. We use large ranges of column densities and kinetic temperature in order to probe a consistent fraction of the space of parameters, corresponding to a variety of physical conditions as wide as possible. We set the value of the parameter ∆ based on the equivalent linewidth that we measured on average spectra using bigaussian functions to model the line profiles. The following steps were followed for each average spectrum using a Python algorithm: 1. We applied sigma-clipping to the spectra of all transitions of 12 CO and 13 CO, using a threshold of 3σ. 2. We measured the specific intensity of each transition on the velocity range corresponding to each spectral feature (e.g. the wings of the shocked clump are treated separately from the core of the line). Then, in each element of the model grid we computed the following quantity, corresponding to a combined reduced χ 2 statistical estimator: are respectively the specific intensity of 12 CO and 13 CO returned by a given RADEX model, s obs the observational intensity and σ the uncertainty associated with the measurements of the intensity s. 3. The physical conditions for each average spectrum were deduced from the minimization of the quantity χ 2 , i.e. we localized the minimum element of the resulting χ 2 RADEX grid and infer the quantities (N, n H 2 , T ) from the corresponding RADEX input.
The errorbars on column density and kinetic temperature were estimated from the uncertainties on the flux and more importantly from the uncertainty on the isotopic ratio 12 CO/ 13 CO.
Results. The results of our LVG analysis are presented in Tab. 7, in which the minimum and maximum boundaries are given for the CO column density and gas mass for the shocked clump, shocked knot, ring-like structure, cloudlet and ambient cloud. As previously, we measured the molecular mass assuming that the H 2 -to-12 CO abundance ratio is equal to 10 4 . The χ 2minimization was successfully attempted for the major part of our analysis. We present the χ 2 diagrams for each region in Fig.  12, where the first and second term of Eq. 12 are independently represented by two sets of filled contours. For most structures, we were not able to determine a precise measurement of the local density from the χ 2 -minimization, as the variation of χ 2 with respect to the choice of the input n H 2 does not strongly favor any LVG model. In all cases we observe that χ 2 increases significantly for densities n H 2 < 10 4 cm −3 , hence our analysis suggests that (i.) the 12 CO lines are thermalized; (ii.) the local density is greater than 10 4 cm −3 across the field of observations.

Discussion
The gas masses measured with the LVG approach are systematically lower than the masses obtained using population diagrams corrected for optical depth (Tab. 6), except for the highvelocity wings where the LVG estimate is higher by a factor ∼ 3 and for the cloudlet where both methods yield the same result. g) shocked knot Fig. 12: χ 2 -minimization diagrams of the comparison between 1×200×200 LVG grids modeled with sets of parameters (n H 2 , N CO , T ) and observations of 12 CO and 13 CO lines towards the boxes corresponding to the shocked clump (a, b, c), ring-like structure (d), cloudlet (e), ambient cloud (f) and shocked knot (g) defined in Tab. A.1. The shades of green represent the areas of minimum χ 2 corresponding to the minimization of the specific intensity of the 12 CO(3-2), 12 CO(2-1) and 12 CO(1-0) lines (first term introduced in Eq. 12), and the shades of pink represent the areas of minimum χ 2 corresponding to the minimization of 13 CO lines (second term introduced in Eq. 12. The filled contour levels are log 10 (χ 2 ) = [0.6, 1, 1.5, 2], with the first contour level highlighted by a solid black line. Hatched areas represent the intersection between the previously described sets of contours, defined by log 10 (χ 2 RADEX ) ≤ [0.5, 0.45, 0.45, 0.5, 0.75, 0.75, 0.5] (Eq. 12) respectively for the boxes a, b, c, d, e, f and g.
Deviation from a single excitation temperature is not sufficient to account for this discrepancy, since the LVG models that fit our data do not predict a disagreement higher than ∼1 K for the excitation temperatures of 12 CO and 13 CO. This discrepancy is mainly due to the fact that the LVG method is applied to average spectra, whereas the LTE approach is applied pixel-per-pixel and channel-per-channel. The channel-per-channel measures of optical depth are particularly different from the average measures towards the shocked clump where the optical depth varies strongly between the center of the line and the high-velocity wings. The highest estimates of the mass are obtained using the CO-to-H 2 conversion factor, based on the emission of 12 CO J=1-0 only. In comparison to previous measures of the molecular mass in the extended G region: 'cloud G' from the 12 CO J=3-2 line. Their measure included a larger field, and they lacked sufficient data to correctly estimate the optical depth.
Our new measurements of the mass are crucial for the interpretation of the interaction of CRs with the ISM. The extended G region is likely to be the main target of interaction with CRs, hence the source of TeV γ-ray emission in IC443. Our results put constraints on the amount of molecular mass that is avail- Table 7: Summary of the results and parameters used in the LVG analysis. We show ranges of values corresponding to the uncertainty interval. N CO is the column density of 12 CO, T kin the kinetic temperature of the gas, n the local density of H 2 , log 10 (χ 2 ) gives the minimum value of χ 2 within the minimization diagrams presented in Fig. 12, M is the mass measured using the [H 2 ]/[ 12 CO] ratio. N grid is the range of column densities parsed across the grid of models, T grid the range of kinetic temperatures, T b and ∆ are respectively the background temperature and the width of the emission line used by RADEX for all models. We indicate the ratio M LVG /M LTE of the estimate of the mass with the LVG approach with respect to the results obtained with population diagrams (Sect. able to interact via bremsstrahlung and pion decay mechanisms. Torres et al. (2010) showed that the characteristics of the γ-ray spectra in IC443 suggest the interaction with two distinct molecular structures: i. a lower mass cloud (∼350 M ) at distance of 4 pc from the SNR, and a higher mass cloud (∼4000 M ) at a distance of 10 pc. Their results were obtained from the analysis of the γ-ray spectra of the SNR in a larger field, yet they are consistent with out findings if we consider that the mass of the ambient gas that we measured could contribute to the distant component. As it has been suggested by Lee et al. (2012), either the shocked clump or cloudlet and ring-like structure could correspond to the closer component, while the distant component might correspond to the ambient cloud. Within errorbars, the mass we measured for these structures could fit with this scenario.

Protostar candidates
In order to characterize the local star formation over the extent of the molecular structures characterized in the previous section, we studied the distribution of optical, infrared and near-infrared point sources in the extended G region. We aimed to check if these infrared point sources can be identified as protostars and if so, to constrain their evolutionary stage based on their infrared fluxes. Finally, we aimed to study their spatial distribution and their association with the molecular clumps found in the extended G region.

Origin of the data
Our field was fully observed both by 2MASS and WISE. 2MASS operates between 1.235 µm and 2.159 µm, and WISE between 3.4 µm and 22 µm. The exact photometric parameters for each band of 2MASS and WISE are given in Tab. 8. To recover point sources from these two catalogues, we used the NASA/IPAC infrared science archive to obtain all entries in a 10 arcmin sized square box around the field center α [J2000] =6 h 16 m 37.5 s , δ [J2000] =+22 • 35 00 , corresponding to the same field that was mapped with IRAM-30m and APEX. A total of 487 point sources in the 2MASS All-Sky Point Source Cat- alog (PSC) and 515 point sources in the AllWISE Source Catalog were found using this query within the extended G region. Both catalogues provide position coordinates, photometric measurements for each band and their uncertainties, signal-to-noise ratio, as well as several flags specifying contamination by extended emission, quality of the PSF profile-fit and other possible sources of bias.

Selection of relevant IR point sources
In order to reject false positives, we applied several selection criteria to our primary catalogues of point sources detected by 2MASS and WISE: 1. We required a complete detection in the WISE bands W 1 , W 2 and 2MASS bands J, H and K (i.e. the measurement is not an upper limit). 2. We selected only the sources characterized by a signal-tonoise ratio greater than 2 for the photometric bands W 1 , W 2 , W 3 , J, H and K.
Article number, page 20 of 36 P. Dell'Ova et al.: Interstellar anatomy of the TeV gamma-ray peak in the IC443 supernova remnant 3. We reject the WISE point sources that were flagged for confusion and/or contamination of the photometric bands by image artifacts. 4. We rejected the 2MASS point sources that were flagged for low quality photometric measurements.
After this selection, 214/515 point sources remain for the AllWISE Source Catalog and 328/487 for the 2MASS PSC. A total of 99 point sources were detected both by 2MASS and WISE. AllWISE point sources marked with a value of ext_flg flag different than '0' are indicated as extended sources. Either their morphology is not consistent with the point spread function of any band or they are spatially associated with a known extended source of the 2MASS Extended Source Catalog (XSC). 2MASS point sources with gal_contam = 0 are also sources that fall within the elliptical profile of a known extended source. A search in the 2MASS XSC catalog shows that 17 extended sources are found in our 10 ×10 field. The shocked clump is particularly crowded, suggesting that the extended emission of bright knot of shocked material are detected by the survey in this area. The photometric measurements are likely to be contaminated by this extended emission, hence the identification of these sources as protostar candidates is uncertain. Nonetheless we do not completely rule out these entries in the catalog since this extended emission might be produced by outflows.

Color-color filtering
In order to select only the point sources that could possibly be young stellar objects, we applied a color-color criteria to our near-infrared point source catalogues. To identify the nature of a protostar candidate detected by 2MASS, we compared the relative flux in the J, K and H photometric bands. We used the following empirical color criteria (Xu et al. 2011), based on the idea that protostars have an infrared excess in the 1.235-2.159 µm range that determines their position in the JHK color-color space and that is directly related to their evolutionary stage (Lada & Adams 1992): Within the color-color diagram, this system of equations defines the color-color domains which mark out the different types of sources (see Fig. 13). This method allows to filter the sample of point sources and to produce a subset of different types of candidate young stellar objects (YSO): -Classical T Tauri stars (CTTS).
Similarly, we used the following color-color criteria (Koenig & Leisawitz 2014) to characterize the point sources from the relative flux in the three first bands W 1 , W 2 and W 3 . Fischer et al. (2016) proposed a slightly modified version of these filtering criteria that excludes an area of the color-color diagram that would be considered as Class II in the original diagram. That area is instead interpreted as shock emission, together with all the point sources that are beyond the left branches of the Class II and Class I areas. In this version of the diagram, there are two well-defined domains of color-color space: 1. A first region is defined by the following system of 4 equations that constrain the infrared excess observed in Class I young stellar objects: 2. A second, adjacent region is defined by the following system of 5 equations that constrain the infrared excess observed in Class II young stellar objects: In the color-color diagram represented in Fig. 13, this system of 9 equations allows to distinguish two samples from the rest of the catalogue and enables to sort the YSO candidates into two distinct expected evolutionary stage based on their infrared excess: -Class I protostar.
Additionnaly to the shocked emission area of the color-color diagram, Fischer et al. (2016) added the labels 'Polycyclic Aromatic Hydrocarbon (PAH) emission', 'AGB stars', 'Tr. (Transition) disks' and 'Star-forming galaxies' to differents parts of the diagram based on the expected emission of these objects in the WISE photometric bands. We defined arbitrary branches to separate the areas corresponding to each label and identified point sources falling into one of these area.
After the color-color filtering of the catalogues, the total amount of remaining point sources is 79/328 for 2MASS and 65/214 for WISE. 9 point sources are detected both by 2MASS and WISE. 1 point source is detected by Gaia and detected as a Class II protostar by WISE. The uncertainty due to extended emission and their identification in the color-color space is the following: 2MASS. 12.7% of the protostar candidates found in the 2MASS PSC are contaminated by extended emission. 23 CTTS candidates, 8 HAeBe candidates and 48 other YSO candidates are identified based on their JHK photometric measurements. WISE. 76.9% of the protostar candidates found in the AllWISE catalog are contaminated by extended emission towards the crowded shocked clump. 53 class I protostar candidates and 12 class II protostar candidates are identified in our field.

Spectral index α of protostar candidates in the range 3.4-12 µm
We measured the spectral index α = d log(λF λ )/d logλ of our protostar candidates detected by 2MASS and WISE to derive an identification based on their infrared spectral energy distribution (SED) (Adams et al. 1987). We used the classification system of Greene et al. (1994) to attribute an evolutionary stages to each point source and rule out the sources characterized by a flat SED: - -Class II: −1.6 ≤ α ≤ −0.3 -Class III: α < −1.6 The photometric fluxes of the bands W 1 , W 2 and W 3 were used to compute the slope of the SED in the range 3.4-12 µm for each protostar candidate found precedently and given in Tab Every single point sources classified as Class I by our measurement of α were also identified as Class I using the color-color diagram (ID 1-4 and ID 16-34). On the contrary, there is only partial agreement between the two methods for the identification of Class II protostar candidates: excluding Flat-SED sources, 75% of the identifications were confirmed by both approaches .
No Class III were identified in our sample.

Gaia point sources
We used Gaia data (see Tab. 8 for exact photometric parameters) to complete our point source census with optical sources and search for multiple detection by Gaia, WISE and 2MASS. We used the ESA Gaia science archive to obtain all entries in a 10 arcmin sized square box around the field center. A total of 468 point sources were found within our 10 ×10 field of observations. In order to estimate the amount of optical point sources that might be evolved stars physically associated with the SNR, we also used a distance criterion to filter the result. Assuming that the distance of IC443 is ∼1.9 kpc, we rejected the point sources with relative distance greater than 2500 pc or lower than 1500 pc, based on the lower and upper bound on the confidence interval for the estimated distance determined by parallax measurements Bailer- Jones et al. (2018). There are 16 point sources in the interval 1500-2500 pc. We repeated the same process with a different distance criterion in order to assess the variability of the result with respect to the amplitude of distance interval applied as a filter. We found a total of 40 optical point sources which relative distance is greater than 1250 pc and lower than 2750 pc. We checked if our samples of optical point sources is correlated with candidate protostars found in our WISE and 2MASS census: A greater amount of point sources is found on the edges of the 24µm bright knot that is spatially associated with the eastern edge of the quiescent molecular cloudlet. Excluding the shocked clump, there is an anti-correlation between both 2MASS and WISE sources and the MIPS-24µm flux map, in particular in the vicinity of the cloudlet. We suppose that the absence of optical point sources within this region is due to the extinction caused by this massive dark cloud. No IR point sources are found at the center of the gap within the ring-like structure.
Interestingly, the few point sources that were not flagged for contamination by extended sources are primarily organized in two clusters in the eastern and north-eastern part of the field (group 1 and group 2, Fig. 14). All these candidates were positively identified by both methods as Class I and Class II protostars. The Class I protostar candidates of group 2 are well associated with Spitzer-IRAC filamentary structures at 4.5 µm, at the edge of a faint structure detected by Spitzer-MIPS at 24 µm. This bright filament can be either part of the SNR shock or unresolved outflows. Most IR sources with protostar signatures are spatially correlated with the molecular shell in IC443G, suggesting enhanced star formation along the rims. Given the protostellar collapse phase timescale of ∼10 5 yr (Lefloch & Lazareff 1994) we dot not establish any causal relation between the SNR shocks and the formation of Class I and Class II protostars. Based on our color-color census for the 65 protostar candidates detected by WISE, the Class II / Class I ratio is equal to 0.23. Following Dunham et al. (2015), we assumed a Class II duration of 2 Myr and used the standard method to infer the age of a stellar population from this ratio (Wilking et al. 1989, Evans et al. 2009). We obtain an age of ∼500 kyr for the YSO population found in the extended G region, which is not consistent with the scenario in which the formation of these YSOs was triggered by interaction of the SNR with its environment. It is likely that the collapse of molecular clouds in the vicinity of the high-mass progenitor was triggered by the compression driven by stellar winds (Xu et al. 2011). To find evidence of enhanced star formation in the SNR shocks would require to detect deeply embedded prestellar cores, which is not attainable with the data at our disposal. In the future, high-resolution (∼1 ) sub-mm and mm observations will be required to disentangle the kinematics of large-scale shock structure and outflows, and to eventually uncover prestellar cores in the extended G region (e.g. as Motte et al. (2018) did for the W43-MM1 star forming region). Still, our sample of protostar candidates provides a unique opportunity to study star formation in an environment that is subject to intense γ-ray flux and shock dynamical feedback. The protons accelerated by young protostars might also provide a source of fresh mid-energy CRs (up to 10-13 GeV and 26-37 GeV respectively for jet acceleration and protostellar surface shock acceleration in lower-mass stars, Padovani et al. 2016). With jet velocities up to 1000 km s −1 , high-mass protostellar shocks could contribute to the TeV γ-ray peak in the extended G region.
Point sources identified as star-forming galaxies from their location in the color-color diagram display a significant correlation with the spiral-shaped, bright and extended infrared features mapped by Spitzer-MIPS at 24 µm. Most of these sources might as well be identified as PAH emission, as our arbitrary criterion of selection puts a significant uncertainty on the cut between the two categories of sources, in particular along the W 1 -W 2 axis (Koenig et al. 2012). In fact, it is more likely that these point sources are associated with PAH, stellar or protostellar emission. With our current criterion, only 2 sources identified as PAH emission are found and located in the same area as point sources that were labelled 'star-forming galaxies'. If most of these detections are actually arising from PAH emission, this result would support the assumption that the bright infrared extended emission arise from a dark dust lane that is spatially associated with the molecular cloudlet detected by our CO observations.
In comparison to previous work on the star formation in IC443: -Based on 2MASS point source catalogues, Xu et al. (2011) found a total of 1666 YSO candidates, 154 CTTS and 419 HAeBe stars in a search circle around IC443 within a 25 radius. Their candidates were mostly concentrated around the CO molecular shell, in particular towards i. the clump C; ii. the extended G region where their YSO candidates are spatially correlated with the shocked clump. They proposed that the formation of these YSOs have been triggered by the stellar wind of the IC443 progenitor. -Based on 2MASS and WISE point source catalogues, Su et al. (2014) found a total of 98 YSO candidates. They proposed a sample of 62 YSO candidates concentrated along the boundary of the radio shell. They also proposed that the formation of these protostars is likely to have been triggered by the stellar winds of the SNR progenitor. In contrast with our results and the findings of Xu et al. (2011), their distribution does not show a strong correlation with the molecular structures in the extended G region.

Summary
transitions of both 12 CO and 13 CO, as well as the two first rotational transitions of C 18 O obtained with the IRAM-30m and APEX telescope. These extensive maps allow to probe the position of the TeV γ-ray peak in IC443G and its surroundings with unprecedented spectral and angular resolution at millimeter and sub-millimeter wavelengths. 2. We proposed a description of the morphology and kinematics of the extended G region based on the definition of four main molecular structures: (i.) a shocked clump; (ii.) a quiescent, "ring-like" structure that might be spurious; (iii.) a quiescent cloudlet that is spatially connected to the ring-like structure and (iv.) a shocked knot. 3. The comparison of our data histograms with modified LTE radiative transfer models revealed enhanced 12 CO J=2-1/1-0 and J=3-2/2-1 ratios towards the shocked clump, in particular within the high-velocity wings of the emission lines that are tracing warm CO (40 -120 K). We provide a rough estimate of the isotopic ratio based on the best match between modified LTE models and the measured intensity of the 12 CO, 13 CO and C 18 O J=1-0 line. 4. We measured the mass of the molecular gas in the extended G region using pixel-per-pixel, channel-per-channel population diagrams corrected for optical depth. We obtained a total molecular mass estimate of ∼230 M , ∼90 M , ∼210 M and ∼4 M respectively for the cloudlet, ring-like structure, shocked clump and shocked knot. We also measured a mass of ∼1100 M for the ambient gas in the [-6.5, -1.5] km s −1 LSR range. The estimate of the mass depends on the adopted 12 CO/ 13 CO isotopic ratio, but it is established that a total molecular mass 0.9-3.1×10 3 M is available to interact with CRs via pion-decay in the extended G region. 5. We proposed a second estimate of the mass using the LVG assumption with a grid of RADEX models. The χ 2minimization of the grid of models with respect to data in each structure yields results that are in partial agreement with the previous measurements, with an estimate that is systematically lower than the LTE estimate, except for the mass measurements of the high-velocity wings towards the shocked clump, which are higher by a factor ∼3 . 6. We studied the spectral energy distribution of infrared and optical point sources using the 2MASS, WISE and Gaia catalogues of point sources. Using color-color diagrams we determined a sample of protostar candidates in the extended G region. 144 protostar candidates are found in the region, and 16 stars are found in the Gaia census but it is worth to note that a large fraction of these point sources are contaminated by extended emission (12.7% for 2MASS, and 76.9% for WISE). Based on the spatial distribution of these candidates, we propose three sites of possible star formation in IC443G (groups 1 and 2 in Fig. 14, and the shocked clump). These protostars might constitute a source of fresh CRs.
= 8π 3 µ 2 3hc If we introduce the partition function we get a relation linking n J to the total density n: Where E J and E J+1 are the level energies. So that: Thus the total column density N can be obtained by computing the following integral over the line of sight: Using the definition of the equivalent width W = τdv and the variable dν = ν J dv/c, we finally get the relation between column density and equivalent width for a given excitation temperature: Where Λ = 3h/8π 3 µ 2 , h is the planck constant and ν i j the frequency of the transition. Applied to the transitions J=1-0 and J=2-1, this relation can be written for two different values of the equivalent width W corresponding to each line. As the total column density is constant in each equation, the ratio of these two relations yields: If we define the temperature of the transition T 0 = hν i j /k and use the variable u = T 0 /T ex then we have: and, finally we get the following relation which permits to determine the temperature from a measure of the ratio of the two equivalent widths: Where F(u) = e −3u cosh(u).
Article number, page 28 of 36 P. Dell'Ova et al.: Interstellar anatomy of the TeV gamma-ray peak in the IC443 supernova remnant      Fig. 6, for the following lines: 12 CO (in black), 13 CO (in blue) and C 18 O (in red). First row presents the (1-0) transition, second row presents the (2-1) transition and third row presents the (3-2) transition for: A, the cloudlet; B, the ring-like structure; C, the shocked clump and D, the shocked knot. Spectral cubes were resampled to allow direct comparison between the different spectra. Spatial resolutions of all transitions were modified to the nominal resolution of C 18 O(2-1), θ = 30.2 . Spectral resolution were set to 0.5 km s −1 for IRAM-30m data, and 0.25 km s −1 for APEX data. On both panels, the LSR of IC443 is indicated with a vertical dashed line (at −4.5 km s −1 ). The black data points with error bars represent the measurements of the optical depth τ corresponding to each channel of these spectra resampled with a spectral resolution of 2 km s −1 , using the [ 12 CO]/[ 13 CO] ratio as a tracer, and red data points represent the measurements of the optical depth obtained using the [ 12 CO]/[C 18 O] ratio. The adopted values for the isotopic ratios are based on the results of section 4.2.1 and indicated in Tab. 6.