Galactic Cosmic-ray Scattering due to Intermittent Structures

Cosmic rays (CRs) with energies $\ll$ TeV comprise a significant component of the interstellar medium (ISM). Major uncertainties in CR behavior on observable scales (much larger than CR gyroradii) stem from how magnetic fluctuations scatter CRs in pitch angle. Traditional first-principles models, which assume these magnetic fluctuations are weak and uniformly scatter CRs in a homogeneous ISM, struggle to reproduce basic observables such as the dependence of CR residence times and scattering rates on rigidity. We therefore explore a new category of"patchy"CR scattering models, wherein CRs are predominantly scattered by intermittent strong scattering structures with small volume-filling factors. These models produce the observed rigidity dependence with a simple size distribution constraint, such that larger scattering structures are rarer but can scatter a wider range of CR energies. To reproduce the empirically-inferred CR scattering rates, the mean free path between scattering structures must be $\ell_{\rm mfp} \sim 10$ pc at GeV energies. We derive constraints on the sizes, internal properties, mass/volume-filling factors, and the number density any such structures would need to be both physically and observationally consistent. We consider a range of candidate structures, both large-scale (e.g. H II regions) and small-scale (e.g. intermittent turbulent structures, perhaps even associated with radio plasma scattering) and show that while many macroscopic candidates can be immediately ruled out as the primary CR scattering sites, many smaller structures remain viable and merit further theoretical study. We discuss future observational constraints that could test these models.

The vast majority of the literature studying low-energy CR propagation and dynamics has focused on simple, phenomenological prescriptions for the effective CR transport rates within the ISM, typically parameterized with an effective diffusion coefficient κeff, or streaming speed vst,eff.Using a variety of methods, classic studies constrained these coefficients by comparing detailed models of CR propagation in a Galactic background to observed CR spectra in the Solar system (including many CR species and a range of CR energies, ratios of primary-to-secondary CRs, radioactive and isotopic abundances, as well as the CR anisotropy on the sky), and/or to diffuse γ-ray observations (Strong & Moskalenko 2001;Vladimirov et al. 2012;Blasi & Amato 2012;Gaggero et al. 2015;Guo et al. 2016;Jóhannesson et al. 2016;Cummings et al. 2016a;Korsmeier & Cuoco 2016;Evoli et al. 2017).The "effective" coefficients inferred by such studies represent, by definition, some weighted average in the ISM between CR sources (e.g.supernova remnants in ⋆ ibutsky@stanford.edu† NASA Hubble Fellow the Milky Way) and the Solar system.Together, the existing observations still only constrain the CR scattering physics in ISM conditions broadly similar to the Solar neighborhood.Since phenomenological models do not explain how such scattering rates actually arise or break many of the degeneracies between CR propagation models, it is by no means clear how to extrapolate their findings to distinct environments (e.g. the Galactic center, CGM, or IGM), galaxy types (e.g.dwarf, starburst, or high-redshift), or CR acceleration cites (e.g.AGN and quasars).Moreover, different choices for such extrapolation can lead to orders-of-magnitude different results in predicted galaxy properties (see references above and Hopkins et al. 2021b;Butsky et al. 2020Butsky et al. , 2023)).
In order to predict how CRs propagate on macro scales in different galactic environments, we first need to understand CR physics on micro scales.The gyroradii, rgyro, of low-energy CRs are extremely small compared to the "macroscopic" scales of the ISM/CGM (e.g.rgyro ∼ 0.1 au, for CRs at the ∼ GeV peak of the spectrum in the diffuse ISM).This means that CRs cannot simply escape their acceleration sites around supernovae at bulk speeds vcr = βcrc, but instead travel along magnetic field lines on gyro orbits with a characteristic gyro frequency Ω ∼ vcr/rgyro, and pitch angle µ ≡ vcr • b relative to the magnetic field direction b ≡ B/|B|.When CRs encounter magnetic field fluctuations, δB, they are scattered in pitch angle, producing some effective pitch-angle scattering rate, νeff.When averaged over large spatial and temporal scales, νeff leads to bulk CR transport that can be parameterized by some effective diffusion coefficient κeff ∼ v 2 cr /νeff and/or streaming speed vst,eff ∼ κeff|∇Pcr|/Pcr.As we summarize in Section 2, traditional scattering models typically assume that the magnetic fluctuations that scatter CRs are weak (|δB|/|B| ≪ 1) and uniformly distributed throughout the volume-filling ISM, but differ from each other in their proposed origin of the magnetic fluctuations.
Constraining CR scattering theories is extremely difficult since it is simply not possible to directly resolve the relevant gyroresonant scales in either ISM observations or in numerical simulations that also include macroscale ISM processes.Even idealized particle-in- CR Source (e.g.SNe)

