Leaked GeV CRs from a Broken Shell: Explaining 9 Years of Fermi-LAT Data of SNR W28

We investigate the morphology of the W28 complex region in the 10–200 GeV regime. The test-statistic (TS) map of this field is shown in Figure 1, where all 3FGL catalog sources except for 3FGL J1801.3–2326e (the northeastern part of W28), 3FGL J1800.8–2402, 3FGL J1758.8–2402, and 3FGL J1758.8–2346 are subtracted. The peak detection significance is ∼20σ and is in the northeastern part.

Zoom In Zoom Out Reset image size

In the source model, we proceeded to replace 3FGL J1800.8–2402 and 3FGL J1758.8–2402 with three point sources at the positions of 240 A, B, and C, respectively. In order to revise the morphology of 3FGL J1801.3–2326e (originally a uniform disk with a 039 radius in the northeastern part), we followed the scheme of a likelihood ratio test adopted by Yeung et al. (2016, 2017). Since its centroid is exactly at the center of 3FGL J1801.3–2326e, we varied only the radius of the extension and left the center unchanged. We assigned it a simple power law (PL). The −ln(likelihood) of different radii in 10–200 GeV is tabulated in Table 1, and the most likely radius is determined to be 0345. We therefore adopted this morphology for the northeastern part of W28 in the subsequent analysis, and we refer to it as Fermi J1801.4–2326 to avoid confusion with 3FGL J1801.3–2326e.

We furthermore re-created the TS map with all sources in the revised model except for 3FGL J1758.8–2346 (the western clump) subtracted, and this is presented in Figure 2. We thus revised the position of 3FGL J1758.8–2346 to be (26947917, −23820494)_J2000, which is the centroid determined on this map, and we renamed it "Source-W." Here, we finalized the source model for the spectral analysis.

Zoom In Zoom Out Reset image size

Motivated by the LAT spectra of 240 B and C presented in Section 2.3, we also investigated their morphology in lower energy bands. We created TS maps in which all sources except 240 B and C are subtracted (Figure 3) in 1–50 GeV (left) and 0.5–1 GeV (right). We adopted PSF3 data to optimize the spatial resolution.

Zoom In Zoom Out Reset image size

In 1–50 GeV, the excess shown in the map is coincident with 240 B and C, but the excesses from these two spatial components are hardly resolved. In 0.5–1 GeV, the excesses from 240 B and C, if any, are buried below the PSF wing of a brighter blob (the southern blob) whose centroid is at (26973584, −24268803)_J2000. The detection significance of this southern blob is ≳30σ.

The diffuse background that is mainly caused by the CR sea interacting with the interstellar medium is particularly high at <10 GeV band. As seen in Figure 4, the southern blob is located inside a bright background region, therefore it is also possible that the southern blob is caused by an imperfect background reduction.

Zoom In Zoom Out Reset image size

For a CC SN with a minimum progenitor mass (8M_⊙), the progenitor wind bubble can be neglected. Therefore, our SNR can directly expand into the ICM from the beginning. Our modeling results are only sensitive to the pre-SN environment (e.g., the progenitor wind bubble) but not to the ejecta mass, therefore a type Ia SNR scenario would deliver results similar to those shown below.

In our work, we adopt the analytical solutions for the ejecta-dominated stage and the Sedov–Taylor stage from Chevalier (1982), Nadezhin (1985), and Ostriker & McKee (1988), Bisnovatyi-Kogan & Silich (1995), and Ptuskin & Zirakashvili (2005), respectively. When the shock is further slowed, we adopt the analytical solution for the pressure-driven snowplow (PDS) stage derived by Cioffi et al. (1988). As shown in Figure 7, inside a homogeneous ICM with density ${n}_{\mathrm{ISM}}=6\,{\rm{H}}\,{\mathrm{cm}}^{-3}$, the SNR spends its first ∼650 years (∼2.5 pc) in the ejecta-dominated stage. At ∼14 kyr (∼9.3 pc), it completes the Sedov–Taylor stage and enters the pressure-driven snowplow stage. When it finally reaches 13 pc at 37 kyr, the present shock velocity is 110 km s⁻¹. The escape energy of the run-away CRs, which is marked as blue lines in Figure 7, is derived from the CR acceleration theory below.

Zoom In Zoom Out Reset image size

