Regions of Suppressed Diffusion around Supernova Remnants?

The recent discovery of the so-called TeV halos has attracted much attention. The morphology of the emission requires that the region is characterized by severe suppression of the diffusion coefficient. This finding raises many questions as to its origin: (1) is the suppressed diffusion to be attributed to instabilities induced by the same radiating particles? (2) or does it actually show that the diffusion coefficient is small throughout the disk of the Galaxy? In both cases, one would expect that the surroundings of supernova remnants (SNRs) should also show evidence of a reduced diffusion coefficient since most remnants are located in the disk and are expected to be sites of effective particle acceleration. Should we expect the existence of regions of extended γ-ray emission from these regions as well? Here, we investigate the transport of cosmic rays (CRs) that have escaped from SNRs in order to assess the viability of the idea of having a cocoon of suppressed diffusion around them. A comparison of our results with the γ-ray emission from the regions around HB9 and W28 does not provide solid evidence of reduced diffusivity. However, if indeed the phenomenon of reduced diffusivity occurs around SNRs surrounded by molecular clouds, our calculations show that the effects on the grammage of Galactic CRs can be significant.


INTRODUCTION
The recent discovery of the so-called TeV halos (see e.g., Abeysekara et al. 2017;Sudoh et al. 2019), regions of spatially extended γ-ray emission at energies TeV around selected pulsar wind nebulae (PWNe), has raised a series of questions that might have a deep impact on our conventional picture of the transport of cosmic rays (CRs).The morphology of these halos is such that the emission is best interpreted as inverse Compton scattering (ICS) of the electrons and positrons originated from the PWN, propagating diffusively in a region of size ∼ 50 pc.The diffusion coefficient in the region appears to be suppressed by a factor ∼ 100 − 1000 with respect to the value deduced for the ISM using the standard indicators, such as the secondary to primary ratio (see e.g., Aguilar et al. 2016) and the abundance of unstable isotopes (see e.g., Strong et al. 2007).A detailed description of the transport of electron-positron pairs in the surroundings of the Geminga pulsar confirms the need for a reduced diffusion coefficient (Schroer et al. 2023).
This discovery raised a series of questions with potentially crucial implications for the standard model of the origin of Galactic CRs: 1) Is the diffusion coefficient small only around PWNe or everywhere in the Galactic disc? 2) A small diffusion coefficient in the TeV energy range implies the existence of a large level of turbulence, or a small coherence length of the turbulence (or both), but is this turbulence extrinsic or rather caused by CRs through the excitation of streaming instabilities? 3) What are the implications of the suppressed diffusion coefficient around sources for the grammage (the mass traversed by CRs during their propagation) and the positron excess?
These questions are not all independent from each other: for instance, the answer to the third question depends very sensibly upon whether the diffusion coefficient is only small around PWNe or also around SNRs, and it also depends on whether the diffusion coefficient is small everywhere in the disc or only in cocoons of limited size around sources.Moreover, the possibility for CRs to excite streaming instability and lead to reduced diffusivity depends on whether the particles are protons from a SNR or rather a quasi-neutral electron-positron beam from PWNe, or even a spatially chargeseparated beam of very high energy electrons or positrons from PWNe (Olmi & Bucciantini 2019;Olmi et al. 2024).
Here we focus on some of these questions: we investigate the information that can be gathered from γ-ray observations of the regions surrounding selected SNRs, in the attempt to infer whether there is evidence of suppressed diffusivity there too, as would be expected if the diffusion coefficient is small everywhere in the Galactic disc.Moreover, since the self-generation of turbulence induced by CR streaming is related to the luminosity of the source, studying SNRs with different intrinsic energy content should lead to different levels of turbulence in the surrounding medium, which should in turn reflect in different γ-ray emission.
In order to start addressing these issues, we consider two approaches to CR diffusion.The first is a simple phenomenological one already used by several authors before (for instance Aharonian et al. (2008); Nava & Gabici (2013); Bao et al. (2019); Bao & Chen (2021)), in which diffusion is assumed to occur in one or three spatial dimensions, at a level that can be calibrated to fit the data.In other words the normalisation and energy dependence of the diffusion coefficient are parameters of the problem.This approach is appropriate to describe transport in a disc where the coherence length is small or there is an enhanced level of turbulence, so that spatial diffusion is slow.The approach may also be used to check the possibility that some set of data may hint at small diffusion coefficient, independent of the origin of the suppressed diffusion.
The second approach is one in which the diffusion coefficient is an output of the calculation and the level of turbulence is due to the excitation of a CR streaming instability.Since the gradient in the number density of non-thermal particles is the largest near the sources, the propagation of CRs in these regions can be strongly affected by their self-generated turbulence, as previously pointed out by several authors (see e.g., Ptuskin et al. 2008;D'Angelo et al. 2016;Nava et al. 2016;D'Angelo et al. 2018;Nava et al. 2019;Recchia et al. 2022).
In the latter case, the assumption is that the pre-existing magnetic field is sufficiently ordered around the source (namely its coherence length is sufficiently large) so that transport can be, in first approximation, considered as one-dimensional.If the field around the source were disordered on scales smaller than the size of the source to start with, the problem would be 3D and the gradients would become too small to allow for self-generation, immediately outside the source region.Moreover, as discussed in previous literature, the process of wave generation induced by particle streaming becomes inefficient (even in one dimension) for particles with energies TeV.
It is clear that a reduced diffusion coefficient around sources (or in the disc) leads to potentially important effects in terms of grammage accumulated by CRs in the Galaxy (D'Angelo et al. 2016;Recchia et al. 2022).
In the present work we focus on several new aspects: 1) we first make an attempt at framing the discussion of the phenomenon of reduced diffusivity in the context of current open issues related to the transport of CRs in the Galaxy (for instance the positron excess and alternative models of CR transport in which grammage is accumulated near sources) and in the light of the discovery of the so-called TeV halos (see discussion in Sec. 2).2) We discuss the requirements to be imposed on clouds in the regions around SNRs in order to identify the effects of reduced diffusivity.This is an all but trivial process, as detailed in Sec. 3, and only few known SNRs remain viable after imposing these constraints.Here we focus on the two SNRs HB9 and W28 and the MCs around them. 3) We assess the role of molecular clouds (MCs) around SNRs in terms of the total grammage that CRs can accumulate around sources, a topic that has recently attracted much attention in the aftermath of alternative models of CR transport (Cowsik et al. 2014;Lipari 2017).We carry out these calculations including the effects of different types of damping.We show that the presence of clouds and the onset of self-confinement may lead to the accumulation of grammage near sources that significantly exceeds naive expectations.4) Finally we comment on different models of particle escape from SNRs in terms of the interpretation of the gamma ray emission from near-source regions.
As mentioned above, we apply the two approaches described above to two SNRs , HB9 and W28, for which there is evidence of extended γ-ray emission from MCs situated in the region around them and illuminated by CRs accelerated at the SNRs, a situation that has been suggested to be ideal for the investigation of CR escape from sources (see e.g., Aharonian et al. 2008).The first SNR we focus on, HB9, is known to be an under-luminous source and is chosen because one should expect weak CR induced instabilities to be excited around such an object and nevertheless γ-ray emission has been observed from the region around it.If we can confirm that no suppression of the diffusion coefficient is required to explain such γ-ray emission, then it would suggest that the small diffusivity is not a property of the Galactic disc at large.The second SNR, W28, has long been claimed to sit in a region with possibly reduced diffusion coefficient (see e.g., Gabici et al. 2007Gabici et al. , 2009;;Ohira et al. 2011;Li & Chen 2012), but as we discuss below, the claim is more or less serious depending on whether diffusion is treated in one or three spatial dimensions (see also Gabici et al. 2009).
The article is organised as follows: in §2, we discuss the difference between SNRs and PWNe as sources of high energy particles and we discuss the value of calculations such as the ones illustrated here to determine the impact of reduced diffusivity on CR transport.In §3, we list some observational facts about HB9 and W28 and discuss the reasons for selecting specific remnants and specific clouds around them to investigate the role of non-linear CR transport.In §4, we describe our theoretical approach to the phenomenological model and the one based on self-generation.The results in terms of gamma ray emission from specific clouds in the regions surrounding HB9 and W28 are discussed in §5.The important implications of suppressed diffusion in a region with MCs on the grammage of Galactic CRs are discussed in §6.Our conclusions are drawn in §7.

