2.1. 7-mm Mopra data and analysis

We conducted 20′ × 20′ and 10′ × 10′ observations of the ISM towards several HESS TeV sources as part of our ‘MopraGam’ survey.Footnote ^a The 7-mm observations towards these studied sources were carried out between 2012 and 2015. Table 1 indicates the position and size of the different observations towards the seven regions that were mapped.

Table 1. Mopra 7-mm coverage of the TeV sources studied in this paper. The central position in (RA, Dec) and size is given for all observations undertaken towards the TeV sources.

We used the Mopra spectrometer MOPS in ‘zoom’ mode, which records 16 sub-bands, each consisting of 4 096 channels over a 137.5-MHz bandwidth. The ‘on the fly mapping’ (OTF) Nyquist sampled these regions with a 1′ beam size and a velocity resolution of ∼ 0.2 km/s. We can thus simultaneously observe tracers that can help understand the structure and morphology of MCs. Among the tracers observed in this paper, both the carbon monosulfide and cyanoacetylene transitions CS(1–0) and HC₃N(5 − 4, F = 4 − 3) can be found in dense molecular gas and star-forming regions (Irvine, Goldsmith, & Hjalmarson Reference Irvine, Goldsmith, Hjalmarson, Hollenbach and Thronson1987). Silicon monoxide SiO emission can also be detected inside shocked dense MCs (Schilke et al. Reference Schilke, Walmsley, Pineau des Forets and Flower1997; Gusdorf et al. Reference Gusdorf, Cabrit, Flower and Pineau Des Forêts2008) and is a good signpost to claim physical association between an SNR and a nearby MC (e.g. see Nicholas et al. Reference Nicholas, Rowell, Burton, Walsh, Fukui, Kawamura and Maxted2012). Finally, the 44-GHz methanol class I maser CH₃OH(I) generally indicates nearby active star-forming regions (Voronkov et al. Reference Voronkov, Caswell, Ellingsen, Green and Breen2014).

2.2. CO(1–0) data and analysis

While the Mopra CO(1–0) survey is well under way (see Braiding et al. Reference Braiding2015, Reference Braiding2018), we used for this work the 4-m Nanten CO(1–0) survey (Mizuno & Fukui Reference Mizuno, Fukui, Clemens, Shah and Brainerd2004), as it encompassed all our sources. The Nanten telescope CO(1–0) survey covered the Galactic plane with a 4′ sampling grid, a velocity resolution Δv = 0.625 km/s, and an averaged noise temperature per channel of ∼ 0.4 K.

The CO(1–0) emission probes the more diffuse molecular gas surrounding PWNe. As PWNe expand inside a low-density medium, we seek velocity ranges where the CO(1–0) emission spatially anti-corresponds with the TeV emission as that would support the PWN scenario (see Section 4.2). Molecular gas may also provide sufficient target material to produce hadronic TeV emission. Therefore, probing extended CO(1–0) emission overlapping the TeV emission is also a powerful means to test the hadronic scenario in the vicinity of a CR accelerator.

We arbitrarily selected CO regions based on the prominence of the CO(1–0) integrated emission compared to the rest of the maps and their proximity with the TeV gamma-ray source, as these can be helpful to test the hadronic/leptonic scenario (see Section 4). In the case where the pulsar is offset from the observed TeV gamma-ray source (i.e. HESS J1026–583), we highlight molecular regions which could produce such asymmetry along the pulsar-MC axis. We also derive physical parameters towards some extended regions defined from our CS detections (e.g. see HESS J1809–193 and HESS J1418–609). Assuming that all the gas traced by CO is embedded within this CS region, we are able to derive an upper limit on the averaged density $n_{\textrm{H}_2}$ .

We fitted the CO components at the velocity range of interest, with Gaussian distributions. We then used the X-factor X _CO = 2 × 10²⁰ cm⁻²/(K km/s) to convert the CO integrated intensity into H₂ column density. Bolatto, Wolfire, and Leroy (Reference Bolatto, Wolfire and Leroy2013) have argued that this value is correct to ∼ 30% across the Galactic plane. We then used Eq. D.4 to obtain the total molecular mass $M_{\textrm{H}_2}\left(\textrm{CO}\right)$ of the cloud, accounting for a 20% contribution from helium. We finally assumed a prolate geometry to obtain the H₂ volume occupied by the molecular gas, and the averaged particle density $n_{\textrm{H}_2}\left(\textrm{CO}\right)$ using Eq. D.5 (see Appendix D). The Gaussian fits for the CO(1–0) components towards the regions of the individual TeV source can be found in Tables C.1–C.6 and their derived physical properties in Tables D.1–D.6.

2.3. HI analysis with SGPS/GASS surveys

From the Southern Galactic Plane Survey (SGPS) and GASS Hi surveys (McClure-Griffiths et al. Reference McClure-Griffiths, Dickey, Gaensler, Green, Haverkorn and Strasser2005; McClure-Griffiths et al. Reference McClure-Griffiths2009), we also obtained the atomic Hi column density N _Hi assuming the (optically thin) conversion factor X _Hi = 1.8 × 10¹⁸ cm⁻²/(K km/s) (Dickey & Lockman Reference Dickey and Lockman1990). Hi self-absorption may however occur towards colder Hi regions and, as a result, the optically thin assumption may underestimate N _Hi by a factor of 2 (Fukui et al. Reference Fukui, Torii, Onishi, Yamamoto, Okamoto, Hayakawa, Tachihara and Sano2015). Comparing the total column density $N_{\textrm{H}}\,{=}\,N_{\textrm{H}\textsc{i}}+2N_{\textrm{H}_2}$ to the absorbed X-ray column density from X-ray counterparts can provide further constraints on the distance to the TeV source. The images showcasing the evolution of the total column density as a function of v _lsr and distance towards several TeV sources can be found in Appendix E. Atomic gas should also be accounted for while testing the hadronic scenario towards TeV sources. Consequently, as per the CO(1–0) emission, we derived the atomic masses $M_{\textrm{H}_I}$ towards the selected CO regions. These masses can be found in Tables D.1–D.6.

3.1. HESS J1809-193

HESS J1809–193 is a bright and extended TeV source whose position is coincident with several potential CR accelerators (Aharonian et al. Reference Aharonian2007). HAWC has also recently detected the TeV source 2HWC J1809–190, associated with HESS J1809–193 (Abeysekara et al. Reference Abeysekara2017). Araya (Reference Araya2018) recently found extended GeV emission associated with HESS J1809–193.

ASCA (Bamba et al. Reference Bamba, Ueno, Koyama and Yamauchi2003) and Suzaku observations (Anada et al. Reference Anada, Bamba, Ebisawa and Dotani2010) revealed non-thermal X-rays likely associated with the pulsar PSR J1809–1917 (shown as a cyan diamond in Figure 1 and pink in Figure 2), with spin down power $\dot{E}_{\textrm{SD}}=1.8\times10^{36}\,{\rm {erg\,s^{-1}}}$ , a characteristic age τ _c = 51 kyr and a dispersion measure distance d ∼ 3.7 kpc (Cordes et al. Reference Cordes, Lazio, Chatterjee, Arzoumanian and Chernoff2002). The presence of two SNR shells G011.0–0.0 and G011.1+0.1 (shown as blue circles in Figure 1) both observed at 330 and 1 465 MHz (Brogan et al. Reference Brogan, Devine, Lazio, Kassim, Tam, Brisken, Dyer and Roberts2004; Castelletti, Giacani, & Petriella Reference Castelletti, Giacani and Petriella2016) adds more complexity to the picture. Additionally, the ∼2 kyr old millisecond pulsar PSR J1811–1925 with spin down energy $\dot{E}_{\textrm{SD}}=6.4\times 10^{36} \,{\rm {erg\,s^{-1}}}$ (Torii et al. Reference Torii, Tsunemi, Dotani, Mitsuda, Kawai, Kinugasa, Saito and Shibata1999) and its progenitor SNR G011.2–0.3, located at d ∼ 4.4 to 7 kpc, are also positioned adjacent to HESS J1809–193 (see Figures 1 and 2) and thus might also contribute to the TeV emission. Finally, the anomalous X-ray magnetar XTE J1810–197 with period P = 5.54 s and period derivative $\dot{P}=1.8\times10^{-12}\,{\rm {s\,s^{-1}}}$ is located ∼ 0.35° south of HESS J1809–193.