Observer
Figure 1.Cartoon illustrating the difference between traditional models of sub-TeV CR pitch-angle scattering in the ISM (Section 2; top) and the "patchy" models proposed here (Section 3; bottom).Low-energy CRs are accelerated at some sites (e.g.supernovae) and confined to gyro orbits with small radii, rgyro, and some pitch angle, µ, along magnetic field lines, B, through the ISM, until they are observed."Traditional" models (Section 2) assume CR pitch angles are scattered constantly by ubiquitous (volume-filling factor ∼ 1) small-amplitude (|δB| ≪ |B|) magnetic fluctuations (with wavelength λ here) uniformly distributed though a medium that is homogeneous up to scales of the CR deflection length ℓ mfp ∼ 10 pc R δs GV , giving rise to some effective mean scattering rate ⟨ν eff ⟩ shown.Reproducing the observed dependence of CR scattering rates and residence times on rigidity therefore requires a specific ratio of median/typical volume-filling magnetic fluctuations of different scales λ.The "patchy" model (Section 3; bottom) assumes that CR scattering is dominated by strong scattering in intermittent "patches" or structures (gray ovals), which have a small volume-filling factor (with size ℓ S , and number density ns).Larger patches, which have an effectively larger "CR scattering optical depth" and scatter both low-and high-energy CRs, are rarer than smaller patches, which only scatter low-energy CRs.An appropriate distribution of smaller and larger patches can therefore produce the observed dependence of CR scattering on rigidity.cell type simulations of CR scattering in an ISM "patch" that is only modestly larger than the CR gyroradii (so obviously unable to span the full dynamic range of conditions) have only just become possible in the last few years (Bai et al. 2019;Bai 2022;Holcomb & Spitkovsky 2019;Plotnikov et al. 2021;Bambic et al. 2021;Ji et al. 2022;Ji & Hopkins 2022).Because of these challenges, existing CR scattering models remain wildly uncertain, even in the "typical" ISM.For example, when applied in galaxy simulations, state-ofthe-art scattering models that are calibrated to reproduce existing observations predict CR scattering rates that differ by as much as ten orders of magnitude (at ∼ GeV energies) in the ISM and predict qualitatively different scalings with properties like magnetic field strength and turbulence (Hopkins et al. 2021c).Additionally, multiple recent studies have pointed out that existing scattering models struggle to even qualitatively capture the correct dependence of CR scattering rate on rigidity at sub-TeV energies (Hopkins et al. 2022b;Kempski & Quataert 2022).
In this paper, we are therefore motivated to propose a novel picture for CR scattering, which may resolve some of these challenges.In Section 2, we summarize some of the central challenges facing "conventional" models in the recent literature.In Section 3, we present a new theoretical framework for "patchy" CR scattering, which is qualitatively distinct from traditional, continuous scattering models and derive a number of criteria any such model must obey to reproduce observations.We discuss a variety of candidate scattering structures, both macroscopic and small-scale in Section 4. We show which candidate scattering structures can be immediately ruled out and discuss the potential connections to intermittent structures in turbulence.In Section 5 we discuss the observational implications for these model categories.We summarize and conclude our results in Section 6.

THE PROBLEM WITH SIMPLE, HOMOGENEOUS THEORIES OF COSMIC-RAY SCATTERING
Empirical constraints on CR propagation in the Galaxy typically infer an effective diffusivity of the form Dxx ∼ κeff = D0 βcr (Rcr/Rcr,0) δs , equivalent to an angle-averaged CR pitch-angle scattering rate, where βcr = vcr/c is the CR velocity, RGV is the CR rigidity in GV, and typical values of the fit parameters correspond to ν0 ∼ 10 −9 s −1 , and 0.5 ≲ δs ≲ 0.7 (e.g.De La Torre Luque et al. 2021, and references in Section 1).This is roughly equivalent to an isotropic diffusion coefficient κeff ∼ 10 29 cm2 s −1 , or effective streaming velocity vst,eff ∼ 300 km s −1 , assuming a gradient scale length of 1 kpc.The effective CR scattering rate translates into a characteristic mean free (or deflection) time, tmfp ∼ 1/⟨ν⟩eff, or mean free (or deflection) path, ℓmfp ∼ vcr/⟨ν⟩eff ∼ (c /ν0) R δs GV , between O(1) deflections in µ.Regardless of the details, an extremely robust observational result is δs > 0: the CR residence times and scattering rates (which scale ∝ ⟨ν⟩eff ∝ R −δs GV ) must decrease with increasing CR rigidity at energies ∼ GeV-TeV, in order to reproduce any of the observed trends, e.g.secondary-to-primary or radioactive isotopic species ratios.
Historical theories of CR scattering (heuristically illustrated in Fig. 1) generally assume that the observed CR pitch-angle scattering is the cumulative result of a large number of uncorrelated, small perturbations to µ generated by encounters with a very large number of independent, small magnetic field fluctuations, δB.These theories also assume that the magnetic field fluctuations occur throughout an ISM that is statistically homogeneous on spatial scales from the gyro scale up to, or larger than, the CR deflection length.This is a questionable starting assumption considering the CR deflection length ℓmfp ∼ 10 pc R δs GV (as large as hundreds of pc for ∼ TeV CRs) is larger than the size scales of much ISM structure.Nonetheless, integrating over an ensemble of perturbations, this leads to the classic predicted scattering rate ⟨ν⟩eff ∼ (vcr/λ) |δB(λ)| 2 /|B| 2 (e.g.Völk 1973), where |δB(λ)| ≪ |B| represents some typical (e.g.root-mean-square) fluctuation amplitude with wavelength λ. 1 This applies equally well to both "self-confinement" theories (in which δB is sourced at λ ∼ rgyro by CR streaming instabilities; Kulsrud & Pearce 1969;Wentzel 1969;Skilling 1975) and "extrinsic turbulence" theories (in which δB is sourced by a turbulent cascade from vastly-larger ISM driving scales ;Jokipii 1966;Völk 1973).
However, as discussed extensively in Hopkins et al. (2022b) and Kempski & Quataert (2022), the most commonly-invoked selfconfinement and extrinsic turbulence theories based on the above assumptions do not reproduce the locally-observed CR spectra at sub-TeV energies.For example, putting in typical ISM values for the relevant parameters gives orders-of-magnitude different normalization (ν0) from that observed (a point already made in Chandran 2000a; Yan & Lazarian 2002;Hopkins et al. 2021c;Chan et al. 2019;Fornieri et al. 2021).However, even allowing for arbitrary renormalization, the CR spectra predicted by traditional CR scattering models will not have the correct shape if one assumes typical scaling parameters (also noted in Yan & Lazarian 2004;Fornieri et al. 2021).Most importantly, these models predict δs ≤ 0, i.e. longer CR residence times for higher-energy CRs (opposite the observed behavior).
The problem with assuming CR scattering is "continuous", is that the only way to reproduce the observed dependence ⟨ν⟩eff ∝ R −δs GV is to invoke some connection between the dominant wavelength of scattering modes and RGV (for example, assuming gyroresonant scattering λ ∼ rgyro ∝ RGV), and a specific power-law spectrum of fluctuations δB(λ). 2 But in self-confinement theory, the only stable steady-state solution is one where all CRs either stream at approximately the MHD Alfvén speed, or free-stream (unconfined) at c, in either case, clearly independent of rigidity (δs = 0).Any solutions out of equilibrium "collapse" to these states on a very rapid timescale (≲ Myr) -an issue that has been noted going back at least to Skilling (1971).In extrinsic turbulence models, it is common to make the phenomenological comparison to Kolmogorov or "Kraichnan"-type spectra δB(λ) ∝ λ 1/4−1/3 , which appear at first to give reasonable estimates for δs.The problem is that this assumes the turbulence is both undamped and isotropic down to scales at least as small as the gyroresonant wavelength, which cannot be true for sub-TeV CRs, where the gyro scales are much smaller than the Alfvén and Kolmogorov/damping/dissipation scales of the turbulence.In that regime, turbulence is highly anisotropic, and the parallel structure necessary for efficient CR scattering are suppressed (Goldreich & Sridhar 1995). 3The prediction is then that scattering is necessarily dominated by larger-scale modes, λ ≫ rgyro, which are independent of rgyro, so δs ≤ 0, a point also noted in Chandran (2000a); Fornieri et al. (2021).
There might be solutions to these problems involving, for example, alternative, volume-filling sources of the modes δB that differ from standard self-confinement or extrinsic turbulence theories (see Hopkins et al. 2022b for discussion).But thus far, there does not appear to be an example of such a model that has actually been shown to reproduce the observed CR spectra.