DIFFERENCES AND ANALOGIES BETWEEN PWNE AND SNRS AS PARTICLE ACCELERATORS
In order to understand the role of self-generation around sources and whether this phenomenon can be responsible for reduced diffusivity, it is important to identify analogies and differences between the two main classes of sources where these considerations are usually applied to.
Multi-frequency observations of SNRs and PWNe leave no doubt that both classes of sources act as effective particle accelerators.However the way particles are energized, the types of particles that get accelerated and the ways the particles escape from the sources are very different in the two cases, as we discuss below.
In SNRs, particle acceleration mainly occurs at the forward shock, where the ISM plasma is slowed down and heated up (see Blasi 2013Blasi , 2019, for reviews), for reviews).In an ideal situation, the accelerated particles are all advected toward downstream, both during the ejecta-dominated and the Sedov-Taylor (ST) phases, with the possible exception of the particles at the maximum energy.The time dependence of the maximum energy depends on many factors, such as the strength of the amplified magnetic field, as due to CR streaming instability excited upstream of the shock (Bell 2004;Amato et al. 2008).The integration over time of the escape flux eventually results in an approximately power-law, typically hard, spectrum (see e.g., Caprioli et al. 2010;Schure & Bell 2014;Cristofari et al. 2020).The rest of the accelerated particles are liberated into the ISM after the shock has slowed down and the flow has become subsonic.In fact, it is believed that escape should take place at the beginning of the radiative phase (see e.g., Cristofari et al. 2020), when instabilities are hard to sustain and particles may be free to leave.However, it is fair to say that a proof of this latter point is missing and there is no observational fact in its support, with the exception of the lack of non-thermal activity from objects that are in the radiative phase.By the time of escape, particles trapped in the downstream have been confined there for tens of thousands of years, and have suffered substantial adiabatic and, for leptons, radiative energy losses (Cristofari et al. 2021).
In a less ideal situation, for instance if the shock is not perfectly spherical and/or is broken at some locations, there may be leakage of the whole spectrum (namely at all energies) from the downstream region at any time.This ambiguity is part of the difficulty in modelling the escape of accelerated particles from a SNR (see discussion in Cristofari et al. 2021).
As discussed by Cristofari et al. (2020), in this basic picture, the high energy particles liberated by a SNR are mainly protons (and nuclei), while the electron spectrum is typically limited to energies 10 TeV, because of energy losses inside the acceleration region and during transport downstream of the shock.This is important in terms of formation of a TeV halo around SNRs, as we discuss below.
Regarding PWNe, particle acceleration is expected to take place mainly at the termination shock, where the pulsar wind is decelerated and heated.It is not yet known whether the process of particle acceleration is solely diffusive shock acceleration (DSA), for which the termination shock is ultrarelativistic and is usually not expected to behave as a good accelerator (Giacinti & Kirk 2018) , or magnetic reconnection (Lu et al. 2021) or a combination of the two.The spectrum of the accelerated particles, which, as inferred from multi-frequency observations, typically has the shape of a broken power-law with a change of slope at lorentz factor of the particles ∼ 10 5 , suggests that the process is not as simple as DSA alone.
It is important to realize that the observed TeV halos surround pulsars that, due to their birth kick motion, have already left their parent SNR.This is of crucial importance in that the accelerated particles may now leave the bow shock nebula that is formed due to the interaction of the pulsar wind with the ISM, through the comet-like tail that is typical of bow shock nebulae.The particles that are accelerated at the termination shock are expected to be electron-positron pairs, with an energy that can in principle be as high as the potential drop of the pulsar.This fact alone, the dominance of leptons and their escape from the acceleration region, leads to the conclusion that PWNe are more effective in injecting leptons in the surrounding region than SNRs, which in turn may explain why so far we have observed TeV halos around PWNe.
The observed TeV halos show a roughly spherical shape.The size and morphology of the γ-ray emission are best explained if the diffusion coefficient in the region around PWNe is small, typically 100−1000 times smaller than inferred in the whole CR propagation volume (see e.g., Abeysekara et al. 2017;Schroer et al. 2023).
Three options can be considered: 1) the diffusion coefficient is small everywhere through the disc of the Galaxy, perhaps due to the larger amount of extrinsic turbulence there; 2) the diffusion coefficient is small around sources, because of CR induced streaming instability in the regions where currents are stronger; 3) the diffusion coefficient is small in some locations, for some reason, and the PWN produces a TeV halo if it happens to cross one such regions.Case 3) would also apply to SNRs with the caveat that electrons are more efficient radiators, and as stressed above, PWNe are better sources of high energy leptons than SNRs.The occurrence rate of extended TeV emission around sources, as in case 3), would reflect the probability of finding such sources in the areas of low diffusivity, and at present cannot be estimated in any reliable way.
In case 1), since both PWNe and SNRs mainly reside in the disc, there should be evidence of γ-ray emission around both classes of sources, although PWNe would shine due to leptonic ICS emission while SNRs require target gas (for instance in the form of MCs) to produce γ-ray emission.
In case 2), the situation requires a deeper analysis: CRs escaping SNRs are known to be able to produce waves through the excitation of streaming instability (see e.g., Ptuskin et al. 2008), although the phenomenon is limited by different mechanisms of damping, depending on the environment (D'Angelo et al. 2016;Nava et al. 2016;D'Angelo et al. 2018;Nava et al. 2019;Recchia et al. 2022).The current of CRs escaping a young SNR may even be sufficient to excite a non resonant streaming instability, as discussed by Schroer et al. (2021Schroer et al. ( , 2022)).
On the other hand, the excitation of streaming instability by the pairs escaping a PWN did not receive equal attention: if the pairs are liberated into the ISM without charge separation, the total current that they carry vanishes and the non resonant branch (Bell 2004) of the instability cannot be excited.The resonant streaming instability is excited even in a situation of zero net current (Evoli et al. 2018) but the resulting level of self-generated turbulence is insufficient to explain the strong suppression of diffusion coefficient, unless the transverse cross section of the flux tube is assumed to be very small (Mukhopadhyay & Linden 2022), which is clearly at odds with the quasi-spherical morphology of the emission.Some evidence that electrons and positrons at energies comparable with the potential drop of the pulsar may be injected at different locations, thereby making the net current locally non-zero, was presented by (Olmi & Bucciantini 2019).The effect of the resulting current for the excitation of the non-resonant streaming instability is currently being investigated Olmi et al. (2024).This phenomenon might be related to another non-thermal unexpected occurrence, namely the long and narrow bright X-ray filaments seen around some PWNe (see for instance Hui & Becker (2007)).This latter phenomenon might be somehow related to the production of TeV halos at later stages.
In the following, we will develop diffusion models in cases 1) and 2) and apply them to SNRs as exemplified by HB9 and W28.In these cases, there is evidence of extended γ-ray emission which, when projected, seems to lie outside (but near) the remnants, and the γ-ray emission appears to be coincident with large MCs, not in direct contact with the remnants' shocks.

