Turbulent Reacceleration of Streaming Cosmic Rays

Subsonic, compressive turbulence transfers energy to cosmic rays (CRs), a process known as non-resonant reacceleration. It is often invoked to explain observed ratios of primary to secondary CRs at $\sim \rm GeV$ energies, assuming wholly diffusive CR transport. However, such estimates ignore the impact of CR self-confinement and streaming. We study these issues in stirring box magnetohydrodynamic (MHD) simulations using Athena++, with field-aligned diffusive and streaming CR transport. For diffusion only, we find CR reacceleration rates in good agreement with analytic predictions. When streaming is included, reacceleration rates depend on plasma $\beta$. Due to streaming-modified phase shifts between CR and gas variables, they are slower than canonical reacceleration rates in low-$\beta$ environments like the interstellar medium (ISM) but remain unchanged in high-$\beta$ environments like the intracluster medium (ICM). We also quantify the streaming energy loss rate in our simulations. For sub-Alfv\'{e}nic turbulence, it is resolution-dependent (hence unconverged in large scale simulations) and heavily suppressed -- by an order of magnitude -- compared to the isotropic loss rate $v_{A} \cdot \nabla P_{\rm CR} / P_{\rm CR} \sim v_{A}/L_{0}$, due to misalignment between the mean field and isotropic CR gradients. Counterintuitively, and unlike acceleration efficiencies, CR losses are almost independent of magnetic field strength over $\beta \sim 1-100$ and are, therefore, not the primary factor behind lower acceleration rates when streaming is included. While this paper is primarily concerned with how turbulence affects CRs, in a follow-up paper (Bustard and Oh, in prep), we consider how CRs affect turbulence by diverting energy from the MHD cascade, altering the pathway to gas heating and steepening the turbulent power spectrum.


INTRODUCTION
Cosmic rays (CRs) are an important non-thermal component of galaxies and their surroundings. In the Milky Way interstellar medium (ISM), their collective energy density is comparable to that in thermal gas, magnetic fields, and turbulence (Boulares & Cox 1990), and they are believed to play a significant role in ionizing molecular clouds (e.g., Dalgarno 2006), driving galactic outflows (Ipavich 1975;Breitschwerdt et al. 1991;Salem & Bryan 2014;Ruszkowski et al. 2017;Buck et al. 2020;Bustard et al. 2020;Hopkins et al. 2021), and impacting pressure balance, stability, and fragmentation of the ISM (Parker 1966;Heintz et al. 2020). Their long residence times (compared to the light crossing time) in the disk, inferred from ratios of their spallation products, suggest that CRs are frequently scattered by magnetic perturbations, forcing them to undergo a random walk along magnetic field lines.
Such frequent scatterings drive the galactic CR distribution to almost perfect isotropy, making it difficult to connect CRs back to their sources; however, detections of their secondary by-products give us clues to their origins and transport throughout the Galaxy. The "standard paradigm", implemented in most phenomenological propagation models (e.g. Strong & Moskalenko 1998;Evoli et al. 2017), is one in which galactic CRs are predominantly created by diffusive shock acceleration (DSA) at supernova remnant shock fronts, followed by energy-dependent, diffusive propagation. This diffusive picture is physically motivated for high-energy CRs above a few hundred GeV, for which we believe CRs are confined by scattering off hydromagnetic waves created by an external turbulent cascade (Chandran 2000;Yan & Lazarian 2004).
For E 300 GeV, however, the dominant confinement mode appears to be self-confinement, where CRs pitch angle scatter off Alfvén waves that the CRs generate themselves through a resonant streaming instability (Wentzel 1968;Kulsrud & Pearce 1969). In this energy range, the resulting transport is a mixture of fieldaligned streaming down the local CR pressure gradient at the Alfvén speed v A,i = B/ √ 4πρ i and field-aligned diffusion that arises from wave damping (Section 2.2) and the subsequently reduced pitch angle scattering rate (see recent reviews by e.g. Zweibel 2013Zweibel , 2017Amato & Blasi 2018). Self-confinement is not only supported by a growing number of theoretical studies but also by a break in the primary spectrum at 300 GeV, which aligns with the theoretically expected transition from self-confinement to extrinsic turbulence (Blasi et al. 2012;Aguilar et al. 2015;Evoli et al. 2018;Amato & Blasi 2018;Evoli et al. 2019;Kempski & Quataert 2022). While this model is increasingly implemented in both galaxy evolution and CR propagation models, there are still many unknown outcomes of diffusive vs streaming transport, both for the interpretation of observables and for our understanding of CR influence on galaxy evolution. The one we concern ourselves with in this paper is turbulent reacceleration.
CRs in a turbulent background can gain energy from both resonant and non-resonant interactions through second order Fermi acceleration (Brunetti & Lazarian 2007;Lynn et al. 2013). The particle mean free path λ mfp relative to the lengthscale l of the turbulent eddy it interacts with determines the relevant regime. In the low scattering rate limit (λ mfp >> l), the key interactions are resonances between particles and magnetic compressions due to fast or slow magnetosonic waves 1 . Although strong damping generally precludes the turbulent cascade from extending down to the CR gyro-scale, CRs with parallel velocities similar to the wave speed can extract energy from it, a process known as transit-time damping (similar to Landau damping). Thus, weakly scattered CRs undergo resonant acceleration.
However, we have argued that low-energy E < 300 GeV CRs in our Galaxy are self-confined, with small mean free paths. This is particularly true of the ∼GeV CRs, which carry most of the energy and have λ mfp ∼pc scale mean free paths in ∼ µG magnetic fields. Since they are strongly scattered (λ mfp << l), the strong coupling with the gas means that they undergo compression and rarefaction with the thermal fluid. If the CRs were completely locked to the fluid, these constitute adiabatic reversible processes, and CRs experience no net energy gain. However, CR diffusion out of overdense regions breaks this symmetry, so that there is net energy transfer from the gas to the CRs (Ptuskin 1988). Heuristically, the CRs gain energy during compression, but diffuse out without giving this energy back. This non-resonant reacceleration is the focus of this paper. In Galactic propagation models, reacceleration is typically included as a default, but it is in increasing tension with data, including both secondaryto-primary ratios as a function of energy (Strong et al. 2007;Gabici et al. 2019), and a growing wealth of multiwavelength data, specifically synchrotron measurements (Trotta et al. 2011;Di Bernardo et al. 2013;Orlando & Strong 2013;Gabici et al. 2019).
The primary goal of this paper is to evaluate the effect of Alfvénic CR streaming, which is characteristic of all self-confined CRs. Since CRs only stream out of local CR maxima, CR streaming also breaks the symmetry between compression and rarefaction, but its effects are far less known. The impact of CR streaming on turbulent reacceleration has only recently been discussed analytically (Hopkins et al. 2022) and never been quantitatively studied in MHD simulations. We find that CR streaming reduces the efficiency of turbulent reacceleration in low β environments, potentially explaining why significant reacceleration appears to be disfavored in models. Somewhat counterintuitively, though, CR streaming losses are not responsible for lower acceleration rates: we don't find a correlation between CR streaming loss rates and reaceleration rates. Interestingly, loss rates are quite different from isotropic loss rates v A · ∇P CR /P CR ∼ v A /L 0 , due to misalignment between the mean field and isotropic CR gradients, and this effect is resolution dependent and may be incorrectly captured in low-resolution galaxy evolution simulations.
Our paper is outlined as follows: First, we provide background on turbulent reacceleration and present the previously-derived growth rate for purely diffusive CRs subject to long-wavelength, subsonic, isothermal turbulence ( §2). We then extend this treatment to include additional self-confinement (streaming) terms, and we derive a simple modification to the canonical reacceleration rate of purely diffusive CRs. In §3, we numerically validate these growth rates using a two-fluid CR implementation in the Athena++ code, and in §4 we quantify the efficiency of streaming energy loss in turbulence with varying plasma β. In §5, we discuss the implications of these results for studies of the ISM, intracluster medium (ICM), and circumgalactic medium (CGM), aimed at both the galaxy evolution and CR propagation communities. We conclude in §6.
While in this paper we focus on how turbulence affects CRs, in a follow-up paper (Bustard and Oh, in prep), we study the back-reaction of CRs on turbulence. This has largely been neglected, although it is clear that CRs are absorbing energy from the flow. We study this both analytically and in exploratory simulations. We find that in many physically plausible scenarios relevant to the ISM and CGM (where CR energy densities are expected to be significant), CRs can absorb a large fraction or even most of the turbulent driving energy at large scales, significantly steepening the turbulent energy spectrum or even wiping out small scale compressive motions.

