Cosmic ray feedback in the FIRE simulations: constraining cosmic ray propagation with GeV gamma ray emission

TL;DR

This paper implements cosmic ray feedback in FIRE simulations, compares different propagation models, and constrains them using GeV gamma-ray observations, revealing key insights into cosmic ray escape and galaxy emission properties.

Contribution

It introduces a novel simulation approach combining CR diffusion and streaming, and constrains CR propagation models using gamma-ray data across galaxy types.

Findings

Constant isotropic diffusion models match gamma-ray observations.

Advection-only and streaming-only models overpredict gamma-ray luminosity.

CRs escape low-density galaxies before significant losses, starbursts are CR calorimeters.

Abstract

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\star$ starburst, and $L\star$ 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 $\gamma$-ray emission from nearby and starburst galaxies. We reproduce the $\gamma$-ray observations of dwarf and $L\star$ galaxies with constant isotropic diffusion coefficient $\kappa \sim 3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$. Advection-only and streaming-only models produce order-of-magnitude too large $\gamma$-ray luminosities in dwarf and $L\star$ galaxies. We show that in models that match the $\gamma$-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 $\gamma$-ray emissivities. Models where CRs are ``trapped'' in the star-forming disk have lower star formation efficiency, but these models are ruled out by $\gamma$-ray observations. For models with constant $\kappa$ that match the $\gamma$-ray observations, CRs form extended halos with scale heights of several kpc to several tens of kpc.

Tables

Table 1. Table 1: Simulation parameters. M vir subscript 𝑀 vir M_{\rm vir} is the virial mass; c 𝑐 c is the halo concentration; M ∗ , disk subscript 𝑀 disk M_{*,{\rm disk}} is the mass of stellar disk; M ∗ , bulge subscript 𝑀 bulge M_{*,{\rm bulge}} is the mass of stellar bulge; M g , disk subscript 𝑀 g disk M_{\rm g,disk} is the mass of gas disk; M g , halo subscript 𝑀 g halo M_{\rm g,halo} is the mass of gas halo; d ∗ , disk subscript 𝑑 disk d_{\rm*,disk} is the stellar disk radial scale length; h ∗ , disk subscript ℎ disk h_{\rm*,disk} is the thickness of stellar disk; d g , disk subscript 𝑑 g disk d_{\rm g,disk} is the gas disk radial scale length; m b subscript 𝑚 b m_{\rm b} is the gas particle mass.

Name	$M_{\mathrm{vir}}$	$R_{\mathrm{vir}}$	c	$M_{\ast ,\mathrm{disk}}$	$M_{\ast ,\mathrm{bulge}}$	$M_{g,\mathrm{disk}}$	$M_{g,\mathrm{halo}}$	$d_{\ast ,\mathrm{disk}}$	$h_{\ast ,\mathrm{disk}}$	$d_{g,\mathrm{disk}}$	$m_b$
	[${10}^{10}{M}_{\odot }$]	$\left[10\mathrm{kpc}\right]$		[${10}^{10}{M}_{\odot }$]	[${10}^{10}{M}_{\odot }$]	[${10}^{10}{M}_{\odot }$]	[${10}^{10}{M}_{\odot }$]	[kpc]	[kpc]	[kpc]	[${10}^3{M}_{\odot }$]
Dwarf	2.9	63	15	0.019	0.0014	0.1	0.01	1.0	0.2	5	3.3
Starburst	21	121	11	0.57	0.14	4.0	1.0	1.0	0.2	15.0	20.0
$L\star$ Galaxy	150	234	12	4.7	1.5	0.9	0.1	3.2	0.24	6.4	2.6

Table 2. Table 2: Different propagation models of CRs. Each column gives the name of our simulation models, while rows list the physics/parameters of the propagation model. The “MHD” column row indicates whether magnetic fields are included. The “Streaming” column indicates whether CR streaming is considered. κ 𝜅 \kappa gives the isotropic/parallel CR diffusion coefficient (CRs will diffuse isotropically if MHD is off, while CRs will diffuse along magnetic fields if MHD is on). c ~ ~ 𝑐 \tilde{c} is the “reduced speed of light” in the two-moment method (see §  2.2.3 ).

	Hydro	MHD	Advec-	$\kappa =3e27$	$\kappa =3e28$	$\kappa =3e29$	MHD	MHD	MHD
	no CR	no CR	-tion					$\kappa =3e28$	$\kappa =3e28$
							Streaming		Streaming
MHD	Off	On	Off	Off	Off	Off	On	On	On
Streaming	Off	Off	Off	Off	Off	Off	On	Off	On
$\kappa$   [${\mathrm{cm}}^2/s$]	-	-	0	$3\times {10}^{27}$	$3\times {10}^{28}$	$3\times {10}^{29}$	0	$3\times {10}^{28}$	$3\times {10}^{28}$
$\tilde{c}$   [km/s]	-	-	-	500	1000	2000	1000	1000	1000

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\right)=0,$$

$$\frac{\partial \rho \mathbf{v}}{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\otimes \mathbf{v}+P_T\mathbb{I}-\mathbf{B}\otimes \mathbf{B}\right)=0,$$

$$\frac{\partial \rho e}{\partial t}+\nabla \cdot \left[(\rho e+P_T)\mathbf{v}-(\mathbf{v}\cdot \mathbf{B})\mathbf{B}\right]$$

$$=P_{\mathrm{cr}}\nabla \cdot \mathbf{v}+{\Gamma }_{\mathrm{st}}+S_{\mathrm{g}}-{\Gamma }_{\mathrm{g}},$$

$$\frac{\partial \mathbf{B}}{\partial t}+\nabla \cdot \left(\mathbf{v}\otimes \mathbf{B}-\mathbf{B}\otimes \mathbf{v}\right)=0,$$

$$\frac{\partial e_{\mathrm{cr}}}{\partial t}+\nabla \cdot {\mathbf{F}}_{\mathrm{cr}}=\mathbf{v}\cdot \nabla P_{\mathrm{cr}}-{\Gamma }_{\mathrm{st}}+S_{\mathrm{cr}}-{\Gamma }_{\mathrm{cr}},$$

$$\displaystyle\frac{DE_{\rm cr}}{Dt}=-\int_{\Omega}d^{3}{\bf x}\,{\Bigl{\{}}$$

$$\frac{1}{{\tilde{c}}^2}\left[\frac{\partial {\tilde{\mathbf{F}}}_{\mathrm{cr}}}{\partial t}+\nabla \cdot \left(\mathbf{v}\otimes {\tilde{\mathbf{F}}}_{\mathrm{cr}}\right)\right]+{\nabla }_{\Vert }{P}_{\mathrm{cr}}$$

$${\kappa }^{\ast }=\kappa +\frac{v_{\mathrm{st}}(e_{\mathrm{cr}}+P_{\mathrm{cr}})}{|\hat{\mathbf{B}}\cdot \nabla e_{\mathrm{cr}}|},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\left(V\mathbf{U}\right)}_i=-\sum _j{\mathbf{F}}_{\mathrm{diff},\mathrm{ij}}^{\ast }\cdot {\mathbf{A}}_{i,j},$$

$${\Gamma }_{\mathrm{cr}}={\tilde{\Lambda }}_{\mathrm{cr}}{e}_{\mathrm{cr}}{n}_{\mathrm{n}}=({\tilde{\Lambda }}_{\mathrm{cr},\mathrm{had}}+{\tilde{\Lambda }}_{\mathrm{cr},\mathrm{Cou}})e_{\mathrm{cr}}{n}_{\mathrm{n}}$$

$$=$$

$$S_{\mathrm{gas}}=$$

$$\left(\frac{n_{\mathrm{n}}}{\mathrm{c}{\mathrm{m}}^{-3}}\right)\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}{\mathrm{s}}^{-1}$$

$$L_{\gamma }(\ge \mathrm{GeV})\approx \sum _i\frac{1}{3}{\beta }_{\pi }{\tilde{\Lambda }}_{\mathrm{cr},\mathrm{had}}{e}_{\mathrm{cr}}{n}_{\mathrm{n}}\Delta V_i,$$

$$\frac{{\dot{E}}_{\mathrm{cr},\mathrm{SN}}}{[erg/s]}={ϵ}_{\mathrm{cr}}\frac{E_{\mathrm{SN}}}{[\mathrm{erg}]}\frac{[{\mathrm{M}}_{\odot }]}{<m_{\ast }>}\xi (m_{\ast }>8{\mathrm{M}}_{\odot })\frac{\mathrm{SFR}}{3.2\times 10^7[{\mathrm{M}}_{\odot }/yr]},$$

$$L_{\gamma }(\ge \mathrm{GeV}{)}^{\mathrm{calor}}\approx 6.7\times 10^{39}\frac{\mathrm{SFR}}{[{\mathrm{M}}_{\odot }/yr]}erg/s.$$

$$L_{\mathrm{SF}}^{\mathrm{constant}}\approx 3.5\times 10^{43}\frac{\mathrm{SFR}}{[{\mathrm{M}}_{\odot }/yr]}erg/s,$$

$$\frac{\partial e_{\mathrm{cr}}}{\partial t}+\frac{\kappa }{({\gamma }_{\mathrm{cr}}-1){\tilde{c}}^2}\frac{{\partial }^2{e}_{\mathrm{cr}}}{\partial t^2}=\kappa {\nabla }^2{e}_{\mathrm{cr}},$$

$$e_{\mathrm{cr}}(\mathbf{r},t)=$$

$${\omega }_{\mathbf{k}}^{\pm }=\frac{-i\pm \sqrt{4{\kappa }^2{k}^2/({\gamma }_{\mathrm{cr}}-1)/{\tilde{c}}^2-1}}{2\kappa /({\gamma }_{\mathrm{cr}}-1)/{\tilde{c}}^2}.$$

$$a_{\mathbf{k}}^{\pm }=\frac{1}{(2\pi {)}^D}\frac{\exp (-{\sigma }^2{k}^2/2)}{1-{\omega }_{\mathbf{k}}^{\pm }/{\omega }_{\mathbf{k}}^{\mp }},$$

$$e_{\mathrm{CR}}(\mathbf{r},t,\tilde{c}\to \infty )\to \frac{\exp \left[-\frac{r^2}{2(2\kappa t+{\sigma }^2)}\right]}{[2\pi (2\kappa t+{\sigma }^2){]}^{3/2}}.$$

$$\frac{1}{{\tilde{c}}^2}{\left[\frac{\partial {\tilde{\mathbf{F}}}_{\mathrm{cr}}}{\partial t}+\nabla \cdot (\mathbf{v}\otimes {\tilde{\mathbf{F}}}_{\mathrm{cr}})\right]}_{\Vert }+{\nabla }_{\Vert }{P}_{\mathrm{cr}}=-\frac{({\gamma }_{\mathrm{cr}}-1)}{{\kappa }^{\ast }}{\tilde{\mathbf{F}}}_{\mathrm{cr},\Vert }$$