Cosmic-ray acceleration in the SNR is well known as the first-order Fermi acceleration, which normally provides a PL CR spectrum at the shock with an index Γ_CR ∼ −2.0 and a cutoff at the escape energy E_max. To explain the super-TeV CRs observed in many young SNRs, Bell (2004) and Zirakashvili & Ptuskin (2008) have developed the acceleration theory of nonresonant streaming instability, in which the magnetic turbulence at the shock upstream is quickly amplified by the CR streaming, and is finally able to boost the escape energy up to hundreds of TeV in a young SNR.

In a strong shock, only CRs with an energy higher than E_max are able to escape from the upstream and propagate to the nearby environment with a flux J_out, while other CRs that are trapped at the shock follow a PL spectrum with an index Γ_CR ∼ −2.0. For a young SNR with a high shock velocity ${v}_{\mathrm{SNR}}\gtrsim 1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$, only these super-TeV CRs are expected to be released. These early escaped super-TeV CRs could explain the TeV emission from the MCs around intermediate-age SNRs well, e.g., SNR W28 (Gabici et al. 2010), and even around young SNRs, e.g., SNR HESS J1731–347 (Cui et al. 2016). With the Fermi-LAT observational study by Abdo et al. (2010), the discovered GeV emission from W28-North, which peaks at ≲0.3 GeV, indicates alternative sources of CRs. Sticking to the run-away CR explanation, one would need to gradually decrease the escape energy to ≲1 GeV, mainly by decreasing the shock velocity, decreasing the magnetic field in the upstream, and/or increasing the density of the neutrals in the upstream.

When the shock velocity v_SNR and the density of the nearby circumstellar medium are obtained, the CR acceleration processes can be calculated through the acceleration theory of nonresonant CR streaming from Zirakashvili & Ptuskin (2008). An analytical approximation of this theory derived by Zirakashvili & Ptuskin (2008) is adopted in our work. It provides the escape CR flux J_out, the CR density at shock n₀, and most importantly, the escape energy E_max; see more details in Zirakashvili & Ptuskin (2008) and Cui et al. (2016). A magnetic field of B₀ = 5 μG and an initial magnetic fluctuation (not amplified by the CR streaming yet) of B_b = 7%B₀ in the ICM are used to calculate the acceleration process.

Knowing the run-away CR spectrum J_out at any given time, we can integrate the entire SNR history as well as the whole surface of the SNR, and eventually, we are able to derive the run-away CR density analytically at any location with a given diffusion coefficient; see details in Section 2.4.1 of Cui et al. (2016). By lowering the escape energy to 10 GeV, which is presented in our damping model (see Figure 8 and text below), a significant amount of GeV CRs can also be released via this run-away process during the late stage of the SNR evolution.

Zoom In Zoom Out Reset image size

The damping of the magnetic waves upstream of the SNR, which is due to the presence of neutral atoms, becomes important in intermediate-age/old SNRs (Zirakashvili & Ptuskin 2017). Hence, in addition to the shock-MC encounter part of the SNR, W28-North, the other parts of this intermediate-age SNR are also likely to release GeV CRs. Both the damping model and the non-damping model are explored in our work. In the non-damping model, the escape energy follows the analytical solution of the nonresonant CR streaming theory mentioned above. In the damping model, this same analytical solution only works during the early stage of the SNR, when the shock velocity drops to ∼1000 km s⁻¹ (R_SNR = 4 pc). The estimate of the escape energy in partially ionized medium of ${E}_{\max }={v}_{\mathrm{SNR}}^{3}{n}_{{\rm{H}}}^{0.5}{n}_{{\rm{n}}}^{-1}$ by O'C Drury et al. (1996) is adopted, which is in accordance with the simulation work by Zirakashvili & Ptuskin (2017). This number ∼1000 km s⁻¹ also gives a smooth transition from the analytical solution by Zirakashvili & Ptuskin (2008) to the estimate by O'C Drury et al. (1996). Here n_n is the number density of neutrals, which is set to be 0.3 cm⁻³ in our damping model.