BACKGROUND AND ANALYTIC ARGUMENTS
Following Ptuskin (1988), let's define three scales: L 0 is the outer eddy scale set by large-scale driving motions; L 1 is some smaller "cut-off" scale, which could be a viscous or damping scale, though here we will associate it with the width of a shock front in the medium 2 ; and L turb ∈ [L 1 , L 0 ] represents a general eddy scale within the turbulent cascade. Ptuskin (1988) derived the nonresonant acceleration rate from an ensemble of random acoustic waves and weak shocks propagating in compressive, subsonic 3 turbulence. We will quickly summarize those results in two illuminating limits (see also e.g. Lynn et al. 2013), first assuming as in Ptuskin (1988) that CR transport is purely diffusive, before considering the effects of streaming.

Non-Resonant Reacceleration of Purely Diffusive Cosmic Rays
Let us first consider acceleration operating on a single scale, the turbulent outer scale L 0 , with turbulent velocity v and CR diffusion coefficient κ. An important insight from Ptuskin (1988) is that in subsonic turbulence, the lifetime of a compression in subsonic turbulence is not the eddy turnover time τ eddy ∼ L 0 /v but instead the sound (or compressive wave) crossing time τ sc ∼ L 0 /v ph , which is shorter than τ eddy . Here, v ph ∼ (P tot /ρ) 1/2 ∼ [P g + P B + P CR )/ρ] 1/2 is the compressible wave phase velocity, and is simply given by the gas sound speed c s ∼ (P gas /ρ) 1/2 if P g P B , P CR . The two relevant timescales are therefore the sound crossing 2 Even in subsonic turbulence, non-linear steepening will still produce weak shocks. 3 For transonic or supersonic turbulence, strong shocks can additionally accelerate CRs through the first order Fermi mechanism, but then the resulting acceleration rate and CR spectrum are qualitatively different from second order Fermi acceleration.
time τ sc ∼ L 0 /v ph , and the diffusion time t diff ∼ L 2 0 /κ across an eddy.
First, consider the case where t diff t sc , i.e. κ >> v ph L 0 .
In this regime, quickly diffusing particles see an effectively static velocity field, and the derivation of the momentum diffusion coefficient, D pp , for a CR with momentum p follows that of the standard second order Fermi argument: D pp ∼ (∆p) 2 /τ scatter ∼ p 2 v 2 /(c 2 τ scatter ) ∼ p 2 v 2 /κ. The energy growth time, accurate to within a factor of a few, is then In the opposite limit, for which κ << v ph L 0 , CRs advect with the fluid and behave quasi-adiabatically during each compression and rarefaction: δp/p ∼ δP CR /P CR ∼ δρ/ρ ∼ v/v ph . Through successive compressions and expansions in a turbulent eddy of size L 0 , the CRs undergo a random walk in momentum space until they manage to diffuse out of the eddy on a timescale L 2 0 /κ. In this case, . Within a factor of a few, the growth time in this limit is Consider κ ∼ v ph L 0 , i.e. the case of maximally efficient acceleration. In this case, both Equations 1 and 2 reduce Consider now the contribution of multiple eddies in a turbulent cascade. If we sum contributions from a variety of scales, the general expression for the momentum diffusion coefficient is (Ptuskin 1988): The alert reader will notice that this minimum growth time (Equation 3) coincides with the cascade time of compressible fast modes of scale l in MHD turbulence, τcas ∼ lv ph /v 2 l (Nazarenko & Schekochihin 2011). This arises because both Fermi II acceleration and fast mode cascade timescales can be understood as the outcome of stochastic random walks. In Fermi II acceleration, tgrow ∼ τscatter(∆E/E) −2 ∼ τscatter(v/c) −2 ∼ v ph L/v 2 , assuming that τscatter ∼ κ/c 2 ∼ v ph L/c 2 . In fast modes, two wave packets of scale l interact on a wave crossing time τ ph ∼ l/v ph τ NL ∼ l/v l , where the non-linear steepening time τ NL ∼ l/v is set by the non-linear advection term v · ∇v in the Euler equation. Thus, each interaction results in a small velocity change ∆v/v ∼ τ ph /τ NL ∼ v/v ph . The cascade time is then set by the number of interactions required for ∆v ∼ v, i.e.
where W 1D (k) is the 1D turbulent power spectrum, normalized such that: What is an appropriate expression for W 1D (k)? Recall that for hydrodynamic turbulence, a standard Hodge-Helmholtz decomposition usually shows that the compressive component of the velocity field is Burgers-  Ptuskin (1988) for the full form of the power-law W 1D (k), including prefactors), and using t grow ∼ p 2 /D pp gives, in the regime v ph L 1 << κ << v ph L 0 (Ptuskin 1988): At scales outside of the power-law spectrum, i.e. at scales larger than L 0 or smaller than the cut-off scale L 1 , the reacceleration time reverts to Equations 1 and 2, respectively. The minimum reacceleration time of Equation 6 is t grow ∼ (9/2)v ph L 0 /v 2 , motivating one to write This is very similar to the resonant reacceleration rate in balanced turbulence (Skilling 1975;Thornbury & Drury 2014;Zweibel 2017): D pp ∼ 1/9(p 2 v 2 A /κ), which is essentially what is implemented in Galprop (Strong & Moskalenko 1998) and other CR propagation codes. The diffusion coefficient is commonly taken to be a power-law in particle rigidity κ ∝ R δ , where the diffusion powerlaw index, δ, is related to the wave spectrum power-law index. In that case, there are additional pre-factors related to the wave energy and power-law index.
Note that for W 1D (k) ∝ k −2 , the integrand in Equation 4 is proportional to 1/(v 2 ph + κ 2 k 2 ), i.e. it diminishes rapidly at high k. Unlike resonant acceleration, the smallest scales do not dominate acceleration. This is good news, as it implies that non-resonant acceleration can be captured in MHD simulations. However, there is still resolution dependence in the κ v ph L 0 regime (see Appendix): as resolution increases, we find a slow increase in acceleration rates in our simulations. In fact, as seen for the analytic L 0 /L 1 = 1000 case in Fig 7, Equation 6 gives a growth time which is almost constant for κ < L 0 v ph , i.e. it scales much more weakly with κ than t grow ∝ 1/κ, as obtained in the single eddy approximation (Equation 2). Heuristically, we can understand this as follows. For κ < L 0 v ph , the acceleration is maximized at some smaller scale l, where κ ∼ lv ph . However, for Burgers turbulence, the velocity at scale l is v l ∝ l 1/2 , so that the minimum growth time t grow ∼ v ph l/v 2 l (when κ ∼ lv ph ) is independent of l. In other words, as long as there is sufficient dynamic range, acceleration in Burgers turbulence is self-similar, with eddies at some scale l ∼ κ/v ph providing the dominant contribution, and other eddies in the κ v ph l or κ v ph l regimes providing subdominant contributions, with the growth rate roughly independent of l, and hence of κ.
Another possibility is that the compressive component has a Kraichnan spectrum W 1D (k) ∝ k −3/2 , as might be expected for compressive fast modes (Cho & Lazarian 2003). The case of resonant acceleration with a Kraichnan spectrum has been studied analytically, and it is a leading model for explaining giant radio halos in galaxy clusters (Brunetti & Lazarian 2007;Miniati 2015). However, our simulations in this paper have very limited dynamical range and do not resolve the fast mode cascade (which in any case is still very uncertain; see e.g. Kowal & Lazarian (2010), who find a Burgers-like fast mode spectrum). We therefore will compare our simulations to Equation 6 appropriate for a Burgers-like spectrum.
Equations 1 -6 are derived assuming isotropic CR diffusion, but cross-field diffusion is much less efficient than field-aligned diffusion. Accounting for anisotropic diffusion changes the above equations (Chandran & Maron 2004;Lynn et al. 2013) in the κ >> v ph L 0 regime (Equation 1), because CRs undergo a 1D random walk along magnetic field lines and are likely to return to the same eddy multiple times before the eddy is randomized in the turbulent flow. The characteristic timescale for turbulence-particle interactions is no longer the diffusion time across an eddy. Rather, it is the decorrelation time t corr ∼ L 0 /v of the eddy, which is longer: t corr /t diff ∼ κ /(vL 0 ) 1, when κ v ph L 0 . Over this time period, a typical CR can diffuse a distance L B ∼ κ t corr , over which time it will interact with N ∼ L B /L 0 ∼ κ || t corr /L 0 ∼ κ || /vL 0 eddies. Since over a period t corr ∼ (κ /vL 0 )t diff ∼ N 2 t diff the CR spends time in N eddies, on average it spends t eddy ∼ N t diff per eddy, i.e. it scatters with each eddy N times, so that it acquires an rms momentum boost (δp) eddy ∼ N (δp) iso . Thus, the momentum diffusion coefficient D pp ∼ (δp) 2 eddy /t eddy ∼ N 2 (δp) iso /(N t diffuse ) ∼ N D pp,iso , i.e. the reacceleration rate increases by a fac-tor of N compared to the momentum diffusion coefficient D pp,iso that assumes isotropic spatial diffusion. Effectively, anisotropic diffusion increases the coherence of CR-turbulent scattering (Chandran & Maron 2004): instead of momentum changes of order (δp) iso on a timescale t diff , the CR makes larger momentum changes (δp) eddy ∼ N (δp) iso on longer timescales t eddy ∼ N t diff . This is similar to reducing the opacity in radiative transfer, so that a photon has a longer mean free path and mean free time. Just as this reduces the photon escape time, the CR reacceleration time is reduced for a longer mean path in momentum space: t grow ∼ p 2 /D pp ∝ N −1 , i.e.: This means that acceleration can still be reasonably efficient off the "sweet spot" κ ∼ v ph L 0 for κ v ph L 0 , since t grow ∝ κ 1/2 instead of t grow ∝ κ. On the other hand, there is little difference in acceleration rates between anisotropic or isotropic diffusion in the κ << v ph L 0 regime. To summarize: in the single eddy limit, assuming isotropic scattering, acceleration rates have a reasonably sharp peak at κ ∼ v ph L 0 , where t grow ∼ v ph L 2 /v 2 is minimized, and t grow ∝ κ, t grow ∝ κ −1 in the κ/v ph L 0 1, κ/v ph L 0 1 regimes respectively. However, this dependence on κ is modified by two effects. Firstly, a hierarchy of eddies contribute to reacceleration, and for κ v ph L 0 , there exists a smaller eddy of scale l, such that κ ∼ v l l. For Burger's turbulence with sufficient dynamic range, t grow becomes almost independent 5 , and hence we often still make use of the 'single eddy' formulae when comparing simulations with analytics. But in reality (or in simulations with much higher dynamic range), we expect a much milder dependence of acceleration rates on κ. of κ. Secondly, anisotropic diffusion increases the coherence of CR acceleration, so that t grow ∝ κ 1/2 increases more weakly with κ in the κ v ph L 0 regime. The net result is a characteristic minimum growth time t grow ∼ v ph L 0 /v 2 which depends only weakly on κ, in a range of several dex around the 'sweet spot' κ ∼ v ph L 0 .