Figure 2. Three colours image showing the MIPSGAL 24 μm and GLIMPSE 8 μm and 4.6 μm in red, green, and blue, respectively, towards HESS J1809–193 overlaid by the HESS TeV gamma-ray counts in dashed white contours and CS(1–0) integrated intensity between v _lsr = 25 to 38 km/s contours (0.6 K) in solid white. The SNRs are shown in red dashed circles, while the pulsars position are indicated in pink diamonds. The black dashed ellipses labelled ‘HC’ indicates the positions of HC₃N(5–4,F = 4–3) detections, while the black solid circles labelled ‘CH’ and ‘S’, respectively, indicate CH₃OH and SiO(1–0,v = 0) detections. The spectra of these regions can be found in Figure 3. The white dotted lines represent the extent of the infrared dark clouds IRDC G010.71–00.16a–h and IRDC G011.11–00.11a–e. (see colour version online)

Based on the kinematic distance of PSR J1809–1917 (d ∼ 3.7 kpc), we here focus on components at v _lsr = 10 − 38 km/s (see blue and pink regions in Figure 3), although our Nanten CO(1–0) data have revealed several other components in the line of sight (see CO(1–0) emission in Figure 3).

Figure 3. Averaged CS(1–0), CO(1–0), SiO(1–0, v = 0), CH₃OH(I), and HC₃N(5–4, F = 4–3) spectra towards the emission from the selected regions in Figures 1 and 2 towards HESS J1809–193. The solid red lines represent the Gaussian fit of the emission whose parameters are shown in Table C.1. The two red vertical lines indicate the pulsar PSR J1809–1917 dispersion measure distance converted to kinematic velocity. The pink and cyan regions represent the velocity range for the CS(1–0) and CO(1–0) integrated intensity maps displayed in Figure 1.

At this velocity range, we remark that the molecular gas generally overlaps the TeV emission shown in black contours. Notably, at v _lsr = 25 to 38 km/s (d ∼ 3.7 kpc), prominent CO emission is found south and east of the TeV emission, while the prominent CO emission spatially overlaps the TeV emission at v _lsr = 10 − 25 km/s (d ∼ 2.7 kpc). Our Mopra CS(1–0) integrated intensity maps (see Figure 1 panels b, d, and f) however provide a clear view of the dense gas.

At v _lsr = 25 to 38 km/s (Figure 1, panels a and b), the MCs in the regions labelled ‘1’ and ‘2’ and located east and south of SNR G011.0+00.0, respectively, appear very extended. From our CO and CS analyses, the masses derived in region ‘1’ attain $M_{\textrm{H}_2}\left(\textrm{CO}\right)=8.1\times10^{4}{\textrm{M}_{\odot}}$ and $M_{\textrm{H}_2}\left(\textrm{CS}\right)=3.2\times10^{4}{\textrm{M}_{\odot}}$ , while we obtain $M_{\textrm{H}_2}\left(\textrm{CO}\right)=2.3\times10^{5}{\textrm{M}_{\odot}}$ and $M_{\textrm{H}_2}\left(\textrm{CS}\right)=3.2\times10^{4}{\textrm{M}_{\odot}}$ in region ‘2’. A significant fraction of the molecular gas in region ‘1’ and ‘2’ is therefore concentrated in clumps. The CS(1–0) emission in region ‘1’ also appears to anti-correspond with the ASCA X-ray diffuse emission shown in blue, supposedly produced by high energy electrons from the pulsar PSR J1809–1317 (see Anada et al. Reference Anada, Bamba, Ebisawa and Dotani2010).

Interestingly, we have found embedded dense filaments in region ‘1’ from C³⁴S(1–0) and HC₃N(5–4, F = 4–3) detections (see Figure F.1 in Appendix F). In fact, the positions of the HC₃N detections labelled ‘HC1 to HC3’ in Figure 2 coincide with the Spitzer infrared (IR) dark cloud IRDC G011.11–00.11a–e (Parsons, Thompson, & Chrysostomou Reference Parsons, Thompson and Chrysostomou2009, see dotted white lines in Figure 2), confirming the molecular gas is foreground to the IR emission.

We have also identified two weak but broad (FWHM ∼4 km/s) SiO(1–0) features in region ‘1’ (see Figures 2 and 3 for spectra), labelled ‘S1’ and ‘S2’. The absence of overlapping IR continuum emission indicates the lack of active star-forming regions and could indicate a possible interaction between the MC and a non-star-forming shock (Schilke et al. Reference Schilke, Walmsley, Pineau des Forets and Flower1997; Gusdorf et al. Reference Gusdorf, Cabrit, Flower and Pineau Des Forêts2008), which could here come from the adjacent SNR G011.0–0.0.

Based on the potential interaction between the SNR G011.0–0.0 and the dense MC at v _lsr ∼ 30 km/s, we thus suggest an alternate SNR distance d ∼ 3.7 kpc compared to the distance d ∼ 3.0 kpc claimed by Castelletti et al. (Reference Castelletti, Giacani and Petriella2016). If SNR G011.0–0.0 indeed interacted with the MC at d ∼ 3.7 kpc, it would then be located at the pulsar PSR J1809–1917 distance and quite likely be its progenitor SNR. Using eq. 3.33a from Cioffi, McKee, and Bertschinger (Reference Cioffi, McKee and Bertschinger1988), we derive that the ambient density required to match the small projected radius r_SNR ∼ 5 pc with the pulsar characteristic age τ _c = 51 kyr is n _amb ∼ 370 cm⁻³. The averaged density $n_{\textrm{H}_2}\left(\textrm{CO}\right)=440\,{\rm{cm}}^{-3}$ found towards region ‘1’ appears somewhat consistent with this estimated averaged density.

In region ‘2’, we have also detected extended HC₃N(5–4, F = 4–3) emission and C³⁴S(1–0) emission overlapping the IR dark cloud IRDC G010.71–00.16a–h (Parsons et al. Reference Parsons, Thompson and Chrysostomou2009), highlighting another dense region. The morphology of the IR dark cloud and our HC₃N detection appear elliptic and embed the anomalous X-ray magnetar XTE J1810–197 (see Gotthelf et al. Reference Gotthelf, Halpern, Buxton and Bailyn2004 and references therein), suggesting their potential physical association. Additionally, a prominent SiO(1–0) detection (see region ‘S3’ in Figure 2) was also found inside this dark cloud. This MC may thus be disrupted by another shock, perhaps caused by the progenitor SNR of XTE J1810–197.

At v _lsr = 10 − 25 km/s (Figure 1 panel d), we observe several CS(1–0) components inside the region here labelled ‘3’. We notably find that the CS(1–0) prominent emission corresponds with the JCMT CO(3–2) peaks found by Castelletti et al. (Reference Castelletti, Giacani and Petriella2016) and the Hii regions (purple contours and blue crosses in Figure 1 panel f, respectively).