SOME OBSERVATIONAL FACTS ABOUT HB9 AND W28
Here we provide some observational facts about SNRs HB9 and W28.Gamma-ray emission appears to arise from several MCs around these SNRs, but only some of this emission can plausibly be used to infer information about the transport of escaping particles in the near-source regions.
There are two MCs in the region around HB9, one (called "R1") that appears to be in touch with the SNR shock and the other (called "R2") that is located at an angular distance from the explosion center of ∼ 1.4 • (corresponding to ∼ 15 pc) (Oka & Ishizaki 2022).
Around W28, there are four MCs, namely, MCs N (which is interacting with the shock directly) A, B, and C (see e.g., Aharonian et al. 2008;Cui et al. 2018).The γ-ray emission of HESS J1800−24 C (corresponding to MC C) is suspected to be contaminated by the contribution from the SNR candidate G5.71−0.08 (Brogan et al. 2006;Aharonian et al. 2008).Hence here we exclude from our analysis the contributions from both MCs N and C.
The γ-ray emission from the surroundings of both SNRs share some similarities in terms of spectrum.A recent Fermi-LAT data analysis of HB9 (Oka & Ishizaki 2022) found that the γ-ray spectra from MCs R1 and R2 are somewhat harder (with power-law indices Γ = 1.84 ± 0.18 and Γ = 1.84 ± 0.14) than that at the SNR shell (Γ = 2.55 ± 0.10).Likewise, for W28, MCs A and B show much harder γ-ray spectra than MC N (which is interacting with the shock directly).The proton spectral indices required to fit the γ-ray spectra of MCs N, A, and B are 2.7, 2.2, and 2.4, respectively (see Fig. 6 in Cui et al. 2018).Both instances can be considered as indications of a possible energy-dependent diffusion of particles to the locations of the MCs.

HB9
HB9 is a thermal composite (or mixed-morphology) SNR, with an estimated age t age ≈ 7 kyr (Leahy & Aschenbach 1995;Iqbal et al. 2009).A radio continuum shell at 408 MHz is found to have a good spatial correspondence with HI and CO clouds (Sezer et al. 2019).The distance to the SNR is measured to be d ≈ 600 pc (Sezer et al. 2019), which is roughly consistent with the estimate (800 pc) in Leahy & Tian (2007).The SNR has a large angular radius of 1 • (∼ 10 pc at 600 pc), with its shock suggested to be interacting directly with the nearby MC R1 as flagged by the detection of 1720 MHz OH maser (Frail et al. 1996).The density of the X-ray-emitting gas is estimated to be ∼ 1.4f −1/2 fill (d/0.6 kpc) −1/2 cm −3 (where f fill is the filling factor) in the eastern region and ∼ 0.9f −1/2 fill (d/0.6 kpc) −1/2 cm −3 in the western region (Sezer et al. 2019).Based on cloud evaporation models (White & Long 1991), Leahy & Tian (2007) suggested that the preshock intercloud density is ∼ (3-10) × 10 −3 cm −3 , with the ratio between the cloud mass m MC and the intercloud medium mass being ∼ 20-50.Therefore the mean ISM density is ∼ 0.2-0.5 cm −3 , which is consistent with the estimate in Sezer et al. (2019).In order to fit the radius of the SNR ∼ 10.5 pc, we adopt the supernova ejecta mass M ej = 3 M ⊙ , the explosion energy E SN = 3 × 10 50 erg, and the mean ambient ISM density n ISM = 0.3 cm −3 .It follows that the velocity of the shock was 2.1 × 10 3 km s −1 at the beginning time of the ST phase, t ST = 1.5 kyr, and is 6.4 × 10 2 km s −1 at the present time (Truelove & McKee 1999) 1 .
Assuming simple dynamical evolution of the SNR in a homogeneous medium, the maximum energy of the accelerated particles in HB9 can be estimated following Equation ( 21) in Bell et al. (2013) or Equation ( 27) in Cardillo et al. (2015).The maximum energy was ∼ 20 TeV at the beginning of the ST phase and is ∼ 5 TeV at present.These estimates can be appreciably lower if damping happens to play an important role.This may be the case especially at late times, when the non resonant hybrid instability can no longer be excited, and ion-neutral damping may be in action, making the maximum energy of the accelerated particles appreciably lower than estimated above (see for instance Fig. 2 of Ptuskin & Zirakashvili (2003)).
Since MC R1 is observed to have been reached by the SNR shock, the γ-ray emission from this location cannot provide information on the CR diffusion in the near-source region.Hence below we focus on MC R2, from which however at present there is no evidence of TeV γ-ray emission.It is not clear whether this may be due to a too low maximum energy or to the fact that the > 10 TeV protons escape too quickly from the circumsource region.