For intermediate-age/old SNRs, the CR acceleration efficiency η is expected to decrease at a low Mach number shock (M_s ≲ 50); see, e.g., the hydrodynamic estimation of the acceleration efficiency by Voelk et al. (1984). We note that the acceleration efficiency η used in this work is the ratio between the energy flux of run-away CRs and the kinetic energy flux of incoming gas onto the shock. In our work, both the sound speed and Alvén speed in ICM are set to be 15 km s⁻¹ (Chevalier 2005). Here we adopt a relationship of $\eta ={\eta }_{0}\exp (-{({M}_{0}/{M}_{{\rm{s}}})}^{{{\rm{\Gamma }}}_{\eta }})$, where η₀ is the acceleration efficiency when the shock velocity is high. This is set to be 0.04 and 0.03 for the damping model and non-damping model, respectively. By choosing Γ_η = 1.2–1.4 and M₀ = 4–5 in our models, this relationship is roughly consistent with the simulation results of 0° < θ < 45° by Caprioli & Spitkovsky (2014), where θ is the angle between the normal of the shock and the magnetic field. When θ ≳ 45°, the acceleration efficiency drops significantly. These acceleration parameters adopted in our model are also constrained by the limitation of the total CR energy; see more information in Section 3.3.5.

The shock-MC encounter time is set at 25 kyr (R_SNR ∼ 11. 5 pc) after the SN. Following the encounter, part (χ) of the SNR shell located at W28-North is rapidly slowed down by the dense MC clump, and only very little gas is shocked (the estimated X-ray emitting mass at W28-North is only ∼1 M_⊙ according to Zhou et al. 2014) before the shock is fully stalled. In fact, the shocked gas mostly comes from the transition region between the diffuse gas and the dense clump. After the shock-MC encounter, a significantly slowed shock will continue to propagate inside the MC clump, but is unable to shock the dense gas there (Gabici & Aharonian 2016). Owing to the high density of the neutral gas, both the magnetic turbulence in the upstream and the downstream will be significantly damped, leaving the CRs confined at and behind this part of the shock, with a total number ∼_χN ready to be released, where N is the total number of CRs trapped inside the SNR. During the SNR evolution history, CRs are continually carried downstream into the interior of the SNR, with a flux J_in ∼ n₀v_SNR/4. Some of these trapped CRs inside the SNR may re-enter the acceleration region and even escape from the upstream. We note that the CRs (E ≲ E_max) confined inside the acceleration region are almost instantly released after the shock-MC encounter, see, e.g., Ohira et al. (2011), while the GeV CRs filling in the inner region of the SNR need to be carried into MC-N by gas flow/turbulence. Assuming a magnetic field of >10 μG inside the SNR, which is consistent with the simulation results by Zirakashvili & Ptuskin (2017), the mean Bohm diffusion distance of a 1 TeV proton is merely ≲0.3 pc after 10 kyr. After the CRs cross the stalled shock and propagate into MC-N and beyond, they will enter an environment with a much higher diffusion coefficient, which is 10% of the Galactic diffusion coefficient in our model.

Simulation work by Zirakashvili & Ptuskin (2012, 2017) has shown that the spatial distribution of the CRs inside the SNR is dependent on the SNR age and CR diffusion coefficient (CR energy). The intermediate-age/old SNRs tend to show a more homogeneous CR distribution inside the SNR, while the CRs inside young SNRs are mostly confined immediately behind the shock. The GeV CRs are more likely to "feel" the gas compression and are more concentrated immediately behind the shock, while super-TeV CRs fill the inside of the SNR more homogeneously. See also the radius profiles of CR pressure and gas density inside different SNRs by Zirakashvili & Ptuskin (2012).

In addition to the CRs in the acceleration region that are immediately released into MC-N after the shock-MC encounter, the GeV CRs located in the inner region of the SNR, which "feel" the gas compression, can also be carried into MC-N by the downstream flows/turbulence. In our model, the CRs released from W28-North (including both the CRs from acceleration region and the CRs carried by flow) are arbitrarily chosen as an instantaneous process. Hence, this release time (12 kyr ago) actually functions like an averaged release time, and the total energy of leaked CRs from W28-North is normalized later in our modeling.

Similar to the spectrum of the CRs confined inside the SNR in the simulation work by Zirakashvili & Ptuskin (2017), the spectrum in our model (also the spectrum of leaked CRs), as seen in Figure 8, is arbitrarily made of two components: first, a low-energy part with a PL index of Γ_CR = −2.0 and a cutoff energy of E_max, which represents the CR spectrum in/near the acceleration region, and second, a high-energy tail extending beyond E_max with a PL index of Γ_CR = −2.3( − 2.5) and a cutoff energy of 10(5)TeV for the damping(non-damping) model. In our models, χ represents the percentage of the SNR shell that is stalled by the MC-N at 25 kyr. After the MC-shock encounter, part (χ) of the shock starts to release all its GeV CRs, while the other 1 − χ GeV CRs will mostly remain behind their local shock. χ is a free parameter here, but it is constrained by the limitation of the total CR energy. By choosing χ = 10%, for instance, we would derive a total CR energy in the SNR of $\sim 4 \% /11 \% \,{{ \mathcal E }}_{51}$ 12 kyr ago in our damping/non-damping model.