$$\frac{1}{{\tilde{c}}^2}{\mathbb{D}}_t{\tilde{\mathbf{F}}}_{\mathrm{cr}}+{\nabla }_{\Vert }{P}_{\mathrm{cr}}$$

$$e_{\mathrm{cr}}(x,t)=2+4v_{\mathrm{st}}t-|x|,$$

$$x_{\mathrm{m}}=\sqrt{{\left(1+\frac{4}{3}{v}_{\mathrm{st}}t\right)}^2+\frac{8v_{\mathrm{st}}t}{3}-1}.$$

Full text

1 Introduction

Cosmic rays (CRs) are charged particles moving with relativistic speeds, mainly generated through shock acceleration of supernova remnants (SNRs) (Bell, 2004) (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 (Strong et al., 2007; Zweibel, 2013; Grenier et al., 2015), and transfer energy to gas via Coulombic and hadronic interactions (Mannheim & Schlickeiser, 1994; Enßlin et al., 2007; Guo & Oh, 2008). Their energy density is comparable to thermal and magnetic energies in the solar neighborhood (Ginzburg & Ptuskin, 1985; Boulares & Cox, 1990), so CRs are believed to be dynamically important in galaxy evolution.

The impacts of CRs on galaxy evolution have been studied with analytic models (Ipavich, 1975; Breitschwerdt et al., 1991, 1993; Zirakashvili et al., 1996; Socrates et al., 2008; Everett et al., 2008; Dorfi & Breitschwerdt, 2012; Mao & Ostriker, 2018) and idealized and cosmological simulations (e.g. Jubelgas et al. 2008; Uhlig et al. 2012; Booth et al. 2013; Wiener et al. 2013; Hanasz et al. 2013; Salem & Bryan 2014; Salem et al. 2014; Chen et al. 2016; Simpson et al. 2016; Girichidis et al. 2016; Pakmor et al. 2016; Salem et al. 2016; Wiener et al. 2017; Ruszkowski et al. 2017; Butsky & Quinn 2018; Farber et al. 2018; Jacob et al. 2018; Girichidis et al. 2018). 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 (Parker, 1992; Hanasz et al., 2009; Kulpa-Dybeł et al., 2011, 2015).

Despite its importance, the details of CR propagation are uncertain. The most popular CR propagation models are self confinement and extrinsic turbulence (Zweibel 2013, 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 (Zweibel, 2013). These two mechanisms confine and isotropise the CR distribution explaining the remarkably low CR anisotropy observed from the Earth (see, e.g. Hillas & Ouldridge 1975) and the long residence time ($>10^{7}\;{\rm yr}$) inferred from the ratios between stable primary and secondary nuclei (Strong et al., 2007). Their long confinement time and small anisotropy imply that CRs have short mean free paths ($\sim{\rm 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. Zweibel 2017, 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 (Ginzburg & Ptuskin, 1976). 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. GALPROP111 https://galprop.stanford.edu/, which attempt to synthesize observational constraints on the MW gas density distribution, CR abundances and spectra, $\gamma$-ray and radio emission, and theoretical models for e.g. galactic winds and diffusive re-acceleration (Strong & Moskalenko, 1998, 2001). These models are commonly used to constrain the “effective” isotropic-equivalent diffusion coefficient of CRs averaged over the whole MW disk (e.g. Strong et al., 2007; Trotta et al., 2011). 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 $\gamma$-ray emission can serve as an independent constraint on CR propagation (Ackermann et al., 2012a; Strong et al., 2000, 2004). High energy ($>{\rm GeV}$) CRs collide with nuclei in the interstellar medium (ISM) and produce pions, which decay into ${\rm GeV}$ $\gamma$-rays. Since pionic $\gamma$-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 AMS Collaboration et al. 2002), CR distribution and propagation can be constrained with high energy $\gamma$-ray observations.

Recently, $\gamma$-ray emission was observed from Local Group (Abdo et al., 2010a, b, c) and starburst galaxies (Acero et al., 2009; VERITAS Collaboration et al., 2009; Abdo et al., 2010d; H. E. S. S. Collaboration et al., 2018), which can be used as a probe of CRs beyond the solar system and the Milky Way (MW). Abdo et al. (2010b) found a correlation between $\gamma$-ray emission and star formation rate (SFR) with a steeper than linear drop at low SFRs ($E_{\gamma}\propto{\rm SFR}^{1.4\pm 0.3}$; summarized in Ackermann et al. 2012b).

To explain this correlation, Lacki et al. (2011) (hereafter L11) 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 (Lacki et al., 2010), the model was used to estimate pionic $\gamma$-ray luminosities of galaxies. They found that in order to explain the correlation between $\gamma$-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 Thompson et al. 2007; see also Abramowski et al. 2012; Yoast-Hull et al. 2013; Wang & Fields 2018), 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 (Kennicutt, 1998), 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. Lopez et al. 2018).

In this study, we investigate the impact of CRs on dwarf, sub-$L\star$ starburst, and $L_{\star}$ 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 (similar to Jiang & Oh, 2018). 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 (Hopkins et al., 2018c, b).222http://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 (Hopkins et al., 2018c), galactic winds (Muratov et al., 2015; Anglés-Alcázar et al., 2017; Hafen et al., 2019), cored central dark matter profiles (Oñorbe et al., 2015; Chan et al., 2015; Wetzel et al., 2016; Fitts et al., 2017), the mass-metallicity relation (Ma et al., 2016) and spatial distribution of gas and metals within galaxies and the CGM (Faucher-Giguère et al., 2015, 2016; Ma et al., 2017; Muratov et al., 2017; Hafen et al., 2017), typical galaxy star formation rates and histories (Sparre et al., 2017), the Kennicutt-Schmidt law (Orr et al., 2018), and galactic magnetic field structure (Su et al., 2018).

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 $\sim$GeV $\gamma$-ray emission from galaxies (and compare our findings with previous theoretical studies). § 2 and § 3 discuss numerical methods, simulated physics, and initial conditions. In § 4.1.1, we investigate how CRs and their propagation influence galactic properties. In § 4.2 we calculate the $\gamma$-ray emission from CRs in our simulations and compare with observational data. In § 5, 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 § 6.

2 Method

2.1 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 Hopkins et al. (2018c), so we only briefly review them here. Our simulations use the GIZMO333http://www.tapir.caltech.edu/p̃hopkins/Site/GIZMO code (Hopkins, 2015a) in its mesh-free Lagrangian finite mass (MFM) mode for (magneto)-hydrodynamics; extensive implementation details and tests of the MHD scheme are presented in Hopkins (2016); Hopkins & Raives (2016). GIZMO uses an updated version of the PM+Tree algorithm from Gadget-3 (Springel, 2005) to calculate gravity and adopts fully conservative adaptive gravitational softening for gas (Price & Monaghan, 2007). 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 Faucher-Giguère et al. (2009) model (accounting for local self-shielding). Stars form in locally self-gravitating, self-shielding, thermally Jeans-unstable gas444We 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 Hopkins et al. 2018c. at densities $n_{\rm H}\geq 100\,{\rm cm^{-3}}$. Once formed, we calculate the energy, momentum, mass and metal return for each star according to the STARBURST99 stellar population synthesis model (Leitherer et al., 1999), for a Kroupa (2002) 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 Hopkins et al. (2018c, b, a).

2.2 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 $\gamma_{\rm cr}=4/3$) in a “single bin” approximation.555One can think of this as evolving only the CR energy density at the energies $\gtrsim$ 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 $\rho$, velocity ${\bf v}$, magnetic field ${\bf B}$, and specific energy $e$, are extended with the equation for the CR energy density $e_{\rm cr}$ (McKenzie & Voelk, 1982):

[TABLE]

where $P_{\rm cr}=(\gamma_{\rm cr}-1)\,e_{\rm cr}$ is the CR pressure; $P_{T}$ is the total (thermal+magnetic+CR) pressure; $\Gamma_{\rm st}=-{\bf v}_{\rm st}\cdot\nabla P_{\rm cr}$ is the CR “streaming loss term” discussed below; $S_{\rm g}$ and $S_{\rm cr}$ are gas and CR source terms (e.g. injection); $\Gamma_{\rm g}$ and $\Gamma_{\rm cr}$ are gas and CR sink/loss (or “cooling”) terms; and ${\bf v}_{\rm st}$ is the CR streaming velocity. ${\bf F}_{\rm cr}$ is the CR energy flux, which can be written ${\bf F}_{\rm cr}=(e_{\rm cr}+P_{\rm cr})({\bf v}+{\bf v}_{\rm st})+{\bf F}_{\rm di}$ where the first term represents advection and streaming, whereas the second term is a diffusive-like flux (e.g. given by ${\bf F}_{\rm di}=-\kappa\hat{\bf B}\otimes\hat{\bf B}\cdot\nabla e_{\rm cr}$ in the “pure diffusion” or “zeroth moment” approximation, but we explicitly evolve this; see §2.2.3).

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})_{\rm eff}=\partial P/\partial\rho=(c_{s}^{2})_{\rm gas}+\gamma_{\rm cr}\,P_{\rm cr}/\rho$, but no other modifications to the MHD method is required.

2.2.1 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 ${\bf v}={\bf v}_{\rm 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. Uhlig et al. 2012):

[TABLE]

where $D/Dt=\partial/\partial t+{\bf v}\cdot\nabla$ is the Lagrangian derivative co-moving with the gas, and $E_{\rm cr}^{i}=\int_{\Omega_{i}}e_{\rm cr}\,d^{3}{\bf x}$ is the conserved total CR energy in the finite-volume domain $\Omega_{i}$ belonging to element $i$. Here $\tilde{\bf F}_{\rm cr}\equiv{\bf F}_{\rm cr}-{\bf v}\,(e_{\rm cr}+P_{\rm cr})={\bf v}_{\rm st}\,(e_{\rm cr}+P_{\rm cr})+{\bf F}_{\rm di}$. Pure advection with the gas is automatically handled in this description. In cosmological simulations, the Hubble flow is included in $\nabla\cdot{\bf v}$.