W28
SNR W28, a prototypical thermal composite SNR, is one of the best-known middle-age remnants, located in a star-forming region which hosts HII regions M8 (at a distance ∼ 2 kpc, Tothill et al. 2002) and M20 (at ∼ 1.7 kpc, Lynds & Oneil 1985).NANTEN 12 CO (J=1-0) data reveal that MCs HESS J1800−240 A, B, and C to the south of the remnant appear coincident with the γ-ray emission from the TeV source HESS J1900−240 (Aharonian et al. 2008).MC C is considered to be correlated with another SNR, G5.7−0.1 (see e.g., Joubert et al. 2016), and hence we do not consider it here.
We assume that they are at the same distance from the shock.The distance of W28 is suggested to be ∼ 1.9 kpc based on HI observations (Velázquez et al. 2002), and its angular radius is ∼ 0.35 • (∼ 12 pc).Radio and near-infrared observations indicate that the ISM around W28 is made up of mainly low-density inter-clump medium which takes up to 90% of the volume and dense clumps (Reach et al. 2005).The mean density of the ISM is ∼ 5 cm −3 .The age of W28 is estimated to be ∼ 40 kyr based on the cooling of the associated neutron star (Yakovlev & Pethick 2004).We assume E SN = 10 51 erg and M ej = 6M ⊙ .
The SNR entered pressure-driven-snowplow (PDS) phase (in which all CRs are assumed to be released into the ISM) at time t PDS ∼ 5 kyr after the explosion, corresponding to a radius of ∼ 7 pc (Cioffi et al. 1988;Leahy & Williams 2017).MCs A and B are located ∼ 0.45 • from the SNR center, corresponding to a distance ∼ 16 pc.The masses of the two clouds are 6 × 104 M ⊙ and 4 × 10 4 M ⊙ , respectively (Aharonian et al. 2008;Gabici et al. 2010).
The transport of non-thermal particles in the near-source region is dominated by diffusion and advection, two processes that are properly described by the Parker equation: where f (p, x, t) is the isotropic part of the CR distribution function in phase space, u( x, t) the advection velocity, and D(p, x, t) the diffusion coefficient.The function Q(p, x, t) describes the injection, in terms of particles leaving a SNR.In the following, we will illustrate our results in terms of two models for diffusion, one of phenomenological nature, in which D is parametrized, and the other in which D is calculated in a non-linear way, as a result of self-generation of waves.
In both cases, the generation of γ-ray radiation is assumed to be due to the inelastic interactions of protons with the gas in a MC located at given distance from the SNR shock.In general, the flux of γ-rays of hadronic origin at photon energy ε can be calculated as (Koldobskiy et al. 2021): where n MC is the number density of the gas in the MC, dσ γ (ε, p)/dε is the differential cross section of production of a photon at energy ε by a proton with momentum p (Koldobskiy et al. 2021).

Phenomenological diffusion coefficient
The phenomenological approach to the description of CR transport around the SNRs HB9 and W28 is based upon the assumption that the diffusion coefficient around these sources is spatially constant and numerically estimated to be a fraction of the Galactic diffusion coefficient D Gal , as estimated from the secondary/primary ratios (see e.g., Ptuskin et al. 2009).The exact numerical value of D Gal is not very important in this context, since it is only used as a reference with respect to which the reduced diffusivity manifests itself.Thus, we assume D(p) = D Gal /κ, with the suppression factor κ > 1, p the proton momentum, and for the Galactic diffusion coefficient we use D Gal = 3.6×10 28 [pc/(1 GeV/c)] 1/3 cm 2 s −1 , aware that deviations from this trend are likely, especially at low energies (Evoli et al. 2020).This simple approach is adopted here following the same logic as in the past literature on the topic (Aharonian et al. 2008;Nava & Gabici 2013;Bao et al. 2019;Bao & Chen 2021), namely to have a quick and approximate assessment of the possible suppression of diffusivity in the near-source regions.
The diffusion equation around the source in this phenomenological approach reads: where the ∇ 2 operator is written in 1D (i.e., the particles only diffuse along the z axis) or 3D (i.e., the particles diffuse spherically symmetrically in the radial coordinate r) depending on the situation at hand and the spatial coordinate is x = z in 1D and x = r in 3D.The injection term Q will be written as Q 1D and Q 3D in the two cases that we investigated , respectively.As discussed in §3, the age of W28 is such that at present it is expected to be in the PDS phase.Hence we assume that the release of the accelerated particles occurred impulsively at the beginning of the PDS phase (t PDS ∼ 5 kyr for W28).
For HB9, since it is still in its ST phase, we assume that the injection of particles occurs continuously in time but at a constant rate.In other words, we neglect the change in time of the velocity and radius of the forward shock.We do so to avoid making additional assumptions on the environment in which the explosion took place.
The geometrical center of the problem, x = 0, is taken as the location of the injection of particles, so that Q(p, x, t) can be written as in the 1D case and in the 3D case, where t esc is the time the particle escapes, ρ the radius of the cylinder coordinate, Ṅesc (p, t esc ) the injection rate of protons, H the Heaviside step function, and ρ Inj the radial spread of the injection region (transverse size of the flux tube in the one dimensional case).We write the injection rate as where p cut is the cutoff momentum, and Q SNR is the injection constant of the SNR which can be obtained using the normalization condition: Here ξ is the efficiency of conversion of the SNR kinetic energy to CR energy (typically ξ ∼ 0.1).In the 1D case, the choice of the SNR radius at the time of particle escape (ρ Inj ) is influenced by several uncertainties: it depends upon the the detailed gas distribution in the region around the SNR and the type of supernova event we deal with.Technically, the uncertainties are consolidated into chosen different values of ρ Inj , and the implications of the different choices will be discussed.
As discussed above, it is possible that escape of accelerated particles occurs only from the upstream region of the shock.In this case the spectrum of the escaping particles may be limited to a narrow neighborhood of the maximum energy at that given time, E max (t), while the spectrum in the nearsource region from the overlap of previous times is expected to show a low energy cut at E max (t).Occasionally in the calculations that follow we will adopt this recipe to mimic the scenario in which escape occurs in this way.
For HB9, which is currently in the ST phase, we calculate the mass of MCs based on an average tube radius of R tube = 8 pc.This average is determined from its present radius of 10 pc and its initial radius of 5 pc at the start of the ST phase.We adopt these values as benchmark while we will warn the reader when other possibilities will be considered.
In contrast, W28, due to its high mean ISM density of approximately 5 cm −3 (Reach et al. 2005), enters its PDS phase early at about 5 kyr and has a shock radius of roughly 7 pc.Given its current age of ∼ 40 kyr, we assume that all particles escaped at the same time t PDS when the radius of the SNR was ρ Inj = 7 pc.
The maximum energy of particles accelerated at the shocks in the two selected SNRs depends on numerous ingredients that are poorly known (such as possibility to excite the non-resonant instability at some time during the evolution of the remnant, role of damping and conditions in the circum-source region).We check that it is likely that there were the conditions for the excitation of the non-resonant instability at the beginning of the ST phase, though for some choices of the parameters the condition is satisfied only marginally.For the sake of the present calculations, given the large uncertainties in the maximum energy and its dependence on time, we assume here that the spectrum of the CRs has a fixed exponential cutoff at 20 TeV for HB9 and 100 TeV for W28.
In the 1D situation, the diffusion problem has a cylindrical symmetry and if the diffusion coefficient is assumed to be spatially constant, the solution can be written as where R d = 4D ISM (t age − t esc ) is the distance that particles can diffuse away from the source at the given time t age .