At 25 kyr after the SN, the energy of TeV CRs trapped inside the SNR only accounts for ≲7% of the energy of the total trapped CRs in our models. Unlike the GeV CRs, all of the diffuse super-TeV CRs trapped in the entire SNR (rather than 10% of the SNR) are able to be gradually released through the leaking hole of W28-North. In our models, the run-away CRs from the early SNR stages dominate the higher energy γ-ray emission at the MC clumps. Therefore, for simplicity, part (χ) of the TeV CRs is bound together with the GeV CRs and released instantaneously from W28-North.

The spatial distribution of Galactic sea CRs, especially the GeV sea CRs, could be quite inhomogeneous. Naturally, the GeV CRs are likely to be concentrated near the star-forming regions, e.g., the spiral arms. Through studying the Fermi-LAT data and the gas density in the entire Galaxy, Acero et al. (2016) have derived a rough radius profile of CR density/spectrum in the Galactic plane. Considering that SNR W28 lies inside the Galactic plane with a distance to the Galactic center of ∼5 kpc, using this radius profile, one can obtain the spectrum of its local CR sea (CR spectrum index Γ_CR,sea ∼ −2.55 to −2.72, energy density ${U}_{\mathrm{CR},\mathrm{sea}}\sim 1.1\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$), which is similar to the one detected on Earth. However, in our models, we explore scenarios with CR sea densities lower than the averaged local one mentioned above, in order to better fit the ≲1 GeV observations (especially for 240A).

In the reproduction of the γ-ray emission at MC-N, MC-A, B, and C, the diffusion coefficient for the entire space outside the SNR is fixed as D(E) = 10²⁷(E/10 GeV)^δ cm² s⁻¹, δ = 0.5, mostly based on the requirement of fitting the GeV data of W28-North and the TeV data of all the γ-ray sources. Thus the fine-tuning of the distances between these MC clumps and the CR sources is performed for a better fitting result, see Table 4. The three-dimensional distances between the MC clumps and the CR sources in our model satisfy the triangle relationship among the three objects: W28-North, SNR center, and each MC clump.

Table 4.  The Distances (pc) between MCs and CR Sources

	MC-N (5 M4a)	MC-A (4.3 M4)	MC-B (6 M4)	MC-C (2 M4)
Damping
SNR center	13	35	31	27
W28-North	0∼1	37	29	28
Non-damping
SNR center	13	35	28	27
W28-North	0∼1	33	26	25

Note.

^aHere ${M}_{4}={10}^{4}{M}_{\odot }$.

Download table as:  ASCIITypeset image

For GeV CRs, their spatial distribution are sensitive to the distances to the CR sources. The distances of MC-A, B, and C to W28-North are about 25–35 pc in our model, and 25 pc is approximately the projected distance between W28-North and MC-B,C. Our model faces difficulties to efficiently bring ≲10 GeV CRs from the SNR to MC-A and MC-B, which have shown γ-ray spectra peaks at 1.4 GeV and 2.8 GeV, respectively. Adding the CR sea could help moving the peaks of the reproduced γ-ray emission to lower energy; see the dotted and solid lines in Figure 9.

In the case of 240 A, when we place MC-A much closer to the SNR, e.g., 20 pc (the same effect as using a higher diffusion coefficient in the entire space), we could move the peak to a lower energy, but the reproduced 10–100 GeV flux will be too high to explain the observation as well. This overproduction of 10–100 GeV flux also occurs for MC-B and C. By modifying the diffusion coefficient, one solution for this situation is to introduce a diffusion coefficient with a smaller PL index (compared to 0.5 for a Galactic one), which allows the 1–100 GeV CRs to diffuse relative faster. Another solution would be introducing a diffusion coefficient spectrum wth BKPLs: the PL index of the diffusion coefficient of 1–100 GeV CRs can be smaller than 0.5. Interestingly, the derived diffusion coefficients inside a Kolmogorov turbulence do show this type of a BKPLs, see, e.g., the numerical simulation results by Casse et al. (2002) and Fatuzzo et al. (2010). However, a PL index of the diffusion coefficient smaller than 0.5 would require a softer spectrum of the leaked GeV CRs (Γ_CR < −2.0) to explain the GeV spectrum of W28-North. Furthermore, an anisotropic or/and inhomogeneous diffusion environment could also help our modeling, e.g., fast-diffusing tubes based on the large magnetic structure. Overall, these alternative diffusion scenarios, which may need a total modification of our model, including the acceleration and releasing of the GeV–TeV CRs, will be left to future work.