The ${P}_{\rm cr}(\nabla\cdot{\bf 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.666To ensure manifest energy conservation, this is solved when the mesh positions are updated. Calculating the volume changes $\Delta V_{i}=\int dt\,\int_{\Omega_{i}}\,d^{3}{\bf x}\,(\nabla\cdot{\bf v})$ with the kernel-weighted divergence of the fluid velocity field (which is the exact discrete change in the domain volume as defined in Hopkins 2015a), we have $\Delta E_{\rm cr}=-P_{\rm cr}\,\Delta 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 $\Gamma_{\rm st}=-{\bf v}_{\rm st}\cdot\nabla{P}_{\rm 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. Wentzel 1968; Kulsrud & Pearce 1969) which are damped and thermalize their energy effectively instantly (compared to our simulation timescales).777With the streaming velocity defined below, the streaming loss term can be written $DE_{\rm cr}/Dt=-E_{\rm cr}/\tau_{\rm st}$ with $\tau_{\rm st}^{-1}=(\gamma_{\rm cr}-1)\,|\hat{\bf B}\cdot\hat{\nabla}e_{\rm cr}|^{2}\,|v_{\rm st}\,\nabla e_{\rm cr}|/e_{\rm cr}$. When this is updated the resulting energy lost $\Delta E_{\rm cr}=\int dt\,\tau_{\rm st}^{-1}\,E_{\rm cr}$ is added to the gas thermal energy.

Finally, the $\int_{\Omega}d^{3}{\bf x}\,\nabla\cdot\tilde{\bf F}_{\rm 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, $\int_{\partial\Omega}d{\bf A}\cdot\tilde{\bf F}_{\rm 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 Hopkins 2015a for details).

We explicitly evolve the conserved quantities $E_{\rm cr}$ and total gas energy $E_{\rm gas}$ which are exchanged (either between gas elements or one another), ensuring manifest total energy conservation.888Because 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 $\lesssim 20\%$ of the injection), and negligible compared to physical non-conservative terms (e.g. collisional/radiative losses, injection).

2.2.2 The Streaming Velocity

CRs stream at some speed $v_{\rm 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. ${\bf v}_{\rm st}\equiv-v_{\rm st}\,\hat{\bf B}\,(\hat{\bf B}\cdot\hat{\nabla}P_{\rm cr})$ where $\hat{\nabla}P_{\rm cr}=\hat{\nabla}e_{\rm cr}=(\nabla P_{\rm cr})/|\nabla P_{\rm cr}|=(\nabla e_{\rm cr})/|\nabla e_{\rm 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_{\rm A}\,(=B/\sqrt{4\pi\rho}\,$) in the low-$\beta$ limit (see Skilling 1971; Holman et al. 1979, or more recently Kulsrud 2005; Yan & Lazarian 2008; Enßlin et al. 2011), so we adopt a fiducial value $v_{\rm st}=v_{\rm A}$.

But in the weak-field, plasma $\beta\gg 1$, regime, the streaming velocity can be boosted by significant wave damping (see discussion in Enßlin et al., 2011; Wiener et al., 2013; Ruszkowski et al., 2017), so we have also tested various streaming speeds in Appendix A. Although the streaming velocity can in principle exceed $v_{A}$ by a large factor, Wiener et al. (2013) and Ruszkowski et al. (2017) argued that the streaming loss $\Gamma_{\rm st}$ should be still limited by $\Gamma_{\rm st}=-{\bf v}_{\rm A}\cdot\nabla P_{\rm cr}$, because this term is mediated by the excitation of Alfven waves. Therefore, regardless of streaming velocities, we set the streaming loss to $-{\bf v}_{\rm A}\cdot\nabla P_{\rm cr}$. When streaming is disabled we simply eliminate terms relevant to streaming.

2.2.3 Diffusive Transport Terms: Two-moment Method

It is common in the literature to treat ${\bf F}_{\rm di}$ in the “zero-th moment” expansion, i.e. approximate it as an anisotropic scalar diffusion with ${\bf F}_{\rm di}=-\kappa\,\hat{\bf B}\,(\hat{\bf B}\cdot\nabla e_{\rm cr})$, where $\kappa$ 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: $\Delta t<C_{\rm cour}\,\Delta x^{2}/\kappa$, where $\Delta x$ is the resolution and $C_{\rm cour}$ the Courant factor) which becomes very small; (2) it implies unphysical super-luminal CR transport when the gradient-scale length $e_{\rm cr}/|\nabla e_{\rm cr}|$ becomes smaller than $\kappa/c\sim 3\,{\rm pc}\,(\kappa/3\times 10^{29}\,{\rm 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 (Sharma et al., 2010); and (5) it encounters the usual difficulties with anisotropic diffusion operators in moving-mesh codes described in Hopkins (2017) (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 Jiang & Oh (2018) and Thomas & Pfrommer (2019) (although the concept and use in CR dynamics are well-established; see Snodin et al. 2006 for examples). Rather than set $\tilde{\bf F}_{\rm cr}=-\kappa\,\nabla e_{\rm cr}$, we explicitly evolve the flux equation:

[TABLE]

where $\tilde{\bf F}_{\rm cr}$ is the CR flux measured in the frame comoving with the fluid, $\nabla_{\|}P_{\rm cr}\equiv\hat{\bf B}\otimes\hat{\bf B}\cdot\nabla P_{\rm cr}$, $\tilde{c}$ is the (reduced) speed of light, and $\kappa^{*}$ is the composite parallel (magnetic field-aligned) diffusion coefficient in the rest frame of the fluid,

[TABLE]

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 Hopkins (2017) with a few modifications.

We solve a general evolution equation of conserved quantities $\left(V\bf{U}\right)_{i}$ of cell i (e.g. CR energy) with a volume $V$ by summing over all adjacent cells j:

[TABLE]

where V is the cell volume and ${\bf F}^{*}_{\rm diff,ij}$ is the interface value of the flux and ${\bf A}_{i,j}$ is the effective face area defined in Hopkins (2015b), in the following steps:

We calculate all relevant coefficients, using the standard gradient estimator in GIZMO for MFM to estimate gradients, e.g. $\nabla P_{\rm cr}$ and $\nabla e_{\rm cr}$, as described in Hopkins (2015b); 2. 2.

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; 3. 3.

We compute the interface face value of the flux ${\bf F}^{*}_{\rm diff,ij}$ by solving the Riemann problem (RP) through the Harten et al. (1983) (HLL) method, using a MINMOD slope limiter, in an operator split manner from the pure MHD. 4. 4.

Finally, the source term in Eq. 3 (not considered in Hopkins 2017) is added implicitly to ensure stability.

We differ from Hopkins (2017) 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 $\tilde{c}$ (since we choose $\tilde{c}$ to be faster than other physical processes).

Unlike Jiang & Oh (2018), we do not modify the momentum transfer from CRs to gas, i.e. the second line of Eq. 2.2, 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 $\sim(1\,{\rm pc})(\kappa/10^{29}{\rm cm^{2}/s})$ are smaller than the resolved length-scales, so the strong coupling assumption is probably relevant in most of the situations.999 The formulation in Jiang & Oh (2018) 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 Jiang & Oh 2018).

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. 3 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 $\mathcal{O}(v/c)$. We discuss some of the subtle differences that can arise in Eq. 3 as a result, in Appendix B.5, 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 $\kappa$ and $v_{\rm st}$ on physical grounds (see e.g. Zweibel, 2017; Thomas & Pfrommer, 2019).

For now, if we ignore streaming, we see that in steady-state and/or when $\tilde{c}$ is large, or $\Delta t\gg\kappa/\tilde{c}^{2}$ (or on spatial scales $\gg\kappa/\tilde{c}$), this equation becomes ${\bf F}_{\rm cr}\approx-\kappa\,\hat{\bf B}\,(\hat{\bf B}\cdot\nabla e_{\rm cr})$, and we recover the usual diffusion equation (see Appendix B.6 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 $\tilde{c}$, and makes the timestep criterion $\Delta t<C_{\rm cour}\,\Delta x/\tilde{c}$,101010We adopt $C_{\rm cour}=0.4$ throughout, and have validated stability (as expected) for this value. $\Delta x$ in the Courant condition is defined in the same manner as Hopkins et al. (2018c) as the local mean inter-particle separation (i.e. the equivalent of the grid spacing in a regular-grid code), $\Delta x\equiv(m/\rho)^{1/3}$. which is only first-order, instead of quadratic, in $\Delta x$.

For true micro-physical CR motion, however, $\tilde{c}\approx 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 $\tilde{c}\ll c$ (namely the “reduced speed of light” approximation), since galaxy properties should still converge, regardless of $\tilde{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 $\tilde{c}\sim 500-2000\,{\rm km\,s^{-1}}$ are sufficient to give converged results, e.g. SFR and $\gamma$-ray emission (see Appendix B.6).

In Appendix B.3, we compare the results using the simpler pure-diffusion (zeroth-moment) approximation: we then simply assume ${\bf F}_{\rm di}\rightarrow-\kappa\,\hat{\bf B}\,(\hat{\bf B}\cdot\nabla e_{\rm cr})$ and solve the anisotropic diffusion equation (with the stricter Courant condition) as described in Hopkins (2017). This is equivalent to adopting $\tilde{c}\rightarrow\infty$, in our Eq. 3. For the same $\kappa_{\ast}$, this gives nearly-identical results to our default Eq. 3 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” $\tilde{c}$ (Levermore, 1984), while the “pure diffusion” approximation is akin to flux-limited diffusion (without the limiter).

2.2.4 The Diffusion Coefficient

The only remaining unspecified parameter in the CR treatment is the effective diffusion coefficient $\kappa$111111 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 (Wiener et al., 2013; Zweibel, 2013). 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 (Enßlin, 2003; Enßlin et al., 2007).

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, Trotta et al. (2011) constrained the isotropically-averaged diffusivity $\kappa$ to be $\sim 6\times 10^{28}{\rm cm^{2}/s}$ to within a factor of a few, at $\sim$GeV energy with GALPROP, using the measured energy spectra and abundances of nuclei species in CRs, and adopting a flat halo diffusion model (Ginzburg & Ptuskin 1976; 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 $\kappa$ and the CR halo height $z_{\rm h}$ (typically 1-10 kpc), out of which CRs can freely propagate (see Figure 3 in Trotta et al. 2011 or Figure 10 in Linden et al. 2010; this issue was also discussed in Ginzburg & Ptuskin 1976). Even in this model, it is possible to match the observational data with a significantly larger $\kappa$ (up to factors of several) if a larger halo size is adopted.121212For example, Fig. 10 in Linden et al. (2010) shows isotropically-averaged $\kappa\sim 3\times 10^{28}{\rm cm^{2}/s}$ with $z_{\rm h}\sim 3\;{\rm kpc}$ but $\kappa\sim 10^{29}{\rm cm^{2}/s}$ with a larger $z_{\rm h}\sim 5\;{\rm kpc})$.

There are other substantial uncertainties in the estimates of $\kappa$, as these empirical constraints usually neglect e.g. local variations in $\kappa$ 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 $\kappa$ 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 $\kappa$, which is a common approach in the literature (e.g. Booth et al., 2013; Salem & Bryan, 2014; Pakmor et al., 2016; Pfrommer et al., 2017b; Wiener et al., 2017), and test a wide range of $\kappa$.

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 $\kappa$ everywhere, even at large heights above the disk. It is therefore likely that our simulations will require a larger $\kappa$ 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 $\kappa$ would be factor $\gtrsim 3$ larger.

2.2.5 Sources & Injection

We assume CR injection from SNe (including Type Ia and Type II), with a fixed fraction $\epsilon_{\rm cr}$ ($=0.1$, as our default value) of the initial ejecta energy ($\Delta E_{\rm cr}=\epsilon_{\rm cr}\,E_{\rm SNe}$ with $E_{\rm SNe}\approx 10^{51}\,{\rm 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 Hopkins et al. (2018b); we therefore reduce the coupled energy by $1-\epsilon_{\rm cr}$ and inject the remaining $\epsilon_{\rm 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 (${\bf F}_{\rm cr}\rightarrow{\bf F}_{\rm cr}+\Delta{\bf F}_{\rm cr}$ with $\Delta{\bf F}_{\rm cr}=\Delta e_{\rm cr}\,\tilde{c}\,\hat{\bf r}$ from the source, where $\hat{\bf r}$ is a unit vector pointing outwards from the source). The injection is therefore operator-split and solved discretely (associated with each SNe).

2.2.6 Hadronic & Coulomb Losses (“Cooling”)

We adopt the estimate for combined hadronic ($\tilde{\Lambda}_{\rm cr,had}$) plus Coulomb ($\tilde{\Lambda}_{\rm cr,\,Cou}$) losses, $\Gamma_{\rm cr}$, from Völk et al. (1996) and Ensslin et al. (1997) as synthesized and updated in Guo & Oh (2008):

[TABLE]

where $n_{\rm n}$ is the number density of nucleons and $x_{e}$ is the number of free electrons per nucleon. Following Guo & Oh (2008) we assume $\sim 1/6$ of the hadronic losses and all Coulomb losses are thermalized, adding a volumetric gas heating term

[TABLE]

The remaining CR losses are assumed to escape in the form of $\gamma$-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 ${\bf F}_{\rm cr}\rightarrow{\bf F}_{\rm cr}(1-\tilde{\Lambda}_{\rm cr}\,n_{\rm n}\Delta 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 Hopkins et al. (2018c).

2.2.7 “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 $\hat{\bf B}$ wherever it appears above (representing projection of motion along field lines) with $\hat{\nabla}P_{\rm cr}$.

3 Simulation setup

3.1 Initial conditions

We study the impact of CRs on three characteristic types of galaxies, dwarf (Dwarf), sub-$L\star$ starburst (Starburst) and $L\star$ ($L\star$** Galaxy**) galaxies, whose details are listed in Table 1.

All of the runs have exponential stellar and gas disks with scale radii $d_{\rm*,disk}$ and $d_{\rm g,disk}$ respectively. We also include small stellar bulges with Hernquist profiles (Hernquist, 1990) 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 Bullock et al. (2001), 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\star$ Galaxy to be 0.1, 0.3, and 1.0 $Z_{\odot}$ 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\star$** 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\star$** Galaxy** and 300 Myr in Dwarf. In the runs with magnetic fields we start with a seed magnetic field with $10^{-2}\mu$G uniformly (over all gas particles) pointing along the direction of disk angular momentum. The magnetic fields rapidly amplify to $\sim\mu$G in dense gas and develop toroidal morphologies with significant turbulent structure, by around a hundred Myr (see Su et al. 2017; Su et al. 2018). 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 ($\sim 0.1\;{\rm g/cm^{2}}$) and SFR ($\sim 5\;{\rm M}_{\odot}/{\rm 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 B.8) that global quantities of interest are robust at our default resolution indicated in Table 1, and that main qualitative effects of CRs on galaxies can be captured at this resolution.

3.2 Cosmic Ray Propagation Models

We consider several different CR propagation models, and a range of diffusion coefficients. All models are listed in Table 2. In particular, we consider diffusion coefficients up to ten times higher and lower than the common inferred isotropically-averaged MW values, sampling a range $\kappa=3\times 10^{27}{\rm cm^{2}/s}-3\times 10^{29}{\rm 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 (see e.g. Strong & Moskalenko, 1998; Jubelgas et al., 2008; Lacki et al., 2010).

We apply the newly developed two-moment method (§ 2.2.3) to both streaming and diffusion with a reduced speed of light, $\tilde{c}$. In Appendix B.6 we test different choices for this parameter and demonstrate that physical properties, e.g. SFR or $\gamma$-ray emission, are not affected by the choice of $\tilde{c}$ as long as it is equal to or larger than the values listed in Table 2.

4 Results

4.1 Distribution of cosmic rays and the effects on galactic properties

4.1.1 Dwarf and $L\star$ galaxies

We begin with a brief overview of the evolution of the gas, CRs and magnetic fields in our $L\star$** 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 (Su et al., 2018). 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 5 and the related text). However, there are also other CR energy gain and loss processes, which we will investigate in §  5.1.2. But in all runs, the CR source distributions are much more concentrated than CR densities, as CRs move from their “birthplace”.

Figure 1 shows the distribution of the CR energy density in a 60kpc$\times$60kpc slice centered on our simulated $L\star$** 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 2, 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 ($\kappa=0$ but $v_{\rm st}\neq 0$) are similar to cases with very low diffusivity ($\kappa\lesssim 3\times 10^{27}\,{\rm cm^{2}\,s^{-1}}$), where CRs are largely confined to the galaxy. From Figure 1 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 ($\kappa\gtrsim 10^{27}\,{\rm 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_{\rm 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_{\rm st}\sim 10\,{\rm km\,s^{-1}}(B/\mu{\rm G})(n_{\rm ISM}/0.1\,{\rm cm^{-3}})^{-1/2}$ (see Su et al. 2018), giving a transport time through a gas halo with a radius $\ell$ of $t_{\rm st}\sim 100\,\,{\rm Myr\;}(\ell/{\rm kpc})(v_{\rm st}/10\,{\rm km\,s^{-1}})^{-1}$, while the corresponding diffusive transport velocity/time is $v_{\rm di}\sim 330\,{\rm km\,s^{-1}}(\kappa/10^{29}{\rm cm^{2}\,s^{-1}})(\ell/{\rm kpc})^{-1}$ and $t_{\rm di}\sim 3\,{\rm Myr}\,(\ell/{\rm kpc})^{2}(\kappa/10^{29}{\rm cm^{2}\,s^{-1}})^{-1}$ . Thus even for quite low $\kappa$, the diffusive flux dominates transport on sub-kpc scales.

But because the diffusion time scales as $\sim\ell^{2}$, if the CRs establish a smooth profile with scale length $\gtrsim 1-10\,$kpc, then on the larger scales the diffusion time eventually could become larger than the streaming transport time, i.e. outside a scale $\ell\sim 30\,{\rm kpc}\,(\kappa/10^{29}{\rm cm^{2}\,s^{-1}})(v_{\rm st}/{\rm 10\,km\,s^{-1}})^{-1}$.

We quantify the above observations with Figure 3, which shows the cumulative distribution of CR energy in the Dwarf and $L\star$** 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_{\rm cr,1/2}$ such that the sphere with $r_{\rm cr,1/2}$ encloses one half of the total CR energy. In the $L\star$** Galaxy**, we find $r_{\rm cr,1/2}$ is around 3 kpc in run “${\rm\kappa=3e27}$”, but it increases to around 10 kpc in run “${\rm\kappa=3e28}$” and 30 kpc in run “${\rm\kappa=3e29}$.” Similarly, the scale-height of the CR energy distribution also increases with increasing $\kappa$. Trends of the CR scale radius with $\kappa$ can be understood with a diffusion model where the CR injection time ($E_{\rm cr}/\dot{E}_{\rm inj}$) is comparable to the CR diffusion escape time ($t_{\rm di}\sim r_{\rm cr,1/2}^{2}/\kappa$) where $\dot{E}_{\rm inj}$ is the CR energy injection rate. Assuming a similar injection time, we find that $r_{\rm cr,1/2}$ is roughly proportional to $\sqrt{\kappa}$, so a faster diffusion leads to a more extended CR distribution.

We show the CR and gas mid-plane densities for our $L\star$** Galaxy** and Dwarf runs in Figure 4. For both galaxy types, CR density profiles are significantly “flatter” (more extended and less centrally-concentrated) with higher $\kappa$. 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 “${\rm\kappa=3e27}$” because of the smaller adiabatic energy gain (a highly non-linear effect), which we will discuss in the next section and Figure 13. The gas midplane density depends rather weakly on CRs. Low-$\kappa$ (or streaming/advection-only) runs have slightly higher midplane densities while higher-$\kappa$ 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-$\kappa$ 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\star$ galaxy model with that observed near the solar circle ($e_{\rm cr}\sim 1\;{\rm eV/cm^{3}}$, see e.g. Grenier et al. 2015), but we must recall that the $L\star$ 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 $\sim 8\,$kpc from its center are lower than the $\sim 1\,{\rm cm^{-3}}$ in the solar neighborhood (Moskalenko et al., 2002; Cox, 2005).

Nevertheless, the model has a stellar mass, gas mass, and SFR similar to the MW. Our runs with isotropic $\kappa\leq 3\times 10^{27}\,{\rm cm^{2}\,s^{-1}}$ produce a mid-plane $e_{\rm cr}$ at $\sim 8\,$kpc from the galaxy center which is high relative to the observed value, those with isotropic $\kappa\geq 3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ are lower, while those in-between are reasonably consistent. Turning on MHD (making the diffusion anisotropic) increases $e_{\rm cr}$ by factors $\sim 2-3$, consistent with the isotropically-averaged $\kappa$ being lower by a similar factor (as expected), so values $\kappa\sim 1-3\times 10^{29}\,{\rm 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 $\gamma$-ray emission in § 4.2.

Fig. 5 shows cumulative SF histories: akin to the disk midplane-pressure effects above, CR runs with very low $\kappa$ suppress SF by modest factors $\sim 1.5-2$, an effect which vanishes at higher $\kappa$. Smaller variations ($\sim 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 Su et al. (2017); Su et al. (2018), but the stochasticity during the early evolution can wash out such effect.

4.1.2 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 $\sim 5\,{\rm M}_{\odot}\,{\rm yr^{-1}}$ with highest central gas densities $\sim 100\,{\rm cm^{-3}}$ or (edge-on) surface densities $\sim 0.1\,{\rm g\,cm^{-2}}$ ($\sim 500\,{\rm{\rm M}_{\odot}\,pc^{-2}})$ (compare Fig. 6 and Weiß et al. 2001), which last $\sim 10\,$Myr. Between starbursts the galaxy has lower SFR and gas densities, with correspondingly longer hadronic loss times, so CRs can escape more easily.

Fig. 6 shows that CR energy densities during starburst phases are around 100 ${\rm eV/cm^{3}}$ with slow transport, similar to the value inferred from modeling the observed $\gamma$-ray spectra of e.g. M82 (L11 and Yoast-Hull et al. 2016). 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 ($\pi G\Sigma_{\rm g}\Sigma_{\rm g}\sim 10^{3-4}{\rm eV/cm^{3}}$), which is primarily comparable to the kinetic energy density in these galaxies (with turbulent velocity dispersions similar to those observed, $\sim 50-100\,{\rm km\,s^{-1}}$). Our findings are therefore consistent with earlier claims by Lacki et al. (2010), L11, 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.

4.2 Pionic $\gamma$-ray emission as a measure of CR propagation

Owing to the lack of direct measurements of primary and secondary CRs at low ($\sim$ GeV) energies from extra galactic sources, pionic $\gamma$-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 $\gamma$-rays. While there is a substantial amount of pionic $\gamma$-ray emission with energy < 1 GeV, it is difficult to isolate it observationally owing to contamination by leptonic emission. For $\gamma$-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 Pfrommer et al. 2017b). We will ignore additional potential channels of $>1$ GeV $\gamma$-ray production, e.g. pulsars or dark matter annihilation. Hence, in the following, we assume all $>1$ GeV $\gamma$-rays are pionic (if there is substantial pulsar contamination, the pionic $\gamma$-ray emission is lower, and higher diffusivities $\kappa$ are required).

The $>1$ GeV $\gamma$-ray luminosity for $\gamma$-rays $L_{\gamma}(\geq{\rm GeV})$ can be calculated as:

[TABLE]

where we sum over gas particle i with volume $\Delta V_{i}$. First, the most of the hadronic loss ($\tilde{\Lambda}_{\rm cr,had}\,e_{\rm cr}\,n_{\rm n}$ in Eq. 6) is responsible for the pion production. Second, only one third of the pions ($\pi^{0}$) produce $\gamma$-rays. Third, $\beta_{\pi}(\approx 0.7)$ is the fraction of the pionic $\gamma$-rays with energy above GeV(Lacki et al., 2011), which is calculated with the GALPROP pionic cross sections (Moskalenko & Strong, 1998; Strong & Moskalenko, 1998; Strong et al., 2000) built on Dermer (1986), 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 $\gamma$-ray emission will be relatively weak, compared to the expectations from the CR injection. We follow Thompson et al. (2007) and L11 and quantify this by comparing the pionic $\gamma$-ray luminosity $L_{\gamma}$ (above $>1$ GeV) with the bolometric “star formation luminosity” $L_{\rm 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_{\gamma}/L_{\rm SF}$ is approximately constant.

The value of $L_{\gamma}$ 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:

[TABLE]

where $\epsilon_{\rm cr}(=0.1)$ is the fraction of SNe energy going into CRs (constant by assumption in our simulations), $E_{\rm SN}(=10^{51}\;{\rm erg})$ is the energy from one supernova (also constant by assumption), $\xi(m_{*}>8{\rm M}_{\odot})(=0.0037)$ is fraction of stars that end as supernovae, and $<m_{*}>(=0.4\;{\rm M}_{\odot})$ is the mean stellar mass, both calculated for the same (Kroupa, 2002) 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. $\dot{E}_{\rm cr,SN}=e_{\rm cr}\tilde{\Lambda}_{\rm cr}$, then the pionic $\gamma$-ray luminosity is

[TABLE]

The corresponding “star formation luminosity” $L_{\rm SF}$ for a constant SFR assuming again the same Kroupa (2002) IMF adopted in our simulations is

[TABLE]

where the prefactor is calculated with STARBURST99131313http://www.stsci.edu/science/starburst99/docs/default.htm (Leitherer et al., 1999; Vázquez & Leitherer, 2005; Leitherer et al., 2010; Leitherer et al., 2014).141414L11 adopted a different conversion factor because they assumed a Salpeter IMF (Salpeter, 1955) (following Kennicutt 1998). But the ratio $L_{\gamma}/L_{\rm sf}$ in L11 is only higher by $\sim 10\%$ since the SNe rate is also adjusted accordingly. So in the constant-SFR, calorimetric limit, we would expect $L_{\gamma}/L_{\rm SF}\approx 2\times 10^{-4}$ (Thompson et al., 2007).

4.2.1 Pionic $\gamma$-ray emission

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

Fig. 8 quantifies the distribution of the $\gamma$-ray luminosity within spheres of increasing radii for all of our runs with cosmic rays. Consistent with the discussion above, $\gamma$-ray emission is much weaker for large $\kappa$.151515As with the CR energy density, we note that sometimes the runs with finite-but-low $\kappa$ exhibit slightly higher $L_{\gamma}$ even than the “Advection” runs owing to non-linear effects discussed in § 5.2. For our $L\star$** Galaxy**, half of the $\gamma$-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. 9 breaks down the distribution of CR energy as a function of local gas density, which (since $L_{\gamma}\propto e_{\rm cr}\,n_{\rm gas}$) effectively determines $L_{\gamma}$. With low-$\kappa$ (or with advection/streaming only), CRs reside longer in the high-density regions where they are injected. If the density (on some scale $\ell$ of the cloud or disk) is larger than $\sim 10\,{\rm cm^{-3}}\,(\ell/{\rm kpc})^{-2}\,(\kappa/10^{29}\,{\rm 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_{\gamma}$ is lower at higher $\kappa$, even when the SFR (injection rate) is higher, because the bulk of the CR energy is at lower $n_{\rm n}<0.1\,{\rm cm^{-3}}$.

4.2.2 $L_{\gamma}/L_{\rm SF}$ ratio and scalings

Fig. 10 shows the time evolution of $L_{\gamma}/L_{\rm SF}$. With a lower $\kappa$ (“Advection” or “${\rm\kappa=3e27}$”), the galaxies are closer to the calorimetric limit (as expected).161616Note there are periods where $L_{\gamma}/L_{\rm SF}$ exceeds calorimetric – this is perfectly allowed. Usually it occurs because of short-timescale drops in the SFR and hence $L_{\rm SF}$, where the CRs take somewhat longer to decay so $L_{\gamma}$ is $\sim$ 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 “${\rm\kappa=3e29}$” runs are lower than the calorimetric limit by more than an order of magnitude in our Dwarf and $L\star$ runs.

Turning around our argument from § 4.2.1 above, if CRs are injected in a structure of size $\ell$ and gas density $n$ with an isotropically-averaged diffusivity $\kappa\lesssim 10^{28}\,{\rm cm^{2}\,s^{-1}}\,(\ell/{\rm kpc})^{2}\,(n/{\rm cm^{-3}})$, the collisional loss time becomes shorter than the escape time, so we expect near-calorimetric $L_{\gamma}/L_{\rm SF}$. On the other hand, at larger $\kappa\gg 10^{28}\,{\rm 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_{\gamma}/L_{\rm SF}$) should scale $\propto\kappa^{-1}$.

In § 5.2, we will show that adiabatic processes are of secondary importance relative to diffusion in reducing $L_{\gamma}$ in runs with high-$\kappa$, but in runs with low-$\kappa$ (e.g. “${\rm\kappa=3e27}$”), they boost $L_{\gamma}$ considerably.

As expected, anisotropic diffusion tends to suppress the isotropically-averaged $\kappa$ by factors $\sim 1.5-3$, and correspondingly increase $L_{\gamma}/L_{\rm SF}$. Streaming slightly increases the escape and dissipates CR energy through streaming loss, so decreases $L_{\gamma}/L_{\rm SF}$, but the effect is very small (and streaming alone produces near-calorimetric results). This is because (as discussed in § 4.1.1 above) the streaming escape time is much longer than the diffusive escape time, even for relatively low $\kappa$, 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 (Farber et al., 2018), which could significantly reduce $\gamma$ 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_{\gamma}/L_{\rm SF}$ , which can be seen in Fig. 10.

4.2.3 Comparison to observations

We now compare the simulations to observational estimates of $L_{\gamma}/L_{\rm sf}$ as a function of either central gas surface density in galaxies ($\Sigma_{\rm g}$) or SFR, as compiled in L11. Most of the observed data is described in L11, but we also include the SFR of the SMC (0.036 ${\rm{\rm M}_{\odot}/yr}$ from Wilke et al. 2004). We add two extra starburst galaxies in Figs. 11-12 (NGC 1068 and NGC 4945; which are listed in Table 2 of Lacki et al. 2011 but not in their figures). The SFRs of starburst galaxies (NGC1068, NGC253, NGC4945, M82) are obtained with the Kennicutt (1998) IR to SFR conversion formula, assuming the Kroupa IMF (Their IR luminosities are also listed in Table 2 of Lacki et al. 2011).

$\gamma$-ray observations of nearby galaxies are limited in spatial extent due to energy resolution and contamination from the diffuse backgrounds and foregrounds (Abdo et al., 2010c, a, b). Hence, we only consider $\gamma$-ray emission within 3 kpc for Dwarf, matched approximately to that used for the SMC. This choice reduces the $\gamma$-ray luminosity by a factor of two compared to using an infinitely large aperture. For Starburst and $L\star$** Galaxy**, we take 10 kpc apertures (which only reduces $L_{\gamma}$ by tens of percent compared to an infinitely-large aperture), matched to those used for e.g. M31, NGC1068, and M82. See Fig. 8 for how this scales with size.

We measure the gas surface densities $\Sigma_{\rm g}$ (averaged over viewing angles) of Starburst within 250 pc, $L\star$** 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\star$171717This choice is smaller than the optical radii that Kennicutt (1998) and L11 used, but the gas surface densities are similar in both choices. and about equal in Starburst (similar to the choice in Kennicutt 1998 and L11).

Figs. 11-12 compare our simulations with the observations (compare to Figure 2 of L11). As expected based on our discussion above, $L_{\gamma}/L_{\rm SF}$ is high and close to the calorimetric limit for Dwarf and $L\star$** Galaxy** with slow CR transport, i.e. for “Advection”, “${\rm\kappa=3e27}$”, and “MHD Streaming”. These values are clearly well above the observationally inferred $L_{\gamma}/L_{\rm SF}$.

With larger diffusion coefficients, $L_{\gamma}/L_{\rm SF}$ decreases as expected. For isotropic diffusion, the observations in dwarf and $L\star$ galaxies appear to require $\kappa_{\rm isotropic}\sim 10^{29}\,{\rm cm^{2}\,s^{-1}}$. Out of the options tested $\kappa=3\times 10^{29}\,{\rm 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_{\gamma}$ is somewhat larger owing to suppressed isotropic-averaged diffusivity, as discussed above, so values of the parallel diffusivity $\kappa_{\|}>10^{29}\,{\rm cm^{2}\,s^{-1}}$ are favored.

For galaxies with high gas surface densities and SFRs, i.e. Starburst, we found $\kappa$ has to be less than $3\times 10^{29}{\rm 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. Ackermann et al. (2012b) and Rojas-Bravo & Araya (2016), considered $L_{\gamma}$ with $0.1{\rm GeV}<E_{\gamma}<100{\rm GeV}$, instead of $E_{\gamma}>1{\rm GeV}$. Therefore, in Appendix C, we compare $L_{\rm 0.1-100\;GeV}$ from simulations with observations, and we find that the same high diffusion simulations provide the best match to the observed $\gamma$ ray emission from galaxies.

It is interesting to compare the results from the simple leaky-box model of L11, as well as more detailed models of CR transport in the MW, with our findings. L11 predicts $L_{\gamma}/L_{\rm 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 Lacki et al. 2010 for details), with an isotropically-averaged $\kappa=3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$. The broad contours of their prediction for $L_{\gamma}/L_{\rm SF}$ as a function of $\Sigma_{\rm g}$ or SFR are similar to our simulations, suggesting – as they argued – that CR escape is required to reproduce the observed trend of $L_{\gamma}/L_{\rm SF}$.

In the MW, much more detailed propagation models have been tested (see e.g. Trotta et al. 2011 and reference in § 1). We again caution that our “$L\star$” model is not an exact MW analogue, since it has higher gas surface density and lacks a central gas deficiency like the MW (Moskalenko et al., 2002), both of which could affect $L_{\gamma}$.

Note that at “face value”, both MW and L11 constraints might appear to favor slightly-lower $\kappa\sim 3-6\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$ compared to the best-match here ($\kappa\sim 3\times 10^{29}\,{\rm 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 L11 models assume relatively small halos out of which the CRs escape instantly, while we assume a constant $\kappa$ everywhere, meaning that our effective halo size is large ($\sim 10-30\,$kpc). Recall (§ 1 and § 2.2.4), the inferred $\kappa$ in the observations increases with the halo height. (2) The gas in the simulations is clumpy where CRs are injected, slightly increasing $L_{\gamma}$ (Boettcher et al., 2013), compared to the smooth mass profiles assumed in those studies (requiring larger $\kappa$ by a factor $\sim 1.5-2$). (3) L11 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 $\kappa$ (what is nominally constrained by the L11 study, for example) is a factor $\sim 2-3$ lower than the parallel $\kappa$. 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. Trotta et al. 2011, and references in § 1.

5 Discussion

5.1 Comparisons to previous studies

5.1.1 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 $\kappa\lesssim 10^{28}\,{\rm cm^{2}\,s^{-1}}$, consistent with many previous findings, e.g. Booth et al. (2013); Salem & Bryan (2014); Pfrommer et al. (2017a). However, such slow transport severely violates constraints from observed $\gamma$-ray emission, and at best results in modest SFR suppression (factor $\sim 1.5-2$). For larger transport speeds required to reproduce the observed $\gamma$-ray emission, CRs have only a weak effect on SF.

Interestingly, Jubelgas et al. (2008) 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 ($\propto{\rm SFR}\propto\rho^{1.5}$ in their model) is balanced by collisional losses ($\propto\rho$) locally (like in the calorimetric limit), in an isothermal-like ISM, so in their models the CR energy density is proportional to $\rho^{1/2}$ while thermal energy densities are proportional to $\rho$: as a result, CR energy was always sub-dominant to thermal energy in their models at gas densities $n>0.2\,{\rm 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 $\kappa$ can be significant even at moderate ISM densities, providing mild suppression of the star formation (similar arguments were presented in Socrates et al. 2008; Booth et al. 2013). For our favorite, large $\kappa$ values CRs escape from the ISM, resulting in practically no effect on the star formation in both our Dwarf and $L\star$** 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 $\kappa$, 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(Hopkins et al., 2019).

5.1.2 $\gamma$-ray emission

Our results are in line with L11: when matching the observed $\gamma$-ray emission, starburst galaxies (with effective isotropic diffusivities $\kappa<3\times 10^{29}{\rm cm^{2}/s}$) are nearly proton calorimeters, while galaxies with lower gas surface density or SFRs (with $\kappa\sim 3\times 10^{28-29}{\rm cm^{2}/s}$) are not proton calorimeters (most CR protons escape).

Salem et al. (2016) also studied hadronic $\gamma$-ray emission with simulations of MW-mass galaxies and argued for isotropically-averaged coefficients $\kappa\sim 3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$; but they only considered the $\gamma$-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 $\gamma$-ray luminosity significantly exceeded the CR injection rate at lower $\kappa$. Moreover, as noted by Jacob et al. (2018), 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 $\kappa$.

Recently, Pfrommer et al. (2017b) also investigated $\gamma$-ray emission with idealized galaxy simulations, assuming CR transport via either advection-only or advection+anisotropic diffusion with $\kappa=10^{28}\,{\rm cm^{2}\,s^{-1}}$. They argued they could (a) reproduce the FIR-$\gamma$-ray correlation and (b) explain the low $L_{\gamma}$ in non-starburst galaxies primarily by adiabatic losses.

But there are several caveats:

(1) Their favored model still over-predicted $L_{\gamma}/L_{\rm 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_{\gamma}/L_{\rm SF}$ is larger than the SMC, LMC, MW, and M33 (not shown therein, but see L11).

They claimed to match the observed FIR-$L_{\gamma}$ correlation, only if an empirical FIR-SFR conversion relation (Kennicutt, 1998) is assumed. However, as they acknowledged, this conversion relation over-predicts $L_{\rm FIR}$ in dwarfs, due to much lower dust opacity/absorption/reddening. Their FIR-$L_{\gamma}$ relation might deviate from observations after taking this correction into account.181818Because both the $L_{\gamma}/L_{\rm SF}$ and $L_{\rm FIR}/{\rm SFR}$ ratios drop in dwarfs, a roughly-linear $L_{\rm FIR}-L_{\gamma}$ relation can still maintain in our simulations (directly related to the “conspiracy” which maintains the FIR-radio correlation; for discussion see e.g. Bell 2003; Lacki et al. 2010).

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

They suggested their over-prediction of $L_{\gamma}/L_{\rm 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-$\kappa$ runs. A possible explanation for the discrepancy is that the observed low $\gamma$-ray luminosities require high diffusion coefficients $\kappa\sim 3\times 10^{28-29}\,{\rm cm^{2}\,s^{-1}}$ as favored by our study here and the modern MW constraints (Trotta et al., 2011).

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

5.2 CR Energetics and the Importance of Different Gain/Loss Terms

Fig. 13 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 $\sim 2$).

Collisional losses are important loss terms (within the galaxies) – and we have already discussed these extensively as they are the origin of the $\gamma$-ray emission. Since they scale $\propto e_{\rm cr}\,n_{\rm n}$ they decrease with “faster” CR transport (higher $\kappa$) as CRs reach lower-density gas faster.191919The (weak) exception to this rule is the $\kappa=3\times 10^{27}\,{\rm cm^{2}\,s^{-1}}$ run in our $L\star$ 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_{\rm cr}$ in dense gas. Streaming losses are comparatively small.

The “adiabatic” term $\dot{E}_{\rm 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. 14 shows the contribution to the total $\dot{E}_{\rm Ad}$ from gas with different densities $n$ or at different galacto-centric radii $n$. For CRs at low ambient $n$ (or large $r$), $\dot{E}_{\rm Ad}$ tends to be a loss term (i.e. CRs are expanding or losing energy in rarefactions). For CRs in high ambient $n_{\rm ISM}$ and small $r$, it tends to be a gain (CRs are being compressed in converging flows). Recall, the “adiabatic” term is defined (Eq. 2) by $-\int d^{3}{\bf x}P_{\rm cr}\,(\nabla\cdot{\bf v})=-(\gamma_{\rm cr}-1)\,\int dE_{\rm cr}\,\nabla\cdot{\bf v}$ (where ${\bf v}$ is the gas velocity). So, combined with Fig. 14, 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. 9).

So, unsurprisingly in Fig. 13, our runs with the most efficient CR escape to large-$r$ and low-$n$ (all the highest $\kappa\sim 3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ runs, or most of the Dwarf runs with even intermediate $\kappa$) show net $\dot{E}_{\rm 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\star$ runs with lower isotropically-averaged $\kappa$) show net $\dot{E}_{\rm Ad}>0$.202020Interestingly, if CRs do not preferentially stay in regions where the adiabatic term is mostly negative or positive, then the $\dot{E}_{\rm 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-$\kappa$ ($\sim 10^{28}\,{\rm cm^{2}\,s^{-1}}$), the qualitative behaviors of $\dot{E}_{\rm Ad}$ in both dwarf and MW-mass systems in Pfrommer et al. (2017b) are similar to what we find here.

However, for the reasons discussed in § 5.1.2, our results do not support their conclusion that adiabatic losses are the dominant factor for the low $L_{\gamma}/L_{\rm SF}$ in dwarfs. At very low $\kappa$, $\dot{E}_{\rm Ad}$ is primarily a gain term. But even at higher $\kappa$ where $\dot{E}_{\rm Ad}<0$ is a loss term, it is insufficient (in itself) to explain the very low $L_{\gamma}/L_{\rm 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. 13). 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_{\gamma}$.

To explicitly show that CR leakage is significant, in Fig. 15 we plot the fraction of CR energy that escapes the central galactic region, $E_{\rm esc}/E_{\rm source}$. We define $E_{\rm source}$ as the total cumulative CR energy input within the central 6 kpc (10 kpc) and $E_{\rm esc}$ as the CR energy that leaves this central region in Dwarf ($L\star$ galaxy)212121In practice, we calculate $E_{\rm 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 8. To calculate $E_{\rm 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 $\kappa$), most of the CR energy indeed escapes the central region, where most $\gamma$ rays are produced. This shows that CR escape is the main reason for reduced $\gamma$ ray emission in low-gas-density galaxies.

6 Conclusions

We explore the effects of CRs on galaxies, in high-resolution, idealized (non-cosmological) (magneto-)hydrodynamic simulations of dwarf, $L\star$, and sub-$L\star$ 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 $\gamma$-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 $\sim 3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ and $\sim$ pc resolution. 2. 2.

The CR “transport parameters”, in particular, the effective diffusivity $\kappa\equiv|{\bf F}_{\rm cr}|/|\nabla_{\|}e_{\rm cr}|$ (which can, microphysically, arise from a combination of streaming and diffusion), have a significant impact on galaxy properties and predicted $\gamma$-ray emission. With very slow propagation ($\kappa\lesssim 10^{28}{\rm 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 $\sim 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 $\gamma$-ray luminosities than the observed for dwarf or MW-like systems. At higher $\kappa\gtrsim 10^{28}\,{\rm 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. 3. 3.

The extent of the CR halo, and correspondingly the extent of the pionic $\gamma$-ray emission, increase with $\kappa$ as expected. For e.g. our $L\star$ galaxy, half the CR energy is located within 10 kpc (30 kpc) for $\kappa=3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$ ($\kappa=3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$). Correspondingly only about $\sim 50\%$ of the $\gamma$-rays are emitted from the central few kpc. 4. 4.

In our sub-$L\star$ starburst galaxies, the CR energy density reaches $\sim 10^{2-3}\,{\rm eV\,cm^{-3}}$ throughout the burst and is larger than thermal or magnetic pressure in the ISM (for any $\kappa$), 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 L11. 5. 5.

We constrain the average CR propagation speed/diffusivity with $\gamma$-ray ($>$ GeV) emission from galaxies. The observed $L_{\gamma}-{\rm SFR}$ relation requires isotropically-averaged diffusivities $\kappa\sim 3\times 10^{29}{\rm cm^{2}/s}$ in dwarf and $L\star$ galaxies, and $\kappa\lesssim 3\times 10^{29}{\rm cm^{2}/s}$ in sub-$L\star$ 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 $\kappa<3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$,then CRs escape galaxies too slowly and produce $\gamma$-ray luminosity close to the calorimetric limit. This over-predicts the observed $\gamma$-ray luminosities in dwarfs (e.g. the SMC, LMC, M33) and $L\star$ systems (M31, the MW) by an order of magnitude or more. However, for faster transport parameters (effective $\kappa\sim 3\times 10^{28-29}\,{\rm cm^{2}\,s^{-1}}$), CRs escape the dense regions rapidly and the $\gamma$-ray luminosity (which scales $\propto e_{\rm cr}\,n_{\rm gas}$) is reduced (especially in dwarf galaxies), predicting $\gamma$-ray luminosities in good agreement with those observed as a function of either gas surface density or SFR (see Figs. 11 and 12). 6. 6.

Given the transport parameters required to reproduce the observed $L_{\gamma}/L_{\rm 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\star$) starburst models are (approximate) proton calorimeters. 7. 7.

CR streaming at trans-Alfvenic speeds is relatively slow and cannot alone reduce $L_{\gamma}/L_{\rm SF}$ significantly below the calorimetric limit in our models (as required by observations), even if we allow modestly super-Alfvenic streaming (with $\sim 4\,v_{\rm A}$; see Appendix A ). For our favored effective $\kappa$, the equivalent streaming speed (using the fact that the CR flux is similar for a diffusivity $\kappa$ or streaming speed $v_{\rm st}\sim\kappa\,|\nabla P_{\rm cr}|/P_{\rm cr}$) is $\sim 10-100$ times the Alfven speed. 8. 8.

“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_{\gamma}/L_{\rm 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_{\gamma}$), 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 $\gamma$-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 (Booth et al., 2013; Salem et al., 2016). 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 (Muratov et al., 2015; Anglés-Alcázar et al., 2017; Hafen et al., 2019) 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 (Condon, 1992). 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 (van der Kruit, 1971, 1973) with galaxy simulations in a manner similar to our analysis of the connection between SFR and $\gamma$-ray emission. However, as mentioned in Lacki et al. (2010), 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 $\kappa$. 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. Jokipii 1966; Enßlin 2003). It would be interesting to investigate galaxy evolution and CR observables in studies where the CR transport coefficients vary dynamically and locally (see e.g. Farber et al., 2018), or with recently-developed models which attempt to actually predict the coefficients self-consistently (Thomas & Pfrommer, 2019).

7 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 (Turk et al., 2011), matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), scipy (Jones et al., 2001), and NASA’s Astrophysics Data System.

Appendix A A comparison of different streaming parameters

In § 2.2.2, we discussed uncertainties in the parameters describing CR streaming: both the streaming speed $v_{\rm st}$ and magnitude of the streaming loss term $\Gamma_{\rm st}=-{\bf v}_{\rm A}\cdot\nabla P_{\rm cr}$ owing to excited and thermalized high-frequency Alfven waves (independent of streaming speed; see §2.2.2). Here we explore these more thoroughly. We consider four model variations:

“$v_{\rm st}=v_{\rm A}$”: This is our default choice in the main text, with streaming speed equal to the Alfven speed. 2. 2.

“$v_{\rm st}=3v_{\rm A}$”: Here $v_{\rm st}=3\,v_{\rm A}$, one of the super Alfvenic speeds considered in Ruszkowski et al. (2017). 3. 3.

“$v_{\rm st}=4v_{\rm A}$”: Here $v_{\rm st}=4v_{\rm A}$, another super Alfvenic speed considered in Ruszkowski et al. (2017) 4. 4.

“$v_{\rm st}=v_{\rm A+c_{\rm s}}$”: Here $v_{\rm st}=\sqrt{v_{\rm A}^{2}+c_{\rm 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 $\beta\sim 1$ is simply assumed).

Fig. 16 summarizes the results: the effects of this choice are much smaller than the variations of e.g. $\kappa$ discussed in the main text. Most importantly, we see no difference in averaged $L_{\gamma}/L_{\rm 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 ${\bf v}_{\rm A}\cdot\nabla P_{\rm cr}$ and continue to increase $v_{\rm st}$ without limit, it will eventually become “fast enough.” In fact, our Eq. 4 shows that, given the manner in which we approximate CR transport numerically, $\kappa$ and $v_{\rm st}$ are formally degenerate if we replace $\kappa$ with $v_{\rm st}\rightarrow\kappa\,|\hat{\bf B}\cdot\nabla e_{\rm cr}|/(e_{\rm cr}+P_{\rm cr})$. From this, using a typical CR gradient scale length $\sim 1-10\,$kpc in our galaxies and favored $\kappa\sim 3\times 10^{29}\,{\rm cm^{-2}\,s^{-1}}$, we see that $v_{\rm st}\sim 100-1000\,{\rm km\,s^{-1}}\gg v_{\rm A}$ (and $\gg c_{s}$ in most of the ISM) is required.

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

Appendix B Numerical Tests

Our CR transport implementation is described in § 2.2.3. 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” $\tilde{c}$), comparison of our two-moment implementation to zeroth-moment “pure diffusion” solvers, and numerical resolution studies.

B.1 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 Sod (1978) shocktube presented in Pfrommer et al. (2006). A 3D box of of dimensions $64\times 16\times 16$ is full of gas (adiabatic index $7/5$) and CRs. Half the box has initial $(\rho,\,v_{x}/c_{s},\,P_{\rm gas},\,P_{\rm 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. 17 compares with analytic solutions from Pfrommer et al. (2006) 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 Hopkins 2015a).

B.2 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 ($\tilde{\bf F}_{\rm cr}\neq 0$). The corresponding “two moment” equation is

[TABLE]

which is a telegraph equation.

We initialize a 1D Gaussian distribution in $e_{\rm cr}^{0}\equiv e_{\rm cr}(r,\,t=0)$, centered at $r=0$, with total energy $E_{\rm cr}=1$, width $\sigma=0.5$ kpc at the center of a $5\,$-kpc cube with $2048$ resolution elements, set constant $\kappa=3\times 10^{27}\,{\rm cm^{2}\,s^{-1}}$, and set $\tilde{c}=100\,{\rm km\,s^{-1}}$ or $500\,{\rm km\,s^{-1}}$, and evolve the system for $5\,$Myr. 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:

[TABLE]

where $\kappa\omega_{\bf k}^{2}/(\gamma_{\rm cr}-1)/\tilde{c}^{2}+i\,\omega_{\bf k}=\kappa\,k^{2}$ and $\omega^{\pm}_{\bf k}$ are two roots of the previous equation,

[TABLE]

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

[TABLE]

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

This problem is entirely scale-free, and we can transform to solutions with any other value of $\kappa$ via suitable re-scaling. As $\tilde{c}^{2}\,t/\kappa\rightarrow\infty$, the solutions progressively approach the solution of the pure diffusion equation222222Specifically, at this limit, the $a^{-}_{\bf k}$ term becomes exponentially small and the $a^{+}_{\bf 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:

[TABLE]

Fig. 18 shows the results of our simulation for varying $\tilde{c}$ at fixed $\kappa$ and $t$ (i.e. varying the dimensionless parameters $\tilde{c}^{2}\,t/\kappa$ and $r\,\tilde{c}/\kappa$ which determine the behavior of the problem)232323We 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. 18 also gives us a practical estimate – for typical units and spatial scales of our simulations – of the rate at which solutions with lower $\tilde{c}$ converge to the solution with $\tilde{c}=c$ (the speed of light). In Fig. 18, the test with $\tilde{c}=500\,{\rm km\,s^{-1}}$ is already effectively indistinguishable from the “pure diffusion” solution. Even at $\tilde{c}=100\,{\rm 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.

B.3 Comparison of Zeroth and Two Moment Approximations in Galaxy Simulations

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

B.4 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 $\kappa=3\times 10^{28}{\rm cm^{2}/s}$ and a fixed B field pointing in the x direction. We set $\tilde{c}=1000\,{\rm km/s}$, initial ${\bf F}_{\rm cr}=0$, and evolve over 5 Myr with the zeroth moment (“pure diffusion”) and two moment schemes. Fig. 20 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 ${\bf F}_{\rm cr}$ perpendicular to the B field; see §B.5). We have also checked that the total CR energy is conserved (with deviations smaller than one in a thousand) in both schemes.

B.5 Test of a variation in the CR flux equation (Eq. 3)

Eq. 3 governs the CR flux evolution, and obviously approaches the desired anisotropic diffusion+streaming equation (with $\tilde{\bf F}_{\rm cr}\rightarrow-\kappa_{\ast}\,\nabla_{\|}e_{\rm cr}$) when $\tilde{c}\rightarrow\infty$ is large, or in local quasi-steady-state ($\partial\tilde{\bf F}_{\rm cr}/\partial t\rightarrow 0$), or on spatial scales $>\kappa/\tilde{c}$ (timescales $>\kappa/\tilde{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 direction242424We thank the anonymous referee pointing out this potential issue.. If we divide the flux into parallel and perpendicular components $\tilde{\bf F}_{\rm cr}=\tilde{\bf F}_{\rm cr,\|}+\tilde{\bf F}_{\rm cr,\perp}$, then we see the perpendicular component obeys $\tilde{c}^{-2}\,\partial\tilde{\bf F}_{\rm cr,\perp}/\partial t=-[(\gamma_{\rm cr}-1)/\kappa_{\ast}]\,\tilde{\bf F}_{\rm cr,\perp}$ (taking ${\bf v}=0$, for simplicity, though this does not change our conclusions below). So an initial $\tilde{\bf F}_{\rm cr,\perp}$ can formally propagate but will be exponentially damped on a timescale $\sim\kappa_{\ast}/\tilde{c}^{2}$. For $\kappa^{*}$ needed to match the $\gamma$ ray luminosity of nearby galaxies ($\sim 10^{29}$) and $\tilde{c}\sim 1000\,{\rm km/s}$, this timescale is $\sim 0.2$ Myr, much shorter than the hadronic interaction time in the MW’s ISM ($n_{\rm ISM}\sim 1\,{\rm cm^{-3}}$), so our $\gamma$ ray constraint will not be affected with or without $\tilde{\bf F}_{\rm cr,\perp}$ (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 $\tilde{\bf F}_{\rm cr,\perp}=0$ every timestep, and evolve only the parallel component of Eq. 3:

[TABLE]

where ${\bf X}_{\|}\equiv(\hat{\bf B}\otimes\hat{\bf B})\cdot{\bf X}$. We have implemented Eq. 17 and compare it (in Fig. 21) to our default Eq. 3 in an anisotropic diffusion test (as § B.4) with an initial (intentionally erroneous) pure-perpendicular flux $\tilde{\bf F}_{\rm cr}=\tilde{\bf F}_{\rm cr,\perp}=(2000\,{\rm km/s})e_{\rm cr}\hat{y}$. This shows that the “original” Eq. 3 does formally allow CRs to propagate perpendicular to the B field, although the effect is rapidly damped, while the “modified” Eq. 17 does not.

Fig. 22 compares two galaxy simulations (Dwarf and $L\star$** galaxy**, with the latter at one level lower resolution than the fiducial main-text case) run with MHD (anisotropic diffusion), $\kappa_{\|}$=3e28, and either Eq. 3 or Eq. 17. 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. 3 and Eq. 17 appear on scales $\ll\kappa/\tilde{c}$; but § B.3 showed that taking $\tilde{c}\rightarrow\infty$ does not change our conclusions.

Of course, CRs can physically propagate across mean magnetic fields due to unresolved (micro-scale) magnetic turbulence (Zweibel, 2017; Farber et al., 2018), so a small $\tilde{\bf F}_{\rm cr,\perp}$ 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 $\kappa$.

It is worth noting that several recent studies including this one and Zweibel (2017), Jiang & Oh (2018), and Thomas & Pfrommer (2019) 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 ($\kappa_{\ast}$), all of these can be written (after some algebra) in the form:

[TABLE]

where the differences are contained in the operator $\mathbb{D}_{t}$. For our default Eq. 3, $\mathbb{D}_{t}{\tilde{\mathbf{F}}_{\rm cr}}=\partial{\tilde{\mathbf{F}}_{\rm cr}}/\partial t+\nabla\cdot({\bf v}\otimes{\tilde{\mathbf{F}}_{\rm cr}})$. For the zeroth-moment or “pure diffusion+streaming” equation in § B.3, $\mathbb{D}_{t}{\tilde{\mathbf{F}}_{\rm cr}}=\mathbf{0}$. For Jiang & Oh (2018), $\mathbb{D}_{t}{\tilde{\mathbf{F}}_{\rm cr}}=[\partial\{{\tilde{\mathbf{F}}_{\rm cr}}+{\bf v}\,(e_{\rm cr}+P_{\rm cr})\}/\partial t]_{\|}$, and for Thomas & Pfrommer (2019), $\mathbb{D}_{t}\tilde{\mathbf{F}}_{\rm cr}=\hat{\tilde{\mathbf{F}}}_{\rm cr}\,[\partial|{\tilde{\mathbf{F}}_{\rm cr}}|/\partial t+\nabla\cdot({\bf v}\,|{\tilde{\mathbf{F}}_{\rm cr}}|)+{\tilde{\mathbf{F}}_{\rm cr}}\cdot\{(\hat{\tilde{\mathbf{F}}}_{\rm cr}\cdot\nabla){\bf v}\}]=\partial{\tilde{\mathbf{F}}_{\rm cr}}/\partial t+\nabla\cdot({\bf v}\otimes\tilde{\mathbf{F}}_{\rm cr})+({\tilde{\mathbf{F}}_{\rm cr}}\cdot\nabla)\,({\bf v}_{\|}-{\bf v}_{\perp})$, with $\tilde{\bf F}_{\rm cr,\perp}=\mathbf{0}.$252525For Jiang & Oh (2018), we have formally taken the limit of the diffusion/interaction tensor $\mathbb{K}\equiv\kappa_{\ast}\,(\hat{\bf B}\otimes\hat{\bf B})+\kappa_{\perp}\,(\mathbb{I}-\hat{\bf B}\otimes\hat{\bf B})$ ($\mathbb{I}$ is the identity tensor) as $\kappa_{\perp}\rightarrow 0$. For Thomas & Pfrommer (2019), we use $\tilde{\mathbf{F}}_{\rm cr}=|\tilde{\mathbf{F}}_{\rm cr}|\,\hat{\tilde{\mathbf{F}}}_{\rm cr}=f_{\rm cr}\,\hat{\bf B}$ in their notation, and replace their scattering term $\partial f/\partial t|_{\rm scatt}$ (which explicitly attempts to account for dynamically-evolving, forward-and-backward propagating gyro-scale Alfven waves) with $-[(\gamma_{\rm cr}-1)/\kappa_{\ast}]\tilde{\bf F}_{\rm cr}$ for a locally-stationary $\kappa_{\ast}$.

What is important for our purposes here is to note that all of these expressions differ only inside the term suppressed by $\tilde{c}^{-2}$. This means they all relax to the same diffusion+streaming equation in quasi-steady state and/or as $\tilde{c}\rightarrow\infty$, and differ only on scales below the CR mean free path $\sim\kappa/\tilde{c}$. To the extent that our galaxy-scale results are converged with respect to the value of $\tilde{c}$, and are not changed if we take $\tilde{c}\rightarrow\infty$ (§ B.3), these differences in $\mathbb{D}_{t}\tilde{\bf F}_{\rm cr}$ cannot alter our conclusions. Moreover, if we adopt a “reduced speed of light” $\tilde{c}\ll c$, then by definition the flux equation does not exactly represent reality on scales $\ll\kappa/\tilde{c}$, regardless of the form of $\mathbb{D}_{t}$ (i.e. the regime where these expressions differ is exactly where they all become inexact when $\tilde{c}\ll 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 $\mathbb{D}_{t}\tilde{\bf F}_{\rm cr}$ from Jiang & Oh (2018) or Thomas & Pfrommer (2019) as defined above: the results are similar to those in Fig. 22.

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

We now examine the effect of varying $\tilde{c}$ on a full simulation (our $L\star$ model). For $\tilde{c}\gtrsim 500\,{\rm 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. 23 shows there is a small impact of $\tilde{c}$ on the SFR (we find the same for all other galaxy properties, not shown here). For $L_{\gamma}/L_{\rm SF}$, we find almost no impact of $\tilde{c}$ in simulations with $\kappa<10^{29}{\rm cm^{2}/s}$ for any values $\tilde{c}\gtrsim 500\,{\rm km\,s^{-1}}$.

Because the effective “advective velocity” of CRs under pure diffusion is $\sim\kappa/\ell$ (where $\ell$ is some gradient scale-length), at much larger diffusion coefficients, e.g. our $\kappa=3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$, where most of the CRs escape diffusively, the value of $L_{\gamma}$ is slightly larger for $\tilde{c}=500\,{\rm km\,s^{-1}}$ compared to much-higher $\tilde{c}$ (because escape is slightly slower). But once $\tilde{c}\geq 1000\,{\rm km\,s^{-1}}$, we see no detectable difference. Moreover by the latter half of the time we run for, the differences even for $\tilde{c}=500\,{\rm km\,s^{-1}}$ (compared to $\tilde{c}\sim 4000\,{\rm km\,s^{-1}}$) at $\kappa=3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ are factors of $\sim 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 $\tilde{c}$ within the range $\tilde{c}\sim 500-4000\,{\rm km\,s^{-1}}$, even at $\kappa=3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$.

B.7 Idealized Streaming Test

Here we test our streaming implementation with an initial 1D symmetric triangular CR profile with $e_{\rm cr}=2-|x|$ (and ${\bf F}_{\rm cr}={\bf v}_{\rm str}e_{\rm cr}$), where we disable the gas motion and CR streaming loss ($v_{\rm st}\cdot\nabla P_{\rm cr}$), but CRs can move across gas particles with a constant $v_{\rm st}=1\;{\rm km/s}$. In the numerical test, we consider evenly spaced 2048 particles over 10 kpc, the reduced speed of light $\tilde{c}=1000\,\rm{km/s}$, and a small diffusion coefficient $3\times 10^{22}\;{\rm cm^{2}/s}$.262626We include a very small diffusion coefficient to avoid potential (numerical) overflows in Eq. 3, i.e. when the CR pressure profile is flat or, equivalently, $\nabla P_{\rm cr}$ is huge. Here, we also consider a variant test in which we relax the limiter on CR pressure gradient272727The 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).282828Because 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 Jiang & Oh (2018), has a flat CR distribution between $x=\pm x_{\rm m}$ and two inclined distributions for $|x|>x_{\rm m}\;\&\;|x|<1$ with

[TABLE]

and $x_{\rm m}$ is determined from the energy conservation:

[TABLE]

Fig. 24 shows the analytic solution agrees well with our numerical solution at $t=10\;{\rm Myr}$, although the “original” run has a round top rather than a plateau mainly because of the limiter on the CR pressure gradient.

B.8 Resolution study

In Fig. 25, we show the properties of $L\star$** Galaxy** runs, at three different resolution levels. Our baseline is the fiducial resolution listed in Table 1, but we compare runs with 10x and 100x poorer mass resolution. Even at 10x poorer resolution, we find very similar SFRs, $L_{\gamma}/L_{\rm 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.

Appendix C Comparison with $L_{\rm 0.1-100GeV}$ observations

For completeness, we also compare our simulations with the broader, 0.1-100 GeV, observed energy range of $\gamma$ rays from Ackermann et al. (2012b) and Rojas-Bravo & Araya (2016).

In the main text, we estimate the $L_{\rm>1GeV}$ following the L11 approach with the hadronic loss rate taken from Guo & Oh (2008) (see Eq. 8 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 $\gamma$ rays and enable direct comparison to both the observations and simulations from Pfrommer et al. (2017b), we follow Pfrommer et al. (2017b) to calculate $L_{\rm 0.1-100GeV}$.

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_{p}c)=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 $\gamma$-ray emissivity from pion decay for the energy band 0.1-100 GeV with Eq. 6 in Pfrommer et al. (2017b), and then integrate to get $L_{\rm 0.1-100GeV}$. While the secondary IC emission from high energy CR electrons also contributes to $\gamma$-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 $\gamma$-ray luminosity. Thus, we simply neglect its contribution in Fig.26.

Fig. 26 shows that $\kappa\sim 1{\rm e}29-3{\rm e}29\;{\rm cm^{2}/s}$ provides the best match to the observations for Dwarf and L$\star$ Galaxy runs, consistent with our conclusions in the main text.