1 Introduction

Cosmic-ray acceleration in supernova remnants (SNRs) is the most promising mechanism for accelerating Galactic cosmic rays, which mainly comprise relativistic protons with energies less than 1015.5 eV. Considerable efforts have been devoted to theoretical works to elucidate the details of particle acceleration (e.g., Bell 1978; Blandford and Ostriker 1978). Recent progress in X-ray and \(\gamma \)-ray observations has allowed us to explore the origin of cosmic rays. Young SNRs, which ages around 2000 yr, are of particular interest because they emit higher-energy X-rays and \(\gamma \)-rays than SNRs younger than 1000 yr (e.g., G1.9+0.3 at \(156\pm 11\) yr, Carlton et al. 2011) or those older than several tens of thousands of years (e.g., W44 at ∼20,000 yr, Wolszczan et al. 1991). In the present article, we focus on the young SNRs RX J1713.7−3946 (hereafter RXJ1713) and RCW 86, which are characterized by the emission of bright X-rays and TeV \(\gamma \)-rays.

The origin of \(\gamma \)-rays from young SNRs is currently under debate, and two mechanisms have been suggested to explain them. One is a leptonic process in which cosmic-ray electrons collide with low-energy photons, boosting them into the \(\gamma \)-ray regime via the inverse Compton process. The other is a hadronic process in which cosmic-ray protons collide with interstellar protons to produce neutral pions that decay into two \(\gamma \)-rays. It is important to verify whether the \(\gamma \)-rays are of hadronic origin for establishing whether SNRs are a major source of cosmic rays in the Galaxy.

In previous works that considered the hadronic model for the \(\gamma \)-ray emissivity of an SNR, the number density of targeted interstellar protons has been used (e.g., to estimate the cosmic-ray energy from the \(\gamma \)-ray luminosity). However, strongly shocked matter (e.g., shock-ionized plasma traced by optical lines and free–free radio/X-ray continuum and shocked molecular hydrogen observed as near-infrared lines), which occupies only a small portion of the interstellar medium (ISM), has mainly been considered as the target gas. Neutral molecular and atomic hydrogen gas has not been considered to play a major role in controlling the radiative properties of SNRs, except for some pioneering works (e.g., Aharonian et al. 1994), while the interaction of SNRs with the ISM has been a subject of interest in some other works (e.g., Dubner and Arnal 1988; Dubner et al. 1998, 1999, 2002a,b, 2004; Tatematsu et al. 1987, 1990a,b; Seta et al. 1998, 2004; Arikawa et al. 1999). In this article, we focus on the bulk neutral gas that dominates the mass the ISM. Recent theoretical studies have revealed that the associated ISM may play a major role in determining the radiative properties of SNRs, and clarifying clarify the role of inhomogeneities in the ISM in regulating the high-energy radiations is important (e.g., Inoue et al. 2012a; Celli et al. 2019). Here we review methods used to identify the ISM associated with a SNR, a theoretical model that incorporates shock–cloud interactions, and their theoretical implications for the origin of very-high-energy \(\gamma \)-rays, with a focus on the SNRs RXJ1713 and RCW 86.

2 The Interstellar Medium (ISM)

The neutral ISM consists of the warm neutral medium (WNM) and the cool neutral medium (CNM) (e.g., Draine 2011). The WNM has a density of 0.1–10 cm−3 and a spin temperature of 300–10000 K, while the CNM has a density of 10–103 cm−3 and a spin temperature of 10–300 K (e.g., McKee and Hollenbach 1980; Heiles and Troland 2003). The masses of the two phases are comparable. The CNM comprises clumps with a volume filling factor of a few percent, while the WNM occupies a large volume of inter clump space (Draine 2011; Fukui et al. 2018). The CO clouds are formed in the CNM clumps, and their volume filling factors are also very small (∼3–4%, e.g., Inoue et al. 2012b; Tachihara et al. 2018). This picture is significantly different from the conventional assumptions made for (one-dimensional) models, which do not appropriately represent these two-phase inhomogeneities.

The 21 cm Hi emission and 2.6 mm CO emission are good tracers of the ISM associated with a SNR. Comparing a high-resolution X-ray image with the CO and Hi distributions provides a powerful method for identifying the interaction between an SNR and the ISM. Because electrons cannot penetrate the dense ISM clumps, we find a small-scale anti correlation between the X-rays and dense ISM clumps at sub-pc scales. Calculations of the ISM column density can be performed using a CO-to-H2 conversion factor (for a review, see Bolatto et al. 2013) and the Hi 21 cm emission, where corrections for the Hi optical depth are required for the dense CNM.

Although the scheme described above works successfully in general, one caveat is that the intermediate density regime between 100 and 1000 cm−3 may not be detectable in the CO emission because of the low CO abundance caused by photo dissociation. In addition, the Hi emission may become optically thick in this density regime, requiring a correction for saturation. Hi is observed to have self-absorption dips; thus, some Hi gas is optically thick (e.g., Sato and Fukui 1978).

The existence of dark gas, which cannot be probed by either CO or Hi emission, was suggested by EGRET \(\gamma \)-ray observations (Grenier et al. 2005) as well as from submillimeter observations by Planck (e.g., Planck Collaboration 2011). Two possible explanations exist for the physical entity of the dark gas: CO-dark H\(_{\mathrm{2}}\) gas and optically thick Hi gas. The concept of CO-dark H\(_{\mathrm{2}}\) gas is based on numerical simulations of the evolution of the H\(_{\mathrm{2}}\)/Hi gas conducted by Wolfire et al. (2010), who showed that the self-shielding of H\(_{\mathrm{2}}\) against UV is strong and a phase of H\(_{\mathrm{2}}\) without an observable abundance of CO can exist outside a CO cloud. This gas, which comprises H\(_{\mathrm{2}}\) without CO, is a reasonable model provided that Hi is converted mainly into H\(_{\mathrm{2}}\). However, Fukui et al. (2014, 2015) used the Planck/IRAS dust emission as a proxy for the hydrogen column density and found that the Hi emission becomes optically thick under the usual Hi gas condition, in which the Hi column density is assumed to be proportional to the dust optical depth. The saturated Hi intensity becomes weaker than in the optically thin case, resulting in a smaller Hi column density than under the conventional optically thin approximation. The column density corrected for the optical depth becomes larger than the optically thin Hi column density by a factor of 1.5–2; this increased column density is large enough to explain the dark gas (Okamoto et al. 2017; Hayashi et al. 2019a,b; Mizuno et al. 2020). The optically thick Hi raises the possibility of Hi being dominant compared with H\(_{\mathrm{2}}\), which requires theoretical justification. Simulations by Inoue et al. (2012b) using a gas-evolution code similar to that employed by Wolfire et al. (2010) presented another evolutionary model of Hi flows showing Hi dominant outside the CO clouds. An important difference exists between the initial conditions in these two simulations. Inoue et al. (2012b) adopted a lower Hi density, which includes both the CNM and WNM. This is a realistic Hi distribution compared with dense Hi comprising only the CNM, which was the initial condition assumed by Wolfire et al. (2010). The lower density leads to slower Hi–H\(_{\mathrm{2}}\) conversion than the high-density initial condition and produces a rich Hi envelope surrounding the CO gas instead of CO-dark H\(_{\mathrm{2}}\) gas. This lends theoretical support to optically thick Hi as an alternative. A possible new probe for this density regime is Ci emission at submillimeter wavelengths, which can serve as another useful tool (Tachihara et al. 2018).

