But What About... Cosmic Rays, Magnetic Fields, Conduction, & Viscosity in Galaxy Formation

TL;DR

This study uses high-resolution cosmological simulations to investigate the roles of magnetic fields, conduction, viscosity, and cosmic rays in galaxy formation, finding cosmic rays can significantly influence star formation in larger galaxies at lower redshifts.

Contribution

It provides a comprehensive analysis of cosmic ray effects across different galaxy masses and redshifts, incorporating explicit magnetic and plasma physics in simulations.

Findings

Magnetic fields, conduction, and viscosity have minimal effects on bulk galaxy properties.

Cosmic rays can suppress star formation by factors of 2-4 in massive galaxies at low redshift.

Cosmic rays build up in halos, supporting cool gas and affecting galaxy evolution.

Abstract

We present a suite of high-resolution cosmological simulations, using the FIRE-2 feedback physics together with explicit treatment of magnetic fields, anisotropic conduction and viscosity, and cosmic rays (CRs) injected by supernovae (including anisotropic diffusion, streaming, adiabatic, hadronic and Coulomb losses). We survey systems from ultra-faint dwarf ($M_{\ast}\sim 10^{4}\,M_{\odot}$, $M_{\rm halo}\sim 10^{9}\,M_{\odot}$) through Milky Way masses, systematically vary CR parameters (e.g. the diffusion coefficient $\kappa$ and streaming velocity), and study an ensemble of galaxy properties (masses, star formation histories, mass profiles, phase structure, morphologies). We confirm previous conclusions that magnetic fields, conduction, and viscosity on resolved ($\gtrsim 1\,$pc) scales have small effects on bulk galaxy properties. CRs have relatively weak effects on all galaxy properties studied in dwarfs ($M_{\ast} \ll 10^{10}\,M_{\odot}$, $M_{\rm halo} \lesssim 10^{11}\,M_{\odot}$), or at high redshifts ($z\gtrsim 1-2$), for any physically-reasonable parameters. However at higher masses ($M_{\rm halo} \gtrsim 10^{11}\,M_{\odot}$) and $z\lesssim 1-2$, CRs can suppress star formation by factors $\sim 2-4$, given relatively high effective diffusion coefficients $\kappa \gtrsim 3\times10^{29}\,{\rm cm^{2}\,s^{-1}}$. At lower $\kappa$, CRs take too long to escape dense star-forming gas and lose energy to hadronic collisions, producing negligible effects on galaxies and violating empirical constraints from $\gamma$-ray emission. But around $\kappa\sim 3\times10^{29}\,{\rm cm^{2}\,s^{-1}}$, CRs escape the galaxy and build up a CR-pressure-dominated halo which supports dense, cool ($T\ll 10^{6}$ K) gas that would otherwise rain onto the galaxy. CR heating (from collisional and streaming losses) is never dominant.

Tables

Table 1. Table 1: Zoom-in simulation volumes run to z = 0 𝑧 0 z=0 (see Hopkins et al. 2018b for details). All units are physical.

Halo/stellar properties listed refer only to the original “target” halo around which the high-resolution volume is centered: these volumes can reach up to $\sim {\left(1-10\mathrm{Mpc}\right)}^3$ comoving, so there are actually several hundred resolved galaxies in total. (1) Simulation Name: Designation used throughout this paper. (2) $M_{\mathrm{halo}}^{\mathrm{vir}}$: Virial mass (following Bryan & Norman, 1998) of the “target” halo at $z=0$. (3) $M_{\ast }^{\mathrm{MHD}+}$: Stellar mass of the central galaxy at $z=0$, in our non-CR, but otherwise full-physics (“MHD+”) run. (4) $M_{\ast }^{\mathrm{CR}+}$: Stellar mass of the central galaxy at $z=0$, in our “default” (observationally-favored) CR+ ($\kappa =3e29$) run. (5) $m_{i,\mathrm{\hspace{0.17em}1000}}$: Mass resolution: the baryonic (gas or star) particle/element mass, in units of $1000M_{☉}$. The DM particle mass is always larger by the universal ratio, a factor $\approx 5$. (6) ${⟨{ϵ}_{\mathrm{gas}}⟩}^{\mathrm{sf}}$: Spatial resolution: the gravitational force softening (Plummer-equivalent) at the mean density of star formation (gas softenings are adaptive and match the hydrodynamic resolution, so this varies), in the MHD+ run. Time resolution reaches $\sim 10-100$yr and density resolution $\sim {10}^3-{10}^4{\mathrm{cm}}^{-3}$. (7) Additional notes.

Table 2. Table 2: Additional high-redshift, massive-halo simulations ( M halo vir ≳ 10 12 ​ M ☉ greater-than-or-equivalent-to superscript subscript 𝑀 halo vir superscript 10 12 subscript 𝑀 ☉ M_{\rm halo}^{\rm vir}\gtrsim 10^{12}\,M_{\sun} ). All units are physical.

Halo/stellar properties (as Table 1) of simulations which feature halos that reach $\sim {10}^{12}{M}_{☉}$ at various redshifts $z$ (i.e. are more massive at $z\sim 0$). Lacking AGN feedback, we do not evolve these halos past masses $\sim {10}^{12-13}{M}_{☉}$. Columns show: (1) Simulation name. (2) $z_{12}$: Redshift at which the halo virial mass reaches $\sim {10}^{12}{M}_{☉}$. (3) $z_f$: Lowest redshift to which the simulation is evolved. (4) $M_{\ast }^{\mathrm{MHD}+}\left(z_f\right)$: Stellar mass (in $M_{☉}$) of the primary galaxy at $z_f$, in the MHD+ run. (5) $M_{\ast }^{\mathrm{CR}+}\left(z_f\right)$: Stellar mass (in $M_{☉}$) of the primary galaxy at $z_f$, in the CR+($\kappa =3e29$) run. (6) $m_{i,\mathrm{\hspace{0.17em}1000}}$: Mass resolution in $1000M_{☉}$. (7) Reference (paper in which this IC first appeared). (8) Additional notes.

Table 3. Table 3: Default physics “suites” in this paper

Description of the basic physics in the simulation sets described here. This notation is used throughout the paper.

Equations

$${\kappa }_{\mathrm{cond}}$$

$${\eta }_{\mathrm{visc}}$$

$$-2$$

$$\frac{\partial e_{\mathrm{cr}}}{\partial t}+\nabla \cdot {\mathbf{F}}_{\mathrm{cr}}$$

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

$$e_{\mathrm{cr}}$$

$$\frac{|\nabla P_{\mathrm{cr}}|}{|\rho \nabla \Phi |}\sim \frac{{\dot{E}}_{\mathrm{cr}}G{R}_s}{3\tilde{\kappa }{f}_{\mathrm{gas}}{V}_{\max }^4}\sim \frac{0.5\alpha {ϵ}_{\mathrm{cr},0.1}}{f_{\mathrm{gas},0.1}{\tilde{\kappa }}_{29}(1+z{)}^{3/2}}\left(\frac{M_{\ast }}{f_b{M}_{\mathrm{halo}}}\right)$$

$$\frac{{\dot{E}}_{\mathrm{loss}}}{{\dot{E}}_{\mathrm{cr}}}\approx \frac{\ln (r_{\max }/r_{\min })}{4\pi \tilde{\kappa }{m}_p{t}_0{n}_0}\frac{f_{\mathrm{gas}}{V}_{\max }^2}{G}\sim \frac{0.15}{{\tilde{\kappa }}_{29}}\left(\frac{{\Sigma }_{\mathrm{gas}}{R}_{\mathrm{gas}}}{0.01\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}\mathrm{kpc}}\right)$$

$$\frac{{\dot{E}}_{\mathrm{loss}}}{{\dot{E}}_{\mathrm{cr}}}$$

$$\frac{F_{\gamma }}{F_{\mathrm{SF}}}$$

$$r_{\mathrm{stream}}\sim 30\mathrm{kpc}\frac{{\tilde{\kappa }}_{29}}{v_{A,10}}\sim 30\mathrm{kpc}\frac{{\beta }_{100}^{1/2}{\tilde{\kappa }}_{29}}{M_{\mathrm{vir},12}^{1/3}(1+z{)}^{1/2}}$$

$$\frac{{\dot{E}}_{\mathrm{loss}}^{\mathrm{stream}}}{{\dot{E}}_{\mathrm{cr}}}(r\ll r_{\mathrm{stream}})=\frac{v_{A}r}{3\tilde{\kappa }}=\frac{1}{3}\left(\frac{r}{r_{\mathrm{stream}}}\right)$$

$$\frac{{\dot{E}}_{\mathrm{cr}}}{L_{\mathrm{cool}}^{\mathrm{tot}}}\sim \frac{0.16\alpha {ϵ}_{\mathrm{cr},0.1}(1+z{)}^{3/2}}{{\Lambda }_{22}(f_{\mathrm{gas}}/f_b)}\left(\frac{M_{\ast }}{f_b{M}_{\mathrm{halo}}}\right)\left(\frac{10^{-4}\mathrm{c}{\mathrm{m}}^{-3}}{⟨n{⟩}_{\mathrm{cool}}}\right)$$

$$\frac{{\dot{e}}_{\mathrm{heat}}^{\mathrm{stream}}}{{\dot{e}}_{\mathrm{cool}}}\sim \frac{0.02\alpha {ϵ}_{\mathrm{cr},0.1}}{{\Lambda }_{22}(f_{\mathrm{gas}}/f_b{)}^2(1+z{)}^{3/2}}\left(\frac{M_{\ast }}{f_b{M}_{\mathrm{halo}}}\right)\left(\frac{r}{R_{\mathrm{vir}}}\right)$$