Owing to the long distances from W28-North to MC-A, B, and C (≳20 pc), 1–100 GeV CRs leaked from W28-North are not able to efficiently reach these MC clumps at the present time with the diffusion coefficient used in our model. Hence, the CR sea is expected to dominate the ≲10 GeV γ-ray emission at these three MC clumps. The density of the CR sea (5 kpc from Galactic center) observed by Acero et al. (2016) contradicts our Fermi-LAT discovery; the non-detection of ≲1 GeV γ-rays at MC-A. These observed MC clumps cover a distance to Earth ranging from 2 to 4 kpc, see Aharonian et al. (2008), but placing them farther from Earth (e.g., 4 kpc instead of 2 kpc; power sources other than W28 could be introduced as well to explain the higher energy band) only leads to a situation in which the MC clumps are placed closer to the Galactic center (3–4 kpc), where a CR sea density higher by about three times than is detected on Earth is suggested by Acero et al. (2016). Noticeably, the diffuse background in the Fermi analysis is very important at this low-energy band, and we cannot rule out that our non-detection of ≲1 GeV γ-ray at MC-A could be the result of an excessive background reduction; see our background in Figure 4.

In summary, if the low <1 GeV flux at 240 A is true, an inhomogeneous distribution of the <10 GeV CR sea becomes necessary. The CR sea density adopted in Figure 9 is half of the averaged local one, except for 240 A, for which we adopt a much lower one of 14% of the averaged local one.

As described in the introduction section, considering the ionization evidence found at MC-N, we believe the CRs with energy down to <1 GeV have already leaked through the broken shell at W28-North and dominate the GeV emission there. By lowering the escape energy, run-away CRs can also contribute significantly to the GeV emission at W28-North. Nonetheless, in the damping theory, the ionization ratio in the ICM, which functions as a key factor on the escape energy, is unknown. With an ICM density of ≲10 cm⁻³, an ICM temperature of T ≳ 10⁴ K, which is just above the ionization temperature of H, is needed to confine the dense clumps and to support the cloud against gravitational collapse (Blitz 1993; Chevalier 1999). Therefore, in our damping model, the neutral density of the upstream ICM is set as low as 0.3 cm⁻³, and the final escape energy of the 100 km s⁻³ shock at the present time is 10 GeV. To further decrease the escape energy to 1 GeV and allow the run-away CRs (from shocks other than W28-North) to dominate the <1 GeV γ-ray emission at MC-N, one requires a neutral density ∼3 cm⁻³ in the ICM, following the estimation by O'C Drury et al. (1996). This requirement slightly contradicts the observed H i density upper limits of 1.5–2 cm⁻³ by Velázquez et al. (2002).

In our SNR evolution model of W28, we have only explored a CC SN scenario that delivers a homogeneous pre-SN environment. This scenario has successfully explained the observed shock velocity at the present time and the measured gas density near the SNR. Here we discuss another three hypotheses that have more massive progenitors and more complex pre-SN wind environments.

• 1.  
  A large pre-SN wind bubble, radius R_b > 13 pc. Owing to the low shock velocity at the present and the lack of nonthermal X-ray, SNR W28 probably does not still evolve inside the pre-SN wind bubble, which has a typical density of ∼0.01 cm⁻³.
• 2.  
  A medium pre-SN wind bubble, radius 13 pc > R_b ≳ 10 pc. If SNR W28 were to have encountered the pre-SN bubble shell relative recently, we would more likely see a shell structure of thermal X-rays rather than the observed features: a diffused one in the center and a bright one at W28-North.