3 RX J1713.7−3946

3.1 Distance determination: nonthermal X-rays and the interacting ISM

Except for several nearby SNRs with historical records, determining distances to SNRs in the Milky Way is not straightforward, because most of SNRs lie in the Galactic plane and are heavily obscured. Although the relation between radio surface brightness and angular diameters (\(\Sigma \)–D relation) has been used to estimate the distances to SNRs, the \(\Sigma \)–D relation does not provide an accurate measure of distance because of its large scatter (e.g., Pavlovic et al. 2014). Shock-excited masers—e.g., OH (1720 MHz) masers—are very useful for identifying the shocked interface, and the distance to SNRs can be obtained from the radial velocities of the masers, although this method can only be used for middle-aged SNRs that have interacted with the ISM for a long time (e.g., Chen and Jiang 2013). The ISM, including CO clouds and Hi gas, associated with an SNR also provides a kinematic velocity, using which the distance can be calculated by comparison with a kinematic model of the Galaxy. The uncertainty is mainly due to the random cloud motions, which are on the order of 10 km s−1.

The nonthermal X-ray SNR RXJ1713 is a case for which two distances that differ by a factor of 6 had been debated. It was first suggested that this SNR is located at a distance of 1 kpc, assuming a typical average Hi density for the X-ray absorption (Koyama et al. 1997), and it was identified as the SNR of 393AD from Chinese historical records (Wang et al. 1997). Slane et al. (1999) claimed that CO clouds at −69 and −94 km s−1 are associated with the SNR, from 8.8 arcmin resolution maps obtained with the CfA 1.2-m telescope, and they derived a large kinematical distance (∼6 kpc). These authors argued that the Hi density toward RXJ1713 is lower than the typical Galactic average owing to an accidental coincidence with a hole in the ISM, while a three-times-larger distance implies an unusually large radius of ∼50 pc for an age of \({\sim}10^{4}\) yrs. The 6 kpc distance has been questioned by a subsequent high-resolution CO survey, which found another CO component at −10 km s−1 to be associated with this SNR at the 2.6 arcmin resolution CO maps obtained with the NANTEN 4-m telescope, as shown in Fig. 1 (Fukui et al. 2003; Moriguchi et al. 2005). The observational signature of this association is the large-scale correlation (on scales of a few pc) between the CO map and the ROSAT X-ray survey, which was confirmed in more detail by the Suzaku data from Sano et al. (2010, 2013) (see Sect. 3.2). As the low velocity CO gas is fainter by an order of magnitude than the −69 and −94 km s−1 clouds, it was ignored by Slane et al. (1999). The smaller distance is now widely accepted, and it has been established that RXJ1713 is a young SNR aged 1600 yr, which was historically recorded (Aharonian et al. 2004; Cassam-Chenaï et al. 2004; Tanaka et al. 2008, 2020; Sano et al. 2015a; Tsuji and Uchiyama 2016). This shows that high resolution and sensitivity are crucial in a comparison with the ISM, raising questions about the distance determinations in previous low-resolution studies.

Fig. 1
figure 1

Map of the ROSAT X-ray flux superposed on the NANTEN CO(\(J\) = 1–0) intensity contours toward RX J1713.7−3946. The CO velocity range is integrated from −11 to −3 km s−1, which corresponds to the velocity component interacting with the supernova remnant (SNR). The lowest contour level and contour intervals are 4 K km s−1. From Fukui et al. (2003) with permission

Full size image

3.2 Detailed correspondence between X-rays and ISM clumps

Fukui et al. (2012) showed that RXJ1713 is associated with molecular and atomic gas, each having masses of \({\sim}10^{4}\) \(M_{\odot }\). The ISM distribution is shell-like, which an 8 pc radius, and the molecular cloud cores are inside the shell. The shell is consistent with a core-collapse SNR, for which the stellar winds from the high-mass progenitor evacuated the surrounding ISM prior to the SN explosion.

A comparison of the Suzaku X-ray data with the CO image revealed details of the interaction (Sano et al. 2010). Figure 2 shows that the CO peaks A–E, G, I, and L are anticorrelated with the X-ray peaks,Footnote 1 indicating that the dense cores hinder cosmic-ray electrons from penetrating them, and the X-rays are limb-brightened on the surfaces of the CO cores (see also Maxted et al. 2012). The cores in the shell were formed via gravitational instability over a timescale of Myr prior to the SN explosion, and three of them (CO peaks A, C, and D) are forming young stars, as shown by the infrared point sources (Sano et al. 2010). In fact, peak C shows obvious signs of active star formation, including bipolar outflows and an IRAS point source (see Fig. 2a).

Fig. 2
figure 2

(a) Suzaku 5–10 keV X-ray image toward the northwestern shell of RX J1713.7−3946. The superposed contours represents the 12CO(\(J\) = 2–1) integrated intensity obtained with NANTEN2. The black crosses indicate the positions of IRAS (Infrared Astronomical Satellite) point sources. The area enclosed by the black box is shown enlarged in Fig. 2(b). (b) Enlarged view of the \(\mathit{Suzaku}\) X-ray image toward CO peak C superposed on the 12CO(\(J\) = 4–3) core. The black contours indicate the total integrated intensity. The blue and red contours represent the blue- and red-shifted outflow components, respectively (from Sano et al. 2010, reproduced by permission of the AAS)

