Cosmic ray CR feedback plays a vital role in shaping the formation and evolution of galaxies through their interaction with magnetohydrodynamic waves. In the CR selfconfinement scenario the waves are generated by the CR gyroresonant instabilities via CR streaming or CR pressure anisotropy and saturate by balancing wave damping. The resulting effective particle scattering rate by the waves nueff critically sets the coupling between the CRs and background gas but the efficiency of CR feedback is yet poorly constrained. We employ 1D kinetic simulations under the MagnetohydrodynamicParticleInCell MHDPIC framework with the adaptive deltaf method to quantify nueff for the saturated state of the CR pressure anisotropy instability CRPAI with ionneutral friction. We drive CR pressure anisotropy by expandingcompressing box mimicking background evolution of magnetic field strength and the CR pressure anisotropy eventually reaches a quasisteady state by balancing quasilinear diffusion. At the saturated state we measure nueff and the CR pressure anisotropy level establishing a calibrated scaling relation with environmental parameters. The scaling relation is consistent with quasilinear theory and can be incorporated to CR fluid models in either the singlefluid or pbyp treatments. Our results serve as a basis towards accurately calibrating the subgrid physics in macroscopic studies of CR feedback and transport.
['Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Astrophysics of Galaxies']
Kinetic simulations of the cosmic ray pressure anisotropy instability cosmic ray scattering rate in the saturated state
{'No Header': 'Draft version September 16, 2024 Typeset using L A T E X twocolumn style in AASTeX631', 'Kinetic simulations of the cosmic ray pressure anisotropy instability: cosmic ray scattering rate in the saturated state': 'Xiaochen Sun, 1, 2 Xue-Ning Bai , 1, 3 and Xihui Zhao 1  \n<!-- image -->  \n1 Institute for Advanced Study, Tsinghua University, Beijing 100084, China 2 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA 3 Department of Astronomy, Tsinghua University, Beijing 100084, China', 'ABSTRACT': 'Cosmic ray (CR) feedback plays a vital role in shaping the formation and evolution of galaxies through their interaction with magnetohydrodynamic waves. In the CR self-confinement scenario, the waves are generated by the CR gyro-resonant instabilities via CR streaming or CR pressure anisotropy, and saturate by balancing wave damping. The resulting effective particle scattering rate by the waves, ν eff , critically sets the coupling between the CRs and background gas, but the efficiency of CR feedback is yet poorly constrained. We employ 1D kinetic simulations under the Magnetohydrodynamic-ParticleIn-Cell (MHD-PIC) framework with the adaptive δf method to quantify ν eff for the saturated state of the CR pressure anisotropy instability (CRPAI) with ion-neutral friction. We drive CR pressure anisotropy by expanding/compressing box, mimicking background evolution of magnetic field strength, and the CR pressure anisotropy eventually reaches a quasi-steady state by balancing quasi-linear diffusion. At the saturated state, we measure ν eff and the CR pressure anisotropy level, establishing a calibrated scaling relation with environmental parameters. The scaling relation is consistent with quasilinear theory and can be incorporated to CR fluid models, in either the single-fluid or p -byp treatments. Our results serve as a basis towards accurately calibrating the subgrid physics in macroscopic studies of CR feedback and transport.  \nKeywords: Plasma astrophysics (1261) - Alfven waves (23) - Magnetohydrodynamics (1964) - Cosmic rays(329)', '1. INTRODUCTION': "Cosmic rays (CRs) are (trans-)relativistic charged particles pervading in space. With relatively low number density ( ∼ 10 -9 cm -3 ), they possess an energy density ( ≥ 1eV cm -3 ) comparable to the internal, kinetic and magnetic energy density of gas in the Galaxy (see reviews, e.g. Blasi 2013; Zweibel 2013; Grenier et al. 2015). Consequently, CRs are expected to play a dynamically important role, known as CR feedback (see reviews, e.g., Ptuskin 1997; Ferri'ere 2001; Zweibel 2017; Naab & Ostriker 2017). At a more fundamental  \nCorresponding author: Xiaochen Sun, Xue-Ning Bai  \nsun.xiaochen@princeton.edu  \nxbai@tsinghua.edu.cn  \nlevel, the interaction between the CRs and the (ionized) gas is primarily through fluctuations (waves) in background electromagnetic field (e.g., Jokipii 1966; Voelk 1975), as opposed to direct collisions (e.g., Radin et al. 1974; Ferrando et al. 1988; Schlickeiser 1989). Based on the source of such waves/fluctuations, CR feedback/transport are classified into two categories: fluctuations generated by CRs themselves, known as 'selfconfinement' (e.g., Wentzel 1969; Skilling 1971; Holman et al. 1979), and fluctuations cascaded from extrinsic turbulence (e.g., Kolmogorov 1941; Iroshnikov 1964; Kraichnan 1965; Goldreich & Sridhar 1995). As CR scattering and transport are primarily sensitive to electromagnetic fluctuations at their gyration scale, it is expected that self-confinement governs low-energy CRs (from sub-GeV to ∼ 10 2 GeV, Blasi et al. 2012; Evoli  \net al. 2018), while extrinsic turbulence likely dominates the transport of higher-energy CRs.  \nAs the CR energy density typically peaks around ∼ GeV (e.g., Adriani et al. 2011, 2013; Consolandi 2014), self-confinement is expected to be the dominant process driving the force of CR feedback. At microphysical level, this is mediated by the CR gyro-resonant instabilities (e.g., Lerche 1967; Wentzel 1968; Kulsrud & Pearce 1969; Ginzburg et al. 1973; Berezinskii et al. 1990). Here, magnetohydrodynamic (MHD) waves in the background gas can be destabilized when the CRs become weakly anisotropic with respect to the bulk gas (by more than about U A /c , where U A is the Alfv'en speed) through gyro-resonant interactions between the CRs and the waves. The classic flavor is the CR streaming instability (CRSI, e.g. Wentzel 1974; Skilling 1975a,b,c), which occurs when the speed of CR bulk motion relative to background gas exceeds U A . This situation naturally arises when the CRs escape from the source. With the gyro-resonant interactions, the growth of MHD waves are fed by the free energy from CR streaming, which in turn reduce the bulk CR motion towards the Alfv'en wave speed (Adkins & Schekochihin 2018). Another flavor of the CR gyro-resonant instability, which is poorly studied, is the CR pressure anisotropy instability (CRPAI, e.g. Lazarian & Beresnyak 2006; Yan & Lazarian 2011; Zweibel 2020), which arises when the level of CR pressure anisotropy exceeds ∼ U A /c . This situation naturally occurs when the bulk gas undergoes compression/expansion and/or shear, where CR anisotropy is driven by a change in background field strength ˙ B/B under the conservation of magnetic moment. The interaction between CRs and electromagnetic fluctuations tends to isotropize the CR momentum distribution to within the ∼ U A /c level in the gas co-moving frame.  \nThe coupling between the CRs and the background gas is reflected in the effective scattering rate ν eff , which is determined by the amplitude of waves. In the selfconfinement regime, the wave amplitudes are set the balance between driving (from the CR gyro-resonant instabilities) and various wave damping mechanisms, which eventually determines the efficiency of CR feedback. These include the ion-neutral collisions (Kulsrud & Pearce 1969; Soler et al. 2016), the non-linear Landau damping (Lee & Volk 1973), the linear Landau damping (Foote & Kulsrud 1979; Wiener et al. 2018), the extrinsic turbulence (Farmer & Goldreich 2004; Lazarian 2016; Xu & Lazarian 2022), and the influence of charged dust grains (Squire et al. 2021). Which of the mechanisms would dominate depends on the local environment, with ion-neutral damping dominating in the neutral medium, while non-linear Landau damping is more prevalent in  \ndilute gas (Zweibel 2013; Armillotta et al. 2021, 2022, 2024).  \nAt macroscopic level, galaxy simulations often model CRs as a fluid and employ approximate prescriptions for the effective scattering rate. Modeled as fluids instead of individual kinetic particles, the CR population is characterized by its energy density, momentum flux, and pressure (McKenzie & Voelk 1982; Guo & Oh 2008; Pfrommer et al. 2017; Jiang & Oh 2018; Thomas & Pfrommer 2019; Chan et al. 2019). Early simulations at galactic scales typically adopt highly simplified treatment (e.g., constant diffusion or streaming) of CR scattering, generally finding significant dynamical consequences especially in heating and in driving large-scale galactic outflows (e.g., Guo & Oh 2008; Dorfi & Breitschwerdt 2012; Uhlig et al. 2012; Hanasz et al. 2013; Girichidis et al. 2016; Wiener et al. 2017b; Jacob et al. 2018; Butsky & Quinn 2018; Chan et al. 2019; Dubois et al. 2019; Dashyan & Dubois 2020), but the results are highly sensitive to the CR prescriptions. More recent works start to incorporate more realistic prescriptions of the CR scattering rates based on the CRSI and various damping mechanisms (Armillotta et al. 2021, 2024), emphasizing the importance of resolving the multi-phase medium, while current state-of-the-art models are still unable to reproduce all observational constraints from the Galaxy (Hopkins et al. 2021).  \nThe main motivation of this work is two-fold. First, we note that the scattering rates in the self-confinement regime, as calculated from balance between wave growth and damping, are usually obtained under the quasilinear theory (QLT) as in the aforementioned simulations. However, it involves a number of approximations that likely yield results that deviate from the reality, and kinetic simulations are needed to clarify and calibrate these results (e.g. Bai 2022). Second, nearly all calculations so far focus on the CRSI (e.g. Holcomb & Spitkovsky 2019; Bai et al. 2019; Haggerty et al. 2019; Lemmerz et al. 2024), which indeed likely dominates CR scattering and heating in the general situations (Zweibel 2020), but the role CRPAI have not been systematically investigated. In this work, we conduct kinetic simulations of the CRPAI, and by balancing wave growth and damping, we aim to calibrate the CR scattering rates from the CRPAI from first principles.  \nOur simulations employ the magnetohydrodynamicparticle-in-cell (MHD-PIC) method (Zachary & Cohen 1986; Lucek & Bell 2000; Bai et al. 2015), which treats CRs kinetically as particles while treat the background plasma as gas described by MHD. Compared to the conventional PIC methods, it avoids resolving the microscopic scales of the background plasma, thus becomes  \nhighly advantageous in simulating the CR gyro-resonant instabilities. Using the ATHENA MHD code (Stone et al. 2008; Bai et al. 2015), Bai et al. (2019) developed the framework for simulating the CRSI, which captures the linear growth over a wide range of wavelengths and quasi-linear evolution over a broad range of CR energies around the peak of the energy distribution (which can be taken to be ∼ GeV). This is facilitated by the implementation of the δf method (e.g. Parker & Lee 1993; Kunz et al. 2014) to suppress Poisson noise, allowing one to use a reasonable number of particles to properly sample a weakly anisotropic particle distribution. Subsequent works of Plotnikov et al. (2021); Bambic et al. (2021) further incorporated a simple prescription of ion-neutral damping. In Bai (2022), we introduced the streaming box framework, where CRs move along an imposed CR gradient, allowing the CRSI to be continuously driven, balanced by damping. This framework thus enables the measurement of the CR scattering rates at the saturated states from first principles.  \nIn this work, we migrate the MHD-PIC simulation to the ATHENA++ MHD code (Stone et al. 2020; Sun & Bai 2023), which is more flexible and computationally efficient. We have already tested that the code well captures the linear growth of the CRPAI. In preparation for this work, we have also developed the expanding/compressing box (Grappin et al. 1993; Sironi & Narayan 2015; Bott et al. 2021) framework for MHDPIC (Sun & Bai 2023), which continuously drives the CR pressure anisotropy to trigger the CRPAI. This approach may mimic a local patch in compressible turbulence, or a mixing layer between two gas phases, or perhaps also including background shear, where a background changing field drives CR pressure anisotropy. In this work, we employ the expanding/compressing box to drive the CRPAI, balanced by ion-neutral damping to achieve a steady state. This allows us to quantitatively measure the effective CR scattering rate in the alternative scenario of CRPAI, establishing their dependence on local medium properties. This study thus fills a major gap in understanding the microphysics of CR feedback.  \nThis paper is structured as follows: Section 2 outlines the formulation and numerical setup. In particular, we introduce the adaptive δf method in Section 2.2. We calculate the expected CR scattering rates under our simulation setup based on QLT in Section 3. The simulation results are presented in Section 4, which also offers quantitative calibration of the CR scattering rates. We discuss the implications of the results towards realistic systems in Section 5. In Section 6, we conclude and discuss prospects for future research.", '2. FORMULATION AND NUMERICAL METHODS': 'In this section, we briefly review the MHD-PIC formulation in the expanding box framework (Sun & Bai 2023), and the treatment of ion-neutral damping (Plotnikov et al. 2021; Bambic et al. 2021; Bai 2022). We next propose a novel adaptive δf method to reduce noise in our simulations, before describing our simulation setup.', '2.1. The governing equations': "The MHD-PIC method treats CRs as super-particles, following the standard approach of PIC, and solves MHD equations for the thermal gas (here specifically just for the thermal ion plasma, see later in this subsection). Each simulation particle carries both position information x and normalized momentum information ( p /m ). The CR kinetic equations are given as:  \nd x = v ,  \nv ≡ ( p /m ) γ , γ ≡ √ 1 + ( p /m ) 2 / C 2 .  \nd t (1) d( p /m ) d t = ( q mc ) ( c E + v × B ) , (2) (3)  \nHere, q/ ( mc ), v , γ and C represent the CR charge-tomass ratio, particle velocity, Lorentz factor and the numerical speed of light, respectively. The CR population meanwhile exerts a Lorentz force on the MHD gas as CR backreaction. The MHD equations with the CR backreaction read,  \n∂ t ρ i + ∇· ( ρ u i ) = 0 , ∂ t ( ρ i u i ) + ∇· ( ρ i u T i u i -B T B + P ) = -ν IN ρ i u i -( Q CR c c E + j CR c × B ) , ∂ t B = -∇× ( c E ) , c E = -u i × B , (4)  \nHere, the variables include ρ i for ion gas density, P ≡ ( ρ i c 2 s + B 2 / 2) I for total pressure with I being the identity tensor, c s for the isothermal sound speed, and u i for the ion gas velocity. The unit for magnetic field is normalized such that magnetic permeability equals one, thus absorbing factors of 1 / √ 4 π . We bypass the CR backreaction on the MHD gas energy in this work by applying an isothermal equation of state (see the reason in Section 3.2). The cosmic ray charge density Q CR and the current density j CR are both functions of the phase space distribution function f ( t, x , p /m ) for the  \nCR population:  \nQ CR c ≡ ( q mc ) m ∫ f d 3 p = ( q mc ) ρ CR , j CR c ≡ ( q mc ) m ∫ v f d 3 p , (5)  \nHere, m represents the mass of an individual CR particle, and ρ CR denotes the mass density of the CR ensemble.  \nGiven that the Lorentz force is approximately proportional to ρ CR , the intensity of the CR backreaction is evaluated through the mass density ratio between CR and thermal ionized gas, ρ CR /ρ i . With the CR cyclotron frequency being Ω ≡ ( q/ ( mc )) | B | , the characteristic rate for CR feedback is ( ρ CR /ρ i )Ω.  \nWe consider the ion-neutral friction as the wave damping mechanism in this work. Ion-neutral damping prevails as the primary damping mechanism in the cold dense medium (Zweibel 2013; Armillotta et al. 2021, 2022; Hopkins et al. 2022), where the ionization fraction ρ i /ρ neu ≪ 1. Under typical ISM conditions, the ion-neutral damping can be treated such that we essentially only simulate the ion gas subjecting to a frictional force -ν IN ρ i u i on top of the background static neutrals (Plotnikov et al. 2021), as included in Equation 4. Here ν IN represents the collision frequency of neutrals with ions. This approach remains valid only in the short-wavelength regime, where the Alfv'en wave frequency ω ( k ) = kU A for the relevant wavenumber k exceeds ν IN (Reville et al. 2007; Plotnikov et al. 2021). Here, U A ≡ B/ √ ρ i stands for the Alfv'en velocity. For the ∼ GeV CRs we focus on in this work,the resonant wavenumber k res ∼ Ω m/p , varies around ∼ AU -1 , much shorter than ν IN /U A in typical ISM conditions (Plotnikov et al. 2021), justifying our approach. We also note that the induction equation retains its ideal MHD form.  \nIn this work, we drive CR pressure anisotropy by mimicking a local uniform patch of gas among the macroscopic system (e.g., ISM/CGM, etc.) undergoing expansion/compression. This is achieved under the expanding box framework (see derivations in Grappin et al. 1993; Sironi & Narayan 2015; Sun & Bai 2023), and CR pressure anisotropy is developed owing to the conservation of the magnetic moment of individual CR particles. This framework has been implemented in the most general manner in the ATHENA++ MHD-PIC module as described in Sun & Bai (2023), where background expansion rate can be flexible in any of the three directions. In the simulations presented here, the background magnetic field lies along the x -direction. We let the background expand or compression at the rate a ( t ) in the y - and z -directions and at a rate of a 2 ( t ) in the x -direction. This configuration enables MHD perturbations to propagate freely  \nalong the background magnetic field line (see Appendix of Sun & Bai 2023).  \nWe reformulate the MHD equations and the CR equations of motion in comoving coordinates x ' of the expanding box while expressing the MHD gas quantities and CR momentum in the lab frame. The equations we eventually solve in this study are as follows:  \n∂ t ρ i + ∇ ' · ( ρ u i ) = -4 ρ i ˙ a a , (6)  \n= -ν IN ρ i u i -( Q CR c c E + j CR c × B  \n∂ t ( ρ i u i ) + ∇ ' · ( ρ i u T i u i -B T B + P ) -4 ∂ t a a ρ i u i -ρ i u i · D ,  \n) (7)  \n∇ ' · B = 0 , (8)  \n∂ t B -∇ ' × ( c E ) = -4 ∂ t a a B + D · B , (9)  \nc E = -u i × B , (10)  \nA · d x ' d t = v , (11)  \nd( p /m ) d t + D · p m = q mc ( c E + v × B ) , (12)  \nwhere,  \n∇  \n' = ( ∂ a 2 ∂x ' , ∂ a∂y ' , ∂ a∂z ' ) , D = diag ( 2 ∂ t a a , ∂ t a a , ∂ t a a ) A = diag ( a 2 , a, a ) .  \n,  \nWe note that even we use an isothermal equation of state, the background expansion or compression brings about a thermodynamic impact on the background MHD gas, following an adiabatic polytropic process (Squire et al. 2020), c s ∝ a -4 / 3 ( t ). The expansion (compression) dynamically reduces (amplifies) the background magnetic field at a rate of ˙ B/B = -2˙ a/a , and the time evolution of the background field strength is  \nB g ( t ) = B g (0) exp ( -2 ˙ a a t ) .  \nConsequently, the CR cyclotron frequency Ω slowly varies over time. Both the gas density and CR density in the comoving frame changes at a rate of -4˙ a/a with time. On the other hand, the mass density ratio ρ CR /ρ i , and the Alfv'en speed U A remains constant throughout.", '2.2. Adaptive δf method': 'Inherent to a particle-based representation of the CR distribution function f , the MHD-PIC method is subject to Poisson noise, which represents a major issue  \nin the simulation of the CR gyro-resonant instabilities (Bai et al. 2019). Simply increasing the particle number N becomes unfeasible as the noise level only reduces as N -1 / 2 which is very inefficient. This difficulty can be overcome by the δf method (e.g. Parker & Lee 1993; Hu & Krommes 1994; Denton & Kotschenreuther 1995; Kunz et al. 2014; Bai et al. 2019). In the δf method, the phase space distribution f ( t, x , p /m ) is divided into two components: f 0 ( t, x , p /m ) as a known analytical form, and δf ( t, x , p /m ), representing the deviation δf ≡ f -f 0 . Consequently, the CR backreaction also splits into two contributions: one associated with f 0 , which can be evaluated analytically, and the other from δf . The simulation particles, along with statistical noise, only contribute to the backreaction from the δf part. As a result, when counting the CR backreaction (Equation 5) from δf , the contribution from each simulation particle is changed by multiplying a weight factor w ,  \nw = 1 -f 0 ( t, x ( t ) , p ( t )) f ( t, x ( t ) , p /m ( t )) . (13)  \nTraditionally, f 0 is often fixed as the initial CR distribution (e.g., Bai et al. 2019; Bai 2022) or as the CR equilibrium distribution after adiabatic evolution (e.g., Sun & Bai 2023). With the expanding box, however, the bulk anisotropy level in f keeps changing due to both the ˙ B/B driving and quasi-linear evolution. This would render the δf method less effective in mitigating statistical noise. Therefore, we adaptively adjust f 0 towards matching the evolving distribution f , so that the statistical noise remains at a relatively low level.  \nBy design, the adjustment of f 0 can be arbitrary. In our case, we initialize the CR particles with an isotropic κ distribution  \nf iso ( p ) = ρ CR m ( πκp 2 0 ) 1 . 5 G ( κ +1) G ( κ -0 . 5) ( 1 + p 2 κp 2 0 ) -κ -1 , (14)  \nwhere G (), κ and p 0 denote the Gamma function, the power index, and the peak momentum, respectively. As the background expands/compresses, the CR distribution becomes anisotropic:  \nf aniso ( t, p ) = ρ CR ξ 4 ma 4 ( πκp 2 0 ) 1 . 5 G ( κ +1) G ( κ -0 . 5) ( 1 + p 2 κp 2 0 ( ξ 4 µ 2 -ξ 2 µ 2 + ξ 2 ) ) -κ -1 , (15)  \nwhere µ and ξ represent the CR pitch angle and anisotropy level. The factor a 4 in the denominator accounts for the reduction in lab frame density resulting  \nfrom the background expansion (see previous subsection). The pitch angle anisotropy is significant for particles whose p ≳ p 0 but becomes negligible when p ≪ p 0 , as the κ -distribution yields f ( p ≪ p 0 , µ ) ∼ const, which changes little in response to expansion/compression. Note that initial κ distribution corresponds to a = 1 and ξ = 1. Under adiabatic evolution without MHD perturbations, the CR distribution f should precisely follow Equation 15 with ξ = a ( t ) (see Appendix of Sun & Bai 2023). However, in our simulations, the MHD waves resulting from the CRPAI tend to isotropize the CR distribution due to quasi-linear diffusion, reducing the anisotropy level at certain energy ranges. Additionally, during scattering with MHD waves, the bulk CR population can lose or gain energy, and p 0 may shift. We thus parameterize the above form of the anisotropic κ distribution (Equation 15) as the ansatz for f 0 , and dynamically estimate ξ and p 0 , by taking the generalized method of moments for Equation 15 after treating f ( t, p ) as f aniso :  \nξ = 2 π ∫ d 3 x d 3 p fp √ 1 -µ 2 ∫ d 3 x d 3 p fp | µ | ,  \np 0 = πκ ( κ -1) G ( κ -0 . 5) 2 G ( κ +1) ∫ d 3 x d 3 p fp √ ξ 4 µ 2 -ξ 2 µ 2 + ξ 2 ∫ d 3 x d 3 p f . (17)  \n(16) √  \nNote that estimating ξ only requires the information about f , but the estimation on p 0 is made convenient when ξ is already known. We also mention that while ρ CR and κ can also be adaptive free parameters, in this work, they are taken as fixed. We expect the CR distribution to be largely uniform in our simulations and thus ρ CR is known from the CR mass conservation in the comoving frame. This is also why the integration above is conducted over space which also serves to substantially reduce noise. We also find changes in the parameter κ in our simulations to be minor.  \nThe evaluation of ξ and p 0 requires integration over f , which is again done by the δf method using f 0 from the previous time-step, denoted as f 0 , old . The result is  \nξ new = 2 π × ( √ πp 0 , old G ( κ +1) N CR 2 ξ old √ κ ( κ -1) G ( κ -0 . 5) + ∫ d 3 x d 3 p δf old p √ 1 -µ 2 ) / ( p 0 , old G ( κ +1) N CR ξ 2 old √ πκ ( κ -1) G ( κ -0 . 5) + ∫ d 3 x d 3 p δf old p | µ | ) , (18)  \nN CR ≡ ∫ d 3 x ρ CR ma 4 .  \nand, we utilize ξ new to calculate p 0 , new :  \np 0 , new = [ √ πκ ( κ -1) G ( κ -0 . 5) 2 G ( κ +1) × ∫ d 3 x d 3 p δf old p √ ξ 4 new µ 2 -ξ 2 new µ 2 + ξ 2 new + p 0 , old N CR 2 × { α 2 + sinh -1 √ α 2 -1 √ α 2 -1 , α ≡ ξ new ξ old > 1 α 2 + sin -1 √ 1 -α 2 √ 1 -α 2 , α ≤ 1 ] /(∫ d 3 x d 3 p δf old + N CR ) . (19)  \nWe have validated the above implementation by observing the adiabatic evolution of f in an expanding box in Appendix A.1. Our implementation for the adaptive δf method attains a similar signal-to-noise level in the original δf method but with more than four times the number of simulation particles (Appendix A.2). The computational cost for adaptively fitting f 0 is not negligible. We thus fit once after a fixed time interval ∆ T adapt .', '2.3. Numerical setup': "Our simulation setup closely follows that of Bai et al. (2019) and the CRPAI simulations of Sun & Bai (2023), which we briefly describe here. We initialize an isotropic homogeneous κ distribution for CRs in a 1D expanding box (Equation (14)) along the x -direction and with a constant collision frequency for ion-neutral damping ν IN . The background gas (thermal ions) has uniform density ρ 0 and uniform magnetic field B 0 along ˆ x . The expansion rates along the three directions (ˆ x, ˆ y, ˆ z ) are set to ( a 2 ( t ) , a ( t ) , a ( t )) as mentioned earlier, with a ( t ) = exp(˙ at ). The background gas is initialized with a series of forward and background-traveling, left- and rightpolarized Alfv'en waves with random phases, with initial wave intensity spectrum as | k | I ( k ) /B 2 0 = A 2 , where k is the wavenumber, so that the wave energy is equally distributed in logarithmic k space, and the normalization is such that the total wave energy integrated over k is a fraction 2 A 2 ln 10 of the background field energy. For the CRPAI, the most unstable wavenumber is given by k 0 ∼ (Ω m ) /p 0 (Sun & Bai 2023). Our simulation box size L x is thus chosen to accommodate multiple most unstable wavelengths in the comoving frame of the expanding box. In the simulations, the numerical units include the initial mass density ρ 0 = 1, the initial intensity of the background magnetic field B g ( t = 0) = B 0 = 1, the Alfv'en velocity of thermal ions U A ≡ B 0 / √ ρ 0 = 1, the initial cyclotron frequency Ω 0 ≡ ( q/ ( mc )) B 0 = 1, and the initial ion inertial length λ i = Ω 0 /U A = 1.  \nTo integrate the governing equations (Equation 6 to Equation 12), we employ the two-stage van Leer time integrator (Stone & Gardiner 2009) in ATHENA++ (Stone  \net al. 2020; Sun & Bai 2023). We use the Roe solver (Roe 1981) with third-order reconstruction to solve the MHDequations, and the Boris pusher (Boris et al. 1972) to integrate CR particles. Interpolation of fluid quantities and CR backreaction deposition both follow the triangular-shaped cloud scheme (Birdsall & Langdon 2004). To enhance the numerical accuracy in CR backreaction, we enable the adaptive δf method (see Section 2.2) with ∆ T adapt = 500Ω -1 0 , and divide the CR population into eight momentum bins (e.g. Bai et al. 2019) which span the momentum range from p 0 / 500 to 500 p 0 in a logarithmic uniform manner. Periodic boundary condition is applied to all fluid quantities and CR particles, except that the gyro phase of particles are randomized upon box-crossing to facilitate quasi-linear evolution (Bai et al. 2019).  \nGiven the extreme level of scale separation in realistic ISM or similar conditions 1 (e.g., U A /c ∼ 10 -4 , ρ CR /ρ 0 ≲ 10 -6 ), it is impractical to adopt such values in our simulation parameters. Instead, we choose parameters in the same regime as realistic conditions, and the results should can be rescaled for real applications. Our default simulation parameters include: C = 200 U A , p 0 = m C , κ = 1 . 25 (so that f ( p ) ∝ p -4 . 5 at large p ), ρ CR = 4 × 10 -5 ρ 0 , ν IN = 2 × 10 -5 Ω 0 , ˙ a = ± 5 × 10 -7 Ω 0 . Note that our choice of ˙ a and ν IN is such that ˙ a ≪ ν IN ≪ ω ( k 0 ) ≪ Ω 0 . This ensures that compression/expansion is adiabatic, and that ν IN is much lower than the frequency of the Alfv'en modes in our simulation box (Plotnikov et al. 2021). Moreover, the anticipated CR anisotropy level ∣ ∣ ξ 2 -1 ∣ ∣ ∼ ( ν IN / Ω 0 ) ( ρ 0 /ρ CR ) ( U A /c ) is kept low at ∼ 10 -2 (see Section 3.1). Our simulation box size is chosen to be L x = 1 . 92 × 10 5 U A / Ω 0 , corresponding to about 150 most unstable wavelengths ( ≈ 2 πp 0 /m Ω 0 ), with a resolution of ∆ x = 5 U A / Ω 0 (about 250 cells per most unstable wavelength). In each cell, we initialize 64 simulation particles per momentum bin. Similar as our earlier works, we adopt initial wave amplitudes with A = 10 -3 / 3 to seed instability growth. Our simulations typically run for a duration of 4 . 5 × 10 5 Ω -1 0 .  \nWe further conduct a brief parameter survey, where we reduce each of the parameters ( C , ρ CR , ν IN , ˙ a ) by half while keeping the other physical parameters fixed to the fiducial values. Additionally, we carry out highresolution simulations with the fiducial physical parameters but a reduced MHD cell size of ∆ x = 2 . 5 U A / Ω 0 , and using 48 simulation particles per momentum bin per  \ncell. For the run with C being half of the fiducial value, the most unstable wavelength is also halved. To achieve equivalent resolution of the most unstable wave, we use the same MHD cell size and number of particles per cell as in the high-resolution run.", '3. THEORETICAL FRAMEWORK': 'Under our simulation setting, the expanding/compressing box continuously drives CR anisotropy, leading to the development of the CRPAI. The wave growth is countered by the ion-neutral damping, reaching a balance in steady state. In the meantime, the waves attempt to isotropize the CRs, balancing the driving force. In this section, we estimate the steady state anisotropy level, wave amplitudes, and the resulting CR scattering rates by quasi-linear theory (QLT). The result is to be compared to and calibrated by our simulation results to be presented afterwards, which can eventually be incorporated into CR (magneto-)hydrodynamics as subgrid models for CR scattering.', '3.1. Prediction on CR anisotropy level and effective scattering rate at the saturated state': "The saturated state is characterized by two balance relations. First, the waves reach a stable amplitude through the competition between CRPAI and damping mechanisms. We extend the analytical study of CRPAI (Sun & Bai 2023) from the non-relativistic regime to the relativistic regime in the expanding box, 2  \nΓ growth ( k ) ≈-ρ CR 2 ρ 0 Ω 0 [ 1 ± ( ξ 2 -1 ) Ω 0 kU A ] Q 2 ( k ) ξ 2 , ≈∓ ρ CR 2 ρ 0 Ω 0 ( ξ 2 -1 ) Ω 0 kU A Q 2 ( k ) ξ 2 , (20)  \nwhere  \nQ 2 ( k ) ≡ √ π κ 1 . 5 G ( κ +1) G ( κ -0 . 5) m Ω 0 kp 0 ( 1 + m 2 Ω 2 0 κk 2 p 2 0 ) -κ . (21)  \nThe linear growth rate Γ growth depends on the wavenumber k and the ± sign denotes the polarization direction. The fastest growth is achieved at k 0 ∼ Ω 0 m/p 0 regardless of the propagation direction. We have also dropped 1 compared to ( ξ 2 -1)Ω 0 /kU A in the bracket of Equation 20. This is because growth is the most prominent  \nnear k ∼ k 0 , making Ω 0 /kU A ∼ C /U A ≫ 1, thus even a weak level of anisotropy would make the second term in the bracket dominate over 1. It should also be noted that this estimation on growth rate assumes that the CR distribution function has the form Equation 15, with a uniform anisotropy level at all momenta characterized by ξ , which is not necessarily true as the CR particles evolve. Therefore, this estimate should only be taken as a proxy.  \nThis growth is to be balanced by the wave damping rate. With our approach of ion-neutral damping (Plotnikov et al. 2021), the damping rate is largely constant for all wavelengths we consider ( ν IN ≪ kU A ),  \nΓ damp ( k ) ≈ ν IN 2 . (22)  \nBalancing growth and damping Γ growth = Γ damp , we can estimate the CR anisotropy level 3 as  \n∣ ∣ ξ 2 -1 ∣ ∣ ( p ) ≈ kU A Ω 0 ν IN Ω 0 ρ 0 ρ CR Q -1 2 ( Ω 0 m p ) . (23)  \nwhich, interestingly, depends only on ν IN but not ˙ B/B . This is the consequence of linear damping, and is similar to the case of CR streaming instability, where the expected streaming speed does not depend on CR pressure gradient under linear damping (e.g. Wiener et al. 2013; Bai 2022). Here, we further approximate ξ to be an explicit function of p (instead of being constant by definition), by substituting k with Ω 0 m/p for the resonance condition and omitting the dependence on the particle pitch angle µ .  \nFurthermore, when treating the CR population as a single fluid, we can then re-express the CR anisotropy level in terms of CR pressure P CR ,  \nξ 2 = P CR , ⊥ P CR , ∥ ≡ ∫ d 3 p fp 2 ( 1 -µ 2 ) 2 ∫ d 3 p fp 2 µ 2 ≈ 1 ± 1 . 78 U A c ν IN Ω 0 ρ 0 ρ CR , (24)  \nwhere ⊥ and ∥ denote the components perpendicular and parallel to the background magnetic field, respectively. Here the factor '1.78' is derived from Q -1 2 (Ω 0 m/p 0 ), where the CRs at the peak of the distribution are taken to be representative of the entire population, as commonly assumed.  \nThe second balance is on the CR distribution, where driving on CRs by ˙ B/B should be balanced by CR scattering, which can be expressed under the quasilinear theory. Using the Fokker-Plank equation, we can  \nwrite (e.g. Jokipii 1966; Kulsrud & Pearce 1969; Skilling 1975c; Schlickeiser 2002)  \n∂ t f + ˙ B 2 B [ p ( 1 + µ 2 ) ∂ p f + µ ( 1 -µ 2 ) ∂ µ f ] = ∂ µ [( D µµ ∂ µ f +sgn( µ ) D pµ ∂ p f )] + 1 p 2 ∂ p [ p 2 (sgn ( µ ) D pµ ∂ µ f + D pp ∂ p f ) ] , (25)  \nwhere,  \nD µµ = ν ( 1 -µ 2 ) 2 , D pµ = ν ( 1 -µ 2 ) 2 pU A v , D pp = ν ( 1 -µ 2 ) 2 ( pU A v ) 2 . (26)  \nThe term proportional to ˙ B/B corresponds to the anisotropy driving, ( D · p ) · ∂f/∂ p (see Equation 12). Note that it not only affects the pitch angle, but also the total momentum. The two terms on the right-hand side represent quasi-linear diffusion (QLD) in the gas comoving frame, resulting from the pitch angle diffusion and CR momentum diffusion when CRs scatter with the Alfv'en waves. The diffusion coefficients, D µµ , D pµ , and D pp , are all characterized by a single parameter, the scattering rate ν ( p, µ ) (e.g. Jokipii 1966; Kulsrud & Pearce 1969; Skilling 1975c; Schlickeiser 2002) which, under quasi-linear theory, is given by  \nν ( p, µ ) = π Ω γ Ω m pµ I ( k res = Ω m pµ ) /B 2 g . (27)  \nIn the saturated state of the CR population, we anticipate ∂ t ln ( P CR , ⊥ /P CR , ∥ ) = 0. We can integrate the above over pitch angle µ . The driving term yields,  \ndriving = ˙ B 2 B ( ∫ d µ [( 1 -µ 4 ) p∂ p f -( 1 -6 µ 2 +5 µ 4 ) f ] ∫ d µ (1 -µ 2 ) f -∫ d µ [( µ 2 + µ 4 ) p∂ p f -( 3 µ 2 -5 µ 4 ) f ] ∫ d µµ 2 f ) .  \nAssuming the CR distribution is close to isotropy (see Section 2.3), P CR , ⊥ ≈ P CR , ∥ , we may ignore the dependence of f on µ as a high-order term and simply use the isotropic distribution function, f ( p, µ ) ≈ f iso ( p ), which leads to  \ndriving ≈ -˙ B 5 B ∂ ln f iso ∂ ln p (28)  \nTo simplify the right-hand side of Equation 25, we firstly eliminate the momentum diffusion terms (inside  \nthe second square bracket) as they are of higher order compared to the pitch angle diffusion by a factor of U A /c . We proceed by applying a partial integration,  \nQLD = -9 2 f iso ∫ d µ ν ( µ -µ 3 ) 2 ( ∂ µ f +sgn( µ ) pU A v ∂ p f ) .  \nTo integrate over µ , we introduce the following assumptions: the slight anisotropic distribution takes the form of f ( p, µ ) = f ( p √ ξ 4 µ 2 -ξ 2 µ 2 + ξ 2 ) ; the anisotropy level ( ξ 2 -1 ) ≪ 1; we simplify the scattering rate as an effective rate independent of µ ,  \nν eff ( p ) ≡ ∫ 1 -1 ( 1 -µ 2 ) ν ( p, µ ) d µ ∫ 1 -1 (1 -µ 2 ) d µ . (29)  \nUnder these assumptions, the contribution from QLD can be written to a form similar to Equation 28,  \nQLD ≈ -ν eff ( 3 5 ( ξ 2 -1 ) + 9 8 U A c ) ∂ ln f iso ∂ ln p . (30)  \nBy inserting the wave saturation condition (Ω 0 /k ∼ p/m in Equation 20) into the CR anisotropy balance condition (driving = QLD), we build up the relation between the effective scattering rate and environmental parameters for each CR momentum bin,  \nν eff ( p ) ≈ 1 3 Q 2 (Ω 0 m/p ) ∣ ∣ ∣ ∣ ∣ ˙ B B ∣ ∣ ∣ ∣ ∣ c U A Ω 0 ν IN ρ CR ρ 0 . (31)  \nThe momentum-dependent factor Q 2 peaks around 0 . 56 (for κ = 1 . 25, the exact value depends on κ ). In the conventional single-fluid treatment of the CRs, peak growth rate is often employed when calculating the balance between wave growth and damping, which yields  \nν eff ≈ 0 . 187 ∣ ∣ ∣ ∣ ∣ ˙ B B ∣ ∣ ∣ ∣ ∣ c U A Ω 0 ν IN ρ CR ρ 0 . (32)  \nCombining Equation 27 and Equation 31, we can further estimate the saturation level of the CRPAI as:  \n( δB B g ) 2 ∼ kI ( k ) B 2 g ∼ √ 1 + Ω 2 / ( kc ) 2 3 π Q 2 ( k ) ∣ ∣ ∣ ˙ B/B ∣ ∣ ∣ Ω 0 c U A Ω 0 ν IN ρ CR ρ 0 ≪ 1 , (33)  \nThe theoretical calculations above have been significantly simplified under a number of approximations, and are intended to provide order-of-magnitude estimate on how the saturated state responds to the environment. Weanticipate the results to provide the appropriate scaling relations, but the pre-factors can be subject to major uncertainties, which we aim to calibrate through our simulations (Section 4).", '3.2. Implication of effective scattering rate in CR (magneto-)hydrodynamics': 'At macroscopic scales (e.g., in the Galaxy), the cosmic rays are typically modeled as fluid, known as CR (magneto-)hydrodynamics, characterized by the CR energy density E CR and the CR energy flux density F CR ,  \nE CR ( t, x ) ≡ ∫ d 3 p f ( t, x , p ) √ c 2 + p 2 c ≡ ∫ dln p E CR ( t, x , p ) , F CR ( t, x ) ≡ ∫ d 3 p f ( t, x , p ) p c 2 ≡ ∫ dln p F CR ( t, x , p ) , (34)  \nwhere we can either approximate CRs as a single fluid or multiple fluids distributed across different energy (momentum) bins. They are coupled with background gas through wave-particle interaction, described by the scattering rates.  \nIn our simulation setting with expanding box, the outcome of the CRPAI in balance with ion-neutral damping only affects E CR , as F CR = 0 in our setups. By integrating the Fokker-Plank equation (Equation 25) over the CR population, and applying similar techniques as described in Section 3.1, the CR fluid energy equation reads  \n∂ t E CR -8 3 ˙ B B E CR = U A c ν eff E CR ( 1 2 ( ξ 2 -1 ) + 4 3 U A c ) . (35)  \nThe term 8 ˙ B/ (3 B ) E CR corresponds to the external driving through adiabatic compression/expansion, and the right-hand side corresponds to the quasi-linear diffusion. Note the different pre-factors compared to those in Section 3.1 due to additional weighting in the moment equations. The above equation can also be naturally re-expressed in the p by p treatment, by substituting E CR ( t, x ) with E CR ( t, x , p ).  \nIn CR (magneto-)hydrodynamics, the anisotropy level ξ 2 and ν eff are expected to be user-specified parameters. Our estimates in Equations 23 & 31 can be considered as sub-grid prescriptions under the assumption that ionneutral damping dominates. Under such prescriptions, the QLD effect in CR energy density equation (the righthand side of Equation 35) can be reformulated in terms of environmental parameters as,  \n∂ t E CR -8 3 ˙ B B E CR ∼ -1 6 ˙ B B U A c E CR ( p ) , (36)  \nwhere we expect the ( ξ 2 -1) term to dominate over the U A /c term on the right hand side of Equation 35. We  \nsee that wave scattering on CRs (which is irreversible) always reduces the adiabatic cooling/heating on CR energy density from external driving, but it is only a minor effect of the order ( U A /c ) / 16. Therefore, for the pure CRPAI case, this effect is largely negligible on CR energy density evolution in our simulation setup (Yan & Lazarian 2011; Zweibel 2020).  \nThis work aims to calibrate ξ 2 and ν eff in the CR energy equation (Equation 35) using kinetic simulations. The CR fluid quantities ( P CR , ⊥ , P CR , ∥ , and E CR ) can be measured by definition. The effective scattering rate of CRSI has been previously measured through the steadystate condition in the CR flux equation (Bai 2022). However, in our simulation, due to the absence of an equilibrium state for the CR energy (Equation 36), we cannot perform a similar measurement of ν eff through Equation 35. Instead, we measure ν by definition (Section 4.2) and then numerically integrate the FokkerPlanck equation (Equation 25) for the QLD effect on CR energy density. By comparing Equation 35 with the numerical integration of Equation 25, we ultimately obtain ν eff (Section 4.3).  \nNote that as CR streaming is not present in our problem setup, we do not include the equation for the CR energy flux. In the more general case with CR streaming, the scattering rate will also affect the evolution of the CR energy flux, and ν eff is likely determined by the combined outcome of the CRSI and CRPAI. This is left for our future work.', '4. SIMULATION RESULTS': 'In this section, we illustrate the instability growth and the establishment of steady CR anisotropy in our simulations. The linear growth of CRPAI has been demonstrated in Sun & Bai (2023). Here, our primary focus lies on the subsequent evolution leading towards saturation. The final outcome is the effective scattering rate ν eff at the saturated state. To compute ν eff , we assess the CR scattering rate by definition and subsequently average it over the CR pitch angle and energy (momentum magnitude). Given Equation 31 is highly approximate, we want to test its scaling on environmental parameters and calibrate the coefficients, by comparing different runs at the end of this section.', '4.1. Fiducial run': 'We start by showing in Figure 1 the time evolution of the wave energy density for all our simulation runs. The initial seed waves are first damped by ion-neutral damping. In the meantime, as the CR anisotropy gradually develops owing to ˙ B/B , the CRPAI starts to overcome damping, eventually leading to wave growth.  \nFigure 1. Time evolution of the wave energy density around k 0 , ∫ 5 k 0 k 0 / 5 I ( k )d k . The solid lines represent the compressing box, while the dashed lines denote the expanding box. Line colors distinguish simulations with varying parameters and resolutions. Lines from the runs sharing similar saturated states largely overlap.  \n<!-- image -->  \nGiven our fiducial parameters, the turning point is expected to occur when the CR anisotropy level reaches ˙ at ∼ ( U A / C )[1 + ( ν IN / Ω 0 )( ρ 0 /ρ CR ) Q -1 2 (Ω 0 m/p )], yielding t ∼ 1 . 5 × 10 4 Ω -1 0 for the fastest growing mode in the fiducial case. The actual turn-over when summing over all waves starts later around t ∼ 4 × 10 4 Ω -1 0 as it takes longer for waves at other wavelength (which grow slower) to catch up. After a few e-folding time, the waves are sufficiently grown to efficiently scatter and isotropize the CR particles. This will reduce the CR anisotropy and hence the wave growth, eventually leading to the saturation of the wave amplitudes, where wave growth due to CRPAI balances ion-neutral damping.', '4.1.1. The wave spectrum': 'To look more closely into the processes, we decompose the waves into the Fourier space, according to polarization and propagation directions (Bai et al. 2019), and plot the wave energy spectra I ( k ) of the fiducial runs at saturated state in Figure 2. In our setup, the expansion ( ˙ B/B < 0) yields an oblate CR distribution ( P CR , ⊥ /P CR , ∥ > 1), where the CRPAI amplifies lefthanded waves propagating in opposite directions and damps the right-handed counterparts, while the compression ( ˙ B/B > 0) leads the waves to saturate at the opposite polarization direction due to the prolate CR distribution ( P CR , ⊥ /P CR , ∥ < 1). The dependence of wave intensity on k are consistent with theoretical expectations of the CRPAI. The wave growth rate peaks at  \n<!-- image -->  \nFigure 2. Wave intensity spectra, kI ( k ), during the saturated state of the fiducial runs at t = 4 × 10 5 Ω -1 0 . The top panel corresponds to the expanding box (˙ a = 5 × 10 -7 Ω 0 ), while the bottom panel corresponds to the compressing box (˙ a = -5 × 10 -7 Ω 0 ). Solid lines represent forward propagating waves, while dashed lines denote backward propagating ones. The left- and right-handed branches are distinguished by black and red markings. Lines in the same color (polarization) largely overlap. The grey dotted lines refer to the theoretical estimate of the saturated wave intensity (Equation 33), peaking at the most unstable wave number k 0 ∼ 5 × 10 -3 Ω 0 /U A for the fiducial runs. The accompanying animation, covering the time span from t = 0 to 4 . 5 × 10 5 Ω -1 0 , is available in the online version of the journal.  \n<!-- image -->  \nthe resonant wave number of particles with p = p 0 given by k 0 ∼ (Ω m ) /p 0 = 5 × 10 -3 Ω 0 /U A (Equation 20). Consequently, after an initial increase in the CR anisotropy level, the waves around k 0 start to grow, and first reach the saturated amplitude (due to QLD, to be discussed next). Waves at other wave numbers require more time to grow and are expected to reach saturation slower with lower amplitudes.  \nBy the end of the simulations we observe that waves spanning from k 0 / 5 to 5 k 0 have approximately reached  \n<!-- image -->  \n(a)  \nFigure 3. The CR distribution function ( f/ ⟨ f ⟩ µ -1 ) for the fiducial runs in the expanding box (top) and the compressing box (bottom) at various evolution snapshots. Dashed lines represent particle momenta resonating with the same wave characterized by a constant k = Ω 0 m/ ( pµ ). The animation spanning from t = 10 3 to 4 . 5 × 10 5 Ω -1 0 is accessible in the online version of the journal.  \n<!-- image -->  \nFigure 4. The CR pressure anisotropy level as a function of CR energy (momentum) at different snapshots of the fiducial runs. The solid lines represent the case in compressing box, while the dashed lines denote the expanding box case. Transparency indicates the simulation run time. The black dotted line corresponds to the theoretical prediction given by Equation 23.  \n<!-- image -->  \nsaturation in the simulations (thus in Figure 1 we actually show the evolution of wave intensity only within this wavenumber range). We also observe from Figure 2 that the wave intensity in this spectral range approximately agree with the estimation from QLT (Equation 33), with small deviation by up to a factor of ∼ 2. Note that the final wave intensity slightly differs between the expanding box and the compressing box. We attribute this to the reduction/enhancement of background magnetic field due to box expansion/compression, which leads to order unity difference in B g within simulation time.', '4.1.2. Anisotropy in the CR distribution function': 'Accompanied by the saturation of wave energy is the asymptotic saturation of the CR distribution f ( p, µ ). Its evolution for the fiducial runs is illustrated in Figure 3, which is normalized by the isotropic distribution ⟨ f ⟩ µ ( p ) averaged over pitch angle. Initially, MHD waves are too weak to scatter the CRs, and hence f ( p, µ ) evolves adiabatically due to expansion/compression approximately following the anisotropic κ distribution with ξ ∼ a ( t ) (Equation 15): CRs with p ≳ p 0 concentrate toward µ ∼ 0/ µ = ± 1 in expanding/compressing box, while CRs with p ≪ p 0 in the κ distribution do not rapidly develop significant anisotropy (see Section 2.2). Following the subsequent development of the CRPAI, waves around k 0 grow fastest and effectively scatter particles whose momentum around p 0 . More precisely, scattering is the most effective for particles satisfying Ω m/ ( pµ ) ∼ k 0 (see dashed lines in Figure 15), leading to pitch angle diffusion and hence isotropization.  \nThe isotropization gradually propagates towards CRs at lower and higher energies as waves in other wavelengths grow. However, even within the momentum range where scattering is the most effective (around p = p 0 ), the particles are not fully isotropized, leaving a mild anisotropy, reflecting the balance with the box expanding/compressing. The CR anisotropy for particles with p ∈ (0 . 5 p 0 , 5 p 0 ) ultimately reaches a quasi-steady state by the end of the simulations in both the expanding and compressing boxes. Finally, we note that the CR energy distribution keeps evolving following Equations 35.  \nWe quantify the CR anisotropy level through its pressure anisotropy, P CR , ⊥ /P CR , ∥ -1, based on Equation 24. In the momentum-by-momentum treatment, we present the temporal evolution of CR anisotropy in Figure 4. In the expanding box, the CR pressure perpendicular to the background magnetic field consistently exceeds the parallel component, while the reverse holds in the compressing box. Owing to the inefficient driving of anisotropy towards p ≪ p 0 (consequence of the κ distribution, see footnote 4), we primarily focus on particles with p ≳ 0 . 2 p 0 where driving is effective. Initially, anisotropy levels across all CR momenta increase with time. Subsequently, quasi-linear diffusion (QLD) induced by amplified waves reduces the anisotropy, ultimately maintaining a steady anisotropic profile. The saturated anisotropy level within the momentum range (0 . 5 p 0 , 5 p 0 ) is broadly consistent with the trend predicted by Equation 23, albeit higher than the theoretical value by order unity. We attribute this deviation to the omission of pitch angle dependence in deriving Equation 23, which is the standard approach in QLT but is prone to error, and a similar situation has been discussed in Bai (2022). Lower-energy CRs in simulations fall below the predicted anisotropy level, which is likely owing to less efficient anisotropy driving. The anisotropy level in the expanding box simulation is systematically higher than that in compressing box.', '4.2. The pitch angle scattering rate and comparison with quasi-linear theory': "This subsection introduces the methodology to quantify the pitch angle scattering rate ν ( p, µ ) and subsequently compares the obtained measurements with quasi-linear theory (QLT) as a diagnostics.  \nIn our setup, two independent mechanisms contribute to the evolution of CR momenta. The expansion/compression of the box stretches the CR momentum while conserving a 2 p ∥ and ap ⊥ (See the appendix of Sun & Bai 2023). Simultaneously, MHD waves interact with CRs, inducing diffusion in both pitch angle µ  \nand momentum magnitude p , characterized by the diffusion coefficients D µµ , D pµ , and D pp , connected to the scattering rate ν ( p, µ ) (Equation 26).  \nWe quantify the diffusion coefficients by tracing simulation particles (Bambic et al. 2021) over some time interval ∆ t = t 2 -t 1 . Since the expanding/compressing box setup itself leads to changes in particle momentum and pitch angle, we counter such effects by considering particle evolution to effectively follow a 'stretch-diffusestretch' fashion. For an arbitrary simulation particle k , its momentum p 1 ,k at time t 1 is first stretched for ∆ t/ 2 to p ' 1 ,k due to box expansion/compression only. Subsequently, the waves scatter p ' 1 ,k to p ' 2 ,k , and finally CR undergoes another stretching for ∆ t/ 2 to p 2 ,k at time t 2 . Given the measured values of p 1 ,k and p 2 ,k , we obtain  \np ' 1 ,k = ( p ' x 1 ,k , p ' y 1 ,k , p ' z 1 ,k ) = ( p x 1 ,k e -˙ a ∆ t , p y 1 ,k e -˙ a ∆ t/ 2 , p z 1 ,k e -˙ a ∆ t/ 2 ) , p ' 2 ,k = ( p ' x 2 ,k , p ' y 2 ,k , p ' z 2 ,k ) = ( p x 2 ,k e ˙ a ∆ t , p y 2 ,k e ˙ a ∆ t/ 2 , p z 2 ,k e ˙ a ∆ t/ 2 ) .  \nWe then consider the particle evolution from p ' 1 ,k to p ' 2 ,k to be due to pure wave scattering, acting on the particle with momentum p ∗ = ( p ' 1 ,k + p ' 2 ,k ) / 2. We choose the measurement time interval ∆ t to be 5 × 10 3 Ω -1 0 to ensure p ' 1 ,k and p ' 2 ,k do not differ significantly. The pitch angle diffusion coefficient is thus determined by definition,  \nD µµ ( µ ∗ , p ∗ ) = 〈   p ' 2 ,k · b g ∣ ∣ ∣ p ' 2 ,k ∣ ∣ ∣ -p ' 1 ,k · b g ∣ ∣ ∣ p ' 1 ,k ∣ ∣ ∣   〉 p ∗ × 2 ∆ t , (37)  \nwhere ⟨⟩ p ∗ denotes averaging over all simulation particles at the midpoint state p ∗ . The measurements for D pµ and D pp follow the same procedure.  \nWith the pitch angle diffusion coefficient D µµ , the scattering rate, ν ( µ, p ), is systematically measured as a function of both pitch angle µ and momentum magnitude p . Figure 5a illustrates ν ( µ, p ) at the saturated state of the fiducial run in expanding box. The scattering rate exhibits a peak around the momentum p 0 , corresponding to particles scattering with the most unstable wave at k 0 , while it diminishes at µ = 0 due to the lack of resonant waves with wave numbers k ≫ k 0 (see Figure 2a). In the context of CR streaming instability, it causes the well-known 90 · pitch angle crossing problem, where this crossing is necessary to reduce the drift anisotropy (Felice & Kulsrud 2001; Holcomb & Spitkovsky 2019; Bai et al. 2019; Zeng et al. 2024). Given the symmetry with respect to 90 · pitch angle in  \nFigure 5. 5a: The CR scattering rate ν measured by definition, as a function of CR momentum p and CR pitch angle µ , corresponding to the fiducial run in the expanding box, at time around ∼ 4 . 5 × 10 5 Ω -1 0 . 5b: The comparison between D µµ and D pp measured by definition. The quasi-linear theory predicts D pp = D µµ ( pU A /v ) 2 (see Equation 25).  \n<!-- image -->  \nFigure 6. The comparison between the pitch angle scattering rate ν ( p, µ ) as a function of µ for different momenta measured by definition (dashed lines) and the one calculated from the wave spectrum (Equation 27, solid lines), for the fiducial run in the expanding box at time around ∼ 4 . 5 × 10 5 Ω -1 0 .  \n<!-- image -->  \nthe pressure anisotropy problem considered here, this crossing does not affect the isotropization process and hence we do not discuss it further.  \nThe measured value of ν remains largely consistent regardless of whether it was made through D µµ , D pµ , or D pp . For example, there exists a pre-factor difference of ( pU A /v ) 2 when calculating from D µµ and D pp (Equation 26), which is verified in Figure 5b for the fiducial run in the expanding box. For the most part, the measured ratio D µµ ( pU A /v ) 2 /D pp stays around 1 as expected. There are deviations near both µ = 0 and ± 1. This is because the diffusion coefficients theoreti-  \nat both µ = 0 and ± 1 (Equation 26), thus the corresponding measurements are considered less faithful. There is also larger numerical noise in the region p < p -0 . 5 0 due to lower wave amplitudes at highk .  \nAdditionally, we compare ν measured from simulation particles and evaluated through the wave intensity based on QLD (Equation 27) in Figure 6. Besides deviations near the µ = 0 (QLD fails) and ± 1 (measurement of D µµ is subject to error) regions, these two independent measurements exhibit consistency, thereby supporting the predictions of QLT. Note that due to the substantial noise in wave intensity, we will employ ν obtained from tracing particles in the subsequent calculations.", '4.3. Anisotropy level and effective scattering rate at the saturated state': 'The final outcomes of our study is to yield the CR pressure anisotropy level, P CR , ⊥ /P CR , ∥ -1, and the effective scattering rate, ν eff . They are expected to be given either as a single fluid, or on p -byp bases, as a function of environmental parameters ( ˙ B/B , ν IN , etc.). We study all our simulation results and compare them with the QLT results derived in Section 3.', '4.3.1. Momentum-by-momentum results': "The relationship between P CR , ⊥ /P CR , ∥ -1 at the saturated state and environmental parameters is depicted in Figure 7. The trends of the CR anisotropy level with respect to CR momentum (energy) remain largely similar after varying environmental parameters (see Section 4.1). The qualitative dependence on magnitude aligns with theoretical expectations (Equation 23). The expansion/compression rate does not significantly impact CRs at the saturated state (within the momentum range (0 . 5 p 0 , 5 p 0 )), as expected with QLT. Reducing C /U A or ρ CR /ρ 0 , which weakens CRPAI, indeed increases the anisotropy level again consistent with Equation 23, while we find the anisotropy level is more sensitive to C /U A than ρ CR /ρ 0 . A weak damping case is counterbalanced by a slow instability growth rate with a smaller anisotropy level, in agreement with results seen in Figure 7. Within identical environmental parameters, the expansion consistently induces a more effective anisotropy compared to compression, by a factor of less than 2. We also see that doubling the resolution leads to identical results as illustrated in Figure 7, thus validating numerical convergence in our simulations.  \nThe effective scattering rate ν eff , normalized with environmental parameters and as a function of the CR momentum, is illustrated in Figure 8. The effective scattering rate encapsulates the wave scattering on CR transport, by taking moments of Equation 25. We numerically integrate the Fokker-Planck equation for ν eff ,  \nFigure 7. The CR pressure anisotropy level, similar to Figure 4, but for different simulation runs at time 4 . 5 × 10 5 Ω -1 0 . The solid lines represent the compressing box runs, while the dashed lines denote the expanding box cases. The simulations with varying parameters and resolutions are given in different colors. The animation covering from t = 10 3 to 4 . 5 × 10 5 Ω -1 0 for all runs in different parameters is accessible in the online version of the journal.  \n<!-- image -->  \nleveraging the known terms in Equation 25 from the simulations. One taking moments method is to retain the dependence of ν eff on momentum (energy) p , namely the p -byp treatment (Bai 2022). This treatment divides the CR population into distinct momentum (energy) groups, with group i representing a CR fluid composed of CR particles within the momentum range from p i, min to p i, max . Through the definition given by Equation 24 & 34, we have the energy density and the anisotropy level for the group i ,  \nE CR ,i = ∫ p i, max p i, min dln p E CR ( p ) , ξ 2 i = ∫ p i, max p i, min d p ∫ 1 -1 d µfp 4 ( 1 -µ 2 ) 2 ∫ p i, max p i, min d p ∫ 1 -1 d µfp 4 µ 2 .  \nThe QLD effect on E CR ,i can be numerically computed from the right-hand side of Equation 25,  \nQLD on group i = ∫ p i, max p i, min dln p ∫ 1 -1 d µp 3 γ C 2 × ∂ ∂ ln p [ p 2 f × ( sgn( µ ) D pµ ∂ µ ln f + D pp ∂ ln p ln f p )] ,  \nwhere both f and the diffusion coefficients are directly measured in the simulations. By comparing the above numerical integration with the CR energy equation  \n(Equation 35), we can obtain ν eff ( p ),  \nν eff ( p i ) = QLD on group i U A C E CR ,i ( 1 2 | ξ 2 i -1 | + 4 U A 3 C ) , (38)  \nwhere p i should be understood as √ p i, min p i, max as the mean particle momentum in the i th bin. Here numerical differentiation are applied to the logarithms of f and p to mitigate potential truncation errors. Also, to further minimize errors, we perform one partial integration in computing 'QLD on group i '.  \nAfter normalized by environmental parameters, the pre-factors of ν eff ( p ) across all simulations closely follow a same curve (Figure 8), which peak around p 0 where resonant waves are the most pronounced. Following early discussions, we consider the scattering rates measured over range of (0 . 5 p 0 , 5 p 0 ) to be reliable. The curves of ν eff in the expanding box largely overlap with the theoretical QLT prediction (Equation 31). On the other hand, in the compressing box, the peak value of ν eff exceeds the QLT prediction by a factor of ∼ 2, achieved at a CR momentum of around 2 p 0 , and the scattering rate is generally higher than QLT prediction by a factor of 2-3 for particles with p ≳ p 0 . This difference is induced by the expanding/compressing box setup which contributes to additional wave damping/growth at the rate of ˙ a (which is much less than our adopted ν IN ), and also influence the peak momentum of CRs and the CR cyclotron frequency Ω. Moreover, in theoretical calculations (Section 3.1), we have assumed that the anisotropy level well exceeds U A /c , which does not entirely hold in simulations, particularly in the compressing box runs (Figure 7). In both expanding and compressing boxes, the main 'outlier' is when we double the parameter Ω 0 /ν IN , yielding pre-factors typically smaller than those from other runs by a factor of ≲ 2. We attribute this deviation to non-negligible numerical dissipation. 4 We have also verified that the results over the range of (0 . 5 p 0 , 5 p 0 ) largely converge in high-resolution runs, while the high-resolution runs yield larger values of ν eff towards low p , where their resonant shortwavelength waves are better resolved.", '4.3.2. Single-fluid results': 'More common in the formulation of CR hydrodynamics is to treat the CRs as a single fluid, effectively integrating over the momentum (energy). We determine ν eff  \nFigure 8. The effective scattering rate ν eff normalized by environmental parameters, depicted as a function of CR momentum (energy) p . The solid lines represent the compressing box runs, while the dashed lines denote the expanding box cases. The simulations with varying parameters and resolutions are distinguished by line colors, albeit with substantial overlap. The black dotted line illustrates the estimation of quasi-linear theory (QLT) on ν eff (Equation 31), in the momentum by momentum treatment.  \n<!-- image -->  \nfor the single fluid in the same way, by taking moments for Equation 25 and calculating the pressure anisotropy for the whole CR population (Equation 24). By default, we set p max = 10 p 0 and p min = 0 . 1 p 0 , which encompasses the majority of CR particles in the κ distribution. However, we also integrate CRs over a short range, (0 . 5 p 0 , 5 p 0 ), where CR ansiotropy reaches a quasi-steady state. We plot the CR anisotropy level and effective scattering normalized by the environmental parameters for all our simulation runs in Figure 9 & 10.  \nIn Figure 9, we see that the mean CR anisotropy | P CR, ⊥ /P CR, ∥ -1 | for (0 . 1 p 0 , 10 p 0 ) particles varies around 10 ( U A /c ) ( ν IN / Ω 0 ) ( ρ 0 /ρ CR ), and their corresponding effective scattering rate in Figure 10 varies around 0 . 07 ∣ ∣ ∣ ˙ B/B ∣ ∣ ∣ ( c/U A ) (Ω 0 /ν IN ) ( ρ CR /ρ 0 ). There is a factor of up to 4 difference between the expanding and compressing box cases, but both deviate substantially from single-fluid QLT predictions. However, in the compressing box cases, when averaged over the range (0 . 5 p 0 , 5 p 0 ), ν eff largely overlaps with the single-fluid QLT prediction, and the anisotropy level is only twice that predicted by the theory. In the expanding box cases, the results still deviate from the QLT predictions when averaged over the same range, but get closer than averaging over (0 . 1 p 0 , 10 p 0 ), with  \nFigure 9. The total pressure anisotropy level after normalized by environmental parameters, in the single fluid treatment. The hollow markers regards the majority of CRs, within 0 . 1 p 0 < p < 10 p 0 , as a CR single fluid, while the filled markers refer to the integration over the CRs which reaches steady pressure anisotropy (0 . 5 p 0 < p < 5 p 0 ). The red circles correspond to the results in the compressing box and the black squares indicate those in the expanding box. The grey dashed line refers to the theoretical prediction, | P CR, ⊥ /P CR, ∥ -1 | ∼ 1 . 78 ( U A /c ) ( ν IN / Ω 0 ) ( ρ 0 /ρ CR ) (Equation 23).  \n<!-- image -->  \n| P CR, ⊥ /P CR, ∥ -1 | ∼ 8 ( U A /c ) ( ν IN / Ω 0 ) ( ρ 0 /ρ CR ) and ν eff ∼ 0 . 09 ∣ ∣ ∣ ˙ B/B ∣ ∣ ∣ ( c/U A ) (Ω 0 /ν IN ) ( ρ CR /ρ 0 ). We can see that over any momentum range considered, the CR pressure anisotropy is primarily dominated by highenergy CRs due to the nature of second moments. However, the standard single-fluid QLT estimation relies on CRs around p 0 , by substituting Q 2 (Ω 0 m/p 0 ) for the whole CR population. The anisotropy level (and effective scattering rate) of high-energy CRs is significantly higher (and lower) than that of CRs around p 0 (Figures 4 & 7). As the average range gets closer to p 0 , both | P CR, ⊥ /P CR, ∥ -1 | and ν eff better converge to the value for p 0 in the p -byp treatment, approaching the theoretical predictions.  \nIn both Figure 9 and 10, the main outlier again is from the run with ν IN reduced by half. We again attribute this to the higher fraction of numerical wave damping as  \nFigure 10. The effective scattering rate ν eff after normalized by environmental parameters, in the single fluid treatment. The markers share the same meaning in Figure 9. The grey dashed line refers to the theoretical prediction, ν eff ≈ 0 . 187 ∣ ∣ ∣ ˙ B/B ∣ ∣ ∣ ( c/U A ) (Ω 0 /ν IN ) (Equation 31).  \n<!-- image -->  \ndiscussed earlier, and the results from this run would be brought closer to the typical values obtained from other runs should numerical dissipation be mitigated. Therefore, all our simulations generally exhibit consistent results over parameter variations, confirming the scaling from QLT. There are modest-to-strong deviations in the normalization factor, which depend on momentum range considered. A convergence in numerical resolution further validates our simulation results.  \nFinally, the total cooling/heating on the CR energy density, after plugging in ξ 2 and ν eff to Equation 35, yields  \n∼ -0 . 3 ˙ B B U A c E CR ,  \nin both expanding boxes and compressing box cases, twice than the QLT prediction (Equation 36). Therefore, its contribution to irreversible CR heating/cooling raises to about (1 / 8) U A /c (instead of (1 / 16) U A /c ) compared to adiabatic heating/cooling, again being very minor to the overall energy. We find this value is insensitive to variations in environmental parameters, nor to the range of momentum involved (within the range of  \n(0 . 1 p 0 , 10 p 0 ) studied here), thanks to the cancellation of such dependencies in ξ 2 and ν eff .', '5. DISCUSSION: APPLICATION OF CRPAI': "Given the dependency of ν eff on environmental parameters, we estimate the value of ν eff in realistic systems and compare the importance of CRPAI with that of the Cosmic Ray Streaming Instability (CRSI), which has been widely considered in modeling CR feedback in the galactic context. The effective scattering rate for the CRSI with ion-neutral damping scales with the environmental parameters as ∼ 0 . 054 ( c/L gal ) ( c/U A ) (Ω 0 /ν IN ) ( ρ CR /ρ 0 ) (Bai 2022), where L gal represents the typical height of the galactic disk scale, also the typical length scale of the CR gradient. Comparing the scaling relations (Equation 32), the ratio of ν eff between the pure CRPAI (Figure 10) and pure CRSI cases is 5 ∣ ∣ ∣ ˙ B/B ∣ ∣ ∣ ( L gal /c ).  \nIn a general sense, the compression/expansion timescale for the CRPAI can be estimated as the eddy turnover time in turbulence. 5 Let the largest eddy size in the ISM be L gal , and the eddy turnover timescale for Alfv'enic turbulence can be estimated as ( L gal /U A ) ∼ ∣ ∣ ∣ B/ ˙ B ∣ ∣ ∣ . Therefore, the CR scattering rate resulting from the CRPAI is smaller than that of CRSI by the order of U A /c . A similar conclusion was drawn by Zweibel (2020) who focused on CR heating. Setting typical values in the galactic context for the cold neutral medium (CNM) 6 , one finds ν eff for CRSI to be around once per year (e.g., Chan et al. 2019; Zweibel 2020; Hopkins et al. 2021), while the effective scattering rate of CRPAI is on the order of 10 -3 per year. Assuming only the CRPAI operates (no CR streaming), we find that the anisotropy level of ∼ GeV CRs to be as low as 10 -5 under such parameters.  \nThe CRPAI can be triggered more effectively by smaller eddies. Given an eddy size l , one can approximately write | δB | ∼ B 0 ( l/L gal ) 1 / 4 ∼ 1 / 3 (1 / 3 for the slow magnetosonic mode and 1 / 4 for the fast magnetosonic mode, see details in Goldreich & Sridhar 1995; Cho & Lazarian 2003), where B 0 represents the mean magnetic field strength. The corresponding variation timescale can be estimated as l/U A , with magnetic field change rate experienced by CRs given by ∣ ∣ ∣ ˙ B/B ∣ ∣ ∣ ∼  \n| δB/B 0 | / ( l/U A ) ∼ ( U A /L gal ) ( l/L gal ) -3 / 4 ∼-2 / 3 , which increases as the eddy size l decreases. Consider a modest scenario taking l ∼ 1pc while L gal ∼ 1kpc, the CR scattering rate resulting from the CRPAI may become comparable to that of CRSI in such cases. The transittime damping in turbulence (e.g. Fisk 1976; Eilek 1979; Achterberg 1981; Miller et al. 1996; Schlickeiser & Miller 1998), which primarily acts on the CR momenta aligned with the background magnetic field, can further enhance CRanisotropy and thereby increase the efficiency of CRPAI (Brunetti & Lazarian 2007).  \nAnother possible scenario where the CRPAI may actively manifest itself is around the mixing layer of multiphase gas (Begelman & Fabian 1990; Klein et al. 1994). Multiphase gas is ubiquitous in astrophysical environments such as the shell of bubbles from stellar explosions (e.g., El-Badry et al. 2019; Lancaster et al. 2021a,b, 2024), the circumgalactic medium (CGM) (e.g., Gronke & Oh 2018; Chen & Oh 2024), galactic winds (Fielding & Bryan 2022), and cosmic filaments (e.g., Mandelker et al. 2020). Mixing layers that separate and mix up different phases at the boundary are dynamic, turbulent and usually rapidly radiating energy away(e.g. Ji et al. 2019; Fielding et al. 2020; Tan et al. 2021), with the inflow gas experiencing cooling contraction and rapid magnetic field amplification (Zhao & Bai 2023; Das & Gronke 2024). While the dynamical effect of the CRs in the mixing layers is yet to be understood, one might approximately estimate that the magnetic field changing timescale | ˙ B/B | -1 as experienced by CRs co-moving with the gas can be as short as ∼ 0 . 01 Myr 7 . In comparison to the largest eddy turnover time in ISM, approximately 10Myr, the multiphase gas mixing may efficiently induce CR anisotropy, coupling the CRs through the CRPAI, potentially leading to a ν eff up to ∼ once per year. A comparison with the CRSI is less straight-  \nd since the level of CR gradient across the mixing layer is unclear, but as significant wave damping is expected in the cold phase, the CR distribution tends to be relatively flat across the phase boundaries (Armillotta et al. 2024).  \nAnother potential applicability of the CRPAI is the CR bottleneck (e.g., Skilling 1971; Wiener et al. 2017a), where CRs accumulate in front of a gas density bump (a dip in the CR streaming speed U A ), leading to a staircase-like CR pressure structure devoid of gradients on both sides of the bottleneck. Such 'meso-scale' structures could be ubiquitous and play a fundamental role in CR feedback (e.g., Bruggen & Scannapieco 2020; Bustard & Zweibel 2021; Tsung et al. 2022; Quataert et al. 2022), yet our exploration and understanding of such structures is still at very early stages. The CR streaming is expected to diminish on the two sides of the bottleneck, and they can be subject to the CRPAI in the presence of ˙ B/B . While it is expected that heating by CRPAI is not necessarily efficient (Bruggen & Scannapieco 2020), it is yet to explore how the additional coupling enabled by CRPAI affect the properties of the CR bottleneck, and towards global scales.", '6. SUMMARY AND OUTLOOK': "In this study, we employ kinetic (MHD-PIC) simulations to investigate the saturated state of CRPAI in the presence of ion-neutral damping. The CR anisotropy is continuously driven by the expanding box, which mimics gas expansion and compression in turbulence (and potentially shear motion) and is characterized by the magnetic field change rate (the expansion rate) ˙ B/B . Meanwhile, the ion-neutral damping is treated in the short-wavelength limit generally applicable in the neutral medium of the ISM, where the ionized gas motion is damped at a rate ν IN / 2. We achieve the saturated state self-consistently, where the growth rate of the CRPAI (Equation 20) is balanced by the ion-neutral damping rate, while the CR anisotropy level is maintained between driving (by expanding box) and quasi-linear diffusion (QLD, Equation 25). The kinetic simulation at the saturated state enables us to calibrate the scattering rate ν eff and the equilibrium pressure anisotropy level P CR , ⊥ /P CR , ∥ -1 from first principles, which are essential for understanding the role of CR feedback at macroscopic (e.g., galactic) scales. These properties scale with local environmental conditions (Equation 23 & 31), including the mass density ratio between CRs and thermal ionized gas ρ CR /ρ 0 , the Alfv'en speed U A , among others. We also vary simulation parameters to validate the scaling relation with environmental parameters.  \nOur main findings are as follows:  \n- · We detailed the formulation of the quasi-linear theory (QLT) under the aforementioned balancing relations, which yields the momentum-bymomentum estimates of the CR scattering rates (Equation 27) and anisotropy level (Equation 23) based on the CRPAI as a function of environmental parameters ( ˙ B/B , ν IN , etc.).\n- · We have verified that the CR pitch angle and momentum diffusion coefficients D µµ , D pp can be accurately reproduced from the wave intensities based on quasi-linear diffusion, lending support to the general validity of QLT.\n- · For those CRs that attain the saturated state (within the momentum range (0 . 5 p 0 , 5 p 0 )) in our simulations, the CR pressure anisotropy P CR , ⊥ /P CR , ∥ -1 largely follows the the momentum dependence predicted by QLT, with highenergy CRs exhibiting greater anisotropy (Figure 4). This pressure anisotropy also scales with the environmental parameters. However, the value of anisotropy in the single-fluid treatment exceeds the QLT prediction by a factor of 2 ∼ 4 (Figure 9).\n- · With the p -byp treatment, the measured CR scattering rates ν eff ( p ) (Figure 8) are consistent with the anticipated trend and the dependence on environmental parameters as predicted by QLT (Equation 27), except that the results with compressing box yields ν eff ( p ) that are higher by a factor of about 2 for p ≳ p 0 particles. In the single fluid treatment for CRs, the values of ν eff are consistent with QLT scaling, though with a factor of up to a few difference in the normalization depending on the range of momentum considered (Figure 10).\n- · The contribution of irreversible CR heating/cooling from the CRPAI is about a factor of 2 higher than QLT prediction (Section 3.1), ∼ -0 . 3 ( ˙ B/B ) ( U A /c ) E CR , which represents minor contributions compared to adiabatic heating/cooling.\n- · In a broader context, the scattering between CRs and MHD waves resulting from CRPAI is notably weaker compared to that of CRSI, by about a factor of U A /c on global scales. Nevertheless, we speculate that the CRPAI could play a more dominant role in smaller-scale structures, including the multiphase mixing layers and the CR bottleneck.  \nOur MHD-PIC simulations on the saturation of the CRPAI, together with that of the CRSI (Bai 2022), provide valuable insights into the microphysics about CR  \nfeedback from the first principles. The ultimate goal of such studies is to provide calibrated CR scattering rates from the CR gyro-resonant instabilities that serve as subgrid physics for CR (magneto-)hydrodynamic simulations at macroscopic scales. So far, we have studied the pure CRSI case driven by imposed CR pressure gradient, and the pure CRPAI case driven by ˙ B/B . The methodology laid out in these two works pave the way for us to further explore the most general case, where both a CR gradient and ˙ B/B are present, which may lead to marked differences from the cases studied so far. In the meantime, future work should explore multidimensional effects (Zeng et al. 2024), as well as the case with other wave damping mechanisms, especially the non-linear Landau damping that dominates in the more volume-filling warm medium (e.g., Armillotta et al. 2021). Finally, with the adaptive δf method, it becomes  \npossible to explore 'meso-scale' phenomena such as the CR bottleneck (e.g., Skilling 1971; Wiener et al. 2017a; Tsung et al. 2022; Quataert et al. 2022) fully kinetically, which would be a major effort to directly link CR feedback from microscopic to macroscopic scales.  \nWe thank Eve Ostriker, Eliot Quataert, Ellen Zweibel  \n1  \n- and Alexandre Lazarian for useful discussions, particu2\n- larly during the Aspen Center for Physics program on 3\n- 'Cosmic Ray Feedback in Galaxies and Galaxy Clusters'. 4\n- This work is supported by National Science Foundation 5\n- of China under grant No. 12325304, and in part by 6\n- the National Science Foundation under Grant No. NSF 7\n- PHY-2210452. Numerical simulations are conducted in 8  \n9  \nthe Orion clusters at Department of Astronomy, Ts-  \n- inghua University and TianHe-1 (A) at National Super10\n- computer Center in Tianjin, China. 11", 'REFERENCES': 'Achterberg, A. 1981, A&A, 97, 259  \nAdkins, T., & Schekochihin, A. A. 2018, Journal of Plasma  \nPhysics, 84, 905840107, doi: 10.1017/S0022377818000089  \nAdriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al.  \n2011, Science, 332, 69, doi: 10.1126/science.1199172  \n- -. 2013, ApJ, 765, 91, doi: 10.1088/0004-637X/765/2/91 Armillotta, L., Ostriker, E. C., & Jiang, Y.-F. 2021, ApJ, 922, 11, doi: 10.3847/1538-4357/ac1db2\n- -. 2022, ApJ, 929, 170, doi: 10.3847/1538-4357/ac5fa9 Armillotta, L., Ostriker, E. C., Kim, C.-G., & Jiang, Y.-F. 2024, arXiv e-prints, arXiv:2401.04169,  \ndoi: 10.48550/arXiv.2401.04169  \nBai, X.-N. 2022, ApJ, 928, 112,  \ndoi: 10.3847/1538-4357/ac56e1  \nBai, X.-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015, ApJ, 809, 55, doi: 10.1088/0004-637X/809/1/55  \nBai, X.-N., Ostriker, E. C., Plotnikov, I., & Stone, J. M. 2019, ApJ, 876, 60, doi: 10.3847/1538-4357/ab1648  \nBambic, C. J., Bai, X.-N., & Ostriker, E. C. 2021, ApJ, 920, 141, doi: 10.3847/1538-4357/ac0ce7  \nBegelman, M. C., & Fabian, A. C. 1990, MNRAS, 244, 26P  \nBerezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin,  \nV. S. 1990, Astrophysics of cosmic rays (North-Holland)  \nBirdsall, C. K., & Langdon, A. B. 2004, Plasma physics via  \ncomputer simulation (CRC press)  \nBlasi, P. 2013, A&A Rv, 21, 70,  \ndoi: 10.1007/s00159-013-0070-7  \nBlasi, P., Amato, E., & Serpico, P. D. 2012, PhRvL, 109, 061101, doi: 10.1103/PhysRevLett.109.061101  \nBoris, J., Shanny, R., of Naval Research, U. S. O., & Laboratory, N. R. 1972, Proceedings: Fourth Conference on Numerical Simulation of Plasmas, November 2, 3, 1970 (Naval Research Laboratory). https: //books.google.com.hk/books?id=zqxSAQAACAAJ  \nBott, A. F. A., Arzamasskiy, L., Kunz, M. W., Quataert, E., & Squire, J. 2021, ApJL, 922, L35,  \ndoi: 10.3847/2041-8213/ac37c2  \nBruggen, M., & Scannapieco, E. 2020, ApJ, 905, 19,  \ndoi: 10.3847/1538-4357/abc00f  \nBrunetti, G., & Lazarian, A. 2007, MNRAS, 378, 245, doi: 10.1111/j.1365-2966.2007.11771.x Bustard, C., & Zweibel, E. G. 2021, ApJ, 913, 106,  \ndoi: 10.3847/1538-4357/abf64c  \nButsky, I. S., & Quinn, T. R. 2018, ApJ, 868, 108, doi: 10.3847/1538-4357/aaeac2  \nChan, T. K., Kereˇs, D., Hopkins, P. F., et al. 2019, MNRAS, 488, 3716, doi: 10.1093/mnras/stz1895  \n- Chen, Z., & Oh, S. P. 2024, MNRAS, 530, 4032, doi: 10.1093/mnras/stae1113\n- Cho, J., & Lazarian, A. 2003, MNRAS, 345, 325, doi: 10.1046/j.1365-8711.2003.06941.x  \nConsolandi, C. 2014, arXiv e-prints, arXiv:1402.0467, doi: 10.48550/arXiv.1402.0467  \nDas, H. K., & Gronke, M. 2024, MNRAS, 527, 991, doi: 10.1093/mnras/stad3125  \nDashyan, G., & Dubois, Y. 2020, A&A, 638, A123, doi: 10.1051/0004-6361/201936339', 'A. BENCHMARK TESTS ON FOR THE ADAPTIVE δF METHOD': 'We conduct two benchmark tests to validate the accuracy and performance enhancement of the adaptive δf method (Section 2.2). These benchmark problems closely match the simulation setup in this study, wherein the expanding box continuously drives the CR anisotropy, and the fitting method for f 0 remains consistent with Equation 18 & 19. However, the parameters in the subsequent simulations are simplified for benchmark purposes and do not hold astrophysical meanings.', 'A.1. Accuracy test: fitting method for f 0 in the expanding box': 'We simulate the adiabatic evolution of the CR population in the expanding box to validate the fitting method outlined by Equation 18 & 19. In the absence of isotropization from waves and solely driven by box expansion/compression, each CR particle conserves pµa 2 ( t ) and pa ( t ) √ 1 -µ 2 (Equation 12, also see the Appendix of Sun & Bai 2023). Consequently, the CR distribution f , starting from an initial isotropic κ distribution (Equation 14), should precisely conform to the anisotropic κ distribution (Equation 15) with ξ = a ( t ). In this benchmark problem, we utilize Equations 18 and 19 to dynamically fit the CR anisotropy parameter ξ and the peak momentum p 0 . If the fitting method performs correctly, the fitted ξ should equal to a ( t ), while the fitted p 0 should remain constant.  \nThe simulation setup largely follows the fiducial runs in the expanding/compressing boxes, albeit slightly simplified for test purposes. We initialize CRs in a uniform isotropic κ distribution (Equation 14) with p 0 = 200 U A and κ = 1 . 25. The simulation particles are segregated into eight momentum bins, with 2 . 4576 × 10 6 particles allocated to each bin, matching the total number of simulation particles in the fiducial runs (Section 2.3). To prevent isotropization, we set the initial wave amplitudes to zero and disable the CR backreaction. The background field and the expanding box setup are identical to those in Section 2.3, while the expansion/compression rate is elevated to ˙ a = ± 10 -2 Ω 0 . The simulation spans 10 2 Ω -1 0 , where the CR anisotropy level becomes significant by the end of the simulation, with f 0 updated every ∆ T adapt = 1Ω -1 0 .  \nThe fitting results for ξ and p 0 are depicted in Figure 11. The fitted ξ closely tracks the expansion rate a ( t ) = exp(˙ at ), while the fitted p 0 remains constant around the initial value of 200 U A , in both the expanding box and the compressing box. We note that there is a systematic deviation in the fitted p 0 less than 1% relative to the anticipated values. This small deviation is related to the initialized CR distribution, where initial simulation particles only span over a finite range of the κ distribution ( p 0 / 500 to 500 p 0 , see Section 2.3) 8 . The fitting error can be further minimized by covering over a larger range of the κ distribution or increasing the fitting frequency (decreasing ∆ T adapt ˙ a ).', 'A.2. Performance test: signal-to-noise improvement in CPPAI': 'We compare the wave spectra triggered by CRPAI in the expanding box, both with the adaptive δf method or the traditional δf method, to assess the effectiveness. The adaptive δf method can alleviate the statistical noise in the CR backreaction, particularly in scenarios where the CR distribution f is less predictable.  \nThe simulation setup closely resembles that in Section 2.3, albeit with simplifications tailored for testing purposes. In comparison to the fiducial runs, we eliminate the ion-neutral damping and augment the box expansion rate to ˙ a = +5 × 10 -5 Ω 0 , to expedite both the wave isotropization and the anisotropy driving. Additionally, for a significant wave intensity within a limited runtime, we set the initial wave amplitude as A = 10 -2 and enhance the CR mass ratio to ρ CR /ρ 0 = 10 -4 . We set the 1D simulation box size to L x = 9 . 6 × 10 4 U A / Ω 0 and the resolution to 10 U A / Ω 0 . All other parameters and setups, such as the initial CR distribution and the number of momentum bins for initial simulation particles, remain identical to those of the fiducial run in the expanding box.  \nFigure 11. Fitting results for the CR anisotropy parameter ξ and the peak momentum p 0 in the test of the adaptive δf method using expanding box (blue) and compressing box (red). The box expansion rate ˙ a is ± 10 -2 Ω 0 . The markers in the top two panels represent the fitting results obtained from Equation 18 and Equation 19, while the solid lines depict the expected values. The bottom panel illustrates the relative deviation of the fitted ξ from the expected values.  \n<!-- image -->  \nWe carry out three parallel runs by varying the choice f 0 in the δf method and particles per cell per momentum bin: the adaptive f 0 with four particles per cell (ppc) per momentum bin, the adiabatic f 0 with four ppc per momentum bin, and the adiabatic f 0 with 16 ppc per momentum bin. The adaptive f 0 fits ξ and p 0 every ∆ T adapt = 5Ω -1 0 (see Section 2.2). The adiabatic f 0 maintains ξ = a ( t ) and a constant p 0 , representing CR adiabatic evolution without wave isotropization. We illustrate the wave spectra measured at t = 3 × 10 4 Ω -1 0 in Figure 12. Comparing the spectra from the two adiabatic f 0 runs with different particles per cell, the statistical noise from the CR backreaction primarily affects the highk regime, and the run with 16ppc has significantly smaller noise level. The adaptive δf method even  \nmore effectively reduce the noise at highk . With only 4ppc, the noise level is even smaller than the run with adiabatic f 0 but using four times less particles.  \nThe effectiveness of the adaptive δf method depends on the specific problem at hand. In this benchmark, it significantly reduces computational costs by over four times compared to the traditional δf method, where f 0 remains the one evolving adiabatically from the initial CR distribution. However, when f is predictable, such as when isotropization from waves dominates, or when quasi-linear diffusion is negligible and CR evolution is adiabatic, the adaptive δf method offers no signifcant improvement in signal-to-noise ratio. In our study, the saturation of CRPAI, where f deviates from both isotropic distribution and adiabatic distribution, the adaptive δf method proves essential. Note that the anisotropy level in the actual simulations varies with p , and further improvement may be considered by allowing ξ to be a parameterized function of p , but this can be left for future work.  \nFigure 12. Wave intensity spectra, kI ( k ), in the benchmark problem with varying f 0 and numbers of simulation particles, at time t = 3 × 10 4 Ω -1 0 . The color and line styles correspond to those in Figure 2. Lines of the same color (representing polarization) substantially overlap.  \n<!-- image -->'}