• 3.  
  A small pre-SN wind bubble, radius R_b ≲ 10 pc. If SNR W28 has already swept the pre-SN bubble shell and is expanding in the ICM before it enters the radiation-dominated stage, a similar shock velocity (v_SNR ∼ 100 km s⁻¹) at the present time is expected (adopting the same SN energy ${{ \mathcal E }}_{\mathrm{ej}}={{ \mathcal E }}_{51}$) because the final shock velocity during the Sedov stage is mostly sensitive to the total swept-up mass, but not to the radius distribution of the gas; see, e.g., the scenario 8 M_⊙ and scenario 15 M_⊙ in Cui et al. (2016). In this new SNR scenario, which is not explored in this work, the GeV CR leaking process that occurs after the sweeping of the pre-SN bubble shell is expected to be similar to that in our scenario, but the release of super-TeV CRs that mostly occurs in the early stage of the SNR evolution will be different.

Since our SNR evolution history is fixed, the key normalization factor on the total accelerated CR energy becomes the acceleration efficiency, see Table 5. Adapting the η_esc = 0.04/0.03 in our damping/non-damping model, at the shock-MC encounter time (25 kyr), the SNR has already released most of its run-away CRs with a total energy about $9.1 \% /6.6 \% \,{{ \mathcal E }}_{51}$. After parts (χ) of the SNR have encountered MC-N (25 kyr), about $0.4 \% /1.1 \% \,{{ \mathcal E }}_{51}$ CRs are released from W28-North in the damping/non-damping model. If we follow the concept by Zirakashvili & Ptuskin (2012, 2017) that most of the accelerated CRs have escaped the SNR through run-away CRs for a intermediate-age/old SNR, then our total CR energy inside the SNR at 25 kpc should be lower than the one of run-away CRs, which leads to a χ ≳ 4.3%/16.7% in the damping/non-damping model. After the shock-MC collision, only the other 1 − χ part of the shell still releases run-away CRs. χ is set ato be 10%/15% in our damping/non-damping model.

Table 5.  Parameters in the Two Models

Modelsa	χ	η0	M0	Γη	Emaxb (TeV)	ECR,runc (${{ \mathcal E }}_{51}$)	ECR,leakd (${{ \mathcal E }}_{51}$)
Damping	10%	0.04	4	1.4	0.012	9.2%	0.4 %
Non-damping	15%	0.03	5	1.2	0.72	6.7%	1.1%

Notes.

^aThe two models share the same SNR evolution history and diffusion environment. ^bThe escape energy at the present time (37 kyr). ^cThe total energy of run-away CRs at the present time (37 kyr). ^dThe total energy of CRs leaked through W28-North at the shock-MC encounter time (25 kyr).

Download table as:  ASCIITypeset image

Based on the X-ray ionization time at W28-North by Zhou et al. (2014), we adopt an instantaneously releasing time (12 kyr ago from present) to set the CRs free at W28-North. In fact, as described in Section 3.2.3, only CRs from the acceleration region are released instantaneously, while the GeV/TeV CRs from the inner region are gradually injected into MC-N advectively/diffusively. After the shock-MC encounter, the tip of MC-N, which is buried inside the SNR, will keep being blown by a downstream wind (flow) with a velocity of v_flow, where v_flow is slightly lower than the current shock velocity (at 25 kyr, v_flow ≲ 150 km s⁻¹). However, this advection-dominated injection of GeV CRs by the downstream flow will decrease with time, because first, the CR density gradually decreases from the region immediately behind the shock to the inner region of the SNR; see e.g., the simulation result of intermediate-age/old SNRs by Zirakashvili & Ptuskin (2012), whose CR density is halved at 10% R_SNR downstream from the shock front. Second, the flows that are falling onto MC-N will decrease when the gas pressure behind the stalled shock reaches a new balance and the eddies stirred by MC-N decay in several crossing times t_c = L_ed/v_flow, where L_ed is the length of the eddies and should be comparable to the size of the MC-N tip that is embedded in the SNR (${t}_{{\rm{c}}}\sim {10}^{4}\,$ year when L_ed ∼ 1 pc). To fully unveil the leaking process at W28-North, an MHD simulation is required in future work, to obtain the flow/turbulence details as well as the diffusion environments in/around this shock-MC encounter region.

We note that this shock-MC encounter at W28-North could be a more complex process than our simple "colliding onto a flat wall" model, considering that the three-dimensional structure of this MC-N could be very complex. This complex encounter process could last for quite a long period, e.g., with an expanding velocity $\sim 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$, SNR W28 would need ≳10⁴ year to fully swallow MC-N, whose size is ≳1 pc.