Full size image

Sano et al. (2013) conducted a comprehensive study of the dense cores probed by CO as well as by cold Hi, and found that the X-ray intensities are well correlated with the masses of the ISM cores, as shown in Fig. 3. Sano et al. (2013) showed that more than 80% of the X-ray peaks, brighter than \(5 \times 10^{-4}\) counts s−1 pixel−1, are associated with the ISM cores, suggesting a tight close connection between X-ray emission and the masses of the interacting ISM cores (Fig. 4). This suggests that the distribution of the dense gas had a significant effect on regulating the X-ray distribution in the recent 1000 yr, assuming a shock velocity of ∼4000 km s−1 (c.f., Tsuji and Uchiyama 2016; Tanaka et al. 2020), which has stimulated a detailed theoretical study of the interaction between a shock front and clumpy ISM.Footnote 2

Fig. 3
figure 3

Correlation between the peak X-ray flux and the clump mass interacting with the SNR shock. The linear regression obtained by least-squares fitting is shown by the solid line, for which the correlation coefficient is ∼0.85 on the double-logarithmic scale. The labels indicate the clump name: A, B, C, D+DW, E+I, G+GE, L, \(\mathrm{O}{+}\mathrm{O}_{\mathrm{b}}{+}\mathrm{O}_{\mathrm{SW}}\), R, U, W, and Z for CO clumps and SE-rim for an Hi clump. From Sano et al. (2013), reproduced by permission of the AAS

Full size image
Fig. 4
figure 4

Schematic illustration of the distributions of molecular (CO) clumps (open crosses), an atomic (Hi) clump (open circle), and X-rays (shaded partial or full circles) superposed on the SNR shell boundary from the \(\mathit{Suzaku}\) X-rays (gray contours). The black open crosses (CO clumps C, E, I, L, and OSW) indicate those fully surrounded by X-rays. From Sano et al. (2013), reproduced by permission of the AAS

Full size image

Most recently, Sano et al. (2020b) obtained high-resolution CO observations (with ∼0.02 pc resolution) toward the northwestern shell of RXJ1713 using the Atacama Large Millimeter/submillimeter Array (ALMA). The authors found dozens of tiny clumps, with typical radii of ∼0.03–0.05 pc (referred to as molecular cloudlets). Figure 5 shows the velocity channel distributions of ALMA CO intensities superposed on Chandra synchrotron X-ray contours. The molecular cloudlets are located not only along the synchrotron X-ray filaments, but also in the vicinity of X-ray hotspots that exhibit flux variations on timescales of a few months or years (see Uchiyama et al. 2007; Higurashi et al. 2020). Because the CO profiles of the cloudlets have narrow widths (see the bottom panels in Fig. 5), the authors concluded that these clumpy structures were formed before the supernova explosion.

Fig. 5
figure 5

Top panels: Velocity channel distributions of Atacama Large Millimeter/submillimeter Array (ALMA) 12CO(\(J\) = 1–0) toward the northwestern shell of RXJ1713. The superposed contours represent Chandra synchrotron X-ray intensities. The yellow crosses represent the positions of X-ray hotspots identified by Higurashi et al. (2020). Typical CO cloudlets—named CL1–6—are also indicated. Bottom panels: CO profiles of CL1–6 (black lines) and Gaussian curves (red lines) fitted using the least-squares method. From Sano et al. (2020b), reproduced by permission of the AAS

Full size image

3.3 Simulations of shock-cloud interactions

Magneto-hydrodynamical (MHD) numerical simulations of a supernova shock propagating in a clumpy ISM show that both the turbulence and the magnetic field are amplified around the dense cores, offering a theoretical basis for the ISM–X-ray correspondence (Inoue et al. 2009, 2012a).

Figures 6(a) and 6(b) show the density distributions at 0 and 1508 yr after since the onset of the shock interaction. The low-density gas (with number density of ∼1 cm−3)—i.e., the WNM—is disturbed significantly by the shock acceleration, while the high-density gas (with number densities exceeding 100 cm−3)—i.e., the CNM—is not much accelerated. In spite of the interaction, the shock front propagates almost freely at its initial velocity of 3000 km s−1 owing to the low filling factor of the CNM cores. Figure 6(c) shows the distribution of the magnetic field which indicates a highly entangled field configuration. The shock front, which was initially planar, is deformed when it hits a dense core and becomes entangled around the core. This produces a turbulent velocity field in the WNM, and the magnetic field becomes turbulent and is amplified to 100 μG or higher from its initial value of 6 μG (Inoue et al. 2009). The amplified magnetic field surrounds the dense cores, as shown in Fig. 6(b), which explains why the synchrotron X-rays are rim-brightened around the dense cores. The simulations show that the magnetic field is amplified around the dense core in a layer of sub pc-scale thickness, where the nonthermal X-rays are enhanced.

Fig. 6
figure 6

(a) Structure of the interstellar medium generated by the thermal instability. The colored image shows the number density of atomic hydrogen, and the black lines represent magnetic field lines. (b, c) Results at \(t = 1508\) yr after the injection of the parallel shock: (b) number density of atomic hydrogen and (c) the magnetic field strength. From Inoue et al. (2009), reproduced by permission of the AAS

Full size image

Figure 7 shows a schematic of the shock–cloud interaction model of Inoue et al. (2012a). Before the supernova explosion of a high-mass progenitor, ambient interstellar gas such as diffuse Hi is completely evacuated by a strong stellar wind. The wind creates a wind bubble with a density of ∼0.01 cm−3 (e.g., Weaver et al. 1977). By contrast, dense molecular clouds, with densities exceeding ∼1000 cm−3 can survive in the wind. Consequently, an inhomogeneous gas distribution with a density difference of five orders of magnitude is formed by the high-mass progenitor, before the supernova explosion occurs. The primary forward shocks from the SN explosion propagate through the wind bubble with clumpy clouds, where particle acceleration operates. Shocks transmitted through the gas cloud are stalled, whereas shock waves in the intercloud medium are not decelerated. These velocity differences are also observed in measurements of the proper motions of the forward shocks (e.g., Tsuji and Uchiyama 2016; Tanaka et al. 2020). The velocity difference generates multiple eddies around the shocked clouds, which enhance the magnetic field strength up to ∼1 mG. Finally, we observe that the shocked clouds are limb-brightened by the synchrotron X-rays. The high magnetic field strength causes short-timescale variability of the synchrotron X-rays, as discovered by Uchiyama et al. (2007). The shocked gas clouds also act as targets for cosmic-ray protons and produce hadronic \(\gamma \)-rays (see Sects. 3.4 and 3.5).