Self-generated diffusion
In the region surrounding a source we can naturally expect CR energy density in excess of the mean Galactic value, as well as stronger CR gradients that result from the presence of a source.In these conditions, as previously discussed by (Ptuskin et al. 2008;Malkov et al. 2013;D'Angelo et al. 2016;Nava et al. 2016;D'Angelo et al. 2018;Nava et al. 2019;Recchia et al. 2022), the transport of CRs may be dominated by self-generated perturbations, mainly in the form of Alfvén waves moving parallel to the regular Galactic magnetic field in which the source is assumed to be located.Since the gradient has a well defined sign, most Alfvén waves move away from the source.This results in an advective term in the transport equation with an advection velocity that equals the Alfvén speed.The main novelty of the present calculations is the fact that we account for the possible presence of dense MCs located in the regions around specific sources and with properties that may lead to the identification of reduced diffusivity.Moreover we investigate the effect of the MCs in terms of grammage accumulated by CRs while propagating away from the sources.
It is important to realize that the process of self-generation is intrinsically non-linear: the level of perturbations that is induced by this phenomenon is related to the local density of CRs at a given location, but on the other hand the CR density is larger when the diffusion coefficient is smaller, namely when the self-generation is more effective.The rate of particle injection per unit volume, and hence the spatial extent of the source region are also important parameters of the problem.
In the below we will assume that the transport of CRs in a region of size similar to the coherence length L c around the source can be considered approximately one dimensional (Ptuskin et al. 2008).Following Beck et al. (2016), we adopt a value L c ≈ 220 pc for the Galactic magnetic field B 0 .When the particles reach a distance from the source that is appreciably larger than L c , the transport is expected to change over to a 3D phenomenon (Giacinti et al. 2013).When this happens the CR density and its gradient become smaller and the role of self-generation quickly becomes unimportant in terms of confinement on regions of size L c .To be more precise, we should say that the local gradients become unimportant, while it is possible that the gradient on Galactic scales may drive the turbulence responsible for the transport of the bulk of CRs on Galactic scales (Blasi et al. 2012), at least up to energies TeV.
As mentioned above, the calculations are carried out in the so-called flux tube approximation, namely we assume that the CRs remain trapped in a region whose cross section is the transverse size of the source that injects them into the ISM.Recently, hybrid particle-in-cell simulations (Schroer et al. 2021(Schroer et al. , 2022) ) were used to show that the reduced diffusion inside the flux tube results in an enhanced pressure gradient that leads the tube to inflate.We neglect this aspect here since it is expected to be the most important for young SNRs.
The 1D diffusion-advection equation reads Since in weakly perturbed fields u is +v A (v A := B 0 / √ 4πm i n i ) for z > 0 and −v A for z < 0, du/dz = 2v A δ(z).The diffusion coefficient D(p, z, t) can be calculated as where r L is the Larmor radius of the particles and the dimensionless spectral power per unit ln(k), F (k, z, t), is calculated at the resonant wave number k = 1/r L (p).The evolution of F (k, z, t) is described by the equation: where the growth rate of self-generated turbulence is given by (Skilling 1971) and Γ D = Γ NLL + Γ FG + Γ IND is the damping term, where Γ FG = 2kv A / √ kL c is the damping rate from interactions with pre-existing MHD turbulence (Farmer & Goldreich 2004), is the rate of non-linear Landau damping (see e.g., Wiener et al. 2013;Nava et al. 2019, and references therein), Γ IND is the ion-neutral damping rate, and k B the Boltzmann constant.The Γ IND can be obtained by solving the equation (Zweibel & Shull 1982) where ω k = kB 0 / √ 4πn i m i , ω is the complex frequency of the wave, ν the ion-neutral momentum transfer ratio and ǫ = n i /n n is the ratio of ion density to neutral density.With ω = ω R + iΓ IND /2, Equation 14 can be written as The equation s are solved using boundary conditions that F | z=Lc/2 = F Gal and a free escape boundary condition f | z=Lc/2 = 0 for the particles, where F Gal = r L c/[3D Gal (E)] is the Galactic turbulence spectral power.The injection term Q(p, z, t) is the one described in Equations 4.

RESULTS
Here we illustrate the results of our calculations of CR transport around sources using phenomenological models and self-generation models and their implications in terms of γ-ray emission from the MCs located in the region s around two SNRs, HB9 and W28.

