Cosmic-ray current-driven instabilities -- revisiting environmental conditions

The growth of magneto-hydrodynamic fluctuations relevant to cosmic ray confinement in and near their sources, and the effects of local plasma conditions is revisited. We consider cases where cosmic rays penetrate a medium which may contain a fraction of neutral particles, and explore the possible effects of high-order cosmic-ray anisotropies. An algorithm for calculating the dispersion relation for arbitrary distributions, and anisotropies is presented, and a general solution for power-law cosmic-ray distributions is provided. Implications for the resulting instabilities near to strong Galactic cosmic-ray sources are discussed. We argue that cosmic-ray streaming in weakly ionised plasmas eliminates the need for the existence of an evanescent band in the dispersion relation, a conclusion which may be confirmed by gamma-ray observations. The necessity for additional multi-scale numerical simulations is highlighted, as understanding the non-linear behaviour is crucial.


INTRODUCTION
Young supernova remnants (SNR) are sources of high-energy particles, revealed by their non-thermal x-ray and γ-ray emission (e.g. Berezhko et al. 2002;Hinton & Hofmann 2009). If they are also sources of cosmic rays (CRs), as is likely the case, then the escape and transport in their local environment should hint at the processes determining confinement, and offer a potential probe of the time-history and conditions near the source. Such information is required to understand their contribution to the total Galactic CR population.
Studies on CR escape from sources, dating back to the early works of Kulsrud & Cesarsky (1971), Skilling (1971) and others, implicitly assume that CR transport, at least at low energies, is dominated by scattering on resonantlyexcited Alfvén waves. Close to efficient CR sources, this is however unlikely to be the case. For example, on scales relevant to the energetic particles emitting the >TeV γ-rays now detected from several young SNRs, estimates for the growth-time of resonantly excited Alfvén modes can exceed the lifetime of the SNR in question. The present understanding is that the non-linear development of a much faster nonresonant short-wavelength instability provides the required scattering and confinement (Bell 2004;Bell et al. 2013).
The concept of magnetic field amplification, and self-E-mail: brian.reville@mpi-hd.mpg.de confinement underpins all current approaches to model acceleration of CRs to the highest energies in SNRs. It will likewise be a requirement of any conceivable alternative source. The model put forward by  and Bell et al. (2013) establishes an upper-limit to the maximum achievable CR energy; directly related to the ability of escaping high-energy particles to amplify the magnetic field (non-resonantly) to a level that inhibits further escape. In this model, the system is continuously leaking a fraction of particles at the highest energies, their mean free path growing with distance from the shock. The picture bears similarities with the free escape boundary approach commonly used as a closure approximation in non-linear steady state models (e.g. Ellison & Eichler 1984;Reville et al. 2009).
This slightly artificial picture is conceptually simple in that it conveniently separates the energetic particles into two populations; one being transported in a highly developed and likely complex field topology close to the shock, and another streaming freely into an ambient medium. The former regime has been the focus of numerous numerical investigations (e.g. Bell 2004;Reville et al. 2008;Gargaté et al. 2010;Rogachevskii et al. 2012;Bai et al. 2015;Marret et al. 2021) although the inherent multi-dimensional, non-linear and multi-scale nature of the problem still poses serious challenges. For the latter scenario, capturing the transport and self-regulated scattering in the immediate surrounding medium is an equally challenging problem, as an accurate model demands connecting the global structure of the accelerator with the local large scale turbulent magnetic field, neither of which are completely understood. Ultimately a unified theory that self-consistently captures the complete global picture bridging the two regimes is desired, but a clear understanding of the underlying assumptions is an essential first step.
Using test particle simulations in an isotropic Kolmogorov turbulent field, Giacinti et al. (2012) demonstrated that in the process of diffusing away from a source, before CRs have diffused over several correlation lengths Lc of the large scale field, the transport is highly filamentary. The filaments develop due to the anisotropic nature of particle diffusion in magnetised plasmas (e.g. Isenberg & Jokipii 1979). Giacalone (2017) has shown that similar filaments may also develop on smaller scales close to the shock, for a given field realisation. Efforts to incorporate this phenomenon into the standard picture have been applied, although for practical purposes in a reduced one-dimensional picture Malkov et al. 2013) while other have included the role of neutral particles in damping/stabilising the self-excited modes (Nava et al. 2016;D'Angelo et al. 2018;Brahimi et al. 2020). These works have exclusively worked within the confines of the standard quasi-linear theory framework discussed previously. Bell et al. (2013) on the other hand, made a first attempt to apply advances in our understanding of non-resonant excitation of magnetic fields into a more self-consistent picture of escape from SNRs. This same approach has also been applied by Blasi et al. (2015) to the escape of CRs from galaxies.
In this work, we revisit the linear analysis of unstable modes, developing a more generic framework for treating the linear perturbations, with the aim of providing fresh insight into the problem, while also laying the groundwork for future numerical investigations. The outline of the paper is as follows. In the next section we describe the general framework, including the key assumptions. Section 3 details application of the new results to ionised plasmas, using new results from Cas A to motivate the parameter choices. In section 4, the combined effects of neutrals and cosmic-rays are investigated. Due to the lengthy nature of the calculation, the derivation of the cosmic-ray response is left to the appendix.

COSMIC-RAY ANISOTROPY AND GENERALISED MHD INSTABILITY
We consider the generic case of transport in the neighbourhood of a strong CR source. A magneto-hydrodynamic description for the background plasma is adopted, which for simplicity we take to be composed of electrons and protons only. Combining the single fluid equations, while accounting for the presence of cosmic-rays and neutral hydrogen, the momentum equation for the charge carrying thermal fluid reads (Bell 2004;Reville et al. 2007) The last term (ncrqcrui × B) arises due to the fact that the return current is drawn relative to the background fluid, and it is this term that ultimately determines the neutral streaming velocity in the fully ionised case. νin is the momentum exchange rate due to elastic charge-exchange collisions with neutral hydrogen atoms. A derivation of this equation is provided in Appendix A, where previous kinetic approaches are also reviewed. Adopting ideal MHD, Faraday's equation reads The above equations must be solved together with momentum conservation for the neutral fluid where from momentum conservation, we have ρnνni = ρiνin.
The following expression for the ion-neutral collision frequency, valid in the range 10 2 K T < 10 6 K (Kulsrud & Cesarsky 1971) is used where Ti is the ion temperature in eV. For lower temperatures T ≈ 10 2 K, we use the results of Osterbrock (1961) which, since we neglect Helium contributions throughout this work, is simply We take this value to hold at lower temperatures, since the same physical arguments presented in Osterbrock (1961) apply. Note however, for such small temperature the weakly collisional regime may come in to play, which requires a different analysis. Bell et al. (2020) have recently demonstrated that a related instability, the magneto-collisional instability, occurs in this regime. It is also known that the environments of SNRs have significant amounts of dust, as do many parts of the ISM. A framework for adding this additional charged species to the system of equations and a detailed linear analysis has been performed by Squire et al. (2020). We will not consider this additional complication here.
We proceed with a standard linear analysis, looking at circularly polarised modes propagating parallel to a mean background field with average CR current along the mean magnetic field direction: j0 B0. Following Bell (2004), we make the substitution j ⊥ /j0 = σB ⊥ /B0 i.e. σ describes the response of the CR current to small transverse fluctuations in the magnetic field. Working in the local fluid rest frame ( ui = 0), and excluding any possible drift between neutrals and the ionised fluid to lowest order, the general dispersion relation is where vA = B0/ √ 4πρi is the Alfvén velocity. The cosmic-ray drift velocity vcr = j0/ncrqcr is introduced here, as well as an ionisation ratio χ = ρi/ρn = νni/νin. 1 For convenience, we have also separated out the ω dependent terms as σ(ω, k) = σ1(k) − (ω/vcrk)σ2(k) (see Appendix B).
A generalised expression for σ is given in the Appendix where an expansion of the cosmic-ray anisotropy to arbitrary order in Legendre polynomials is presented. Other basis functions (e.g. Chebyshev) are equally straightforward to derive, however, here we aim to develop a framework to compare to future simulations which will build on the work of  who used the related spherical harmonic functions as a basis. The penultimate expression for σ, which involves an integral over the magnitude of momentum, is given in a form that allows for straightforward numerical calculation which can be applied to an arbitrary distribution. In the present work, we focus on a fixed power-law distribution at all orders of the expansion f ∝ p −s Θ(p; p1, p2), where Θ(x; a, b) is a top hat function, equal to unity for a < x < b and zero otherwise. In this situation, the integrals, although cumbersome, can still be expressed in terms of standard functions, and can be shown to reduce to previous results in the appropriate limits. e.g. first order expansions with s = 4 were originally presented in Bell (2004), while thermal/composition effects, and alternative power-laws have been explored in Reville et al. (2007). We revisit and extend these results in the following.
Higher order anisotropies (to second order in Legendre polynomials) have been presented in Bykov et al. (2011), although there the emphasis was on the firehose instability. Using the asymptotic krg,2 1 values given in the appendix, where rg,2 = rg(p2) is the gyro-radius of the maximum energy cosmic-rays, independent of s the usual firehose expression is found where P ,⊥ are the parallel,perpendicular components (w.r.t. the guide field) of the cosmic-ray pressure tensor (Bykov et al. 2012;Scott et al. 2017). If the pressure anisotropy ∆Pcr = P − P ⊥ is large enough the system is firehose unstable. We will see that for what concerns the acceleration of the highest energy particles, this regime is generally too slow to be of any influence, at least for the cases we consider. We may also compare our result to that of the resonant ion-cyclotron instability, valid in the weak damping approximation. As we will show below, this case applies in the limit vcr/c UB/Ucr, where UB and Ucr are the energy densities of the magnetic field and CRs respectively. In this limit the first term on the rhs of equation 6 dominates and the real part of the dispersion relation corresponds to left or right propagating Alfvén waves Re(ω) = ±kvA. Using the asymptotic expressions given in equations (C8) and (C9) of the Appendix, we recover the familiar linear growth rate expression for the resonant streaming instability where Ω0 = eB0/mpc is the non-relativistic gyro-frequency, and ncr(> p) = 4π ∞ p f p 2 dp. Note that only waves with ω/k and vcr of the same sign have positive growth. We see here also the neutral streaming speed has the correct form, (e.g. Melrose & Wentzel 1970), namely that steady streaming is proportional to the Alfvén velocity. If the self-generated waves reduce the streaming velocity to the threshold condition vcr = svA/3, the next order term in the anisotropy can be included. Using the results in the appendix B5, the pressure anisotropy instability recently investigated by Zweibel (2020) is recovered. The instabilities are not mutually exclusive, but rather limiting forms of the resonantly driven instability for vanishing first and second order anisotropies respectively.
Interest in cosmic-ray streaming has been renewed in recent years, where significant efforts have been made to implement the process in large galaxy-scale simulations (Jiang & Oh 2018;Thomas & Pfrommer 2019;Chan et al. 2019). In this regard, it is necessary to highlight the range of validity of the streaming instability result above. The previously mentioned limit on applicability of the weak-damping approximation corresponds to a cosmic-ray energy density of Ucr 10 (vcr/vA) −1 n 1/2 gasBµG keV cm −3 . This limit has a simple physical interpretation, being equivalent to the condition that the jcr × B force on the background fluid does not dominate over the magnetic tension (a necessary condition for Alfvén wave propagation) on any scale in the resonant range. The above limit is not always guaranteed to be the case, particularly in the neighbourhoods of strong CR sources, where both Ucr and vcr/vA are enhanced relative to Galactic averages. When this condition is violated, the short-wavelength non-resonant branch first investigated by Bell (2004) dominates. However, as we show below, the resonant modes also become increasingly sensitive to the nature and details of the anisotropy.
With the goal of generalising the prescription for parallel propagating modes in magnetised plasmas containing a population of drifting CRs in mind, it will prove convenient to introduce the dimensionless variables (ω,ν) = (ω, νin) × rg1/vA andk = krg1, where rg1 is the gyro-radius of the lowest energy cosmic rays in the background mean field. We keep the treatment general, such that vcr can take any value, but will specialise later to the case of shock precursors.
In dimensionless form, equation (6) readŝ where we have introduced the dimensionless quantity In this notation, the condition for Bell-instability to operate is ζ1 > 1. If this condition is not satisfied, the magnetic tension dominates at krg,1 = 1, and the non-resonant branch is effectively stabilised. The weak-damping result given above generally applies in this case, provided the anisotropy is small. If the anisotropy is large, specifically if ∆Pcr/Pcr > vcr/c, the leading order terms are comparable. The effect of increased anisotropy will be shown in the next section. The second form for ζ1 given in (9) above is close 2 to that used in Bell (2004) although the assumption that vcr = u sh is no longer implicit. Note that setting UB to the total magnetic energy (turbulent + mean field) in this expression, the saturation condition given in Bell (2004) corresponds to ζ1 = 1. An alternative approach introduced by Bell et al. (2013), which is more applicable to CR escape upstream of a shock, is to define an escaping flux j pescc = ηesceρtotu 3 sh , where ηesc represents the fraction of the total mass energy density processes by the shock being carried away by escaping particles, andpesc the average momentum of these escaping particles. Setting this to p1, we find where MA,n = u sh /vA,n is the Mach number of the shock.
Here we introduce the reduced Alfvén velocity vA,n = B/ 4π(ρi + ρn) in anticipation of shocks that encounter an incompletely ionised medium in the far upstream. This reduces to the expression in Bell et al. (2013) in the limit ρn = 0. We explore the new generalised dispersion relation below for a selection of specific relevant cases.

STRONGLY IONISED PLASMASχ → ∞
Observations of several SNRs in Hα provide unambiguous evidence that SNR shocks interact with partially ionised plasmas. As previously pointed out by , the strong heating in the non-linear phase of CR-driven magnetic field amplification is likely to ionise most material before it reaches the shock, subject to the obvious condition that the driving is strong enough to reach this stage of nonlinear development. The apparent anti-correlation between synchrotron X-ray rims, a tracer of strong magnetic field amplification, and Hα emission in some cases (Winkler & Long 1997) indicate that this might be the case. Hence, the completely ionised scenario likely applies to the developed precursors of strong SNR shocks. Here we apply the new generalised expressions, motivated by recent observational results which better constrain the cosmic-ray parameters.

Cas A:
Recent observations by Ahnen et al. (2017) and in particular the combined Fermi-VERITAS analysis of Abeysekara et al. (2020) provide robust constraints on the CR parameters in Cas A. While a full numerical model that captures self-consistent acceleration and field growth is desirable, here we simply aim to use the revised dispersion relation to explore the various modes that mediate the scattering. We focus on the pure hadronic model of Abeysekara et al. (2020), which indicates a total CR energy ECR ≈ 1.7 × 10 50 erg, with best fit above the pion bump exhibiting a power-law with exponent s = 4.17 and exponential cut-off at pc = 17 TeV. Adopting for the total available energy of the SNR a value of ESNR = 2 − 5 × 10 51 erg (Chevalier & Oishi 2003), the hadronic model implies a total (i.e. intergrated over the SNR lifetime) CR acceleration efficiency ≈ 5%.
To determine the dispersion of the linear modes, we require two additional pieces of information. The cosmic-ray energy density in the precursor, and limits on the higherorder CR anisotropy. For the latter, we will restrict ourselves to second order, i.e. pressure anisotropy. For a shock radius of 2.5 pc, adopting the total CR energy in Cas A given by Abeysekara et al. (2020), the average energy density is Ucr ≈ 60 keV cm −3 . This is to be compared with the mean energy density of the shock processed material U sh ≡ 1 2 ρiu 2 sh ≈ 180 keV cm −3 . In the latter estimate, we have adopted the same gas density upstream of the shock ρi = 2.34 × 10 −24 g cm −3 as used in Abeysekara et al. (2020). In terms of hydrodynamic feedback, this gives η ≡ Pcr/ρ0u 2 sh ≈ 0.2. Comparable hydrodynamic efficiencies for other young remnants have been claimed previously in the literature, although here we will adopt it as an upper limit, as it can not yet be ruled out that most of the energy resides deep in the SNR interior.
With this picture in mind, we can repeat the estimates of Bell et al. (2013), assuming the non-resonant mode dominates, and that confinement demands a minimum of N ≈ 5 growth times (tNR = 2rg1/vAζ1). In a quasi-steady state, we expect the escaping flux to match the average accelerating flux close to the cut-off (Drury et al. 1999), where v · n is the projection of the particle velocity vector onto the shock normal, ∆p the average momentum change per cycle, and ∆u the velocity jump across the shock. Assuming the cosmic-rays populate a power-law distribution, extending without break to low energies, This integral can be expressed as a combination of Hypergeometric functions, but it suffices to know that it is of order unity for any realistic choice of parameters. Using this value for the escaping flux of CRs upstream, the confinement condition tNR < tSNR/N thus becomes pmax mpc where ωpi is the background ion plasma frequency upstream of the shock. For η = 0.2, and s = 4.17, this gives Emax ≈ 50 TeV, reasonably close to the observed cut-off, suggesting the efficiency of η ≈ 0.2 is justified. Regarding the pressure anisotropy, the standard diffusion approximation predicts ∆Pcr/Pcr = 2(u sh /c) 2 (Jokipii & Williams 1992). However, as shown by Luo & Melrose (2009), the non-resonant short wavelength instability of Bell (2004) leads to an enhancement of the anisotropy, and requires resonant scattering to counteract it. It is thus desirable to explore the relevant time-scales involved using the updated constraints in this object.
We first consider the case relevant to confined particles that occupy the precursor in the immediate upstream of the shock. As shown above, the acceleration efficiency, as inferred by the joint Fermi/VERITAS observations, is close to 20%, unless of course the CRs are mostly confined deep inside the SNR. We thus explore the possible parameter space, Figure 1. Real (left) and imaginary (right) components of the phase speeds for waves in a fully ionised background, normalised to the ambient Alfvén velocity, as a function of pressure anisotropy ∆ P = (∆Pcr/Pcr)/(2(u sh /c) 2 ) (i.e. the fractional pressure anisotropy relative to the diffusion approximation result). Dashed curves are for modes rotating in the same sense as the driving cosmic-rays, while solid lines in the opposite sense. We use values expected for Cas A: vcr = u sh = 5, 000 km s −1 , ρ i = 2.34 × 10 −24 g cm −3 and B 0 = 3 µG. We take p 2 /p 1 = 10 4.5 , such that for E min = 1GeV, Emax ≈ 30 TeV. We have taken a CR efficiency of Pcr/ρ i u 2 sh = 0.2 (see text). Note the cut-off in the real part of the frequency below kr g2 = 10 −4.5 due to onset of the firehose instability. η/η 0 = 10 0 η/η 0 = 10 −0.5 η/η 0 = 10 −1.0 η/η 0 = 10 −1.5 η/η 0 = 10 −2.0 η/η 0 = 10 −2.5 Figure 2. Real component of the phase speeds (left) and growth rates (right) in a fully ionised background and growth rate, as a function of cosmic-ray acceleration efficiency η = Pcr/ρ i u 2 sh . Pressure anisotropy is set to the diffusion approximation ∆ P = 2(u sh /c) 2 , and η 0 = 0.2 (see text). All other parameters are as in Fig.1. In the lower plot, for the conditions given, the relevant time scale for growth is r g1 /v A ≈ 10 6 (p 1 /mpc) s. We recall that the maximum growth rate is predicted to be Im(ωr g1 /v A = ζ 1 /2), where in these plots ζ 1 ≈ 1200η/η 0 keeping fixed the inferred spectral index of s = 4.17. In figures 1 and 2 we plot the phase velocity and growth rates for the full dispersion relation, truncating the expansion at second order. In figure 1 the CR efficiency is kept fixed at the upper limit inferred from observations, and we explore the effect of increasing anisotropy. The first thing to note is that the anisotropy has no effect on the maximum growth rate of the Bell instability This results from the non-resonant nature of the instability, for which only the total current matters. The CR current is kept fixed for all plots in Fig.1. Second, and contrary to the standard picture, only at extremely short wavelengths (or small ζ1) does the usual Alfvén branch emerge. At all other wave-lengths, the modes are highly super-Alfvénic and also dispersive. There is thus no unique wave-frame, and naively, one might anticipate this to have an effect on the transport in the precursor, and as a consequence the acceleration. The dispersion free (ω/k = constant) behaviour that is commonly assumed strictly speaking only applies to s = 4 power-law spectra, or in the ζ1 1 limit, where the waves are Alfvénic. The imaginary components are shown on the right, where we similarly scale by k to highlight the different growth rates for left and right circularly-polarised modes. The blue lines in Figures 1 and 2 show the same data. We do not explore the effect of wave dispersion further here, as the extremely short time-scales associated to the non-resonant waves is likely to dominate the field growth and subsequent particle transport (e.g. Bell et al. 2013).
In figure 2, we show the dependence of the phase velocity and growth rates on the cosmic-ray efficiency η relative to the reference value η0 = 0.2. This demonstrates the sensitivity of the growth rate, and phase speeds to the assumed cosmic-ray efficiency, as well as the effect of a radial decrease in ζ1. As expected, provided ζ1 > 1, scatterers in the resonant band are not adequately described by Alfvén waves, as is often assumed to be the case. Thus, in contrast to the standard model, the range over which the scattering fluctuations are non-Alfvénic may extend over a significantly larger volume than previously considered. Recently, several groups have performed kinetic simulations of the resonant ion-cyclotron instability in the weak-damping limit (e.g Holcomb & Spitkovsky 2019; Bai et al. 2019b). We will explore this process, and its impact on the surrounding environment, in a separate numerical study.
For Cas A parameters, the reference time rg1/vA is on the order of a few weeks, while the inverse lifetime of Cas A is rg1/vAtSNR ≈ 10 −4 . This coincidentally is close to the growth rate for waves resonant with particles close to the inferred cut-off in Cas A, provided acceleration remains efficient. Long wavelength, non-resonant firehose unstable fluctuations however appear to be too slow to be of significance, although an enhancement by almost an order of magnitude could be achieved if the anisotropy is greatly enhanced relative to its diffusive limit (see Fig 1). The existence of such strong anisotropy would need to be properly motivated.

COSMIC-RAY PROPAGATION IN PARTIALLY/WEAKLY IONISED PLASMAS
The linear theory of parallel propagating plane-waves in incompletely ionised magnetised plasmas is well known (Kulsrud & Pearce 1969;Zweibel & Shull 1982;Tagger et al. 1995;Xu et al. 2016), as is its effect on diffusive shock acceleration (Drury et al. 1996). Interest around this topic has grown in recent years, primarily in the context of escape from SNRs (Bykov et al. 2013;Nava et al. 2019;Inoue 2019), but also regarding the penetration of CRs into molecular clouds including their role as an ionising agent, as well as their transport in the multi-phase ISM (Padovani et al. 2009;Phan et al. 2018;Brahimi et al. 2020). Complexities associated to the field topology and its impact on CR observable such as the Boron to Carbon ratio have also been explored (D'Angelo et al. 2018). The analysis here follows the approach of Reville et al. (2007). The standard QLT approach of applying the diffusion approximation of Skilling (1971): ∂f /∂µ ≈ λ mfp ∂f0/∂x (x here is distance along the mean field) only holds if gradients are long relative to the mean free path λ mfp . This may not be satisfied for dynamically evolving systems, and hence we keep ζ1 as a free a parameter. Our aim here is partly to bring attention to areas where more numerical efforts are required to address this problem self-consistently.
Incompletely (χ ∼ 1) and weakly ionised gases without cosmic rays For incompletely ionised (χ > 1) gases, the linear analysis is well known, and we refer the reader to Zweibel & Shull (1982) for a clear presentation. The key limiting results can be read off Fig. 3. At high frequencies, ω νin, inter-species collisions do not much affect the oscillations. This results in slowly damped Alfvén waves. At low frequencies ω νni, the ions and neutrals are collisionally coupled, which results in the two species oscillating as a single fluid with increased inertia. This has the effect of reducing the Alfvén velocity ω/k = vA,n, and since the coupling is more effective at lower frequencies, the damping naturally decreases rapidly (∝ k 2 ) as one moves to longer wavelengths. The upper blue line in the right hand panel corresponds to damping rates of nonpropagating perturbations in the neutral fluid, should they exist.
The case of weakly ionised gases (χ 1) reveals a number of additional complexities, and is of course more relevant to CRs impacting on dense molecular material. In the absence of CRs, and subject ot the condition χ < 1/8, waves in the so-called evanescent range 2νni/vA,n < k < νin/2vA do not propagate (Zweibel & Shull 1982). For sharp resonance (k = 1/rg) this equates to cosmic-ray energies in the approximate range with reference values given in subscripts. For many relevant parameters, this falls in an energy range pertinent for GeV-TeV γ-ray observations. Tagger et al. (1995) have pointed out that the wave propagation can be viewed in two ways (see also Soler et al. 2013). In the generally standard approach of taking coherent plane waves with real k, no solutions for the real part of ω are found in this wave band. However, taking the opposite approach by considering a wavepacket at fixed real ω and examine how a wave grows/decays as it propagates, Tagger et al. (1995) show that waves in this band do in fact propagate, damping as they do so. This latter approach is also well-suited to shock precursor studies, where wave amplitudes are expected to develop in-homogeneously as they approach the shock. A comparison of the two approaches is shown in Fig. 3, where the lower panels are for an ionisation factor χ = 10 −2 . As expected, the two methods agree at high and low frequencies, but perturbations excited with Fourier components in the evanescent band can clearly be seen to propagate only in the latter approach. Implications are discussed in section 5 Weakly ionised gases with cosmic rays Here, we focus exclusively on gases in which χ < 1/8, i.e. the threshold condition for the evanescent band to exist (Zweibel & Shull 1982). The inclusion of CRs into the calculations should be considered with some care. It is convenient to continue working with the cosmic-ray and magnetic energy 10 −2 10 −1 10 0 10 −2 10 −1 10 0 10 1 10 2  . Real (ωr) and imaginary (ω i ) components of the frequency showing phase velocity and damping rates in the absence of cosmic rays for χ = 1 (top) and χ = 10 −2 (bottom). Both plots assume a fixed ntot = n i + nn = 10 cm −3 , a temperature of 1 eV, a magnetic field of 5 µG, while r g1 corresponds to that of a 100 GeV proton in this field. In these units, we find for χ = 10 −2 , the evanescent band occurs in the range 0.3 < kr g1 < 0.8, as shown.
densities, as both can in principle be inferred from observations. However, one should make sure to check in extremely weakly ionised cases, that the MHD approximation holds, i.e. ncr ni. The short-wavelength ζ1 > 1 case was previously considered by Reville et al. (2008). The present results extend those results to arbitrary wavelengths, and upper cut-offs in the distribution. As discussed in the appendix, the contribution from the end points of the distribution must be correctly accounted for to ensure consistent results.
The damping of scatterers due to ion-neutral collisions in and near CR accelerators has been explored previously in the context of spectral breaks (Malkov et al. 2011), cosmic ray escape from SNR (e.g. Xu et al. 2016;Brahimi et al. 2020), and either the penetration through (Cesarsky & Volk 1978;Morlino & Gabici 2015; Inoue 2019) or exclusion from (e.g. Skilling & Strong 1976) molecular clouds. Such studies are particularly interesting when they occur in the neigh-bourhood of a CR source. As discussed above, the existence of an evanescent band is an artefact of the chosen approach, and may not exist in a turbulent plasma, provided the MHD cascade continuously replenishes the wave-band in question. Here we show that even a minute presence of CR anisotropy also causes waves in the evanescent band to propagate.
In Figure 4 we consider gas conditions representative of the warm and cold neutral medium phases (WNM and CNM respectively) of the ISM, and the effect of varying CR anisotropy. The details can be read from the caption. Note that we adopted a modest value of Ucr/UB = 10 in these plots, and explored the effect of increased anisotropy. The dashed curves correspond to Ucr/UB = 100 to show an example of ζ1 > 1 limit. This again shows how waves deviate from standard Alfvén branch in the strongly driven regime. While not shown here, it is easily shown that provided ζ1 1, the imaginary part of the frequency matches the expectation from a linear addition of the damping and For the WNM we adopt the parameters n i + nn = 0.5 cm −3 , T = 5, 000K and χ −1 = 50, while for the CNM we take n i + nn = 50 cm −3 , T = 50K and χ = 10 −3 . In both cases the magnetic field is set to 5 µG. Note that for the two dashed lines ζ 1 > 1.
growth rates (with an appropriate choice for the Alfvén velocity at long wavelengths for the latter). This also turns out to be a reasonable approximation to the growth/decay in the ζ1 > 1 limit, i.e. Im(ω) ≈ Im(ω)no CRs + G(k)Im(ω) no neutrals , where and with a linear fit in between. This renormalizing factor is necessary to account for the increased inertia of the plasma. The agreement is typically good to within a few percent although slightly larger variations may be found close to the evanescent band that occurs in the absence of CRs. It is well known that for a fully ionised plasma, the maximum growth rate for field-aligned CR driven modes is given by the expression in equation (14), regardless of the value of ζ1 (e.g. Bell 2013). However, there is a fundamental difference between the physical processes underlying the instability in the two limits of ζ1 greater than or less than unity. The presence of neutrals can expose this difference.
It is easily verified that the magnitude of the ion-neutral damping rate is always less than the oscillation frequency kvAG. If ζ1 > 1, the maximum growth rate is also kmaxvAG, where kmax = ζ1/2rg1 is the wavenumber of the fastest growing mode. This is merely confirmation of the proof in Reville et al. (2007), that the non-resonant instability can not be stabilised by neutral damping. The opposite is true for resonant modes. If ζ1 < 1, and assuming the distribution is peaked at an energy corresponding to gyroradius rg1, by definition kmax = 1/rg1 and γmax < GkvA and hence can be stabilised.
Two possibilities exist. If rg1 < ζ1vA/νin (c.f. the lower value in the inequality (15)) growth is still possible, however for small ζ1 this may correspond to very low energies. On the other hand, if rg1 > ζ1vA/νin net damping is expected at short wavelength, but since the damping rate at long wavelengths scales as k 2 , and resonantly driven modes scale as k s−3 (neglecting pressure anisotropy effects), higher energy particles, if present, can still produce net growth. Thus it is possible to define a critical energy, below which damping must dominate.
Taking for simplicity s = 4, we see this equates to an energy ζ −1 1 times the upper energy limit in the inequality (15), with all energies below this being unable to self-generate scattering waves. While it is tempting to speculate on the impact of this effect for escape from SNR, or indeed penetration into a molecular cloud, we will defer such discussion until reliable numerical simulations can explore the full multi-scale aspects.

DISCUSSION
The stability of parallel propagating transverse MHD fluctuations in the presence of an external current of cosmic rays has been explored. The solution given in equation (6), and in particular the method for determining the cosmic-ray response provided in the appendix, generalises many of the existing results in the literature. The approach can readily be applied to other plasma conditions, where parallel propagating modes are of interest. A limitation of the analytic results presented here, is the use of a single power-law dependence at all orders of anisotropy. For example, in the context of shock precursors, one might anticipate the anisotropy at large momenta to have a different spatial profile to that at low momenta, which might for example enhance the growth rate of the fire-hose instability. A follow-up study of cosmic-ray anisotropy in multi-scale self-generated fields will exploit the generic form presented here.
The standard quasi-linear theory approach to studying cosmic-ray transport and self-confinement applies in the limit of ζ1 < 1, or Ucr/UB c/vcr. In the Galactic disk, we expect Ucr ≈ UB and vcr ≈ vA c, meaning this condition is generally satisfied. However, near to sources, where both Ucr and vcr are greatly enhanced, this is generally not the case. Similar conclusions may apply for the low density environment at the periphery of the Galactic halo (e.g. Blasi & Amato 2019). We plan to bridge the two regimes in an upcoming numerical study.
The transport of cosmic-rays in the ISM, and the related topic of their escape from sources, has recently received significant interest, with particular emphasis on the different phases of the ISM. The role of ion-neutral damping features heavily in these studies. That different conceptual approaches to the problem produce opposite conclusions regarding propagation of Alfvénic fluctuations at intermediate frequencies in weakly ionised plasmas (νni < ω < νin) is already known (Tagger et al. 1995;Soler et al. 2013). However, studies to date have focussed exclusively on only one such approach, which, motivated by the new results, we now argue is in fact less likely to apply.
Turbulent MHD cascades can be conceptualised as interactions between colliding wave-packets (Kraichnan 1965). To mediate efficient exchange of energy between Fourier modes, the semi-classical conservation laws for three-wave couplings: k1 + k2 = k3 and ω(k1) + ω(k2) = ω(k3) should be satisfied. Here ω(k) is the 3D generalisation of the dispersion relation in this paper 3 . It is well-known that for incompressible MHD, coupling requires either vanishing k1 · vA or k2 · vA, i.e. one of the Fourier modes must be nonpropagating and in a plane perpendicular to the guide field. This generally results in an anisotropic turbulence spectrum (see Zhou et al. 2004, for a review). If the waves are however dispersive, which will always be the case in the neighbourhood of an hypothetical evanescent band, it is straightforward to show that this is no longer a requirement, provided ω = 0. Thus, based on the dispersion curves shown in section 4, there is no physical justification to exclude waves with k3 falling within the evanescent band, and allowing them to propagate with phase velocity between the two Alfvén velocity limits, as shown in Fig 3. The existence of a region of non-propagating fluctuations would appear to be a highly idealised concept. The dispersive nature of the waves however may still cause modifications to the resulting turbulence spectrum and may yet leave an imprint on energetic particle transport. To our knowledge, there has been no thorough investigation of this effect, which may reveal itself in the γ-ray emission from CRs in molecular clouds.
To this end, recent studies of the diffuse γ-ray emission from giant molecular clouds with Fermi-LAT (e.g. Aharonian et al. 2020) hint at a dependence of the Galactic cosmicray density with Galacto-centric radius, although could also correlate with nearby sources. The lack of any clear trend in the spectra however hint at the absence of unique energy-dependent behaviour/trends. Surveys that look at higher energies with instruments such as HAWC (Lauer 2015), LHAASO (Bai et al. 2019a) and the anticipated southern hemisphere counterpart SWGO (Albert et al. 2019), may provide a different picture. The analysis presented in this work is important for the transport at such energies, and future observations in the VHE to UHE γ-ray regimes may reveal new insight into the role of cosmic-ray feedback on molecular clouds. Unfortunately, we currently lack predictive capability.
Finally, we note that attention was restricted in the present study to relativistic particles. The results could however be extended to non-relativistic particles without difficulty. Since low energy cosmic rays may play an important role in the ionisation of dense molecular clouds, the ability of self-induced scattering to inhibit penetration into the cores of such clouds warrants further investigation.

APPENDIX A: DISPERSION RELATION -KINETIC VS FLUID APPROACH
For completeness, we compare different approaches to deriving the linear dispersion relation for circularly-polarized modes propagating parallel to a background field (B0 = B0x). Following standard text-books, we explore first order small perturbations as a spectrum of periodic transverse waves (e.g. Ichimaru 1973;Krall & Trivelpiece 1973). For a uniform background, the zeroth order total current and charge density must vanish. Taking all first order quantities to be of the form ξ(x, t) = ξ k e i(kx−ωt) , with k · δE k = 0 Maxwell's equations can be combined to give where is the perturbed current associated to species s, with qs its charge. The summation is taken over all species in the plasma.
To derive a dispersion relation, it is necessary to describe the response of the current to a fluctuating magnetic field. This is determined from the Vlasov equation, using the zero-th order helical trajectories along the mean field(see for example Krall & Trivelpiece 1973, chapter 8), and gives Here fs(p , p 2 ⊥ ) is the gyrotropic zeroth order equilibrium solution of sepcies s. The i0 + term ensures the integration is taken along the Landau contour, and ωg,s = qsB0/γmsc is the relativistic gyro frequency. Introducing the complex quantities j ⊥ = jy + i jz (and similarly for δE), where = ±1 determines the polarisation, the current takes the rather simple form This expression is found in most textbooks, however, it proves useful to write it out explicitly here, as we will use it in the next section. Applying the same notation to Eq (A1) and substituting in the expression for the current, one recovers the familiar expression where the plasma susceptibility has been introduced. This is the starting point for kinetic investigations of circularly polarised linear modes that propagate parallel to a guide field. As is standard, an equilibrium condition must be considered that satisfies the general constraints of zero net charge and current. When a current is provided by streaming cosmic rays, a compensating return current must be drawn by the background plasma. Since the cosmic rays are generally small by way of number, one approach is to simply assume the thermal electrons establish a slow drift relative to the background ions that provides the required return current (e.g. Achterberg 1983;Reville et al. 2006;Luo & Melrose 2009). An alternative scenario involves only a sub-population of electrons are drawn to provide the return current (Amato & Blasi 2009;Zweibel & Everett 2010). We note however, that such a configuration is also two-stream unstable, and for realistic scenarios, is likely to relax rapidly to the former. In the work of Amato & Blasi (2009) the two approaches are shown to ultimately give equivalent results, at least in the low frequency regime of interest. In this regime, it is of course expected that the results approach the MHD limit, which is insensitive to the microphysics of the thermal gas.
In the fluid approach, one need only consider the bulk properties of all quantities, although a general treatment still requires a kinetic treatment of the CRs. The fluid approach has the tremendous advantage that it allows easy inclusion of collisions, something that would be undesirably cumbersome in a kinetic approach.
Starting from the fluid equations for the thermal species, which for simplicity we limit to be electrons, ions and neutrals, for each species we have where νsα is the interspecies momentum exchange rate, and we have assumed isotropic pressure for each species. In the following, we consider only collisions between the ions and neutrals.
We proceed by summing over the charged components in the usual fashion, to give where ρ = nimi + neme ≈ ρi and u = (nimiui + nemeue)/ρ ≈ ui . In the absence of CRs, the last two terms equate to 0 and jtot respectively. However, as in the kinetic description above, we must not neglect the CR contribution. Thus, from quasineutrality and Ampere's law, we have s=i,e nsqs = −ncrqcr and c 4π Subsituting into the above, we find which in the ideal MHD limit, E = −u × B/c reduces to Equation (1). The kinetic elements of this problem are introduced via jcr, which we describe in the next section.

APPENDIX B: THE COSMIC-RAY MAGNETIC RESPONSE
Following Bell (2004), we wish to express the response of the cosmic-ray current to fluctuations in the magnetic field. As discussed and summarised above, this equivalent procedure for electric fields can be found in most plasma text books. However, we can easily connect the two methods. Using the complex notation from the previous section, it follows from Faradays' law, which on substitution into (A4) gives Note we have dropped the subscript for species, as we are only concerned with CRs. We seek a general formalism to determine the response for a given Legendre polynomial expansion of the distribution function. In this case it is convenient to first transform to spherical momentum coordinates where we have introduced µ0 = ω kv − ωg kv , with k > 0. Expansion of the distribution in a Legendre polynomial series f (p, µ) = f (p)P (µ), leads to where from Rodrigues' formula one can show with generalised binomial coefficients a n = 2 n +n−1 2 .
The integration over pitch angle is contained in the term Recall that j ⊥ /j0 = σB ⊥ /B0 and therefore This format is convenient when the series is to be truncated after a finite number of terms, i.e. setting f >Lmax = 0. Consider for example Lmax = 3, the first few An terms are For the integral over µ in In, we note that if |µ0| < 1, the integral has a simple pole at µ = µ0, and hence in general we may write with P denoting the Cauchy principal value. The In satisfy the recurrence relation In+1 = µ0In + 3 2 1 + cos nπ (n + 1)(n + 3) and hence in computing the pitch angle integration it is necessary to specify only I0: This formalism is still completely general (we have not made any assumptions about ω in high/low frequency limits, etc.). It is particularly convenient for numerical calculations, as the final integral over momentum can be performed using standard methods. The approach may thus be applied to other scenarios where field aligned modes are of interest. A similar result is found using Chebyshev polynomials for the angular basis functions, but since one of the motivations is to provide a general framework for comparison with future simulations using an expansion in the related Spherical Harmonic basis functions , for brevity, we leave this extension to the reader.
Before proceeding to derive expressions for power-law distributions of relativistic particles, we make a brief digression to demonstrate how the above technique could be applied to other cases, such as a slow drifting Maxwellian distributions. To compare with previous results, we use the previous expression for the susceptibility, Equation (A6), which in spherical coordinates reads where ωp = (4πnq 2 /m) 1/2 is the plasma frequency, and Z = µ0|v=v th . This agrees with the standard result for the long wavelength limit of the plasma dispersion function (e.g. Ichimaru 1973). While the approach is not as straightforward, the algorithm is simple and trivial to implement numerically. Most importantly, it may prove useful for situations in which the moments of the distribution function do not reduce to standard integrals.

APPENDIX C: SOLVING FOR POWER-LAW COSMIC RAY DISTRIBUTIONS
To simplify the analysis that follows (although not strictly necessary at this stage), we consider only modes with ω ωg (which can be checked after the fact) making µ0 = − ωg/kv ≡ − λ. Here we have defined the normalised wavelength λ = ωg/kv = (krg) −1 . In this limit We seek a closed form solution to develop astrophysical applications. To this end, we follow Bell (2004), taking a power-law distribution f (p) = φ p −s Θ(p; p1, p2) where Θ(p; p1, p2) = H(p − p1) − H(p − p2), and H(p) = p −∞ δ(x)dx is the Heaviside step function. The assumption that all orders in the expansion have the same power-law shape is an obvious short-coming, as in non-equilibrium systems the anisotropy may well have a non-trivial momentum dependence, but we defer any detailed discussion to future numerical experiments, and a full non-linear theory is beyond the scope of this paper.
In the following, we consider only situations where p1 mc such that we can approximate v/c = 1. We note that In their final results, Bell (2004) and Reville et al. (2007) ignored the contribution from the end point delta-functions, as they play no role in the short wavelength limit. However, they are essential in the long wavelength limit, and neglecting them would result in the unphysical result of lim k→0 σ = 1, i.e. cosmic rays would not follow the field lines on large scales! Substituting into An, and changing variable to λ = λ1p1/p = λ2p2/p, we find Lmax =n a n φ λ 2 1 δ(λ − λ1) − λ 2 2 δ(λ − λ2) − (n + s)λ Θ(λ; λ2, λ1) = λ s−1 (p1λ1) s Θ(λ; λ2, λ1)Ãn Likewise, changing variable in the momentum integral it follows that where we have introduced the bulk CR velocity vcr = jcr/qncr. From here on, we restrict our attention to positive k (and λ) and define the integral which also satisfies a useful recursion relation: Using this allows us to write Finally, for convenience in the main part of the text, we will split up the linear susceptibility into frequency dependent/independent parts, Lmax =n a n φ (n + s)K s−3 n − λ s−2 In .

C1 Evaluating K m n
From the definition of In it follows that K m n contains only terms which are either integrals of powers of λ or the logarithmic function via: where m ≥ 0. The latter has solution for general m where Gauss' Hypergeometric function appears via We note the limiting forms of Fm are where the long wavelength limit is only accurate for m > 1. This completes our description of the necessary equations to generate the full solution for an arbitrary order expansion. We conclude by writing out the first few terms: Below we consider some specific integer examples which have special regular solutions.

C3 Moments of the distribution function
To help visualise the above results, normalised by the unhelpful φ terms, it is convenient to introduce more familiar fluid concepts. Using the ortho-normality of the Legendre polynomials, one can easily recover the cosmic-ray density, current and pressure tensor, being respectively: ncr = d 3 p f = 4π dp p 2 f0 jcr = q d 3 p vxf = 4π 3 q dp p 2 vf1 T µν cr = d 3 p µ v ν f = 4π 3 dp p 3 v f0δ µν + f2 5 diag (2, −1, −1) From the latter, we identify the isotropic and anisotropic pressures Pcr = 4π 3 dp p 3 vf0 , ∆Pcr ≡ T 11 cr − T 22 cr = 4π 5 dp p 3 vf2 .
For power-law distributions as we have used in the previous section, we can re-express the first few φn in terms of these fluid parameters: C4 Susceptibilities -example to second order in f and limiting cases Using these expressions, we provide the required scalar susceptibilities, i.e. the response of the CR current to linear magnetic perturbations. The relative magnitudes of the different moments will depend on the physical scenario under investigation. For example, in the diffusion approximation, it is generally found that φ +1 ≈ (vcr/c)φ (e.g. Jokipii & Williams 1992).
As we have only considered the first three fluid moments, we consider only the terms in the Legendre expansion up to f2, Im(σ 1 ) Im(σ 1 ) Figure C1. σ 1 plotted for different power-law indices s. In all cases p 2 /p 1 = 10 4 . For the real part, 1 − Reσ is plotted, to emphasise how it approaches zero as k → 0, as it must. where expressions for K m n are found in equations (C3) above. The real and imaginary parts of σ1 are plotted in Figure C1, where we have taken ∆Pcr > 0 in all cases. A crucial test of the solution is that the long wavelength limit (krg2 1), the cosmic rays must follow the magnetic fields, i.e. j ⊥ /j = B ⊥ /B , or equvalently, σ → 1. We observe that this is indeed the case. σ2 is dominated by the leading two terms in equation C5, and plotting for different parameters is un-necessary, although graphically σ2 look very similar to the left column in Figure C1. This is also to be expected, since it is this term which stabilises the resonant streaming instability when ω/vcrk approaches unity, which in turn determines the neutral streaming speed of the cosmic rays. Strong anisotropy of course changes this (see for example Zweibel 2020).