Fig. 7
figure 7

Schematic of the shock–cloud interaction model (see the text). From Inoue et al. (2012a), reproduced by permission of the AAS

Full size image

3.4 The clumpy gas distribution and the \(\gamma \)-ray spectrum

Many previous theoretical works on gamma-ray production via the hadronic process have assumed uniform density or radial density gradients, with a typical ambient density of ∼1 cm−3 (e.g., Ellison et al. 2010; Lee et al. 2013). Even though a low-density environment such as a wind-swept cavity and non-uniformity have been discussed to explain the observed \(\gamma \)-ray spectra, potentially energy dependent propagation effects in the multi phase ISM have not been considered (cf., Berezhko and Völk 2010). One reason why the multi phase ISM has not generally been considered is because MHD effects in SNRs are mostly negligible for low-density contrasts \(\precsim 10^{2}\) (Berezhko et al. 2013). The real ISM environment surrounding a SNR, however, contains preexisting inhomogeneities with density contrasts of ∼105 in the ISM, the importance of which was first emphasized by Cox and Smith (1974) and McKee and Ostriker (1977).

Conventional \(\gamma \)-ray spectra in the hadronic and leptonic scenarios for a uniform (or small density-contrast) ISM are shown in Fig. 8 (Abdo et al. 2011). The hadronic spectrum is relatively flat, with a spectral index \(\Gamma \) of ∼2.0, while the leptonic spectrum is hard in the GeV band (\(\Gamma {\simeq }1.5\)). Therefore, on average, over several SNRs—the shape of the \(\gamma \)-ray spectra—should allow discrimination between the two scenarios if we consider a simple one-zone assumption. However, such assumptions do not hold for individual SNRs, as has been discussed in the community, even for the typical young SNR Cassiopeia A (e.g., Atoyan et al. 2000; Helder and Vink 2008). RXJ1713 is also a representative examples: the inclusion of the clumpy ISM has a significant effect and significantly changes the hadronic spectrum as discussed below (Zirakashvili and Aharonian 2010; Inoue et al. 2012a; Gabici and Aharonian 2014; Celli et al. 2019).

Fig. 8
figure 8

Spectral energy distributions (SEDs) of RX J1713.7−3946 in \(\gamma \)-rays. The green-shaded areas show the uncertainty bands obtained from maximum-likelihood fits, assuming a power-law spectrum extending from 0.5 to 400 GeV. The gray-shaded areas represent systematic uncertainties in the model fits. The solid and dashed curves represent various models (Berezhko and Völk 2006; Porter et al. 2006; Ellison et al. 2010; Zirakashvili and Aharonian 2010). The upper panel shows the hadronic models, whereas the bottom panel displays the leptonic models. From Abdo et al. (2011), reproduced by permission of the AAS

Full size image

The hadronic \(\gamma \)-ray scenario requires a high density of target protons, whereas the large volume of low-density gas is required for efficient particle acceleration via diffusive shock acceleration (DSA). Because the maximum energy of the accelerated particles is proportional to the shock speed, it is assumed to be that a low-density medium (typically less than ∼1 cm−3) and the large volume is needed to efficient acceleration. Under the assumption of a uniform ISM, an ISM density greater than 1 cm−3 is required to produce the \(\gamma \)-rays observed from RXJ1713 (Ellison et al. 2010). For such high densities, strong thermal X-rays due to shock heating are expected. This is contradicted the purely nonthermal X-rays observed from RXJ1713. Although some theoretical studies have avoided this contradiction by considering a wind-bubble model and/or a thermal nonequilibrium state between the electrons and protons downstream (e.g., Morlino et al. 2009; Berezhko and Völk 2010), the absence of bright thermal X-ray emission has been used to exclude hadronic \(\gamma \)-rays and support leptonic \(\gamma \)-rays (Ellison et al. 2012).

However, the lack of thermal X-ray emission in the presence of high gas densities can be explained naturally by considering a realistic clumpy ISM distribution. Most of the volume of the interclump medium is of low density, which allows the DSA to work, whereas a large fraction of the mass of the ISM is contained in dense cores, which can work as a dense, massive, target material, enabling the cosmic-ray protons to produce hadronic \(\gamma \)-rays via the p–p reaction. In the dense cores, the shock fronts are decelerated significantly, with no heating, and thermal X-rays are suppressed, as indicated by the equation 1 (Inoue et al. 2012a):

$$\begin{aligned} k_{\mathrm{B}} T_{\mathrm{c}} = 2 \times 10^{-4}\; \biggl( \frac{v_{\mathrm{sh,d}}}{3000\;\mathrm{km\;s^{-1}}}\biggr)^{2}\; \biggl(\frac{n_{\mathrm{d}}}{0.01\;\mathrm{cm^{-3}}}\biggr) \\ \biggl(\frac{n_{\mathrm{c}}}{10^{3}\;\mathrm{cm^{-3}}}\biggr)^{-1} \;( \mathrm{keV}), \end{aligned}$$
(1)

where \(k_{\mathrm{B}}T_{\mathrm{c}}\) is the maximum temperature of the protons, \(v_{\mathrm{sh,d}}\) is the shock velocity in the intercloud or diffuse region, \(n_{\mathrm{d}}\) is the number density of the intercloud or diffuse gas, and \(n_{\mathrm{c}}\) is the number density in the dense clouds. Because the density differences between \(n_{\mathrm{d}}\) and \(n_{\mathrm{c}}\) are on the order of 5, no strong thermal X-rays are expected from RXJ1713 under the clumpy ISM distribution. The recently discovered thermal X-ray line emission from RXJ1713 originates from the supernova ejecta (Katsuda et al. 2015), and hence it is consistent with the shock–cloud interaction model.