Phenomenological models of CR transport
The phenomenological models allow us to infer some generic properties of CR transport in the nearsource regions without delving with the complex aspects involved in self-generation.In particular, these models can provide an easy flag for appreciable suppression of diffusivity in the near-source regions (see for instance Gabici et al. (2007); Nava & Gabici (2013)).
The γ-ray emission expected from the MC (R2) near the HB9 SNR is shown in Figure 1, for the 1D (left panel) and 3D (right panel) case, respectively.The different curves refer to different distances and masses of the R2 MC, so as to illustrate the effect of the uncertainty in these two quantities.One can see that for the ranges of mass and distance that are allowed by observations, a suppression of the diffusion coefficient with respect to the Galactic one is required, for both the 1D and 3D cases.The suppression factor κ ranges between ∼ 10 and ∼ 200.
It is important to stress that the main need for a suppression in the diffusion coefficient comes from the spectral shape of the γ-ray emission in the low energy range: the downward inflection of the spectrum requires that the low energy particles do not reach the cloud.Clearly this effect is more pronounced, namely requires more suppression of the diffusion coefficient, when the distance from the SNR shock to the cloud d cl is larger.
These results would suggest that for the HB9 SNR there is evidence of suppressed diffusivity and that such evidence is independent of the assumed dimensionality of the problem.Given the importance of the statement, it is important to check whether other factors may mimic the same phenomenology.One possibility is related to the way particles escape the remnant: as discussed above, the general, and somewhat idealized, picture of particle escape is that CRs leave the remnant from upstream of the shock in such a way that only particles with energy in a narrow neighborhood of the maximum energy at that given time can escape.This means that at a given time only particles with energy > E max (t) can fill the medium surrounding the SNR.From the phenomenological point of view we can mimic this picture by repeating the previous calculation but assuming that the spectrum of escaping particles has a minimum energy and check how low such minimum energy should be in order to fit the observed γ-ray emission.
The results of this exercise are summarized in Fig. 2, where we allowed the minimum energy to be 1 TeV and 100 GeV and the diffusion coefficient is the same as assumed for the Galaxy at large.The suppression of the γ-ray emission at low energies here is not due to confinement but rather to the fact that only particles with energies above some minimum value have escaped the remnant at the time corresponding to the age of HB9.The low energy behaviour of the γ-ray emission reflects the energy dependence of the differential cross section of pion production ∝ E −1 π .This exercise, based on a phenomenological description of the diffusion of particles in the nearsource region, suggests that there is a degeneracy between the model of particle escape and the phenomenon of reduced diffusivity: the observed γ-ray emission is both compatible with a maximum energy in the hundred GeV range at the age of HB9 or with the escape of all particles from the shock and a suppressed diffusivity around the remnant, although the physical reason for the low energy suppression in the γ-ray emission, as discussed above, is completely different in the two scenarios.
The issue of whether it is reasonable to have a maximum energy in the ∼ 100 GeV range for a SNR having the age of HB9 is strongly related to the damping mechanisms that are at work in the acceleration region.For instance, in Fig. 2 of Ptuskin & Zirakashvili (2003) one can read the maximum energy as a function of time in the presence of damping.If only non-linear Landau damping is at work, it is not feasible to decrease the maximum energy to such low values.If ionneutral damping near the shock dominates, a maximum energy below ∼ 1 TeV seems somewhat too low for HB9, but it is hard to rule it out, given the uncertainties on the local density and the ionization fraction in the acceleration region.Hence, unfortunately some level of degeneracy between the two scenarios discussed above cannot be excluded.
The γ-ray emission expected in our phenomenological models of CR transport around the W28 SNR is shown in the left (right) panel of Figure 3 for MC A (B).Our predictions are compared with the γ-ray observations of Fermi-LAT and H.E.S.S.We use the nominal values of the mass of the two MCs.
In the 1D models, one can see that the observations are well described by adopting the Galactic diffusion coefficient with no need for suppression.Note that the acceleration efficiency ξ required in these calculations is extremely low, of order 0.3%.We also find that in the 1D situation, d cl can hardly affect the spectra because for z ≪ R d ∼ 4.0 × 10 2 (E/1 TeV) 1/6 pc, the distribution function f (p, z, t) ∝ 1/R d , but is independent on d cl .This might suggest that not the whole volume of the clouds is actually illuminated by the CRs accelerated in W28, so that more reasonable values of the efficiency can be adopted.In any case, larger and more common values of the acceleration efficiency would strengthen the evidence that no suppression in the diffusion coefficient is required in the case of 1D transport around W28, consistent with the earlier results by Nava & Gabici (2013).
In the 3D case of CR transport, as expected, the requirement for suppressed diffusivity becomes evident: for values of the distance d cl between the cloud and the remnant that are accepted in the literature and for the nominal values of the masses of the two clouds, the diffusion coefficient is required to be suppressed by a factor ∼ 50 − 200.This is a generic feature due to the fact that in the 3D case the CR density drops faster with the distance and hence a smaller diffusion coefficient is required to obtain the same γ-ray emission.
As discussed earlier in this manuscript, this latter scenario is of limited interest in that, in order for the diffusion to be treated as three dimensional, the coherence scale of the turbulent field is required to be much smaller than the size of the SNR, which by itself suggests that the diffusion coefficient is smaller than the Galactic value.Notice also that, from a purely observational point of view, the data on the low energy γ-ray emission from W28 are not unambiguously suggestive of a flux reduction, namely there is no clear need for a suppression of diffusivity.

Models of self-generated transport
The excitation of streaming instability near the source leads to a reduction of the diffusion coefficient that in turn increases the CR density (CR confinement).The solution of the transport equation is made difficult by the fact that the diffusion coefficient entering the diffusion equation depends on the distribution function of particles, making the problem non-linear.Moreover, the spectrum of self-generated perturbations evolves in time trying to reach a balance between the excitation of the instability (growth) and damping.The latter depends rather critically on the level of ionization of the background plasma.
Below we illustrate the effects of this phenomenon in terms of the CR confinement and γ-ray emission in the regions around HB9 and W28.We focus on a 1D problem in which CRs propagate non-linearly along the local flux tube determined by the direction of the Galactic magnetic field in the region where the source is located.The 3D case is less interesting in terms of self-generation, since the small coherence scale of the field necessary to justify the assumption of 3D transport already implies that the local diffusion coefficient is smaller than in the Galaxy at large.
For the sake of simplicity and to limit the number of parameters, we assume that the transverse sizes of the SNR and of the cloud are the same (ρ Inj ) and that the cloud is a cylinder extending from a distance z 1 to a distance z 2 away from the source (see Table 2 for the numerical values of these parameters as well as the cloud mass).
The spectrum of γ-ray emission from MC R2 around HB9 in the self-generated case is shown in Figure 4 for the case of fully ionized and partially ionized plasma in the near-source region.The left (right) panel refers to a coherence length L c = 100 pc (300 pc).The injection spectrum adopted for the computation is ∝ E −2.4 all the way down to small energies.As discussed in Sec.5.1, it is possible that this is not the case and that escape from upstream of the shock is concentrated in a narrow energy region around the maximum energy at that given time.This picture results in an 10 0 10 1 10 2 10 3 10 4 Energy [GeV] . γ-ray emission from MC R2 near HB9 in the case of a spectrum of particles with a low energy cut at E < 100 GeV.The three curves show the case of Galactic diffusion coefficient, the case of selfgeneration with R tube = 10 pc, (n i , n n ) = (0.3, 0) cm −3 and the case of self-generation with R tube = 10 pc, (n i , n n ) = (0.05, 0.25) cm −3 .The power-law index in the momentum space α = 4.4 is used.
effective low energy cutoff in the spectrum injected over time, at the maximum energy reached at the age of HB9.We will comment below on how this scenario affects self-generation.
In the partially ionized case, ion-neutral damping leads to substantial suppression of Alfvén waves produced at high wavenumbers, responsible for resonances with low energy CRs.As a consequence, the diffusion coefficient of the CRs with energies 1 TeV becomes large and inhibits self-confinement.In the case in which the gas is fully ionized, the main channel of damping is non-linear Landau damping, but its effect is not sufficient to prevent self-confinement.This can be seen very clearly in Figure 4, which shows a downturn of the γ-ray emission at energies 100 GeV, corresponding to the fact that the diffusion coefficient is lower and the low energy particles can hardly reach the R2 MC.This interpretation can be easily validated by looking at the γ-ray emission obtained for a Galactic diffusion coefficient (red curve s) and a Galactic diffusion coefficient suppressed by hand by a factor 200 (green curve s).The latter shows a similar downturn in the γ-ray emission at energies 100 GeV.The main difference between the cases with L c = 100 pc and L c = 300 pc is that in the latter the confinement extends to higher energies, as expected.
Taken at face value, this finding would be clear evidence that even in the case of an under-luminous SNR such as HB9 there is suppression of diffusivity around the source.Unfortunately, as discussed in Sec.5.1, the low energy γ-ray flux reduction suggesting a smaller diffusivity can also be caused by the differential escape of CRs from the upstream region of the shock, if the maximum energy of the accelerated particles for a SNR with the age of HB9 is as low as ∼ 100 GeV.Although this condition requires a role of ion-neutral damping at the shock that is probably larger than what one would expect (Ptuskin & Zirakashvili 2003), the scenario cannot be ruled out given the uncertainties in the environmental parameters.
Figure 6.γ-ray emission from MC A (left panels) and MC B (right panels), modelled in the self-generated turbulence case.For all curves we have considered the convection.The injection term is impulsive in time at the beginning of ST phase.In the first row we assume the free escape boundary at 100 pc, in the second row we assume the free escape boundary at 300 pc.In all rows we assume that the total proton energy is 4 × 10 48 erg.
In order to show more clearly the effect of this degeneracy, in Fig. 5 we show three cases of CR transport in a situation in which the spectrum of escaping particles is characterized by a low energy cut at 100 GeV: the first case corresponds to assuming that the diffusion coefficient is the same as in the Galaxy at large; in the second case the self-generation is active but there is no neutral gas; the third case shows the result of self-generation with ion-neutral damping taken into account.The plot clearly shows that the observed γ-ray emission simply requires a small renormalization of the cloud mass in the three cases, within the uncertainties on such parameter.Hence the evidence for self-generation around HB9 can, at present, only be considered as circumstantial.
In Figure 6, we show the γ-ray fluxes from MC A (left panels) and MC B (right panels) around W28 for different values of the n n /n i ratio, where the total density is fixed to 5cm −3 .One can find that the effect of self-generation is only visible in the unlikely case of high density of fully ionized gas (red curve), and the resulting γ-ray emission is not well described by this situation.In all cases with neutral gas the ion-neutral damping is so severe that the diffusion coefficient at all proton energies of relevance here remains close to the Galactic diffusion coefficient.In fact, the low energy suppression in all the curves that refer to the cases with neutral gas occurs basically at the same energies obtained by using the Galactic diffusion coefficient (see the top panels in Figure 3).At high energies, small differences with respect to the cases shown in Figure 3 are due to the different boundary condition s: in the self-generation cases we assume free escape at some distance from the SNR, while in the phenomenological models the boundary escape is at infinity.