We observe that the CS(1–0) emission averaged over the grid regions (Figure 1 panels e–f) exhibits considerable variation of the peak velocity ranging between v _lsr = 10 − 22 km/s. For example, the two peaks at v _lsr ∼ 12 and 18 km/s inside ‘boxes 9 and 10’ merge to a single peaked emission with v _lsr ∼ 15 km/s in ‘box 15’. This dense molecular region appears to host several Hii regions (see ‘H1–H3’ in Figure 1f). Consequently, the two spectral components may actually probe the same molecular gas, disrupted by the driving motion forces from various Hii regions. We also note that the MC anti-corresponds with the two SNRs. The disrupted gas could consequently be caused by one of these SNRs. As a result, we cannot rule out the SNR G011.0–0.0 distance d ∼3.0 kpc suggested by Castelletti et al. (Reference Castelletti, Giacani and Petriella2016).

3.2. HESS J1026–583

The TeV source HESS J1026–583 was discovered from energy dependent morphology studies towards HESS J1023–591 (Abramowski et al. Reference Abramowski2011). The latter source is thought to be powered by the colliding winds from Wolf–Rayet stars within the massive stellar cluster Westerlund 2 at $d=5.4^{+1.1}_{-1.4}$ kpc (Furukawa et al. Reference Furukawa, Dawson, Ohama, Kawamura, Mizuno, Onishi and Fukui2009). Acero et al. (Reference Acero2013) catalogued HESS J1026–583 as a PWN candidate based on the detection of a nearby radio quiet gamma-ray pulsar PSR J1028–5810 (Ray et al. Reference Ray2011) responsible for the GeV emission towards 3FGL J1028.4–5819 (shown as a red circle in Figure 4). PSR J1028–5819 has a spin down power $\dot{E}_{\textrm{SD}}=8.3\times10^{35}\,{\rm {erg\,s^{-1}}}$ , characteristic age τ = 89 kyr, and a dispersion measure distance d = 2.3 ± 0.7 kpc. However, HESS J1026–583 shows a hard VHE spectral index Γ _γ = 1.94. It also does not exhibit any X-rays that are spatially coincident with the TeV source. Additionally, diffuse GeV gamma-ray emission has been detected with Fermi-LAT towards Westerlund 2 (Yang, de Oña Wilhelmi, & Aharonian Reference Yang, de Oña Wilhelmi and Aharonian2017). The authors have argued for a hadronic origin based on its 200 pc extension of the 1–250 GeV emission. Consequently, a clear identification remains to be seen. Because of its proximity to HESS J1023–575, several ISM features have already extensively been studied (Dame Reference Dame2007; Fukui et al. Reference Fukui2009; Furukawa et al. Reference Furukawa, Dawson, Ohama, Kawamura, Mizuno, Onishi and Fukui2009; Furukawa et al. Reference Furukawa2014; Hawkes et al. Reference Hawkes2014).

Figure 4. Nanten CO(1–0) and Mopra CS(1–0) emission between v _lsr = −23 to 13 km/s and v _lsr = 0 to 20 km/s towards HESS J1026–582 and HESS J1023–575 whose TeV gamma-ray counts are shown in black contours. The position of the pulsars PSR J1028–5819, PSR J1023–5746, and PSR J1019–5749 are indicated as cyan diamonds. The GeV emission 3FGL J1028–5819 is shown as a red circle. The cyan ellipses indicate the selected regions (labelled A and B) from our CO analysis, while the black circles (labelled 1 to 7 in panels b and d) show the position of selected CS(1–0) regions. The location of WR 21a is shown as a blue cross in panel d, while the purple square and cross indicate the position of the Hii region GAL 284.65–00.48 and the reflection nebula GN 10.23.6, respectively. Panel e is a two-colour image showing the Nanten CO(1–0) integrated intensity at v _lsr = 2 to 6 km/s (red) and v _lsr = 6 to 16 km/s (green) overlaid by the HESS TeV contours in white towards HESS J1026–582. The cyan dashed ellipse represents the possible molecular ring structure (see Section 3.2). Panel f illustrates the averaged CO(1–0) and CS(1–0) emission from the selected regions. The red lines indicate the fit used to model the emission (see Table C.2 for fit parameters). The two red vertical lines show the dispersion measure distance of the pulsar PSR J1018–5819. The blue and pink regions indicate the velocity range shown in panels a to d.

Figure 4 shows the Nanten CO(1–0) integrated intensity at v _lsr = −23 to −13 km/s (panel a) and 0 to 20 km/s (panel c), inferring a distance d ∼ 2.3 kpc and d ∼ 4.9 kpc, respectively, positioning the detections along the Carina arm (see Figure A.3 in Appendix B). At v _lsr = −23 to −13 km/s, we note that the CO emission does not overlap the TeV emission, nor the pulsar’s position. In fact, the molecular region with prominent CO emission, which we labelled ‘A’, is located east of the pulsar position, and could support the crushed PWN scenario (see Blondin et al. Reference Blondin, Chevalier and Frierson2001). Assuming a kinematic distance d ∼ 2.3 kpc, we obtain a total mass $M_{\textrm{H}_2}\left(\textrm{CO}\right)=9.7\times10^{3}{\textrm{M}_{\odot}}$ . Lastly, from our CS survey, we have also identified two unresolved CS(1–0) features labelled ‘1’ and ‘2’ in Figure 4 panel b.

At v _lsr = 0 to 20 km/s (d ∼ 4.9 kpc, see Figure 4 panel c), we observe prominent CO emission at v _lsr ∼ 4 km/s, labelled ‘B’, spatially coincident with the HESS J1026–582 TeV peak. The MC in region ‘B’ with total mass $M_{\textrm{H}_2}\left(\textrm{CO}\right)=4.9\times10^{4}{\textrm{M}_{\odot}}$ appears next to a possible shell like structure (see dashed purple and cyan ellipses in Figure 4 panels c and e, respectively) which overlaps the HESS J1023–575 TeV emission. The molecular region is however located at the tangent of the Sagittarius arm (see Figure A.3); thus, it is possible that the various MCs found in Figure 4 panel c may be unrelated. From Figure 4 panel d, it nonetheless appears that the CO(1–0) emission at v _lsr = 2 to 6 km/s fills in some missing segment of the shell-like structure and might form a ring (indicated as a dashed-cyan ellipse) whose centre is located south-west of HESS J1026–582 observed at v _lsr = 6 to 16 km/s. The molecular structure at v _lsr = 6 to 16 km/s could therefore be physically connected to the MC in region ‘B’.

Among the CS(1–0) features found at v _lsr = 0 to 20 km/s, the gas clumps labelled ‘5’ to ‘7’, embedded in the molecular gas in region ‘B’ form a partial shell structure overlapping the TeV peak emission with a combined mass $M_{\textrm{H}_2}\left(\textrm{CO}\right)\sim5.0\times10^{3}{\textrm{M}_{\odot}}$ , very similar to the CO mass obtained in region ‘B’. No massive stars have been catalogued towards the centre of this structure. However, we note the presence of the Hii region GAL 284.65–00.48 towards region ‘5’ and the reflective nebula GV 10.23.6 at the centre of the shell (see purple box and cross in Figure 4 panel d, respectively).

3.3. HESS J1119–614