Inoue et al. (2012a) studied \(\gamma \)-ray production in the shock–cloud interaction in RXJ1713 based on MHD numerical simulations, and they showed that a hard spectrum (\(\Gamma _{\mathrm{GeV}} = 1.5\) in the GeV band) like that of leptonic \(\gamma \)-rays is well reproduced in the hadronic scenario with a clumpy ISM. Within the SNR shell, the dense cores are exposed to cosmic-ray protons, and the p–p interactions inside the cores produce \(\gamma \)-rays. The p–p reactions thus lead to a spatial correspondence between the \(\gamma \)-rays and the ISM mass.

The most essential argument of Inoue et al. (2012a) is taking into account the penetration depths of cosmic-ray protons into dense cores, which are limited by the amplified turbulent magnetic field around the dense cores (Inoue et al. 2012a: see also Zirakashvili and Aharonian 2010). This tends to dilute the \(\gamma \)-ray–ISM correspondence. Equation 2 gives an expression for the penetration depth as a function of the magnetic field from Inoue et al. (2012a),

$$\begin{aligned} l_{\mathrm{pd}} = 0.1\; \eta ^{0.5}\; (E / 10\;\mathrm{TeV})^{0.5}\; (B / 100\;\upmu \mbox{G})^{-0.5} \\ (t_{\mathrm{age}} / 1000\;\mathrm{yr})^{0.5} \ (\mathrm{pc}), \end{aligned}$$
(2)

where \(E\) is the particle energy, \(t_{\mathrm{age}}\) is the age of the SNR, and the cosmic-ray energy spectrum is assumed to have the form \(N(E)dE \propto E^{-p} dE\) above the critical energy for \(\pi ^{0}\) creation and below the maximum energy of the accelerated nuclei. Here, \(\eta \) is the degree of randomness of the turbulent magnetic field, which is ∼1 around the cores (e.g., Uchiyama et al. 2007). This indicates that the cosmic-ray protons can penetrate on the order of 0.1 pc into the surface layers of the cores in; this is smaller than the typical size of 1 pc of the CO cores in RXJ1713. Consequently, low-energy cosmic-ray protons cannot penetrate deeply into the dense regions of the cores, and the \(\gamma \)-ray spectrum is modified accordingly (Inoue et al. 2012a; Gabici and Aharonian 2014; Celli et al. 2019; Inoue 2019). The ISM mass that interacts with the cosmic rays is proportional to the cosmic-ray energy \(E\): \(M(E) \propto l_{\mathrm{pd}}(E) \propto E^{1/2}\), if we assume a cosmic-ray distribution with the spectral index \(p = 2\) for high-energy nuclei, as in conventional DSA theory (Inoue et al. 2012a). The increase in \(M\) with \(E\) produces a hard \(\gamma \)-ray spectrum (\(\Gamma { \approx } 1.5\)), similar to the leptonic spectrum produced in a uniform low-density ISM. Figure 9 shows the \(\gamma \)-ray SED of RXJ1713, which can be well reproduced by the hadronic scenario when the penetration depths of the cosmic-ray protons are considered (Celli et al. 2019). In addition, the penetration depth is expected to cause local anticorrelations between the \(\gamma \)-rays and the ISM density on the 0.1-pc scales, which may become directly observable in the future using high-resolution \(\gamma \)-ray observations with the Cherenkov Telescope Array (CTA).

Fig. 9
figure 9

Hadronic models for the \(\gamma \)-ray SEDs of RX J1713.7−3946. The hadronic models (solid lines) refer to a configuration with a magnetic field strength inside the gas clump of \(B_{\mathrm{c}}\) = 1 μG (black) or to \(B_{\mathrm{c}}\) = 10 μ (blue). The field in the skin of the gas clump is fixed at \(B_{\mathrm{s}}\) = 100 μG in both models. The bottom pane shows the residuals for the data/model. From Celli et al. (2019) with permission

Full size image

3.5 Correspondence between the ISM and \(\gamma \)-rays

3.5.1 ISM–\(\gamma \)-ray correspondence

In early works (Aharonian et al. 2006; Fukui 2008), only CO clouds were compared with the \(\gamma \)-ray distributions. Fukui et al. (2012) made a comprehensive comparison between the ISM (including both molecular and atomic hydrogen) and the TeV \(\gamma \)-rays observed using H.E.S.S. (Aharonian et al. 2007) aiming to test the spatial correspondence. A new feature of this work was that the Hi emission was included in the comparison, and the amount of the Hi gas was found to be comparable with the molecular mass. Figure 10 shows that the azimuthal distributions of the two quantities have good spatial correspondence at an effective resolution of 3 pc, which is limited by the H.E.S.S. observations. These results lend support to the hadronic scenario (see below for further discussion). The 3-pc resolution is considerably larger than the typical penetration depth of ∼0.1 pc and the effect of the exclusion of low-energy cosmic rays from the dense region is not significant in the plot. The good spatial correspondence between the total ISM protons and \(\gamma \)-rays shows for the first time that the necessary condition for the hadronic origin of the \(\gamma \)-rays, which was not clear in Aharonian et al. (2006), is satisfied. Subsequent detailed comparative studies among the ISM, gamma-rays, and synchrotron X-rays supported this result by separating the hadronic and leptonic gamma-rays quantitively (Fukui et al. 2021).

Fig. 10
figure 10

(a) Distribution of the total proton column density \(N_{\mathrm{p}}\)(H2+Hi) toward RX J1713.7−3946 superposed on the TeV \(\gamma \)-ray contours. The lowest contour and the contour intervals are 20 and 10 smoothed counts, respectively. (b) Azimuthal distributions of the proton column densities of molecular hydrogen \(N_{\mathrm{p}}\)(H2), atomic hydrogen \(N_{\mathrm{p}}\)(Hi) and \(N_{\mathrm{p}}\)(H2+Hi), and \(\gamma \)-rays between the two elliptical rings shown in (a). The same plot for the area inside the inner ring is also shown on the right-hand side of (b). From Fukui et al. (2012), reproduced by permission of the AAS

Full size image

3.5.2 Cosmic-ray energy

The total cosmic-ray energy \(W_{\mathrm{tot}}\) is estimated from equation 3 below (from Aharonian et al. 2006) to be \({\sim}10^{48}\) erg, which corresponds to 0.1% of the total energy of a typical supernova explosion.

