{'fraction_non_alphanumeric': 0.03747498635619429, 'fraction_numerical': 0.015038505851676673, 'mean_word_length': 4.071340713407134, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-linear Cosmic Ray propagation close to the acceleration site

Recent advances on gamma-ray observations from SuperNova Remnants and Molecular Clouds offer the possibility to study in detail the properties of the propagation of escaping Cosmic Rays (CR). However, a complete theory for CR transport outside the acceleration site has not been developed yet. Two physical processes are thought to be relevant to regulate the transport: the growth of waves caused by streaming instability, and possible wave damping mechanisms that reduce the growth of the turbulence. Only a few attempts have been made so far to incorporate these mechanisms in the theory of CR diffusion. In this work we present recent advances in this subject. In particular, we show results obtained by solving the coupled equations for the diffusion of CRs and the evolution of Alfven waves. We discuss the importance of streaming instabilities and wave damping in different ISM phases.


Introduction
We study the non-linear diffusion of a population of Cosmic Rays (CRs) after they escape from the acceleration site (for example, a supernova remnant shock). As a zero order approximation the problem can be described by an impulsive injection of CRs of a given energy at a given location in the Galaxy. We study the transport of CRs taking into account the growth of Alfvén waves caused by particle-wave interactions, and their back-reaction on the CR diffusion. We also consider possible linear wave damping mechanisms that limit the growth of the waves. The problem is described by two differential equations that govern the evolution of the CR pressure and of the wave energy density as a function of time and distance from the acceleration site. We present numerical solutions of these two coupled equations for different values of the physical parameters involved in the problem.

The model
In our model the transport of CRs is assumed to be regulated by the resonant scattering off Alfvén waves, i.e. a CR of energy E resonates with waves of wave number k = 1/r L (E) where r L is the particle Larmor radius (see e.g. [1]). The (normalized) energy density I(k) of Alfvén waves is defined as: where B 0 is the ambient magnetic field and δ B the amplitude of the magnetic field fluctuations. According to quasi-linear theory, CRs diffuse along the magnetic field lines with a diffusion coefficient equal to [2]: which can be expressed as the ratio between the Bohm diffusion coefficient D B (E) and the energy density of resonant waves I(k = 1/r L ). Formally, quasi-linear theory is valid for δ B/B 0 ≪ 1. In this limit the diffusion perpendicular to the field lines can be neglected, being suppressed by a factor of (δ B k /B 0 ) 4 ≡ I(k) 2 with respect to the one parallel to B 0 (see e.g. [3,4]). This implies that under the condition that δ B/B 0 remains small, the problem is one-dimensional. We consider the streaming of CR to be the main source of Alfvénic turbulence. For CRs that stream along the direction of B 0 (that we consider aligned along the coordinate z) the growth rate of Alfvén waves Γ CR is proportional to the product between the Alfvén speed V A and the gradient of the pressure of resonant CRs, and can be defined as (see e.g. [3,5]): The sign indicates that only waves traveling in the direction of the streaming are excited. For dimensional consistency, the pressure of CRs P CR has to be normalized to the magnetic energy density B 2 0 /8π. On scales smaller than the magnetic field coherence length, a flux tube characterized by a constant magnetic field strength B 0 and directed along the z-axis can be considered. Two coupled equations (for the evolution of CRs and waves along the flux tube) must be solved: the diffusion coefficient of CRs of energy E (Eq. 2.2), indeed, depends on the energy density of resonant waves I(k), and in turn the growth rate of these waves (Eq. 2.3) depends on the gradient of the pressure of resonant CRs. The two coupled equations then read: where the left side in both expressions is the time derivative computed along the characteristic of excited waves: The advective terms V A ∂ P CR /∂ z and V A ∂ I/∂ z are neglected in the following since we verified that they play little role in the situation under examination. The last two terms in Eq. 2.5 account for possible mechanisms of wave damping, operating at a rate Γ d , and for the injection Q of turbulence from an external source (i.e. other than CR streaming). In this work we consider only linear damping mechanisms, for which the damping rate is independent of space and time. The term representing the external source of turbulence can be set to Q = 2Γ d I 0 , so that when streaming instability is not relevant the level of the background turbulence is at a constant level I = I 0 .
As already pointed out by [6], the approach described above decouples the process of acceleration of particles, which operates, for example, at a SNR shock, from the particle escape from the acceleration site. Though such a separation might seem artificial, the problem defined above can still be useful to describe the escape of particles form a dead accelerator, in which the acceleration mechanism does not operate any more, or operates at a much reduced efficiency [6]. This situation would probably apply to the case of old SNRs. On the other hand, [7] suggested that Equations. 2.4 and 2.5 could be also used to describe an intermediate phase of CR propagation in which the CRs have left the source but are not yet completely mixed with the CR background. For the case of a SNR shock, the equations above could thus be applied to those CRs characterized by a diffusion length large enough to decouple them from the shock region. Typically, this happens during the Sedov phase to the highest energy particles confined at the SNR shock when the diffusion length D B /u s gradually increases with time up to a value larger than some fraction χ of the SNR shock radius R s , where u s is the shock speed and χ ≈ 0.05...0.1 [8,9].
In both scenarios, the initial conditions for Eqs 2.4 and 2.5 can be set as follows: where a is the spatial scale of the region filled by CRs at the time of their escape from the source. Following [6], we introduce the quantity Π, defined as: where: To understand the physical meaning of Π, consider the initial setup of the problem, in which CRs are localized in a region of size a. The CR pressure within a is P 0 CR , and then Φ CR = P 0 CR a. The initial diffusion coefficient outside the region of size a is equal to D B /I 0 . The problem of transport is characterised by three relevant timescales: i) the growth time τ CR , ii) the time it takes the CR cloud to spread due to diffusion τ di f f ≈ a 2 /D ≈ a 2 I 0 /D B , and iii) the damping time τ d = 1/2Γ d . In order to have a significant growth of waves due to CR streaming, the timescale for growth must be shorter than the two other timescales: τ g < min(τ di f f , τ d ). In terms of the parameter Π this conditions reads: It is evident, then, that the parameter Π regulates the effectiveness of CR-induced growth of waves. In the absence of a damping term for waves, Π > 1 is a necessary condition for streaming instability to be relevant, while in the presence of efficient wave damping, a more stringent condition on Π applies.
CRs play no role in the generation of Alfvén waves, and Equation 2.4 can be solved analytically: .11 is referred to as test-particle solution. When wave growth cannot be neglected, the solution deviates from this analytic test-particle solution. In the extreme scenario Π >> max(1, τ di f f /τ d ) waves grow so quickly that the the diffusive term in Eq. 2.4 becomes negligible when compared to advection term V A ∂ P CR /∂ z, and the advection term cannot be neglected anymore. This describes a situation in which CRs are "locked" to waves and move with them at a velocity equal to V A [10]. An identical result was found by [11] in a study of the escape of ≈ MeV CRs from sources, and also suggested by [12].