Djannati-Ataï (Reference Djannati-Ataï, Marandon and Chaves2009) and Abdalla et al. (Reference Abdalla2018a) reported the detection of the TeV source HESS J1119–614 (see solid black contours in Figure 5). It is thought to be associated with either the PWN, powered by the radio pulsar PSR J1119–6127 (see cyan diamond in Figure 5), with spin down period P = 400 ms (Camilo et al. Reference Camilo, Kaspi, Lyne, Manchester, Bell, D’Amico, McKay and Crawford2000) and characteristic age τ _c = 1.9 kyr (Weltevrede, Johnston, & Espinoza Reference Weltevrede, Johnston and Espinoza2011); or its progenitor SNR G292.2–0.5 (see dashed blue circle in Figure 5). The 3″ × 3″ PWN has been resolved in X-ray with Chandra (Gonzalez & Safi-Harb Reference Gonzalez and Safi-Harb2005; Safi-Harb & Kumar Reference Safi-Harb and Kumar2008). The shell of progenitor SNR has been observed both in the radio band with ATCA (Crawford et al. Reference Crawford, Gaensler, Kaspi, Manchester, Camilo, Lyne and Pivovaroff2001), and in X-ray with ASCA between 0.4 and 10 keV (Pivovaroff et al. Reference Pivovaroff, Kaspi, Camilo, Gaensler and Crawford2001). Follow-up studies with Chandra and XMM-Newton (Kumar et al. Reference Kumar, Safi-Harb and Gonzalez2012; Ng et al. Reference Ng, Kaspi, Ho, Weltevrede, Bogdanov, Shannon and Gonzalez2012) revealed additional information about the nature of the X-ray emission. Indeed, from the four regions studied by Kumar et al. (Reference Kumar, Safi-Harb and Gonzalez2012) (here labelled ‘X1’ to ‘X4’ in Figure 5), the three regions ‘X1’ to ‘X3’ shown as blue ellipses have a non-thermal component. However, the X-ray emission in region ‘X4’ shown as a blue-pink ellipse is prominent and is thought to be of thermal origin.

Figure 5. Nanten CO(1–0) emission between v _lsr = −20 to −10 km/s (panel a) and v _lsr = 20 to 40 km/s (panel b) towards HESS J1119–614 whose TeV gamma-ray emission is shown as solid black contours. The progenitor SNR G292.2–0.5 is delimited by the blue dashed circle. The solid blue and dashed blue-pink ellipses labelled ‘X1’ to ‘X4’ (see text) highlight bright X-ray regions studied by Kumar, Safi-Harb, & Gonzalez (Reference Kumar, Safi-Harb and Gonzalez2012) with XMM-Newton and Chandra. The pulsar PSR J1119–6127 (P1)’s position is indicated as a cyan diamond. The cyan ellipses show the selected regions (labelled A to D) for our CO analysis. Panels c and d show the variation of the averaged CO(1–0) spectra over the black grid of boxes shown in panel (b). The cyan and pink regions indicate the velocity ranges mapped in panels a and b. Panel e shows the averaged CO(1–0) emission from the selected regions towards HESS J1119–164. The red lines indicate the fit used to parametrise the emission. The fit parameters are displayed in Table C.5. The pink and cyan regions show the velocity range used for the above integrated intensity maps.

Caswell, McClure-Griffiths, & Cheung (Reference Caswell, McClure-Griffiths and Cheung2004) estimated a distance d ∼ 8.4 kpc for the pulsar and its progenitor SNR, based on the Hi and magnetic field studies. However, reconciling the characteristic age τ _c = 1.9kyr of PSR J1119–6127 and the ∼ 25 pc diameter of the progenitor SNR requires a low-density medium (Ng et al. Reference Ng, Kaspi, Ho, Weltevrede, Bogdanov, Shannon and Gonzalez2012) and a massive progenitor star. Gonzalez & SafiHarb (Reference Gonzalez and Safi-Harb2005) argued for an SNR distance at d = 3.6 to 6.3 kpc by modelling the X-ray spectrum. Our gas study thus aims to provide additional constraints on the distance and how it could affect the nature of the TeV source.

No extended CS emission have been found towards this region. From our Nanten CO data, several components have been detected along the line of sight at v _lsr ∼ −30 km/s, v _lsr ∼ −10 km/s (near/far distance d ∼ 2.6/4.5 kpc, near Carina-arm, see Figure A.3), v _lsr ∼ 20 km/s (distance d ∼ 8.6 kpc, far Carina-arm), v _lsr ∼ 30 km/s (kinematic distance d ∼ 9.7 kpc, far Carina-arm).

At v _lsr = −20 to −10 km/s, we observe two MCs (labelled ‘A’ and ‘B’ in Figure 5 panel a) positioned north-east and west to the SNR, respectively, and with respective masses $M_{\textrm{H}_2}\left(\textrm{CO}\right)=2.3\times10^{4}{\textrm{M}_{\odot}}$ and $M_{\textrm{H}_2}\left(\textrm{CO}\right)=7.1\times10^{4}{\textrm{M}_{\odot}}$ . We also note that the molecular ISM anti-corresponds with all X-ray regions.

At v _lsr ∼ 20 to 40 km/s we observe that the bulk of the molecular gas is found in two regions labelled ‘C’ and ‘D’, with masses $M_{\textrm{H}_2}\left(\textrm{CO}\right)=1.3\times10^{4}{\textrm{M}_{\odot}}$ and $M_{\textrm{H}_2}\left(\textrm{CO}\right)=2.3\times10^{5}{\textrm{M}_{\odot}}$ , respectively. The morphology of the gas overlaps the TeV gamma-ray detection. It also corresponds with the thermal X-ray in region ‘X1’ both at v _lsr ∼ 20 km/s and v _lsr ∼ 30 km/s (see Figure 5 panel d), highlighting potential SNR-MC interaction. The molecular gas also anti-corresponds with the other X-ray regions ‘X2’ to ‘X4’, thought to have non-thermal components likely produced by the high energy electrons from both the SNR and the PWN. Therefore, we argue that the morphology of the CO(1–0) both at v _lsr ∼ 20 km/s and v _lsr ∼ 30 km/s appears consistent with the X-ray results, inferring a source kinematic distance d = 8.6 − 9.7 kpc. We finally remark that these two components may be physically connected and highlight motion caused by the progenitor star.

We also compare the column densities, shown in Figure E.1, towards the four X-ray regions. Towards the regions ‘X1’ and ‘X2’, we match the X-ray modelled column density N _h (grey regions in Figure E.1 top panels) at v _lsr = 22 − 30 km/s, also suggesting a distance d = 8.6 − 9.6 kpc, positioning the source in the far-Carina arm. Towards the X-ray regions ‘X3’ and ‘X4’ however, our column densities do not reach the X-ray one, or yield to unreasonable large distance d ∼ 15.5 kpc. Clumps unresolved by the Nanten CO(1–0) may account for the missing components. Our column density studies thus provide a lower-limit kinematic distance d > 8.6 kpc, consistent with our morphological studies and with Caswell et al. (Reference Caswell, McClure-Griffiths and Cheung2004). To explain both the large SNR radius and the enhanced thermal X-ray emission in region ‘X4’, we argue that the SNR once propagated in a low-density medium (Ng et al. Reference Ng, Kaspi, Ho, Weltevrede, Bogdanov, Shannon and Gonzalez2012), until it reached the denser gas in region ‘D’.

3.4. Kookaburra and rabbit