$$\begin{aligned} {W_{\mathrm{tot}} \sim 1-3 \times 10^{50} \biggl( \frac{d}{1\, \mathrm{kpc}} \biggr)^{2} \biggl( \frac{n}{1\, \mathrm{cm^{-3}}}\biggr)^{-1}} \ {\mathrm{(erg)},} \end{aligned}$$
(3)

where \(d\) is the 1-kpc distance to the source and \(n\) is the average density of interstellar protons associated with the SNR. This value is unlikely to represent the total cosmic-ray energy involved in RXJ1713. The number of dense cores in the shell is around 20, each with a typical radius of 1 pc (Sano et al. 2013). For the volume of the shell, which has an 8-pc radius, the volume filling factor of the cores is estimated to be ∼10%, implying that about 1/10 of the cosmic-ray protons are interacting with the ISM protons. This means that the total energy of cosmic rays is 10 times larger than the estimated value of \(W_{\mathrm{tot}}\). The time evolution must be taken into account further, and we expect the maximum energy of the \(\gamma \)-rays to shift to lower values with 10-fold increase in the total energy. The total cosmic-ray energy in the middle aged SNRs such as W44 and W28 is estimated to be 1049 erg in the hadronic scenario (e.g., Yoshiike et al. 2013, 2017). In addition, it is likely that escaping cosmic rays can be more significant than accelerated cosmic-rays inside the SNR. In W44, Uchiyama et al. (2012) showed the nearby clouds outside the SNR are \(\gamma \)-ray bright, and may involve a significant amount of cosmic rays on the order of (0.3–3)\(\times 10^{50}\) erg. This value corresponds to ∼10% of the typical kinematic energy released by a supernova explosion (∼1051 erg). Therefore, SNRs such as RXJ1713 may substantially contribute to the Galactic cosmic-ray energy budget because on average 10% of the total explosion energy is needed to sustain the Galactic cosmic-ray flux (see Gabici 2013 and references therein).

3.5.3 Hadronic vs. leptonic scenarios

The major process responsible for \(\gamma \)-ray production in RXJ1713 is still under debate. The similar shell-like distributions in both X-rays and \(\gamma \)-rays make two scenarios possible. In the leptonic scenario, cosmic-ray electrons are responsible for both \(\gamma \)-ray and X-ray production. In the hadronic scenario, the shell-like ISM distribution produces shell-like \(\gamma \)-rays and the shock–cloud interaction produces a shell-like distribution owing to the enhanced nonthermal X-rays around the dense cores. A possible difference between the two is that the shock–cloud interaction causes spatial offsets on the order of 0.1–0.4 pc between the cores and the \(\gamma \)-ray peaks, which are still below the current limit of \(\gamma \)-ray resolution.

Broad-band fitting of the \(\gamma \)-ray spectrum can be reconciled with the both scenarios according to a recent analysis by H.E.S.S. Collaboration (2018). They argued that both the leptonic and hadronic scenarios can explain the \(\gamma \)-ray observations obtained by H.E.S.S., based on a spectral analysis of the \(\gamma \)-ray observations from 2003 to 2018. In this work, the H.E.S.S. \(\gamma \)-rays and Suzaku X-rays were combined and the broad-band energy spectrum from 1 keV to 10 TeV was presented. The spatial resolution was conservatively set to 3.2 pc, which was adopted for the Suzaku analysis (Tanaka et al. 2008), while the new H.E.S.S. data achieved a nominal resolution of 0.6 pc. In the leptonic scenario, the cosmic-ray electrons are responsible for the \(\gamma \)-ray production. The electrons are cooled by synchrotron energy loss, so the magnetic field has to be small (∼10 μG) for the leptonic scenario to work. In shock–cloud interactions, the field may become as strong as 0.1–1 mG locally, although it may not be so strong over a large volume. H.E.S.S. Collaboration (2018) did not make a comparison with the ISM, which remains to be done in the future.

4 RCW 86

4.1 X-rays and the ISM

The Type Ia SNR RCW 86 is located at (\(l\), \(b\)) = (315.4, −2.3) and it exhibits both thermal and nonthermal X-rays (Fig. 11a). The distance has been determined to be ∼2.5 kpc using several methods (e.g., Westerlund 1969; Rosado et al. 1996; Helder et al. 2013), and its ∼1800 yr age is based on a possible counterpart to the historical supernova event SN 185 recorded in Chinese literature in 185 AD (Stephenson and Green 2002; Smith 1997; Dickel et al. 2001; Vink et al. 2006). The ISM in the velocity range from −46 to −28 km s−1 is associated with the SNR, as found from the spatial correlation with the X-rays (Ajello et al. 2016; Sano et al. 2017a). RCW 86 is close to a CO cloud, that is part of the supershell GS 314.8−0.1−34 identified in CO emission by Matsunaga et al. (2001). Sano et al. (2017a) compared the X-rays with CO and Hi for two energy bands, 0.5–1.0 and 2.0–4.5 keV, where the low energy band is dominated by thermal X-rays and the high energy band by nonthermal X-rays (e.g., Rho et al. 2002; Ajello et al. 2016, see also Fig. 11a).

Fig. 11
figure 11

(a) RGB image of RCW 86 observed by \(\mathit{XMM}\)-\(\mathit{Newton}\). Red, green, and blue represent the X-ray energy bands, 0.5–1.0, 1.0–2.0, and 2.0–4.5 keV, respectively. The white solid lines outline the region observed. The contours represent the radio continuum centered at a frequency 843 MHz. (b) RGB image of RCW 86 and its surroundings. Here, red, blue, and green represent the XMM-Newton X-rays (0.5–4.5 keV), NANTEN2 12CO(\(J\) = 2–1), and the ATCA & Parkes Hi map, respectively. The velocity ranges of CO and Hi are from −46.0 to −28.0 km s−1. The contours indicate the Hi integrated intensity. From Sano et al. (2017a) with permission

Full size image

The ISM surrounding the SNR is dominated by Hi. The overall distribution of both Hi and CO is shell like, and the X-rays are distributed within the cavity (Fig. 11b). In addition, extended weak Hi emission is found inside the shell. According to Sano et al. (2017a), the thermal X-ray filaments with a typical thickness of 1 pc show tight correlations with Hi and CO (Fig. 12), while the nonthermal X-rays are associated with low-density regions of Hi, with densities around \(4 \times 10^{21}\) cm−2. The ISM density in RCW 86 is lower than that in RXJ1713, which is consistent with a Type Ia SNR. Only a CO cloud with a size of 5 pc × 3 pc interacts with the shock, as verified by the high temperature estimated from the high ratio of the two CO transitions, CO \(J\) = 2–1 and 1–0. No heating sources such as IRAS / AKARI infrared point sources or high-mass stars exist in these regions, except for the shock front of RCW 86 (Sano et al. 2017a).