Non-Resonant Reacceleration of Streaming
Cosmic Rays 5 It is important to appreciate that in our simulations, the limited dynamic range means that this effect does not really kick inacceleration rates are closer to the 'single eddy' approximation, with tgrow ∝ κ −1 in the κ v ph L 0 regime While CRs in the ∼ GeV energy range are close to isotropic, even a small amount of drift anisotropy can excite magnetic perturbations through the resonant streaming instability (Wentzel 1968;Kulsrud & Pearce 1969;Skilling 1971). When the instability acts, forward traveling (relative to the CR drift) Alfvén waves are most efficiently excited, with backward waves quickly damped. In the absence of wave damping, CRs pitch angle scatter off these waves until they isotropize in the wave frame, and thus advect along the magnetic field at the local Alfvén speed. In the presence of wave damping, however, a steady-state balance between growth and damping gives a finite scattering rate which dictates that CRs diffuse relative to the wave frame. Since the streaming instability growth rate declines with increasing CR momenta, the diffusivity is also energy-dependent and generally rises with increasing CR energy. The resulting fluid CR transport is then a mixture of diffusion and streaming, as long as scattering by streaming generated waves dominates over scattering by extrinsic turbulence, a cut-off predicted to occur around 300 GeV (Blasi et al. 2012;Kempski & Quataert 2022), coincident with an observed change in the proton spectral index ). For energies above 300 GeV, CR propagation is thought to be purely diffusive.
While diffusion and streaming are similar in that they unlock CRs from the gas, the unique behavior of streaming fundamentally alters CR-wave interactions and, therefore, the reacceleration rate for CRs of any energy E ∼ < 300 GeV. Let's consider a single compression. In the no transport case where CRs are perfectly locked to the gas, all energy gained via compression will be lost via rarefaction. Crucially, CR diffusion introduces a π/2 phase shift between CR pressure and gas density perturbations. This "drag" against CRs provides a frictional force on compressive motions, giving rise to an acceleration that damps the wave. This phenomenon, akin to a damped simple harmonic oscillator, is known as Ptuskin damping (Ptuskin 1981) and is responsible for the nonresonant transfer of wave energy to CR energy (see §3.5 for additional discussion).
The π/2 phase shift between CR and gas perturbations requires a CR flux F CR ∝ ∇P CR , such that the CR restoring force in response to perturbations is proportional to velocity (e.g., see §2.1.1 of Tsung et al. (2022) for details), which is the requirement for either damping or driving. On the other hand, the CR flux with pure streaming is F CR ∝ P CR , resulting in a CR restoring force proportional to displacement and therefore unable to create a π/2 phase shift necessary for Ptuskin damping and reacceleration -as in a simple harmonic oscillator, energy is conserved. With pure streaming and no diffusion, there is then no energy transfer from waves to CRs.
Another way to understand this is in terms of phase space transport. Spatial diffusion is tied to momentum diffusion (as can be seen from the relationship between the spatial (κ) and momentum (D pp ) diffusion coefficients, e.g. in equation 4), while spatial advection is tied to momentum advection (as is true, for instance, in adiabatic compression or expansion). Since there is more phase space at higher momenta, CRs diffusing in phase space inevitably diffuse to higher energy. By contrast, in CR advection, there is no stochasticity: CR evolution is adiabatic in the Alfvén wave frame. This can be seen in the equation for the distribution function f (x, p) (Skilling 1971): where w = v + v A is the net velocity of Alfvén waves in the lab frame. Since ∇ · w = 0, i.e. there is no net converging or diverging flow for Alfvén waves over a box which is homogeneous on large scales, Df /Dt = 0and a time-stationary distribution function means no acceleration is taking place. All energy changes are reversible, since the energy gained due to converging Alfvén waves during compression is returned to the plasma as CRs stream out of the overdensity, resulting in diverging Alfvén waves. Advection only produces a net energy change when ∇ · w = 0 -for instance, when there is net compression or expansion of the fluid ∇ · v = 0, or there is a net change in Alfven speed, due to net gradients in B-field strength or gas density ∇ · v A = 0. The latter happens, for instance, when CRs stream in a stratified medium, resulting in CR energy loss at a rate v A · ∇P c . So what is the role of streaming in CR-wave interactions? Streaming introduces two additional effects. First, the advective transport modifies the perturbed acceleration from the phase-shifted CR force. One can derive the acceleration rateu in the diffusion only and then streaming-modified cases (Begelman & Zweibel 1994;Tsung et al. 2022): where u is the velocity, c c ∼ (P CR /ρ) 1/2 , and Γ ∼u/u is the damping/growth rate of gas motions. Note that only the streaming-dominated case allows for growth (Γ > 0). If there were a background CR gradient much larger than the perturbed gradients due to gas motions, growth can occur in the diffusion-only case when the diffusion time is shorter than the sound crossing time across a CR scale length L c (Drury & Falle 1986); however, that is not the case for the unstratified medium we assume here. In the streaming case, though, unstable growth is possible, regardless of L c , if β 1. This case was studied in detail in Tsung et al. (2022). For now, we will operate in the β ≥ 1 regime. Note also that Equation 11, as written, assumes an isothermal equation of state where gas heating by CRs is neglected. In non-isothermal gas, Equation 11 is multiplied by 1 ± (γ g − 1)v A /c s (Begelman & Zweibel 1994).
The ± in Equation 11 refers to gas motion in same (+) or opposite (-) direction relative to the direction of cosmic ray drift. Overdensities are created by converging gas flows, but CRs stream out of them in the opposite direction. These opposing flows reduce the net amount of CR compression in the Alfven wave frame, and hence reduce CR energization. Thus, the '-' sign is appropriate, and CRs take energy from gas motions at a new, reduced rate modified by the factor 1 − v A cs . The new CR reacceleration time relative to the pure diffusion growth time, obtained by dividing Equation 10 by Equation 11, is then Note that the growth rate now depends on plasma β, which for an isothermal (γ g = 1) gas is β = 2c 2 s /v 2 A . For an adiabatic γ g = 5/3 gas, the growth rate is increased by a factor of (1+2v A /3c s ), still resulting in a significant decrease in reacceleration when v A ∼ c s . As β increases, the rate converges to the pure diffusion rate. Note also that Equation 12 is only appropriate in the regime Second, wave excitation by the streaming instability drains energy from the CR population at a rate H = v A · ∇P CR (Zweibel 2017). This energy is subsequently transferred to the gas by wave damping. For an isothermal equation of state, the gas does not gain energy, but the energy sink for the CRs remains and competes with compressional heating to decrease the CR reacceleration rate. As we'll see in §4, however, reacceleration rates are not predominantly stunted by streaming energy loss; instead, when reacceleration rates are slowest (the low-β regime), CR energy loss rates are also slowest. The major correction to growth rates appears to come from the modified phase-shifted CR force.
3. MHD SIMULATIONS Motivated by the above considerations, the agenda for numerical simulations is straightforward: • Study growth times in the pure diffusion case as a function of κ, and compare to Ptuskin's predictions (Equations 1-3, 6).
• Study the effect of anisotropic diffusion, which is expected to reduce t grow → t grow /( κ/vL 0 ) in the κ v ph L 0 regime.
• Study the combined effect of CR streaming and diffusion on turbulent reacceleration, and check our new analytic expectations (Equations 12, 13), particularly as a function of plasma β. When streaming is important, we expect acceleration to be inefficient at low β.
• Quantify the fractional turbulent dissipation going into CRs, and the non-linear saturation of CR reacceleration.
• Study the resolution dependence of both CR reacceleration and streaming energy loss. Are these well captured in zoom-in cosmological or even fully cosmological simulations of galaxy formation? We show additional convergence tests and discuss implications in Appendix A.