The two TeV sources HESS J1418–609 and HESS J1420–607 have been classified as PWNe based on their spatial coincidence with the X-ray (Roberts & Romani Reference Roberts and Romani1998) and GeV gamma-ray counterparts (Acero et al. Reference Acero2013). Ng, Roberts, & Romani (Reference Ng, Roberts and Romani2005) indicated that two diffuse non-thermal X-ray sources were associated with the pulsar PSR J1420–6048 (labelled P1 in Figure 6), a 68.2-ms period pulsar with spin down power $\dot{E}_{\textrm{SD}}=1.0\times10^{37}\,{\rm {erg\,s^{-1}}}$ , characteristic age τ _c = 13 kyr and dispersion measure distance d ∼ 5.6 kpc; and the 108-ms radio-quiet gamma-ray pulsar PSR J1418–6058 (labelled P2), with a characteristic age τ _c = 1.6 kyr. It has also been noted that P2’s location is offset (∼ 8.4′) from the HESS J1418 − 09 TeV peak.

Figure 6. Nanten CO(1–0) and Mopra CS(1–0) emission between v _lsr = −65 to −55 km/s, −52 to −42 km/s, and −10 to 0 km/s towards HESS J1420–607 and HESS J1418–609 shown in black contours. The black dashed box indicates our Mopra 7-mm coverage. The position of the pulsars PSR J1420–6048 and PSR J1418–6058 (labelled P1 and P2 here) are shown as cyan diamonds while the nearest SNR G312.4–0.04 is indicated as a blue dashed circle. Our CO regions ‘A’ and ‘B’ are shown in cyan in panels c and e. The blue ellipses on panel c and the black ellipses on the right panels indicate the position of the CS regions (‘1’ to ‘15’). The position of the CH₃OH(I) detections labelled ‘CH1’ to ‘CH3’ are shown in green circles.

Roberts et al. (Reference Roberts, Brogan, Gaensler, Hessels, Ng and Romani2005) argue that the expected cometary shape of the PWN produced by a pulsar with a high space velocity would explain the offset position of the gamma-ray emission with respect to the pulsar. The distance to P2 has not yet been constrained. Ng et al. (Reference Ng, Roberts and Romani2005) have argued a PWN distance d = 2 to 5 kpc, while Wang (Reference Wang2011) has claimed a much smaller distance d = 1.4 − 1.9 kpc. Thus, by looking at the gas distribution in various velocity ranges, we aim to highlight the distance which would support the PWN scenario.

From our CO(1–0) observations, we have detected several molecular complexes along the line of sight. We mostly focus at the velocity ranges v _lsr = −65 to −55 km/s, (d ∼ 5.6 kpc), ∼50 km/s (Scutum Crux arm, distance d ∼ 3.5 kpc), and ∼ −5 km/s (local arm, d ∼ 0.1 kpc), respectively, as shown in Figure 6. In Figure 7, we also note a CO component at v _lsr ∼ −25 km/s. As opposed to the other velocity ranges, no prominent CS detections have been found at v _lsr ∼ −25 km/s in this region, and thus we suggest that the MC is located on the far distance.

Figure 7. Averaged CS(1–0), CH₃OH(I) and CO(1–0) emission from the regions towards HESS J1418–609 and HESS J1420–607 shown in Figure 6. The red lines indicate the fit used to model the emission at the velocity ranges where the regions are shown in Figure 6 (see Table C.4 for fit parameters). The red vertical lines show the dispersion measure distance of the pulsar PSR J1420–6048 (P1 in Figure 6). The pink, cyan and brown regions show the velocity ranges for the integrated intensity maps in Figure 6.

At v _lsr = −65 to −55 km/s (see Figure 6a), we note the bulk of the CO(1–0) emission is located towards the west and north side of HESS J1418–609. We also remark that the molecular emission shows little overlap with any of the TeV sources at this velocity range. We observed extended CS(1–0) emission, labelled ‘1’ to ‘5’, north of HESS J1418–609. Assuming a distance d = 5.6 kpc, the combined mass of these clumps attains $M_{\textrm{H}_2}\left(\textrm{CS}\right)=1.1\times10^{4}{\textrm{M}_{\odot}}$ . Notably, the bulk of the CS(1–0) emission appears to wrap around the TeV emission (as expected for leptonic IC emission).

At v _lsr = −52 to −42 km/s, we note that the CO(1–0) emission overlaps the two TeV sources, and peaks west of HESS J1018–609 and south east of HESS J1420–607. However, due to the small velocity separation between the components and the components at v _lsr = −65 to −55 km/s, it is somewhat difficult to accurately describe the morphology of the diffuse molecular gas at these velocities. We particularly observe that the prominent CO emission in region ‘A’, with mass $M_{\textrm{H}_2}\left(\textrm{CO}\right)=3.5\times10^{4}{\textrm{M}_{\odot}}$ and averaged density $n_{\textrm{H}_2}\left(\textrm{CO}\right)=5.6\times10^{2}\,{\rm {cm}}^{-3}$ , anti-corresponds with the TeV source HESS J1418–609. Our CS(1–0) results highlight a filamentary structure (see regions ‘6’ to ‘10’ in Figure 6d.), with a combined CS mass $M_{\textrm{H}_2}\left(\textrm{CS}\right)=5.2\times10^{3}{\textrm{M}_{\odot}}$ (CO mass $M_{\textrm{H}_2}\left(\textrm{CO}\right)=3.1\times10^{4}{\textrm{M}_{\odot}}$ ) which crosses HESS J1420–607. Interestingly, the dense gas in regions ‘6’ to ‘8’ also appears in a shell-like arrangement centred towards the south-east of HESS J1420–607. CS(1–0) clumps, labelled ‘11’ to ‘13’ are also found towards HESS J1418–609. The regions ‘11’ and ‘12’ are in fact coincident with the TeV peak emission.

At v _lsr = −10 to 0 km/s (bottom panels), we found prominent CO emission south of HESS J1418–609, east of HESS J1420–607, and particularly between the two pulsars. Notably, the latter molecular region (labelled ‘B’) with mass $M_{\textrm{H}_2}\left(\textrm{CO}\right)=6.8{\textrm{M}_{\odot}}$ nests an extended dense clumps (see region ‘14’) with mass $M_{\textrm{H}_2}\left(\textrm{CS}\right)=2.7{\textrm{M}_{\odot}}$ and with averaged density reaching $n_{\textrm{H}_2}\left(\textrm{CS}\right)=3.8\times10^{4}\,{\rm{cm}}^{-3}$ .

In Figure E.2, we compare the column density from CO and Hi measurements with the column density $N_{\textrm{H}}=2.7\pm0.2\times10^{22}\,{\rm{cm}}^{-2}$ towards the pulsar PSR J1418–6058 derived from X-ray measurements (Kishishita et al. Reference Kishishita, Bamba, Uchiyama, Tanaka and Takahashi2012). Assuming that all components, except that at v _lsr ∼ 25 km/s, are in the near distance, we observe in Figure E.2 that the column density value from X-ray measurements is matched at v _lsr ≤ −50 km/s, which infers a distance range d ∼ 3.5 kpc. As a result, this column density X-ray study appears to favour a PWN location at d ∼ 3.5 kpc (v _lsr = −52 to −42 km/s), while our morphological study of the molecular ISM favours a PWN distance at d ∼5.6 kpc (v _lsr = −65 to −55 km/s), as it anti-corresponds with the TeV gamma-ray emission (see Figure 6).

3.5. HESS J1303-631