$$\displaystyle\frac{\hat{\tilde{\bf F}}_{\rm cr}}{\tilde{c}^{2}}\,{\Bigl{[}}\frac{\partial|\tilde{\bf F}_{\rm cr}|}{\partial t}$$

$$+{\nabla }_{\Vert }{P}_{\mathrm{cr}}=-\frac{({\gamma }_{\mathrm{cr}}-1)}{{\kappa }_{\ast }}{\tilde{\mathbf{F}}}_{\mathrm{cr}}$$

Full text

1 Introduction

Galaxy and star formation are intrinsically multi-physics processes that involve a competition between gravity, collisionless dynamics, fluid dynamics and turbulence, radiation-matter coupling and chemistry, relativistic particles, magnetic fields, and more. Many of these processes enter most dramatically via “feedback” from massive stars, whereby the radiation, winds, and explosions from massive stars dramatically alter subsequent generations of star and galaxy formation.

In recent years, there has been tremendous progress in modeling and understanding the effects of multi-phase gas in the interstellar and circum/inter-galactic medium (ISM, CGM, and IGM), radiative cooling, turbulence, and self-gravity, and how these processes couple to stellar feedback. For example, it is now well-established that without feedback, these processes lead to runaway fragmentation and gravitational collapse that turns dense gas into stars on a single gravitational free-fall time (on both giant molecular cloud [GMC] and galactic scales; Tasker 2011; Hopkins et al. 2011; Dobbs et al. 2011; Harper-Clark & Murray 2011), and transforms most of the baryons in the Universe into stars (Katz et al., 1996; Somerville & Primack, 1999; Cole et al., 2000; Springel & Hernquist, 2003; Kereš et al., 2009), both in stark contrast to observations. Moreover direct and indirect effects of feedback are ubiquitously observed in outflows and enrichment of the CGM and IGM (Martin, 1999; Heckman et al., 2000; Pettini et al., 2003; Songaila, 2005; Martin et al., 2010; Sato et al., 2009; Weiner et al., 2009; Steidel et al., 2010). Of course, feedback (even restricting just to feedback from stars) comes in a variety of forms, including radiative (ionization, photo-heating and radiation pressure), mechanical (thermal and kinetic energy associated with supernovae Types Ia & II, stellar mass-loss, and protostellar jets), the injection of magnetic fields, and acceleration of cosmic rays. Numerical simulations have begun to directly resolve the relevant scales of some of these processes and therefore have begun to explicitly treat some of these feedback channels and their interactions on ISM and galactic or inter-galactic scales (e.g. Tasker, 2011; Hopkins et al., 2011, 2012; Wise et al., 2012; Kannan et al., 2014; Agertz et al., 2013; Roškar et al., 2014).

One example is the suite of simulations studied in the “Feedback In Realistic Environments” (FIRE)111See the FIRE project website:

http://fire.northwestern.edu

For additional movies and images of FIRE simulations, see:

http://www.tapir.caltech.edu/~phopkins/Site/animations/ project (Hopkins et al., 2014). These simulations have been used extensively in recent years to explore the interplay between a multi-phase ISM and CGM and both radiative (Hopkins et al., 2018a) and mechanical feedback (Hopkins et al., 2018c) processes from stars. These processes alone, coupled to the physics of radiative cooling, gravity, and star formation, appear to explain a wide variety of observed phenomena in galaxies, including their abundances (Ma et al., 2016; Escala et al., 2018), star formation “main sequence” and fluctuations (Sparre et al., 2017), satellite mass functions (Wetzel et al., 2016; Garrison-Kimmel et al., 2018), color distributions (Feldmann et al., 2016), stellar (Wheeler et al., 2015, 2017) and gas-phase (El-Badry et al., 2018b, a) kinematics, radial gradients and internal thick/thin disk structure (Ma et al., 2017a, b; Bonaca et al., 2017), stellar halos (Sanderson et al., 2017; El-Badry et al., 2018c), multi-phase fast outflows (Muratov et al., 2015, 2017), dark matter profiles (Oñorbe et al., 2015; Chan et al., 2015), and more.

However, essentially all of the conclusions above were based on simulations that treated the gas in the hydrodynamic limit. The detailed plasma physics of the ISM and CGM/IGM is of course quite complex and still a subject of active research. But it is well-established that magnetic fields are ubiquitous and important for local particle transport, that anisotropic transport processes (e.g. conduction and viscosity) can become significant in hot, tenuous gas, and that the ISM and IGM contain a spectrum of relativistic charged particles (cosmic rays [CRs]). In a very broad sense at least in the local Solar-neighborhood ISM, magnetic field and CR energy densities are order-of-magnitude comparable to thermal and turbulent energy densities (Ginzburg & Ptuskin, 1985; Boulares & Cox, 1990), suggesting they may not be negligible dynamically.

Magnetic fields have been studied and discussed in the galactic and extra-galactic context for decades (for reviews, see Beck et al., 1996; Beck, 2009). They can, in principle, slow star formation in dense gas (Piontek & Ostriker, 2005, 2007; Wang & Abel, 2009; Beck et al., 2012; Pakmor & Springel, 2013; Kim & Ostriker, 2015), or alter fluid mixing instabilities (Jun et al., 1995; McCourt et al., 2015; Armillotta et al., 2017) and the evolution of SNe remnants (Jun & Norman, 1996b, a; Jun & Jones, 1999; Thompson, 2000; Kim & Ostriker, 2015), and of course determine the actual dynamics of anisotropic transport. Anisotropic thermal conduction and viscosity in hot gas have also been studied extensively in the past (albeit not quite as widely), and it has been widely suggested that both could be important for plasma heating and dynamics on galaxy cluster scales (Reynolds et al., 2005; Sijacki & Springel, 2006; Markevitch & Vikhlinin, 2007; Sharma et al., 2009, 2010; Parrish et al., 2012; Choi & Stone, 2012; Armillotta et al., 2017) or, again, in SNe remnants (see references above), or for the mixing/survival of cool clouds in hot galactic outflows or the CGM (Brüggen & Scannapieco, 2016; Armillotta et al., 2017). However, it is not clear if these processes are particularly important for galaxy properties. In fact, most studies in the past have argued the effects of these physics in dwarfs and Milky Way (MW)-mass ($\sim L_{\ast}$) galaxies are relatively small – perhaps not surprising since the magnetic dynamo in super-sonic turbulence appears to saturate with magnetic fields always sub-dominant to turbulence (effectively, passively-amplified; see Federrath et al., 2014; Su et al., 2018a; Squire & Hopkins, 2017; Colbrook et al., 2017; Rieder & Teyssier, 2017; Banda-Barragán et al., 2018; Martin-Alvarez et al., 2018) and conduction/viscosity depend strongly on gas temperature and are weaker in the cooler gas of sub-$L_{\ast}$ galaxy halos.222Throughout this manuscript, we use “halo” to refer to the extended circum-galactic gas, stars, and/or dark matter extending from outside the “luminous” galaxy to around the virial radius, i.e. from $\sim 10-30\,$kpc to $\sim 200-400\,$kpc. This is standard notation in many fields (e.g. cosmology, CGM/IGM, galaxy formation), but we emphasize that this is quite different from what is usually called “halo” gas in much of the Galactic CR and magnetic field literature, where “halo” is often used to refer to gas within a scale height $\lesssim 1-10\,$kpc of the disk midplane out to $\sim 10\,$kpc (i.e. the “thick disk” in standard galaxy/CGM notation; see Haverkorn 2015; Han 2017). In Su et al. (2017), we attempted to study the effect of magnetic fields, anisotropic conduction, and viscosity in addition to the physics described above in FIRE simulations, and concluded the effects were minimal. However, the simulations in that paper were mostly non-cosmological (although it did include two cosmological cases), so might not capture all the important effects in the CGM. Moreover, a swathe of recent work in plasma physics has argued that conductivity and viscosity of dilute plasmas might be self-limiting under exactly the relevant conditions of the ISM and CGM (see Kunz et al., 2014; Riquelme et al., 2016; Roberg-Clark et al., 2018; Squire et al., 2017b; Komarov et al., 2018, and references therein), and these effects were not accounted for in previous studies (although they generally act to weaken the conductivity and viscosity).

The situation with CRs is much less clear. In the MW (and, it is widely believed, most dwarf and $\sim L_{\ast}$ or star-forming galaxies), the CR pressure and energy density (and correspondingly, effects on both gas dynamics and heating/cooling rates of gas via hadronic or Coulomb collisions or excitation of Alfvén waves in various plasma “streaming instabilities”; Mannheim & Schlickeiser 1994; Enßlin et al. 2007; Guo & Oh 2008) are dominated by mildly-relativistic $\sim$ GeV protons accelerated primarily in supernova remnants (with $\sim 10\%$ of the SNe ejecta energy ultimately in CRs; Bell 2004). In more massive galaxies, which host supermassive black holes (BHs) and have little star formation, most of the CR production appears to be associated with AGN jets and “bubbles.” CRs and their influence on galaxy evolution have been a subject of interest in both analytic (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 numerical simulation (Jubelgas et al., 2008; Uhlig et al., 2012; Booth et al., 2013; Salem & Bryan, 2014; Ruszkowski et al., 2017; Farber et al., 2018a) studies for decades – with an explosion of work in recent years. This work has argued that CRs could, in principle, drive galactic outflows (Simpson et al., 2016; Girichidis et al., 2016; Pakmor et al., 2016; Wiener et al., 2017), suppress star formation in low (or high) mass galaxies (Hanasz et al., 2013; Chen et al., 2016; Jacob et al., 2018), provide additional pressure to “thicken” galactic gas disks (Wiener et al., 2013b; Salem et al., 2014), alter the phase structure of the CGM (Salem et al., 2016; Girichidis et al., 2018; Butsky & Quinn, 2018), “open up” magnetic field lines or otherwise alter the galactic dynamo (Parker, 1992; Hanasz et al., 2009; Kulpa-Dybeł et al., 2011, 2015), and more.

However, a number of major uncertainties and limitations remain in this field. First, the actual CR transport processes, and their coupling to the gas, remain deeply uncertain (owing to the extremely complicated plasma processes involved) – the physics that gives rise to some “effective diffusivity” and/or CR streaming is still debated (see Strong et al., 2007; Zweibel, 2013; Grenier et al., 2015), and there is no widely-accepted a priori model which predicts the relevant transport coefficients in the way of, say, Spitzer-Braginskii conductivity and viscosity. There are some empirical constraints from e.g. $\gamma$-ray emission in nearby galaxies or more detailed products (e.g. spallation) in the MW, but critically any inferred constraint on the “effective diffusion coefficient” or “streaming speed” of CRs is strongly model-dependent (as it depends on the density distribution the CRs propagate through, the magnetic field configuration, etc.). Really, one must forward-model these constraints in any galaxy model, to test whether the adopted CR transport assumptions are consistent with the observations. Second, almost all previous studies of CRs on galaxies either focused on (a) idealized “patches” of the ISM or CGM, ignoring the global dynamics of accretion, outflows, star formation, etc., or (b) galaxy simulations with (intentionally) highly-simplified models for the turbulent, multi-phase ISM, star formation, stellar feedback from supernovae, stellar winds, radiation, and more. But these details are critical for determining the balance of CR heating and cooling, how CRs will be trapped or escape galaxies, whether CRs will influence outflows or gas in the CGM “lofted up” by other processes, and whether CRs ultimately matter compared to the order-of-magnitude larger energy input in mechanical (thermal+kinetic) form in SNe. To give an extreme example: almost anything will have a large effect relative to a “baseline” model which includes weak or no stellar feedback. It is much less clear whether the inclusion of CR physics will “matter” once mechanical and radiative feedback from stars is already accounted for.

Working towards this goal, Chan et al. (2018) performed and presented the first simulations combining the specific physics from the FIRE simulations, described above, with explicit CR injection and transport, accounting for advection and fully-anisotropic streaming and diffusion, as well as hadronic and Coulomb collisional and streaming (Alfvén) losses. They systematically varied the transport coefficients and treatment, and compared with observational constraints, to argue that – at least given this particular physics set and treatment of CR transport – the observations required diffusivities $\gtrsim 10^{29}\,{\rm cm^{2}\,s^{-1}}$, and that within the allowed range of diffusivities, the effects on galaxy star formation rates and gas density distributions were modest. However, these simulations were restricted to non-cosmological, isolated galaxies, representative of just a couple of $z=0$ galaxy types (e.g. one dwarf, one MW-mass system). Predicting the consequences for galactic winds or the CGM and therefore long-term galaxy evolution (e.g. stellar masses, etc.) requires cosmological simulations.

In this paper, we therefore introduce and explore a new, high-resolution suite of FIRE simulations, with $>150$ fully-cosmological simulations spanning halo masses from ultra-faint dwarfs through MW-mass systems at a range of redshifts (reaching $\sim$ pc-scale resolution), and systematically exploring all of the physics above. Specifically, we compare our standard physics assuming hydrodynamics, to simulations with explicit MHD and anisotropic conduction and viscosity as in Su et al. (2017), and simulations with all of the above plus explicit treatment of cosmic rays as in Chan et al. (2018). We moreover systematically survey the treatment of CR transport physics and coefficients, and compare with observations where possible to constrain the allowed range of assumptions. Our intention here is to identify which physics might have an influence on bulk galaxy properties (e.g. SFRs, stellar masses, morphologies), and where uncertain parameters exist (e.g. CR diffusivities), what range of those parameters is allowed and how the effects (if any) on galaxies depend on them within the allowed range. We also limit our study to dwarf and $\sim L_{\ast}$ galaxies where it is widely believed that SNe dominate the CR injection. In companion papers (e.g. Su et al. 2018d, Su et al., in prep.) we will study the complementary role of AGN injecting CRs in much more massive galaxies, and in other companion papers (e.g. Ji et al. 2019 and Chan et al., in prep) we will study the (potentially much larger) effects of CRs on the CGM around galaxies and the origin and properties of the weak, CR-driven outflows.

In § 2 we review the numerical methods and describe the simulation suite. Before analyzing the simulations, § 3 presents a simple analytic model for the effects and equilibrium distribution of CRs, given our assumptions in the simulations, which allows us to predict and estimate (with surprising accuracy) many of the scalings we will observe in the cosmological simulations. § 4 briefly presents the key results from the simulations, which we discuss and analyze in more detail – attempting to break down the effects of different physics on different scales – in § 5. Finally, we summarize and conclude in § 6.

2 Methods

The simulations in this paper were run with the multi-physics code GIZMO333A public version of GIZMO is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html (Hopkins, 2015), in its meshless finite-mass MFM mode. This is a mesh-free, finite-volume Lagrangian Godunov method which provides adaptive spatial resolution together with conservation of mass, energy, momentum, and angular momentum, and the ability to accurately capture shocks and fluid mixing instabilities (combining advantages of both grid-based and smoothed-particle hydrodynamics methods). We solve the equations of ideal magneto-hydrodynamics (MHD), as described and tested in detail in Hopkins & Raives (2016); Hopkins (2016), with fully-anisotropic Spitzer-Braginskii conduction and viscosity and other diffusion operators implemented as described in Hopkins (2017). Gravity is solved for gas and collisionless (stars and dark matter) species with fully-adaptive Lagrangian force softening (so hydrodynamic and force resolutions are consistently matched).

Our simulations are fully-cosmological “zoom-in” runs, with a high-resolution Lagrangian region identified surrounding a $z=0$ “primary” halo of interest in a large cosmological box (see Oñorbe et al., 2014). Tables 1-2 list the specific volumes run, and the properties of each “primary” halo.444For the MUSIC (Hahn & Abel, 2011) files necessary to generate all ICs here, see:

http://www.tapir.caltech.edu/~phopkins/publicICs We note that the high-resolution volumes reach as large as $\sim(10\,{\rm Mpc})^{3}$ in the largest runs, so there are many galaxies present (our set has hundreds of galaxies with $>100$ star particles each). Also, although we focus on $\sim 30$ zoom-in volumes, because we systematically vary the physics and CR parameters like the diffusion coefficient, our default simulation set includes well over 100 full-physics high-resolution simulations. However, to simplify our analysis and presentation, avoid ambiguities in galaxy matching and separating systematic differences between satellite and field galaxies, and to focus on the best-resolved galaxies possible, we focus only on the most massive “primary” galaxies in each box. We note though that a brief comparison indicates that our conclusions appear to apply to all galaxies in the box. Unless otherwise specified (e.g. Table 2), all are run to $z=0$.

2.1 Default FIRE-2 (“Hydro+”) Physics

All our simulations here include the physics of cooling, star formation, and stellar feedback from the FIRE-2 version of the Feedback in Realistic Environments (FIRE) project, described in detail in Hopkins et al. (2018b), but briefly summarized here.

Gas cooling is followed from $T=10-10^{10}\,$K including free-free, Compton, metal-line, molecular, fine-structure and dust collisional, and we also follow gas heating from photo-electric and photo-ionization by both local sources and a uniform meta-galactic background including the effect of self-shielding.555As detailed in Hopkins et al. (2018b) Appendix B, in our runs that do not include explicit CR transport (“Hydro+” and “MHD+”), the cooling/ionization tables do assume a uniform MW-like CR background ($\sim 1\,{\rm eV\,cm^{-3}}$) for gas at densities $>0.01\,{\rm cm^{-3}}$. This is generally negligible for heating, but is important for e.g. the small ionized fraction in GMCs. In our runs with explicit CR transport, these terms are replaced with the explicitly-evolved CR background and collisional+streaming heating rates described below. We explicitly follow 11 different abundances, including explicit treatment of turbulent diffusion of metals and passive scalars as in Colbrook et al. (2017); Escala et al. (2018). Gas is turned into stars using a sink-particle prescription: gas which is locally self-gravitating at the resolution scale following Hopkins et al. 2013c, self-shielding/molecular following Krumholz & Gnedin 2011, Jeans unstable, and denser than $n_{\rm crit}>1000\,{\rm cm^{-3}}$ is converted into star particles on a free-fall time. Star particles are then treated as single-age stellar populations with all IMF-averaged feedback properties calculated from STARBURST99 (Leitherer et al., 1999) assuming a Kroupa (2001) IMF. We then explicitly treat feedback from SNe (both Types Ia and II), stellar mass loss (O/B and AGB mass-loss), and radiation (photo-ionization and photo-electric heating and UV/optical/IR radiation pressure), with implementations at the resolution-scale described in Hopkins et al. (2018b), Hopkins et al. (2018c), and Hopkins et al. (2018a).

The simulations labeled “Hydro+” in this paper include all of the physics above, but do not include magnetic fields, physical conduction or viscosity, or explicit treatment of cosmic rays.

2.2 Magnetic Fields, Conduction, & Viscosity (“MHD+”)

Our simulations labeled “MHD+” in this paper include all of the “Hydro+” physics (e.g. radiative cooling, star formation, stellar feedback), but add magnetic fields and physical, fully-anisotropic conduction and viscosity. These are described in Su et al. (2017) but we briefly summarize here.

As noted above, for magnetic fields we solve the equations of ideal MHD, as described in Hopkins & Raives (2016); Hopkins (2016). We can optionally include ambipolar diffusion, the Hall effect, and Ohmic resistivity as in Hopkins 2017, but these are completely negligible at all resolved scales in our simulations here.666We have, in fact, run one m10q and one m12i run from $z=0.1$ to $z=0$ with the full non-ideal MHD terms enabled and the coefficients calculated following Zhu et al. 2014, to confirm that they are negligible corrections to the induction equation at the scales resolved here (usually by several orders of magnitude, except for ambipolar diffusion, which can be a $\sim 10\%$-level correction to the induction equation in the dense, neutral gas in GMCs, but this is far smaller than other physical uncertainties).

For conduction and viscosity we include the physical scalings for fully-anisotropic Spitzer-Braginskii transport, with the conductive heat flux $\kappa_{\rm cond}\,\hat{\bf B}\,(\hat{\bf B}\cdot\nabla T)$, where

[TABLE]

(Spitzer & Härm, 1953; Braginskii, 1965), where $m_{e}$ and $e$ are the electron mass and charge, $k_{B}$ the Boltzmann constant, $f_{i}$ is the ionized fraction, $\ln{\Lambda_{c}}\sim 37$ is a Coulomb logarithm and $\ell_{e}=3^{3/2}\,(k_{B}\,T)^{2}/4\,n_{e}\,\pi^{1/2}\,e^{4}\,\ln{\Lambda_{c}}$ is the electron deflection length (Sarazin, 1988), $\ell_{T}\equiv T/|\hat{\bf B}\cdot\nabla T|$ is the parallel temperature gradient scale-length, and the plasma $\beta\equiv P_{\rm thermal}/P_{\rm magnetic}$ is the usual ratio of thermal-to-magnetic pressure. Note that the $\ell_{e}/\ell_{T}$ term ensures proper behavior in the saturated limit (and a smooth transition between un-saturated and saturated limits, e.g. Cowie & McKee 1977), while the $\beta$ term accounts for micro-scale plasma instabilities (e.g. the Whistler instability) limiting the heat flux in high-$\beta$ plasmas (see Komarov et al., 2018). For viscosity we modify the momentum and energy equations with the addition of the viscous stress tensor $\Pi\equiv-3\,\eta_{\rm visc}\,(\hat{\bf B}\otimes\hat{\bf B}-\mathbb{I}/3)\,(\hat{\bf B}\otimes\hat{\bf B}-\mathbb{I}/3):(\nabla\otimes{\bf v})$ (where $\otimes$ is the outer product, $\mathbb{I}$ the identity matrix, and $:$ the double-dot-product), with

[TABLE]

where like the above $m_{i}$, $\ell_{i}$, and $Z_{i}\,e$ are the mean ion mass, deflection-length, and charge, and $\ell_{|v|}\equiv|{\bf v}|/|\hat{\bf B}\otimes\hat{\bf B}:\nabla\otimes{\bf v}|$ is the parallel velocity gradient scale-length. The upper and lower allowed values of the anisotropic stress are limited (capped) according to the latter expression, to account again for plasma instabilities (the mirror and firehose at positive and negative anisotropy, respectively) limiting the flux (see Squire et al., 2017c, a, b).

These numerical methods have been extensively tested and discussed in previous papers, to which we refer for details (e.g. Hopkins & Raives, 2016; Hopkins, 2016, 2017; Hopkins & Conroy, 2017; Su et al., 2017, 2018a; Lee et al., 2017; Colbrook et al., 2017; Seligman et al., 2018). We note that, given the mass resolution and Lagrangian nature of the code, the Field (1965) length $\lambda_{\rm F}$ (approximately, the scale below which thermal conduction is faster than cooling) is resolved ($\Delta x<\lambda_{\rm F}$) in the simulations here in gas hotter than $T\gtrsim 2\times 10^{5}\,{\rm K}\,(n/0.01\,{\rm cm^{-3}})^{0.4}$.

2.3 Cosmic Rays (“CR+”)

Our simulations labeled “CR+” in this paper include all of the “MHD+” physics, and add our “full physics” treatment of CRs. The CR physics is described in detail in Chan et al. (2018) but we again summarize here: we include injection in SNe shocks, fully-anisotropic CR transport with streaming and advection/diffusion, CR losses (hadronic and Coulomb, adiabatic, streaming), and self-consistent CR-gas coupling.

CRs are treated as an ultra-relativistic fluid777The fluid limit adopted here is motivated by the fact that $\sim$ GeV CR gyro radii $r_{\rm L}\ll{\rm au}$ are always much smaller than any resolved spatial scales, and CR “deflection lengths” or scattering lengths $\sim\kappa/c\sim$ pc are also always smaller than e.g. the CR pressure gradient scale lengths. However, we emphasize that fundamental questions remain about the validity of other assumptions (like the form of the CR distribution function) when the CR diffusivity is large (scattering rates are low), particularly when $|\kappa_{\|}\nabla_{\|}P_{\rm cr}|\gg v_{A}\,(e_{\rm cr}+P_{\rm cr})$, which occurs in our favored models at galacto-centric radii $\ll 30\,$kpc, as discussed below. (adiabatic index $\gamma_{\rm cr}=4/3$) in the “single bin” approximation, which we can think of either as evolving only the CR energy density ($e_{\rm cr}$) at $\sim$ GeV energies that dominate the CR pressure ($P_{\rm cr}\equiv(\gamma_{\rm cr}-1)\,e_{\rm cr}$), or equivalently assuming a universal CR energy spectral shape. CR pressure contributes to the total pressure and effective sound speed in the Riemann problem for the gas equations-of-motion according to the local strong-coupling approximation, i.e. $P=P_{\rm gas}+P_{\rm cr}$ and $c_{s,\,{\rm eff}}^{2}=\partial P/\partial\rho=c_{s}^{2}+\gamma_{\rm cr}\,P_{\rm cr}/\rho$. Integrating over the CR distribution function and energy spectrum, we evolve the CR energy density as McKenzie & Voelk (1982):

[TABLE]

where $S_{\rm cr}$ and $\Gamma_{\rm cr}$ are source and sink terms; $\langle{\bf v}_{\rm cr}\rangle\equiv{\bf v}_{\rm gas}+{\bf v}_{\rm stream}$ is the bulk CR advection velocity, de-composed into the gas velocity ${\bf v}_{\rm gas}$ and the “streaming velocity” ${\bf v}_{\rm stream}$; and ${\bf F}_{\rm cr}$ is the lab-frame or total CR energy flux which can be de-composed into ${\bf F}_{\rm cr}\equiv\langle{\bf v}_{\rm cr}\rangle\,(e_{\rm cr}+P_{\rm cr})+{\bf F}_{\rm di}\equiv{\bf v}_{\rm gas}\,(e_{\rm cr}+P_{\rm cr})+\tilde{\bf F}_{\rm cr}$, where $\tilde{\bf F}_{\rm cr}$ is the flux in the fluid frame (with ${\bf F}_{\rm di}$ the “diffusive” flux).

For $S_{\rm cr}$, we assume CR injection in SNe shocks, with a fixed fraction $\epsilon_{\rm cr}$ of the initial ejecta kinetic energy ($\approx 10^{51}\,{\rm erg}$) of every SNe going into CRs. In our default simulations, $\epsilon_{\rm cr}=0.1$ is adopted. This is coupled directly to gas in the immediate vicinity of each explosion, alongside the thermal and kinetic energy, mass, and metals, as described in Hopkins et al. (2018c). For $\Gamma_{\rm cr}$, we follow Guo & Oh (2008) and account for both hadronic/catastropic and Coulomb losses with $\Gamma_{\rm cr}=e_{\rm cr}\,n_{\rm n}\,(1+0.28\,x_{e})/t_{0}\,n_{0}$, where $n_{\rm n}$ is the nucleon number density, $t_{0}\,n_{0}\equiv 1.72\times 10^{15}\,{\rm s\,cm^{-3}}$, and $x_{e}$ is the number of free electrons per nucleon. A fraction $\sim 1/6$ of the hadronic products and all Coulomb losses thermalize, and contribute a volumetric gas heating term $Q_{\rm gas}=e_{\rm cr}\,n_{\rm n}\,(0.17+0.28\,x_{e})/t_{0}\,n_{0}$.

As shown in Chan et al. (2018) the remaining terms in Eq. 4 can be decomposed (in Lagrangian form) into simple advection (automatically handled in our Lagrangian formulation) and adiabatic “PdV work” terms (solved in an exactly-conservative manner with our usual MFM solver), the $\tilde{\bf F}_{\rm cr}$ term (see below), and a “streaming loss” term888As discussed at length in Chan et al. (2018) and studied in the Appendices therein, we have considered simulations where the “streaming losses” scale either as ${\bf v}_{\rm stream}\cdot\nabla P_{\rm cr}$ or ${\bf v}_{A}\cdot\nabla P_{\rm cr}$. The latter (streaming losses limited to the Alfvén speed, as they arise from damped Alfvén waves) is our default choice, even if $v_{\rm stream}>v_{A}$. However we do not find this choice significantly alters any of our conclusions in this paper. ${\bf v}_{A}\cdot\nabla P_{\rm cr}$ which is negative definite and represents energy loss via streaming instabilities that excite high-frequency Alfvén waves (frequency of order the CR gyro frequency, well below our resolution; Wentzel 1968; Kulsrud & Pearce 1969) that damp and thermalize almost instantaneously, so the energy lost via this term is added to the gas thermal energy each timestep. The streaming velocity always points down the CR pressure gradient, projected along the magnetic field, so ${\bf v}_{\rm stream}=-v_{\rm stream}\,\hat{\bf B}\,(\hat{\bf B}\cdot\hat{\nabla}P_{\rm cr})$.

It is widely argued that micro-scale instabilities regulate the streaming speed to of order the Alfvén speed (Skilling, 1971; Kulsrud, 2005; Yan & Lazarian, 2008; Enßlin et al., 2011), although super-Alfvénic streaming can easily emerge in self-confinement models for CR transport (Wentzel, 1968; Holman et al., 1979; Achterberg, 1981; Wiener et al., 2013a; Lazarian, 2016); moreover in partially-neutral gas (ionized fraction $f_{\rm ion}<1$) there is an ambiguity about whether the appropriate Alfvén speed is the ideal-MHD Alfvén speed $v_{A}\equiv|\bf B|/(4\pi\,\rho)^{1/2}$ or the (larger) ion Alfvén speed $v_{A}^{\rm ion}\equiv|{\bf B}|/({4\pi\,\rho_{\rm ion}})^{1/2}\sim f_{\rm ion}^{-1/2}\,v_{A}$ (Skilling, 1975; Zweibel, 2013; Farber et al., 2018b). So we simply adopt the ad-hoc $v_{\rm stream}\approx 3\,v_{A}$ as our default, although we vary this widely below setting $v_{\rm stream}=v_{A}$, disabling streaming entirely, or allowing highly super-sonic/Alfvénic streaming with $v_{\rm stream}=3\,(v_{A}^{2}+c_{s}^{2})^{1/2}$ (several times the fastest possible MHD wavespeed), and show (both here and in Chan et al. 2018) it has little effect on our conclusions.

Following Chan et al. (2018) we treat $\tilde{\bf F}_{\rm cr}$ using a two-moment scheme (similar to other recent implementations by e.g. Jiang & Oh 2018 and Thomas & Pfrommer 2018), solving

[TABLE]

where $\nabla_{\|}P_{\rm cr}=(\hat{\bf B}\otimes\hat{\bf B})\cdot(\nabla P_{\rm cr})=(\gamma_{\rm cr}-1)\,\hat{\bf B}\,(\hat{\bf B}\cdot\nabla e_{\rm cr})$ is the parallel derivative of the CR pressure tensor, $\tilde{c}$ is the maximum allowed CR “free streaming” speed,999Following Chan et al. (2018) we note that $\tilde{c}$ is a nuisance parameter (the simulations evolve to identical solutions independent of $\tilde{c}$, so long as it is faster than other bulk flow speeds in the problem), so rather than adopt the microphysical $\tilde{c}=c$ (speed of light), we adopt $\tilde{c}=1000\,{\rm km\,s^{-1}}$ by default (but show below that our solutions are independent of $\tilde{c}$ over a wide range). and $\kappa_{\ast}\equiv\kappa_{\|}+\gamma_{\rm cr}\,v_{\rm st}\,P_{\rm cr}/|\nabla_{\|}P_{\rm cr}|$ is the effective parallel diffusivity (we are implicitly taking the perpendicular $\kappa_{\bot}=0$). We note this reduces to the simpler pure-anisotropic diffusion+streaming equation in steady state and/or on large spatial and time-scales ($\tilde{\bf F}_{\rm cr}\rightarrow-\kappa_{\ast}\nabla_{\|}e_{\rm cr}=\kappa_{\|}\nabla_{\|}e_{\rm cr}+{\bf v}_{\rm stream}\,(e_{\rm cr}+P_{\rm cr})$), but unlike a pure-diffusion equation (where one forces $\tilde{\bf F}_{\rm cr}$ to always be exactly $-\kappa_{\ast}\nabla_{\|}e_{\rm cr})$) it correctly handles the transition between streaming and diffusion and prevents un-physical super-luminal CR transport.101010Super-luminal CR transport would occur in the “pure-diffusion” approximation for $\tilde{\bf F}_{\rm cr}$ wherever the resolution scale $\Delta x\lesssim\kappa/c\sim 3\,{\rm pc}\,(\kappa/3\times 10^{29}\,{\rm cm^{2}\,s^{-1}})$. The simulations in this paper routinely reach this or better spatial resolution, so this distinction is important.

As discussed in Chan et al. (2018), if the “streaming loss” term is limited or “capped” to scale with the Alfvén speed ($\sim v_{A}\nabla P_{\rm cr}$; see above), and streaming is super-Alfvénic ($v_{\rm stream}\gg v_{A}$), then only the “effective” diffusivity $\kappa_{\ast}$ (which can arise from a combination of microphysical diffusion and/or streaming) – as compared to $\kappa_{\|}$ or $v_{\rm stream}$ individually – enters the large-scale dynamics. This effective $\kappa_{\ast}$, or equivalent CR transport speed $v_{\rm cr,\,eff}\sim\kappa_{\ast}\,|\nabla P_{\rm cr}|/P_{\rm cr}$, is what we actually constrain in our study here.

There are many approximations in this description, and the effective “diffusion coefficient” $\kappa_{\ast}$ for CRs on these (energy, spatial, and time) scales remains both theoretically and observationally uncertain. Therefore, we treat $\kappa$ as a constant but vary it systematically in a parameter survey, with values motivated by the comparison with observational constraints in Chan et al. (2018).

3 A Theoretical Toy Model for CRs

In this section, we develop a simple theoretical “toy model” for CRs, which provides considerable insight into the phenomena we see in the simulations.

Assume a galaxy has a quasi-steady SFR $\dot{M}_{\ast}$, so the associated CR injection rate is $\dot{E}_{\rm cr}=\epsilon_{\rm cr}\,\epsilon_{\rm SNe}\,\dot{M}_{\ast}$, where $\epsilon_{\rm SNe}\sim(10^{51}\,{\rm erg}/70\,M_{\odot})$ is the time-and-IMF-averaged energetic yield per solar mass of star formation from SNe (for the IMF here including Type-II and prompt Ia events). Also assume that (since we are primarily interested in large scales) the CR injection is concentrated on relatively small scales, so it can be approximated as point-like, and assume – for now – that diffusion with constant effective coefficient $\tilde{\kappa}\sim\langle|\hat{\bf B}\cdot\hat{\nabla}e_{\rm cr}|^{2}\rangle\,\kappa_{\ast}\sim\kappa_{\ast}/3$ dominates the transport (e.g. collisional and streaming losses are negligible; we will return to these below). This has the trivial equilibrium solution:111111Eq. 6 is still approximately valid (up to an order-unity constant) in most cases for a non-constant $\tilde{\kappa}$ or effective transport/streaming speed $v_{\rm cr}$ (replacing $\tilde{\kappa}\rightarrow v_{\rm cr}\,r$). This follows from the fact that the radial CR flux (in a spherical hydrostatic system with no losses) is ${\bf F}_{\rm cr}\rightarrow-\bar{\kappa}\,(\partial e_{\rm cr}/\partial r)\,\hat{r}\sim-\bar{\kappa}\,(e_{\rm cr}/r)\,\hat{r}\sim v_{\rm cr}\,e_{\rm cr}\,\hat{r}$.

[TABLE]

This gives rise to the CR pressure gradient $\nabla P_{\rm cr}=-e_{\rm cr}\,\hat{r}/(3\,r)$.

Now also assume, for simplicity, the (diffuse) gas and DM are in an isothermal sphere (the detailed profile shape is not important) with gas fraction $f_{\rm gas}$ and circular velocity $V_{\rm c}\sim V_{\rm max}$, with some characteristic halo scale radius $R_{s}\sim R_{\rm vir}/c$. The ratio of the CR pressure gradient to the gravitational force at $R_{s}$ is ${|\nabla P_{\rm cr}|}/{|\rho\,\nabla\Phi|}\sim\dot{E}_{\rm cr}\,G\,R_{s}/(3\,\tilde{\kappa}\,f_{\rm gas}\,V_{\rm max}^{4})$ (this is identical up to an order-unity constant if we assume the gas is in e.g. an NFW profile or a Mestel disk). Now recall, $\dot{E}_{\rm cr}\propto\dot{M}_{\ast}$. For star-forming (sub-$L_{\ast}$) galaxies, $\dot{M}_{\ast}=\alpha\,M_{\ast}/t_{\rm Hubble}$ where $t_{\rm Hubble}(z)$ is the Hubble time at redshift $z$ and $\alpha\approx 1-2$ is (in a time-and-sample averaged sense) very weakly dependent on galaxy mass (see e.g. Mitra et al., 2017). Using this and the canonical scaling relations121212For Eq. 7, we approximate the system around the radii of interest as an isothermal sphere with $V_{c}\approx V_{\rm max}$, gas $\rho=f_{\rm gas}\,V_{c}^{2}/(4\pi\,G\,r^{2})$, $|\nabla\Phi|=G\,M_{\rm enc}/r^{2}=V_{c}^{2}/r$, and $V_{\rm max}\approx(2\,G\,M_{\rm vir}/R_{\rm vir})^{1/2}$ (a reasonable approximation for NFW-type halos with concentrations $\sim 1-20$), and $\dot{E}_{\rm cr}=\epsilon_{\rm cr}\,u_{\rm SNe}\,\dot{M}_{\ast}=\epsilon_{\rm cr}\,u_{\rm SNe}\,(\alpha\,M_{\ast}/t_{\rm Hubble}[z])$ (with $u_{\rm SNe}\approx 10^{51}\,{\rm erg}/70\,M_{\sun}$). For simplicity we adopt redshift scalings for a matter-dominated Universe: $M_{\rm halo}=M_{\rm vir}\equiv(4\pi/3)\,\Delta_{c}\,\rho_{c}\,R_{\rm vir}^{3}$ with $\Delta_{c}\approx 180$ (and $\rho_{c}=3\,H_{0}^{2}\,(1+z)^{3}/8\pi\,G$ the critical density) and $t_{\rm Hubble}[z]\approx 14\,(1+z)^{-3/2}\,{\rm Gyr}$, and evaluate everything at a halo scale radius $r\sim R_{s}\sim 0.1\,R_{\rm vir}$, but the choice of cosmology or exact radius do not qualitatively alter our conclusions. for halo $R_{s}$, $R_{\rm vir}$, and $M_{\rm halo}=M_{\rm vir}$, we can re-write this ratio as:

[TABLE]

where $\epsilon_{{\rm cr},\,0.1}\equiv\epsilon_{\rm cr}/0.1$, $f_{{\rm gas},\,0.1}\equiv f_{\rm gas}/0.1$, $\tilde{\kappa}_{29}\equiv\tilde{\kappa}/10^{29}\,{\rm cm^{2}\,s^{-1}}$, and $f_{b}\equiv\Omega_{b}/\Omega_{m}$ is the universal baryon fraction, so ${M_{\ast}}/{f_{b}\,M_{\rm halo}}$ is the stellar mass relative to what would be obtained converting all baryons into stars.

This leads immediately to the prediction that CRs can have a large effect in relatively massive (intermediate, LMC-like through MW-mass; $M_{\rm halo}\sim 10^{11-12}\,M_{\sun}$) halos, at low-to-moderate redshifts ($z\lesssim 1-2$), while their effect is limited at high redshifts or in low-mass halos: $M_{\ast}/M_{\rm halo}$ is a strongly-increasing function of mass (roughly, $M_{\ast}\propto M_{\rm halo}^{2}$ at low masses). So for a MW mass halo, CR pressure is approximately able to balance gravitational forces for $\tilde{\kappa}_{29}\sim 1$, while for a true dwarf halo like m10q, the CR pressure will be order-of-magnitude too small, because the SFR (hence the CR injection rate, which is proportional to the CR energy density/pressure/pressure gradient in steady state) is several orders of magnitude smaller in such a tiny dwarf. We also show below that in dwarfs, CRs escape efficiently from both the galaxy and CGM, while SNe cool less efficiently, all of which make CRs relatively less-dominant compared to mechanical stellar feedback.

Note also that our derivation above, at face value, would imply that a lower diffusion coefficient would give stronger effects of CR pressure (because the CRs are more “bottled up” so the steady-state pressure is larger). However, we neglected collisional and streaming losses: we show below that if $\tilde{\kappa}_{29}\ll 1$, these quickly dominate and prevent CRs from having any significant effects.

3.1 Collisional Losses at Low-$\kappa$ & Observational Constraints

We neglected collisional losses above. Assuming ionized gas (the difference is small, this just determines the contribution from Coulomb terms) these scale as $\dot{e}_{\rm loss}=e_{\rm cr}\,n/(t_{0}\,n_{0})$ with $t_{0}\,n_{0}\sim 1.7\times 10^{15}\,{\rm s\,cm^{-3}}$ (§ 2.3). If we take the same isothermal sphere model above, and calculate the steady-state volume-integrated CR loss rate $\dot{E}_{\rm loss}$, we obtain:

[TABLE]

where in the latter equality we use $f_{\rm gas}\,V_{\rm max}^{2}/G=M_{\rm gas}/R_{\rm gas}\approx\pi\,\Sigma_{\rm gas}\,R_{\rm gas}$ (and $\ln{(r_{\rm max}/r_{\rm min})}\sim 5$). Using the observational facts that (sub-$L_{\ast}$, star-forming) galactic disks have approximately constant effective $\Sigma_{\rm gas}\sim 10-30\,M_{\odot}\,{\rm pc^{-2}}$ and $R_{\rm gas}\sim 2\,R_{\ast}\sim 10\,{\rm kpc}\,(M_{\ast}/10^{11}\,M_{\odot})^{1/3}$ (e.g. Courteau et al., 2007), we can equivalently write:

[TABLE]

Of course, $\dot{E}_{\rm loss}$ cannot exceed $\dot{E}_{\rm cr}$ in steady-state: this gives an upper bound to the expression above at the“calometric limit” ($\dot{E}_{\rm loss}=\dot{E}_{\rm cr}$) at which point all input CR energy is lost to collisions.

The ratio $\dot{E}_{\rm loss}/\dot{E}_{\rm cr}$ is directly related to the observed ratio of GeV $\gamma$-ray flux or luminosity ($F_{\gamma}\propto L_{\rm gamma}\propto\dot{E}_{\rm loss}$) to bolometric flux from massive stars ($F_{\rm SF}\propto L_{\rm SF}\propto\dot{M}_{\ast}\propto\dot{E}_{\rm cr}$). Using the values given in § 4 below for conversion factors between collisional losses and $\gamma$-ray luminosity or flux, Eq. 8 becomes:

[TABLE]

with the calorimetric limit “capping” this at $F_{\gamma}/F_{\rm SF}\sim 2\times 10^{-4}$.

Note that the geometry of the gas is not especially important here. If we assume instead the gas is in a thin, exponential disk with scale height $H/R\ll 1$, and effective radius $R_{e}$, but the CRs diffuse approximately spherically (as a random walk), then we obtain, a result which differs only by an order-unity constant from the isothermal-sphere scaling above even as $H/R\rightarrow 0$ (because even though the disk occupies a vanishingly small volume, for fixed surface density $\Sigma$, the three-dimensional density $n$ in the disk, hence the CR loss rate while within it, increases inversely with $H/R$).

Three immediate consequences follow from this. First, essentially all observed galaxies with SFRs below $\sim 10\,{\rm M_{\odot}\,yr^{-1}}$ (including the MW, LMC, SMC, M31, and M33) have observed $F_{\gamma}/F_{\rm SF}\leq 10^{-5}$ (Lacki et al., 2011), so we predict that a diffusion coefficient $\tilde{\kappa}_{29}\gtrsim 1$ is required to match the observations. The same $\tilde{\kappa}$ is required to reproduce the canonical Milky Way constraints – most recent studies agree that for a MW with a diffuse gaseous halo of scale length $\gtrsim 10\,$kpc (appropriate for the simulations here, since the galaxies have extended halos and the diffusivity is, by assumption, constant) an effective, isotropically-averaged diffusivity $\tilde{\kappa}_{29}\gtrsim 1$ is needed (Blasi & Amato, 2012; Vladimirov et al., 2012; Gaggero et al., 2015; Guo et al., 2016; Jóhannesson et al., 2016; Cummings et al., 2016; Korsmeier & Cuoco, 2016; Evoli et al., 2017; Amato & Blasi, 2018).131313For the MW, we note that the observations (e.g. secondary-to-primary ratios and the like) do not really constrain the “diffusion coefficient” or “residence time” (these are model-dependent inferences), but rather the effective column density or “grammage” $X_{s}\equiv\int_{\rm CR\,path}\,n_{\rm gas}\,d\ell_{\rm CR}=\int_{\rm CR\,path}\,n_{\rm gas}\,c\,dt$ integrated over the path of individual CRs from their source locations to the Earth, with $X_{s}\approx 3\times 10^{24}\,{\rm cm^{-2}}$ measured. Repeating our calculation above for either an isothermal sphere or thin-disk gas distribution, it is straightforward to show that the grammage (integrated to infinity, as opposed to the solar circle) is directly related to the hadronic losses as $\dot{E}_{\rm loss}/\dot{E}_{\rm cr}=X_{s}^{\infty}/(2\,t_{0}\,n_{0}\,c)\sim 0.01\,(X_{s}^{\infty}/3\times 10^{24}\,{\rm cm^{-2}})$. Thus the constraints from matching the direct MW observations are essentially equivalent to matching the observationally “inferred” $F_{\gamma}/F_{\rm SF}\sim(0.3-1)\times 10^{-5}$ in the MW.

Second, we see why low diffusion coefficients cannot be invoked to increase the CR pressure (as noted above) – if one lowers $\tilde{\kappa}_{29}\ll 1$, then not only will the model fail to reproduce the observations, but the collisional losses will quickly dominate. If $\dot{E}_{\rm loss}\gtrsim\dot{E}_{\rm cr}$, it means that all CR energy is rapidly lost in the ISM, so there is no steady-state, high-pressure CR halo (which requires they escape the galaxy in the first place). In small dwarfs with low densities (low $\Sigma_{\rm gas}\,R_{\rm gas}$) one may be able to make $\tilde{\kappa}$ slightly lower without losing all of the CR energy, but we show below this still leads to large violations of the observational constraints.

Third, this means we do not expect CRs to (at least locally) have strong effects in extreme starburst-type galaxies, at either low or high redshift (where they are more common), owing to the very high $\Sigma_{\rm gas}$ observed in such systems, which should lead to efficient losses (see Chan et al. 2018 for more discussion).

3.2 Streaming & The Critical Radius

Now consider the streaming terms explicitly. If the streaming velocity is $\sim v_{A}\sim\beta^{-1/2}\,c_{s}$, and $\beta$ is approximately constant (as one might expect from e.g. a transsonic turbulent dynamo), then in an isothermal halo $v_{A}$ is also constant, so the radial streaming “flux” is $\sim v_{A}\,e_{\rm cr}$. The diffusive flux is $\sim\tilde{\kappa}\,\nabla e_{\rm cr}\sim\tilde{\kappa}\,e_{\rm cr}/r$, so the streaming dominates at $r\gtrsim r_{\rm stream}\sim\tilde{\kappa}/v_{A}$ (note we could have derived this instead in terms of the transport timescales to reach a given radius, and would reach the identical conclusion). This gives:

[TABLE]

where $v_{A,\,10}\equiv v_{A}/10\,{\rm km\,s^{-1}}$, $\beta_{100}\equiv\beta/100$, and $M_{{\rm vir},\,12}\equiv M_{\rm halo}/10^{12}\,M_{\odot}$.

Beyond this radius, if we continue to neglect losses, the solution would become that for a steady-state wind with constant velocity, i.e. $e_{\rm cr}(r>r_{\rm stream})\sim{\dot{E}_{\rm cr}}/({4\pi\,v_{A}\,r^{2}})$, so it falls more steeply compared to the diffusion-dominated case ($e_{\rm cr}\propto r^{-1}$). This has several important consequences. First (1) if trans-Alfvénic streaming dominates transport, then the streaming transport timescale $\sim r/v_{\rm stream}$ is always of the same order as the streaming loss timescale $\sim e_{\rm cr}/\dot{e}^{\rm stream}_{\rm loss}\sim e_{\rm cr}/|{\bf v}_{\rm stream}\cdot\nabla P_{\rm cr}|\sim(3/2)\,r/v_{\rm stream}$. So the CRs will, by definition, lose energy to Alfvénic and ultimately thermal energy as they stream at the same rate they propagate out, at this radius. This means the energy density must eventually decay, further accelerating streaming losses. So (2) a non-negligible fraction of the CR energy is thermalized within a factor of a few of this radius. And (3) since the CR pressure drops more rapidly, the CRs eventually provide small pressure support at $r\gg r_{\rm stream}$, even if they dominate over gravitational pressure at $r\lesssim r_{\rm stream}$ according to our arguments above.

Thus, in general, the dominant role of CRs is predicted to be confined to an “inner CGM” CR “halo” at $\lesssim 30-100\,{\rm kpc}$. For fixed $\beta$ and $\tilde{\kappa}$, lower $v_{A}$ means that this radius extends further in smaller halos, but (as we argued above), the CR injection and energy density also declines rapidly in smaller halos, so this simply means that the (relatively small) CR energy density more efficiently escapes, rather than being thermalized, in smaller halos. In massive halos, $r_{\rm stream}$ becomes comparable to the halo scale-radius.

Streaming losses still occur in the inner parts of the halo, even if streaming does not dominate the transport. Integrating the streaming loss rate $={\bf v}_{\rm stream}\cdot\nabla P_{\rm cr}$ within a maximum radius $r$, we obtain:

[TABLE]

so this is always a small (but not negligible) fraction of $\dot{E}_{\rm cr}$ within $r<r_{\rm stream}$. However, from the scaling of $r_{\rm stream}$, we see that streaming imposes yet another reason why CRs will not have a large effect at very low diffusion coefficients. If $\tilde{\kappa}$ is too low, then even if the gas density is so low that collisional losses are negligible, the CRs will lose much of their energy to Alfvénic damping (“streaming losses”) at very small $r$. Effectively, if $r_{\rm stream}$ is smaller than the gas disk effective radius, it means that the CR energy is thermalized before it can escape the star-forming disk.

Note that if we make $v_{\rm stream}$ much larger and allow the “streaming loss” term to increase with $\sim v_{\rm stream}\nabla P_{\rm cr}$ (instead of limiting it at $\sim v_{A}\nabla P_{\rm cr}$), then although streaming moves CRs faster, it also means streaming losses occur much faster and closer to the galaxy, where they can be radiated away more efficiently. In Eq. 11, for example, if we set $v_{\rm stream}^{2}=v_{A}^{2}+c_{s}^{2}$, then for $\beta\gg 1$ the expression becomes $r_{\rm stream}\sim 3\,{\rm kpc}\,\tilde{\kappa}_{29}\,M_{{\rm vir},\,12}^{-1/3}\,(1+z)^{-1/2}$. Thus, especially if combined with a lower diffusion coefficient, this particular super-Alfvénic streaming model means that streaming losses would thermalize most of the CR energy within the galaxy and ISM, where it would be efficiently radiated away.

Note that because $v_{A}$ and $r_{\rm stream}$ depend directly on $|{\bf B}|$ or $\beta$, if magnetic fields are stronger (weaker) it shifts our predictions for CR transport accordingly. For CRs, this is effectively degenerate with how we treat the scaling of $v_{\rm stream}$ with $v_{A}$: our default model ($v_{\rm st}\approx 3\,v_{A}$) is akin to a model with $v_{\rm st}\approx v_{A}$ but $10$ times lower $\beta$. Our experiments with different $v_{\rm st}$ therefore give some insight into how our predictions depend on the ultimate strength of magnetic fields.

3.3 What about CR Heating?

We have argued for the importance of the CR pressure/adiabatic terms above, and discussed CR losses. But can the CRs also be important as a thermal heating mechanism?

If we assume $100\%$ of the CR energy is thermalized (of course an upper limit, since some CRs escape, some energy is lost doing adiabatic work, and for hadronic losses $5/6$ of the energy goes into products like $\gamma$-rays which escape rather than thermalize), we can compare this to the cooling luminosity $L_{\rm cool}^{\rm tot}=\int\,\Lambda\,n^{2}\,d^{3}{\bf x}$, where $\Lambda\sim\Lambda_{22}\,10^{-22}\,{\rm erg\,s^{-1}\,cm^{3}}$ is roughly constant over the temperature range of interest. Using this and the various scalings above (for fully-ionized gas), we obtain:

[TABLE]

where $\langle n\rangle_{\rm cool}$ is the cooling-luminosity-weighted density (i.e. density where most of the cooling occurs). We immediately see, again, that CR heating cannot have a large effect in dwarf galaxies, owing to their very low $M_{\ast}/M_{\rm halo}$ (and correspondingly low SFRs and CR energy production), even if the CRs couple at very low CGM densities. In massive halos, at ISM densities ($\langle n\rangle_{\rm cool}\sim 1$), we see that there is no possible way CR heating can compete with radiative cooling. Moreover, the CR energy injection rate is an order-of-magnitude smaller than that from mechanical energy in SNe shocks.141414CR heating can be non-negligible in cold gas with $T\ll 10^{4}\,$K, where $\Lambda$ is much smaller, and the gas is strongly self-shielded so photo-electric heating is negligible, provided the local CR energy density is high. However this has a negligible effect on the dynamics of the cold gas (it is mostly important for accurate ionization fraction calculations), or on the total cooling budget of the ISM which is dominated by warmer gas cooling down to these low temperatures.

In the very lowest-density CGM near the virial radius (around the mean gas density of the halo, $\langle n\rangle_{\rm cool}\sim 10^{-4}\,(1+z)^{3}\,{\rm cm^{-3}}$), Eq. 13 suggests CR heating could become significant, but that assumes all the CR energy couples in the least-dense gas just inside $R_{\rm vir}$ (and ignores gas outside/inside). More accurately, if we assume streaming losses dominate (with the equilibrium $e_{\rm cr}$ for streaming at large $r>r_{\rm stream}$, defined above),151515In Eq. 14, because we assume streaming dominates the transport ($r>r_{\rm stream}$), the steady-state CR energy $e_{\rm cr}\propto\dot{E}_{\rm cr}/(v_{\rm stream}\,r^{2})$ while the streaming losses scale as $\sim v_{A}\,\nabla P_{\rm cr}\propto(v_{A}/v_{\rm stream})\,\dot{E}_{\rm cr}/r^{3}$. For $v_{\rm stream}\sim v_{A}$ (assumed in Eq. 14), the dependence on $v_{A}$ vanishes, while for super-Alfvénic streaming, the heating rate from streaming losses is reduced by a factor $\sim v_{A}/v_{\rm stream}$. Likewise, for streaming inside of $r<r_{\rm stream}$ where diffusion dominates transport, the heating rate is reduced by a factor $\sim v_{A}\,r/\kappa\sim r/r_{\rm stream}<1$. and that these are instantly thermalized, and that the gas is in an isothermal sphere, we can then calculate the ratio of the local thermal heating rate from CRs to the cooling rate:

[TABLE]

So the CR thermal heating is unlikely to be relevant at any radius (at least for the halos of interest here).

3.4 (Lack of) Effects Interior to Star-Forming Galactic Disks

The model above immediately implies that CRs have very weak effects within the star-forming galaxy disk at any mass scale. In order to avoid losing all the energy to collisional losses, it is required that the CR diffusion time ($\sim 0.03\,{\rm Myr}\,(L/100\,{\rm pc})^{2}\,\tilde{\kappa}_{29}$ for diffusion on scale $L$) is much faster than dynamical times in the disk. This essentially means CRs “diffuse out” of any locally dynamically interesting region of the disk (e.g. a GMC or strong shock) well before they can do interesting adiabatic work on that region. As a result, the CRs form (as assumed here) a quasi-spherical profile.

For a MW-like galaxy, our Eq. 6 predicts a CR energy density at the solar circle of $e_{\rm cr}\sim 1\,{\rm eV\,cm^{-3}}\,(\dot{M}_{\ast}/{\rm M_{\odot}\,yr^{-1}})\,\epsilon_{{\rm cr},\,0.1}\,(r/8\,{\rm kpc})^{-1}\,\tilde{\kappa}_{29}^{-1}$, more or less exactly the canonical value, and in order-of-magnitude equipartition with other disk-midplane energy densities. However, because of rapid diffusion, the CR pressure gradients are necessarily weak in the disk. If we assume a vertically-exponential gas disk balancing gravity in vertical hydrostatic equilibrium (disk midplane pressure $P_{\rm mid}\approx\pi\,G\,\Sigma_{\rm gas}^{2}$), then the ratio of the vertical pressure gradients $|\partial P_{\rm cr}/\partial z|/|\partial P_{\rm mid}^{\rm other}/\partial z|$ (where $\partial P_{\rm mid}^{\rm other}/\partial z$ is the gradient of thermal/magnetic/turbulent pressure in vertical hydrostatic equilibrium) scales as $\sim(P_{\rm cr}/P_{\rm mid}^{\rm other})\,(H/R)^{2}$ – i.e. the CR pressure forces (or gradients) are sub-dominant by a factor of $\sim(R/H)^{2}\gg 100$ in the MW midplane. The difference is even more dramatic if we consider still smaller-scale sub-structure (e.g. turbulent substructure in the ISM or clumps/cores in GMCs, where the relevant turbulent or magnetic/thermal support terms have structure on sub-pc scales).

3.5 Summary: The “Sweet Spot” for CRs

This toy model illustrates that, although for realistic (or observationally allowed) parameters we do not expect CRs to be dynamically dominant in the evolution of small dwarf galaxies, there is a potential “sweet spot” of galaxy parameter space in which CRs (from SNe) might be quite important for intermediate (LMC) through massive (MW-mass) galaxies ($M_{\rm halo}\sim 10^{11-12}\,M_{\sun}$) at low-to-intermediate redshifts ($z\lesssim 1-2$), via the creation of an extended CR halo in the inner CGM with pressure sufficient to maintain virial equilibrium and therefore support gas which might otherwise accrete onto the galaxy.

However even this requires some “sweet spot” in the cosmic ray transport parameter space. If $\tilde{\kappa}$ is too low, the CRs are trapped and collisional+streaming losses dissipate all their energy rapidly (in contradiction to all present observational constraints for Local Group galaxies). If $\tilde{\kappa}$ is too high (which may be observationally allowed), or is not constant but rises very rapidly outside of the galaxy (certainly allowed observationally), CRs will simply “free stream” out of the CGM without building up a significant pressure gradient or thermalizing their energy – although this may require extremely high $\tilde{\kappa}$. If the CR injection fraction $\epsilon_{\rm cr}\ll 0.1$, there is simply not enough energy in CRs to have an effect at any mass scale, while if it is too large ($\epsilon_{\rm cr}\sim 1$) it would violate direct observational constraints.

That is not to say this “sweet spot” requires very special fine-tuning. In fact, the characteristic parameters in the “sweet spot”: $\tilde{\kappa}\sim 10^{29}-10^{30}\,{\rm cm^{2}\,s^{-1}}$, $v_{\rm stream}\sim v_{A}$, $\epsilon_{\rm cr}\sim 0.1$ are both theoretically predicted (or at least plausible) and observationally allowed. These parameters are also in good agreement with the ”sweet spot” values identified in idealized galaxy simulations in Chan et al. (2018). And we stress that each of these can be varied by a factor of several (as we have done) without radically altering our conclusions. However, caution is needed, as these remain deeply uncertain.

4 Results

We now present the results of our simulations. Across all the parameter surveys described below, we study $>150$ high-resolution full-physics cosmological zoom-in simulations run to $z=0$. We discuss the implications of these results, and use them to test our simple model expectations, in § 5 below.

4.1 Case Studies

Figs. 1, 2, & 3, summarize a number of basic properties of the simulated galaxies: the “archeological” star formation history (distribution of star formation times); galaxy stellar mass and metallicity versus redshift; dark matter and baryonic mass profiles, stellar effective radii, and circular velocity curves at $z=0$. A number of other properties can be directly inferred from this (or are trivially related to information here), including e.g.: the “burstiness” of star formation, the distribution of stellar ages, the evolution of the mass-metallicity relation and stellar mass-halo mass relation, the existence of “cusps” or “cores” in the dark matter profile, baryon content of the halo, baryonic mass concentration, etc. We explicitly note the $z=0$ stellar mass, stellar effective radius, and metallicity. These figures intentionally parallel our previous comparison of feedback physics (adding or removing SNe, radiation, etc.), numerical methods and resolution in Hopkins et al. (2018b) – additional details of how each quantity is measured are given in that paper.

We show this for each of the “primary” (most well-resolved) galaxies in the high-resolution region of a number of representative simulations from Table 1. For each, we compare the different physics variations in Table 3: (1) our default FIRE-2 physics “Hydro+” run, (2) the FIRE-2 + MHD + anisotropic conduction and viscosity “MHD+” run, (3) the FIRE-2 + MHD + anisotropic conduction and viscosity + full CR physics “CR+” run with a “lower” parallel diffusion coefficient $\kappa\equiv 3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$, and (4) a “CR+” run with higher parallel diffusion coefficient $\kappa\equiv 3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$.

Simulations from Table 1 not shown explicitly in Figs. 1-3, as well as the (many) less-massive galaxies in each box, are omitted for brevity, but exhibit very similar behavior to those shown at similar masses. Likewise, Appendix A shows additional simulations with varied resolution: the qualitative conclusions are similar to the survey in Figs. 1-3.

Qualitatively, we see that “MHD+” and “Hydro+” runs are very similar in all respects and “CR+($\kappa=3e28$)” runs generally differ by only a small amount (with a couple exceptions). However “CR+($\kappa=3e29$)” runs produce suppressed SFRs and stellar masses in massive ($M_{\rm halo}\gtrsim 10^{11}\,M_{\sun}$) systems at low redshifts ($z\lesssim 1-2$).

4.2 $\gamma$-Ray Emission

Fig. 4 compares the $\gamma$-ray emission predicted in our simulations to observational constraints. With the exception of the MW, where more detailed constraints from spallation and other measurements exist (see § 3), $\gamma$-ray emission represents one of the most direct observational constraints available on the CR energy density in nearby galaxies. This was studied in detail in Chan et al. (2018), in non-cosmological simulations, so we extend that comparison here with the addition of our cosmological runs. Many of the $\gamma$-ray observations (and the equivalent constraint for the MW, essentially the measured grammage or “residence time”) are collected in Lacki et al. (2011) (note that more recent studies, e.g. Tang et al. 2014; Griffin et al. 2016; Fu et al. 2017; Wojaczyński & Niedźwiecki 2017; Wang & Fields 2018; Lopez et al. 2018, find consistent results). We compare directly to these constraints in Fig. 4, mimicking the observations as best as possible (see Chan et al. 2018 for additional details). Specifically, we compute the predicted luminosity $L_{\gamma}$ in $\sim$ GeV $\gamma$-rays produced by hadronic CR collisions. To do so, we follow Guo & Oh (2008); Chan et al. (2018) and take the exact hadronic loss rate computed in-code, assume $5/6$ of the losses go to pions with a branching ratio of $1/3$ to $\pi^{0}$, which decay to $\gamma$-rays with a spectrum that gives $\sim 70\%$ of the energy at $>1\,$GeV, and integrate this inside a $\sim 5\,$kpc aperture for dwarfs ($M_{\rm halo}<3\times 10^{11}\,M_{\sun}$) or $\sim 10\,$kpc aperture for more massive galaxies (similar to the effective areas used for observations, although this has a relatively small effect). We also compute the central161616“Central” radius is defined here following Lacki et al. (2011) as a projected radius $\sim 2\,$kpc for dwarfs, and $\sim 4\,$kpc for $L_{\ast}$ galaxies. However using a constant $\sim 1\,$kpc or half the effective radius gives qualitatively similar results. sightline-averaged gas surface density $\Sigma_{\rm central}$, and the luminosity from young stars $L_{\rm SF}$ (computed with STARBURST99 convolving all star particles $<100\,$Myr old with their appropriate ages and metallicities). This defines the ratio $L_{\gamma}/L_{\rm SF}$, versus $\Sigma_{\rm central}$, as measured in Lacki et al. (2011).

If all of the CR energy injected by SNe were lost collisionally, in steady-state with a time-constant SFR, this would produce a steady-state $L_{\gamma}/L_{\rm SF}\sim 2\times 10^{-4}$, which we label as the “calorimetric limit.” Of course, galaxies can violate this in transient events (or by a small systematic amount if e.g. star formation is non-constant in time). The ratio of $L_{\gamma}/L_{\rm SF}$ to calorimetric gives, approximately, the fraction of CR energy lost collisionally. We compare this for our “default” suite from Figs. 1-3, as well as a detailed surveys varying $\kappa$ (discussed below) and the streaming speed. For the streaming-speed study, we specifically re-run simulations m10q, m11b, m11q, m11d, m11h, m11f, m11g, m12f, m12i, and m12m, all with $\kappa=3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$, either with our default streaming speed, or with an arbitrarily much larger speed ($=3\,(c_{s}^{2}+v_{A}^{2})^{1/2}$, chosen ad hoc to be faster by a factor of a few than the largest MHD wavespeed).

We see that increasing $\kappa$ leads to more-efficient CR escape from the dense galactic gas, lowering $L_{\gamma}$ (as expected). Increased streaming speed produces a similar but much weaker effect, for the values we consider. The observations appear to strongly rule out $\kappa\ll 10^{29}\,{\rm cm^{2}\,s^{-1}}$: reproducing them requires $\kappa\sim 3-30\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$. In this regime, galaxies with dense nuclei (e.g. starbursts and bulge centers) are approximate proton calorimeters, but less-extreme systems (essentially all dwarfs and much of the volume of $\sim L_{\ast}$ galaxies) see the large majority ($>90\%$) of the CR energy escape the ISM without producing $\gamma$-rays.

4.3 Magnetic and CR Energies, Pressures, and Heating Rates

Fig. 5 shows the gas thermal ($=n\,k_{B}\,T$), magnetic ($=|{\bf B}|^{2}/8\pi$), CR ($=(\gamma_{\rm cr}-1)\,e_{\rm cr}$), and kinetic ($=\rho\,|{\bf v}|^{2}/2$) pressures (equivalent to energy densities), as a function of radial distance $r$ from the galaxy center at $z=0$ (in our standard physics suite). We focus on a representative sub-set of our ensemble, but the qualitative results shown here are generic to the ensemble of simulations. Each property is computed for each resolution element, exactly as they are determined in-code (self-consistently), and we show the $5-95\%$ inclusion contour of these properties at each galacto-centric radius. Because of our Lagrangian numerical method, this is effectively a gas-mass-weighted distribution of these values. To contrast, we therefore also plot the volume-averaged values in spherical shells. Note that where most of the gas is in a thin disk (e.g. in m12i at $r\lesssim 10\,$kpc), this means the spherically-volume-weighted average value of certain quantities (if they are concentrated in the thin disk) will be lower than its midplane value (closer to the mass-weighted average value) by a factor $\sim H/R$ (the ratio of the disk scale-height $H$ to radius $R$). We compare these to the “gravitational” or approximate local virial energy density, $\equiv\rho\,V_{c}^{2}/2\equiv\rho\,G\,M_{\rm enc}(<r)/2\,r$. Using a different definition based on the total potential gives quite similar results for our purposes here.

Note that the kinetic energy density here is defined as $\rho\,|{\bf v}|^{2}$, with ${\bf v}$ defined relative to the mean velocity of the whole galaxy, i.e. this includes rotation. This is done in part for simplicity because the dwarfs do not have strong rotation and separating rotation vs. dispersion (even in simulations) is in general quite challenging (see El-Badry et al., 2018b). It also allows us to immediately see (whether primarily in rotation or dispersion) where the gas is primarily “held up” by kinetic energy.

We see (discussed below) that magnetic pressure is almost always sub-dominant ($\beta\equiv P_{\rm thermal}/P_{\rm magnetic}\gg 1$), while CR pressure can dominate and maintain virial equilibrium in massive halos, especially in our high-$\kappa$ runs, at radii outside the galaxy from $\sim 30-200\,$kpc.

Fig. 6 considers a similar comparison of the radial profiles of the gas cooling rate, computed in-code at $z=0$ in each radial annulus using our full cooling function, to the heating rate from the CR “streaming loss” (gyro-resonant Alfvén-wave heating) and “collisional” (thermalized hadronic+Compton) terms. These are always sub-dominant to the gas cooling rates (dominated by gas at $\sim 10^{4}-10^{6}\,$K, near the peak of the cooling curve).

4.4 Internal Galaxy Properties: Metallicities, Star Formation Rates, ISM Phase Structure, Morphologies, Angular Momentum

We next survey a more detailed set of internal galaxy properties.

Fig. 7 shows each of the “primary” galaxies in our suite (Table 1) on the stellar mass-halo mass and stellar mass-stellar metallicity relations at $z=0$. These re-affirm the trends suggested in Figs. 1-3: our “MHD+” and “CR+($\kappa=3e28$)” runs do not deviate significantly from “Hydro+,” while the “CR+($\kappa=3e29$)” runs show systematically lower stellar mass in massive halos ($M_{\rm halo}\gtrsim 10^{10.5-11}\,M_{\odot}$) by a factor $\sim 2-3$, but all runs move along a single (relatively tight) mass-metallicity relation. In addition to the “primary” galaxies, we randomly select $\sim 50$ additional non-satellite galaxies within the box171717These galaxies are chosen within the high-resolution “zoom-in” region ($<1\%$ contamination by low-resolution dark matter particles inside their virial radii), but $>500$ kpc away from the “primary” and outside the virial radius of any more massive halo. We select from the boxes m12i, m12f, m12m (each with resolution $m_{i,\,1000}=7$) because these have relatively large zoom-in regions. surrounding our m12 halos, to compare with the “primary” dwarfs. Although the dwarfs in the larger m12 boxes are simulated at lower resolution ($m_{i,\,1000}=7$) compared to the smaller m09 or m10 boxes around a single dwarf (reaching $m_{i,\,1000}=0.25$), the results are consistent.

In Fig. 7 we also include observational estimates for guidance, but stress that we are not attempting to compare rigorously (e.g. we are not matching what is actually measured, comparing scatter, etc.) – such comparisons are presented (for non-CR runs) in Ma et al. (2016); Hopkins et al. (2018b); Ma et al. (2018a); Garrison-Kimmel et al. (2018); Wheeler et al. (2018). All our galaxies lie within the approximate $\sim 5-95\%$ scatter inferred around the median $M_{\ast}-M_{\rm halo}$ relation at all masses shown (Garrison-Kimmel et al., 2017a), although interestingly at the highest masses ($M_{\rm halo}\gtrsim 10^{11.5}\,M_{\odot}$) the low-$\kappa$ CR and Hydro/MHD runs appear to be systematically high in $M_{\ast}$ while the high-$\kappa$ CR runs are systematically low (bracketing the observed median). In metallicity the simulations agree extremely well with observations at $M_{\ast}\gtrsim 10^{8}\,M_{\odot}$ but appear to fall more steeply in very faint dwarfs compared to the Local Group satellites in Kirby et al. (2013): this is explored in detail in Wheeler et al. (2018).

Fig. 8 shows the locations of the example galaxies from Figs. 1-3 on the observed Schmidt-Kennicutt relation, at $z=0$. We specifically follow Kennicutt (1998) and measure the SFR (defined as the average over the last $<10\,$Myr) and total neutral (HI+H2) gas mass inside a projected circular aperture enclosing $90\%$ of the starlight (approximately equivalent to the RC2 sizes used therein, for MW-mass systems, but this allows us to include low-surface-brightness dwarfs), to define $\Sigma_{\rm SFR}$ and $\Sigma_{\rm gas}$. We also show the $\sim 90\%$ inclusion contour of observed galaxies compiled from Kennicutt et al. (2007); Bigiel et al. (2008); Genzel et al. (2010), although this is again intended only for reference, as we are not attempting a rigorous comparison with observations (such a comparison is presented in Orr et al. 2018, for those interested) but only to discern systematic effects between different runs. Briefly, we see no apparent offset in the relation between any of our physics suites.

Fig. 9 compares the distribution of gas phases (densities and temperatures) in the interstellar medium. We plot the distribution of mass as a function of density, in the different traditional ISM phases (molecular, cold neutral, warm neutral, warm ionized, hot ionized), where for simplicity we define “ISM” as gas within $r<10\,$kpc (the exact threshold makes little difference). For simplicity, we only compare “MHD+” and “CR+($\kappa=3e29$)”, as “Hydro+” and “MHD+” are nearly-identical (see Su et al. 2017 for more detailed comparison), and “CR+($\kappa=3e28$)” is as well. Even with higher-$\kappa$, the differences are very subtle (tens of percent shift in cold-to-warm-neutral medium mass). In a temperature-density diagram, this subtle shift is nearly undetectable. We see similar effects examining the phase structure of outflows, specifically, but defer a detailed study of this to future work.

Figs. 10, 11, 12, 13 shows the visual morphologies of the galaxies at $z=0$ in mock HST images of starlight, as well as the gas morphology in different phases and on different spatial scales. The stellar images are mock $u/g/r$ composite ray-tracing images determined using STARBURST99 to compute the age-and metallicity-dependent spectrum of each star (the same assumptions used in-code) and adopting a constant dust-to-metals ratio to use the gas and metals distribution determined in-code to attenuate and extinct the light; the gas images are volume-renderings with iso-temperature contours centered on broad (log-normal) temperature bands with dispersion $\sim 0.5$ dex around “cold,” “cool” or “warm,” and “hot” gas (see Hopkins et al. 2014 for details of the rendering). For dwarfs, we see no major difference with any physics change. For massive halos, the high-$\kappa$ CR runs are later-type, consistent with their lower mass, and we see more warm (as compared to hot) gas in the CGM (beyond the disk). The CGM properties will be studied in greater detail in future work.

We have also examined the morphology of the magnetic fields, specifically, but find they are highly tangled on all scales here, in both “MHD+” and “CR+” runs, consistent with our more detailed studies in Su et al. (2017) and Su et al. (2018a).

Fig. 14 shows the “quantitative kinematic morphology” (angular momentum distribution) of stars (weighted by stellar mass or visual luminosity) and neutral gas. In all of these, we again compare our ensemble of galaxies from Table 1 across our “core physics variations” from Table 3. Again, the variations with physics are weak, and (where present) consistent with the morphological changes described above.

4.5 Variations in CR Transport Physics

We have widely-varied the CR transport coefficients $\kappa$ and $v_{\rm stream}$ in a subset of our simulations. The effect on $\gamma$-ray emission is shown in Fig. 4.

Fig. 15 shows the effects on galaxy properties (using the same style as Fig. 1) of varying $\kappa$ more widely and densely in our galaxies m10q, m11b, m11q, m12i. This confirms that CRs have weak effects for $\kappa\lesssim 10^{29}\,{\rm cm^{2}\,s^{-1}}$ and maximal effect in massive halos at $z\lesssim 1-2$ for $\kappa\sim 3-30\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$.

Fig. 16 varies the fraction $\epsilon_{\rm CR}$ of SNe energy which is injected into CRs, focusing on a low-resolution version of run m12i with $\kappa=3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ (since MW-mass halos with this diffusivity are where we see the most dramatic CR effects). As expected, increasing the CR energy input increases their effect, but the effect is quite weak (sub-linear).181818Throughout, we have assumed CRs are injected at the sites of SNe: because the efficiency of CR acceleration scales with the shock velocity, models (and SN remnant observations) generally assign most of the acceleration to the “fastest” SNe shocks with $v\gtrsim 1000\,{\rm km\,s^{-1}}$ and “swept up” ISM mass $\sim M_{\rm ejecta}\sim 1-10\,M_{\odot}$, which is always un-resolved in our simulations. However one might imagine that if CRs are injected primarily outside the galaxy from e.g. structure-formation or superbubble-CGM shocks, then this would allow one to match the observations with lower $\kappa$ – but only if we dramatically lowered the injection from SNe themselves (otherwise this would simply increase $L_{\gamma}$). So this is unlikely to change our conclusions here, but could be important for CRs from AGN (where jet termination shocks reach into the CGM).

Fig. 17 compares m10q, m11q, m12i in a survey of basic CR physics: comparing our default “CR+” implementation (§ 2) to runs (1) without MHD (where lacking ${\bf B}$ directions we assume CR diffusion/streaming, and conduction/viscosity, are isotropic, i.e. take the projection tensor $\hat{\bf B}\otimes\hat{\bf B}\rightarrow\mathbb{I}$); (2) without collisional (hadronic or Coulomb) CR losses; (3) without streaming (setting $v_{\rm stream}\rightarrow 0$); and (4) allowing super-sonic and super-Alfvénic streaming by setting the streaming speed to a few times the fastest MHD wavespeed $v_{\rm stream}\rightarrow 3\,(c_{s}^{2}+v_{A}^{2})^{1/2}$. These effects are discussed in detail below but their effects are generally small compared to including/excluding CRs at all.

Note that for Fig. 4 we noted we have actually run a large ensemble of simulations with lower $\kappa=3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$ and the larger $v_{\rm stream}\rightarrow 3\,(c_{s}^{2}+v_{A}^{2})^{1/2}$: we have also compared the resulting galaxy properties but omit them for brevity as the difference owing to faster streaming is very small. We have also run a large suite with $\kappa=3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ varying $v_{\rm stream}=(1-3)\,v_{A}$ more modestly, and find (unsurprisingly) essentially no effect.

Additional purely-numerical tests are given in Chan et al. (2018) and Appendix A. This includes modifications of the CR pressure tensor, form of the flux equation, and other detailed assumptions that stem from any two-moment expansion of the CR transport equations. We do not see large effects from these on our predictions, but of course can only survey a limited range of possibilities.

4.6 High-Redshift Galaxies

Figs. 18-19 repeats our earlier exercises comparing galaxy properties and magnetic+CR energy densities (from Figs. 1 and 5) for the suite of halos in Table 2, which reach large masses ($M_{\rm halo}>10^{12}\,M_{\sun}$) at increasingly higher redshifts $z\sim 1-10$. These simulations are not run to $z=0$, hence we analyze them separately.

We discuss the results in detail below, but briefly, we see no appreciable effects of CRs or MHD above $z\gtrsim 1-2$.

5 Discussion

We now explore the implications of the results presented in § 4.

5.1 Magnetic Fields, Conduction, & Viscosity

In Su et al. (2017), we used similar FIRE-2 simulations to study the effects of magnetic fields and anisotropic Spitzer-Braginskii conduction and viscosity on galaxy properties. The study there included much more detailed measurements of properties like the gas phase distributions of the ISM and CGM, outflow properties, turbulence and energy balance in the ISM, magnetic field amplification, and more. There, we concluded that there was no appreciable systematic effect on any global galaxy properties from these physics.

The “MHD+” simulations here improve on those studied in Su et al. (2017) in three significant ways. (1) Our mass resolution is an order-of-magnitude better, allowing us to much better-resolve the Field length and other small-scale effects. (2) The simulation suite here is an order-of-magnitude larger, and all fully-cosmological, allowing us to assess and improve the statistics and avoid uncertainties owing to inevitable run-to-run stochastic variations in galaxy properties (discussed therein or in Keller et al. 2018; Genel et al. 2018). (3) Our treatment of anisotropic conduction & viscosity is more accurate. Specifically, a large body of recent work in the plasma physics literature (both theoretical and experimental) has shown that the parallel transport coefficients for heat and momentum (e.g. $\kappa_{\rm cond}$ and $\eta_{\rm visc}$) are strongly self-limited by micro-scale plasma instabilities (e.g. the Whistler, firehose, and mirror instabilities) in high-$\beta$ plasmas like the ISM and CGM (Kunz et al., 2014; Komarov et al., 2016, 2014; Riquelme et al., 2016; Santos-Lima et al., 2016; Roberg-Clark et al., 2016, 2018; Tong et al., 2018; Squire et al., 2017c, a, b; Komarov et al., 2018). Our treatments in § 2.2 include these effects, while our previous work did not.

Despite this (or perhaps because of it), we confirm the conclusions of Su et al. (2017): in every property we measure, the “MHD+” runs do not appear to systematically differ significantly from the “Hydro+” runs. This is perhaps not surprising: the only physical change of (1)-(3) above is (3), which has the effect of uniformly decreasing the magnitude of conduction and viscosity at high-$\beta$ (e.g. the ISM and CGM). Because of this, we will not discuss these variations further (for much more detailed discussion of why the results do not change, we refer to Su et al. 2017).

Our goal here is to understand effects of magnetic fields on galaxies, not to consider a detailed comparison with observations. However, Guszejnov et al. (2019) directly compared the values of ${\bf B}$ predicted in our MW-like m12i simulation in atomic/molecular clouds to Zeeman observations from Crutcher et al. (2010) and found remarkably good agreement at all observed densities ($n\gtrsim 10\,{\rm cm^{-3}}$). In the ionized/warm ISM and thick disk/inner halo, observations are much more uncertain and model-dependent (see Han, 2017), e.g. assuming gas density profiles and equipartition between CR and magnetic energy to infer $|{\bf B}|$ from rotation measures (RMs) or synchrotron emission ($I_{s}$). But preliminary estimates from our simulations show the combination of large “clumping factors” (e.g. $\langle|{\bf B}|^{4}\rangle/\langle|{\bf B}|\rangle^{4}$), violation of equipartition, and more extended gaseous halos can produce similar RM and $I_{s}$ to observations, often with much lower median ${\bf B}$ compared to simpler models. More detailed forward-modeling is clearly needed. In the extended CGM/halo, observations are lacking, but comparing our Fig. 5 to other cosmological simulations of MW and dwarf galaxies including stellar feedback shows that we predict similar, or somewhat higher, values of $|{\bf B}|$ as a function of either baryon density ($n/\langle n\rangle$) or galacto-centric distance $r/r_{\rm vir}$ (Vazza et al., 2014; Marinacci et al., 2015; Martin-Alvarez et al., 2018; Katz et al., 2019). Older simulations not including stellar feedback show substantially lower $|{\bf B}|$ (Dubois & Teyssier, 2008; Dolag et al., 2008; Donnert et al., 2018), consistent with recent studies of the role of turbulence and winds in field amplification (Su et al., 2018a; Martin-Alvarez et al., 2018).

5.2 Cosmic Rays in Dwarf Galaxies

In essentially every galaxy property we examine, the effects of CR physics on dwarf galaxies ($M_{\ast}\ll 10^{10}\,M_{\odot}$, $M_{\rm halo}\lesssim 10^{11}\,M_{\odot}$) are relatively small. Above in § 3, we argued that this is generically expected: in equilibrium for realistic star-forming galaxies on the observed “main sequence” of star formation, the ratio of CR pressure forces on gas in the halo, compared to gravitational forces, scales $\propto M_{\ast}/M_{\rm halo}$. Of course, other forms of feedback scale this way as well. But with other forms of stellar feedback in place (notably SNe), CR pressure is relatively inefficient at re-accelerating winds or stalling accretion on scales of order the galaxy and halo scale radius. One key point is that on small (ISM) scales, in dwarfs, with low metallicities and densities, SNe cool relatively inefficiently and can convert a large fraction of their ejecta energy into work and/or thermal energy which is not immediately radiated (see Hopkins et al. 2018c). Since the CR energy is only $\sim 10\%$ of the mechanical energy, by construction, this means mechanical energy will always dominate, if it is not efficiently radiated (as it would be in dense gas, in more massive galaxies).

A second key point follows directly from the observations. In dwarf galaxies, almost all of the CRs escape efficiently from the ISM: but unlike massive galaxies where there is a large CGM “hot halo” extending to $>300\,$kpc which can “confine” the CRs and allow them to build up a substantial pressure profile, in dwarf galaxies there is no such hot halo (it cannot be built up from virial shocks, and outflows are much less confined as well so simply escape) and the size of the virial radius is much smaller ($\ll 100\,$kpc). So CRs in dwarfs simply escape from their CGM as well (and even if they stay, they have relatively little mass to support and prevent from accreting).

We do see some effects of CRs, but these are small compared to the effects of different treatments of SNe and stellar radiation (see Hopkins et al. 2018b, c, a), and appear primarily only at lower $\kappa$ than allowed observationally. Table 1 and Figs. 1-3 (also 7, 15) show that at ultra-faint masses (e.g. m09 and m10v or m10q at very early times, with $M_{\ast}\lesssim 3\times 10^{5}\,M_{\odot}$), runs with CRs have slightly increased stellar mass by $\sim 0.1-0.2$ dex. At somewhat higher masses CRs can slightly suppress SF in dwarfs, for $\kappa$ within the range $\kappa\approx 10^{28}-10^{29}\,{\rm cm^{2}\,s^{-1}}$. This effect is maximal for m10q, but even there is only a factor of $\sim 2$ (much smaller than many radiative or mechanical feedback effects in such small galaxies, which are only turning $\ll 1\%$ of their baryons into stars). By slightly higher masses (m11b, m11q at $M_{\ast}\sim 10^{7.5}-10^{8.5}\,M_{\sun}$) the effect weakens to $\sim 0.1$ dex or nothing at all ($M_{\ast}\gtrsim 10^{9}\,M_{\sun}$ in m11v, m11c, see Table 1). Figs.  8, 9, 10 & 14, show that even where the effect on stellar masses is maximal it has no systematic qualitative effect on the visual/stellar or gas morphology, phase structure, or kinematics either within the galaxy itself or the inner CGM, nor on the star formation efficiency of the galaxy.

Figs. 5-6 illustrate why this is the case, in the context of the arguments above. In all cases, CR “heating” is vastly sub-dominant to gas cooling (as expected – recall the virial temperatures of these halos are near the peak of the cooling curve). Moreover, for higher-$\kappa$ ($\gtrsim 10^{29}\,{\rm cm^{2}\,s^{-1}}$), even with zero losses, the CR pressure in the galaxy and CGM is always sub-dominant to thermal pressure (SFRs are simply too low to support an energetically-dominant CR halo). By making $\kappa$ lower and trapping CRs one can build up more CR pressure close to the galaxy (keeping warm gas “puffy” and non-star-forming), but for $\kappa\lesssim 10^{28}\,{\rm cm^{2}\,s^{-1}}$ diffusion is too slow and much of the CR energy is lost collisionally. Even at the “sweet spot” where the CR pressure is maximized (about $\kappa\sim 3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$), CR pressure only becomes comparable to thermal pressure (never strongly dominant) in the more massive dwarfs.

Because these effects are somewhat fine-tuned, and CRs are never strongly dominant, they are also sensitive to other details of the physical treatment: Fig. 17 shows that if we allow modestly super-Alfvénic streaming (increasing both the streaming speed and “streaming loss” term in these runs), the combination of enhanced losses and faster escape from the galaxy (extending but lowering the CR pressure) eliminates the (already small) effect of CRs within galaxies. Likewise, isotropic transport (effectively allowing slightly faster diffusion) allows the CRs to escape and weakens the effects. These conclusions are discussed and demonstrated in more detail in Chan et al. (2018) (and also consistent across our more extensive studies for e.g. Fig. 4).

Finally, and perhaps most important, Fig. 4 shows that the relatively low $\kappa$ required to produce an effect on dwarfs leads to a factor $\sim 10$ over-prediction of the observed $\gamma$-ray luminosity in low-surface-density galaxies. Recall, the observed points include the LMC, SMC, and M33, analogous to several of our simulated dwarfs. In order to match these, Fig. 4 shows $\kappa>10^{29}\,{\rm cm^{2}\,s^{-1}}$ is required, at which point the effect of CRs on most internal dwarf galaxy properties vanishes.

Briefly, we note that our conclusions here appear to contradict some claims in the literature that CRs can have a strong effect on SF in dwarf galaxies (e.g. Jubelgas et al., 2008; Booth et al., 2013; Chen et al., 2016). However, to our knowledge, all such claims either (a) adopted low diffusion coefficients, $\kappa\ll 10^{29}\,{\rm cm^{2}\,s^{-1}}$, without comparing to the constraints from $\gamma$-ray fluxes which we argue prohibit such coefficients (in several such studies, collisional losses were not included at all, which allows CRs to artificially “build up” at low $\kappa$ when their energy should be lost); or (b) used idealized (non-cosmological) simulations, or simulations with very weak (or non-existent) stellar feedback from other sources (e.g. mechanical SNe and/or radiative feedback), such that the SFR and ratio $M_{\ast}/M_{\rm halo}$ was much higher than observed (which, according to the scalings in § 3, would allow CRs to have a large effect, but directly violates observations of galaxy stellar masses and SFRs).

We do stress that CRs could still have an effect in the CGM, or ICM further away from the galaxy, at least at intermediate mass scales. This will be studied in detail in future work, but briefly, we note that even in the outer CGM or IGM out to $\sim 4\,R_{\rm vir}$, we see no obvious systematic effect of CRs in very low-mass halos (e.g. $M_{\rm halo}\lesssim 10^{10}\,M_{\sun}$, i.e. m10q and smaller galaxies). However we do see some effects on the velocity field of gas even at surprisingly small halo mass scales (down to $M_{\rm halo}\sim 4\times 10^{10}\,M_{\sun}$, corresponding to m11b). In these intermediate-mass systems, Fig. 5 shows the CR pressure is not completely negligible around $\sim R_{\rm vir}$ (although it is not dominant), so this is plausible. But since our primary focus here is galaxy properties, we defer a more detailed investigation to future work.

In summary, for any observationally-allowed CR parameters explored here, the effects on any galaxy property studied are small (much smaller than effects of e.g. mechanical or radiative feedback).

5.3 Cosmic Rays in Intermediate and Milky Way-Mass Galaxies

As we look at progressively more massive low-redshift galaxies, above $M_{\rm halo}\gtrsim 10^{11}\,M_{\odot}$ (at $z=0$), however, we see that CRs can have significant effects, as predicted in § 3.

However, with a lower diffusion coefficient, $\kappa\lesssim 10^{29}\,{\rm cm^{2}\,s^{-1}}$, the effects on galaxy properties (in e.g. Figs. 2-3, 7-9, 11-14, etc.) are very weak – typically $\sim 0.1\,$ dex or so in SFR and stellar mass, at most. The effects at CRs at low-$\kappa$ are even weaker with super-Alfvénic streaming, as seen with dwarfs (see also Chan et al. 2018). The reason is obvious in Fig. 4: for $\kappa\lesssim 10^{29}\,{\rm cm^{2}\,s^{-1}}$, the massive galaxies become proton calorimeters, i.e. lose most of their CR energy to collisions within the galaxy – also as predicted in § 3. Note that the massive galaxies have higher central densities compared to the dwarfs (almost all have $\Sigma_{\rm central}\gtrsim 10^{-2}\,{\rm g\,cm^{-2}}$), so it requires larger $\kappa$ for CRs to escape without losing most of their energy. This is also obvious in Fig. 5 – for these lower $\kappa$ values, the CR pressure/energy density outside the galaxy is order-of-magnitude below the predicted value if the CRs diffused without losses.

Fig. 15 shows that, as a result, the “sweet spot” where CRs have maximal effect occurs at $\kappa\sim 3-30\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$.191919Note that although the low-resolution ($m_{i,\,1000}=56$) m12i shown in Fig. 15 shows increasing effects of CRs going from $\kappa\sim 3\times 10^{29}$ to $\kappa\sim 3\times 10^{30}$, the low-resolution version of this simulation is a relatively high-density galaxy with $\Sigma_{\rm central}\sim 10^{-1.2}\,{\rm g\,cm^{-2}}$ in Fig. 4, so the higher $\kappa$ improves escape. We have run both higher resolution ($m_{i,\,1000}=7$) m12i (where the galaxy is systematically less dense), and m11f with $\kappa\sim 3\times 10^{30}$, to limited redshift ($z\sim 1$), and find in these less-dense ($\Sigma_{\rm central}\sim 10^{-(1.7-2.1)}\,{\rm g\,cm^{-2}}$) systems that the CR effects become weaker at this still-higher $\kappa$, as predicted if they escape “too” efficiently. For these $\kappa$, Figs. 2-3 show that as the the galaxies approach stellar masses $M_{\ast}\sim 10^{10}\,M_{\ast}$, or halo mass $\sim 10^{11}\,M_{\sun}$ (SFRs $\gtrsim 0.1-1\,{\rm M_{\sun}\,yr^{-1}}$) the CRs begin to have an effect suppressing SF. The suppressed SFRs in the “CR+($\kappa=3e29$)” runs tend to flatten once this mass scale is reached, while SFRs in the “Hydro+”, “MHD+”, and “CR+($\kappa=3e28$)” continue to rise, so at the peak of the “Hydro+” run SFRs ($z\sim 0$ for the lower-mass m11 runs, or $z\sim 1$ for the higher-mass m12 runs) the difference in SFR can be as large as factor $\sim 3-10$, although the integrated difference in stellar mass by $z\sim 0$ is usually a more modest factor $\sim 2-3$.

These effects do depend on redshift, as discussed below (§ 5.4) and shown in Figs. 18-19. Essentially all the effects we see from CRs are confined to relatively low redshifts $z\lesssim 1-2$.

Figs. 16-17 (see also Appendix A & the more detailed studies in Chan et al. 2018) show that the generic behaviors described above are not especially sensitive to other details of the CR transport physics, provided similar large $\kappa$ (discussed further in § 5.6). In Fig. 16 we systematically vary the fraction of SNe energy injected as CRs from $\epsilon_{\rm cr}\sim 0.05-0.2$: as expected, more efficient CR production produces stronger SFR suppression, but the effect is highly sub-linear (a factor of $\sim 4$ change in $\epsilon_{\rm cr}$ produces a factor $\sim 1.5$ change in stellar mass or SFR). Consistent with Fig. 5, the effect is primarily a “threshold” effect: once sufficient CRs reach large radii to set up a halo that can pressure-support the cool gas, adding somewhat more produces little effect. Of course eventually if $\epsilon_{\rm cr}$ is too low, the CR halo cannot maintain pressure support, and the effect of CRs should rapidly vanish.

Fig. 4 shows that for these MW-mass systems, the same $\kappa\sim 3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ which produces large effects on their SF histories appears to agree well with the spallation constraints in the MW and observed $\gamma$-ray luminosities in the MW and local galaxies, while the lower $\kappa\lesssim 10^{29}\,{\rm cm^{2}\,s^{-1}}$ (which produces weak effects on the galaxies) predicts excessive $\gamma$-ray flux.

The effects of CRs within the galaxies, even at high-$\kappa$, appear quite weak, as expected (see § 3) – properties like the stellar/visual and/or gas morphology and kinematics within the galactic disks, abundances/metallicities, baryonic and dark matter mass profiles, disk sizes and rotation curves, outflow rates, and star formation efficiencies (location on e.g. the Kennicutt-Schmidt relation) do not appear to be strongly altered at any $\kappa$ or mass scale we study, except insofar as the total galaxy mass and SFR shift up or down (changing e.g. the total gas supply, or total mass in metals produced, or total mass in baryons contributing to $V_{c}$). This is not surprising: in the disk midplane, we confirm the CR pressure gradients are sub-dominant to thermal and turbulent forces. But in the CGM around the galaxies, the CRs appear to have a direct and dramatic effect. We discuss this further in § 5.5 below.

However Figs. 5 and 11-13 demonstrate that the effects of CRs on the CGM around these galaxies, at radii $\sim 10-100\,$kpc, are dramatic. In the cases where CRs suppress SF, they establish a high-pressure CR halo outside of the galactic disk (extending to or even past the virial radius), supporting a large reservoir of gas which is much cooler ($T\ll 10^{6}$ K) than would be required to maintain thermal pressure equilibrium. In contrast in the “MHD+” or “Hydro+” cases denser and/or cooler halo gas cools rapidly, then falls onto the galaxy, leaving a virialized halo of only the “leftover” gas which is more tenuous. In the “CR+” runs around massive galaxies, Fig. 5 shows the CR pressure is dominant over thermal and magnetic pressure outside the disk, all the way to the virial radius (Fig. 6 shows the direct CR heating is negligible). For low-$\kappa$ the CR pressure is suppressed owing to losses, but for the high-$\kappa$ runs the CR pressure profile agrees remarkably well with the analytic predictions in § 3, which also predict accurately the mass scale where this can support enough gas mass to suppress gas inflows and (ultimately) SF in the galaxies at a significant level. We discuss the dynamics of the CRs in the CGM in more detail below and in future work.

5.4 Dependence on Redshift (and Super-$L_{\ast}$ Massive Galaxies)

The galaxies of interest in § 5.3, where CRs have appreciable effects, have $M_{\rm halo}\sim 10^{11-12}\,M_{\sun}$ at $z\sim 0$. It is natural to ask whether galaxies in this mass range at higher redshifts also experience similar effects of CR transport.

First, note that high-redshift dwarfs are represented in our sample already, in Figs. 1-3, etc. These are simply the progenitors of the $z=0$ galaxies. Since we plot their entire evolutionary history, it is straightforward to see in these figures that essentially all of the differences owing to CRs (at any mass scale) manifest only at relatively low redshifts $z\lesssim 1-2$. In most cases, it is quite notable: e.g. for galaxies m11i, m11d, m11h, m11f, m11g, m12i, m12f, m12m, the SFH, stellar mass, and metallicity (as well as all other properties, such as morphology, outflow properties, etc.) in the CR+ runs closely track the MHD+ runs until $z\lesssim 1-2$, where they begin to diverge significantly. For the lower-mass galaxies, this is less surprising: their progenitors at $z\gtrsim 1-2$ are small dwarfs (for which the effects of CRs are weak at $z=0$). But note that some of the more massive galaxies (e.g. m12m and m12f) reach stellar masses $\gtrsim 10^{10}\,M_{\sun}$ and halo masses $\gtrsim 10^{11.5}\,M_{\sun}$ long before they begin to diverge – well above the threshold where we see effects appear at $z=0$. This indicates that there is some additional redshift dependence here.

To explore more massive halos at higher redshifts, i.e. those with $M_{\rm halo}\sim 10^{11-12}\,M_{\sun}$ at $z\sim 1-10$, Figs. 18-19 consider the extended suite of high-redshift, high-mass galaxies from Table 2. These are halos which reach $\sim 10^{12}\,M_{\sun}$ (comparable to our most massive $z=0$ halos) already at various redshifts $z=z_{12}\sim 1-10$. The earliest of these represent progenitors of what will be massive galaxy clusters, by $z\sim 0$. In Fig. 18 we can see essentially no detectable systematic effects of MHD or CRs on these halos at any redshift, except for the lowest-redshift example (m12z1) which begins to show a modest suppression of its SFR and stellar mass only at $z\lesssim 1.5$. We have also examined the other properties in this paper (e.g. morphologies, gas phase distributions, outflow velocities) and similarly find no differences above $z\gtrsim 1-2$ (hence their not being shown here, for brevity).

This is actually predicted by the simple scalings in § 3. At high redshift, SFRs are higher at a given stellar mass (as obvious in Fig. 18), so the CR production/injection rate is also higher (scaling $\propto t_{\rm Hubble}^{-1}\sim(1+z)^{3/2}$). However, the density of the CGM and IGM, and ram pressure of inflowing gas, is much higher (scaling $\sim(1+z)^{3}$). So CRs are unable to maintain a super-virial pressure which can suppress inflow/accretion in this dense gas (as in Fig. 5). We have directly confirmed this, in fact, comparing the CR and virial pressure in Fig. 19 – for the high-redshift massive halos, most of the inflow is under-pressurized relative to what would be needed to maintain virial equilibrium, with or without CRs.202020In the high-redshift halos in Table 2, there is some low-density gas in the CGM for which the CR pressure exceeds or is comparable to that needed for virial equilibrium. This may imply that CRs have an effect in the CGM even where they have little effect on the bulk galaxy properties. However, this under-dense gas represents very little of the gas which accretes onto the galaxy (with or without CRs). Moreover, owing to their high inflow rates, most of the star formation in the high-redshift systems occurs in starbursts with very high SFRs and, correspondingly, very high gas densities within the galaxy (obeying the extension of the observed Schmidt-Kennicutt law in Fig. 8 to higher densities) – we see this directly in Fig. 18 where all the massive high-$z$ systems reach $\dot{M}_{\ast}\gtrsim 100\,M_{\sun}\,{\rm yr^{-1}}$ and the couple most extreme reach $\dot{M}_{\ast}\gtrsim 1000\,M_{\sun}\,{\rm yr^{-1}}$. During these phases, the gas surface densities reach $\gg 0.1\,{\rm g\,cm^{-2}}$ (see Fig. 8 and the more detailed studies in Hopkins et al. 2011; Orr et al. 2018), or $\gg 10^{3}\,M_{\sun}\,{\rm pc^{-2}}$. At these densities, we see in Fig. 4 that essentially all observed galaxies, and all of our simulations (even with extremely high $\tilde{\kappa}>10^{30}\,{\rm cm^{2}\,s^{-1}}$) approach the proton-calorimetric limit – in other words, a substantial fraction of the CR energy is lost to collisions in the dense ISM. We confirm this directly in the simulations in Fig. 18 in their “peak starburst” phases. Interestingly, there is some evidence that away from the starburst peaks, at the lowest SFRs, the CR+ runs exhibit slightly more-suppressed SFRs, but these phases of course contribute negligibly to the total SF and stellar mass (and therefore most other galaxy properties).

Importantly, these simulations also highlight the critical need for some additional feedback beyond the stellar feedback mechanisms (SNe types Ia and II, stellar mass loss from O/B and AGB stars, radiation, cosmic rays accelerated in SNe) and microphysical processes (magnetic fields, conduction, viscosity) studied here, in order to explain the suppression of SF (i.e. “quenching”) in massive ($\gtrsim L_{\ast}$) galaxies. In Su et al. (2018b), we show this explicitly, in non-cosmological simulations of halos with $M_{\rm halo}\sim 10^{12-14}\,M_{\sun}$ at $z\sim 0$ – demonstrating that all the physics studied here cannot solve or even substantially mitigate the “cooling flow” and quenching problems. Here in Figs. 18-19 we are essentially showing the same for the high-redshift progenitors of these massive halos (at $z\sim 2-10$, when their halo masses were in the range $M_{\rm halo}\sim 10^{12-13}\,M_{\sun}$). Not only do these massive halos sustain high SFRs ($\dot{M}_{\ast}\gtrsim 10-100\,M_{\sun}\,{\rm yr^{-1}}$) as long as we run them (including to $z\sim 0$) – i.e. clearly fail to quench – but they also form extremely dense central bulges in their starbursts with (in the most extreme cases) central circular velocities approaching $\sim 1000\,{\rm km\,s^{-1}}$ (well in excess of the most massive galaxies observed). But these are precisely the systems ($M_{\rm halo}\gtrsim 10^{12-13}\,M_{\sun}$ at $z\sim 2$, and $\gtrsim 10^{14}\,M_{\sun}$ at $z\sim 0$, with bulge-dominated $M_{\ast}\gtrsim 2\times 10^{11}\,M_{\sun}$) expected to host extremely massive super-massive black holes and quasars. And indeed, Anglés-Alcázar et al. (2017b) presented preliminary studies of the most extreme two halos here (m12z4 and m12z3) including models for black hole growth (but not AGN feedback), where the black holes reached masses $\gtrsim 10^{8-9}\,M_{\sun}$. So it will clearly be of particular interest in the future to explore the effects of AGN feedback in these systems.

5.5 Where Do CRs Act Most Efficiently? (Effects Inside versus Outside Galaxies)

We noted in § 5.3 above that CRs appear to have very weak effects on instantaneous properties within galaxies, even in the regime (e.g. large-$\kappa$ and high-$M_{\ast}$) where they have a large cosmologically-integrated effect on galaxy masses and SFHs.

5.5.1 Weak Effects Within Galaxies

For example, in Figs. 1-3 & 7, in the systems where CRs suppress SF, the metallicity and central circular velocity are also lower, accordingly, but these are essentially consistent with moving along, not off of, the observed mass-metallicity relation and Tully-Fisher relations, respectively (for more detailed studies of those, see Ma et al., 2016; El-Badry et al., 2018a). From these plots, we see that the galaxies with higher (lower) star formation rates have systematically higher (lower) baryonic masses in their centers, i.e. they appear to be shifting with the “supply” of gas. Fig. 8 shows this more rigorously, demonstrating that the different physics runs do not systematically differ in the normalization of the Schmidt-type relations – i.e. they are not consuming gas faster/slower or more/less “efficiently.” Rather, the galaxies with suppressed star formation have moved along the relation. This is quite different from what occurs if we increase/decrease the number or mechanical energy of supernovae, which systematically moves the relation down/up (for the same gas mass, fewer/more massive stars are required to regulate against gravitational collapse; see Hopkins et al. 2011; Ostriker & Shetty 2011; Faucher-Giguère et al. 2013; Agertz & Kravtsov 2015; Orr et al. 2018). Furthermore, although we defer a detailed study of the effect of CRs on galactic winds to future work, we find that gas outflow rates and velocities in the immediate vicinity of the galaxy are not strongly influenced by CRs – again unlike the case if we were to change the energetics or rate of SNe (see Hopkins et al., 2012, 2013a, 2013b, 2018c). It is possible, of course, that CRs contribute to outflows via “slow” or “gentle” acceleration of material at sub-virial velocities, and this increases at larger and larger radii (discussed below) but they do not directly launch “fast” outflows ($V_{\rm launch}\gtrsim V_{\rm escape}$).

Similarly, Figs. 10-13 show that the galaxy-scale stellar/visual and gas morphologies are only weakly modified even in any of our core suite. Where the stellar masses are strongly suppressed at high-$\kappa$ and high-$M_{\ast}$, the galaxies do tend to have slightly later-type visual morphologies, but this is entirely consistent with their lower masses – they simply resemble the earlier-time versions of their “Hydro+” counterparts (i.e. they have simply evolved less along the mass-morphology sequence; for more details of that see Garrison-Kimmel et al. 2017b). Conversely, some of the low-$\kappa$ runs which produce slightly higher stellar masses exhibit earlier-type morphologies. In the gas within the disks there are some slight differences in the “sharpness” of features in the cool gas in Figs. 10-13, but nothing qualitatively distinct. Fig. 9 shows some non-trivial (but still quantitative, rather than qualitative) effects on the ISM distribution of phases: namely, CRs can support more warm neutral gas within the galaxy. This is not surprising (they both heat and pressure-stabilize this phase of gas, without ionizing it significantly), but it is quantitatively non-negligible for detailed comparisons of ISM phase structure. Also, in Fig. 14 we see the CRs do not radically alter the kinematics of gas or stars (again, except insofar as they suppress the amount of low-redshift star formation; see El-Badry et al. 2018b for a detailed study of how this varies as a function of mass in our “Hydro+” runs).

All of this is expected from § 3 and Figs. 5-6. Within the star-forming galactic disk, CRs can have roughly equipartition-level energy densities, but their large diffusivity means that (as observed) the CR scale-height/length is much larger than the cold star-forming gas scale height (let alone the size of structures like GMCs, cores, etc). This means that CR pressure gradients – which actually determine the forces – are usually order-of-magnitude smaller than the small-scale forces from gravity, magnetic fields, gas thermal pressure, and turbulent ram pressure, in the multi-phase ISM.

5.5.2 (Potentially) Strong Effects in the CGM

On the other hand, we see strong CR effects in the CGM, even in many simulations (e.g. with lower $\kappa$) where the CRs do not have a large effect on galaxy masses.

This is especially evident in Figs. 10-13: for dwarfs with halo masses $M_{\rm halo}\gtrsim 10^{11}\,M_{\sun}$, the “CR+” runs often feature much more prominent cool gas in the CGM. Figs. 1-3 show that the total baryonic mass density on CGM scales ($\sim 10-100\,$kpc) is not drastically modified by CRs, although there are subtle (order-unity) changes evident in many cases (e.g. some of the intermediate-mass m11 halos with lower-$\kappa$ feature enhanced gas density on large scales). The runs where CRs dominate over thermal pressure in the CGM (see Fig. 5) typically feature modest (factor $\sim 2-3$) overall enhancements of gas density at some intermediate range of radii of order the halo scale radius ($R_{s}\sim R_{\rm vir}/10$), with little effect on total gas density out to or beyond the virial radius $\sim R_{\rm vir}$.

As noted above, Fig. 6 shows this enhancement is not driven by CR “heating” via either collisional or streaming processes. Most of the CGM mass in these runs is in “cool” or “warm” CGM phases, where the cooling times are relatively short and the temperature is maintained largely by photo-ionization equilibrium (Ji et al., 2019), hence their low (sub-virial) thermal pressure. Rather, Fig. 5 shows that the “maintenance” of this gas density profile owes to the CRs establishing a quasi-virial-equilibrium pressure profile on large scales. This is qualitatively similar to the findings of e.g. Salem et al. (2014, 2016); Chen et al. (2016) in their cosmological simulations with CRs, despite their adopting different numerical methods and treatment of the CR and star formation/feedback “microphysics,” suggesting the conclusion is robust to these details.

These, plus the weak effects of CRs on essentially all star formation and outflow and internal galaxy properties discussed above, demonstrate that the CRs primarily operate as a “preventive” feedback mechanism, rather than an “ejective” or “responsive” feedback mechanism. Rather than launching strong outflows, or removing gas from the halo, or preventing gas which has already accreted into the galaxy from efficiently forming stars, the primary effect of CRs appears to be preventing gas in the halo from actually accreting rapidly onto the galaxy. These CR effects result in a more subtle re-arrangement of gas mass within the halo and between different phases.212121Note that it is not possible here to completely dis-entangle the role of CRs preventing “new” gas from accreting onto the galaxy at all, versus “lofting” or “gently/slowly accelerating” cool gas in galactic fountains (at velocities $\ll V_{c}$) and preventing it from re-accreting. In an instantaneous sense these are the same thing: suppressing cool gas from falling into the galaxy. Future work following trajectories of individual gas elements will allow us to better disentangle these possibilities. Even when $\gtrsim L_{\ast}$ halos become dominated by a hot, virialized halo gas, cooling of that hot atmosphere or re-accretion of previously ejected gas can provide a substantial gas supply for late time star formation (Kereš et al., 2005; Muratov et al., 2015; Anglés-Alcázar et al., 2017a). Preventing such accretion (or re-accretion) has important consequences for the late time star formation and growth of galaxies.

Note that even a modest (factor $\sim 2$) effect on the gas density profile around the halo scale radius is significant – this is the radius where most of the halo gas mass resides, and most of the total baryon supply is in the halo (especially in the lower-mass galaxies) rather than in stars, so this relatively modest effect can easily account for differences of factors $\sim 2-5$ in star formation rates and stellar masses within the galaxies.

5.6 Cosmic Ray Transport

We now discuss the different detailed CR transport processes modeled here, their effects, and sensitivity to uncertain physics.

5.6.1 The Critical Role of the CR “Effective Diffusivity” (Transport Speed)

In § 5.2 & 5.3, we noted the existence of a “sweet spot” in the CR diffusion coefficient $\kappa$, demonstrated throughout but especially in Fig. 15. As predicted analytically in § 3 and confirmed in the simulations in Fig. 4, at too-low $\kappa$, CRs take too long to escape dense gas, and lose their energy rapidly to collisional processes. Advection alone (e.g. CRs being “carried” out of galaxies in super-bubbles and cold outflows) is sub-dominant even at quite low $\kappa$. Even without collisional losses, it requires a relatively high $\kappa$ for CRs to build up in the CGM before trans-Alfvénic streaming (which is also “lossy”) takes over. On the other hand, at arbitrarily high-$\kappa$, CRs would escape “too efficiently” even from the CGM/extended halo and the steady-state CR pressure ($\propto 1/\kappa$) would become too small to support gas or do any interesting work.

For dwarfs, this “special” value of $\kappa\sim 3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$, but even there the effects of CRs are weak, and moreover the observations of $\gamma$-ray luminosities from dwarfs like the SMC, LMC, and M33 favor higher $\kappa\gtrsim 3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ (where they escape low-density dwarfs more efficiently and do little work). But for massive (MW-mass) systems, interestingly, the “sweet spot” in $\kappa\sim 3-30\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ appears to neatly coincide with the observationally-favored values.

It is also encouraging that the more detailed study of CR transport physics in isolated (non-cosmological) simulations in Chan et al. (2018) identified approximately the same critical $\kappa$ needed to reproduce the $\gamma$-ray observations. In detail we favor slightly higher $\kappa$ here, owing to the larger gaseous halos present in cosmological simulations, which contribute to increasing the probability of CRs re-entering the disk (hence “residence time”) before they escape, but this is expected. We also note that the value we quote here is the parallel diffusivity. If one assumed isotropic CR diffusion, the corresponding isotropically-averaged $\tilde{\kappa}$ would be a factor $\sim 3$ lower.

At both mass scales, we find this “ideal” range of $\kappa$ spans approximately an order-of-magnitude in range, and within this range, the predictions are not especially sensitive to $\kappa$ (i.e. changing $\kappa$ by a factor $\sim 2$ about this “ideal” value produces relatively small effects). So although the physical scalings of CR diffusivity remain poorly-understood, this does not necessarily require particular fine-tuning.

This highlights, however, what is likely the most uncertain-but-important assumption in this work: that CR transport can be approximated with a constant diffusivity (or a fluid-like model at all). CR transport – what gives rise to “diffusive-like” behavior at all – remains deeply theoretically and observationally uncertain, as it is generally believed that diffusive behavior arises non-linearly from the interplay between different plasma instabilities and transport processes along field lines. For this reason we simplified and explored different, constant $\kappa$ models here. If, in nature, the “effective” $\kappa$ simply varies systematically from galaxy-to-galaxy (e.g. with cosmological environment or halo mass) this is not so problematic: it simply amounts to choosing a different “effective $\kappa$” as “most appropriate” for our different simulations. If the effective $\kappa$ varies on small spatial or time scales, this is also not such a concern, as the CRs diffuse sufficiently rapidly that these variations are averaged out on the large CGM scales ($\sim 10-100\,$kpc) on which they act (and indeed, because the transport is anisotropic, local turbulent magnetic field-line variations already mean that the actual local diffusivity is constantly varying by factors of several on small scales).

The biggest cause for concern is if, in reality, there are large systematic variations in effectively diffusivity with some property that varies between galaxies and CGM (e.g. spatial scale, density, magnetic field strength). It is plausible to imagine models where $\kappa\sim 3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ in the ISM (producing the same $\gamma$-ray luminosity), but CRs “de-couple” or “free-stream” or escape much more rapidly in the CGM, dramatically reducing their effects on galaxy formation (see § 6).

5.6.2 (Trans-Alfvénic) CR Streaming

We confirm the conclusion in Chan et al. (2018), that streaming at trans-Alfvénic speeds (as parameterized here) appears to be sub-dominant in transport. For observationally-favored (high) $\kappa$ values, transport is sufficiently fast that streaming at speeds $\sim v_{A}$ can only take over as a dominant transport mechanism, and begin to dissipate a large fraction of the CR energy, outside a radius $\gtrsim 30-100\,$kpc (§ 3 and Fig. 6). Thermalized energy from streaming losses are never, in the simulations here, able to compete with gas cooling losses (Fig. 6). As a result, turning off “streaming” entirely (Fig. 17) produces only minor effects on galaxy properties, around the favored $\kappa$.

Above, we argued that for CRs to have large effects, high “effective” $\kappa_{\ast}$ is needed in order for CRs to escape the dense star-forming regions of galaxies, and these high-$\kappa$ values are also favored by observations of dwarfs and MW-like galaxies (which appear to be well below calorimetric). Because the actual physics which regulates streaming speeds remains uncertain (and streaming losses will not appear in the $\gamma$-ray constraints), one might reasonably wonder whether the same effects might be accomplished by maintaining a relatively low diffusivity (e.g. $\kappa_{\|}=3\times 10^{28}\,{\rm cm^{2}\,s^{-1}}$) but increasing the streaming velocity by a modest factor. This is the experiment performed in Figs. 4 and 17. If the “streaming losses” scale with $v_{\rm stream}$ – as assumed in this test for the sake of illustration – then increasing $v_{\rm stream}$ is not equivalent to increasing $\kappa$. In fact, increasing $v_{\rm stream}$ and decreasing $\kappa$ in this manner leads to weaker, not stronger, effects of CRs on galaxy properties. And while increasing $v_{\rm stream}$ (with $\kappa$ constant) does slightly reduce the predicted $L_{\gamma}/L_{\rm SF}$ (Fig. 4), the effect is minor for the values considered here. There are two fundamental reasons for this. First, if the “streaming loss” term scales with $v_{\rm stream}$, then more of the CR energy is dissipated close to the galaxy (where cooling is efficient). Second, if streaming dominates then the effective diffusivity is $\kappa_{\ast}\sim v_{\rm stream}\,\ell_{\rm cr}$ (see § 2.3; where $\ell_{\rm cr}\equiv P_{\rm cr}/|\nabla_{\|}P_{\rm cr}|$), so for the large $\kappa_{\ast}$ favored we require $v_{\rm stream}\sim 1000\,{\rm km\,s^{-1}}(\kappa_{\ast}/3\times 10^{29}\,{\rm cm^{2}\,s^{-1}})\,(\ell_{\rm cr}/{\rm kpc})^{-1}$), much larger than $v_{A}$ or $c_{s}$.