GRAMMAGE
Although at present there is no compelling evidence of self-generation induced by CRs around SNRs, this phenomenon must occur at some level.The difficulty in observing, as discussed above, may be due to the complexity that is typical of the regions where these accelerators are located: such complexity would lead us to a multiplicity of phenomena taking place which would lead those regions to be excluded by an analysis such as the one presented here.
The presence of regions of reduced diffusivity near sources leads to an increased residence time in these regions, and to the corresponding accumulation of grammage that is to be added to the grammage accumulated by CRs during their transport in the Galaxy.
This phenomenon is particularly important in the presence of dense MCs in the circum-source regions, despite the fact that the time accumulated in such regions remains small compared to the escape time in the Galaxy (D'Angelo et al. 2016;Recchia et al. 2022).This is due to the larger mean density of the target gas that CRs traverse when MCs are present.The grammage accumulated in the absence of MCs was computed by D' Angelo et al. (2016) and was found to be sizeable only for relatively large density of ionized gas and absence of neutral hydrogen, unlikely to be the case in the disc of the Galaxy where SNRs are located.On the other hand, the confinement time is enhanced even if the ionized gas density is small.In this latter case, as discussed below, the accumulated grammage can be appreciable if MCs are present, even though surrounded by relatively rarefied ionized gas.
Here, for simplicity, for the sole purpose of illustrating the effect, we assume that the CR protons are injected impulsively at z = 0 at the edge of a cylinder with a given radius, mimicing the transverse size of the source, depending on the gas density in the surrounding medium.
We solve the 1D diffusion-advection equation in the presence of self-generation, assuming a free escape boundary condition at a distance of 100 pc from the SNR, and we compute the grammage as where ρ gas is the mass density.The grammage accumulated by CRs in the near-SNR region is shown in Figure 7 for n i = 0.01cm −3 (left panel) and for n i = n n = 1cm −3 (right panel).In both cases the cloud is assumed to be such that the MC density times the MC length n M C L M C = 800cm −3 pc (clearly in a one dimensional problem this is the only relevant quantity).The size of the source is set to 30 pc for the low density case (left panel) and 12 pc for the high density case (right panel).
The different curves refer to different values of d cl , the distance between the source and the center of the cloud.As expected, the latter parameter has little effect on the resulting grammage.
In the left panel of Figure 7, the grammage accumulated in the surrounding gas (without clouds) and using the Galactic diffusion coefficient is shown with a purple line.Even at ∼ 10 GeV the corresponding grammage is negligibly small, ∼ 10 −4 g cm −2 , as expected.When self-generation is taken into account, the residence time in the circum-source region increases and the corresponding grammage (black curve) increases correspondingly, but remains negligible.When a MC is present, on the other hand, the grammage becomes much larger, by roughly an amount where L c is the location of the free escape boundary.It is worth stressing that even in the absence of self-generation, in the presence of a MC the grammage can be boosted considerably, and a few g cm −2 can be accumulated at 10 GeV.Clearly, this does not mean that CRs from all SNRs undergo such phenomenon, but it is worth keeping in mind that the presence of MCs may affect the estimate of the grammage accumulated near sources.The situation when neutral gas is present (right panel in Figure 7) is more interesting in that the effect of ion-neutral damping limits self-generation in an interesting manner as described below.
At very low energies, the effect of ion-neutral damping is so severe that the diffusion coefficient increases roughly to the Galactic value, hence in the absence of clouds the grammage overlaps with the purple line.When the cloud is present, the enhanced grammage is solely due to the fact that now the effective density in the region is ∼ n M C L M C /L c ∼ 8cm −3 , so that the grammage increases correspondingly.At very high energies the diffusion coefficient is again comparable with the one obtained with Galactic diffusion, due to the fact that the particles of such high energies escape too quickly for the effect of self-generation to be relevant.At energies between ∼ 100 and ∼ 10 4 GeV, on the other hand, ion-neutral damping is not sufficient to deplete self-generated waves completely, and hence the grammage shows a substantial enhancement, to ∼ 1 − 10 g cm −2 .Clearly, such an enhanced grammage would not be compatible with the observed grammage, as inferred from the B/C ratio and similar indicators, if this were a general phenomenon, around any source of CRs in the Galaxy.
On the other hand, measurements of CR fluxes have reached such a high level of precision that it is conceivable that effects like that described above, even if subdominant, may be observable as deviations or anomalous behaviours of standard indicators of CR transport (for instance the B/C, B/O ratios).
We mention in passing that there are alternative models of CR transport, mainly proposed to accommodate the positron excess in a purely secondary context, that advocate that most grammage might be accumulated in dense cocoons near sources (Cowsik et al. 2014;Lipari 2017).Although such models face daunting problems of self-consistency, especially when unstable secondary nuclei are included in the analysis, they forced us to consider the problem of CR transport near sources with more attention: given the accuracy of the data available at this time, even relatively small deviations from the standard picture of CR transport, such as those induced by CR transport around SNRs in the presence of MCs, might soon be unveiled.

CONCLUSION AND DISCUSSION
There are several reasons for exploring the possibility that CRs excite instabilities around sources, that in turn produce waves responsible for self-confinement.The first such reason is that there is now circumstantial evidence of reduced diffusion coefficients around some PWNe and (less clearly) SNRs, inferred from detailed γ-ray observations.The second reason is that this phenomenon, if sufficiently general, can impact on the grammage that CRs traverse while propagating in the Galaxy, a quantity that we infer from observing secondary to primary ratios in the cosmic radiation.For reasons explained above, SNRs and PWNe are expected to behave in qualitatively different ways and a dedicated investigation should be devoted to each of the two classes of sources.
While this issue was investigated in previous articles (D 'Angelo et al. 2016'Angelo et al. , 2018;;Recchia et al. 2022), here we focused on several new aspects of the problem: 1) we identified a few instances of SNR-MC associations for which the γ-ray emission may flag the suppression of diffusivity and discussed a set of criteria to be adopted for this sake; 2) we analyzed the γ-ray emission from two SNRs, W28 and HB9, that are not expected to be particularly bright; 3) we investigated the effect of self-confinement on the grammage that CRs traverse while escaping the circum-source regions with special emphasis on the presence of MCs.
Extracting relevant physical information from the γ-ray emission associated with MCs around SNRs is particularly complex because of environmental issues: first, there are several cases in which the cloud is hit by the shock front of the SNR, a situation sometimes flagged by maser emission; in this case the acceleration process is severely affected and one should not necessarily expect that the γ-ray emission provides information about the diffusion coefficient in the medium surrounding the SNR.We avoid such situations.
Second, in some cases, due to projection effects, it may be unclear whether a MC is associated with a given SNR.We try to avoid such situations as well, in that the results would be, to say the least, ambiguous.
In the attempt to identify some clean cases, in which the emission from (some of the) clouds is reasonably well associated with the production of pions by CRs diffusing away from a given SNR, we chose one MC around SNR HB9 and two clouds around W28, as targets of our investigation and carried out the calculations in such cases.
A rough assessment of the possible need for a suppression in the diffusion coefficient around sources can be made by using a phenomenological calculation, in which the diffusion coefficient is parametrized as a fraction of a Galactic-like diffusion coefficient (Aharonian et al. 2008;Nava & Gabici 2013;Bao et al. 2019;Bao & Chen 2021).Meanwhile, a non-linear computation that includes self-generation and damping, clarifies whether the conditions are appropriate to expect a suppression of the diffusion coefficient as a result of the excitation of resonant streaming instability.
We find that the case of MC R2 around the SNR HB9 is particularly useful: the diffusion coefficient required to account for the low energy part of the γ-ray emission is about two orders of magnitude smaller than the Galactic one and we show that it can be interpreted as a result of self-generation.Unfortunately this apparently clear signature is made less evident by the degeneracy between this scenario and one based on a different picture of particle escape from the SNR: if the maximum energy for HB9 at the present age is as low as ∼ 100 GeV, and the escape is assumed to be quasimonochromatic, namely the particles escape the remnant from upstream only when their energy is within a narrow neighborhood of the maximum energy at that given time, then the effective spectrum of the particles in the near-source region has a low energy cut which reflects into a bending of the γ-ray spectrum, very similar to the observed one.We show that the γ-ray spectrum is equally well fit in the two scenarios, which leads to the impossibility to discriminate between them.We conclude that the evidence of suppression of diffusivity around HB9 is, at present, only circumstantial.
The case of W28 is even more ambiguous: phenomenological models suggest that the evidence of small diffusion coefficient is limited to a 3D treatment of transport around this source, confirming previous results by Nava & Gabici (2013).This is due to the fact that in three dimensions the density of CRs drops faster and for a given γ-ray emission this implies a smaller diffusion coefficient.On the other hand, a 3D treatment is justified only if the field lines in the region around the source are tangled on scales smaller than the source size.In this situation one can expect a smaller diffusion coefficient even in the absence of self-generation.
In the one-dimensional case, our calculations of non-linear transport around W28 confirm that one should not expect appreciable self-generation, due to the presence of neutral gas in the region around W28.
In conclusion, neither one of the SNR-MC associations considered here provides convincing proof that the diffusion coefficient is appreciably suppressed around these sources.One could argue that, besides being cleaner than other cases, neither SNR is particularly well suited for this purpose: in fact HB9 is known to be an underluminous SNR, so that the density of particles and their spatial gradients near the source could not be large.W28 is a relatively old SNR, hence possible non-linear transport would most likely be due to the particles injected in the past.Perhaps not surprisingly we cannot claim the identification of suppression of diffusivity in either case.
Clearly, this does not imply that self-generation is not at work around SNRs, but rather that the most active regions that might be more suitable for this type of investigations are too complex and we have been unable, so far, to identify plausible candidates.If such regions do exist and the selfgeneration is at work, it is meaningful to ask what the implications would be in terms of grammage accumulated by CRs.While this problem was addressed in earlier studies (D'Angelo et al. 2016(D'Angelo et al. , 2018;;Recchia et al. 2022), here we focused on the possible presence of MCs in the circum-source regions, possibly embedded in a medium that can be either dilute and ionized or denser and quasi-neutral.We find that the presence of MCs can substantially enhance the grammage accumulated near the source.In a realistic situation in which the MC is embedded in a partially ionized medium, the effect of ionneutral damping is to reduce the level of self-generation for low energy particles, while self-generation remains effective for particles in the range 10 2 − 10 4 GeV.Such a type of situation, if sufficiently common, would induce a feature in the same range of energy per nucleon in the secondary to primary ratios.Given the implications that this might have for the global CR transport, this phenomenon would definitely benefit from a dedicated investigation.