HESS J1303–631 was first classified as a ‘dark source’ due to its lack of any counterparts at other wavelengths (Aharonian et al. Reference Aharonian2005; Acero et al. Reference Acero2013). However, energy-dependent morphology of the TeV source (see Abramowski et al. Reference Abramowski2012b) unambiguously highlighted its association with the pulsar PSR J1301–6305 (P1 in Figure 8) with spin down energy $\dot{E}_{\textrm{SD}}=2.6\times10^{36}\,{\rm {erg\,s^{-1}}}$ , a rotation period P = 184 ms, and a characteristic age τ _c = 11 kyr. Follow-up observations with XMM-Newton revealed diffuse X-ray emission towards the pulsar with a power-law spectral index $\Gamma_X=2.0^{+0.6}_{-0.7}$ (Abramowski et al. Reference Abramowski2012b). Acero et al. (Reference Acero2013) detected a GeV counterpart with a gamma-ray spectral index Γ _γ = 1.7. Finally, from the 1.384-GHz ATCA observations, Sushch et al. (Reference Sushch, Oya, Schwanke, Johnston and Dalton2017) recently announced the presence of a plausible SNR radio shell with radius ∼12′ next to the pulsar PSR J1301–6305, although this association appears unlikely.

Figure 8. Nanten CO(1–0) integrated intensity map towards HESS J1303–631 (shown in black contours) between v _slr = −35 to −25 km/s, −25 to −15 km/s, −5 to +5 km/s, and 15 to 35 km/s. The 7-mm map is shown as a black-dashed box, while the position and size of the SNR candidate from Sushch et al. (Reference Sushch, Oya, Schwanke, Johnston and Dalton2017) is indicated by a blue dashed circle. The position of the pulsar PSR J1301–6305 (P1) is shown in cyan diamond. The regions labelled ‘1’ to ‘7’ where CS(1–0) was detected are shown in black circles. The positions of prominent CO detections slightly overlapping the TeV emission are shown in cyan ellipses.

Based on the dispersion measure of this pulsar, Cordes et al. (Reference Cordes, Lazio, Chatterjee, Arzoumanian and Chernoff2002) suggested a distance d ∼ 6.6 kpc, much closer than the previous distance d ∼ 12.6 kpc (Taylor and Cordes Reference Taylor and Cordes1993). From the Nanten CO(1–0) components identified (see Figure 9), we focus on several molecular complexes in the line of sight at v _lsr = −35 to −25 km/s (distance d ∼ 6.6 kpc, Scutum Crux arm), v _lsr = −25 to −15 km/s (distance d ∼ 1.5 kpc, near Sagittarius–Carina arm), v _lsr = −5 to 5 km/s (d ∼ 0.1 kpc, local arm) and v _lsr = 25 to 35 km/s (distance d ∼ 12.6 kpc, far Sagittarius–Carina arm) shown in Figure 8. From our 7-mm CS observations, which cover the north-west part of the SNR towards the TeV source, we have identified several molecular clumps, which we have labelled ‘1 to 7’ (see Tables C.5 and D.4 for their physical parameters) but no extended CS(1–0) emission has been detected.

Figure 9. The averaged CS(1–0) and CO(1–0) emission from the different regions shown in Figure 8 towards HESS J1303–631. The red lines indicate the fit used to parametrise the emission. The fit parameters are displayed in Table C.5. Finally, the red and green vertical lines represent the dispersion measure distance of the pulsar P1 as predicted by Cordes et al. (Reference Cordes, Lazio, Chatterjee, Arzoumanian and Chernoff2002) and Taylor and Cordes (Reference Taylor and Cordes1993), respectively. The pink, cyan, brown, and yellow regions indicate the velocity range of the integrated intensity maps shown in Figure 8.

At all the aforementioned velocity ranges, it appears that some CO(1–0) emission always overlaps the TeV emission. It should be noted that the CO emission peaks south of HESS J1303–631 inside region ‘A’ at v _lsr = −35 to −25 km/s. At the local arm (v _lsr = −5 to 5 km/s), most of the CO(1–0) emission is distributed north east of the TeV source. At v _lsr = 15 to 35 km/s, prominent CO emission is found overlapping HESS J1303–631 (see regions ‘C’ and ‘D’). We also note that the CO(1–0) emission does overlap the position of the SNR candidate represented in blue circle in Figure 8 at v _lsr = −35 to −25 km/s and v _lsr = 15 to 35 km/s.

Interestingly, we have found a CO(1–0) emission dip at v _lsr = −25 to −15 km/s localised towards the SNR candidate. From the position–velocity plots shown in Figure 10, we observe prominent CO emission which in fact appears to surround the SNR position (whose boundaries are shown in red dashed lines in Figure 10 bottom panels) between v _lsr = −22 to −15 km/s. From the CO(1–0) integrated intensity region at this velocity range, we also observe little spatial overlap between the molecular gas and the SNR candidate. Consequently, it may highlight the presence of a putative molecular shell at d ∼ 1.5 kpc coincident with the recently observed SNR candidate. The green cross and ellipse in Figure 10 indicate the position and the expansion speed (v _exp ∼ 4 km/s) of the putative molecular shell surrounding the SNR candidate. From the dashed magenta circular region in Figure 10 centred towards the SNR candidate position, the maximum mass swept by the SNR candidate reaches $M_{\textrm{H}_2}=3.3\times10^{4}{\textrm{M}_{\odot}}$ . We obtain as a upper limit a required kinetic energy $E_{\textrm{kin}}=1/2M_{\textrm{H}_2}v_{\textrm{exp}}^{2}=5.1\times10^{47}$ erg which represents ∼0.05% of the total kinetic energy from powerful O stars stellar winds over 1-Myr time-scale (see Weaver et al. Reference Weaver, McCray, Castor, Shapiro and Moore1977 for detailed study on interstellar bubbles). Consequently, the molecular shell may have been produced by the SNR progenitor star. We thus suggest this new SNR candidate is at distance at d ∼ 1.5 kpc, supporting its non-physical association with PSR J1301–6305 claimed by Sushch et al. (Reference Sushch, Oya, Schwanke, Johnston and Dalton2017).

Figure 10. (panel a) Nanten CO(1–0) integrated intensity between v _lsr = −22 to −15 km/s overlaid by the HESS TeV contours in black. The blue dashed circle indicates the size of the candidate SNR (Sushch et al. Reference Sushch, Oya, Schwanke, Johnston and Dalton2017). The cyan diamond shows the position of the pulsars PSR J1301–6305 (P1). The purple dashed circle indicates the region used to compute the mass of the putative molecular shell (see discussion in Section 4). The green grid of boxes indicates the position of the displayed CO(1–0) spectral lines (panel b). (panels c and d) Galactic longitude-velocity (l, v) and latitude-velocity (b, v) images integrated between l = [304.25°:304.55°] and b = [–0.34°:–0.04°], respectively (shown as red dashed lines in top left panel). The green cross-hair and ellipse show the location of a putative expanding molecular shell while the red dashed lines indicate the boundaries of the candidate SNR.

3.6. HESS J1018–589

HESS J1018–589 was first reported by Abramowski et al. (Reference Abramowski2012a) and actually consists of two distinct sources. The gamma-ray binary 1FGL J1018.6–5856 appears to be responsible for the HESS J1018–589A TeV emission (Abramowski et al. Reference Abramowski2012a), while HESS J1018–589B (shown in dashed black circle in Figure 11) is thought to be a PWN powered by the pulsar PSR J1016–587 (shown by a cyan diamond), with a rotation period P = 107 ms, a spin down energy $\dot{E}_{\textrm{SD}}=2.6\times10^{36}$ erg, and a characteristic age τ _c = 21 kyr (Camilo et al. Reference Camilo2001). It has been suggested that the pulsar, with dispersion measure distance d = 8 kpc, is not associated with the nearby SNR G292–1.8 located at d ∼ 2.9 kpc (Ruiz & May Reference Ruiz and May1986).

