Uchiyama et al. (2012) investigate the GeV gamma-ray emission detected by the Fermi-LAT from molecular clouds in the vicinity of the SNR W44. The CR model is adapted from Ohira et al. (2011). CRs are injected in a burst-like manner. Reduced diffusion solutions are possible with χ taking values down to 0.1 (notice that in this study δ was kept fixed to 0.6). Peron et al. (2020) propose an updated analysis of W44 and its surroundings using 9.7 years of Fermi data. Interestingly, their analysis shows that two clouds lying along the galactic plane have enhanced gamma-ray emission that is incompatible with a contribution of the local CR background. These observations likely necessitate anisotropic CR diffusion from W44. Complementary observations are needed in order to infer constraints on the values of the diffusion coefficients κ s , ‖ and κ s , ⊥ .

Using both GeV and TeV data and an isotropic spherical diffusion model with burst-like injection, Hanabata et al. (2014) propose χ to be in the range 1/20–1/2 to explain the gamma-ray profiles around W28, while κ s is found to have a mild dependence with the CR energy E 0.1 − 0.3 .

These two first studies triggered several others with diverse conclusions: The MAGIC collaboration (MAGIC Collaboration et al., 2023) fits extended emission around the γ Cygni SNR and finds (spherical model, burst-like injection) χ in the range 0.03–0.1; Li et al. (2023) using Fermi data find (spherical model, burst-like injection at an energy-independent time) solutions with χ < 1 and δ in the range of 0.5–1 around G15.4 + 0.1 SNR; using a similar model, Xin and Guo (2023) using Fermi data find results consistent with χ ∼ 1 around DA 530 SNR; Oka and Ishizaki (2022) arrive (spherical model, burst-like injection) to similar conclusions in the case of HB 89 SNR. Finally, Katagiri et al. (2016) used Fermi data to find (spherical model, continuous injection) an enhanced CR energy density by a factor of approximately 10 with respect to local values in the environment of the HB3 SNR, in particular in the W3 source associated with an HII complex in interaction with the SNR.

Aharonian et al. (2019) use Fermi, HESS, and ARGO data to investigate the extended emission around the Westerlund 1 massive star cluster, around the Arches, Quintuplet clusters in the central molecular zone, as well as the Cygnus cocoon. In each case, the radial profile of the gamma-ray emission shows a 1/r decline within 50 pc around the central source. The objects also show an energy spectrum scaling as E − 2.2 without any noticeable cut-off above 10 TeV. Fixing a maximum value of the overall wind mechanical power W max and a maximum acceleration efficiency for protons above 10 TeV to not exceed 0.01, it is possible to derive an upper limit of the diffusion, still above 10 TeV, not exceeding 1 0 28   c m 2 / s , giving χ ≤ 0.01 .Recent investigations propose some re-analysis of some of these clusters and obtain more mitigated results. Aharonian et al. (2022), using HESS observations of Westerlund 1, do not show strong evidence of CRs leaking out of the cluster region. Adopting a distance of 3.9 kpc to fix the CO map around the object, they still obtain a density radial profile compatible with the results by Aharonian et al. (2019). Abeysekara et al. (2021) propose a deep analysis of the Cygnus cocoon combining Fermi and HAWC data. If a GeV gamma-ray profile reproduces a 1/r profile well, a TeV gamma-ray profile may be well explained by a constant CR density produced by a bursting event. Constraints on the diffusion coefficient still conclude a suppression of CR diffusivity with 0.002 ≤ χ but still ≪ 1 .

Gamma-ray halos around evolved pulsars are a new phenomenon detected in the GeV and multi-TeV bands. This shows extended diffuse structures of several pc scales; hence, they are much extended compared to the pulsar wind nebula. López-Coto et al. (2022) propose a review of the gamma-ray halo phenomenon based on the discovery of Geminga and Monogem pulsar halos. Fang (2022) proposes a complement to this list, including two more objects, namely, LHAASO J0621 + 3,755 [see also Liu (2022)] and HESSJ 1831-098. For each of these objects, the gamma-ray profile centred on the pulsar and its nebula (whose extension is small with respect to the halo), a spherical gamma-ray profile either in the GeV or multi-TeV range fits a 1 / r law well. The particle injection is time dependent with an amplitude scaling with the pulsar spin-down power evolution. Typical reduced particle diffusivity can be deduced for this list of objects in the range 5   1 0 − 3 (for LHAASO J0621 + 3,755) to 0.03 (for Monogem); see Fang (2022) and the references therein.