Computational Methods
We ran a suite of simulations using a version of the magnetohydrodynamics (MHD) code Athena++ (Stone et al. 2020) modified to include CRs (Jiang & Oh 2018). CR transport via diffusion and streaming is implemented via a two-moment method developed originally for radiative transfer, and the efficiency and accuracy of this implementation have been extensively tested (Jiang & Oh 2018).
The equations solved are a combination of the ideal MHD equations and CR evolution equations for a mixture of streaming and diffusive transport: where ρ is the gas density, v g is the gas velocity, v st = −v A (B·∇P CR )/|B·∇P CR | is the CR streaming velocity, B is the magnetic field, E is the thermal energy, E CR is the CR energy, F CR is the CR energy flux, and, P B = B 2 /8π, P CR = (γ CR − 1)E CR are the magnetic and CR pressures.
Note that the CR adiabatic index is γ CR = 4/3, and we use an isothermal equation of state, for which the gas adiabatic index is γ g = 1. While gas can be effectively isothermal when radiative cooling is strong, this isothermality assumption is not rigorously appropriate. However, since the reacceleration time for a fixed κ depends on the phase velocity of compressible waves, i.e. the sound speed, we enforce γ g = 1 to keep a constant phase velocity and facilitate an easier comparison to analytic expectations. Our conclusions and applicability to the ISM, CGM, and ICM are quite insensitive to this choice: note that the minimum growth time from Equation 6 scales as v ph ∝ √ γ g and, hence, changes by the small factor of 5/3 modulating between an isothermal and adiabatic equation of state. An alternative implementation, with γ g = 5/3 and including radiative cooling, requires an additional heating term tuned to enforce global thermal balance (rather than allowing for a cooling runaway), but this heating term is often invoked as a substitute or rough parametrization of CR heating itself. As a natural starting place for our study of cosmic ray -turbulence interplay, the isothermality assumption is cleaner. The two-moment method presents itself through the inclusion of a maximum speed of light parameter v m and an interaction coefficient (Equation 10 of Jiang & Oh (2018) . κ is the CR diffusivity. Source terms in the momentum and energy equations depend on this interaction coefficient and encapsulate how CRs exchange momentum and energy with the gas. In the gas thermal energy and CR energy equations, these source/sink terms account for collisionless energy transfer from the CRs to the thermal gas due to wave damping. We will sometimes refer to this collisionless energy transfer as "streaming energy loss." Time-dependent hydromagnetic wave energy is not explicitly tracked here, as streaming instability growth times are generally much shorter than other timescales of interest; waves are assumed to couple CRs to gas unless ∇P CR → 0 (see Thomas & Pfrommer 2019 for an implementation which tracks wave energy).

Generating Turbulence
To generate turbulence, we take advantage of the turbulent stirring module in Athena++, which uses an Ornstein-Uhlenbeck process (Uhlenbeck & Ornstein 1930) to smoothly generate a prescribed mixture of compressive and solenoidal velocity perturbationsv = f shearvshear + (1 − f shear )v compressive over a correlation time (similar methods are employed in e.g. Eswaran & Pope 1988;Federrath et al. 2008Federrath et al. , 2010Lynn et al. 2012). The turbulent reacceleration rate of CRs depends only on the compressive component of turbulence. For this paper, we will focus on purely compressive forcing (f shear = 0); increasingly solenoidal forcing leads to weaker compressions and rarefactions for a given Mach number, therefore decreasing the turbulent reacceleration rate further until it becomes zero with purely solenoidal perturbations. The advantage of purely compressive driving is that the turbulent dynamo (which depends on solenoidal driving) does not operate, so our simulations have roughly constant plasma β. Otherwise, since we need to drive the simulations for many eddy turnover times to see CR turbulent reacceleration, turbulent magnetic field amplification would obscure the plasma β dependence of CR reacceleration that we wish to study. Instead of enforcing a specific turbulent powerlaw over many scales, we use parabolic driving between modes 1 < k < 3, and the resulting turbulent cascade to higher wavenumbers is created organically. For driving, we set the autocorrelation timescale to be t corr = L/c s and drive fluctuations every t drive = 2 × 10 −3 (L/c s ). Our results are not sensitive to these assumptions.
To check simulated growth times vs analytic predictions, we use grids with 64 3 -256 3 cells in a square domain of length 2L, and we measure the outer-scale turbulent eddy (the one with the most power) to be of size ∼ 2L/3. This low resolution gives us only a short turbulence inertial range -dissipation sets in at a scale ≈ 30 cell widths (Federrath et al. 2010(Federrath et al. , 2021) -but in practice, we find that even our lowest resolution runs Table 1. Range of simulation parameters used in §3.3, 3.4. In §3.3, to check analytic growth predictions, we fix β = 2, the initial PCR/Pg = 100, and the grid size to 64 3 . We test convergence with respect to grid size and maximum speed of light, vm, in the Appendix. In §3.4, all parameters are the same, but we vary β between 2 and 200, and we run a set of higher resolution simulations with grid sizes of 128 3 -512 3 . with a 64 3 box give reasonably converged reacceleration rates matching analytic expectations 6 . This allows us to run a large parameter study to verify the scalings of Section 2.

Box Size
In §3.4, we find that we must be a bit more careful with our resolution and simulation box sizes. While reacceleration rates are again converged even at our lowest resolution of 2L/64, streaming energy loss rates in low-β plasmas are sensitive to the amount of magnetic field tangling captured in the simulation; to show this, we vary some of our simulations to have resolutions of 2L/128, 2L/256, and 2L/512. For each simulation, we choose a fiducial maximum speed of light v m = 50c s , and we show convergence with respect to this choice in the Appendix. Table 1 compiles our fiducial parameters. In Section 3.4, all parameters are the same, except we vary the initial β between 2 and 200.

Results: Pure Diffusion
We begin with simulations without streaming, i.e. solely anisotropic diffusion along the local magnetic field direction. Each simulation has β = 2 and low initial CR pressure (P g /P CR = 100) so the phase velocity is approximately the isothermal sound speed, c s . We choose three different forcing amplitudes to give turbu- lent Mach numbers of M ∼ 0.25, 0.35, and 0.5 (measured at late times when the RMS velocity has saturated), and we vary the diffusion coefficient over more than 3 orders of magnitude to test whether we recover the growth time prediction of Equation 6.
The results are plotted in Figure 1. The simulation growth times are calculated by summing the total CR energy in the box, E tot CR , and fitting a line to the log 10 (E tot CR ) curve in the exponential growth phase between 3 and 6 eddy turnover times. This time interval was chosen because it spans many time outputs and occurs after the kinetic energy has saturated but before the CR energy rises past equipartition with the thermal energy, at which point changes to the gas compressibility (changes to v ph ) decrease the growth rate and affect the normalization of κ relative to L 0 v ph . A rolling derivative confirms that these growth times are representative of the exponential growth phase.
Our simulations match analytic expectations (Equation 6) very well, not only reproducing the correct scaling with Mach number (t grow /τ eddy ∝ 1/M) but also the correct growth curve shape in the κ ≤ v ph L 0 regime if one sets L 0 /L 1 ≈ 20. The appropriate value of L 0 /L 1 is determined, in reality, by the characteristic width of a shock front, L 1 , relative to the outer driving scale, L 0 . In simulations, however, L 1 is limited by resolution; our choice of L 0 /L 1 ≈ 20 is motivated by L 1 spanning roughly one cell width in our 64 3 grid. In the Appendix, we show how other choices of L 0 /L 1 , which should correspond to other spatial resolutions, change the growth curves.
In the κ > v ph L 0 limit, we also recover the expected decrease in growth time due to anisotropic rather than isotropic diffusion (Chandran & Maron 2004): t grow → t grow /( κ/vL 0 ). The dashed line in Figure 1 shows t grow /τ eddy in the κ >> v ph L 0 regime when correcting for anisotropic diffusion, and our simulated growth times match this trend very well.

Results: With Streaming
We now include streaming transport and a variety of initial plasma β ∈ [2, 200] to test Equation 12 in the streaming-dominated regime, as well as the predicted collapse of t stream grow → t diff grow in the κ >> c s L 0 regime. We vary the diffusion coefficient between κ = 0.15L 0 c s (the maximal growth case without streaming) and κ = 15L 0 c s .
An important note is that we use the same forcing for each simulation, i.e. = ρv 3 /L 0 is held constant; therefore, for increasing plasma β, the magnetic field back-reacts on the flow less, leading to a slightly higher turbulent Mach number and average velocity divergence. This is a mild effect. Nonetheless, to make consistent comparisons, we run each simulation with and without streaming for each plasma β and focus our analysis on the ratio of growth rates. As in Section 3.3, the initial CR pressure is 1% of the thermal pressure, so CRs don't significantly affect the properties of turbulence.
Our results are shown in Figure 2. The left panel shows growth times for streaming and non-streaming simulations. These are normalized by the outer scale eddy turnover time, which assumes M ∼ 0.5 for each run, but decreased growth times for pure diffusion at higher β reflect the reduced MHD forces on the flow, leading to M > 0.5. The x-axis denotes the evolved plasma β of the simulation; because our forcing has no solenoidal component, magnetic field amplification is relatively inefficient, and plasma β decreases by a factor less than 2 during the time interval of our analysis.
The right panel shows the ratio of reacceleration times for anisotropically diffusing CRs, with and without additional streaming. At low β, typical of the ISM, the growth time is an order of magnitude longer than the growth time with pure diffusion. The discrepancy decreases at higher β but is still a factor of ∼ 2 even at β ∼ 10. The dashed blue line shows the expected ratio in the κ < v ph L 0 regime (Equation 12). Indeed, our results broadly fit with expectations. When streaming is important (κ = 0.15L 0 c s ), energy gains are largely reversible, and growth times are much longer at low β. At higher β, the discrepancy drops, following the predicted scaling of Equation 12 fairly well. For κ = 15L 0 c s , the effects of streaming are lower, following the expectation from Equation 13 for the κ >> v ph L 0 regime (the blue dot-dashed line). Note that growth times in this regime, even without additional streaming modifications, are already much longer than the minimum growth time when κ < L 0 c s . Overall, the main point is clear: CR streaming alters CR-turbulence interactions and significantly decreases reacceleration rates in low-β, ISM-like environments compared to the pure diffusion growth rates first derived in Ptuskin (1988) and frequently assumed in CR propagation models.

Saturation
Now that we've analyzed the linear regime when P CR << P g and compared to analytic growth rates, we move to the nonlinear regime when CRs become dynamically significant. In principle, CRs can back-react on the turbulent flow, changing its energy spectrum and cascade rates, much as for instance magnetic fields alter hydrodynamic turbulence, resulting in MHD turbulence. A fuller exploration of this interesting issue requires higher resolution simulations and is the subject of our follow-up paper (Bustard and Oh, in prep). Given our focus on energetics, however, we can at least study in this work how turbulent kinetic energy is dissipated. In hydrodynamics, the turbulent energy forcing rate˜ ≈ ρv 3 /l ≈const (independent of scale l) is equal to the gas heating rate. In MHD, some fraction of the kinetic energy is used to amplify magnetic fields via the turbulent dynamo. Similarly, in a two-fluid CR-gas system, some fraction of the kinetic energy does not dissipate into heat at small scales, but is instead diverted through the CR population.
SinceṖ CR ∝ P CR /t grow , where t grow is initially independent of P CR , we expect exponential growth in P CR (t), similar to the magnetic turbulent dynamo wherė P B ∝ ωP B (where ω is the fluid vorticity), and magnetic fields initially grow exponentially. For the turbulent dynamo, the fluid vorticity decreases due to magnetic tension from the growing magnetic field, and the dynamo transitions from exponential to linear growth, and finally saturation at roughly equipartition values.
The non-linear saturation of CR turbulent reacceleration is also interesting. Presumably, for fixed driving, as P tot increases, the fluid becomes less compressible due to an increase in v ph , which decreases M ph and increases t grow . While saturation in the turbulent dynamo is due to a decrease in ω = ∇ × v, it is related to a decrease in ∇·v for CR reacceleration. This holds true in our simulations. In Figure 3, we compare simulated growth times for diffusion-only simulations with κ = 0.15L 0 c s and varying P CR /P g > 1. We find, as expected given Equation 2, that due to the rise in v ph due to the increase in CR energy density, the growth time increases secularly: t grow ∝ v 2 ph ∝ P CR /P g . Thus,Ṗ CR ∝ P CR /t grow → constant, and growth transitions from exponential to linear, just as for the magnetic turbulent dynamo. Note that this saturation only occurs in the κ < v ph L regime and once P CR ≥ P g . Otherwise, there is little change in t grow .

SUPPRESSION OF STREAMING ENERGY LOSS
We additionally quantify how magnetic field strength affects the streaming energy loss rate. This can be thought of interchangeably as the gas heating rate due to CR streaming, if the gas equation of state is adiabatic, but here no heating occurs because the gas is isothermal.  (Equation 13). Compared to canonical rates assuming pure diffusion, reacceleration is less efficient in low-β environments when streaming is included. Figure 4 shows the collisionless energy loss time, calculated as t CR,loss = P CR /|v A · ∇P CR | for simulations with varying plasma β. It also shows the naive expectation that, if the CR scale height is approximately the outer eddy scale L 0 , the relative loss time should be L 0 /v A , so that loss times decrease monotonically for stronger B-fields. In fact, that is far from the case. At low β, the loss time is orders of magnitude longer, rising with increasing resolution at low β but seemingly converged with resolution for β > 10. The loss time reaches a minimum near β ∼ 10, which corresponds to an Alfvén Mach number M A = v/v A ∼ 1. Note also that streaming energy loss time is not inversely correlated with the reacceleration time, as it would be if streaming energy loss were the main factor stunting reacceleration. Instead, slow reacceleration in the low-β regime is also accompanied by slow energy loss, pointing to the modified phase-shifted CR force as the dominant correction to reacceleration rates (Equations 11 -13). Why does CR energy loss have this M A dependence? This largely arises due to misalignment between the magnetic field and CR pressure gradient. The streaming loss rate v A · ∇P CR is sensitive to the angle between the magnetic field B and CR pressure gradient ∇P CR . In M A < 1 turbulence, magnetic tension is strong, and field line tangling is suppressed. While compressions (and hence ∇P CR ) can occur in all directions, the mean magnetic field maintains its initial orientation, so that v A · ∇P CR is suppressed. This effect is apparent in Figure 5, which shows slice plots of three 256 3 simulations along the z = 0 axis after 5.8 eddy turnover times, when the turbulence is developed. Each column corresponds to different evolved plasma β, going from strong field (left) to weak field (right). The top row shows density with magnetic field lines overplotted; peaks and troughs are most apparent in the weak field case since magnetic tension is weakest, also allowing the field to tangle more easily. The middle row shows the collisionless energy loss rate, calculated as |v A · ∇P CR |, divided by |v A ∇P CR |. This effectively quantifies the misalignment between the magnetic field and CR pressure gradient. Clearly, the β ∼ 1 simulation shows large regions of misalignment, while the β ∼ 100 simulation has a more isotropic magnetic field and, hence, better alignment between B and ∇P CR .
Eventually, one would like a fitting formula for this loss rate for use in e.g. sub-grid models of CR transport for galaxy evolution simulations; however, that is premature at this stage. Clearly, the degree of misalignment depends on resolution at low β. For β = 2, our lowest resolution 64 3 simulation has energy loss suppressed by a factor of ∼ 50, while in our highest resolution 512 3 simulation, the loss is suppressed by a factor of ∼ 10. Higher resolution allows for a more accurate capture of field line tangling. Future higher resolution simulations are needed to probe convergence and facilitate the development of a sub-grid CR heating model, but we can clearly conclude that low-resolution galaxy evolution simulations, while possibly capturing reacceleration (which is dominated by the outer scale eddies) likely  Figure 3. Top: Two representative simulations starting from PCR/Pg = 1/100 and 1, run out to hundreds of eddy turnover times. Growth is exponential, with increasing growth times and a transition to quasi-linear growth once PCR > Pg. Bottom: Simulated growth time vs PCR/P g for simulations with constant κ < v ph L0 and varying PCR > Pg. We expect (and find) that tgrow ∝ v 2 ph ∝ PCR/Pg, consistent with Equation 2.
cannot capture true loss rates. However, at higher β, the suppression is weaker (approaching a few tenths) for all resolutions we tried, consistent with the field becoming more isotropic, and t CR,loss ∼ L 0 /v A is approximately realized.
An additional interesting effect is that streaming CRs create "bottlenecks" on a timescale of order the Alfvén crossing time. In the streaming dominated regime, P CR ∝ (v A + v) −γ CR along a flux tube. A minimum in (v A + v) (e.g., due to a density spike) creates a situation where CRs would have to stream up their pressure gradient, which is not possible. Instead, CRs form a flat pressure profile where no momentum or energy are transferred to the gas (Skilling 1971;Wiener et al. 2017Wiener et al. , 2019Bustard & Zweibel 2021); multiple regions can take the form of a 'staircase' structure (Tsung et al. 2022). In sub-Alfvénic turbulence, these bottlenecks, which form after an Alfvén crossing time, can form before the flow randomizes during an eddy turnover time, potentially creating ∇P CR ≈ 0 regions where streaming energy loss is suppressed. We see this in the bottom row of Figure  5, which shows the parallel (magnetic field-aligned) CR scale length L CR = P CR /∇ || P CR = P CR /(b · ∇P CR ) relative to the outer eddy scale. In large portions of the simulation box, L CR >> L 0 for β ∼ 1. At higher β, due to the combined effects of greater field line tangling and slower development of bottlenecks, L CR ∼ L 0 for much more of the volume. In 1D stratified media, bottlenecks have only weak effects on the net loss rate, though they greatly increase the spatial and temporal intermittency of energy loss (Tsung et al. 2022). In Bustard and Oh 2022, in prep, we probe CR energy loss / gas heating in greater detail and at higher spatial resolution, specifically how much gas heating occurs via grid-scale dissipation or via streaming energy loss.

Varying the Diffusion Coefficient
As a sanity check, we briefly explore sensitivities of the streaming energy loss time to varying κ, with each simulation run on a 64 3 grid and with initial plasma β varying from 2 to 200. As κ increases, both the magnitude of energy loss and the discrepancy between simulated loss and expectation should decrease. In the top panel of Figure 6, the solid lines show the expected loss times, either L 0 /v A for the streaming-dominated run or (κ/v)/v A for the diffusion dominated runs, the difference owing to the appropriate guess for the CR scale height. Indeed, we see that, as the scale height increases with increasing diffusivity, the loss rate goes down for all β. Energy loss suppression due to misalignment of B and ∇P CR also decreases, since ∇P CR is no longer set by compression but by the speed of field aligned diffusion.

Solenoidal vs Compressive Driving
In hydrodynamic 3D turbulence, the fraction of power in compressive modes is given by (Federrath et al. 2010) where f shear is a stirring parameter we can vary between 0 (purely compressive forcing) and 1 (purely solenoidal forcing). A natural mixture of f shear = 0.5 yields F long /F tot = 1/3. Since CR pressure gradients, which lead to collisionless energy loss, are developed by compressive fluctuations, we expect the loss rate to be a . Each simulation has Ms ≈ 0.5 and streaming-dominated transport, with κ ∼ 0.15L0cs, hence L0/vA = τstream < L 2 0 /κ = τ diff for β 100. At high-β, turbulence easily tangles field lines, and the energy loss time approaches L0/vA. At low-β, when there is a strong guide field, field line tangling is more difficult, and misalignment between the magnetic field and the CR pressure gradient suppresses energy loss by a factor ∼ 0.1. The amount of energy loss is resolution-dependent for low β. the same total driving rate split into different mixtures of compressive and solenoidal driving. Note that, because solenoidal motions more easily amplify magnetic fields than compressive motions, the saturated plasma β shown on the x-axis differs greatly from the initial plasma β, saturating, for example, near β ∼ 10 (M A ∼ 1) for initially super-Alfvenic turbulence. As expected, f shear = 0, corresponding to our purely compressive driving simulations shown in the top panel of Figure  4, gives the largest CR reacceleration rates (∼ 72t eddy for f shear = 0 vs ∼ 170t eddy for f shear = 0.5) and the largest energy loss rates. Going from purely compressive forcing (f shear = 0) to a natural mixture (f shear = 0.5) decreases the energy loss rate by a factor of 5 or more. That this decrease is greater than a factor of 3, assuming F long /F tot = 1/3, can be understood if δP CR /P CR and δv/v are not well-coupled. We speculate that solenoidal motions lead to longer CR scale lengths due to the aforementioned bottleneck effect for streaming CRs; this decreases δP CR /P CR and ∇P CR . This is worth studying in simulations where CRs are dynamically important, but for now, we just emphasize that solenoidal driving decreases both CR reacceleration rates and CR energy loss rates for a given driving rate .

Issues with Reacceleration in the ISM
In propagation models, reacceleration is typically included as a default, and it is an attractive alternative to otherwise "leaky box" models in that it provides a natural fit to low-energy (R ≈ 1 GV) boron-to-carbon data while maintaining the standard paradigm of diffusive propagation, i.e. a single power-law dependence of the diffusion coefficient D(E) or, inversely, the escape path length λ esc ∼ 1/D(E) (Heinbach & Simon 1995). At the same time, canonical reacceleration rates are energetically troubling, as surprisingly high fractions of total CR energization (up to ∼ 50%, comparable to the contribution from diffusive shock acceleration) have been attributed to turbulent reacceleration (Thornbury & Drury 2014;Drury & Strong 2017), and they are increasingly in tension with data. If reacceleration were the dominant acceleration mechanism in the ∼ 1 − 100 GeV range in our Galaxy, we would then see a progres- sive increase in secondary to primary ratio as a function of energy; instead, we see otherwise (Strong et al. 2007;Gabici et al. 2019). Additionally, at this low-energy end, a growing wealth of multiwavelength data, specifically synchrotron measurements, are best fit in models with little or no reacceleration (Trotta et al. 2011;Di Bernardo et al. 2013;Orlando & Strong 2013;Gabici et al. 2019), in tension with canonical reacceleration rates assuming purely diffusive CR transport. These low energy CRs, however, are precisely those which are self-confined, and where the impact of CR streaming must be considered. Hopkins et al. (2022) argues analytically that timescales for CR reacceleration, energy loss, and convection obey the ordering τ conv < τ loss < τ reacc , in which case reacceleration is negligible and convection in a large-scale wind presents a compelling alternative to explain the bump in B/C at low energies. Our simulations specifically demonstrate that non-resonant reacceleration from subsonic, compressive turbulence is ineffective in typical ISM environments with low β, although here we find the suppression is due to phase-shifted CR forces rather than CR streaming losses. Streaming stunts energy gain due to large-scale compressive motions, even when CRs diffuse slowly and would otherwise gain energy from turbulence at the maximal rate. Even if typi-cal ISM CR diffusivities lie in the κ >> v ph L 0 regime 7 , where streaming corrections are smaller, reacceleration is, in any case, less efficient than in the well-trapped (κ < v ph L 0 ) regime. That turbulence is not purely compressive further decreases the efficiency of non-resonant reacceleration since solenoidal motions don't energize CRs. All considered, the neglect of streaming in Galactic transport codes (see Hanasz et al. 2021 for a recent review) overestimates CR energization by second order Fermi acceleration by a large factor in β ∼ 1 plasmas like the warm ISM. What about other reacceleration mechanisms? Resonant second order Fermi acceleration, which relies on the presence of both forward and backward propagating waves, is not possible for self-confined CRs, for which the only waves that are excited are those that co-move with the CRs (Zweibel 2017). Reacceleration by transittime damping (TTD), essentially the magnetized version of Landau damping, is similarly inefficient at these energies: scattering rates from TTD are orders of magnitude lower than the rate of gyroresonant scattering by selfconfinement (e.g. Figure 2 of Yan & Lazarian 2004).
While we here consider only second order Fermi mechanisms, first order Fermi mechanisms, such as DSA or turbulent reconnection (Lazarian et al. 2020), may be efficient in regions of the ISM such as molecular clouds (Gaches et al. 2021) or superbubbles (Vieu et al. 2022) where turbulence is supersonic and therefore cascades into shocks. These first order processes imprint distinctly different spectral signatures than second order Fermi acceleration.

Reacceleration in the CGM and ICM
While streaming considerably stunts reacceleration in low-β galactic environments, non-resonant reacceleration in other environments is still quite plausible. Table  2 lists reasonable values for the phase velocity, outer driving scale, and plasma β in the warm ISM (WIM), CGM, and ICM. 7 Note, however, that this presumes a constant diffusion coefficient. The empirical diffusion coefficients which are used in Galactic propagation models more likely reflect conditions in the halo, where diffusing particles spend most of their time. In quasi-linear theory of self-confinement, where diffusion expresses transport relative to the Alfvén frame, the diffusion coefficient adapts to local conditions (e.g., see Wiener et al. 2017): where v D ∼ O(v A ) is the net drift velocity. This can in fact result in κ/csL 0 ∼ < 1, rather than κ/csL 0 1. Table 2. Typical values for the warm ISM (WIM), CGM, and ICM. Note that β appears to vary quite widely in simulations of the CGM, e.g. van de Voort et al. (2021), where β ∼ 0.01 in localized regions coincident with galactic outflows, but β ∼ 10 − 100 in quiescent regions. Recent observations of a fast radio burst passing through a foreground galaxy halo suggest β > 1 (Prochaska et al. 2019).
WIM CGM ICM v ph (cm/s) 10 6 10 7 10 7 L0 100 pc 1-10 kpc 100 kpc κcrit = v ph L0 (cm 2 /s) 3 × 10 26 3 × 10 28−29 3 × 10 30 β 1 1-10? 100 The CGM is a plausible candidate for efficient nonresonant reacceleration. Collisional loss rates are long in these diffuse galaxy halos, and for typically assumed CR diffusion coefficients of κ ∼ few ×10 28 cm 2 /s, κ < L 0 v ph and CRs may be well-trapped in turbulent eddies of the CGM, prolonging their residence time in the CGM and boosting their energy density. The effect of streaming is tied to the local plasma β and is therefore locationdependent. Observational constraints on plasma β in the CGM are sparse, but observations of a fast radio burst passing through a foreground galaxy halo suggest β > 1 (Prochaska et al. 2019). In simulations, CGM plasma β is often large (e.g., β ∼ 10 − 100 in FIRE simulations; Hopkins et al. 2020), though it can significantly fall in regions affected by galactic winds (van de Voort et al. 2021), with β ∼ 1 most favorable in largescale galactic winds. The efficacy of reacceleration could thus vary spatially, being most efficient in high-β regions, while being stunted in others. Overall, turbulent reacceleration could play a much more significant role in regulating the CR content of the CGM, compared to the ISM.
Reacceleration rates are largely unchanged by streaming losses in the high β environment typical of the ICM, where reacceleration is frequently invoked to explain radio emission in merging galaxy clusters. CR energy densities in the ICM are constrained to be quite low but are still consistent with models of (efficient) turbulent reacceleration (Brunetti & Lazarian 2007. While most theoretical analyses to-date of reacceleration in merging clusters have focused on reacceleration from TTD, we simply remind the reader that non-resonant reacceleration will, at some level, always be present and can have a competitive growth time compared to resonant reacceleration in high-β environments (e.g., see Fig. 6 of Brunetti & Lazarian 2007). This is useful to keep in mind due to the universality and simplicity of non-resonant reacceleration. The resonant case, on the other hand, relies on turbulence cascading down to gyroresonant scales and highly uncertain details of the spectra and damping scales in compressible MHD turbulence (Brunetti & Lazarian 2011;Miniati 2015;Pinzke et al. 2017). For instance, if compressible modes dissipate in weak shocks (Burgers turbulence), as is quite plausible, then particle acceleration rates are too low to explain observed giant radio halos. Even if compressible turbulence is Kraichnan, standard TTD on thermal particles gives problematic damping scales for turbulence, and a reduction of particle mean free path by plasma instabilities (mirror, firehose) might be necessary (Brunetti & Lazarian 2011). By constrast, non-resonant turbulent reacceleration is a comparatively robust, wellunderstand mechanism further validated by our numerical simulations.

Do Galaxy Evolution Simulations Accurately
Capture CR Energy Transfer?
In Appendix A, we briefly review the impact of spatial resolution. We find that turbulent reacceleration is adequately resolved as long as the outer scale of turbulence is sufficiently well resolved (by ∼ 20 cells). Most galaxy evolution simulations refine based on density, meaning they decrease spatial resolution going from the dense ISM to the diffuse halo. Fortunately, outer eddy scales similarly increase going from the ISM to the halo, meaning that simulations likely resolve at least the outer driving scale and, in fact, have a good chance of capturing accurate reacceleration rates. One caveat is turbulent reacceleration in the κ v ph L 0 regime. As the dynamic range of the simulation increases, the acceleration time decreases, due to acceleration by smaller eddies, particularly those of size l where κ ∼ v ph l.
Unfortunately, spatial resolution does affect the influence of CRs on the background gas. When small-scale field line tangling is not well-resolved, streaming energy loss is artificially reduced (Figure 4), especially in our low-β simulations. This has abundant implications for large-scale simulations. For instance, one effect of CRs that has garnered significant attention is their ability to drive multiphase galactic winds (Ipavich 1975;Breitschwerdt et al. 1991;Salem & Bryan 2014;Ruszkowski et al. 2017;Buck et al. 2020;Hopkins et al. 2020;Bustard et al. 2020;Quataert et al. 2022;Bustard & Zweibel 2021;. Indeed, turbulent reacceleration could affect CR wind driving -for instance, by reenergizing CRs in the halo despite strong energy losses in streaming scenarios arising from streaming down steep density gradients ). An unsolved issue with CR acceleration of multi-phase gas is how CRs accelerate cold (T ∼ 10 4 K) dense gas clouds. Do they provide direct acceleration by exerting CR forces on the cold gas, due to steep CR gradients which develop at the cold-hot gas interface (Wiener et al. 2017(Wiener et al. , 2019Brüggen & Scannapieco 2020;Bustard & Zweibel 2021)? Or is direct acceleration inefficient (as might be expected if field lines wrap around the cloud), and CRs first accelerate the background hot gas, which then accelerates the cold gas via mixing-induced momentum transfer (Gronke & Oh 2018? There are hints of the latter process in Bustard & Zweibel (2021); . The relative importance of indirect vs direct acceleration depends on the relative efficacy of CR thermal vs momentum driving, particularly in the background hot medium. Figure 4 suggests that, for M A < 1, the CR energy loss rate is a relatively flat or even rising function of β; a tangled field overcompensates for a lower magnetic field strength, at least until the field is isotropic (M A ∼ 1), at which point energy loss scales inversely with β. Hints of this behavior were recently seen in , where CR heating and, in turn, indirect cloud acceleration were actually more efficient in runs with decreased field strength. A goal of future work, currently premature given the incomplete parameter space we've so far explored, should be to develop a sub-grid model for energy loss / gas heating as a function of resolution, M A , and different turbulent driving modes. If convergence (not clearly seen in Figure 4) can be achieved with higher resolution runs, this may not be far off.

CONCLUSIONS
In this paper, we reviewed the non-resonant reacceleration of CRs in subsonic, compressive turbulence and derived heuristic modifications to reacceleration rates when CRs are self-confined or "streaming." As CRs are believed to be self-confined up to E ∼ < 300GeV CR reacceleration rates are modified throughout this entire energy range. After describing our analytic expectations, we ran a suite of MHD simulations to verify previously derived reacceleration rates for purely diffusive CRs (Ptuskin 1988;Chandran & Maron 2004), test the effects of streaming, and probe the nonlinear regime. Our main findings are as follows: • Our simulations, which show a Burger's-like power spectrum, verify the analytic reacceleration rates derived for a k −2 spectrum (Ptuskin 1988; Equations 1,2), including the expected modifications due to anisotropic field-aligned transport at large κ (Chandran & Maron 2004; Equation 8). In particular, reacceleration rates peak for diffusion coefficients κ ∼ 0.1L 0 v ph , and have scalings consistent with analytic expectations. To our knowledge, this is the first time CR turbulent reacceleration has been shown in MHD simulations with a fluid CR treatment and is an encouraging test of the Athena++ CR module implemented in Jiang & Oh (2018).
• When CR streaming is introduced, the rate of net energy gain can be substantially suppressed compared to the diffusion only case, due to modified phase shifts between CR and gas variables. For κ ∼ < L 0 v ph , the regime where reacceleration is canonically most efficient, growth times in low-β plasmas typical of the ISM are significantly longer than canonical expectations (Equation 12; right panel of Fig 2). By contrast, in high β environments like the CGM and ICM, streaming does not have a significant impact on turbulent reacceleration. In diffusion-dominated regimes (κ > L 0 v ph ), the impact of streaming is milder, but reacceleration is in any case already inefficient. New reacceleration times t grow ∼ p 2 Dpp in κ < v ph L 0 and κ > v ph L 0 regimes, respectively, are given below, assuming a k −2 kinetic energy spectrum (see Equation 4 and the ensuing discussion for how to calculate reacceleration rates with alternative power spectra): L 0 is the outer eddy scale, L 1 is the smaller characteristic scale of shocks in the medium, and f corr encodes corrections due to anisotropic diffusion (Chandran & Maron 2004) and streaming transport assuming an isothermal gas (this work): The β-dependent terms are relevant for self-confined CRs with energy E 300 GeV; at higher energies, CRs are no longer self-confined and these can be dropped.
To diagnose the limitations of lower-resolution galaxy evolution simulations and as a step towards sub-grid modeling of CR energetics and influence, we also determine some sensitivities to resolution.
• Reacceleration rates with pure diffusion are largely insensitive to resolution (the minimum growth time is well-captured even when the outer eddy scale is only resolved by 20 cells), but higher resolution more accurately captures power at small scales, boosting simulated reacceleration rates. This is important in the κ v ph L 0 regime.
• Streaming energy loss v A · ∇P CR is more strongly dependent on Alfvén Mach number M A and simulation resolution. Because streaming energy loss only occurs when the CR pressure gradient is aligned with the magnetic field, not every compression induces streaming energy loss. This misalignment between the magnetic field and ∇P CR is a function of plasma β since it is more difficult for turbulence to tangle magnetic field lines in highly magnetized plasmas. The net result is that relative CR loss rates are much lower than v A /L 0 in low-β plasmas (Figures 4 and 5) and are clearly sensitive to resolution. When field line tangling is resolved, the average misalignment between ∇P CR and the magnetic field decreases. Counterintuitively, due to the countervailing effects of increased CR streaming speeds and decreased field alignment at lower β, our highest resolution results over β ∼ 1 − 100 show that CR energy loss is insensitive to plasma β. This also shows that the reduced acceleration rates at low β is not due to increased CR losses.
An important issue we have not considered in this paper is the effect of density stratification, which results in a background CR gradient. This provides constant CR coupling, pressure support and heating, and if sufficiently strong, can drive a wind. How does CR reacceleration proceed in such a background? One issue is that the density fluctuations created by turbulence can create CR 'bottlenecks' (Skilling 1971;Wiener et al. 2017): small decreases in the Alfvén speed along a magnetic flux tube cause CRs to pile up or "bottleneck." For purely streaming CRs, they readjust to these conditions and create a flat CR pressure upstream of the dip in Alfvén speed. As ∇P CR → 0, CRs no longer excite confin-ing Alfvén waves and instead free-stream at close to the speed of light, no longer transferring energy or momentum to the gas. It is as yet unclear how this stochastic coupling affects CR energization and escape rates from stratified, turbulent media. We have also not considered the case where CR phase shifts and/or heating conspire to reverse the sign of energy transfer, such that CRs give energy to gas motions, rather than vice-versa (Begelman & Zweibel 1994;Tsung et al. 2022). This would occur when β ≤ 0.25, a regime we plan to study in future work.
Finally, although it is clear that CRs can damp compressive turbulence, their back-reaction on the turbulent cascade is relatively unexplored but possibly significant, for instance, in CR-dominated galaxy halos (Ji et al. 2020). We explore this in a series of higher-resolution follow-up simulations (Bustard and Oh 2022, in prep), where we find that CRs, even with very low reacceleration rates, can absorb a significant fraction of large-scale turbulent energy, subsequently modifying the compressive cascade. some care must be taken with the flux term, specifically the value of v m , which must be greater than all other speeds in the system (including the CR propagation speed). We fiducially set v m = 50c s , which gives converged reacceleration rates at not only our fiducial resolution of 2L/∆x = 64 but also at higher resolutions 2L/∆x = 128, 256. While not shown here, we also find that v m = 50c s is necessary to get converged CR loss rates, with lower v m artificially boosting loss rates. Figure 7 shows simulation data with maximum speed of light in the range v m /c s = 10 to 400, each for purely compressive turbulence with turbulent Mach number M ∼ 0.5. Simulated growth times are well converged with respect to v m and show an excellent match to the analytic derivation from Ptuskin (1988), assuming L 0 /L 1 = 20.
This convergence may seem a bit surprising at first since the time-dependent flux term 1/v 2 m ∂F CR /∂t is only small when v A /v m < ∆x/L; however, this is only the criterion for the two-moment equations to effectively collapse to the one-moment equations. Many of the test problems in e.g. Jiang & Oh (2018); Tsung et al. (2021) violate this criteria but show convergence to analytic results. Indeed, one of the most powerful components of the two-moment method is the ability to unlock stability and convergence from the quadratic time-stepping requirement needed for one-moment implementations.
Convergence with respect to spatial resolution proves a bit trickier in our simulations. For instance, in Figure 2, one can see that the growth times for pure diffusion continue to decrease with increasing resolution. The value of L 1 , which we associate with the width of a shock front in the medium, should be set by the spatial resolution. For our 64 3 simulations, the outer scale is resolved by ∼ 20 cells in each direction, making L 0 /L 1 ∼ 20. In real astrophysical plasmas, though, the scale separation is much larger, and the velocity divergence, which ultimately energizes CRs, gets additional contributions from smaller scales. Lines in Figure 7 denote different L 0 /L 1 . Notably, the growth times in the κ < v ph L 0 regime are much shorter; presumably, simulations of higher resolution would match these curves. To test this, we ran a few 128 3 and 256 3 simulations, again using v m /c s = 50. Growth times in the κ < v ph L 0 regime do decrease: The 128 3 runs match the analytic curve with L 0 /L 1 = 40 quite well; the 256 3 runs, for which we have a smaller set of data points, show somewhat shorter growth times than the minimum of the L 0 /L 1 = 80 curve, but the growth time near κ/v ph L 0 ∼ 0.01 returns to the L 0 /L 1 = 80 curve. Overall, we recover the expected trend that higher resolution in the κ < v ph L 0 regime should lead to shorter growth times. Since CRs with κ > v ph L 0 diffuse quickly over small-scale structures, higher resolution makes very little difference in this regime. Analytic, L 0 /L 1 = 20 Analytic, L 0 /L 1 = 40 Analytic, L 0 /L 1 = 80 Analytic, L 0 /L 1 = 1000 Analytic (Anisotropic) Figure 7. Testing growth time convergence with respect to the reduced speed of light vm (relative to the fastest propagation speed, cs), and numerical resolution. Simulations with vm/cs = 50 appear well-converged. Magenta and green points denote simulations with 128 3 and 256 3 resolution, respectively, for which the velocity divergence and hence compressional heating is slightly larger and decreases the reacceleration time in the well-trapped κ < v ph L0 regime. For large κ, CRs diffuse quickly over small-scale structures, which then do not need to be well-resolved to get converged reacceleration rates.