So, if we set $v_{\rm stream}$ to be much larger ($\gtrsim 500-1000\,{\rm km\,s^{-1}}$), and limit streaming losses to scale with the Alfvén speed ($\sim v_{A}\,\nabla P_{\rm cr}$, as these come from the excited Alfvén waves which cannot propagate faster than this), then streaming will become essentially degenerate with high-$\kappa$ diffusion. We see this indirectly via the fact that the maximum “transport speed” of CRs has little effect on the results once it is sufficiently large (Fig. 22). How, micro-physically, these transport speeds arise remains uncertain.

5.6.3 Collisional and Streaming Losses

As shown in Fig. 6, thermalized CR collisional or streaming losses are essentially never important as a source of heating the gas (they are always highly sub-dominant to cooling). However, they can be important as loss mechanism for the CRs themselves. Although it is sometimes (incorrectly) assumed that CRs are “lossless” (or “don’t cool”), the collisional loss timescale for CRs in dense gas on length scale $\sim\ell$ is shorter than the diffusion time if $n\gtrsim 10\,{\rm cm^{-3}}\,\tilde{\kappa}_{29}\,(\ell/{\rm kpc})^{-2}$. And indeed, essentially all observed starburst galaxies, which have very high nuclear gas densities (see Fig. 4) are consistent with being proton-calorimeters (i.e. all or most of the CR energy appears to be lost). For that reason, at either low-$\kappa$ ($\ll 10^{29}\,{\rm cm^{2}\,s^{-1}}$) or high gas surface densities, collisional losses play an important role in limiting CR energetics and effects on galaxies. Otherwise, in “steady state” very low-$\kappa$ would produce more trapping and allow one to artificially build up essentially arbitrarily high CR energy densities. On the other hand, once $\kappa$ is relatively high, in galaxies with central densities comparable to or less than the MW, then most of the CRs escape without decaying (e.g. $L_{\gamma}$ is well below calorimetric in Fig. 4). In these cases, which include most of the cases of greatest interest above, collisional losses necessarily play a minor role, and so de-activating these losses in Fig. 17 produces relatively small effects under these conditions.