Fig. 12
figure 12

Top left: Thermal X-ray image of RCW 86 in the energy band 0.5–1.0 keV. Other panels: Maps of Hi (right panels, North and Southwest) and 12CO(\(J\) = 2–1) (bottom-left panel, East) obtained with ATCA \(\&\) Parkes and NANTEN2 (color scale) superposed on the X-ray contours. From Sano et al. (2017a) with permission

Full size image

Because the thermal X-rays are produced by shock heating of neutral gas with a density less than 100 cm−3, the correlation with the Hi gas is a natural outcome of the shock heating. The critical difference between RCW 86 and RXJ1713 is probably the total mass of Hi gas inside the SNR shell. For RXJ1713, dense molecular clouds with densities greater than \(10^{4}\) cm−3 remain without being swept up by the SNR shocks. The diffuse and intercloud Hi gas were completely evacuated by strong stellar winds, and then creating the very low-density cavity and the swept-up dense wall (e.g., Fukui et al. 2003; Moriguchi et al. 2005; Sano et al. 2010; Maxted et al. 2013). Therefore, we did not find any Hi envelopes around the dense molecular clouds in RXJ1713, and strong thermal X-rays were not detected at the peripheries (see also Sect. 3.4). On the other hand, RCW 86 shows an Hi envelope surrounding the molecular cloud in the east, and a large amount of diffuse Hi gas remains inside the SNR (see Fig. 11b). This difference is explained by the hypothesis that accretion winds (or disk winds) from the progenitor system of RCW 86 are weaker than the stellar winds from the high-mass progenitor of RXJ1713.

Sano et al. (2017a) also studied the association of the ISM with nonthermal X-rays. Figure 13 shows nonthermal X-ray images superposed on the CO and cold Hi contours toward the east and southwest regions. The CO and Hi clumps are limb-brightened in nonthermal X-rays, which are dominated by synchrotron radiation. This indicates that shock–cloud interactions occur in the surface layers of dense clumps, where the magnetic filed is strongly amplified. RCW 86 has no clumpy CO inside the shell, and the gas density seems to be smoother than in RXJ1713. The authors interpreted the lower ISM mass to be related to the older evolutionary state of the region after star formation because RCW 86 is part of a molecular supershell with an age of ∼2 Myr (Matsunaga et al. 2001). By contrast, in RXJ1713 we see ongoing star formation within the SNR shell, as evidenced by the protostellar outflow in CO peak “C” and a few more YSO candidates (Sano et al. 2010).

Fig. 13
figure 13

Nonthermal X-ray images (\(E\): 2.0–4.5 keV) toward the (a) East and (b) Southwest. Superposed contours indicate (a) CO and (b) cold Hi. From Sano et al. (2017a) with permission

Full size image

4.2 The origin of the \(\gamma \)-rays

Figures 14(a) and 14(b) show overlays of the synchrotron X-ray and radio continuum distributions on the TeV \(\gamma \)-ray contours, where the lowest \(\gamma \)-ray contour level of \(\gamma \)-ray corresponds to a significance level of \({\sim}5\sigma \) (see also Fig. 1 of H.E.S.S. Collaboration 2018). In these images, the synchrotron emission appears as a nearly circular shell, whereas the \(\gamma \)-ray emission seems to deviate from this shape, particularly toward the West, where the \(\gamma \)-ray emission appears to be shifted toward the SNR interior. If true, this would argue against the leptonic scenario. However, by employing a quantitative analysis that uses radial profiles in the different regions of interest, H.E.S.S. Collaboration (2018) showed that within statistical errors, there is no significant difference between the TeV \(\gamma \)-rays and the X-ray synchrotron emission. Future \(\gamma \)-ray observations with high sensitivity and high angular resolution are needed to confirm the spatial match or the difference between the two images.

Fig. 14
figure 14

(a–b) Synchrotron radiation images from (a) XMM-Newton X-rays (\(E\): 2–5 keV) and (b) Molonglo Observatory Synthesis Telescope (MOST) 843 MHz radio continuum, both toward RCW 86. The superposed solid contours represent the H.E.S.S. TeV \(\gamma \)-rays. The dashed contours indicate the intensity distributions for each map smoothed with the same point-spread function as the \(\gamma \)-ray image. (c) Hi map of RCW 86 obtained with ATCA \(\&\) Parkes superposed on the TeV \(\gamma \)-ray contours. (d) Distribution of \(N_{\mathrm{p}}\)(H2 + Hi) toward RCW 86 superposed on the TeV \(\gamma \)-ray contours. (e) Azimuthal distributions of \(\gamma \)-rays, \(N_{\mathrm{p}}\)(H2), \(N_{\mathrm{p}}\)(Hi), and \(N_{\mathrm{p}}\)(H\(_{2}+\)Hi), which are averaged over every 30∘ in azimuthal angle within the region enclosed by the white circle in Fig. 14(d). From Sano et al. (2019a), reproduced by permission of the AAS

Full size image

A spatial comparison between the ISM and \(\gamma \)-rays therefore remains essential for understanding the origin of the \(\gamma \)-rays. In Fig. 14(c) the TeV \(\gamma \)-ray distribution is overlaid on the Hi, and the comparison shows good spatial correspondence. Figure 14e shows a plot of the azimuthal distributions of the \(\gamma \)-rays and of the total proton column density \(N_{\mathrm{p}}\)(H\(_{2} + \)Hi), derived using both the Hi and the CO. The \(\gamma \)-rays show good correspondence with the Hi distribution, while no enhancement is seen toward the CO cloud (Figs. 14d and 14e in the azimuth angles from \(-90^{\circ }\) to \(0^{\circ }\)). The poor correlation with the CO is ascribed to the effect of cosmic-ray exclusion due to the small penetration depths (less than 1 pc; Sano et al. 2019a, see Sect. 3.4). More precisely, the radii of the shock-interacting molecular clouds constitute a key difference between RXJ1713 and RCW 86. For RXJ1713, the cloud radius is ∼1 pc or less; hence, the molecular clouds play a significant role as targets of cosmic-ray protons. However, the cloud radius of RCW 86 is more than 5 pc (Sano et al. 2017a). Additionally, the weak disk winds from the progenitor system of RCW 86 were not able to strip off the Hi envelope from the surface of molecular cloud. Therefore, cosmic-ray protons accelerated in RCW 86 can interact only with the low-density Hi gas on the surface of the molecular cloud.