AN ALTERNATIVE: "PATCHY" SCATTERING BY MESO-SCALE STRUCTURES
Here, we consider a more radical alternative: patchy scattering, in which CR scattering rates are high in discrete "scattering structures" and low in-between such structures (illustrated in Fig. 1).
More formally, we drop the assumption of the "classic" models in Section 2 that the fluctuations, δB, from which CRs scatter are homogeneous and uniformly volume filling.In this case, the effective scattering rate, ⟨ν⟩eff, inferred from various observational constraints should not be thought of as a uniform, volume-filling rate, but as an average rate of encountering scattering structures.Equivalently, the observationally-inferred effective deflection length, ℓmfp, would no longer represent the length over which a sufficient number of small deflections are continuously accrued to change µ by O(1), but rather would represent the mean free path between scattering "patches," Importantly, we will show that this reinterpretation frees us from being forced to assume there is a unique one-to-one correspondence between the measured rigidity dependence, δs, and the shape of the power spectrum of magnetic fluctuations on gyro scales, as in Section 2. In the following sections, we discuss the basic constraints that a plausible scattering structure would need to meet.

Geometry and General Constraints
First, consider the basic geometric properties of candidate scattering structures.It is helpful to think of structures in terms of their effective dimensionality, D, equal to the number of dimensions along which the system has a highly elongated axis ratio.Using this definition, D = 2 describes sheet or pancake structures with two long axes (ℓL) and one short axis (ℓS), D = 1 describes filamentary (or tube-like) structures with one long axis and two short axes, and D = 0 describes spherical (or "point-like") structures with all axis lengths of order ℓS.The cross-sectional area of the structure, across a random set of viewing angles, is dominated by the two larger dimensions, , while the relative depth of the structure (as seen by e.g.CRs traversing it) is dominated by the short-axis distance, ℓS.From these definitions, the volume of a scattering structure scales as Vs ∼ AsℓS.