Streaming losses are similarly not dominant if $\kappa_{\ast}$ is sufficiently large, at least out to radii $r_{\rm stream}\sim\kappa/v_{\rm stream}$ where trans-Alfvénic streaming can begin to dominate the transport. Removing streaming losses with similar $v_{\rm stream}$ and $\kappa$ simply removes the losses at large $r$, making the equilibrium CR pressure profile fall somewhat less rapidly. But because, in this regime, it is already falling more rapidly than in the halo center (see Fig. 5), this is a second-order effect. Of course, as discussed above, if we substantially increase the streaming velocity beyond the Alfvén speed, and correspondingly increase the streaming loss rate (and shrink $r_{\rm stream}$ so the losses happen closer to the galaxy), then it can become important limiting the effect of CRs. Of course, as noted above, if we increase $v_{\rm stream}$ but remove or suppress the loss term so it does not scale accordingly, this is just equivalent to increasing $\kappa$.

5.6.4 Isotropic vs. Anisotropic Transport & Magnetic Fields

Fig. 17 considers some experiments where we remove magnetic fields, and assume all transport processes are isotropic (instead of projecting them along field lines). There, and in the other galaxy and CGM properties studied here, we see no radical or qualitative change in behavior. In the dwarfs, the results are very similar; in the MW-mass experiment (m12i), switching to isotropic transport and removing magnetic fields (Fig. 17) leads to a slightly weaker effect from CRs here, largely because the CRs escape more efficiently (but the difference is small compared to removing CRs entirely). In general, as shown in Chan et al. (2018), the leading-order effect of anisotropic transport is to alter the effective (galaxy-averaged) $\kappa$, by a modest (order-unity) factor, with a net effect of somewhat lower diffusion coefficients required in isotropic diffusion case for similar effect.