Figure 1 .
Figure1.Hadronic γ-ray emission from MC R2 around HB9 for the 1D (left panel) and 3D (right panel) phenomenological models.The different curves refer to different combinations of values of the suppression factor for the diffusivity, the distance of the cloud from the SNR, and the mass of MC. ξ = 9% is adopted in all panels.

Figure 2 .
Figure2.Hadronic γ-ray emission from MC R2 around HB9 for the 1D phenomenological model, assuming a minimum energy of the escaping particles as indicated in the legend.Two injection spectra have been adopted, with slope s α E = 2 and 2.4.

Figure 3 .
Figure 3. Phenomenological models for the hadronic γ-ray emission from MC A (left panels) and B (right panels) of the W28 complex.The different curves refer to different values of the suppression factor for the diffusivity, different value of the distance of the cloud from the SNR and different masses of the MC.The top panels illustrate the situation in the 1D models, while the 3D cases are shown in the bottom panels.

Figure 7 .
Figure7.CR grammage accumulated near a SNR for n i = 0.01 cm −3 (left panel) and n i =n n =1 cm −3 (right panel).The radius of the SNR is assumed to be 30 pc in the low density case (left panel) and 12 pc in the high density case (right panel).The MC is approximated to be a cylinder with radius 30 pc and centered at z = z MC .The free escape boundary is assumed to be 100 pc away from the source.

Table 1 .
Fitting Parameters for the linear phenomenological scenario

Table 2 .
Fitting Parameters for the nonlinear self-generated turbulence case