Size and Internal Scattering Constraints
The first, most basic, size constraint is that in order to scatter a CR with some rigidity RGV, the structure must be larger than that CR's gyroradius, where BµG is the magnetic field in microGauss.
Additionally, we assume that there is a local CR pitch-angle scattering rate, νs, inside the scattering structures.It is important to distinguish the CR scattering rate within the structures from the ISM-averaged effective scattering rate, ⟨ν⟩eff, as by definition νs ≫ ⟨ν⟩eff.Therefore, the depth of a scattering structure must be large enough so that the scattering time within it (∼ 1/νs) is shorter than the unscattered CR crossing time of that structure (∼ ℓS/vcr).In other words, the structure must be large enough so that CRs are reliably scattered as they traverse it, ℓS ≳ vcr/νs. (4) Another way of saying this is that the scattering "optical depth" to CRs of some rigidity RGV, τs ∼ νs ℓS/vcr ≳ 1, must exceed unity in order for CRs of that RGV to be strongly scattered.
On the other hand, the structure cannot be so large that CRs are effectively trapped inside of it for much longer than their inferred total residence times in the ISM.Assuming there is a large scattering rate inside the structure, the time CRs spend inside of the scattering structure, ts, is set by the effective diffusivity, κs ≈ v 2 cr /νs, giving ∆ts ∼ ℓ 2 S /κs ∼ (ℓS/vcr) 2 νs. (5) By definition, if this were longer than the ISM-averaged deflection time ∼ ℓmfp/vcr, then this would dominate the total residence time and exceed the limits above, violating the basic assumptions of our framework.Therefore, we require ∆ts < ℓmfp/vcr, or It is only possible to satisfy both equation ( 4) and equation ( 6) if the depth of the scattering structure is significantly smaller than the mean free path between scattering structures ℓS ≪ ℓmfp. ( But this is of course implicit in a "patchy" scenario.

Number Densities, Surface Densities, and Mass or Volume-Filling Factors
The mean free path between patches that can scatter CRs of a given rigidity RGV is set by the cross-sectional area of the structures as well as their relative abundance or "number density" ns, as ℓmfp ∼ 1/(nsAs).Given the observationally-constrained mean free path for CRs of that rigidity, the scattering structures therefore must have a volume-averaged ISM number density of Next, consider the typical surface density (or column density; Σs) of a structure, viewed from a random angle, Σs ∼ Ms/As, in terms of the mass of the structure (Ms) and its cross-sectional area (As).Using the relations above, Σs ∼ Ms/As ∼ ρsℓmfp, where ρs ∼ Msns is the volume-averaged mass density of the scattering structures within the ISM as a whole.Assuming that the total mass contained in the scattering structures is some fraction, fM, of the ISM (with some volume-averaged ISM density ρISM ≈ mHnISM), we can rewrite the expression for the surface density as where ρISM = mH nISM is the volume-averaged ISM density, mH is the mass of a hydrogen atom, and nISM ≈ 1 cm −3 .Assuming that fM < 1, we can place an upper limit on the surface density of scattering structures, where ρs = Ms/Vs is the internal density of a scattering structure, i.e.Σs < 5 × 10 −5 g cm −2 R δs GV or column density < 3 × 10 19 cm −2 R δs GV .Similarly, the volume-filling fraction of the scattering structures, fV , would be These above quantities (ns, Σs, fM, fV ) depend on CR rigidity either explicitly, or implicitly, through ℓmfp ∝ R δs GV .This tells us the bulk properties of the scattering patches that scatter lower-energy CRs (∼MeV; e.g.their number densities, column densities, mass and volume-filling factors) will differ from the bulk properties of the patches that scatter higher-energy CRs (∼TeV).Of course, the sizes of the structures (ℓS, ℓL) that scatter CRs of different RGV may also be distinct, as we might naturally expect from the constraints in Section 3.1.

Sufficiently Weak Scattering Between Structures
A key requirement of the patchy scattering model is that CR scattering not be dominated by scattering in the medium between patches.This means that the diffuse/volume-filling ISM cannot have substantial CR scattering from either extrinsic turbulence or from the CR streaming instability (SI).For the extrinsic turbulence case, this is easy to satisfy, as the more detailed calculations in Section 2 argue that the theoretically-favored extrinsic scattering rates in the warm ISM for sub-TeV CRs are orders of magnitude smaller than the mean observed ⟨ν⟩eff.However, one must still avoid runaway growth of the CR SI between patches, which would self-confine CRs to move no faster than the Alfvén speed, over-confining them especially at higher energies (Hopkins et al. 2022c).
There are a few ways in which the overconfinement problem could be avoided.For example, at energies ≫ GeV, the SI growth rate is proportional to the number density of CRs and so drops rapidly.Thus, for higher energy CRs, the constraint is not so severe, and it may be possible that CRs just around ∼ 1 GeV are self-confined while other physics take over between 1 − 1000 GeV (though this may require some fine-tuning; see Kempski & Quataert 2022).Additionally, some recent MHD-PIC simulations have argued that SI growth rates may be slower than expected from simple quasi-linear expressions (Bai et al. 2019;Holcomb & Spitkovsky 2019).
Alternatively, we can use the weak scattering requirement to place some constraints on the properties of the volume-filling ISM.In steady state, with some large-scale background CR gradient, there would be some net flux of CRs moving away from the galactic center, leading to the growth of SI on some timescales.The patchy scattering model can still hold, so long as the scattering rate due to the saturated SI in the medium between scattering patches is less than the observationally-inferred effective CR scattering rate, νSI < νeff.This is equivalent to requiring that the effective diffusivity due to the SI be larger than the empirically-constrained average diffusion coefficient, where

Key Requirements to Reproduce Observational Scalings and Differences from the Homogeneous Models
Provided a potential scattering structure meets all the above criteria, the key point is that by having discrete structures that are not volume filling, we no longer need to impose a specific scattering rate νs, or a specific distribution of magnetic field fluctuations δB(λ) inside the structure in order to produce the observed dependence of λmfp or ⟨ν⟩eff on CR rigidity (Section 2).Essentially, we have replaced the requirement that the observed scaling ⟨ν⟩eff ∝ R −δ GV reflects a specific power-law scaling of δB(λ) ∝ λ α B , with a different requirement: that the number of scattering structures that can scatter CRs of a given rigidity depends in a specific power-law fashion on that rigidity.
How plausible is this?First, we consider a simplified scenario in which scattering structures all have some internal scattering rate νs, and vary from each other only in size.In this case, the scattering patches would have a wide range of sizes ℓS, spanning a range from the smallest to largest gyroradii of the GeV-to-TeV CRs of interest (∼ 0.1 − 100 au for microGauss fields).Per Section 3.1, the patches that are able to scatter CRs of a given rigidity will have ℓS > ℓS,min ∼ rgyro ∝ RGV.Our number density constraint then becomes: ns(> ℓS) ∝ ℓ D−2−δs S ℓ −D L .Realistically, the scattering rate within the proposed structures will likely have some dependence on CR rigidity in addition to the size of the structure νs ∝ ℓ −α S S R −α R GV .Combined with the requirement that structures can scatter CRs (ℓS > vcr/νs) and the other scalings from Section 3.1-3.1.2,this gives a more generalized constraint on the number density, ns(> ℓS) ∝ ℓ L .In either case, for quasi-3D structures (D = 0) or ℓL ∝ ℓS, reasonable values for αS ≈ 1, αR ≈ 0.5, and 0.3 ≲ δs ≲ 0.7, we obtain ns(> ℓS) ∝ ℓ −αn S with 2 ≲ αn ≲ 3, which is broadly similar to the distribution of sizes of many classes of objects in the ISM including molecular clouds and H II regions, H I filaments, star clusters, stellar wind termination shocks/bubbles/magnetospheres (Guszejnov et al. 2018).For sheet-like structures (D = 2), the scaling is quite similar to the distribution of shock widths seen in supersonic, isothermal turbulence (Squire & Hopkins 2017;Mocz & Burkhart 2019).As advertised, the scaling does not depend sensitively or in the same manner as we discussed in Section 2 on how νs (or implicitly δB) depends on wavelength λ.
Thus while it has proven (surprisingly!) challenging to theoretically construct a power spectrum of magnetic fluctuations that satisfies the observational requirements from Section 2, it appears, at least in principle, straightforward to conceive of models with a size distribution of "patches" that satisfy the relevant requirements ).The blue/green shaded regions with solid boundaries (labeled) show the contours that would produce roughly the correct mean free paths for MeV, GeV, and TeV CRs, obeying all the other constraints in Section 3, for quasi-3D (D = 0) structures.The slope of the blue/green regions is set by fixing the CR mean free path, ℓ mfp ∼ 1/(nsℓ 2 S ), and the maximum allowed size is set by ℓ S > ℓ mfp .An ideal candidate scattering structure, which could explain CR scattering in the ∼ MeV-TeV range, would intercept all of these.The warmer-colored shaded regions with dotted boundaries (labeled) show the approximate location of various known large-scale structures in the ISM including stellar magnetospheres, H II regions, molecular clouds (GMCs), planetary nebulae (PNe), supernova remnants (SNRs), and galactic spiral arms (arms).None of these appear viable: they might scatter CRs, but their abundance is too low to account for most observed CR scattering in the ISM.
(i.e. will produce the same observables) without violating any obvious constraints.