Some of the reason for this is that magnetic fields are highly disordered on most scales, especially in the CGM. This is shown more rigorously in previous studies (e.g. Su et al., 2018a), and will be explored in subsequent work as well (Ji et al., 2019). This is not surprising, as the CGM of $\sim L_{\ast}$ and dwarf galaxies is generically trans-sonically turbulent and has plasma $\beta\gg 1$ (so field lines are essentially passively advected into locally disordered configurations). So typical anisotropic “suppression factors” vary within order-unity values but strong, systematic suppression of CR transport is never really possible. This may differ in more massive halos (e.g. clusters), where the turbulent Mach numbers in the (much hotter) diffuse gas are expected to be much smaller.

Related to this, we see no evidence for CRs systematically “opening up” field lines on large scales. This is also not surprising, as we have already argued the primary role of CRs in the simulations here is not in violently ejecting material from small radii to large, but in quasi-statically maintaining the halo.

5.7 Resolution & Numerical Effects

We have noted in several places above various numerical tests in our simulations. For example, we have run our standard physics set (“Hydro+”, “MHD+”, “CR+($\kappa=3e28$)”, “CR+($\kappa=3e29$)”) at several different resolution levels in boxes m10q (with resolution $m_{i,\,1000}=0.25,\,2.1,\,16$), m11q ($m_{i,\,1000}=0.88,\,7,\,56$), m12i, m12f, and m12m (each with $m_{i,\,1000}=7,\,56,\,450$), i.e. $60$ simulations. Some examples are shown in Appendix A (and in Fig. 7). As shown exhaustively in Hopkins et al. (2018b), the dwarf galaxies m10q and m11q show very weak resolution-dependence in “Hydro+” runs (or “MHD+”, in Su et al. 2017). We find the same for “CR+” runs, reinforcing our conclusions. As studied extensively in Hopkins et al. (2018b), our massive halos m12i, m12f, m12m do show systematic resolution dependence in all the runs. This is of course important and the scalings and reasons for this are discussed in Hopkins et al. (2018b). What is important, for our purposes here, is that our qualitative conclusions about the systematic effects (or lack thereof) of MHD and CRs, are independent of resolution. Like the dwarf runs, at every resolution we find the MHD runs and low-$\kappa$ CR runs have little or no systematic effect on the massive halos, while the high-$\kappa$ CR runs suppress SF significantly, via the same physical mechanisms. The effect is slightly stronger in the higher-resolution runs because the galaxy is overall lower-density and lower-mass even in the “Hydro+” runs (this owes to better resolution of galactic winds “venting” and mixing in the CGM), which means (per Fig. 4) the collisional losses are somewhat less.