The method
In order to solve Eqns. 2.4 and 2.5 it is convenient to perform the following change of coordinates [6]: and re-normalise the parameters as follow: In these notations (and after dropping the advection terms) the equations become: Their solution (i.e., P and W as a function of the variables τ and s) depends only on three parameters: the initial values P 0 CR , W 0 , and Γ ′ 0 . Note that: i) P 0 CR = Π, ii) W 0 = V A a/D 0 , and iii) Γ ′ = Γ ′ 0 (since we consider here only a linear damping constant in both time and space). In these notations, the condition for effective growth of waves is Π > max(1, Γ ′ W 0 ). , as a function of the CR energy E. Different colours refer to two different ISM phases: warm ionised medium (WIM, red), and warm neutral medium (WNM, blue). Right: parameter Π as a function of the CR energy for a WIM and WNM (solid lines). In absence of efficient wave damping, values of Π < 1 correspond to a test particle solutions, since the wave growth rate is inefficient. Conversely, for Π larger that Π max ∼ 10 the growth rate is large and the diffusion coefficient becomes smaller than the Bohm diffusion coefficient D B . In presence of efficient wave damping, both these limits are shifted to higher values, depending on the value of Γ d .
We present examples of numerical solutions to the coupled equations 3. . In order to understand why we chose these values, and to which physical conditions they correspond, it is necessary to specify the values of V A , a, D 0 , D B , and the damping mechanisms. These quantities in general depend on the properties of the ISM, and on the particle energy. We consider two different ISM phases: a warm neutral medium (WNM) and a warm ionised medium (WIM). The typical values of density, temperature, ionisation fraction, and magnetic field strength for these two phases are taken from [13]. In these ISM phases, the two most relevant linear damping mechanisms are the ion-neutral friction (Γ IN d ) and the damping by background MHD turbulence suggested by Farmer & Goldreich [14] (Γ FG d ). The value of Γ FG d as a function of the particle energy is shown in Fig. 1 (left panel, long-dashed lines), for both ISM phases. The ion-neutral damping rate is estimated by numerically solving the equations in [15]. The results are shown in Fig. 1 (left panel, solid lines), for both ISM phases: at low energies The level of coupling between ions and neutrals affects also the Alfvén speed, that is in the range 10 6 − 10 7 cm/s. Since it is useful for our calculations, we show in Fig. 1 the value of the ratio V a /a, as a function of energy (dash-dotted lines). In order to estimate a we consider a ≈ R esc (E), where R esc (E) is the escape radius and depends on the CR energy. We model its dependence on the energy following [16]: we found that its value ranges between a few pc (for the highest energy particles) and ∼20 pc (for ∼GeV particles). Since we are interested in following the CR diffusion from the scale a where they are released up to a distance z 10 3 pc, this roughly corresponds to solve the normalised equations in a range of s from 1 to several hundreds. Moreover, we want to study the CR pressure at different times t in the range 10 3 yr to 10 5 yr. Considering the value of V a /a as a function of energy displayed in Fig. 1, this corresponds to τ = 10 −2 − 1.
In order to identify the relevant values of the parameters Π,W 0 , Γ ′ we proceed as follow: ranges between 10 −9 and 10 −12 s −1 (Fig. 1, left  panel), depending on the particle energy and ISM phase. This corresponds to 1 < Γ ′ 10 4 .
• Π -the values of Π as a function of the CR energy, for the two different ISM is shown in Fig. 1 (right panel, solid lines).
• W 0 -for the background diffusion D 0 we can refer to the average Galactic value D 0 = D ISM = 10 28 (E/10 GeV) 0.5 cm 2 /s. Since V A a ranges between 10 27 cm 2 /s and 10 25 cm 2 /s (at low and high energies, respectively), the range of relevant values for W 0 is 10 −5 < W 0 < 10 −1 .