EXAMPLE CANDIDATE STRUCTURES OR PHYSICAL MECHANISMS
We now consider some different physical mechanisms and/or candidate "scattering structures."

Quasi-Static/Coherent and "Macroscopic" ISM Structures
One possibility is that the "patches" of interest could be associated with some known population of quasi-static or coherent ISM structures that are already known to perturb the magnetic field structure on some scale ℓS.In evaluating whether such a population is viable as the dominant source of CR scattering, it is helpful to place them on a sort of modified "Hillas plot" for CR pitch-angle scattering in the ISM, which we show in Fig. 2. There, we plot ns(> ℓS) versus ℓS, for quasi-3D objects (D = 0, which reasonably describes all the systems we consider in the plot), and show the allowed regions that produce the observed mean free paths, combining all of the constraints from Section 3.1-3.1.24for CR protons with energies of ∼ MeV, GeV, and TeV.We also place some rough estimates of the range of number densities and sizes of various "macroscopic" scattering candidates, including molecular clouds (GMCs), stellar magnetospheres, planetary nebulae (PNe), supernova remnants (SNRs), HII regions, and Galactic spiral arms (rough estimates of number this plot is also valid for any dimensionality of structures, D = 0, 1, 2).We compare the observationally and physically-allowed regions, and provide an example (black shaded region) of a hypothetical model for a distribution of scattering patches, with size ℓ S ∼ rgyro and S , which would satisfy all of the observational constraints without over-predicting the CR scattering rate or violating any obvious physical or observational limits.This could arise from intermittent structures in turbulence (Section 4.3), where a very small volume-filling factor or volumetric probability P V ∼ 10 −7 − 10 −5 of structures with size-scale ∼ ℓ S ∼ au featuring O(1) magnetic field fluctuations would be sufficient to explain the observed CR scattering.Interestingly (Section 4.2), the size scales here are very similar to those inferred for small-scale ISM structures responsible for radio-wave plasma scattering, and the volume-filling factors f V might be consistent as well, but orders-of-magnitude uncertainties in the observed f V prevent us from reliably placing specific examples on this plot (likewise for some proposed physical mechanisms in Section 4.4).
An ideal candidate scattering structure would, in this plot, intersect not just one but all three of the allowed CR "bands" without over-predicting the scattering rate for any energy range.Instead, all of the plotted candidates appear clearly ruled out as the dominant source of CR scattering: at a given ℓS, ns is orders-of-magnitude too low.In other words, the mean free path between structures is too large, or, alternatively, the maximum ISM-mean scattering rate, ⟨ν⟩eff ∼ vcr/ℓmfp ∼ vcr ns As ∼ vcr ns ℓ 2 s , is much smaller than required to account for the observed CR scattering.So while these structures can scatter CRs (we know, in fact, the Heliosphere does so), they cannot produce most of the observed scattering.For this reason structures that have similar sizes but are even rarer in the ISM (e.g.Bok globules, pulsar wind nebulae, globular cluster cores, colliding wind binaries) are also immediately ruled out.
Thus while it is conceivable that such a "macroscopic" ISM population might exist, we are unable to identify an obvious candidate.Furthermore, the probability that macroscopic structures act as the primary source of CR scattering is reduced due to the fact that macroscopic structures are predominantly found within the galactic disk.Such a disk-centric distribution of scattering candidates contrasts with the inferred CR confinement several kiloparsecs outside of the galactic disk.