Chan et al. (2018) also considered extensive tests of the imposed maximum CR free-streaming speed $\tilde{c}$, and showed that as long as this is faster than typical bulk velocities of gas in the simulated systems, it has no effect on the results: Appendix A shows this holds in our simulations so long as $\tilde{c}\gtrsim 500\,{\rm km\,s^{-1}}$. More detailed resolution tests and idealized validation tests of our numerical CR methods, and explicit comparison of the results of two-moment vs. one-moment approximations for the CR transport flux solver, are all presented in Chan et al. (2018): none of these presents obvious numerical concerns here (see also Appendix A).

We stress that extensive numerical tests of almost every other aspect of these simulations (resolution, force softening, hydrodynamic solvers, radiation-hydrodynamics solvers, etc.) are presented in Hopkins et al. (2018b), Hopkins et al. (2018c), and Hopkins et al. (2018a), for our “Hydro+” models. For more detailed resolution and physics tests of the “MHD+” models, we refer to Su et al. (2017, 2018a), and for extensive numerical validation tests of the MHD, conduction, viscosity, and anisotropic transport solvers we refer to Hopkins & Raives (2016); Hopkins (2016, 2017).

6 Conclusions

6.1 Overview

We present and study a suite of $>150$ new high-resolution ($\sim 100-10000\,M_{\sun}$, $\sim 1-10\,$pc, $10-100\,$yr, $\sim 10^{3}-10^{4}\,{\rm cm^{-3}}$) FIRE-2 cosmological zoom-in simulations, with explicit treatment of stellar feedback (SNe Types Ia & II, O/B & AGB mass-loss, photo-heating and radiation pressure), magnetic fields, fully-anisotropic conduction and viscosity (accounting for saturation and their limitation by plasma instabilities at high-$\beta$), and cosmic rays. Our CR treatment accounts for injection in SNe shocks; advection, fully-anisotropic streaming and diffusion; and losses from collisional (hadronic+Coulomb), streaming, and adiabatic processes. We systematically survey different aspects of the CR physics and uncertain transport coefficients. We examine the effects of these physics on a range of galaxy properties including: stellar masses, star formation rates and histories, metallicities and abundances, stellar sizes and baryonic mass profiles, dark matter mass profiles, rotation curves, visual morphologies and kinematics of stars and gas, and gas phase distributions. We survey these properties across a suite of simulations spanning all redshifts, and masses at $z\sim 0-10$ ranging from ultra-faint dwarf ($M_{\ast}\sim 10^{4}\,M_{\sun}$, $M_{\rm halo}\sim 10^{9}\,M_{\sun}$) through Milky Way mass. We summarize our conclusions as follows:

• •

We confirm the growing body of work showing that magnetic fields, physical conduction, and viscosity on resolved scales have little effect on any galaxy property studied (Figs. 1-3 and 7-8). These simulations reach higher resolution (sufficient to resolve the Field length in warm and hot gas with $T\gtrsim 2\times 10^{5}\,{\rm K}\,(n/0.01\,{\rm cm^{-3}})^{0.4}$), and include more detailed treatment of the effect of plasma instabilities on transport coefficients, compared to our previous work (Su et al., 2017), but this only serves to reinforce that conclusion. It is of course possible there are important un-resolved effects which could be important via “sub-grid” effects (e.g. altering the effective cooling rates or stellar initial mass function).

• •

Magnetic fields are highly-tangled, at all mass scales we survey, with or without CR physics (see also Ji et al. 2019). Per Fig. 5, the plasma $\beta\equiv P_{\rm thermal}/P_{\rm magnetic}$ varies enormously in the ISM ($\beta\sim 0.1-100$) but in warm/hot phases is usually large (median $\sim 1-10$ in the diffuse ISM in MW-mass galaxies, and larger $\sim 10-30$ or $\sim 30-300$ in lower-mass $M_{\rm halo}\sim 10^{11}\,M_{\sun}$ and $\sim 10^{10}\,M_{\sun}$ dwarfs, respectively), consistent with many recent studies of amplification (see e.g. Martin-Alvarez et al., 2018, and references therein). It rises in the CGM with galacto-centric distance (to a median $\sim 100-1000$ at $\gtrsim R_{\rm vir}$, with local variations reaching $\sim 10^{4}$).