For the azimuthal angles from \(-180^{\circ }\) to \(-120^{\circ }\), we find \(\gamma \)-ray excesses relative to the total proton column density. We see a hint of the \(\gamma \)-ray excess by a factor of 1.2, which may be ascribed to the contribution from a leptonic \(\gamma \)-ray component, because bright synchrotron radiation is detected toward the SE region. The possible contribution of leptonic \(\gamma \)-rays to the total \(\gamma \)-rays may be as low as 5% if the arguments given above are correct. The total cosmic-ray energy in RCW 86 is estimated to be \({\sim}1 \times 10^{48}\) erg. A similar argument to the case of RXJ1713 is perhaps applicable here, and the energy is a typical one for a young SNR with an age of ∼2000 yr. Considering the lower sensitivity of the current H.E.S.S. data for RCW 86, we do not discuss the \(\gamma \)-ray origin further until more sensitive \(\gamma \)-ray data become available.

5 Summary and future prospects

We have reviewed the ISM associated with young supernova remnants and its connection with the observed X-rays and \(\gamma \)-rays. These high-energy radiations are closely related to cosmic rays that are probably accelerated via DSA. The present focus was on the young, bright TeV \(\gamma \)-ray SNRs RXJ1713 and RCW 86. Shock–cloud interactions in a clumpy ISM play an important role in regulating the X-ray properties as well as the \(\gamma \)-ray production.

The pursuit of the ISM in SNRs is being extended to nearby galaxies, including the Large Magellanic Cloud (LMC) and the Small Magellanic Cloud (SMC). The ISM associated with the SNRs has been identified in many of them, providing promising candidates for further studies on the role of the ISM (Sano et al. 2015b, 2017b,c, 2018, 2019b,c, 2020a; Yamane et al. 2018, 2021; Roper et al. 2018; Alsaberi et al. 2019). The smaller contamination along the line of sight is advantageous for identifying molecular clouds associated with the SNRs and for using optical and infrared datasets without strong stellar absorption. Although the Magellanic SNRs lie at 50–60 times larger distances from us than compared with the distance of RXJ1713 from us (e.g., Hilditch et al. 2005; Pietrzyński et al. 2013), ALMA’s unprecedented sensitivity and spatial resolution have enabled us to resolve spatially the cloud-scale structures of the ISM associated with the Magellanic SNRs.

Figure 15 shows CO mapping results toward the five Magellanic SNRs and a superbubble that are bright in X-rays. The molecular clouds lie nicely along, or embedded within, the X-ray shells, indicating that the shock–cloud interactions have occurred. 30 Doradus C, the brightest X- and \(\gamma \)-ray superbubble in the Local group (see Kavanagh 2020 and references therein), provides one of the best laboratories for studying shock–cloud interactions because of its large apparent diameter and very bright TeV \(\gamma \)-rays (H.E.S.S. Collaboration 2015). Future observations using ALMA, Chandra, and the CTA will allow us to study in detail not only shock–cloud interactions but also the hadronic \(\gamma \)-ray radiation from the Magellanic SNRs.

Fig. 15
figure 15

CO results toward the five Magellanic SNRs N63A (Sano et al. 2019b), N49 (Yamane et al. 2018), RX J0046.5−7308 (Sano et al. 2019c), N132D (Sano et al. 2020a), N103B (Sano et al. 2018), and the Magellanic superbubble 30 Doradus C (Yamane et al. 2021). The CO observations of the LMC SNRs N63A, N49, N132D, N103B, and the LMC superbubble 30 Doradus C were obtained using ALMA, while the SMC SNR RX J0046.5−7308 was observed using the Atacama Submillimeter Telescope Experiment (ASTE). From Sano et al. (2018, 2019b,c, 2020a) and Yamane et al. (2018, 2021), reproduced by permission from the AAS

Full size image

We summarize the main points of this review as follows:

  1. 1.

    Spatial comparisons between X-rays and the ISM provide a powerful approach for determining the distances to SNRs, which is otherwise difficult for SNRs in the Galactic plane because of heavy obscuration. For RXJ1713, the ISM delineates the outer boundary of the nonthermal X-ray shell at pc scales, while X-ray bright spots (or filaments) are anti-correlated with the ISM clouds/cloudlets at small-scales (sub-pc scales), allowing a robust distance estimate. RCW 86 is an easier case, for which distance has been determined accurately owing to its relatively high Galactic latitude.

  2. 2.

    When we consider the ISM surrounding the SNRs, it is important to recognize that the neutral ISM is characterized by highly clumped distributions, with clumps that are much denser than the nominal uniform density of 1 cm−3 assumed in many previous works. The ISM shell in a SNR surrounds a low-density cavity, and the shell comprises a layer that includes tiny, dense clumps with a small volume filling factor. Consequently, an SNR has a large low-density volume, where DSA works, as well as dense clumps, which interact with the cosmic rays and shocks. Magnetohydrodynamical numerical simulations of shock–cloud interactions reveal that the dense clumps cause strong deformations of the shock fronts. The deformations produce highly turbulent velocity distributions, which entangle and amplify the magnetic field to mG levels (Inoue et al. 2012a). In most parts of the low-density cavity, however, the shocks propagate with little deceleration or deformation.

  3. 3.

    The clumpy ISM picture requires modification of the previous scenarios that were based on a homogeneous ISM. The observed \(\gamma \)-ray spectrum of RXJ1713 is hard, according to the Fermi collaboration, and the spectrum has previous been interpreted in terms of the leptonic scenario in a uniform ISM picture. However, we have shown that a the clumpy ISM can explain the hard spectrum in the GeV band (\(\Gamma _{\mathrm{GeV}}\) = 1.5) equally well because the penetration of cosmic rays into dense clumps is inhibited by diffusive scattering due to the turbulent magnetic field. We have discussed the implications of a clumpy ISM for interpreting the observations of several SNRs and have outlined the future prospects for this field.