"Small"-Scale ISM Structure
An alternative possibility is that the proposed scattering patches could be associated with some small-scale ISM stuctures or magnetic field features.
While Fig. 2 focused on relatively large-scale structures, in Fig. 3, we repeat the same exercise, but focus on the smaller end of the size range of ℓS (though note the axis ranges do overlap).For small-scale structures, especially where the effective dimensionality may not be known, we find it useful to focus on the volume-filling factor fV ∼ ns Vs ∼ ns As ℓS ∼ ℓS/ℓmfp (Section 3.1.2),rather than ns specifically.This also lets us compress the vertical dynamic range of the plot and factor out the dependence on D, so this plot is valid for sheet-like or filamentary structures, not just quasi-spherical structures.
In this plot, we also show a hypothetical model that would explain the required rigidity dependence of CR scattering and reproduce Solar-system observables fairly naturally, assuming ℓS ∼ rgyro and then calculating ns(> ℓS) such that ℓmfp scales with rgyro and therefore RGV as observed (as in Section 3.3).Since fV ∼ ℓS/ℓmfp we can rearrange to obtain fV ∼ 10 −7 B −δs µG (ν0/10 −9 s −1 ) (ℓS/0.2 au It is noteworthy that the required volume-filling factors ranging from fV ∼ 10 −7 − 10 −5 for structures with sizes ranging from ∼ 0.01 − 100 au (or ∼ a few ×10 −7 around ∼ 1 au) are intriguingly similar to some estimates of the volume-filling factor of so-called "tiny-scale atomic structures" (TSAS; Heiles 1997;Stanimirović et al. 2003;Stanimirović & Zweibel 2018;McEvoy et al. 2015) in the ISM as well as the volume-filling factor estimated in some models of ISM plasma structures causing "extreme scattering events" (ESEs; Romani et al. 1987;Cordes & Lazio 2001;Bannister et al. 2016;Jow et al. 2023).However, we caution that the ESE filling factor is largely unconstrained and model-dependent (Stanimirović & Zweibel 2018).While the sizes (ℓS) of TSAS structures are broadly agreed to lie in the range plotted in Fig. 3, some other observational estimates of their volume-filling factor are as high as ∼ 10 −2 , much larger than what is needed to explain CR scattering (see e.g.Brogan et al. 2005).
Note these are categories of ISM structures classified by their effects on radio waves: physical explanations for such structures range widely, but often invoke intermittent turbulent structures, which we discuss below.

Connection to Intermittency
A natural category of candidates for the scattering structures suggested by Fig. 3 is intermittent turbulent structure (Zhdankin et al. 2016;Mallet & Schekochihin 2017;Dong et al. 2018).Recent test particle simulations of CRs in intermittent magnetic fields suggest that magnetic structure may enhance CR diffusion, even for a fixed magnetic field power spectrum (Shukurov et al. 2017;Seta et al. 2018).The anomalous CR diffusion due to magnetic field intermittency has qualitatively similar behavior to the more general approach of modeling CR diffusion in Fokker-Plank-like equations with non-Markovian statistics (e.g.Wilk & Włodarczyk 1999;Snodin et al. 2016;Zimbardo & Perri 2020).Below, we place constraints on magnetic intermittency in the context of the patchy CR scattering model.
Recall, in the traditional model, we had ⟨ν⟩eff ∼ (vcr/λ) |δB(λ)| 2 /|B| 2 , with the assumption that scattering was dominated by ubiquitous (volume-filling factor fV ∼ 1) but weak (|δB| ≪ |B|) fluctuations, which must obey specific conditions on their power spectra at scales ≪ 100 au in order to reproduce the observed dependence of CR residence time on rigidity.In the patchy model here, we have ⟨ν⟩eff ∼ vcr/ℓmfp ∼ (vcr/ℓS) fV (using fV ∼ ℓS/ℓmfp from Section 3.1.2).So, per Fig. 3, we instead assume that CR scattering is dominated by regions with strong magnetic fluctuations ("patches") but very low volume-filling factor ( fV ≪ 1).For example, if we were to assume gyroresonant λ ∼ ℓS ∼ rgyro, then in order to reproduce the observed ⟨ν⟩eff ∝ R −δs GV , in the "traditional" models we must have |δB(λ)| 2 ∝ λ 1−δs , while in the patchy model we replace this with the requirement fV (ℓS) ∝ ℓ 1−δs S .While the latter does not (and indeed cannot, mathematically) reduce the number of observational requirements on the model, it does avoid all of the mathematical and physical challenges to the "traditional" models.Specifically, the intermittent scattering model removes the requirement that CR scattering theories produce spectra of the form |δB(λ)| 2 ∝ λ 1−δs at sub-100 au scales in the ISM, as reviewed in Hopkins et al. (2022b) and Kempski & Quataert (2022).
In Fig. 3, we show an example of a hypothetical successful intermittent scattering model (black shaded region), assuming the size of scattering structures scales with the CR gyroradius, ℓS ∼ rgyro and δs ∼ 1/2.In order for a scattering patch with a size scale of O(1) gyroradius to reliably scatter CRs (have a "scattering optical depth" of order unity) at that rigidity, this hypothetical model requires a magnetic fluctuation amplitude O(|δB(ℓS)|/|B|) ∼ 1.So, as heuristically demonstrated in Fig. 3, a model that features intermittent structures with O(|δB(ℓS)|/|B|) ∼ 1 on size scales ℓS ∼ 1 − 1000 au, with small volume-filling factor fV ∼ 10 −6 (ℓS/10 au) 1/2 , would automatically give rise to the "desired" (empirically-inferred) CR scattering rates at ≲ GeV through ≳ TeV energies.The exact scaling, of course, would also depend on the details of how the magnetic field strength scales with the size of scattering patches (e.g.Lemoine 2023;Kempski et al. 2023).
We can also reason about the prevalence of intermittent turbulent structures in terms of the shape of the probability distribution function (PDF) of magnetic fluctuations.Consider the volumetric PDF PV (δB | ℓS) of fluctuations |δB| with a given size scale/wavelength ℓS: to calculate the contribution of different fluctuations with the given ℓS to CR scattering we should integrate over this PDF.If the PDF is Gaussian/normal -i.e.completely nonintermittent -then the contribution to CR scattering will be dominated by the "weak" ±1σ fluctuations in the core of the PDF, giving rise to the usual scattering rate ∝ ⟨|δB| 2 ⟩ ∼ |δB| 2 median .But now, as is commonly parameterized for intermittent systems, consider a PDF with power-law tails in the rare-event (large-δB) regime, dPV /d ln |δB| ∝ |δB| −α P (i.e.dPV /dδB ∝ |δB| −α P −1 ).The critical division between the "patchy" and "traditional" behaviors will then occur at αP = 2.If the PDF of the magnetic fluctuations that scatter CRs falls more steeply (αP > 2), then the "core" of the PDF dominates CR scattering and the behavior will resemble the traditional model.If the PDF falls more slowly/is more shallow (αP < 2), then the contributions to CR scattering will be instead dominated by the largest (non-linear, O(1)) fluctuations in the tails -our "patchy" behavior.

Candidate Mechanisms for Intermittent, Small-Scale ISM Structure
Briefly, we note that both the empirically-observed small-scale scattering structures from Section 4.2 and the theoretical category of "intermittent" structures from Section 4.3 could be related to a variety of distinct microphysical processes in the ISM.Each of these is, in a sense, a "candidate" for the patchy scattering structures.This includes plasma sheets in MHD turbulence (Dong et al. 2022), turbulent boundary/mixing layers (Ji et al. 2019;Yang & Ji 2023), magnetic mirrors and traps (Chandran 2000b;Bustard & Zweibel 2021;Lazarian & Xu 2021;Tharakkal et al. 2023), plasmoid instabilities (Fielding et al. 2023), weak shocks (Kadomtsev & Petviashvili 1973;Makwana & Yan 2020;Kempski & Quataert 2022), regions with strong dust-CR coupling (Squire et al. 2021;Ji et al. 2022), and regions where self-confinement has "run away" (e.g. in "staircase-like" instabilities; Quataert et al. 2022b;Huang & Davis 2022;Tsung et al. 2022), to name a few.CR scattering by rare (i.e.not volume-filling) regions of large field-line curvature, proposed recently by Lemoine (2023) and Kempski et al. (2023), is a promising microphysical origin of the patchy transport model.Small-scale bends of magnetic field lines may be a generic intermittent feature of MHD turbulence and exist even on scales much smaller than the turbulence injection scale.These bends are therefore plausible candidates for the "patches" in our patchy transport model.The scattering rates derived in Lemoine (2023) and Kempski et al. (2023) depend on the volume-filling fractions of the large-curvature patches, broadly consistent with the calculations presented here.In particular, our result that ℓmfp ∼ ℓS/ fV is somewhat similar to equation (3.1.)in Lemoine (2023).
All of these scenarios appear viable in principle, but none has reached the level of theoretical development where we can plot an unambiguous prediction in Fig. 3.However, our hope is that in future studies modeling the structures predicted by such mechanisms, Fig. 3 can prove a useful "figure of merit" to test whether such mechanisms are (or are not) viable candidates for explaining observed CR scattering in the ISM.

OBSERVATIONAL TESTS
A key property of the patchy CR scattering models discussed here is that they, by definition, produce the same observational constraints and phenomenologically-derived CR properties (e.g.CR scattering rate ⟨ν⟩eff, scattering/deflection mean free path ℓmfp, CR residence time, grammage, etc.) as "traditional" continuous CR scattering models.So long as the CR residence time (observationally inferred to be ≳ Myr) is longer than the scattering time ∼ 1/⟨ν⟩eff (which for the example observationally-inferred values quoted in Section 3 is ∼ 30 yr at ∼ 1 GeV), or equivalently so long as the total "path length" traveled along the trajectory of a given CR is larger than one mean free path (equivalent to saying that the grammage exceeds X > ρISM ℓmfp ∼ 10 −5 g cm −2 R −0.5 GV , which is satisfied by several orders of magnitude), this ensures all of the predictions for e.g.CR spectra, primary-to-secondary or radioactive or isotopic ratios, etc., will be identical in a "patchy" model and a "continuous/traditional" model with the same ⟨ν⟩eff.This means, for example, that while one might naively expect a "patchy" model to produce a larger observable CR anisotropy, such a difference would only be present in the population at distances ≲ ℓmfp ≲ 10 pc from the initial CR sources, and the CRs will (by definition) converge to isotropic on a timescale ∼ 1/⟨ν⟩eff ∼ ℓmfp/vcr ∼ 30 yr.Since no SNR exists so close to Earth, we predict no difference in the anisotropy of local CRs.
More detailed tests of e.g.higher-order statistics and correlations will, in general, require a specific model for the origin of the "patchy" CR scattering (e.g.different specific physical scenarios discussed in Section 4).However, we can propose some generic predictions that may conceivably be testable with future instruments.
Consider, for example, a scattering patch of depth ℓS with internal scattering rate νs, and a "CR scattering optical depth" of τs ∼ ℓsνs/vcr (Section 3.1).By definition, this patch can scatter, i.e. temporarily confine, CRs below some critical rigidity Rs, enhancing their density relative to the ambient medium by a factor ∼ τs (just like with multiply-scattered photons).This would, in principle, enhance the γ-ray emission at some energies from the patch.However, that effect alone would be strictly degenerate with local variations in the source density, variations in the volume-filling scattering rate in the "traditional" models, and variations in the background nucleon number density of the ISM.But since, by definition, CRs with rigidities R ≫ Rs are not confined by the patch, their relative density is not enhanced.Therefore, the patches would exhibit a γ-ray spectrum that is more biased to emission from R < Rs compared to the ambient medium.Not only would such variation exist, but the change in spectral shape would be correlated with the sizes and number densities of the scattering structures as we have discussed above (because we require a spectrum of patches that scatter CRs of different energies with different relative rates).
The challenge here is that measuring such an effect would require orders-of-magnitude superior resolution in γ-ray telescopes compared to current state-of-the art instruments like Fermi.Ideally, observing this effect would require the ability to resolve the γ-ray spectra at ∼ 0.1 − 100 GeV energies with spatial resolution ∼ ℓS ∼ au (or at least ≪ 10 3 au) -i.e.angular resolution in the ISM (for typical spatial distances) reaching sub-arcsecond levels, whereas current instruments typically achieve few-degree resolution at these energies.
Alternatively, one could look for a similar effect in radio synchrotron emission from CR electrons with similar energies, where the angular resolution of current instruments is much less limited.But even in this case, the desired angular resolution is still a challenge, and far beyond the scope of current single-dish surveys, requiring interferometry with ≫ km baselines.More problematically, the synchrotron spectral slope at these energies is also strongly influenced by the relative strength of inverse Compton and synchrotron losses, which would be at least partially degenerate with the desired measurement.Still, despite these challenges, there may already be hints of the relative/patchy CR enhancement described above in existing observations, perhaps related to inferences of an excess of low-energy ionizing CRs in GMCs (e.g.Indriolo et al. 2009;Indriolo & McCall 2012;Indriolo et al. 2015;Yang et al. 2014;Baghmanyan et al. 2020), or to the observed "point-sourcelike" excess of ∼ 1 GeV (aka "soft" according to the arguments above) γ-ray emission from the Galactic center (Hooper & Goodenough 2011;Cholis et al. 2015;Lee et al. 2016;Bartels et al. 2016;Ackermann et al. 2017;Hooper et al. 2017).
Another possible set of tests would be to look for the variations caused by discreteness noise in the "number of scattering structures" on scales comparable to ℓmfp around specific CR sources.We caution that the point of interest here is not the emission from the acceleration region or CR source itself (e.g.not emission from SNRs or PWNe), but the secondary emission from the ISM on ∼ pc scales around such a source.Here one must be careful to estimate out to which radii the CR background "seen by" the ISM would be dominated by the local source, versus the collective Galactic background.Moreover, one has the same resolution challenges as noted above, since the most obvious tests would require resolving inhomogeneity on a spatial scale of order ∼ ℓS.

SUMMARY
In light of the theoretical challenges facing "traditional" CR scattering theories, which assume low-energy (∼ MeV-TeV) CRs are scattered by a uniform, volume-filling population of weak magnetic field fluctuations, we propose a novel category of "patchy" CR scattering models, in which CRs are scattered by inhomogeneous/intermittent/punctuated structures with small volume-filling factors.We derive a number of constraints any such structures must obey (e.g.their sizes, internal CR scattering rates/magnetic field fluctuations, number densities, mass and volume-filling factors) in order to be both internally self-consistent and to reproduce existing observational constraints.We show that fundamentally, in this category of models, we can reproduce the observed dependence between CR scattering and rigidity (⟨ν⟩eff ∼ R −δs GV ) by imposing a size distribution of scattering structures: larger "patches" (with greater optical depth to higher-energy CRs) are rarer but scatter a wider range of CR energies, while smaller "patches" are more abundant but only scatter lower-energy CRs.
We consider a variety of observational and physical candidates for such structures.We show that many "macroscopic" quasi-static ISM structures (e.g.GMCs, SNRs, PNe, stellar magnetospheres, HII regions, PWNe) cannot be the dominant scattering regions since the mean-free path between them, as seen by CRs, is too large.However we show that small-scale or intermittent structures in the ISM, with size scales as small as ≲ au and volume-filling factors as small as ∼ 10 −7 , could plausibly explain the observed CR scattering rates from ∼ MeV-TeV energies, while avoiding any obvious theoretical or observational objections.These may even be related to other small-scale ISM structures inferred from radio-wave scattering observations.However, as of yet, there is no single obvious physical mechanism that is clearly predicted to produce the desired scattering rates.We propose a "figure of merit," akin to a "Hillas plot" for CR pitch-angle scattering in the ISM, with which to compare any such future models.
We discuss some direct observational tests of this model category, though we conclude that, for now, the required resolution remains far beyond current capabilities.

Figure 2 .
Figure2.Constraints on the size, ℓ S , and volume number density in the ISM, ns, of possible patchy scattering structures (Section 4.1).The blue/green shaded regions with solid boundaries (labeled) show the contours that would produce roughly the correct mean free paths for MeV, GeV, and TeV CRs, obeying all the other constraints in Section 3, for quasi-3D (D = 0) structures.The slope of the blue/green regions is set by fixing the CR mean free path, ℓ mfp ∼ 1/(nsℓ 2 S ), and the maximum allowed size is set by ℓ S > ℓ mfp .An ideal candidate scattering structure, which could explain CR scattering in the ∼ MeV-TeV range, would intercept all of these.The warmer-colored shaded regions with dotted boundaries (labeled) show the approximate location of various known large-scale structures in the ISM including stellar magnetospheres, H II regions, molecular clouds (GMCs), planetary nebulae (PNe), supernova remnants (SNRs), and galactic spiral arms (arms).None of these appear viable: they might scatter CRs, but their abundance is too low to account for most observed CR scattering in the ISM.

Figure 3 .
Figure3.Similar to Fig.2, we constrain the minor axis, ℓ S , and volumefilling factor, f V , required of the patchy scattering model (Section 4.2-4.4;this plot is also valid for any dimensionality of structures, D = 0, 1, 2).We compare the observationally and physically-allowed regions, and provide an example (black shaded region) of a hypothetical model for a distribution of scattering patches, with size ℓ S ∼ rgyro and f V ∼ ℓ S /ℓ mfp ∼ ℓ Farmer & Goldreich 2004;Zweibel 2017b;Squire et al. 2021én waves, eB ≡ |B| 2 /8π is the magnetic energy density, vA is the Alfvén speed, and ecr is the CR energy density (equation 7 inHopkins et al. 2021c).Following the assumptions in Hopkins et al. 2021c, we can turn the above equation into a rough, order-ofmagnitude lower limit on the required damping rate to sufficiently suppress the CRSI in the volume-filling ISM, Γdamp ≳ 10 −9 s −1 , a rate on the higher end of that obtained via simple estimates in the ionized ISM, though potentially reasonable (e.g.Farmer & Goldreich 2004;Zweibel 2017b;Squire et al. 2021).