• •

CRs have relatively weak effects on the galaxy properties studied, in dwarfs with $M_{\rm halo}\ll 10^{11}\,M_{\sun}$ ($M_{\ast}\ll 10^{10}\,M_{\sun}$), for essentially any physical CR parameters considered, once realistic mechanical and radiative feedback are already included (Figs. 1, 7, 10). This is both predicted and easily understood from basic energetic considerations (§ 3). Previous claims to the contrary have either failed to account for (a) the fact that realistic dwarfs have very low $M_{\ast}/M_{\rm halo}$ (and correspondingly low SFRs, hence SNe and CR energy injection rates), without much “hot halo” gas, (b) realistic supernova and radiative feedback which easily overwhelm the effects of CRs in dwarfs, and/or (c) observational constraints from $\gamma$-ray emission in dwarfs which place upper limits on the collisional loss rate and require that $>90-99\%$ of the CR energy escape without hadronic losses (prohibiting low transport speeds). It remains possible CRs modify the CGM in these dwarfs, or modify processes like ram-pressure stripping in dwarf satellites around massive galaxies (not studied in detail here).

• •

CRs (from SNe) also have relatively weak effects on galaxy properties, at any mass scale, at high redshifts $z\gtrsim 1-2$ (Fig. 18). Although SFRs (and therefore CR production rates) are higher (at a given mass) at high-$z$, the CGM density and ram pressure of dense, cold inflows (carrying most of the mass) is much higher, and so CR pressure is insufficient to strongly suppress those inflows (Fig. 19). Moreover the periods of strongest inflow are often accompanied by dense starbursts within the galaxy where surface densities exceed $\gg 0.1\,{\rm g\,cm^{-3}}$ ($\gg 1000\,M_{\sun}\,{\rm pc^{-2}}$), during which CRs experience strong hadronic losses. Consistent with observed low-redshift starbursts at these densities, they become approximate proton calorimeters.

• •

CRs can have significant effects in massive galaxies ($M_{\rm halo}\gtrsim 10^{11}\,M_{\sun}$, $M_{\ast}\gtrsim 10^{10}\,M_{\sun}$) at relatively low redshifts ($z\lesssim 1-2$), reducing their peak star formation rates by as much as factors $\sim 5$ and $z=0$ stellar masses by factors $\sim 2-3$ (Figs. 2, 3, 7). This in turn significantly reduces the central peak in their rotation curves, and moves the galaxies along the mass-metallicity, Kennicutt-Schmidt, morphology-mass, and other mass-based scaling relations. This maximal effect requires effective diffusivities in the range $\kappa_{\ast}\sim 3\times 10^{29}-3\times 10^{30}\,{\rm cm^{2}\,s^{-1}}$ (which can arise from a combination of microphysical diffusion and/or streaming).

• •

In these systems, the CRs primarily operate on the CGM, and have relatively little direct effect within the galaxies (e.g. Figs. 8, 9, 14). Essentially any model where CRs are sufficiently well-trapped in galaxies to, e.g. slow down SF directly in relative dense ($>1\,{\rm cm^{-3}}$) disk gas, or launch strong outflows (velocities $\gtrsim V_{c}$) is ruled out observationally and results in excessive CR losses (greatly limiting their net effect). We see essentially all galaxy properties move along, not off of, standard scaling relations. This is dramatically different from e.g. increasing the mechanical energy or rate of SNe. Instead, CR feedback is primarily “preventive.” In massive halos at low redshifts, with the appropriate $\kappa_{\ast}$, the CRs establish a quasi-virial-equilibrium pressure profile, which dominates over the gas thermal and magnetic pressure from $\sim 20-200\,$kpc, and supports warm+cool gas ($T\ll 10^{6}\,$K) which would otherwise cool and rain onto the galaxy (Fig. 5 and Ji et al. 2019).

• •

Given present uncertainties, the most important-yet-uncertain parameter determining the effects of CRs (especially in massive galaxies) is the effective diffusivity $\kappa_{\ast}$. If too low, CRs are trapped and (a) lose energy to collisional processes in dense gas (negating their ability to do work, and directly violating observational constraints from spallation and $\gamma$-ray emission), and (b) cannot propagate to large-enough radii to slow accretion in the outer halo. If $\kappa_{\ast}$ is too high, CRs should simply escape and their equilibrium pressure, even in the outer halo, would be too low to do interesting work. The “sweet spot” is approximately an order-of-magnitude in width, and agrees well in the simulations and simple analytic models.

• •

We show that with $\kappa_{\ast}\sim 3-30\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$, cosmological simulations reproduce observed $\gamma$-ray emissivities (and similar constraints from spallation in the MW), at dwarf through Local Group through starburst-galaxy density/star formation rate scales (Fig. 4). This echoes the conclusion from recent non-cosmological simulations (Chan et al., 2018), for the first time in cosmological simulations.

• •

Cosmic ray “heating” of the gas (via energy transfer through collisional or streaming losses) is orders-of-magnitude smaller than gas cooling rates and never important for the gas (Fig. 6). It is however important as a loss mechanism for CRs, especially in dense starburst-type systems (which become proton calorimeters) or the outer halo (where trans-Alfvénic streaming can dominate).

• •

Around the effective diffusivities of interest, anisotropy and streaming at the Alfvén speed play second-order roles in CR transport (Fig. 17). Anisotropy does not have radical obvious effects, but does lead to somewhat suppressed diffusion, especially in dense gas where losses can be large (necessitating larger-$\kappa_{\ast}$ to enable escape). Trans-Alfvénic streaming is too slow to account for the required $\kappa_{\ast}$ in the ISM and inner CGM, although it can dominate transport outside of some radius $\sim\kappa_{\ast}/v_{A}$ (in the outer halo).

• •

We present a simple analytic equilibrium model in § 3, which is able to at least qualitatively predict essentially all of the relevant CR effects here. Particularly in the mass and $\kappa_{\ast}$ range of greatest interest, this analytic model provides a remarkably accurate description of when CRs become important, the favored values of $\kappa_{\ast}$, and the resulting CR pressure profile and equilibrium gas density profile in the CGM of CR-dominated halos.

6.2 Future Work & More Massive Galaxies

This study – despite its length – is far from comprehensive. In future work, we plan to explore a number of properties of these simulations in greater detail, including: their more detailed outflow and gas phase structure, properties of GMCs in the ISM and accretion onto the galactic disk, observable ions in the CGM, properties of resolved satellites around massive galaxies in the modified CR-dominated CGM. There may be more detailed tracers which can distinguish between different models in greater detail, or areas where MHD/conduction/viscosity/CRs have large effects which we have failed to identify here. Moreover, even assuming our simple CR treatment is a reasonable representation, there are a number of basic physical processes meriting further investigation. For example it is unclear how the non-linear thermal instability operates in a stratified halo supported primarily by a CR pressure gradient (nor is it likely this can be fully resolved in the simulations here). It is not obvious how much of the additional cool gas in the halo in the high-$\kappa_{\ast}$, high-$M_{\ast}$ CR runs owes to “pure suppression” of new inflowing material vs. CRs “gently” (and slowly) re-accelerating or “lofting” cool gas which would otherwise recycle in galactic fountains. And it is worth exploring how the virial shock and transition “out of” the CR dominated-regime at large radii occur. These and many additional questions clearly merit exploration in idealized, high-resolution numerical experiments.

Future and parallel work will also explore the role of CRs in more massive ($>L_{\ast}$) galaxies. We predict and find in our simulations that the strength of the effects from CRs sourced via SNe scales as $\propto\dot{M}_{\ast}/M_{\rm halo}$, which is maximized around $\sim L_{\ast}$. In very massive halos (in nature), this declines both because (a) $M_{\ast}/M_{\rm halo}$ drops, and (b) star formation is “quenched.” As a result, our parallel study in Su et al. (2018b) demonstrates that CRs sourced by SNe (including Ia’s) cannot possibly solve the “cooling flow problem” and resist excessive cooling and star formation in very massive ($M_{\rm halo}\sim 10^{14}\,M_{\sun}$) halos. And here we show essentially the same at higher redshifts when these halos initially form most of their stars and dense, compact bulges (when their progenitor halo masses are $\lesssim 10^{12-13}\,M_{\sun}$).

However, in those same massive galaxies, super-massive black holes (SMBHs) and AGN appear to dominate CR production (seen in e.g. jets and “bubbles”) by orders-of-magnitude compared to SNe, and in this paper we only accounted for CRs produced in SNe. In a companion study (Su et al., 2018c) we show that injection of $\sim 10^{43}\,{\rm erg\,s^{-1}}$ in CRs in a $\sim 10^{14}\,M_{\sun}$ halo – modest for the SMBHs and AGN in those systems, but $\sim 100-1000$ times the production rate from SNe assuming $\sim 10\%$ of the SNe energy goes into CRs – can have a dramatic effect on the halo gas and cooling flow in these massive galaxies. Obviously it will be important to revisit these high-redshift halos with an explicit model for AGN feedback.

6.3 Fundamental Physical Uncertainties

We wish to strongly emphasize that our conclusion is only that CRs could be important to massive galaxies at low redshifts, not that they are necessarily. We have parameterized our ignorance by adopting a simple two-moment model with fixed parallel diffusivity $\kappa_{\|}$, but the reality is that deep physical uncertainties remain. It is not clear what actual physics determines the “effective diffusivity” or transport parameters of CRs in all the regimes here, let alone how these parameters should scale as a function of local plasma properties (e.g. magnetic field strength, density, mean-free-path, strength of turbulence, etc.). Although constraints from detailed modeling of the Milky Way CR population and $\gamma$-ray observations of nearby galaxies favor an effective diffusivity $\kappa_{\ast}\gtrsim 10^{29}\,{\rm cm^{2}\,s^{-1}}$, this constraint is (a) only measured in a few $z=0$ galaxies at present, and (b) more importantly, only constrains the effective diffusivity or “residence time” in the relatively dense, star-forming galactic disk gas (where the hadronic collision rate, which scales as $n_{\rm cr}\,n_{\rm nucleon}$, is maximized). It is completely plausible to imagine models where $\kappa_{\ast}$ is $\sim 3\times 10^{29}\,{\rm cm^{2}\,s^{-1}}$ in this gas, but once CRs reach the CGM at $r\gg 10\,$kpc (where the gas has low-density, very weak magnetic fields $\beta\sim 100-10^{4}$, and long mean-free paths $\sim 1-100\,$pc), the transport speed rises dramatically and CRs simply “escape” rather than forming a high-pressure halo. This is, essentially, the classic “leaky box” model. Even if the CRs are confined in the halo on $\sim 100\,$kpc scales, it is not clear whether the local tight-coupling approximation is valid in the tenuous CGM – i.e. can CR “pressure” actually be simply “added” to to the local gas stress tensor in the MHD equations? Or can CRs “slip” or couple with non-negligible lag? Fundamentally, our treatment of CRs as a fluid depends on assuming something about their ability to reach a micro-scale gyrotropic equilibrium distribution function, which may not be valid when scattering rates are low, i.e. transport speeds are very large (see e.g. Holcomb & Spitkovsky 2019). These are critical questions which could completely alter our conclusions (in the examples above, they would make the effects of CRs much weaker), and may fundamentally require PIC-type simulations that can follow explicit kinetic plasma processes to answer.

Observationally, direct constraints on CRs in the regimes of interest are unfortunately very limited. Essentially the only such measure outside the MW is the $\gamma$-ray luminosity discussed extensively here. Recall, most of the CR energy/pressure is in $\sim$ GeV protons. So constraints on the CR electrons (e.g. synchrotron), which tend to dissipate their energy far more efficiently and closer to sources, are not particularly illuminating for this specific question, nor are constraints on the high-energy CR population (which has very different diffusivity and contains negligible CR pressure). However, given that the regime of greatest interest is precisely where we predict CRs should have a large effect on the temperature and density distribution of CGM gas around $\sim L_{\ast}$ galaxies – a topic of tremendous observational progress at present – it is likely that indirect observational constraints from the CGM will represent, in the near future, the best path towards constraining the CR physics of greatest interest here.

Acknowledgments

We thank the anonymous referee for a number of insightful comments and suggestions. Support for PFH and co-authors 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. DK was supported by NSF grant AST-1715101 and the Cottrell Scholar Award from the Research Corporation for Science Advancement. CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522, by NASA through grant 17-ATP17-0067, and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. Numerical calculations were run on the Caltech compute cluster “Wheeler,” allocations from XSEDE TG-AST130039 and TG-AST120025 and PRAC NSF.1713353 supported by the NSF, and NASA HEC SMD-16-7592.

Appendix A Additional Numerical Tests

In this Appendix we present some examples of additional numerical tests, discussed in § 5.7 in the main text.

Figs. 20-21 show resolution tests, for MW-mass halos and varied diffusivity $\kappa$. As noted in § 5.7 we have also extensively varied the resolution in our dwarf halos m10q and m11q, but since (a) those halos exhibit no large effect from either MHD or CRs (at any $\kappa$) and (b) our extensive resolution studies in Hopkins et al. (2018b) (as well as Su et al. 2017 and Chan et al. 2018) have shown the dwarf properties in “Hydro+” and “MHD+” runs are extremely insensitive to resolution, it is not surprising that we find there is no change at any resolution level for these halos (this is also demonstrated in Fig. 7). Therefore we focus on MW mass where both CR physics and resolution have much larger effects.

As discussed in § 5.7, Figs. 20-21 do show a systematic resolution dependence in our massive m12 halos (compare also Fig. 20 to 3). In massive halos, lower-resolution runs produce larger stellar masses and, correspondingly, higher central baryonic densities and enhanced circular velocities. This is all studied in extensive detail in Hopkins et al. (2018b). Critically, for our study here, we see the same systematic effects (or lack thereof) of MHD and CRs at all resolution levels.

In Fig. 22, we systematically compare halos m10q, m11q, m12i, chosen with the values of $\kappa$ where effects on the galaxy are maximal, varying the numerically-imposed maximum CR free-streaming speed $\tilde{c}$. As expected, we see that this has little effect on our conclusions (variations are smaller than stochastic run-to-run variations) for essentially all values $\tilde{c}\gg 200\,{\rm km\,s^{-1}}$, fast enough that CRs can (when they should, according to the diffusivity $\kappa$ and streaming speed $v_{\rm stream}$) diffuse or stream faster than local gas rotation/outflow velocities. Otherwise, CR effects are somewhat weaker, as the CRs can spend more time artificially trapped in dense gas (where collisions sap energy). This is consistent with our more detailed numerical study in Chan et al. (2018).

Chan et al. (2018) also describe and test in detail the CR flux equation (Eq. 5), i.e. the equation for the diffusive CR flux $\tilde{\bf F}_{\rm cr}$. They showed that using alternative forms of this equation, or simply omitting it entirely and solving directly a single 0th-moment (pure diffusion) CR energy equation (setting $\tilde{\bf F}_{\rm cr}\equiv\kappa_{\ast}\,\nabla_{\|}e_{\rm cr}$) give essentially identical results. We confirm this in our cosmological simulations in Fig. 23. We re-run a low and high-mass halo where CRs have a large effect (m11i and m12i), adopting the alternative formulation of the flux equation from Thomas & Pfrommer (2018):

[TABLE]

with $\hat{\tilde{\bf F}}_{\rm cr}\equiv\hat{\bf B}$ by definition. As discussed in Chan et al. (2018), this is essentially the same as our formulation, up to whether $\hat{\tilde{\bf F}}_{\rm cr}$ appears inside or outside of the derivative terms. The Thomas & Pfrommer (2018) formulation arises if we assume the perpendicular fluxes that arise when magnetic fields rotate are instantaneously converted into gyrotropic motion by micro-scale instabilities; this and the assumption of frame in which the motion is gyrotropic produce the small “pseudo-force” correction $\tilde{\bf F}_{\rm cr}\cdot\{(\hat{\tilde{\bf F}}_{\rm cr}\cdot\nabla)\,{\bf v}_{\rm gas}\}$. In practice, the choice of Eq. 15 instead of our default Eq. 5 only produces differences below the CR mean free path, and has no appreciable effect here. Unsurprisingly, the flux equation from Jiang & Oh (2018), which gives behavior “in-between” our Eq. 5 and Eq. 15 above (as shown in Chan et al. 2018) is even more similar to our default.

We have also experimented with a modified CR pressure term in the hydrodynamic equation motivated by Jiang & Oh (2018): replacing $\nabla P_{\rm cr}\rightarrow\nabla_{\bot}\,P_{\rm cr}-(\gamma_{\rm cr}-1)\,({\bf F}-{\bf v}_{\rm st}\,[e_{\rm cr}+P_{\rm cr}])/\kappa_{\ast}\approx\nabla P_{\rm cr}+\tilde{c}^{-2}\,[\partial\tilde{\bf F}/\partial t+\nabla\cdot({\bf v}_{\rm gas}\tilde{\bf F}_{\rm cr})]$, i.e. keeping the perpendicular CR pressure but only the parallel component contributing to CR scattering. Because this correction only enters in the $\sim 1/\tilde{c}^{2}$ terms, and because field lines are tangled, this has only weak effects.

Finally, Fig. 23 also considers two additional variants of the streaming speed: taking $v_{\rm stream}=v_{A}$ (as compared to our default $=3\,v_{A}$), or setting $v_{\rm stream}=0$ but retaining the “streaming loss” term: in all of these the loss term scales as $v_{A}\,\nabla P_{\rm cr}$. Given the much larger variations to the streaming speed considered in the main text (which produce small effects), it is expected that the effect of these variations is very small.