In this short review, we first detail in Section 2 the kinetic models describing cosmic-ray escape from their sources and then present the cosmic-ray cloud (CRC) model. We consider three types of sources: supernova remnants, massive star clusters, and pulsars through their gamma-ray halos. Then, Section 3 describes the potential dynamical effects associated with CR pressure and current effects. The effects due to ionisation are also invoked in this section. Section 4 concludes and proposes some perspectives.

Modified CR propagation is usually invoked in the framework of CR self-generated turbulence; that is, CRs can produce the turbulence that they scatter off. However, specific local properties of a background (or extrinsic) turbulence can also explain, to some extent, the gamma-ray and CR halo phenomena. We now examine these two possibilities in more detail.

The first step for the production of a CR halo comprises the escape of these particles from the accelerator. One could define a CR halo as a cloud of energetic particles confined in the close spatial and temporal environment of a source (the spatial and temporal extension of the cloud need to be properly estimated) but disconnected from the acceleration process [see Malkov et al. (2013) for more discussion]. Let us illustrate the properties of a halo in the case of SNRs. The above definition is, in fact, a bit arbitrary. Indeed, in the diffusive shock acceleration process, likely at the heart of the acceleration of CRs, a more correct approach would be to consider that the most energetic particles have a 3D random walk and a probability to return to the shock, which is dropping to 0, as discussed by Drury (2011). This description, however, requires a minimum level of extrinsic turbulence to allow particles close to the maximum energy to not freely escape; otherwise, the escape process is controlled by the energy flux carried by the highest CR to bootstrap the acceleration process (Bell et al., 2013); (Blandford et al., 2023). However, 3D modelling of the escape process is out of reach of current computational resources. The authors are reduced to adopting a 1D spatial description in which the escape process is controlled by a physical boundary (named the far escape boundary) or even simply setting a maximum energy; see the discussion in Reville et al. (2009). Physically, this boundary may mark the transition between a zone where the presence of CRs triggers a magnetic dynamo and the interstellar medium (Xu and Lazarian, 2022b). Notice that CRs can also escape from downstream if, for instance, the turbulence level drops (Ptuskin and Zirakashvili, 2005.

Several models have attempted to describe the escape process and the transition from the in-source to around-source CR distribution. These models focus on SNRs. In these approaches, the maximum CR energy E max or momentum p max is either fixed or parameterised as a function of time and in some cases, derived from source modelling. We describe some of them below. The main assumptions are resumed in Table 1 for conciseness.

Ohira et al. (2011) consider the case of SNRs close or in interaction with molecular clouds (MCs). They use a 1D spherical model and introduce a parameterised form for p max = p 0 ( t / t Sed ) − α , where p 0 = 1 PeV/c, and t Sed is the Sedov timescale and α = 6.5 . However, this value can be modulated in the case that the SNR interacts with an MC directly because particles can escape due to the damping of magnetic perturbations by ion-neutral collisions. The authors apply their model to four specific SNRs: W51C, W44, W28, and IC443. They derive the effective CR spectrum, which is marked by a succession of energy spectral breaks quoting the population of CRs reaching the inner and outer bound of the MCs. In order to fit the gamma-ray data, a reduced diffusion coefficient is necessary in all cases. Makino et al. (2019) continue with similar modelling but also derive the contribution of MeV protons to the iron 6.4 keV fluorescence line if the SNRs interact with a molecular cloud.

Nava and Gabici (2013) propose a similar model using 2D cylindrical geometry. In this model, the maximum CR energy is also parameterised as by Ohira et al. (2011) but with p 0 = 5 PeV/c. Particles are injected in a magnetic flux tube in the (x, y) plane at z = 0 within a circular region given the SNR radius corresponding to the escape time of particles at a given energy E. The model is 2D because it includes the effect of magnetic field random walk due to magnetic, turbulent, chaotic motions. The diffusion coefficient is parameterised with the two above-defined parameters, χ and δ (fixed to be 0.5). The authors show that because anisotropic diffusion concentrates the CR flux in flux tubes, it can alleviate low χ values (the authors use normal values with χ = 1 ) when it matters to reproduce the gamma-ray signal from MCs surrounding the W28 SNR.

Yan et al. (2012) propose a 1D Cartesian model for the CR precursor ahead, the SNR shock. They accurately derive a value for p max based on a physical model accounting for magnetic amplification by the CRs themselves. They treat CR propagation around the SNR using a 1D spherical approach where the diffusion coefficient is scaled down by a factor χ with respect to its ISM values. The CR distribution is obtained using the solution of a stationary 1D diffusion. Wave growth and damping are only derived in the CR precursor region. The χ factor is constrained by comparing the model with gamma-ray data from W28, leading to a reduced diffusivity with χ ≃ 0.05 .

Fujita et al. (2011) adopt a 1D spherical two-zone model: a shock zone corresponding to a CR precursor and an escape zone where particles at high energy move diffusively with a diffusion coefficient derived from quasi-linear theory (or QLT, see next). The coefficient includes the amplitude of the self-generated waves. In each zone, an equation for the self-generated waves is solved using a quasi-linear expression for the growth rate. The authors consider three types of models: A, B, and C. Model A corresponds to the general case where the amplitude of the waves is derived in parallel to the evolution of the CR distribution; models B and C are used as tests. In model B (C), the diffusion coefficient is fixed to its Bohm value (the amplitude of the wave is fixed). In all cases, there is no explicit calculation for p max , which is likely fixed. The authors find that wave self-generation can strongly limit CR expansion in a zone around 2 R s , where R s is the SNR radius. The extension of the CR distribution in model A lies between the two other models.

Telezhinsky et al. (2012) include a treatment of particle acceleration at the SNR shock in 1D spherical geometry, including hydrodynamics of fluid evolution. The model can hence derive a physically based expression for p max . For this, the authors define different diffusion regions (their D2 model): one region around the shock with radii between 0 and R s h ( 1 + ℓ ) with ℓ = 0.05 where Bohm diffusion κ = κ B = E q B (B is the total magnetic field accounting for magnetic field amplification) is adopted, a transition region where κ = χ κ 0 , and then an outer region beyond 2 R s h where κ = κ 0 . Hence, p max is fixed by comparing the spectra of escaped CRs and diffusion models. The authors consider two scenarios for the interacting material producing gamma rays: either a molecular cloud in the case of a type Ia SNR or a dense wind-blown shell in the case of a core-collapse SNR. In a recent work, Brose et al. (2020) further elaborate on Telezhinsky et al. (2012) calculations, including 1D hydrodynamic solutions for the thermal gas. The calculation is performed in the test-particle limit and so does not account for CR pressure effects. The large-scale magnetic field is fixed up- and downstream of the SNR shock. The CR injection and transport are extended to include the SNR radiative phase. The authors also include an explicit treatment of the wave amplitude, including self-generation and damping, following Brose et al. (2016). They also include an explicit calculation of the CR electron spectrum. Finally, the authors allow for a time-dependent injection of CRs at the SNR forward shock; it scales as a power-law function over time. Figure 1 illustrates the propagation modes of CRs near their sources, which is interesting in terms of the purpose of this review. Figure 1up shows the CR solutions over the whole volume of integration: around the shock and upstream, as defined in Telezhinsky et al. (2012). It illustrates the impact of the CR transport over the maximum energy reachable by the particles. Bohm diffusion maximises particle confinement and leads to higher energies than in the case of self-generated wave solutions. The latter model also produces an energy break around 10 GeV at late evolution timescales where the energy CR spectrum passes from a power law with an index of 2–2.4, consistent with the required injection spectrum for galactic propagation models. Figure 1(down) shows the spectra of escaped CRs at distances of 7.5 pc from the SNR shock in a 2-pc-thick layer at a time of 37 kyr. Four models are displayed: Bohm and self-generated transport with or without the inclusion of a background CR population (not shown in Figure 1). The main impact of background CRs is to compensate for the CR normalisation decrease at low energy because of reacceleration. The self-generated solutions lead to less CR pressure but a wider distribution than the Bohm case, where the escaped particles are more peaked and populate higher energies. This difference in shapes is explained by the fact that the maximum particle energy evolves more rapidly with time in the case of the self-generated wave model.

Yang et al. (2015) also elaborate on Telezhinsky et al. (2012) but go beyond the test-particle limit by including non-linear effects over magnetohydrodynamics of shock solutions following the numerical recipes by Zirakashvili and Ptuskin (2012). The authors propose a calculation of the CR distribution outside the acceleration zone by modulating the escaped CR distribution (so a distribution based on a non-linear acceleration model) by the kernel obtained from the 1D spherical solution of the diffusion equation. This is their model B. Their model A is similar to the approach adopted by Telezhinsky et al. The main trend of the results, regardless of age and distance (up to 100 pc ahead of the SNR), is that model A leads to a moderately enhanced CR component with respect to model B, which includes non-linear effects in the source. The difference is, at most, one order of magnitude in the energy flux.

A recent analysis considered the escape of CRs from SNRs propagating in the Parker spiral produced by a massive progenitor star (Kamijima and Ohira, 2022). Test-particle solutions show that the energy of escaping particle E max is partly controlled by the electrostatic potential difference between the pole and the equator Δ ϕ , with E max ≥ e Δ ϕ . Red supergiant stars can contribute to producing CRs above 10 TeV if Gauss-level magnetic fields can be maintained in the wind. Fast-rotating Wolf–Rayet stars can contribute to producing CRs above 1 PeV. The same study in an ISM magnetic field, when adapted to type Ia SNRs, gives maximum energies of several tens of TeV (Kamijima and Ohira, 2021).

At this stage, an alternative approach combines acceleration, escape, and transport in the ISM. In this model, Petrosian and Chen (2014), instead of adjusting data with model parameters, revert the analysis by using CR data from sources and direct measurements to infer the CR momentum diffusion coefficient and angular scattering frequency. The authors concentrate their analysis on electron spectra (for instance, from SNR RXJ1713.7–3946 or SN 1993J). They conclude that it is likely that complex relationships beyond a simple random walk should link escape and scattering times, possibly invoking the presence of stochastic interaction with compression perturbations, as discussed in Section 2.3.8.

Summarising briefly, the aforementioned models all tend to show that an SNR can inject a CR distribution harder than the background ISM population over relatively long timescales of a few tens of kyr. When self-generated waves are included, the confinement of particles is reduced with respect to the case of prescribed transport [see, e.g., Brose et al. (2020)], that is, by fixing a reduced factor χ . In all cases, MCs can track CRs and are good test targets to probe the escaping CR distribution.

We now turn to the cosmic-ray cloud (CRC) model, which is a different class of model. These are 1D Cartesian. They calculate the propagated CR spectrum along a magnetic flux tube surrounding the mother source and provide a parallel solution of the self-generated wave spectra. In some senses, CRC models are less accurate about the escaping process from the accelerator but propose a more accurate treatment of the propagated particles once these are disconnected from the acceleration process. Notice that assuming a population of particles can be completely disconnected from acceleration processes occurring at the SNR is the main assumption of the CRC model. The assumption somehow hides the physics of escape and, hence, the transition from in-source to near-source propagation. This subject is still open, and the above-described models only treat this aspect in a phenomenological way. In reality, this issue requires a refined modelling of the microphysical processes responsible for this transition [see discussions in Malkov et al. (2013) and Bykov et al. (2017)].

The cosmic-ray cloud model has been developed by Ptuskin et al. (2008) and Malkov et al. (2013). In this model, CR acceleration and escape phases are decoupled, and acceleration is not explicitly treated separately in the construction of the CR injection term Q C R .

The gamma-ray halo phenomenon briefly described in Section 1 has also been interpreted in the framework of CR self-confinement. In that case, self-generated waves are produced by electron–positron pairs injected from the pulsar termination shock and the pulsar wind nebula, although other interpretations are possible [see (Giacinti et al., 2020)]. Evoli et al. (2018) derive the theory of resonant streaming instability in the context of gamma-ray halos. In this model, electron–positron pairs trigger magnetic fluctuations following k r L ∼ 1 . Pairs are injected with a power-law distribution f ( p ) ∝ p − a . The growth rate of the resonant modes can be derived using an alternative form to the one proposed in Equation 4, which implicitly assumes particles are diffusing. The growth rate can be expressed in terms of the drift speed of the pair beam V d : Γ g ≃ π 2 a − 3 a − 2 V d V a n e e n g Ω c p k r l , 0 a − 4 ,

where r l , 0 is the Larmor radius of the particles with the minimum momentum in the distribution. The solution is obtained by solving a similar system to Equations 1, 2. The nonlinear Landau damping process (see Equation 1) is the damping process retained in this work. The electron–positron diffusion coefficient is obtained by solving Equation 3. Mukhopadhyay and Linden (2022) use a similar approach but include supplementary damping processes in the wave equation, namely, ion-neutral damping and turbulent damping. The authors added two more ingredients to the model: a 3D description of the propagation that is still in spherical geometry and the presence of an SNR, which contributes by injecting cosmic rays (electrons and protons) as supplementary sources of self-generated waves. The authors consider two models for their 3D runs. Both have a pulsar injection efficiency of 100% into pairs (or α = 1 ). The “tuned” model has a harder injection distribution with a = 3 instead of a = 3.5 in the fiducial case. This model also has a smaller spin-down time and a smaller background magnetic field strength. Figure 10 presents the solutions of the electron-positron diffusion coefficients as a function of time at 10 TeV (at energies corresponding to the gamma rays probed by HAWC). The solutions clearly show that reducing the diffusion coefficient by several orders of magnitude is challenging in this framework. The maximum reduction occurs close to the pulsar, at 5 pc. The addition of the turbulence generated by the particles injected at the SNR can help in reducing the coefficient further but requires the transfer of much of the SNR mechanical energy into protons, about 40%. Mukhopadhyay and Linden (2022) also propose a 1D model where a diffusion coefficient reduction by more than two orders of magnitude is recovered, although this model predicts diffusion suppression over scales larger than 100 pc. The latter effect partly motivated the above-mentioned 3D developments. From this analysis, it seems difficult to fully account for the gamma-ray halo phenomenon in the framework of self-generated waves or, at least, it requires other sources of turbulence, such as the one generated by the SNR. However, self-generated wave models can be invoked if the electron–positron beam is injected after the pulsar birth if 1D propagation can be sustained, which requires specific properties of the background turbulence (Malik et al., 2023). As stated in Section 2.3.8, another possible explanation for pulsar halos is to have multiple processes contributing: self-generated waves close to the pulsar wind nebula, other sources of turbulence farther away, and specific magnetic geometries. Time-dependent effects should also be important in shaping the halo.

This instability has been considered the main instability driving magnetic field turbulence in the ISM (Morlino, 2018; Blasi, 2019; Marcowith et al., 2021). Through this instability, CRs transfer their momenta into the background gas via the production of Alfvén waves and ultimately participate in heating the ISM (Wiener et al., 2013). Once they can trigger resonant magnetic fluctuations, CRs are rapidly locked to the gas (Skilling, 1971). Their bulk speed is the ambient Alfvén speed, or f V a , with f > 1 being a boost factor fixed by the main wave damping process (Ruszkowski et al., 2017).

The main application of this theory concerns the CR feedback in star formation through the sustaining of galactic winds. In effect, while driving the background gas in motion, CRs extract some gas from the galactic disk, lowering the star formation rate [see the review by Ruszkowski and Pfrommer (2023) and references therein]. This review is, however, restricted to smaller-than-kpc-scale ISM. The interested reader is referred to the above-mentioned review for details about CR-driven physics at larger galactic scales.

Around sources, the pressure gradient set by the onset of the CR streaming instability can have dynamical effects over the ambient gas because ∇ ⃗ P CR is a force term in the gas momentum equation. Commerçon et al. (2019) have used CR-magnetohydrodynamic simulations to investigate the necessary condition under which CRs can become trapped with gas at meso-galactic scales, close to their sources, in practice around SNRs. In the ISM, turbulence imposes some scaling relations, known as the empirical Larson laws (Larson, 1979; Larson, 1981), linking the turbulent velocity σ and the size of the probed turbulent region. Namely, we have σ = σ 1   p c L q , where q = 1 / 3 corresponds to a Kolmogorov turbulence, and σ ( 1   p c ) is the turbulent rms speed at 1 pc. In ionised gas, a recent investigation using the WHAM survey dataset improves the estimated density fluctuations at large scales (Chepurnov and Lazarian, 2010). Improving the lever arm of the data permits a better fit with a Kolmogorov spectrum, with q = 1 / 3 . Lines broadening studies in the partially ionised atomic gas allow, for instance, deriving q ≃ 0.35 and σ ( 1   p c ) ≃ 1 − 1.5 km/s (Roy et al., 2008). In molecular gas, we have q ≃ 0.5 and σ ( 1   p c ) ≃ 1 km/s (Hennebelle and Falgarone, 2012). If CRs scatter some magnetic fluctuations leading to a typical diffusion coefficient κ along the mean magnetic field lines, then, for a region of size L, CRs can be trapped in the gas if t diff = L 2 / κ > t turb = L / σ ( L ) . Hence, the diffusion coefficient must be limited to κ ≤ σ L L . (8)

In that situation, the CR gradient along the magnetic field lines is high enough to produce some feedback effects. This condition is true regardless of the ISM phase. If κ fulfils the condition in Equation 8, then the gas is entrained by the streaming CR motion, and the local sound speed is also modified as the effective pressure now combines the gas and the CR pressure. All macro instabilities that are important for ISM dynamics (Parker–Jeans, Rayleigh–Taylor, Kelvin–Helmoltz, and magneto-rotational) have a modified growth rate partly controlled by κ [see Marcowith (2023)]. In particular, close to SNRs, Commerçon et al. (2019) estimate that the CR gradient sustained by the presence of an SNR must exceed a typical ISM value of P C R , r e f ′ ≃ 5 × 1 0 − 33   d y n e s / c m 2 P CR 1   e V / c m 3 L C R 100   p c − 1 , where the prime denotes the gradient taken along the mean magnetic field direction, and L CR is the CR gradient length, here set with respect to the galactic disc height. The triggering of the resonant modes is, in principle, an efficient way to reduce κ and to increase P CR ′ above P C R , r e f ′ . A rough calculation balances the resonant streaming growth rate with the dominant damping rate and fixes the expected level of turbulence. It is found that for an SNR propagating in the atomic phase of the ISM (WNM, CNM), the ratio P CR ′ / P C R , r e f ′ ≥ 100 for GeV CRs. Hence, one should expect CRs recently released in such an ISM phase to possibly have substantial feedback because of pressure gradient effects. This simple estimate, however, requires that dedicated numerical simulations be confirmed or invalidated.

Schroer et al. (2022) and Schroer et al. (2021) have considered a different scenario than the ones invoked in the framework of the CRC model. In their setup, CRs escaping the SNR can carry a current J CR strong enough to trigger the non-resonant branch of the steaming instability (Bell, 2004). In that case, magnetic perturbations grow due to the return shielding current ported by thermal background electrons. The non-resonant streaming instability grows faster than the resonant one, but non-resonant modes only grow from wavenumbers k > r l − 1 to k max = B 0 J CR 2 ρ V a 2 , with a maximum growth rate Γ max = k max V a reached at k max / 2 [hence the non-resonant character of the instability; see Amato and Blasi (2009)]. Here, B 0 and ρ are the background magnetic field amplitude and gas density, respectively. In this model, CRs are released in a single burst from the SNR. CR distribution follows a power-law distribution with f ( p ) ∝ p − 4 . The condition for triggering the non-resonant branch of the instability is (Zweibel and Everett, 2010; Bell, 2004) U CR U B 0 > c V D , (9)

which translates into k max r L > 1 . This is the necessary condition for this instability to become destabilised. However, even if the non-resonant modes are strongly driven, the instability growth time must be short enough with respect to the residence timescale of the CRs t res , either upstream of a shock front in an SNR or in a flux tube. In the former case, this timescale is set by the CR shock precursor size, ℓ prec ≃ κ / u sh , where u sh is the shock speed, and κ the diffusion coefficient of the most energetic CRs. Hence, we have t res = ℓ res / u sh . In the latter case, the residence time is t res ≃ ℓ c / V D , and the ratio of the flux tube coherence length to the CR drift speed. Hence, the condition is Γ max t res > 1 .

In Schroer et al. (2022) and Schroer et al. (2021), 2D and then 3D hybrid simulations are conducted where CRs are injected on a subpart of the left boundary to mimic the injection by an SNR into the magnetic flux tube. Figure 11 shows the evolution of the 3D morphology of the flux tube using the hybrid approach. The simulation setup respects the condition in Equation 9, even if the CR density is upscaled with respect to a realistic object. A clear inflation due to the lateral CR pressure effect can be noticed. The diffusion coefficient along the flux tube direction converges to a value a few times its Bohm limit. This work shows that while being injected around the SNR, CRs may modify the surrounding ISM gas and magnetic distribution. This effect should be accounted for while dealing with any gamma-ray halo morphology. However, the effective expansion of the flux tube in realistic situations is not clear at this stage. In effect, if only 10% of the mechanical energy is imparted into CRs at the position of the SNR, 90% is still available to inflate the magnetic flux tube. The inflation information propagates at a speed corresponding to the maximum between the local sound and the Alfvén speed. While moving ballistically, the CRs will reach regions of the flux tube far ahead, but as they start to trigger magnetic fluctuations, their mean speed slows to V a or a few times V a , depending the main damping mechanism. Hence, the total effect due to CRs should not exceed the effect produced by the SNR expansion. Still, if CRs can be more efficiently scattered by adding the contribution of the magnetic fluctuations issued from the generation of the non-resonant instability, then one may expect an enhancement of the CR grammage around a source.

Blasi and Amato (2019) propose a similar approach for the case of CRs escaping the galactic disc. CRs escape diffusively the galactic disc with a flux ϕ CR . The flux escaping from the disc is proportional to the disc CR luminosity L CR ∼ 1 0 41 erg/s dispersed over a surface proportional to R d 2 , R d being the disc radius. The flux thus obtained is found to be high enough to sustain the onset of the non-resonant streaming instability based on Equation 9 and a background magnetic field not in excess of approximately 21 0 − 2 μ G. Once the instability grows, magnetic perturbations also grow, and CRs start to propagate diffusively. This induces a CR gradient, which drives the background gas into motion to a speed close to the Alfvén speed. The non-resonant instability can be a source of magnetic field injection into the galactic halo. Finally, CRs can produce diffuse gamma-ray and high-energy neutrino emissions while interacting with circum-galactic gas.

By CR pressure-driven feedback, we mean the induced effects associated with second-order anisotropy in the CR distribution function (Bykov et al., 2014; Zweibel, 2020). Zweibel (2020) proposes an ensatz for the particle distribution function given by F p = F 0 p + ζ p ∂ p F 0 p P 2 μ , (10) where F 0 ( p ) is the isotropic part of the CR distribution and P 2 ( μ ) is the Legendre second-order function as a function of the CR velocity pitch-angle cosine μ . The controlling anisotropy parameter is ζ = 5 12 Δ P CR P C R , where P CR = ( 2 P C R , ⊥ + P C R , ‖ ) and P C R , ⊥ , P C R , ‖ are the perpendicular and parallel CR pressures with respect to the mean magnetic field direction, respectively. The pressure excess is derived from the second term in the RHS of Equation 10. For F 0 ( p ) ∝ p − a , we directly reproduce the pressure-anisotropy-driven instability growth rate [see also Lazarian and Beresnyak (2006)] Γ C R , A k ≃ π 8 a − 2 a − 3 × Ω c p × n C R p > p 1 n g 3 a a − 2 a 2 − 1 c V a ζ − 1 ,

where p 1 = m p Ω c p / k ‖ . The CR scattering frequency in the pressure-driven turbulence can be fixed by balancing the above growth rate with the appropriate damping rate. Zweibel (2020) showed that the CR pressure-driven usually leads to a less efficient ISM heating and CR scattering frequency than the resonant streaming instability by a typical factor of V a / c . This is explained by the fact that magnetic field expansion, which occurs at a typical speed of V a , is balanced CR scattering, whereas, in the case of streaming instability, CRs tend to move along the magnetic field at the speed of light. However, this instability is interesting, as it is a source of CR work in the transversal direction of the magnetic flux tube. This instability can also be relevant close to CR sources where stronger anisotropy should be expected, such as those induced by the CR escape process (Bykov et al., 2017) (In this study, other types of CR pressure instabilities may be triggered, namely, the non-resonant firehose/mirror modes that control the CR anisotropy). Again, as for the streaming instability, it seems important to develop new numerical tools and set-ups capable of integrating different degrees of CR anisotropy to address the effect of triggering both pressure-driven and streaming instabilities (Reville et al., 2021).