Results and Discussion
We numerically solve Eqs. 3.3 and 3.4 and show the results for P(s, τ) vs. s (solid lines in Fig. 2) at three different normalised times τ = 10 −2 , 10 −1 , 1 (corresponding to different colours in the figure). The test particle solution is also shown for comparison (dotted lines). Each row in the figure corresponds to results obtained by varying only one parameter at the time.
The upper row shows the impact of the parameter Π for fixed values of Γ ′ = 10 and W 0 = 10 −4 . Since for Π < 1 the numerical solution coincides with the test particle solution, we show our results for values of Π > 1. When Π = 2, the solution starts to deviate from the test particle case: due to the growth of wave by streaming instability the diffusion is slower, the CR pressure in the vicinity of the accelerator is larger, and the diffusion length is smaller. The effect is more and more evident for increasing values of Π (upper row, middle and right panels).
The middle row shows the effect of wave damping, for fixed Π = 10 and W 0 = 10 −4 . Since Π >> 1, when wave damping is negligible (Γ ′ =0, left panel) the diffusion is suppressed as compared to the test particle case. Larger values of Γ ′ limit the growth of Alfvén waves and the solution approaches the test particle solution. Note that Γ ′ =10 (middle row, central panel) corresponds to a damping mechanism that becomes important at τ > 0.1. For this reason the solution at τ = 0.01 and τ = 0.1 are unaffected by the increased damping rate, and only the solution at late time τ = 1 is modified, and approaches the test particle solution. The right panel shows the case of a even faster damping mechanism, Γ ′ =100, whose effect starts to be relevant at times corresponding to τ > 0.01, explaining why the solution at intermediate and late times approaches the test particle solution, while the solution at early time is unaffected by wave damping.
The bottom row shows calculations performed at Π = 5 and Γ ′ =10, with W 0 varying from 10 −5 to 10 −2 . A large W 0 (left panel) corresponds to slow diffusion and small diffusion length, both in the test particle case and in the numerical solution. The two solutions differ at early and intermediate times because the growth rate is faster than the diffusion rate, and suppresses the diffusion. At later times instead the growth rate decreases and the numerical solution does not appreciably differ from the test particle case. A change in the initial value of the background turbulence has the same effect on both the analytical and numerical solution: the diffusion is faster for decreasing values of W 0 (see bottom row, from left to right), and the CRs can reach larger distances. The deviation from a P CR P P CR CR Figure 2: Normalized CR pressure as a function of the variable s, which represents the distance from the source z normalised to the size of the source a. In each row the results are shown by varying only one of the three parameters (whose value is reported above each panel) and keeping the other two fixed to the values reported in the shaded area. The solid curves show the results from the numerical computation, while dotted curves represent the test particle solution, where the role of streaming instabilities in enhancing the turbulence is neglected. Different colours refer to different normalised times τ = 10 −2 , 10 −1 , 1. test particle solution however remains unchanged, since both the growth rate and the diffusion rate scale as W 0 .
This study outlines the importance of considering, in the physics of CR transport, both the growth of waves due to streaming instability and the role of mechanisms that damp the turbulence. When damping is neglected, even relatively small values of the integrated CR pressure Π = 10 (corresponding, for example, to particles of a few TeV, in a WIM -see Fig. 1) in a environment with small background turbulence W 0 = 10 −4 (relevant for TeV particles when the background diffusion coefficient D 0 is similar to the average Galactic diffusion coefficient) lead to the conclusion that diffusion is strongly suppressed by the efficient growth of waves, and particles cannot travel large distances, even at late times. This situation is shown in Fig. 2 (middle row, left panel) where the numerical computation is performed in absence of a wave damping term: the diffusion is strongly suppressed also at τ = 1, corresponding to old remnants (t age ∼ 10 5 yr).