Cosmic ray feedback in the FIRE simulations: constraining cosmic ray propagation with GeV gamma ray emission
T. K. Chan D. Keres P. F. Hopkins E. Quataert K. -Y. Su C. C. Hayward C. -A. Faucher-Giguere



Cosmic Rays in FIRE-2]Cosmic ray feedback in the FIRE simulations: constraining cosmic ray propagation with GeV gamma ray emission

T. K. Chan et al.]
  T. K.  Chan1 Email: (TKC) tkc004@physics.ucsd.edu,
  D.  Kereš1 Email: (DK) dkeres@physics.ucsd.edu,
  P. F. Hopkins2,
  E.  Quataert3, K.-Y. Su2,C. C.  Hayward4, & C.-A. Faucher-Giguère5 

  1Department of Physics, Center for Astrophysics and Space Sciences,University of California at San Diego,
 9500 Gilman Drive, La Jolla, CA 92093, USA

  2TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA 

  3Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720, USA

  4Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA

  5Department of Physics and Astronomy and CIERA, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA



[
    [
    
=====




We present the implementation and the first results of cosmic ray (CR) feedback in the Feedback In Realistic Environments (FIRE) simulations. We investigate CR feedback in non-cosmological simulations of dwarf, sub-L⋆ starburst, and L⋆ galaxies with different propagation models, including advection, isotropic and anisotropic diffusion, and streaming along field lines with different transport coefficients. We simulate CR diffusion and streaming simultaneously in galaxies with high resolution, using a two moment method. We forward-model and compare to observations of γ-ray emission from nearby and starburst galaxies. We reproduce the γ-ray observations of dwarf and L⋆ galaxies with constant isotropic diffusion coefficient κ∼ 3× 10^29  cm^2 s^-1. Advection-only and streaming-only models produce order-of-magnitude too large γ-ray luminosities in dwarf and L⋆ galaxies. We show that in models that match the γ-ray observations, most CRs escape low-gas-density galaxies (e.g. dwarfs) before significant collisional losses, while starburst galaxies are CR proton calorimeters. While adiabatic losses can be significant, they occur only after CRs escape galaxies, so they are only of secondary importance for γ-ray emissivities. Models where CRs are “trapped” in the star-forming disk have lower star formation efficiency, but these models are ruled out by γ-ray observations. For models with constant κ that match the γ-ray observations, CRs form extended halos with scale heights of several kpc to several tens of kpc.




galaxies: evolution — galaxies:  kinematics and dynamics






§ INTRODUCTION




Cosmic rays (CRs) are charged particles moving with relativistic speeds, mainly generated through shock acceleration of supernova remnants (SNRs)  <cit.> (and possibly also in active galactic nuclei in massive galaxies). Unlike thermal energy, they can propagate through the galactic interstellar medium (ISM) rapidly via advection, diffusion and streaming <cit.>, and transfer energy to gas via Coulombic and hadronic interactions <cit.>. Their energy density is comparable to thermal and magnetic energies in the solar neighborhood <cit.>, so CRs are believed to be dynamically important in galaxy evolution.

The impacts of CRs on galaxy evolution have been studied with analytic models <cit.> and idealized and cosmological simulations (e.g. ). These studies show CRs can drive multiphase winds, reduce star formation rates in low mass galaxies, thicken gaseous disks, and modify the phase structure of the circum-galactic medium (CGM). It has also been suggested that CRs may play an important role in the galactic dynamo <cit.>. 

Despite its importance, the details of CR propagation are uncertain. The most popular CR propagation models are self confinement and extrinsic turbulence (, and reference herein). In the former picture, CRs interact with a series of Alfven waves, which results in random scattering in pitch angles. The waves are then amplified via the streaming instability of CRs, increasing the scattering and trapping CRs in a background medium. These “self confinement” interactions effectively transfer energy from the CRs to thermal plasma. In the extrinsic turbulence model, CRs propagate through random magnetic field lines and are scattered by the background turbulent magnetic fields. This mechanism is especially important for high energy CRs, since there are too few high energy CRs amplifying the Alfven waves and the self-confinement mechanism fails <cit.>. These two mechanisms confine and isotropise the CR distribution explaining the remarkably low CR anisotropy observed from the Earth (see, e.g. ) and the long residence time (>10^7  yr) inferred from the ratios between stable primary and secondary nuclei <cit.>. Their long confinement time and small anisotropy imply that CRs have short mean free paths (∼ pc) and their propagation can therefore be approximated by a random walk, so CRs can be well described as a diffusive fluid, obeying an advection diffusion equation (see, e.g. , for arguments for the CR fluid theory).
 
Most of the studies of CR propagation have focused on an approximate picture of the Milky Way described by the flat halo diffusion model <cit.>. This model consists of a cylindrical gaseous halo with a radius around 20 kpc and a height larger than 1 kpc, and a thinner yet more dense cylindrical internal disk with CR sources. CRs are assumed to diffuse isotropically (averaged over the scale of hundreds of pc) with a spatially constant but energy dependent diffusion coefficient, and “escape” through the halo boundaries to intergalactic space. Extensions of this model are commonly used in numerical CR propagation codes, e.g. GALPROP[ <https://galprop.stanford.edu/>], which attempt to synthesize observational constraints on the MW gas density distribution, CR abundances and spectra, γ-ray and radio emission, and theoretical models for e.g. galactic winds and diffusive re-acceleration <cit.>. These models are commonly used to constrain the “effective” isotropic-equivalent diffusion coefficient of CRs averaged over the whole MW disk <cit.>. However, there are still large uncertainties in the role of gas dynamics and small-scale gas density fluctuations (“clumpiness”), magnetic field geometries on small scales, the spatial and temporal distribution of CR sources, the size and mass distribution of the gaseous galactic halo, and the CR propagation model. To make progress, self-consistent modeling of galaxy evolution that includes CR propagation together with hydrodynamics or magneto-hydrodynamics is required.




In addition to the CR energy density and abundance of nuclei, high energy γ-ray emission can serve as an independent constraint on CR propagation <cit.>. High energy (> GeV) CRs collide with nuclei in the interstellar medium (ISM) and produce pions, which decay into GeV γ-rays. Since pionic γ-ray luminosity is proportional to CR density and most of the energy density of the CRs is at energies around GeV (from the direct CR observations, e.g. in ), CR distribution and propagation can be constrained with high energy γ-ray observations. 


Recently, γ-ray emission was observed from Local Group <cit.> and starburst galaxies <cit.>, which can be used as a probe of CRs beyond the solar system and the Milky Way (MW). <cit.> found a correlation between γ-ray emission and star formation rate (SFR) with a steeper than linear drop at low SFRs (E_γ∝ SFR^1.4± 0.3; summarized in ).


To explain this correlation, <cit.> (hereafter ) constructed one-zone leaky box model of galaxies where a fixed fraction of SN energy is injected as CRs. They assume CRs escape with an energy dependent escape time and that the CR energy density and spectral distribution are in a steady state (the injections and losses are balanced).  Constrained with the observed CR abundances and the far infrared (FIR)-radio correlation <cit.>, the model was used to estimate pionic γ-ray luminosities of galaxies. They found that in order to explain the correlation between γ-ray emission and SFR, in starburst galaxies, most CR protons are required to lose their energy via collisions with the ISM (i.e. that starbursts are “CR proton calorimeters”, as in the earlier calculations of ; see also ), while in dwarf galaxies, most of CR protons should escape. The main drivers of this effect are that SFR drops with decreasing gas surface density <cit.>, and that lower gas densities enable CRs to escape before heavy pionic losses. Subsequent observational studies have reached the same conclusion regarding efficient escape in galaxies like the MW, Andromeda (M31), the Large and Small Magellanic Clouds (LMC, SMC; see, e.g. ).

In this study, we investigate the impact of CRs on dwarf, sub-L⋆ starburst, and L_⋆ galaxies, using idealized simulations of galaxy evolution. We run galaxy simulations with both CR diffusion and streaming with high spatial resolution and diffusivity thanks to the newly developed two-moment method <cit.>. We also couple explicit CR transport and losses to an explicit, local stellar feedback model which time-resolves individual SNe explosions, as well as stellar mass-loss and radiative feedback, which together enable self-consistent generation of galactic winds and a turbulent, multi-phase ISM, critical for understanding CR transport and emission in that same ISM. Specifically, our CR implementation in the code GIZMO is coupled to the FIRE-2 (Feedback In Realistic Environments 2) algorithm for star formation and stellar feedback <cit.>.[<http://fire.northwestern.edu/>] Cosmological simulations with these physics (without explicit CR transport) have been shown to successfully reproduce many observed galaxy properties including stellar masses <cit.>, galactic winds <cit.>, cored central dark matter profiles <cit.>, the mass-metallicity relation <cit.> and spatial distribution of gas and metals within galaxies and the CGM <cit.>, typical galaxy star formation rates and histories <cit.>, the Kennicutt-Schmidt law <cit.>, and galactic magnetic field structure <cit.>. 

In this, the first paper in a series, we introduce our implementation of the CR propagation model (including isotropic and anisotropic diffusion and streaming), simulate galaxies with several CR propagation models, and focus on constraining the model using the observations of ∼GeV  γ-ray emission from galaxies (and compare our findings with previous theoretical studies).  <ref> and  <ref> discuss numerical methods, simulated physics, and initial conditions. In  <ref>, we investigate how CRs and their propagation influence galactic properties. In  <ref> we calculate the γ-ray emission from CRs in our simulations and compare with observational data. In  <ref>, we compare our findings with the previous studies and analyze the relative importances of different CR energy gain and loss processes. We summarize our findings in  <ref>. 



§ METHOD




 §.§ Simulation code



All the physics and numerical details in this study, except for CRs, follow the FIRE-2 version of the FIRE algorithms presented in detail in <cit.>, so we only briefly review them here. Our simulations use the GIZMO[<http://www.tapir.caltech.edu/p̃hopkins/Site/GIZMO>] code <cit.> in its mesh-free Lagrangian finite mass (MFM) mode for (magneto)-hydrodynamics; extensive implementation details and tests of the MHD scheme are presented in <cit.>. GIZMO uses an updated version of the PM+Tree algorithm from Gadget-3 <cit.> to calculate gravity and adopts fully conservative adaptive gravitational softening for gas <cit.>. Gas cooling is calculated with tabulated cooling rates from T=10-10^10 K, including atomic, metal-line, and molecular cooling. While our simulations are idealized and do not include cosmological environments, we do include the present-day  ultraviolet background, from the <cit.> model (accounting for local self-shielding). Stars form in locally self-gravitating, self-shielding, thermally Jeans-unstable gas[We assume the strong coupling limit between gas and CRs, so the effective sound speed in the virial parameter includes both thermal and CR pressure. See Appendix C in .] at densities  n_ H≥ 100  cm^-3. Once formed, we calculate the energy, momentum, mass and metal return for each star according to the STARBURST99 stellar population synthesis model <cit.>, for a <cit.> IMF, accounting for SNe Types Ia & II, O/B and AGB star mass-loss, and radiation (photo-electric and photo-ionization heating and radiation pressure with a five-band approximate radiation-hydrodynamics treatment). For details see <cit.>.





 §.§ Cosmic Rays


The implementation of CR physics in GIZMO includes fully-anisotropic cosmic ray transport with streaming and advection/diffusion, CR cooling (hadronic and Compton, adiabatic, and streaming losses), injection in SNe shocks, and CR-gas coupling. The CRs are treated as an ultra-relativistic fluid (adiabatic index γ_ cr=4/3) in a “single bin” approximation.[One can think of this as evolving only the CR energy density at the energies ≳ GeV, which dominate the CR pressure, and approximating the CR energy spectrum as having a universal shape at all positions.] Integrating over the CR distribution function and spectrum, the usual ideal-MHD equations solved for gas density ρ, velocity v, magnetic field B, and specific energy e, are extended with the equation for the CR energy density e_ cr <cit.>: 

    ∂ρ/∂ t+∇· ( ρ v )=0,
       ∂ρ v/∂ t+∇· ( ρ v⊗ v+P_T𝕀- B⊗ B )=0,
       ∂ρ e/∂ t+∇· [ (ρ e+P_T) v-( v· B) B ]
                  =P_ cr∇· v+Γ_ st+S_ g-Γ_ g,
       ∂ B/∂ t+∇· ( v⊗ B- B⊗ v )=0,
       ∂ e_ cr/∂ t + ∇· F_ cr =  v·∇ P_ cr- Γ_ st+ S_ cr - Γ_ cr,

where P_ cr = (γ_ cr-1) e_ cr is the CR pressure; P_T is the total (thermal+magnetic+CR) pressure; Γ_ st=- v_ st·∇ P_ cr is the CR “streaming loss term” discussed below; S_ g and S_ cr are gas and CR source terms (e.g. injection); Γ_ g and Γ_ cr are gas and CR sink/loss (or “cooling”) terms; and v_ st is the CR streaming velocity. F_ cr is the CR energy flux, which can be written F_ cr = (e_ cr + P_ cr)( v+ v_ st)+  F_ di  where the first term represents advection and streaming, whereas the second term is a diffusive-like flux (e.g. given by F_ di= -κB̂⊗B̂·∇ e_ cr in the “pure diffusion” or “zeroth moment” approximation, but we explicitly evolve this; see <ref>). 

For the gas equations-of-motion, note when solving the Riemann problem between neighboring fluid elements, P_T includes the CR pressure (i.e. we make the local strong-coupling approximation: CRs and gas strongly interact), and the effective sound speed of the two-fluid mixture is modified to (c_s^2)_ eff = ∂ P/∂ρ = (c_s^2)_ gas + γ_ cr P_ cr/ρ, but no other modifications to the MHD method is required. 



  §.§.§ CR Transport: Advection & Streaming


In our method, each mesh-generating point (which defines the gas resolution “elements”) represents a finite-volume domain that moves with the fluid velocity v= v_ gas in a quasi-Lagrangian fashion. After operator-splitting the source/injection and loss/cooling terms, it is convenient to re-write the advection and streaming terms in the following Lagrangian, finite-volume form (see e.g. ): 

    D E_ cr/D t = -∫_Ω d^3 x {   P_ cr   (∇· v) +Γ_ st +  ∇·F̃_ cr}

where D/Dt = ∂/∂ t +  v·∇ is the Lagrangian derivative co-moving with the gas, and E_ cr^i = ∫_Ω_i e_ cr d^3 x is the conserved total CR energy in the finite-volume domain Ω_i belonging to element i. Here F̃_ cr≡ F_ cr -  v (e_ cr+P_ cr) =  v_ st (e_ cr+P_ cr) +  F_ di. Pure advection with the gas is automatically handled in this description. In cosmological simulations, the Hubble flow is included in ∇· v. 

The P_ cr(∇· v) term represents adiabatic changes to the CR energy via compression/expansion (the “PdV work”), which exchanges energy with gas. We will refer to this as the “adiabatic” term throughout.[To ensure manifest energy conservation, this is solved when the mesh positions are updated. Calculating the volume changes Δ V_i = ∫ dt ∫_Ω_i d^3 x (∇· v) with the kernel-weighted divergence of the fluid velocity field (which is the exact discrete change in the domain volume as defined in )], we have Δ E_ cr = -P_ cr Δ V_i. This is removed from the total energy equation after the hydrodynamic Riemann problem is solved to determine the total gas energy update.


The Γ_ st=- v_ st·∇P_ cr term represents “streaming loss”, which transfers energy to gas and is always positive because CRs always stream down the CR pressure gradient (see the next section). As CRs stream, instabilities excite high-frequency Alfven waves (frequency of order the gyro frequency, well below our simulation resolution limits; see e.g. ) which are damped and thermalize their energy effectively instantly (compared to our simulation timescales).[With the streaming velocity defined below, the streaming loss term can be written D E_ cr / D t = - E_ cr / τ_ st with τ_ st^-1 = (γ_ cr-1) |B̂·∇̂ e_ cr|^2 |v_ st ∇ e_ cr|/e_ cr. When this is updated the resulting energy lost Δ E_ cr = ∫ dt τ_ st^-1 E_ cr is added to the gas thermal energy.]

Finally, the ∫_Ω d^3 x ∇·F̃_ cr term does not change the total CR energy, but represents flux of energy between resolution elements, caused by CR streaming and diffusion. This can be transformed via Stokes's law into a surface integral, ∫_∂Ω d A·F̃_ cr, which is then solved via our usual second order-accurate, finite-volume Godunov MFM method (in a manner identical to the hydrodynamic equations, see  for details). 

We explicitly evolve the conserved quantities E_ cr and total gas energy E_ gas which are exchanged (either between gas elements or one another), ensuring manifest total energy conservation.[Because we do not evolve a total energy equation, if we use adaptive timesteps, total energy conservation is formally exact at integration-error level rather than machine-accurate. However we have verified that the errors are typically small (percents-level over hundreds of millions of years evolution, although in the most extreme case we find the cumulative amount over 500 Myr can be ≲ 20% of the injection), and negligible compared to physical non-conservative terms (e.g. collisional/radiative losses, injection).] 





  §.§.§ The Streaming Velocity


CRs stream at some speed v_ st down the local CR phase-space density gradient (which is equivalent in our single-bin approximation to CR pressure or energy density gradient), projected along the magnetic field lines, i.e. v_ st≡ -v_ st B̂ (B̂·∇̂ P_ cr) where ∇̂ P_ cr = ∇̂ e_ cr = (∇ P_ cr)/|∇ P_ cr| = (∇ e_ cr)/|∇ e_ cr| is the direction of the CR pressure/energy density gradient. 

It is generally believed that micro-scale instabilities limit the streaming velocity to Alfven speed v_ A (=B/√(4πρ)) in the low-β limit (see , or more recently ), so we adopt a fiducial value v_ st = v_ A.

But in the weak-field, plasma β≫ 1, regime, the streaming velocity can be boosted by significant wave damping <cit.>, so we have also tested various streaming speeds in Appendix <ref>. Although the streaming velocity can in principle exceed v_A by a large factor, <cit.> and <cit.> argued that the streaming loss Γ_ st should be still limited by Γ_ st=- v_ A·∇ P_ cr, because this term is mediated by the excitation of Alfven waves. Therefore, regardless of streaming velocities, we set the streaming loss to - v_ A·∇ P_ cr. When streaming is disabled we simply eliminate terms relevant to streaming.




  §.§.§ Diffusive Transport Terms: Two-moment Method


It is common in the literature to treat F_ di in the “zero-th moment” expansion, i.e. approximate it as an anisotropic scalar diffusion with F_ di = -κ B̂ (B̂·∇ e_ cr), where κ is the effective diffusion coefficient,  which parameterizes the unresolved CR physics. However at high resolution this is problematic for several reasons: (1) it imposes a quadratic timestep criterion (if evaluated with an explicit scheme: Δ t < C_ cour Δ x^2/κ, where Δ x is the resolution and C_ cour the Courant factor) which becomes very small; (2) it implies unphysical super-luminal CR transport when the gradient-scale length e_ cr/|∇ e_ cr| becomes smaller than κ/c ∼ 3  pc (κ / 3× 10^29  cm^2 s^-1) (resolution often reached for simulations in this paper); (3) it cannot smoothly handle the transition between streaming and diffusion-dominated regimes; (4) it will develop spurious numerical oscillations near extrema when handling streaming <cit.>; and (5) it encounters the usual difficulties with anisotropic diffusion operators in moving-mesh codes described in <cit.> (including e.g. difficulty if CRs are “trapped” at local maxima). 

Hence, we adopt a simple two-moment approximation for CR diffusion and streaming, independently developed for this work but similar in concept to the recently-presented implementations in <cit.> and <cit.> (although the concept and use in CR dynamics are well-established; see  for examples). Rather than set F̃_ cr=-κ ∇ e_ cr, we explicitly evolve the flux equation:

    1/c̃^2 [ ∂F̃_ cr/∂ t + ∇·(  v⊗F̃_ cr) ] + ∇_ P_ cr   = -(γ_ cr-1)/κ^∗ F̃_ cr,

where F̃_ cr is the CR flux measured in the frame comoving with the fluid, ∇_P_ cr≡B̂⊗B̂·∇ P_ cr, c̃ is the (reduced) speed of light, and κ^* is the composite parallel (magnetic field-aligned) diffusion coefficient in the rest frame of the fluid, 

    κ^*=κ+v_ st(e_ cr+P_ cr)/|B̂·∇ e_ cr|,

where the second term includes the CR streaming with the streaming velocity specified above. 


For the numerical implementation of CR energy and flux, we follow the treatment of diffusion operators in MFM, outlined in Section 2 in <cit.> with a few modifications.

We solve a general evolution equation of conserved quantities ( VU )_i of cell i (e.g. CR energy) with a volume V by summing over all adjacent cells j:

    d/d t ( VU )_i=-∑_j F^*_ diff,ij· A_i,j,

where V is the cell volume and F^*_ diff,ij is the interface value of the flux and A_i,j is the effective face area defined in <cit.>, in the following steps:


  
  * We calculate all relevant coefficients, using the standard gradient estimator in GIZMO for MFM to estimate gradients, e.g. ∇ P_ cr and ∇ e_ cr, as described in <cit.>;  
  
  * We estimate the values on the left and right sides of the face from the values of cells i and j through a linear reconstruction and use them to solve the Riemann problem; 
  
  * We compute the interface face value of the flux F^*_ diff,ij by solving the Riemann problem (RP) through the <cit.> (HLL) method, using a MINMOD slope limiter, in an operator split manner from the pure MHD. 

  
  
  * Finally, the source term in Eq. <ref> (not considered in ) is added implicitly to ensure stability.

  
  We differ from <cit.> since (1) we explicitly evolve the CR flux (instead of calculating it from the CR energy gradient) and (2) in the Riemann solver, we consider the fastest wavespeed to be generally c̃  (since we choose c̃ to be faster than other physical processes).
  

Unlike <cit.>, we do not modify the momentum transfer from CRs to gas, i.e. the second line of Eq. <ref>, since we assume the strong coupling limit between gas and CRs throughout the paper. This, however, will over-estimate the momentum transfer from CRs to gas when the strong coupling assumption breaks down, i.e. in regimes where the CR and gas coupling is weak and the CR mean free paths are long. In all of our simulations, the mean free paths ∼ (1  pc)(κ/10^29 cm^2/s) are smaller than the resolved length-scales, so the strong coupling assumption is probably relevant in most of the situations.[ The formulation in <cit.> will also over-estimate the momentum transfer in the weakly coupling regime if the “reduced speed of light” approximation is introduced (see below and  5.2 in ).] 

We stress that while the flux equation can be generically obtained by integrating the first moment of the Vlasov equation (with some model for closure of higher  moments, equation-of-state, and scattering terms), one should not take Eq. <ref> to represent a physical two-moment expansion of the relativistic Vlasov equation for CRs. Doing so requires making a number of additional assumptions about e.g. the CR phase space distribution function, ratio of gyro radii to resolved scales, and order of truncation in terms 𝒪(v/c). We discuss some of the subtle differences that can arise in Eq. <ref> as a result, in Appendix <ref>, but stress that on large and/or long time scales these vanish, and so they have no effect on our conclusions in this paper. For our purposes here, it is better to think of it as a generic two-moment numerical expansion of the anisotropic diffusion + streaming equation which eliminates all of the numerical pathologies (1)-(5) above. In future work, it will be interesting to explore more detailed physically-derived transport models that include these higher-order terms, and attempt to actually predict the coefficients κ and v_ st on physical grounds <cit.>.

For now, if we ignore streaming, we see that in steady-state and/or when c̃ is large, or Δ t ≫κ/c̃^2 (or on spatial scales ≫κ/c̃), this equation becomes F_ cr≈ -κ B̂ (B̂·∇ e_ cr), and we recover the usual diffusion equation (see Appendix <ref> for a comparison between the pure diffusion and two-moment methods). However, the two-moment method smoothly limits the maximum bulk transport velocity of the CRs to c̃, and makes the timestep criterion Δ t < C_ cour Δ x / c̃,[We adopt C_ cour=0.4 throughout, and have validated stability (as expected) for this value. Δ x in the Courant condition is defined in the same manner as <cit.> as the local mean inter-particle separation (i.e. the equivalent of the grid spacing in a regular-grid code), Δ x ≡ (m/ρ)^1/3.] which is only first-order, instead of quadratic, in Δ x. 

For true micro-physical CR motion, however, c̃≈ c, the speed of light, which still requires a impractically small timestep. Fortunately, for our purposes in these simulations – where we only capture bulk CR properties in the fluid limit – it is more convenient to consider c̃≪ c (namely the “reduced speed of light” approximation), since galaxy properties should still converge, regardless of c̃, provided it is set to some value faster than other relevant physical processes, e.g. the CR cooling or the actual bulk flow speeds realized in our simulations. We have experimented extensively with this and find that, for the simulations here, values c̃∼ 500-2000  km s^-1 are sufficient to give converged results, e.g. SFR and γ-ray emission (see Appendix <ref>). 

In Appendix <ref>, we compare the results using the simpler pure-diffusion (zeroth-moment) approximation: we then simply assume F_ di→ -κ B̂ (B̂·∇ e_ cr) and solve the anisotropic diffusion equation (with the stricter Courant condition) as described in <cit.>. This is equivalent to adopting c̃→∞, in our Eq. <ref>. For  the same κ_∗, this gives nearly-identical results to our default Eq. <ref> in our galaxy simulations, demonstrating that the form  of the CR flux  equation is not a significant source of uncertainty here. We also find an excellent agreement between the zeroth- and two-moment methods in a pure diffusion test given a high enough reduced speed of light.

It is worth noting that our CR treatment is akin to the first-moment or “M1”
moment-based method for radiation hydrodynamics (with different closure relations and the scattering terms), with the “reduced speed of light” c̃ <cit.>, while the “pure diffusion” approximation is akin to flux-limited diffusion (without the limiter).










  §.§.§ The Diffusion Coefficient


The only remaining unspecified parameter in the CR treatment is the effective diffusion coefficient κ[ We do not attempt to calculate the diffusion coefficients from microphysics, but treat them as empirical parameters to be varied/constrained.]. However, there is still substantial uncertainty on its value from a theoretical or observational perspective. In the self confinement picture, it depends on wave damping mechanisms, which are currently not well constrained <cit.>. In the extrinsic turbulence picture, CRs are scattered through turbulent magnetic fields, but we have limited knowledge of the small scale magnetic fluctuations and the coupling between magnetic field turbulence and CRs <cit.>.

Fortunately, there are some empirical constraints on the effective diffusion coefficients, i.e., the diffusion coefficients that broadly reproduce observations of cosmic rays in the Milky Way (even though it is possible that the microscopic model of diffusion is not the correct one for cosmic-ray transport).


For example, <cit.> constrained the isotropically-averaged diffusivity κ to be ∼ 6× 10^28 cm^2/s to within a factor of a few, at ∼GeV energy with GALPROP, using the measured energy spectra and abundances of nuclei species in CRs, and adopting a flat halo diffusion model (; see Introduction for a brief description). Implicitly, these abundances depend on the residence time of CRs in the Galaxy, so there is a degeneracy between κ and the CR halo height z_ h (typically 1-10 kpc), out of which CRs can freely propagate (see Figure 3 in  or Figure 10 in ; this issue was also discussed in ). Even in this model, it is possible to match the observational data with a significantly larger κ (up to factors of several) if a larger halo size is adopted.[For example, Fig. 10 in <cit.> shows isotropically-averaged κ∼ 3× 10^28 cm^2/s with z_ h∼ 3  kpc but κ∼ 10^29 cm^2/s with a larger z_ h∼ 5  kpc).]

There are other substantial uncertainties in the estimates of κ, as these empirical constraints usually neglect e.g. local variations in κ or magnetic field structure, the role of advection, halo density profiles (in addition to sizes), small-scale gas density variations (“clumpiness”), and the complicated spatial and temporal distributions of CR sources. The value of κ is even more poorly constrained outside the MW. 

Given these uncertainties, we do not attempt a self-consistent calculation of the diffusion coefficient. Instead, we simply assume a constant κ, which is a common approach in the literature <cit.>, and test a wide range of κ.

Unlike a flat halo or “leaky box”-type diffusion model, where CRs simply freely escape after crossing the boundary of the halo, we assume CR diffusion with constant κ everywhere, even at large heights above the disk. It is therefore likely that our simulations will require a larger κ than the value from a flat halo model with a small halo size.

We will also consider anisotropic CR diffusion with a constant parallel diffusivity. Because the above estimate is isotropically-averaged, if magnetic fields are tangled or toroidal, the equivalent anistropic diffusion coefficient κ would be factor ≳ 3 larger.




  §.§.§ Sources & Injection
 

We assume CR injection from SNe (including Type Ia and Type II), with a fixed fraction ϵ_ cr (=0.1, as our default value) of the initial ejecta energy (Δ E_ cr = ϵ_ cr E_ SNe with E_ SNe≈ 10^51  erg) of every SNe explosion going into CRs. SNe explosions inject thermal and kinetic energy into neighboring gas resolution elements according to the algorithm described in detail in <cit.>; we therefore reduce the coupled energy by 1-ϵ_ cr and inject the remaining ϵ_ cr energy alongside the metals, mass, and thermal+kinetic energy using the same relative “weights” to determine the CR energy assigned to each neighbor. Likewise the CR flux is updated assuming the CRs free-stream at injection (F_ cr→ F_ cr + Δ F_ cr with Δ F_ cr = Δ e_ cr c̃ r̂ from the source, where r̂ is a unit vector pointing outwards from the source). The injection is therefore operator-split and solved discretely (associated with each SNe). 






  §.§.§ Hadronic & Coulomb Losses (“Cooling”)


We adopt the estimate for combined hadronic (Λ̃_ cr,had) plus Coulomb (Λ̃_ cr, Cou) losses, Γ_ cr, from <cit.> and <cit.> as synthesized and updated in <cit.>: 

    Γ_ cr =Λ̃_ cr e_ cr n_ n = (Λ̃_ cr,had + Λ̃_ cr, Cou) e_ cr n_ n
     =     5.8× 10^-16 (1 + 0.28 x_e)( e_ cr/ erg  cm^-3)( n_ n/ cm^-3)  erg  cm^-3s^-1

where n_ n is the number density of nucleons and x_e is the number of free electrons per nucleon. Following <cit.> we assume ∼ 1/6 of the hadronic losses and all Coulomb losses are thermalized, adding a volumetric gas heating term 

    S_ gas =     0.98×10^-16 (1 + 1.7 x_e)( e_ cr/ erg cm^-3)×
       ( n_ n/ cm^-3)  erg  cm^-3s^-1

The remaining CR losses are assumed to escape in the form of γ-rays and other products to which the gas is optically thin.  

Due to the hadronic and Coulomb losses, we have to consider the Boltzmann equation with a weak collision term, instead of the Vlasov equation. Since the collision term affects both CR energy density and flux, in the two moment method, we also update the CR flux as F_ cr→ F_ cr (1- Λ̃_ cr n_ nΔ t).

The loss and heating terms are operator-split and solved together with all other gas heating/cooling terms with our usual fully-implicit cooling scheme described in <cit.>.



  §.§.§ “Isotropic” Runs


By default, we solve the CR equations coupled to the ideal MHD equations, and treat the CR transport (streaming and advection/diffusion) fully anisotropically. However in many of the tests below we consider isotropic CR diffusion without MHD and streaming, so we simply solve the hydrodynamic equations, remove the terms relevant to streaming, and replace B̂ wherever it appears above (representing projection of motion along field lines) with ∇̂ P_ cr.





§ SIMULATION SETUP




 §.§ Initial conditions


We study the impact of CRs on three characteristic types of galaxies, dwarf (Dwarf), sub-L⋆ starburst (Starburst) and L⋆ (L⋆ Galaxy) galaxies, whose details are listed in Table <ref>. 

All of the runs have exponential stellar and gas disks with scale radii d_ *, disk and d_ g, disk respectively. We also include small stellar bulges with Hernquist profiles <cit.> and gas halos with beta profiles (beta=2). The latter enable CRs to diffuse far from the galaxies, since CRs cannot diffuse without the presence of neighboring gas particles in our numerical scheme. 

Halo spin parameters (which determine the rotation of the halo gas and dark matter) are set to be 0.033, close to the median value of simulated halos in <cit.>, and the initial Toomre Q is set to one uniformly in the gas and stellar disks.  We set the metallicity of all star and gas particles in our initial conditions (ICs) in Dwarf, Starburst and L⋆ Galaxy to be 0.1, 0.3, and 1.0 Z_⊙ respectively. Ages of stars present in our ICs are set to > 10 Gyr to avoid excessive SNe from old stellar populations when the simulation begins, which could significantly affect the early evolution of our simulations.

In all L⋆ Galaxy and Dwarf runs, we delay turning on the CRs because of initial instabilities from settling of the ICs and to allow magnetic fields to first amplify to a steady-state strength. We enable CRs after initial evolution of 150 Myr in L⋆ Galaxy and 300 Myr in Dwarf. In the runs with magnetic fields we start with a seed magnetic field with 10^-2μG uniformly (over all gas particles) pointing along the direction of disk angular momentum. The magnetic fields rapidly amplify to ∼μG in dense gas and develop toroidal morphologies with significant turbulent structure, by around a hundred Myr (see ). In the following, we define t=0 at the time when CRs are turned on.

Starburst is designed to mimic dwarf galaxies with high gas surface density (∼ 0.1  g/cm^2) and SFR (∼ 5 / yr) (e.g. M82 or NGC253). We set up a massive gas reservoir with the extended disk and halo such that gas can continuously accrete to the galaxy and trigger intense star formation for an extended period of time. In Starburst runs, we inject CRs immediately at the beginning of the run, since we want to study the transient phenomena (namely, the starburst).

For a subset of our runs we have performed resolution studies and show (see Appendix <ref>) that global quantities of interest are robust at our default resolution indicated in Table <ref>, and that main qualitative effects of CRs on galaxies can be captured at this resolution. 










 §.§ Cosmic Ray Propagation Models



We consider several different CR propagation models, and a range of diffusion coefficients. All models are listed in Table <ref>. In particular, we consider diffusion coefficients up to ten times higher and lower than the common inferred isotropically-averaged MW values,  sampling a range κ = 3 × 10^27 cm^2/s - 3 × 10^29 cm^2/s.

The most complete (and potentially the most realistic) CR propagation model we test includes fully anisotropic diffusion with MHD and streaming. However, given the uncertainties in the magnetic field configuration on small scales as well as uncertainties in the streaming parameters, we evolve a range of simulations with isotropic diffusion without streaming. This model also enables straightforward comparison with other work as it is the most prevalent propagation model in the literature <cit.>.

We apply the newly developed two-moment method ( <ref>) to both streaming and diffusion with a reduced speed of light, c̃. In Appendix <ref> we test different choices for this parameter and demonstrate that physical properties, e.g. SFR or γ-ray emission, are not affected by the choice of c̃ as long as it is equal to or larger than the values listed in Table <ref>.





§ RESULTS



 §.§ Distribution of cosmic rays  and the effects on galactic properties



  §.§.§ Dwarf and L⋆ galaxies


















We begin with a brief overview of the evolution of the gas, CRs and magnetic fields in our L⋆ Galaxy simulations with different CR propagation models. Once initial transients are damped away, the galaxy has a relativity steady, low SFR with weak galactic winds driven by SNe and other stellar feedback processes. The magnetic fields amplify and develop irregular yet roughly toroidal morphology through dynamo action <cit.>. After 150 Myr when we turn on CRs, the galaxy is in approximate steady state. SNe inject a fraction of energy into the surrounding gas as CRs, which is transported via advection, diffusion, or streaming. The timescale for CR hadronic and Coulomb losses is long enough that steady-state CR pressure support can arise within/around the galaxy. The total CR energy at any time roughly follows the CR injection from SNe, which is proportional to the total stellar mass formed (see Figure <ref> and the related text). However, there are also other CR energy gain and loss processes, which we will investigate in   <ref>. But in all runs, the CR source distributions are much more concentrated than CR densities, as CRs move from their “birthplace”. 

Figure <ref> shows the distribution of the CR energy density in a 60kpc×60kpc slice centered on our simulated L⋆ Galaxy. Runs with higher diffusion coefficients result in lower CR energy density at the galaxy center but develop stronger CR pressure away from the disks. The strong CR pressure gradients continue accelerating gas out to a large radius in the radial direction, although stellar feedback without CRs can also drive winds. In Figure <ref>, we also show that galaxies with CR diffusion in general have the smoother CGM structure, and outflowing gas is present further from the galactic centers. 
The study of CR driven outflows, including a comparison with thermally-driven outflows and their effect on the CGM, will be presented in a companion paper (Chan et al., in preparation).  Simulations with only streaming (κ=0 but v_ st 0) are similar to cases with very low diffusivity (κ≲ 3×10^27  cm^2 s^-1), where CRs are largely confined to the galaxy. From Figure <ref> it is clear that combining streaming with diffusion lowers the concentration of CRs in the disk plane and spreads them to larger distances. In almost all of our runs with non-negligible diffusion (κ≳ 10^27  cm^2 s^-1), diffusion makes the CR energy distribution approximately spherical, as opposed to flattened (only the streaming-only and advection-only runs show strongly flattened e_ cr, as the CRs do not efficiently escape the star-forming disk).


These scalings are easy to qualitatively understand. In the ISM, bulk transport speed for streaming is typically, v_ st∼ 10  km s^-1(B/μ G)(n_ ISM/0.1  cm^-3)^-1/2 (see ), giving a transport time through a gas halo with a radius ℓ of t_ st∼ 100   Myr (ℓ/ kpc)(v_ st/10  km s^-1)^-1, while the corresponding diffusive transport velocity/time is v_ di∼ 330  km s^-1(κ/10^29 cm^2 s^-1)(ℓ/ kpc)^-1 and t_ di∼ 3  Myr (ℓ/ kpc)^2(κ/10^29 cm^2 s^-1)^-1 . Thus even for quite low κ, the diffusive flux dominates transport on sub-kpc scales.

But because the diffusion time scales as ∼ℓ^2, if the CRs establish a smooth profile with scale length ≳ 1-10kpc, then on the larger scales the diffusion time eventually could become larger than the streaming transport time, i.e. outside a scale ℓ∼ 30  kpc (κ/10^29 cm^2 s^-1)(v_ st/ 10 km s^-1)^-1. 


We quantify the above observations with Figure <ref>, which shows the cumulative distribution of CR energy in the Dwarf and L⋆ Galaxy runs. CR energy density is most extended vertically in simulations with the largest diffusion coefficients and it is most concentrated in Advection simulations. We define the (3D) CR scale radius r_ cr, 1/2 such that the sphere with r_ cr, 1/2 encloses one half of the total CR energy. In the L⋆ Galaxy, we find r_ cr, 1/2 is around 3 kpc in run “κ=3e27”, but it increases to around 10 kpc in run “κ=3e28” and 30 kpc in run “κ=3e29.” Similarly, the scale-height of the CR energy distribution also increases with increasing κ. Trends of the CR scale radius with κ can be understood with a diffusion model where the CR injection time (E_ cr/Ė_ inj) is comparable to the CR diffusion escape time (t_ di∼ r_ cr,1/2^2/κ) where Ė_ inj is the CR energy injection rate. Assuming a similar injection time, we find that r_ cr,1/2 is roughly proportional to √(κ), so a faster diffusion leads to a more extended CR distribution.


We show the CR and gas mid-plane densities for our L⋆ Galaxy and Dwarf runs in Figure <ref>. For both galaxy types, CR density profiles are significantly “flatter” (more extended and less centrally-concentrated) with higher κ. Consequently, in runs with fast diffusion, CRs have smaller impact on the central region of galaxies, providing less pressure support to the central gas, but they can be more important in the CGM. Interestingly, the “Advection” runs have lower CR central densities than “κ=3e27” because of the smaller adiabatic energy gain (a highly non-linear effect), which we will discuss in the next section and Figure <ref>.
The gas midplane density depends rather weakly on CRs. Low-κ (or streaming/advection-only) runs have slightly higher midplane densities while higher-κ runs have midplane densities similar to “no CR” runs. 
  
This is likely caused by additional pressure support from CRs trapped in the midplane in low-κ runs, which allows gas to reach higher densities before fragmenting and forming stars.

It is interesting to compare the CR energy density in our L⋆ galaxy model with that observed near the solar circle (e_ cr∼ 1  eV/cm^3, see e.g. ), but we must recall that the L⋆ model was not constructed to be an exact MW analogue. For example, it has a more steeply-rising central gas density, without the gas deficiency that appears in the center of the MW (i.e. it does not have a “star-forming ring” and corresponding “hole” in the central few kpc), and the gas densities at ∼ 8kpc from its center are lower than the ∼ 1  cm^-3 in the solar neighborhood <cit.>.

Nevertheless, the model has a stellar mass, gas mass, and SFR similar to the MW. Our runs with isotropic κ≤ 3× 10^27  cm^2 s^-1 produce a mid-plane e_ cr at ∼ 8kpc from the galaxy center which is high relative to the observed value, those with isotropic  κ≥ 3× 10^29  cm^2 s^-1 are lower, while those in-between are reasonably consistent. Turning on MHD (making the diffusion anisotropic) increases e_ cr by factors ∼ 2-3, consistent with the isotropically-averaged κ being lower by a similar factor (as expected), so values κ∼ 1-3×10^29  cm^2 s^-1 are marginally more favored. Given the lack of a detailed match between our models and the MW, stronger constraints on CR propagation come from γ-ray emission in  <ref>. 

Fig. <ref> shows cumulative SF histories: akin to the disk midplane-pressure effects above, CR runs with very low κ suppress SF by modest factors ∼ 1.5-2, an effect which vanishes at higher κ. Smaller variations (∼ 10%-level) are generally dominated by stochastic run-to-run variations. Runs with MHD generally show slightly higher SFRs (all else equal), an effect discussed in detail in <cit.>, but the stochasticity during the early evolution can wash out such effect.









  §.§.§ Starburst galaxy

The Starburst model is designed to reach high SFRs and gas densities, which are transient phenomena since strong stellar feedback after the starburst will disrupt the galaxy and reduce the gas density. Thus, our Starburst run reaches SFR peaks of ∼ 5   yr^-1 with highest central gas densities ∼ 100  cm^-3 or (edge-on) surface densities ∼ 0.1  g cm^-2 (∼ 500  pc^-2) (compare Fig. <ref> and ), which last ∼ 10Myr. Between starbursts the galaxy has lower SFR and gas densities, with correspondingly longer hadronic loss times, so CRs can escape more easily.


Fig. <ref> shows that CR energy densities during starburst phases are around 100 eV/cm^3 with slow transport, similar to the value inferred from modeling the observed γ-ray spectra of e.g. M82 ( and ). Although these are high relative to the MW, and a factor of several higher than the thermal (or magnetic) energy densities, they are lower than the pressure required for hydrostatic balance (π G Σ_ gΣ_ g∼ 10^3-4 eV/cm^3), which is primarily comparable to the kinetic energy density in these galaxies (with turbulent velocity dispersions similar to those observed, ∼ 50-100  km s^-1). Our findings are therefore consistent with earlier claims by <cit.>, , and others who showed that CRs are dynamically unimportant at least in the cores of the starbursts, but they might be more important away from the central dense region.




 §.§ Pionic γ-ray emission as a measure of CR propagation












Owing to the lack of direct measurements of primary and secondary CRs at low (∼ GeV) energies from extra galactic sources, pionic γ-ray emission is one of the few observables that constrain CR propagation outside of the MW. CRs interact with nuclei and produce pions that decay into pionic γ-rays. While there is a substantial amount of pionic γ-ray emission with energy < 1 GeV, it is difficult to isolate it observationally owing to contamination by leptonic emission. For γ-rays with energies > 1 GeV, the leptonic emission is less than one tenth of the pionic emission (for CRs with a spectrum consistent with our default model assumptions; see calculations by ). We will ignore additional potential channels of  > 1 GeV γ-ray production, e.g. pulsars or dark matter annihilation. Hence, in the following, we assume all > 1 GeV γ-rays are pionic (if there is substantial pulsar contamination, the pionic γ-ray emission is lower, and higher diffusivities κ are required).

The > 1 GeV γ-ray luminosity for γ-rays L_γ(≥ GeV) can be calculated as:

    L_γ(≥ GeV)≈∑_i 1/3β_πΛ̃_ cr,had e_ cr n_ nΔ V_i,

where we sum over gas particle i with volume Δ V_i. First, the most of the hadronic loss (Λ̃_ cr,had e_ cr n_ n in Eq. <ref>) is responsible for the pion production. Second, only one third of the pions (π^0) produce γ-rays. Third, β_π(≈ 0.7) is the fraction of the pionic γ-rays with energy above GeV<cit.>, which is calculated with the GALPROP pionic cross sections <cit.> built on <cit.>, assuming the CR spectrum (between 1 GeV and 1 PeV) follows E^-p, where E is the CR proton energy and p (= 2.2) is the spectral index. 

If CRs can propagate fast enough that a significant fraction can leave galaxies without interacting with ISM, the γ-ray emission will be relatively weak, compared to the expectations from the CR injection. We follow <cit.> and  and quantify this by comparing the pionic γ-ray luminosity L_γ (above >1 GeV) with the bolometric “star formation luminosity” L_ SF (UV/optical/IR luminosity ultimately contributed by stellar radiation from massive stars, estimated assuming a time-constant SFR and the same stellar IMF as in our simulations), since the CR injection is proportional to SN injection rate and thus to the SFR. If the SFR is constant and we are in the “proton calorimetric limit” (all CR protons instantly lose their energy to collisional processes, without any other processes influencing their energies or spatial distribution, assuming the same time-constant SFR), then the ratio L_γ/L_ SF is approximately constant. 

The value of L_γ in the calorimetric limit is derived as follows. If the SFR is constant, the SNe rate is dominated by Type-II events, and the CR injection rate is: 

    Ė_ cr,SN/ [erg/s] =ϵ_ crE_ SN/ [erg] []/<m_*>ξ(m_*>8) SFR/ 3.2× 10^7[/yr],

where ϵ_ cr(=0.1) is the fraction of SNe energy going into CRs (constant by assumption in our simulations), E_ SN(=10^51  erg) is the energy from one supernova (also constant by assumption), ξ(m_*>8)(=0.0037) is fraction of stars that end as supernovae, and <m_*>(=0.4 ) is the mean stellar mass, both calculated for the same <cit.> IMF used in our simulations. If this injection rate of CRs is balanced by collisional losses without any other energy gain/loss processes, i.e. Ė_ cr,SN=e_ crΛ̃_ cr, then the pionic γ-ray luminosity is 

    L_γ(≥ GeV)^ calor≈ 6.7× 10^39 SFR/ [/yr]  erg/s .



The corresponding “star formation luminosity” L_ SF for a constant SFR assuming again the same <cit.> IMF adopted in our simulations is

    L_ SF^ constant≈ 3.5× 10^43 SFR/ [/yr] erg/s,

where the prefactor is calculated with STARBURST99[<http://www.stsci.edu/science/starburst99/docs/default.htm>] <cit.>.[ adopted a different conversion factor because they assumed a Salpeter IMF <cit.> (following ). But the ratio L_γ/L_ sf in  is only higher by ∼ 10%  since the SNe rate is also adjusted accordingly.] So in the constant-SFR, calorimetric limit, we would expect L_γ/L_ SF≈ 2× 10^-4 <cit.>.








  §.§.§ Pionic γ-ray emission



Figure <ref> shows the projected pionic γ-ray surface brightness of the L⋆ Galaxy run, for different values of κ. γ-rays mostly originate from galactic disk, i.e. where gas and CR densities are the highest. The γ-ray surface brightness drops by over an order of magnitude a few kpc away from the disk plane. For higher-κ, the emission is dimmer but more spatially extended (reflecting the CR energy distribution). 

Fig. <ref> quantifies the distribution of the γ-ray luminosity within spheres of increasing radii for all of our runs with cosmic rays. Consistent with the discussion above, γ-ray emission is much weaker for large κ.[As with the CR energy density, we note that sometimes the runs with finite-but-low κ exhibit slightly higher L_γ even than the “Advection” runs owing to non-linear effects discussed in  <ref>.] For our L⋆ Galaxy, half of the γ-ray luminosity originates from the inner 5-7 kpc. The half-luminosity radius is smaller for our dwarf galaxy, as expected, since the galaxy itself (gas and stellar) is smaller. 


Fig. <ref> breaks down the distribution of CR energy as a function of local gas density, which (since L_γ∝ e_ cr n_ gas) effectively determines L_γ. With low-κ (or with advection/streaming only), CRs reside longer in the high-density regions where they are injected. If the density (on some scale ℓ of the cloud or disk) is larger than ∼ 10  cm^-3 (ℓ/ kpc)^-2 (κ/10^29  cm^2 s^-1), then the collisional loss time of CRs becomes shorter than the diffusion time, so the CRs decay close to their injection. This means L_γ is lower at higher κ, even when the SFR (injection rate) is higher, because the bulk of the CR energy is at lower n_ n<0.1  cm^-3.









  §.§.§ L_γ/L_ SF ratio and scalings





Fig. <ref> shows the time evolution of L_γ/L_ SF. With a lower κ (“Advection” or “κ=3e27”), the galaxies are closer to the calorimetric limit (as expected).[Note there are periods where L_γ/L_ SF exceeds calorimetric – this is perfectly allowed. Usually it occurs because of short-timescale drops in the SFR and hence L_ SF, where the CRs take somewhat longer to decay so L_γ is ∼ constant. But it can also owe to adiabatic compression increasing CR energy, or the contribution of SNe Ia's, or smoothly declining SF histories, all of which violate the assumptions used to derive the calorimetric limit.] But the “κ=3e29” runs are lower than the calorimetric limit by more than an order of magnitude in our Dwarf and L⋆ runs.
 

Turning around our argument from  <ref> above, if CRs are injected in a structure of size ℓ and gas density n with an isotropically-averaged diffusivity κ≲ 10^28  cm^2 s^-1 (ℓ/ kpc)^2 (n/ cm^-3), the collisional loss time becomes shorter than the escape time, so we expect near-calorimetric L_γ/L_ SF. On the other hand, at larger κ≫ 10^28  cm^2 s^-1 in the limit where CRs do escape, if we assume the galaxy gas structure is otherwise similar, then the ratio of escape time to loss time (proportional to the fraction of CR energy lost in escaping, and therefore L_γ/L_ SF) should scale ∝κ^-1.



In  <ref>, we will show that adiabatic processes are of secondary importance relative to diffusion in reducing L_γ in runs with high-κ, but in runs with low-κ (e.g. “κ=3e27”), they boost L_γ considerably.


As expected, anisotropic diffusion tends to suppress the isotropically-averaged κ by factors ∼ 1.5-3, and correspondingly increase L_γ/L_ SF. Streaming slightly increases the escape and dissipates CR energy through streaming loss, so decreases L_γ/L_ SF, but the effect is very small (and streaming alone produces near-calorimetric results). This is  because (as discussed in  <ref> above) the streaming escape time is much longer than the diffusive escape time, even for relatively low κ, but with the caveat that we do not consider the decoupling between CRs and gas in the cold ISM due to the low ionization fraction and ion-neutral damping <cit.>, which could significantly reduce γ ray emission from dense gas.
 

For the same CR propagation model, galaxies with higher gas densities and larger sizes (effectively larger column densities of dense gas with which CRs must interact to escape) have a larger L_γ/L_ SF , which can be seen in Fig. <ref>. 




  §.§.§ Comparison to observations


We now compare the simulations to observational estimates of L_γ/L_ sf as a function of either central gas surface density in galaxies (Σ_ g) or SFR, as compiled in . Most of the observed data is described in , but we also include the SFR of the SMC (0.036 /yr from ). We add two extra starburst galaxies in Figs. <ref>-<ref> (NGC 1068 and NGC 4945; which are listed in Table 2 of  but not in their figures). The SFRs of starburst galaxies (NGC1068, NGC253, NGC4945, M82) are obtained with the <cit.> IR to SFR conversion formula, assuming the Kroupa IMF (Their IR luminosities are also listed in Table 2 of ). 

γ-ray observations of nearby galaxies are limited in spatial extent due to energy resolution and contamination from the diffuse backgrounds and foregrounds <cit.>. Hence, we only consider γ-ray emission within 3 kpc for Dwarf, matched approximately to that used for the SMC. This choice reduces the γ-ray luminosity by a factor of two compared to using an infinitely large aperture. For Starburst and L⋆ Galaxy, we take 10 kpc apertures (which only reduces L_γ by tens of percent compared to an infinitely-large aperture), matched to those used for e.g. M31, NGC1068, and M82. See Fig. <ref> for how this scales with size.
 

We measure the gas surface densities Σ_ g (averaged over viewing angles) of Starburst within 250 pc, L⋆ Galaxy with 4 kpc, and Dwarf within 2 kpc – chosen to be twice the sizes of the active star-forming region in Dwarf and L⋆[This choice is smaller than the optical radii that <cit.> and  used, but the gas surface densities are similar in both choices.] and about equal in Starburst (similar to the choice in  and ). 




Figs. <ref>-<ref> compare our simulations with the observations (compare to Figure 2 of ).  As expected based on our discussion above, L_γ/L_ SF is high and close to the calorimetric limit for Dwarf  and L⋆ Galaxy with slow CR transport, i.e. for “Advection”, “κ=3e27”, and “MHD Streaming”. These values are clearly well above the observationally inferred L_γ/L_ SF.

With larger diffusion coefficients, L_γ/L_ SF decreases as expected. For isotropic diffusion, the observations in dwarf and L⋆ galaxies appear to require κ_ isotropic∼ 10^29  cm^2 s^-1. Out of the options tested κ = 3× 10^29  cm^2 s^-1 provides the best match but the range of data allows for a slightly lower value as well. 
 
For anisotropic diffusion, L_γ is somewhat larger owing to suppressed isotropic-averaged diffusivity, as discussed above, so values of the parallel diffusivity κ_ > 10^29  cm^2 s^-1 are favored. 

For galaxies with high gas surface densities and SFRs, i.e. Starburst, we found κ has to be less than 3× 10^29 cm^2/s. On the face value, this implies that CR transport is effectively slower in high gas surface density regions or during starburst. However, in our Starburst runs, for the highest diffusion coefficient tested, we did not include MHD and anisotropic diffusion that, depending on the magnetic field configuration, could slow down the transport of CRs out of high density regions. The models also have gas configurations and typical gas densities that are not an excellent match to observed starburst galaxies, and our analysis only includes brief time intervals when central gas density reaches the values similar to observed starbursts.

Many observations, e.g. <cit.> and <cit.>, considered L_γ with 0.1  GeV<E_γ<100 GeV, instead of E_γ>1 GeV. Therefore, in Appendix <ref>, we compare L_ 0.1-100 GeV from simulations with observations, and we find that the same high diffusion simulations provide the best match to the observed γ ray emission from galaxies.

It is interesting to compare the results from the simple leaky-box model of , as well as more detailed models of CR transport in the MW, with our findings.  predicts L_γ/L_ SF as a function of gas surface density by assuming the Kennicutt-Schmidt law and a one-zone leaky box model with CR diffusion (see  for details), with an isotropically-averaged κ = 3×10^28  cm^2 s^-1. The broad contours of their prediction for L_γ/L_ SF as a function of Σ_ g or SFR are similar to our simulations, suggesting – as they argued – that CR escape is required to reproduce the observed trend of L_γ/L_ SF. 

In the MW, much more detailed propagation models have been tested (see e.g.  and reference in  <ref>). We again caution that our “L⋆” model is not an exact MW analogue, since it has higher gas surface density and lacks a central gas deficiency like the MW <cit.>, both of which could affect L_γ. 

Note that at “face value”, both MW and  constraints might appear to favor slightly-lower κ∼ 3-6×10^28  cm^2 s^-1 compared to the best-match here (κ∼ 3×10^29  cm^2 s^-1), but this is a relatively small offset and completely expected if we account for the points below. (1) The MW observations and  models assume relatively small halos out of which the CRs escape instantly, while we assume a constant κ everywhere, meaning that our effective halo size is large (∼ 10-30kpc). Recall ( <ref> and  <ref>), the inferred κ in the observations increases with the halo height. (2) The gas in the simulations is clumpy where CRs are injected, slightly increasing L_γ  <cit.>, compared to the smooth mass profiles assumed in those studies (requiring larger κ by a factor ∼ 1.5-2). (3)  did not consider galactic winds and adiabatic losses/gains in their fiducial models; the MW constraints did not account for galactic winds in a self consistent manner (i.e. they do not consider CR-driven winds and the radial/temporal variations of the winds). In our anisotropic runs, we also find that the isotropically-averaged κ (what is nominally constrained by the  study, for example) is a factor ∼ 2-3 lower than the parallel κ. Accounting for all of these facts, our favored coefficients appear to be consistent with other state-of-the-art constraints on CR propagation in the MW from e.g. , and references in  <ref>.




§ DISCUSSION




 §.§ Comparisons to previous studies




  §.§.§ Suppression of star formation by cosmic rays


In our idealized non-cosmological simulations, we find that SF can be suppressed by CR feedback in simulations with either advection or streaming only, or very low diffusivities κ≲ 10^28  cm^2 s^-1, consistent with many previous findings, e.g. <cit.>. However, such slow transport severely violates constraints from observed γ-ray emission, and at best results in modest SFR suppression (factor ∼ 1.5-2). For larger transport speeds required to reproduce the observed γ-ray emission, CRs have only a weak effect on SF. 

Interestingly, <cit.> found that while CRs reduce SFRs in dwarf galaxies, they have almost no effect in MW mass galaxies. Their conclusion was likely due to their “local equilibrium” assumption, namely that CR injection (∝ SFR∝ρ^1.5 in their model) is balanced by collisional losses (∝ρ) locally (like in the calorimetric limit), in an isothermal-like ISM, so in their models the CR energy density is proportional to ρ^1/2 while thermal energy densities are proportional to ρ: as a result, CR energy was always sub-dominant to thermal energy in their models at gas densities n>0.2  cm^-3. In contrast, in our simulations, CRs can propagate far from their injection sites, 
so local equilibrium is not valid and we find that the ratio of CR pressure to gas thermal or turbulent pressure for low κ can be significant even at moderate ISM densities, providing mild suppression of the star formation
(similar arguments were presented in 
). For our favorite, large κ values CRs escape from the ISM, resulting in practically no effect on the star formation in both our Dwarf and L⋆ Galaxy simulations.

However, we caution that because our simulations are non-cosmological, they do not account for the effect of CRs on the CGM and IGM (the source of fuel for galaxies). As CRs escape the galaxies more efficiently with the favored larger κ, we have shown they have proportionally much higher energy density/pressure in the CGM, which means they could (in principle) be important for the long-term cosmological evolution and accretion onto galaxies. This is likely most important in more massive galaxies that build quasi-hydrostatic halos whose late-time cooling influences galaxy growth. We will explore this in cosmological simulations in future work<cit.>.





  §.§.§ γ-ray emission



Our results are in line with : when matching the observed γ-ray emission, starburst galaxies (with effective isotropic diffusivities κ < 3× 10^29 cm^2/s) are nearly proton calorimeters, while galaxies with lower gas surface density or SFRs (with κ∼ 3×10^28-29 cm^2/s) are not proton calorimeters (most CR protons escape).

<cit.> also studied hadronic γ-ray emission with simulations of MW-mass galaxies and argued for isotropically-averaged coefficients κ∼ 3× 10^28  cm^2 s^-1; but they only considered the γ-ray emission in the CGM and they did not include hadronic/collisional losses in the simulations, which led to some unphysical results. For example, their predicted pionic γ-ray luminosity significantly exceeded the CR injection rate at lower κ. Moreover, as noted by <cit.>, neglecting collisional CR losses allows CRs to build up in dense gas or the disk midplane without being rapidly lost (as they should), which artificially enhances the strength of CR-driven winds. Nevertheless, we broadly agree on the preference for a relatively high κ.

Recently, <cit.> also investigated γ-ray emission with idealized galaxy simulations, assuming CR transport via either advection-only or advection+anisotropic diffusion with κ=10^28  cm^2 s^-1. They argued they could (a) reproduce the FIR-γ-ray correlation and (b) explain the low L_γ in non-starburst galaxies primarily by adiabatic losses. 

But there are several caveats: 

(1) Their favored model still over-predicted L_γ/L_ SF by a factor of a few or more in non starburst galaxies, e.g. dwarf and MW-mass galaxies. For their actual simulated points (see their Fig. 3) without diffusion, the predicted L_γ/L_ SF is larger than the SMC, LMC, MW, and M33 (not shown therein, but see ).

They claimed to match the observed FIR-L_γ correlation, only if an empirical FIR-SFR conversion relation <cit.> is assumed. However, as they acknowledged, this conversion relation over-predicts L_ FIR in dwarfs, due to much lower dust opacity/absorption/reddening. Their FIR-L_γ relation might deviate from observations after taking this correction into account.[Because both the L_γ/L_ SF and L_ FIR/ SFR ratios drop in dwarfs, a roughly-linear L_ FIR-L_γ relation can still maintain in our simulations (directly related to the “conspiracy” which maintains the FIR-radio correlation; for discussion see e.g. ).]

(2) We do not consider the same CR models and the same range of γ ray energy. They consider L_γ(0.1-100  GeV), i.e. including all CRs from 0.1-100GeV, instead of the choice here and in , which is restricted to L_γ(>1  GeV). They also assumed a shallower CR spectrum (∝ E^2.05, as compared to ∝ E^2.2 here and in ). Together with this, our L_γ can differ from theirs by a factor of ∼ 2-3. However, even if these differences are considered, their L_γ/L_ SF are still greater than the observed dwarfs.

They suggested their over-prediction of L_γ/L_ SF might be reconciled with simulations that could resolve the multi-phase ISM, since CRs may preferentially spend time in low density regions, which dominate the volume. We do have the multi-phase ISM here, but predict similar results in our advection-only or low-κ runs. A possible explanation for the discrepancy is that the observed low γ-ray luminosities require high diffusion coefficients κ∼ 3×10^28-29  cm^2 s^-1 as favored by our study here and the modern MW constraints <cit.>.


(3) We will show immediately below that when κ is in the favored range, adiabatic processes are less important than CR transport in reducing L_γ, although if κ is small, adiabatic processes tend to increase L_γ.










 §.§ CR Energetics and the Importance of Different Gain/Loss Terms




Fig. <ref> shows the relative importances of various CR gain/loss terms in our simulations: SNe injection, collisional (hadronic+Coulomb) losses, “streaming losses” (energy loss to excitations of Alfven waves), and “adiabatic” terms (“PdV” work lost pushing gas, or CR energy gain in compression). 

The initial injection from SNe is proportional to the SFR (with a few Myr delay), so it tracks the SFR and varies only by a relatively small amount in our different runs of a given galaxy model (even the highest/lowest SFR runs differ by at most a factor ∼ 2). 

Collisional losses are important loss terms (within the galaxies) – and we have already discussed these extensively as they are the origin of the γ-ray emission. Since they scale ∝ e_ cr n_ n they decrease with “faster” CR transport (higher κ) as CRs reach lower-density gas faster.[The (weak) exception to this rule is the  κ=3×10^27  cm^2 s^-1 run in our L⋆ and Starburst models, where collisional losses are slightly larger than in the corresponding “Advection” runs. This is caused by the slightly stronger adiabatic compression term boosting e_ cr in dense gas.] Streaming losses are comparatively small.







The “adiabatic” term Ė_ Ad can be comparable to injection or collisional terms, but can be a gain or loss process. To better understand where the gains and losses occur, Fig. <ref> shows the contribution to the total Ė_ Ad from gas with different densities n or at different galacto-centric radii n. For CRs at low ambient n (or large r), Ė_ Ad tends to be a loss term (i.e. CRs are expanding or losing energy in rarefactions). For CRs in high ambient n_ ISM and small r, it tends to be a gain (CRs are being compressed in converging flows). Recall, the “adiabatic” term is defined (Eq. <ref>) by -∫ d^3 x P_ cr (∇· v) = -(γ_ cr-1) ∫ d E_ cr ∇· v (where v is the gas velocity). So, combined with Fig. <ref>, this simply means that at high gas densities within galaxies, more of the ISM is collapsing/converging or being compressed (in e.g. shocks), while at low densities outside galaxies, more the gas is expanding in outflows. Whether one or the other term dominates depends on where most of the CR energy resides (shown in Fig. <ref>). 

So, unsurprisingly in Fig. <ref>, our runs with the most efficient CR escape to large-r and low-n (all the highest κ∼ 3×10^29  cm^2 s^-1 runs, or most of the Dwarf runs with even intermediate κ) show net Ė_ Ad < 0, since CRs rapidly migrate to the expanding regions. In contrast, those with the least efficient escape (e.g. all the Starburst runs and the L⋆ runs with lower isotropically-averaged κ) show net Ė_ Ad > 0.[Interestingly, if CRs do not preferentially stay in regions where the adiabatic term is mostly negative or positive, then the Ė_ Ad term will be relatively small, since adiabatic gains compensate adiabatic losses, as occurs in a couple of our “Advection” runs.]

For either advection-only or low-κ (∼ 10^28  cm^2 s^-1), the qualitative behaviors of Ė_ Ad in both dwarf and MW-mass systems in <cit.> are similar to what we find here. 

However, for the reasons discussed in  <ref>, our results do not support their conclusion that adiabatic losses are the dominant factor for the low L_γ/L_ SF in dwarfs. At very low κ, Ė_ Ad is primarily a gain term. But even at higher κ where Ė_ Ad<0 is a loss term, it is insufficient (in itself) to explain the very low L_γ/L_ SF observed in the SMC/LMC without significant CR leakage (the adiabatic+SNe terms are not enough to explain the loss terms in the top right panel of Fig. <ref>). More importantly, the adiabatic losses arise only after the CRs have already escaped the dense gas, i.e. the regions which produce
most of L_γ.

To explicitly show that CR leakage is significant, in Fig. <ref> we plot the fraction of CR energy that escapes the central galactic region, E_ esc/E_ source. We define E_ source as the total cumulative CR energy input within the central 6 kpc (10 kpc) and E_ esc as the CR energy that leaves this central region in Dwarf (L⋆ galaxy)[In practice, we calculate E_ source by summing up all positive CR energy gains within the central region, including SNe, adiabatic, and the small numerical error terms mentioned in footnote <ref>. To calculate E_ esc, we sum up all CR energies outside the central region and compensate for the collisional, streaming, and adiabatic losses.]. Outside of this central region the (hadronic) dissipation time is much longer than 50 Myr. At high CR propagation speeds (e.g. high κ), most of the CR energy indeed escapes the central region, where most γ rays are produced. This shows that CR escape is the main reason for  reduced γ ray emission in low-gas-density galaxies.





§ CONCLUSIONS



We explore the effects of CRs on galaxies, in high-resolution, idealized (non-cosmological) (magneto-)hydrodynamic simulations of dwarf, L⋆, and sub-L⋆ starburst galaxies, using the FIRE-2 treatment of the multi-phase ISM, star formation, and stellar feedback, accounting for CR injection from SNe, collisional (hadronic+Coulomb) losses, and CR transport via diffusion and streaming. We focus on constraining CR propagation models (e.g. diffusion and streaming coefficients) using observations of GeV γ-rays from galaxies. Our main conclusions include:




  * We adopt a newly developed two moment method for CR transport, and show that it is computationally efficient and accurate, allowing us to simulate CR transport simultaneously including diffusion and streaming with diffusivities up to ∼ 3×10^29  cm^2 s^-1 and ∼ pc resolution.



  * The CR “transport parameters”, in particular, the effective diffusivity κ≡ | F_ cr|/|∇_e_ cr| (which can, microphysically, arise from a combination of streaming and diffusion), have a significant impact on galaxy properties and predicted γ-ray emission. With very slow propagation (κ≲ 10^28 cm^2/s), CRs are trapped in the disk and contribute to the mid-plane pressure gradients, so suppress SF (albeit only by modest factors ∼ 1.5-2, if hadronic losses which limit the CR energy density are accounted for). However, these models are ruled out because they produce much larger γ-ray luminosities than the observed for dwarf or MW-like systems. At higher κ≳ 10^28  cm^2 s^-1, CRs form extended halos. This means they have weak effects on gas within the disk, but could help accelerate galactic winds or provide support via pressure gradients in the CGM.


  * The extent of the CR halo, and correspondingly the extent of the pionic γ-ray emission, increase with κ as expected. For e.g. our L⋆ galaxy, half the CR energy is located within 10 kpc (30 kpc) for κ=3×10^28  cm^2 s^-1 (κ=3×10^29  cm^2 s^-1). Correspondingly only about ∼ 50% of the γ-rays are emitted from the central few kpc.




  * In our sub-L⋆ starburst galaxies, the CR energy density reaches ∼ 10^2-3  eV cm^-3 throughout the burst and is larger than thermal or magnetic pressure in the ISM (for any κ), but is still much smaller than the energy density in turbulent motions or that required to maintain hydrostatic equilibrium. This leads to weak CR effects at the central region of starburst, consistent with the results in .




  * We constrain the average CR propagation speed/diffusivity with γ-ray (> GeV) emission from galaxies. The observed L_γ- SFR relation requires isotropically-averaged diffusivities κ∼ 3× 10^29 cm^2/s in dwarf and L⋆ galaxies, and κ≲ 3×10^29 cm^2/s in sub-L⋆ starburst galaxies.

 If CRs are transported only by gas advection, or streaming only at speeds which cannot exceed modestly super-Alfvenic values, or (equivalently) low isotropically-averaged effective diffusivities κ < 3× 10^28  cm^2 s^-1,then CRs escape galaxies too slowly and produce γ-ray luminosity close to the calorimetric limit. This over-predicts the observed γ-ray luminosities in dwarfs (e.g. the SMC, LMC, M33) and L⋆ systems (M31, the MW) by an order of magnitude or more. However, for faster transport parameters (effective κ∼ 3× 10^28-29  cm^2 s^-1), CRs escape the dense regions rapidly and the γ-ray luminosity (which scales ∝ e_ cr n_ gas) is reduced (especially in dwarf galaxies), predicting γ-ray luminosities in good agreement with those observed as a function of either gas surface density or SFR (see Figs. <ref> and <ref>).  




  * Given the transport parameters required to reproduce the observed L_γ/L_ SF, we find most CR protons escape from dwarf galaxies, i.e. low-gas-surface-density systems are not proton calorimeters, while our (sub-L⋆) starburst models are (approximate) proton calorimeters.


  * CR streaming at trans-Alfvenic speeds is relatively slow and cannot alone reduce L_γ/L_ SF significantly below the calorimetric limit in our models (as required by observations), even if we allow modestly super-Alfvenic streaming (with ∼ 4 v_ A; see Appendix <ref> ). For our favored effective κ, the equivalent streaming speed (using the fact that the CR flux is similar for a diffusivity κ or streaming speed v_ st∼κ |∇ P_ cr|/P_ cr) is ∼ 10-100 times the Alfven speed.


  * “Adiabatic” effects on CR energy densities (losses in expansion, or gains in compression) can be comparable to injection or collisional loss terms, but cannot alone reduce L_γ/L_ SF close to the level required by observations of the MW/SMC/LMC/M33. In dense gas within the galaxies, the net effect of these terms is primarily to increase CR energy density (and L_γ), while in low-density gas outside galaxies, it is primarily to decrease the CR energy via expansion in outflows. This means that CR “adiabatic losses” are significant only after CRs already diffuse out of the dense ISM gas (where γ-rays are produced). 



Our study only scratches the surface of the rich phenomena of CRs in galaxies and leaves out many important details. 
For example, it is clearly important to study the effects of CRs in cosmological galaxy simulations, which can treat CRs and magnetic field evolution consistently, explore the effects of CRs on magnetic field amplification, self-consistently generate starburst systems (in e.g. mergers), and (perhaps most importantly) explore the interaction of CRs with inflows and outflows in a “live” CGM/IGM environment. We will explore such cosmological runs in future work (Hopkins et al., in preparation).

Although we briefly mentioned the effects of CRs on galactic winds (which are ubiquitous in these simulations), we have not investigated them here. It has been proposed that CR-driven winds could have very different phase structure (compared to thermally-driven winds) and strongly modify the CGM properties <cit.>. Although extensive literature on this topic exists, detailed study of CR winds in (cosmological) simulations that can already self-consistently drive galactic winds with stellar feedback <cit.> is largely unexplored. Our simulations provide a unique combination of a high-resolution, multi-phase ISM, with explicit treatment of local star formation in self-gravitating substructures, individually time-resolved SNe and their thermal and momentum feedback combined with CR injection and transport. It will therefore be especially interesting to explore the effects of CRs on the development of galactic winds (Chan et al., in preparation).

We also do not study another important indirect CR constraint, the radio emission from synchrotron radiation, which has been observed in many galaxies <cit.>. These observations provide independent constraints on primary CRs, secondary CR electrons from CR protons, and magnetic fields. It is worth also exploring the observed FIR-radio correlation <cit.> with galaxy simulations in a manner similar to our analysis of the connection between SFR and γ-ray emission. However, as mentioned in <cit.>, these correlation requires the consideration of secondary CRs, which we plan to incorporate in the future. 

Because this was an idealized parameter study, we have adapted a simple model with a constant isotropic/parallel diffusivity κ. But in essentially any physical model, this coefficient depends on local properties of the gas and CRs, in a manner which remains deeply uncertain both theoretically and observationally (see e.g. ). It would be interesting to investigate galaxy evolution and CR observables in studies where the CR transport coefficients vary dynamically and locally <cit.>, or with recently-developed models which attempt to actually predict the coefficients self-consistently <cit.>.
 




  


§ ACKNOWLEDGEMENTS

We thank the anonymous referee for the detailed comments that helped to improve this manuscript. We thank Patrick Diamond, Ellen Zweibel, Michael Norman, and Todd Thompson for insightful suggestions and advice, Eve Ostriker for pointing out a typo, and Bili Dong for his help with YT. We would like to thank the Simons Foundation and the participants of the Galactic Superwinds symposia for stimulating discussions. TKC and DK were supported by NSF grant AST-1715101 and the Cottrell Scholar Award from the Research Corporation for Science Advancement. Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NSF Collaborative Research Grant #1715847 and CAREER grant #1455342, and NASA grants NNX15AT06G, JPL 1589742, 17-ATP17-0214. EQ was supported in part by a Simons Investigator Award from the Simons Foundation and by NSF grant AST-1715070. The Flatiron Institute
is supported by the Simons Foundation. CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522, by NASA through grants NNX15AB22G and 17-ATP17-0067, and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. 

The simulation presented here used computational resources granted by the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant no. OCI-1053575, specifically allocation TG-AST120025. Numerical calculations were also run on the Caltech compute cluster “Wheeler”, allocations from XSEDE TG-AST130039 and PRAC NSF.1713353 supported by the NSF, and NASA HEC SMD-16-7592. This work uses data hosted by the Flatiron Institute's FIRE data hub. This work also made use of YT <cit.>, matplotlib <cit.>, numpy <cit.>, scipy <cit.>, and NASA’s Astrophysics Data System.

mn2e







§ A COMPARISON OF DIFFERENT STREAMING PARAMETERS






In  <ref>, we discussed uncertainties in the parameters describing CR streaming: both the streaming speed v_ st and magnitude of the streaming loss term Γ_ st=- v_ A·∇ P_ cr owing to excited and thermalized high-frequency Alfven waves (independent of streaming speed; see <ref>). Here we explore these more thoroughly. We consider four model variations:


  * “v_ st=v_ A”: This is our default choice in the main text, with streaming speed equal to the Alfven speed.


  * “v_ st=3v_ A”: Here v_ st = 3 v_ A, one of the super Alfvenic speeds considered in <cit.>.


  * “v_ st=4v_ A”: Here v_ st = 4 v_ A, another super Alfvenic speed considered in <cit.>


  * “v_ st=v_ A+c_ s”: Here v_ st = √(v_ A^2+c_ s^2), the fastest MHD wavespeed (which has no particular physical motivation but resembles what might be inferred in hydrodynamic models or observations where a plasma β∼1 is simply assumed).


Fig. <ref> summarizes the results: the effects of this choice are much smaller than the variations of e.g. κ discussed in the main text. Most importantly, we see no difference in averaged L_γ/L_ SF (although there are significant fluctuations due to stochasticity), implying that streaming at these speeds – even if modestly super-Alfvenic streaming is allowed – is ineffective at transporting CRs from dense regions. Of course, if we limit the “streaming losses” to scale with v_ A·∇ P_ cr and continue to increase v_ st without limit, it  will eventually become “fast enough.” In fact, our Eq. <ref> shows that, given the manner in  which we approximate CR transport numerically, κ and v_ st are formally degenerate if we replace κ with v_ st→κ |B̂·∇ e_ cr|/(e_ cr+P_ cr). From this, using a typical CR gradient scale length ∼ 1-10kpc in our galaxies and favored κ∼ 3×10^29  cm^-2 s^-1, we see that v_ st∼ 100-1000  km s^-1≫ v_ A (and ≫ c_s in most of the ISM) is required.

We have run the similar tests with our L⋆ galaxy model (not shown) and find qualitatively identical results (with even smaller differences between streaming models).











§ NUMERICAL TESTS




Our CR transport implementation is described in  <ref>. Here we present some numerical tests, including simple code validation problems, tests of the effect of the maximum CR free-streaming speed (or “reduced speed of light” c̃), comparison of our two-moment implementation to zeroth-moment “pure diffusion” solvers, and numerical resolution studies.






 §.§ CR Shocktube Test



We test our code implementation of CR coupling to adiabatic and advective terms and in the MHD Riemann problem using a variation of the <cit.> shocktube presented in <cit.>. A 3D box of of dimensions 64×16×16 is full of gas (adiabatic index 7/5) and CRs. Half the box has initial (ρ, v_x/c_s, P_ gas, P_ cr)=(1,1,1,2) and the other half has (0.25,1.5,0.275,0.275). We also consider another shocktube with (1,1,1,2) and (0.25,10,0.00384,0.00384). In both our mass resolution is 0.004. CRs have no diffusion or streaming (just advection). 

Fig. <ref> compares with analytic solutions from <cit.> at t=5. The agreement is good (despite very small shock broadening and numerical oscillations near discontinuities). The small “bumps” on the left close to x=20 are due to the slope limiter and should converge away at a higher resolution in our MFM method (see the hydrodynamic Sod shocktube test in ).












 §.§ Idealized Diffusion Test



We now test a “diffusion” problem to validate the diffusive terms in our CR transport implementation, i.e. a problem where the gas does not move or respond to CRs (we disable the terms by which the CRs act on the gas), but the CRs are transported (F̃_ cr 0). The corresponding “two moment” equation is 

    ∂ e_ cr/∂ t +κ/(γ_ cr-1)c̃^2∂^2 e_ cr/∂ t^2  = κ∇^2e_ cr,

which is a telegraph equation.

We initialize a 1D Gaussian distribution in e_ cr^0≡ e_ cr(r, t=0), centered at r=0, with total energy E_ cr=1, width σ=0.5 kpc at the center of a 5-kpc cube with 2048 resolution elements, set constant κ=3×10^27  cm^2 s^-1, and set c̃=100  km s^-1 or 500  km s^-1, and evolve the system for 5Myr. We do not include magnetic fields so the diffusion is isotropic. Given the symmetry of the problem, this can be solved exactly via the usual separation of modes, giving solutions of the form:

    e_ cr( r, t) =   ∫ d k  a^±_ k exp{i [ k· r - ω^±_ k t ] },

where κω_ k^2/(γ_ cr-1)/c̃^2 + i ω_ k = κ k^2 and ω^±_ k are two roots of the previous equation,

    ω^±_ k=-i±√(4 κ^2k^2/(γ_ cr-1)/c̃^2  -1)/2 κ/(γ_ cr-1)/c̃^2.


In our “diffusion” test, the initial CR flux is set to zero. Together with the initial CR energy density, we can solve for

    a^±_ k = 1/(2π)^Dexp(-σ^2k^2/2)/1-ω_ k^±/ω_ k^∓,

where D is the dimension of the Gaussian packet. With Eqs. <ref> and <ref>, the time evolution of CR energy density can then be calculated by integrating Eq. <ref> numerically.


This problem is entirely scale-free, and we can transform to solutions with any other value of κ via suitable re-scaling. As c̃^2 t/κ→∞, the solutions progressively approach the solution of the pure diffusion equation[Specifically, at this limit, the a^-_ k term becomes exponentially small and the a^+_ k term approaches the “pure diffusion” solution. Thus the solution to the two moment equation converges to the “pure diffusion” solution insensitive to the initial condition.], which is:

    e_ CR( r, t, c̃→∞) →exp[ -r^2/2(2κ t +σ^2)]/[2π (2κ t+ σ^2)]^3/2.


Fig. <ref> shows the results of our simulation for varying c̃ at fixed κ and t (i.e. varying the dimensionless parameters c̃^2 t/κ and r c̃/κ which determine the behavior of the problem)[We have turned off the HLL flux in the simulation to avoid small numerical diffusion.]. In all cases, the agreement with the exact solution is excellent, with numerical integration errors less than one percent. This validates that our two-moment implementation correctly solves the desired diffusion problem.

Moreover, Fig. <ref> also gives us a practical estimate – for typical units and spatial scales of our simulations – of the rate at which solutions with lower c̃ converge to the solution with c̃ = c (the speed of light). In Fig. <ref>, the test with c̃=500  km s^-1 is already effectively indistinguishable from the “pure diffusion” solution. Even at c̃=100  km s^-1, the solutions differ only by less than 10%. 

Of course, in real problems with bulk gas flows (e.g. galactic rotation), such a slow maximum CR transport speed would mean CRs would lag advection, in an unphysical manner, which motivates our additional tests below.














 §.§ Comparison of Zeroth and Two Moment Approximations in Galaxy Simulations



In Figure <ref>, we compare results of a full galaxy simulation using, instead of our default two-moment expansion (where we explicitly evolve the CR flux F̃_ cr with a finite maximum free-streaming speed c̃, as discussed in  <ref>), the “zeroth moment” or “pure diffusion” method  (for detailed tests of our pure isotropic/anisotropic diffusion solver, we refer to ). In the equations of motion, the latter simply replaces the explicitly-evolved flux vector with the value F̃_ cr = -κ ∇ e_ cr. This is mathematically equivalent to taking c̃→∞, and imposes a number of numerical difficulties discussed in  <ref> (not least of which is a much smaller timestep limit). However, Fig. <ref> shows there are only small differences in SFR and L_γ/L_ SF between these two approaches. We find the same in all other galaxy properties studied here.








 §.§ Anisotropic Diffusion Test


 We next test our scheme with an intrinsically multi-dimensional anisotropic diffusion test. We consider an initial 3D (256 particles on a side) spherically symmetric Gaussian profile with an anisotropic diffusion coefficient κ=3× 10^28 cm^2/s and a fixed B field pointing in the x direction. We set c̃=1000  km/s, initial F_ cr=0, and evolve over 5 Myr with the zeroth moment (“pure diffusion”) and two moment schemes. Fig. <ref> show the resultant CR energy densities through the x or y planes. We found a good agreement between two schemes and that CRs only diffuse along the B field direction (if there is no initial F_ cr perpendicular to the B field; see <ref>). We have also checked that the total CR energy is conserved (with deviations smaller than one in a thousand) in both schemes. 








 §.§ Test of a variation in the CR flux equation (Eq. <ref>)



Eq. <ref> governs the CR flux evolution, and obviously approaches the desired anisotropic diffusion+streaming equation (with F̃_  cr→ -κ_∗ ∇_e_ cr) when c̃→∞ is large, or in local quasi-steady-state (∂F̃_ cr/∂ t→ 0), or on spatial scales > κ/c̃ (timescales >κ/c̃^2). However, on small spatial/time scales when out-of-equilibrium, small fluxes perpendicular to the magnetic fields can arise from either integration error or rapid small-scale fluctuations in the magnetic field direction[We thank the anonymous referee pointing out this potential issue.]. If we divide the flux into parallel and perpendicular components F̃_ cr = F̃_ cr, + F̃_ cr,⊥, then we see the perpendicular component obeys c̃^-2 ∂F̃_ cr,⊥/∂ t = -[(γ_ cr-1)/κ_∗] F̃_ cr,⊥ (taking v=0, for simplicity, though this  does not change our conclusions below). So an initial F̃_ cr,⊥ can formally propagate but will be exponentially damped on a timescale ∼κ_∗/c̃^2. For κ^* needed to match the γ ray luminosity of nearby galaxies (∼ 10^29) and c̃∼ 1000  km/s, this timescale is ∼ 0.2 Myr, much shorter than the hadronic interaction time in the MW's ISM (n_ ISM∼ 1  cm^-3), so our γ ray constraint will not be affected with or without F̃_ cr,⊥ (well within the uncertainties). Likewise since this is much shorter than relevant dynamical times, essentially no effects of CR pressure and other dynamics should be altered.






If desired, we can trivially set F̃_ cr,⊥=0 every timestep, and evolve only the parallel component of Eq. <ref>: 

    1/c̃^2[ ∂F̃_ cr/∂ t + ∇·( v⊗F̃_ cr) ]_ + ∇_ P_ cr  = -(γ_ cr-1)/κ^∗ F̃_ cr,

where X_≡ (B̂⊗B̂)· X. We have implemented Eq. <ref> and compare it (in Fig. <ref>) to our default Eq. <ref> in an anisotropic  diffusion test (as  <ref>) with an initial (intentionally erroneous) pure-perpendicular flux F̃_ cr = F̃_  cr,⊥ = (2000  km/s)e_ crŷ. This shows that the “original” Eq. <ref> does formally allow CRs to propagate perpendicular to the B field, although the effect is rapidly damped, while the “modified” Eq. <ref> does not.
 
Fig. <ref> compares two galaxy simulations (Dwarf and L⋆ galaxy, with the latter at one level lower resolution than the fiducial  main-text case) run with MHD (anisotropic diffusion), κ_=3e28, and either Eq. <ref> or Eq. <ref>. We find essentially no difference in their star  formation rates or gamma-ray luminosities (apart from a small stochastic deviation associated with one slightly-stronger burst in the dwarf run).  
This is expected given the arguments above. Moreover, note that all the differences between Eq. <ref> and Eq. <ref> appear  on  scales ≪κ/c̃; but  <ref> showed that taking c̃→∞ does not change our conclusions.


Of course, CRs can physically propagate across mean magnetic fields due to unresolved (micro-scale) magnetic turbulence <cit.>, so a small F̃_ cr,⊥ might not be unreasonable. But in the main text, we show that even much larger perpendicular diffusivities given  by assuming isotropic diffusion (much  larger than the fluxes that arise from  the numerical or physical effects described above) do not strongly alter our conclusions regarding either the effects of CRs or the observationally-favored effective κ.


It is worth noting that several  recent studies including this one and <cit.>, <cit.>, and <cit.> have adopted slightly-different forms of the CR flux equation (although the CR energy equation is consistent in all cases). If we impose some desired scattering (κ_∗), all of these can be written (after some algebra) in the form:

    1/c̃^2 𝔻_tF̃_ cr + ∇_P_ cr   =  -(γ_ cr-1)/κ_∗ F̃_ cr,

where the differences are contained in the operator 𝔻_t.  For our default Eq. <ref>, 𝔻_t𝐅̃_ cr = ∂𝐅̃_ cr/∂ t + ∇· ( v⊗𝐅̃_ cr).  For the zeroth-moment or “pure diffusion+streaming” equation in  <ref>, 𝔻_t𝐅̃_ cr=0. For <cit.>, 𝔻_t𝐅̃_ cr = [∂{𝐅̃_ cr+ v (e_ cr+P_ cr) }/∂ t]_, and for <cit.>, 𝔻_t𝐅̃_ cr  = 𝐅̂̃̂_ cr [∂|𝐅̃_ cr|/∂ t + ∇·( v |𝐅̃_ cr|) + 𝐅̃_ cr·{(𝐅̂̃̂_ cr·∇) v}]   =  ∂𝐅̃_ cr/∂ t + ∇· ( v⊗𝐅̃_ cr ) + (𝐅̃_ cr·∇) ( v_- v_⊥), with F̃_ cr,⊥=0.[For <cit.>, we have formally taken the limit of the diffusion/interaction tensor 𝕂≡κ_∗ (B̂⊗B̂) + κ_⊥ (𝕀 - B̂⊗B̂) (𝕀 is the identity tensor) as κ_⊥→ 0. For <cit.>, we use 𝐅̃_ cr =|𝐅̃_ cr | 𝐅̂̃̂_ cr =f_ cr B̂ in their notation, and replace their scattering term ∂ f/∂ t|_ scatt (which explicitly attempts to account for dynamically-evolving, forward-and-backward propagating gyro-scale Alfven waves) with -[(γ_ cr-1)/κ_∗] F̃_ cr for a locally-stationary κ_∗.]


What is important for our purposes here is to note that all of these expressions differ only inside the term suppressed by c̃^-2. This means they all relax to the same diffusion+streaming equation in quasi-steady state and/or as c̃→∞, and differ only on scales below the CR mean  free path  ∼κ/c̃. To the extent that our galaxy-scale results are converged with  respect  to the value of c̃, and are not changed if we take c̃→∞ ( <ref>), these differences in 𝔻_tF̃_ cr cannot alter  our conclusions. Moreover, if  we adopt a “reduced speed of light” c̃≪  c, then by definition the flux equation does not exactly represent reality on scales ≪κ/c̃, regardless  of the form of 𝔻_t (i.e. the regime where these expressions differ is exactly where they all  become inexact when c̃≪ c is adopted). We have explicitly verified that our conclusions are robust to these choices of flux equation, running limited galaxy-scale simulations with the forms of  𝔻_tF̃_ cr from <cit.> or <cit.> as defined above: the results are similar to those  in Fig. <ref>.








 §.§ Comparison of Different (Finite) Maximum CR Propagation Speeds in Galaxy Simulations



We now examine the effect of varying c̃ on a full simulation (our L⋆ model). For c̃≳ 500  km s^-1, which is generally faster than the bulk rotation and outflow speeds (at least those containing most of the gas) in the galaxies, Fig. <ref> shows there is a small impact of c̃ on the SFR (we find the same for all other galaxy properties, not shown here). For L_γ/L_ SF, we find almost no impact of c̃ in simulations with κ< 10^29 cm^2/s for any values c̃≳ 500  km s^-1.

Because the effective “advective velocity” of CRs under pure diffusion is ∼κ / ℓ (where ℓ is some gradient scale-length), at much larger diffusion coefficients, e.g. our κ=3×10^29  cm^2 s^-1, where most of the CRs escape diffusively, the value of L_γ is slightly larger for c̃=500  km s^-1 compared to much-higher c̃ (because escape is slightly slower). But once c̃≥ 1000  km s^-1, we see no detectable difference. Moreover by the latter half of the time we run for, the differences even for c̃= 500  km s^-1 (compared to c̃∼ 4000  km s^-1) at κ=3×10^29  cm^2 s^-1 are factors of ∼ 1.4, not large enough to change any of our conclusions.

In tests run for shorter duration and tests of our Dwarf galaxy (not shown), we have also verified similar conclusions, and found that runs with magnetic fields (since these slow down the transport) and finite streaming velocities (since these dominate the transport over diffusion in some regimes) exhibit even weaker dependence on c̃ within the range c̃∼ 500-4000  km s^-1, even at κ=3×10^29  cm^2 s^-1.











 §.§ Idealized Streaming Test


Here we test our streaming implementation with an initial 1D symmetric triangular CR profile with e_ cr = 2-|x| (and F_ cr= v_ stre_ cr), where we disable the gas motion and CR streaming loss (v_ st·∇ P_ cr), but CRs can move across gas particles with  a constant v_ st=1  km/s. In the numerical test, we consider evenly spaced 2048 particles over 10 kpc, the reduced speed of light  c̃=1000  km/s, and a small diffusion coefficient 3× 10^22  cm^2/s.[We include a very small diffusion coefficient to avoid potential (numerical) overflows in Eq. <ref>, i.e. when the CR pressure profile is flat or, equivalently, ∇ P_ cr is huge.] Here, we also consider a variant test  in which we relax the limiter on CR pressure gradient[The limiter is necessary to avoid divergence and other numerical issues in realistic galaxy simulations.] in order to correctly capture the plateau region and also turn off the HLL flux to avoid numerical diffusion smoothing the discontinuity (so  we can  cleanly separate the effects of numerical dissipation, which are resolution-dependent, and the effects of the actual form of the flux equation and its implementation of streaming).[Because the exact solution to this test contains dis-continuous first derivatives, any numerically stable method will introduce some dissipation at the cusps, even at infinite resolution.] We also show the “original” test run with the original limiter and the HLL flux included, to illustrate the version that we used in the main text. 



The analytic solution, calculated in <cit.>, has a flat CR distribution between x = ± x_ m and two inclined distributions for |x| > x_ m &  |x| < 1  with

    e_ cr(x, t) = 2 + 4v_ stt-|x|,

and x_ m is determined from the energy conservation:

    x_ m = √( ( 1+4/3v_ stt  )^2+8v_ stt/3-1).


Fig. <ref> shows the analytic solution agrees well with our numerical solution at t = 10  Myr, although the “original” run has a round top rather than a plateau mainly because of the limiter on the CR pressure gradient.




 §.§ Resolution study


In Fig. <ref>, we show the properties of L⋆ Galaxy runs, at three different resolution levels. Our baseline is the fiducial resolution listed in Table <ref>, but we compare runs with 10x and 100x poorer mass resolution. Even at 10x poorer resolution, we find very similar SFRs, L_γ/L_ SF, and all other galaxy properties studied here; the same is true in the Dwarf runs (not shown). However systematic offsets do begin to appear at 100x poorer resolution. Other more detailed properties (e.g. phase structure of galactic winds) may require much higher resolution - this will be explored in future work.




§ COMPARISON WITH L_ 0.1-100GEV OBSERVATIONS



For completeness, we also compare our simulations with the broader, 0.1-100 GeV, observed energy range of γ rays from <cit.> and <cit.>.

In the main text, we estimate the L_ > 1 GeV following the L11 approach with the hadronic loss rate taken from <cit.> (see Eq. <ref> and the related text), with an implicit assumption that all CR energy is at > 1GeV to simplify the calculation. However, to properly account for the lower energy γ rays and enable direct comparison to both the observations and simulations from <cit.>, we follow <cit.> to calculate  L_ 0.1-100 GeV.


We assume that CRs follow a power-law spectrum with a constant spectral index 2.2 (as used in the main text) and the low momentum cutoff q = P_p/(m_pc)=1, where m_p is the proton mass. The normalization is determined by integrating the spectrum over energy and comparing to the local CR energy density. Then we calculate the energy-integrated γ-ray emissivity from pion decay for the energy band 0.1-100 GeV with Eq. 6 in <cit.>, and then integrate to get L_ 0.1-100 GeV. While the secondary IC emission from high energy CR electrons also contributes to γ-ray luminosity in this range, we find it is relatively unimportant compared to emission from pion decay, since (1) in high B field regions, CR electrons cool preferentially through synchrotron radiation; (2) in low B field regions, IC luminosity is only a small fraction (< 20%) of the γ-ray luminosity. Thus,
we simply neglect its contribution in Fig.<ref>. 

Fig. <ref> shows that κ∼ 1 e29-3 e29  cm^2/s provides the best match to the observations for Dwarf and L⋆ Galaxy runs,
consistent with our conclusions in the main text.