Figure 11. Nanten CO(1–0) integrated intensity map towards HESS J1018–589 between v _lsr = −23 to 10 km/s overlaid by the TeV gamma-ray emission from HESS 1018–589 in solid black contours. The dashed black circle indicates the size and position of HESS J1018–589B. The SNR G284.3–01.8 is shown as a blue dashed circle, while the pulsar PSR J1019–5857 is shown as a cyan diamond. The extended CO region labelled ‘A’ is shown in cyan, the position of the SiO(1–0, v = 0) ‘S1’ are displayed in yellow, and the CH₃OH maser found in the region ‘CH1’ is shown in green. Their respective spectral lines are displayed on the right-hand side. The pink region illustrates the aforementioned velocity range.

From the Nanten CO(1–0) observations shown in Figure 11, we have identified CO emission at v _lsr = −23 to −10 km/s, inferring a near distance d ∼ 2.8 kpc, matching the SNR distance. The CO(1–0) emission appears filamentary north of HESS J1018–589B and partially overlaps the northern rim of SNR G292–1.8. The molecular gas in region ‘A’, with mass attaining $M_{\textrm{H}_2}\left(\textrm{CO}\right)=2.9\times10^{3}{\textrm{M}_{\odot}}$ , shows quite broad emission (Δv ∼ 12 km/s, see Table C.6 and Figure 11). Although no CO emission was revealed at v _lsr ∼ 30 km/s, we note from Figure E.3 significant Hi emission, which may suggest that the ISM surrounding HESS J1018–589B mostly consists of atomic gas.

4.1. TeV emission from CRs

Based on the mass estimates towards molecular regions overlapping the TeV sources, we use eq. 10 from Aharonian (Reference Aharonian1991) to derive the CR enhancement factor k _CR = w _CR/w _⊙ which represents the ratio between w _CR the energy density of CRs towards an MC next to a TeV source, and w _⊙ ∼ 1.0 eV cm⁻³ the energy density found in the solar neighbourhood. We use the combined atomic and molecular mass to obtain the total amount of target material available. Table 2 indicates the k _CR values required for each molecular regions (partially) overlapping the TeV sources to account for the observed TeV fluxes. The k _CR value can then be compared to typical CR enhancement factors predicted in the vicinity of SNRs (and potentially towards PWNe) to check the plausibility of hadronic contribution.

Table 2. Cosmic-ray enhancement factors k _cr derived using eq. 10 from Aharonian (Reference Aharonian1991), required to reproduce the observed TeV emission above 1 TeV $F(\gt 1\,{\rm TeV}) = \int N_0 E_{\gamma}^{-\Gamma}{\rm d}{E}_{\gamma}$ via CR–ISM interaction.

*Aharonian et al. (Reference Aharonian2007), ^†Abramowski et al. (Reference Abramowski2011), ^‡Abdalla et al. (Reference Abdalla2018a), ^±, ^∐ Aharonian et al. (Reference Aharonian2006), Abramowski et al. (Reference Abramowski2012b), *Abramowski et al. (Reference Abramowski2015)

^a We scaled down the HESS J1809–193 photon flux by 16%, 16% and 12% for regions ‘1’, ‘2’, ‘3’, respectively, corresponding to the ratio between the molecular regions and the TeV emission sizes.

^b We scaled down the HESS J1420–607 photon flux by 4%, ∼1%, <1%,, 1%, and 3% for the regions 6 to 10, respectively.

^c We scaled down the HESS J1303–631 photon flux by 44%, 7% and 7% for region ‘A’, ‘B’, and ‘C’, respectively.

4.1.1. Contribution from nearby SNRs?

Nearby SNRs are the most likely viable candidates to produce CRs energy densities up to ∼10³ eV cm⁻³ (see Reynolds Reference Reynolds2008 and references therein). CRs propagate along the magnetic field lines and they scatter from their interaction with magnetic field perturbations (provided the scale of the perturbation roughly equals the CR gyroradius). As the magnetic field turbulence is thought to be enhanced in MCs, we here assume an isotropic diffusion of CRs and electrons as a first-order approximation. If we assume an isotropic diffusion of CRs and neglect energy losses, we can estimate the energy density distribution of CRs at a distance R from the SNR (here assumed as an impulsive source of CRs; see Aharonian & Atoyan Reference Aharonian and Atoyan1996):

(1) $$\begin{align} n\left(E,R,t\right)&=\frac{\eta_{\textrm{pp}}E_\textrm{SNR}}{\left(m_\textrm{p}c^2\right)^{2-\alpha}}\frac{E^{-\alpha}}{\pi^{3/2} R_\textrm{d}^{3}}\exp\left(-\left(\frac{R}{R_\textrm{d}}\right)^2\right) \end{align}

(2) $$\begin{align} R_\textrm{d}&=2\left(\chi D_{10}t_\textrm{age}\sqrt{E/10\textrm{ GeV}}\right)^{1/2}\end{align}

with m _p being the proton mass, D ₁₀ = 10²⁸ cm²s⁻¹ being the diffusion coefficient of 10 GeV CRs, α the spectral index of the proton distribution, and η _pp the ratio of the total SNR energy E _SNR transferred to CRs. R _d(t) represents the diffusion radius travelled by CRs after a time t. The diffusion suppression factor χ accounts for slower diffusion of particles which can be caused, for instance, by streaming instabilities or perturbations caused by shocks (see Nava et al. Reference Nava, Gabici, Marcowith, Morlino and Ptuskin2016; Malkov et al. Reference Malkov, Diamond, Sagdeev, Aharonian and Moskalenko2013, and references therein). Here we use, χ = 0.01 to 1 which matches the slow and fast diffusion coefficient regimes defined by Aharonian & Atoyan (Reference Aharonian and Atoyan1996). Various studies of the CRs interacting with nearby molecular clouds W 28 (Giuliani et al. Reference Giuliani2010; Li & Chen Reference Li and Chen2010; Gabici et al. Reference Gabici, Casanova, Aharonian, Rowell, Boissier, Heydari-Malayeri, Samadi and Valls-Gabaud2010) have, in fact, suggested a suppression factor between χ = 0.01 and 0.1. However, the value of χ remains poorly constrained. From Eq. 1, we can then obtain the total energy density of CRs w _CR(R, t) using the following equation:

(3) $$\begin{equation} w_\textrm{CR}\left(R,t\right)=\int_{\epsilon_0}{n\left(E,R,t\right)E\textrm{d}E} \label{getenergy}\end{equation}

where ε ₀ = m _p c ² is the proton energy at rest. Figure 12 illustrates the range of k _CR produced by SNRs with initial energy E _SNR = 10⁵¹ erg as a function of the SNR age, using a standard proton spectral index α = 2.2 and η _pp = 0.1. From the distance between the SNRs and the surrounding ISM regions (see Table 2), and the age of the SNRs, we can then check whether the required enhancement factors from Table 2 fall within the predicted range from nearby SNRs.

Figure 12. Evolution of the cosmic-ray enhancement factor k _CR range as a function of time at a distance d = 10 pc (red), d = 30 pc (green), and d = 50 pc (blue) away from an impulsive source with initial energy E _SNR = 10⁵¹ erg and initial CR spectral index α = 2.2. A energy-dependent diffusion of CRs (see Section 4) has been applied with a diffusion coefficient at 10 GeV bounded between D ₁₀ = 10²⁵ to 10²⁸ cm² s⁻¹.

4.1.2. CR contribution from PWNe?

Additionally, a few authors (see Amato Reference Amato2014 and references therein) have also argued that high energy hadrons could also be produced inside the pulsar environment, and be responsible for several features inside the PWN (e.g. wisps in the Crab PWN; Gallant Arons Reference Gallant and Arons1994). Providing high energy CRs have indeed been produced inside PWNe and not suffered heavy adiabatic losses, we will thus discuss whether the pulsars considered in this work could also generate the required CR enhancement factors shown in Table 2. In order to model the high CR energy density potentially produced by those pulsars, we account for the evolution of the spin down power $\dot{E}_\textrm{SD}$ as a function of time t, which can be described as follows:

(4) $$\begin{eqnarray} &\dot{E}_\textrm{SD}\left(t\right)=\dot{E}_\textrm{SD}\left(t_\textrm{age}\right)\left(1+\left(n_\textrm{b}-1\right)\frac{\dot{P}\left(t-t_\textrm{age}\right)}{P}\right)^{-\Gamma} \end{eqnarray}

(5) $$\begin{eqnarray} &\Gamma=\frac{n_\textrm{b}+1}{n_\textrm{b}-1} \end{eqnarray}

with t _age being the current age of the PWN, P(t) and $\dot{P}\left(t\right)$ being the current period and period derivative of a pulsar, respectively, at time t, $\dot{E}_\textrm{SD}\left(t\right)$ being the pulsar spin down power at time t, and n _b being the pulsar braking index. In order to obtain the density of CRs at a given radius R from the pulsar, we rewrite Eq. 1 using the source term $S\left(t\right)=\eta_\textrm{pp}\dot{E}_\textrm{SD}\left(t\right)/\left(m_\textrm{p}c^2\right)^{2-\alpha}$ , with η _pp being the fraction of the spin down power transferred to CRs. We thus numerically solve the following equation:

(6) $$ \begin{equation} n\left(E,R,t\right)=\frac{E^{-\alpha}}{\pi^{3/2}}\int_{t_\textrm{age}}^{0}{-\frac{S\left(\xi-t_\textrm{age}\right)}{R_d^3\left(\xi\right)}\exp\left(-\frac{R^2}{R_d^2\left(\xi\right)}\right)\textrm{d}\xi} \label{diffusionpulsar} \end{equation}

with ξ = t _age − t. As per the SNR scenario, we obtain the CR energy density at a given distance R from pulsar using Eq. 3.

Figure 13 illustrates the CR energy density produced by the various pulsars as a function of the diffusion coefficient suppression factor χ. Here, we assumed n _b = 3 for all pulsars, except for PSR J1119–6127 in which we used the measured braking index n _b = 2.68. From the spectral modelling towards several PWNe, Bucciantini, Arons & Amato (Reference Bucciantini, Arons and Amato2011) suggested that the maximum energy rate available to be transferred to CRs must be at most 0.20 of the total spin down power. We thus have used η _pp= 0.20 and the k _CR values shown in Figure 13 are consequently upper limits. In the following subsections, we compare these predictions with the various required k _CR (see Table 2) derived inside the molecular regions.

Figure 13. Predicted energy density k _CR from pulsars at 10 pc (top panel) and at 20 pc (bottom panel) distance as a function of the diffusion coefficient suppression factor χ (see colour version online).

4.2. TeV emission from high energy electrons

As opposed to CRs, high energy electrons suffer heavy radiation losses as they diffuse inside the dense ISM because of the potentially enhanced magnetic field strength (see Crutcher et al. Reference Crutcher, Wandelt, Heiles, Falgarone and Troland2010). In the lack of intense radiation fields, synchrotron losses, with time-scale $\tau_{\textrm{sync}}=2.45\times10^{7}\big(B^2_{\textrm{mg}}E/m_\textrm{e}c^2\big)^{-1}$ yr (B _mg = B/1_mG being the magnetic field strength and E being the total energy of the electron), are likely to dominate over IC losses inside MCs. Although Bremsstrahlung emission, with time-scale τ _brem ∼ 4.6 × 10⁷(n/1 cm⁻³)⁻¹ yr, may contribute inside MCs with densities n _H > a few × 10³ cm⁻³ in the 0.1 to 1 TeV band, we expect the TeV emission above 1 TeV to anti-correspond with the molecular ISM as the IC radiation is likely to dominate above 1 TeV. Table 3 indicates the synchrotron energy loss time-scale τ _sync for electrons with energy E = 10 TeV, required to produce photons with energy E_γ ∼ 1 TeV. To obtain B values inside MCs, we use the Crutcher et al. (Reference Crutcher, Wandelt, Heiles, Falgarone and Troland2010) relation:

(7) $$\begin{eqnarray} & B=10\,\mu\textrm{G} & \textrm{for } n_{\textrm{H}} \lt 300\,\textrm{cm}^{-3} \end{eqnarray}

(8) $$\begin{eqnarray} & B=10\left(\dfrac{n_{\textrm{H}}}{300\textrm{ cm}^{-3}}\right)^{0.6}\,\mu\textrm{G} &\textrm{for } n_{\textrm{H}}\ge300\,\textrm{cm}^{-3}. \end{eqnarray}

We compare these time-scales with the age of nearby high energy sources and the diffusion time-scale τ _diff required for the particles with energy above 10 TeV to fully cross a molecular region. We here define τ _diff = (〈r〉 + d _source)²/6D (E, B, χ), with 〈r〉 being the mean radius of the molecular region and d _source being the distance separating the centre of the high energy source to the centre of the molecular region. To check the leptonic scenario, we first select dense molecular regions which anti-correspond with each TeV source and/or its X-ray counterpart. Here, we preferably use the extended CS regions as most of the observed gas is confined. However, in the cases where no CS(1–0) emission was detected (i.e. HESS J1018–189B, HESS J1303–631, and HESS J1119–164), we use the CO regions. We again note that the densities $n_{\textrm{H}_2}$ from our CO analyses are upper limits. Thus, these synchrotron time-scale values should be, in these cases, used as upper limits. Then, we derive the synchrotron time-scale and compare it to τ _diff required for the high energy electrons to cross the cloud. In the case where the synchrotron time-scale is much smaller than the diffusion time-scale, we expect the high energy electrons to radiate most of their energy and thus become unable to produce X-rays and TeV emission inside the MC. An anti-correlation between the ISM and X-rays would support the physical connection between the TeV source and the ISM as well as the leptonic nature of the TeV source (see Sano et al. Reference Sano2013).

Table 3. The diffusion time-scale t _diff = (〈r〉 + d _source)²/6D(E, B) for a particle with energy E = 10 TeV to cross the listed molecular regions with mean radius 〈r〉 and positioned at d _source from the centre of the listed sources. The age of the source t _age, the CR–ISM interaction time-scale τ _pp for protons, and the synchrotron time scale τ _sync for electrons with energy E = 10 TeV are also displayed as means of comparison.

^a Obtained from eq. 21 from Crutcher et al. (Reference Crutcher, Wandelt, Heiles, Falgarone and Troland2010) (see Section 4.2).

^b $\tau_{\rm pp} = 6 \times 10^{7}/2n_{{\rm H}_{2}}$ where $n_{{\rm H}_{2}}$ can be found in Tables D.1–D.6.

^c $\tau_{sync} = 2.45 \times 10^{7} (B^2_{mG}E/m_{e})^{-1}$ yr (Ginzburg & Syrovatskii Reference Ginzburg and Syrovatskii1964).

^d See eqs. 2 and 3 from Gabici, Aharonian, Blasi (Reference Aharonian2007) to obtain the diffusion time-scale. The results shown assume a diffusion coefficient suppression factor χ = 0.1

^e SNR candidate, see text and Sushch et al. (Reference Sushch, Oya, Schwanke, Johnston and Dalton2017) for further details and Figure 8 for location. We used standard SNR expansion parameters to provide an age estimate of the SNRc if located at this distance (see text).

4.3. Discussion of Individual Sources