Cosmological Simulations of Dwarf Galaxies with Cosmic Ray Feedback

Cosmological Simulations of Dwarf Galaxies with Cosmic Ray Feedback

Jingjing Chen, Greg L. Bryan, Munier Salem
Department of Astronomy, Columbia University, New York, NY 10027, USA

We perform zoom-in cosmological simulations of a suite of dwarf galaxies, examining the impact of cosmic-rays generated by supernovae, including the effect of diffusion. We first look at the effect of varying the uncertain cosmic ray parameters by repeatedly simulating a single galaxy. Then we fix the comic ray model and simulate five dwarf systems with virial masses range from 8 to M. We find that including cosmic ray feedback (with diffusion) consistently leads to disk dominated systems with relatively flat rotation curves and constant star formation rates. In contrast, our purely thermal feedback case results in a hot stellar system and bursty star formation. The CR simulations very well match the observed baryonic Tully-Fisher relation, but have a lower gas fraction than in real systems. We also find that the dark matter cores of the CR feedback galaxies are cuspy, while the purely thermal feedback case results in a substantial core.

galaxies: formation - methods:numerical

1 Introduction

Isolated dwarf galaxies are commonly observed in the local Universe, and provide an interesting way to test models of galaxy formation and evolution. Below a stellar mass limit of about M (or rotation velocity of about 100 km/s), such systems are largely star-forming systems except when they become satellites of larger halos (Geha et al., 2012; Bradford, Geha & Blanton, 2015). Naively, one might expect such systems to be easy to understand, as they have small or non-existent bulges and are free from uncertainties related to AGN feedback. However, detailed comparisons with cosmological models have encountered two kinds of problems.

The first problem has to do with the baryonic content of the halos – simply put, galaxy formation must become increasingly less efficient as the halo mass decreases, for standard CDM to match observations. It has long been recognized that this is implied by the different slopes of the low-mass end of the dark matter halo mass function and the faint end of the luminosity function (White & Frenk, 1991). More recent work has sharpened this mismatch with both semi-analytic models (e.g., Guo et al., 2010; Behroozi, Wechsler & Conroy, 2013; Leauthaud et al., 2012) and simulations (e.g., Stinson et al., 2007; Governato et al., 2010; Sawala et al., 2010). A closely related problem is the ”missing” dwarf satellites around Milky-Way mass halos (Moore et al., 1999; Springel et al., 2008).

A second issue has to do with the inner structure of dwarf galaxies – in particular the relation between the halo mass and the rotational velocity near the core. Dark matter only simulations predict quite cuspy profiles with relatively high dark matter densities in the center, while observations generally find cored profiles (e.g., Moore, 1994; Salucci & Burkert, 2000; de Blok & Bosma, 2002; Walter et al., 2008). A related version of this problem has recently been highlighted among the dwarf satellites as the ”too big to fail” problem, which highlights the paucity of satellites with rotation velocities in the 30-70 km/s range (Boylan-Kolchin, Bullock & Kaplinghat, 2012).

Solutions to both of these problems in an astrophysical context require efficient feedback mechanisms that can blow out a large amount of gas, suppress star formation, bring down the baryonic fraction, and possibly reduce the central dark matter density of lower mass dwarfs. It has long been clear that the low potential depth of dwarf systems makes them particularly susceptible to stellar feedback (Dekel & Silk, 1986); however, modeling feedback from supernovae is numerically challenging (e.g., Navarro & White, 1994).

A variety of numerical feedback mechanisms have been tested in different papers. For example, by suppressing cooling immediately after a star formation event to counteract the artificially enhanced radiative losses (Governato et al., 2007; Guedes et al., 2011), simulated galaxies produced less massive bulges. Supernovae feedback which injects momentum directly in to the surrounding gas has also been shown to drive mass-loaded outflows (e.g., Oppenheimer & Davé, 2008). This has been combined with other physical processes – for example, Sawala et al. (2011) combined supernovae feedback with photoionization from ultraviolet background, and produced dwarf galaxies with high mass-to-light ratios.

More recent simulations have made progresses in reproducing the observed properties using subgrid models with strong feedback (e.g., Di Cintio et al., 2014; Hopkins et al., 2014; Schaye et al., 2015). For example, Sparre et al. (2015) showed that they could reproduce the star formation main sequence relation at . Zoom-in simulations afford more resolution, and have shown much recent promise in reproducing the stellar mass content of dwarf satellites (e.g., Governato et al., 2010, 2012; Hopkins et al., 2014; Oñorbe et al., 2015). Brooks & Zolotov (2014) found that energetic supernovae feedback and subsequent tidal stripping has a significant effect on reducing the dense core of Milky Way mass galaxies

Here we explore the impact of a specific physical feedback mechanism on dwarf galaxies, namely dynamical cosmic ray (CR) feedback. CRs are high-energy, quickly diffusing particles accelerated by shock waves. They have a wide range of energies, with the largest contribution to the energy density coming from those particles with energies around a few GeV. The energy density of CRs is comparable to the magnetic energy and turbulent energy in galaxies (Kulsrud, 2005).

In this paper, we present our cosmological simulations of dwarf galaxies with CR feedback. We use the Enzo code, with an implementation of CR feedback by Salem & Bryan (2014). This adopts a relatively simple two-fluid model that provides an extra pressure term and also allows for diffusion of the CRs (Drury, 1985). In Salem & Bryan (2014) and Salem, Bryan & Corlies (2016), we applied this model to the cosmological simulation of a Milky Way mass galaxy. Here we examine the effect of CR feedback and diffusion on dwarf galaxies with halo masses in the range of 8 to M.

In section 2, we describe the simulation model, the feedback scheme, initial conditions, and the selection of halos for further simulations. In section 3.1, we choose the best set of parameters by comparing the simulation result of a single system with varying CR physics. Then we apply this model to a series of dwarf galaxies, with halo mass spanning from 8 to M in section 3.2, and we evaluate the results by comparing to observed scaling relations. Finally, in Section 4, we discuss our results.

Figure 1: Edge-on views of the same halo modeled with different CR physics. From left to right, the panels show a simulation with no CRs, followed by runs with CR diffusion coefficients of 0, , , cms. From top to bottom, the panels show the projected (density-weighted) gas density, stellar density, and CR energy density. Each panel represents a kpc kpc area.

2 Method

We begin with a brief introduction to our hydrodynamics method, and then describe the two-fluid model used to model the CRs, and finally outline the cosmological initial conditions.

2.1 Numerical code: Enzo

Our simulations use the open-source code Enzo, which is a three-dimensional, Cartesian, grid-based hydrodynamics code that includes adaptive mesh refinement (Bryan & Norman, 1997; O’Shea et al., 2004; Hummels & Bryan, 2012; Bryan et al., 2014). We include a non-equilibrium chemical model following H, H, He, He, He, and include radiative cooling from all the above chemical species as well as metal cooling computed in a lookup table as described in Smith, Sigurdsson & Abel (2008). We include a metagalactic radiative background as computed in Haardt & Madau (2012). The star formation recipe of Enzo is based on the criteria presented in Cen & Ostriker (1992) and described in detail in Bryan et al. (2014). Briefly: when a parcel of dense gas is identified to be Jeans unstable and collapsing, some of its gas may be converted into a star particle with some efficiency. For a more detailed discussion of the star formation (and thermal-only) feedback scheme used in this paper, see Hummels & Bryan (2012) – in this work we adopt an efficiency of star formation per free fall time of 1% and an energetic feedback efficiency of of the stellar rest mass energy.

2.2 Cosmic Ray Feedback Scheme

We use a model for CR feedback and diffusion as implemented, tested and described in Salem & Bryan (2014). The model uses a two-fluid method (Drury, 1985; Drury & Falle, 1986; Jun, Clarke & Norman, 1994). It treats the high-energy particles as an ideal gas with . Anisotropic CR pressure and the dynamical impact of magnetic fields are assumed to be subdominant. The fluid equations for the thermal gas stay as usual, except that the momentum of the gas is affected by the pressure of both the gas and CR components. In addition, we optionally include isotropic diffusion of the cosmic rays. Here we simply highlight the additional equation that describes the evolution of the CR energy density :


where u is the fluid velocity, P is the pressure from CR, is the CR isotropic diffusion coefficient, and is the source term for CR. Except for the isotropic diffusion, the CRs are assumed to be tightly coupled to the gas through a tangled magnetic field (which is not explicitly simulated). While this model is approximate, simulations including magnetic fields and/or anisotropic diffusion have also been shown to drive winds with similar mass-loading factors (e.g., Hanasz et al., 2013; Girichidis et al., 2016). We also do not include any non-adiabatic losses in the CRs; in reality, CRs lose energy primarily through collisions with the thermal plasma; this is generally subdominant to adiabatic cooling in our simulations except for runs with zero (or low) diffusion coefficients, which are not realistic in any case. A more detailed discussion of this point and the shortcomings of the model in general can be found in Salem & Bryan (2014) and Salem, Bryan & Corlies (2016).

When a Type II SNe is formed, it will eject some of its mass and energy back into the surrounding cells. In our model, we assume that some of the energy is in thermal form and some accelerated by blast waves (generally not resolved in the grid). The fraction of energy that is deposited to thermal gas and CR is 0.7 and 0.3 respectively; this CR fraction is somewhat larger than canonical values but within suggested ranges (e.g., Wefel, 1987; Ellison et al., 2010; Kang & Jones, 2006). Moreover, we assume a relatively steep initial mass function (our chosen corresponds to 1 SN per 185 M masses of stars produced), and hence a relatively low SN energy injection rate, so the CR energy produced per solar mass of stars formed is within typical estimates. We adopt these values for consistency with previous work (Hummels & Bryan, 2012; Salem & Bryan, 2014), where we also explored the impact of varying them.

2.3 Initial Conditions and Halo Selection

We are interested in modeling a number of dwarf galaxies forming out of cosmological initial conditions. Therefore, we use MUSIC (Hahn & Abel, 2011) to set up the initial density and velocities fields. This is done using the cosmological parameters from the best fit Planck 2013 results (Planck Collaboration et al., 2014), i.e. H = 67.11 km/s/Mpc, , , . All the simulations start at and stop at the present day.

ID CR physics R M M
H1-0 no CR

R is in unit of kpc. M and M are in unit of M.

Table 1: Simulation parameters for fixed halo, varying CR runs.

First we perform a dark matter only run in a (20 Mpc/h) box and identify halos that are around to M. All of these halos become a little more massive in the final runs, listed in Table 2. To study the impact of CR feedback, we select relatively isolated halos so that the physical properties of the dwarfs are not affected by nearby massive galaxies. In particular, all the selected systems are at least 0.5 Mpc away from any halos that are more massive than 0.1 times of the target halo. Then we trace the dark matter particles currently in the halo to their positions at the beginning of the simulation. The selected halos are re-simulated at high-resolution with baryons and adaptive mesh refinement (AMR) is performed in the cubical region that encloses all the dark matter particles from the beginning. Our root grid is 128 and we use three additional levels in the initial conditions for a particle mass of M M. During the simulation, additional refinement is added whenever the gas (dark matter) mass in a cell exceeds () M. We allow up to 10 levels of refinement, resulting in a smallest cell size of 227 pc.

In order to explore the impact of some of the uncertain parameters in our model, simulations with different CR physics, including a no CR run, and CR runs with diffusion coefficients of , , , and cm s are done on the most massive of the selected halos. The results of these simulations are listed in Table 1. We analyze the results from these runs (see below) and choose the coefficient of cm s to carry on the simulations with the other selected halos. As discussed below, this choice is somewhat arbitrary, but it is consistent with estimates from recent GALPROP models (e.g., Ptuskin et al., 2006; Ackermann et al., 2012) and with observational measurements (e.g., Strong & Moskalenko, 1998; Tabatabaei et al., 2013). It is also found to be consistent with galactic gamma ray emission in Salem, Bryan & Corlies (2016). The properties of all selected halos simulated with this choice of parameters are listed in Table 2.

H1 174.32
H2 161.82
H3 132.10
H4 115.78
H5 118.60

R is in unit of kpc. M and M are in unit of M.

Table 2: Simulated dwarf galaxy properties

3 Results

First, we focus on comparisons of runs with different CR physics, using the most massive of the halos. Second, we discuss the physical properties of the suite of five halos with fixed CR physics. Third, we compare our simulation results to observations, specifically the Baryonic Tully-Fisher relation and HI stellar mass relation. Analysis has been carried out using the yt package (Turk et al., 2011).

3.1 Impact of Varying CR Physics

Figure 2: Circular velocity profiles of the neutral hydrogen gas in the disk from the simulations with non-zero CR diffusion, measured out to a limiting HI surface density of cm.

We choose an isolated halo of around M and study the effect of CR feedback by varying the CR physics and fixing all the other parameters (See Table 1 for the simulation parameters). Figure 1 shows edge-on images of the gas density, stellar density and CR energy density of the dwarf galaxies produced in each run at , focusing on just the center of the halo where the stellar system forms. Since H1-0 and H1-1 are not disk-shaped, we show them in the direction of their total angular momentum within 10 kpc radius.

The run without any CR feedback (H1-0; left-most panels) forms a diffuse, irregular cloud. We find that the energetic thermal feedback is actually remarkably effective in regulating the formation of dense, star-forming clumps. This behavior is quite different than that seen in high-mass galaxies using identical star formation and feedback prescriptions (e.g., Hummels & Bryan, 2012) and is more like the behavior observed in very low-mass galaxies Simpson et al. (e.g., 2013, 2015). We explore this point in more detail in Section 4.1

Moving on to the simulations with cosmic-rays, we see in Figure 1 that H1-1, with no CR diffusion, has an almost spherically symmetric gas distribution. This occurs because when the CRs are tied to the gas, they provide complete pressure support and rotational support is unimportant. Clearly this is unrealistic (and probably violates gamma-ray emission from dwarf irregulars – see Salem, Bryan & Corlies (2016)).

Figure 3: Comparison of the star formation rate history in simulations of the same halo but different CR physics, as noted in the legend.

The other runs, which have non-zero CR diffusion, all produce disk-shaped galaxies, generally becoming thinner as the diffusion coefficient increases. H1, which has the highest diffusion coefficient of cm s, results in the most extended disk. H1-1, which has the lowest non-zero diffusion coefficient, results in a smaller disk which is surrounded by denser gas. The gas density and CR energy density immediately surrounding the disk also decreases as increases. The CRs escape more easily, both because they diffuse out more rapidly and because they drive stronger outflows (Salem & Bryan, 2014), and therefore they provide less pressure support, decreasing the weight of gas they can support in the diffuse (not-rotating) halo.

Figure 4: Edge-on views of the simulated galaxies with different halo masses but a fixed diffusion coefficient ( cm s. The plot show halos H1, H2, H3, H4, H5 from left to right, and gas density, stellar density, and CR energy density from top to bottom. Each panel represents a kpc kpc area at .

Since all the runs with non-zero CR diffusion produce disk galaxies, we can examine the rotational velocity of the gas in the disk. Figure 2 shows the rotational velocities of the HI gas in H1-2, H1-3, and H1. To compute these rotation velocity curves, we look at the actual rotational velocities of the HI gas (rather than a simple dynamical measure such as ). To do this, we adopt cm as the lowest detectable HI surface density in order to make an approximate comparison to observations – the results are not strongly sensitive to this threshold, but disks from runs with low diffusion coefficients start to pick up some CR-pressure supported gas and show falling rotation curves. In particular, in H1-2, the run with a diffusion coefficient of 3 10 cm s is affected by this issue: the circular velocity peak at 160 km/s at about 1.5 kpc radius and falls quickly beyond that radius. As can be seen from its edge-on plot, the disk in H1-2 has a small radius, but its surrounding gas is relatively dense, which makes the HI detectable radius even larger than the disk radius. In that run, outside the disk, the gas is actually supported in part by the pressure of the CRs, which explains why the circular velocity curve drops so quickly after the peak.

From this figure, we see that, at least for the two higher diffusion coefficient runs, the circular velocity curves are relatively flat, in reasonable agreement with observed systems. For example, in the run with the highest diffusion coefficient (H1), the circular velocity is relatively flat at about 140 km/s out to 6-7 kpc.

Next, we examine the impact of CR pressure on the star formation rate. This is shown in Figure 3 for all the runs of halo H1 with different CR physics. For all of our more realistic runs, with non-zero diffusion coefficient, the SFRs are quite similar: an early rapid rise follow by a nearly constant rate that falls slowly from 1 to 0.5 M/yr. The run without diffusion shows a rate which is a factor of two larger, due to the higher gas densities in its core, while the thermal-only feedback is lower by nearly an order of magnitude and shows considerably more variation with time – consistent with the bursty feedback seen in the gas distribution.

3.2 Varying Halo Mass with Fixed CR Physics

In this section, we fix the CR diffusion coefficient to be cms and study the impact on dwarf galaxies with different halo masses, varying the total mass from about M to M. We examine how CR feedback behaves in different scales and also use this suite of runs to compare with the observational scaling relations in Section 3.3.

All of these runs produce dwarf galaxies with extended disk features; Figure 4 shows edge-on images of these galaxies. The halo masses decrease from left to right: H1 to H5 as described in Table 2. There is a fairly clear trend of decreasing disk extent in the stellar distribution, with the larger halos having more extended disks. H2 is a slight outlier in this trend with a truncated disk, particularly in the gas and CR components. In addition, the circumgalactic medium is less dense (with a lower CR energy density) for the lower mass halos.

Figure 5: Measured circular velocity of the neutral hydrogen in the disk for halos H1 to H5. As in Figure 2, we plot the circular velocity only out to a limiting HI surface density of cm.

We plot the rotational velocities of the gas for these five galaxies in Figure 5. Within our adopted HI detection limit, the curves of H3, H4, and H5, stay flat at about 100 km/s, just above 80 km/s, and just below 80 km/s, respectively. The velocities are consistent with their masses, and although we do not make a detailed comparison the curves are similar to observed rotation curves – in particular, they do not posses large bulges with centrally peaked rotation curves, as seen in many simulations without effective feedback. H1 has a higher rotational velocity (not surprising given its larger mass) and shows some mild features at 1 and 5 kpc. On the other hand, the curve of H2 looks quite different: it is highly peaked, and falls sharply at about 4 kpc (again because of partial CR pressure support). As can be seen from the edge-on plot of H2, this dwarf galaxy does have a dense core, larger than the other dwarfs.

Figure 6 shows the star formation histories of these same halos. There is a clear progression of increasing SFR with increasing mass. Despite some bumps and wiggles, the overall histories are remarkably flat, with each profile showing typically a factor of 2-3 variation over 12-13 Gyr. There is a general trend for slowly falling star formation rate with time, particularly evident in H3. This occurs because of the relatively tight self-regulation – feedback acts to decrease the gas density, slowing star formation.

Figure 6: Star formation histories of our simulated halos H1 to H5 (all with fixed CR diffusion coefficient).

3.3 Comparison with Observed Scaling Relations

Dwarf galaxies are less efficient in star formation than normal galaxies, and the baryon fraction of galaxies decreases as the total mass goes down (e.g., McGaugh et al., 2010). In this section, we compare our simulation results with a range of observational probes, including the baryonic Tully-Fisher relation, as well as observed gas-to-star ratios and inferred stellar mass-to-halo mass relations.

We begin with the observed relation between the disk mass (the combined stellar and gas masses) and the flat part of the rotational velocity of the disk – the baryonic Tully-Fisher relation. Figure 7 shows the simulation results from varying halo mass with fixed CR physics overlaid on the fitted line of the Baryonic Tully-Fisher relation from McGaugh et al. (2010). They provide a fit for disk galaxies in the range of to km/s given by M with and . In Figure 7, the solid line represents their best fit, and the grey area represents one standard deviation of . We estimated the scatter in the observed as the average for spiral and gas disk galaxies in Table 2 of McGaugh et al. (2010).

For the simulations, we determined the baryonic mass as the sum of the stellar and gas masses, with the gas mass calculated as times the HI mass. Although the code evolves the HI density, it is underestimated in dense gas due to the lack of self-shielding. To correct for this, the HI density is estimated as follows: when the total H number density cm, we assume that the neutral fraction is nearly unity and take ; otherwise we assume the HI density is zero. This is approximately consistent with radiative transfer calculations (Rahmati et al., 2013) that take into account the metagalactic radiation background (a local radiation background would decrease the neutral gas fraction even lower, and so our bounds here are conservative).

Both the stellar mass and HI mass are calculated within the same HI column density limit we used to determine the radial extent for the rotation curves ( cm) and extend above and below the plane by 2 kpc. The circular velocity is calculated as the average rotation velocity of the gas at a radius of 2.2 scale lengths of the stellar surface density. To determine the scale length, we fit the stellar surface density curve with a straight line from 2 kpc to 10 kpc.

Our simulation results all lie within the grey area and scatter around the observed relation. Although we do not have a significantly large numerical sample to make a detailed comparison, these preliminary results appear to be in good agreement, indicating that the baryonic content of our galaxies is consistent with that observed. We note that we cannot include the purely thermal simulation (H1-0) on this plot as the systems has no systematic rotation.

Figure 7: Comparison between the five dwarf simulations (H1-H5) run with our standard diffusion coefficient and the baryonic Tully-Fisher relation from observations. The solid line is the fit from McGaugh et al. (2010). The red crosses are simulation data. The grey area shows the region of one sigma error in . The dashed line is taken from Komatsu et al. (2009) and indicates the relation between the halo circular velocity and the cosmic baryon fraction of the halo.

The baryonic Tully Fisher relation compares the galaxies’ baryonic mass with a measure of its total mass. We would also like to compare the two main components of baryons: neutral hydrogen and stars. Figure 8 shows the ratio of the HI to stellar components as a function of the stellar mass, all quantities which are directly observable. The observed relation found in Evoli et al. (2011) (their Eq. 5) can be approximated as:


In Figure 8, the solid line shows this relation. We also estimate the observed scatter around this relation by taking a variation of +/- 0.5 log M/M, based on a visual inspection of Figure 1 of Evoli et al. (2011); this is shown as a grey area in our Figure 8. The ratio from the simulations is also plotted, again using the stellar and neutral hydrogen masses within our adopted HI surface density limit, and +/- 2 kpc in vertical extent.

As can be seen from Figure 8, all of our simulation results are below the observed relation, indeed most lie below the estimated scatter in the observations (the grey band). This demonstrates that our simulation produces consistently lower HI fraction than observed. There is some indication that the simulated relation shows a rising gas fraction with decreasing stellar mass; however, the sample size is too small to be sure.

This mismatch in the gas to star ratio, by approximately a factor of five, when combined with the agreement in the baryonic Tully-Fisher relation (which measures the total baryonic content), implies that our galaxies contain the correct amount of baryons, but are overly efficient in their conversion of gas to stars. One alternative way to examine just the stellar content of our systems is based on the dark matter abundance matching measurements in Evoli et al. (2011), which provides a relation between the stellar mass and halo mass under the assumption that there is a monotonic relation between halo mass and stellar mass. Indeed, if we compare our results to their inferred M-M relation, we find that our values are systematically too large by factors of 2-5, consistent with the idea that the division between stellar and gas mass is incorrect in these simulations. Finally, we note that these relations also imply that the simulations disagree with the observed Kenicutt-Schmidt relation between surface gas density and star formation surface density.

Figure 8: Comparison between this work and the relation of HI mass and stellar mass inferred from observations. The solid line is derived in Evoli et al. (2011). The grey area shows an estimated scatter of +/- 0.5 log around the observed relation in the vertical direction. The red crosses are the simulation results.

3.4 Dark Matter Profiles

There are a number of observational indications that dwarf galaxies do not contain the cuspy dark-matter profiles predicted in dark-matter only simulations (e.g., de Blok et al., 2008; Walker & Peñarrubia, 2011; Boylan-Kolchin, Bullock & Kaplinghat, 2011, 2012; Salucci et al., 2012; Klypin et al., 2014). A number of simulations with strong feedback have found that it is possible for the feedback to decrease the dark matter density in the center of halos (e.g., Navarro, Eke & Frenk, 1996; Governato et al., 2010; Pontzen & Governato, 2012; Ogiya & Mori, 2014).

Therefore, we examine our simulation to see the impact of thermal vs. CR feedback on the core properties of the systems. Figure 9 shows the dark matter density profile of our most massive halo in the simulations with different CR physics. Interestingly, all of the CR runs show very similar results, indicating that although they produce different distributions for the gas and stars, the presence of CR feedback leads to similar dark matter core properties. In all cases, the profile is quite steep, with an inner density profile between and (similar to dark-matter only runs). On the other hand, the run without CR diffusion shows a clear core in the dark matter with a core radius of several kpc.

The cored dark matter profile in the no-CR case might be due to the large and rapid motions of the gas, which can disturb the dark matter distribution (Pontzen & Governato, 2012) and is consistent with the observed asymmetric gas distribution and the bursty nature of the feedback in the no-CR run.

Figure 9: Dark matter profiles of simulations with different CR physics.

4 Discussion

We have performed a preliminary investigation into the impact of CRs on the formation and evolution of dwarf galaxies with rotation velocities in the 70-140 km/s range using cosmological simulations. This complements similar work (Salem, Bryan & Hummels, 2014) which looked at a higher mass galaxy (with a rotation velocity above 200 km/s). As in that work, we found that the inclusion of CR feedback had a strong impact on the gas and stellar properties. In both cases, CR feedback led to robust outflows and rotationally supported disks with relatively flat rotation curves, but only if CR diffusion was included. A comparison to observations showed good agreement with the baryonic Tully-Fisher relation, indicating that the total disk mass was reasonable. This mass, converted into a fraction of the total halo mass, is well below the cosmic value, indicating that feedback was important in generating mass-loaded winds. On the other hand, a comparison of the gas content of the simulated galaxies against observations showed that too much gas had been converted to stars. In addition, the gas depletion times of our galaxies (computed as ), are a few Gyr, well below observed values for low-mass isolated dwarfs (Huang et al., 2012).

Our star formation and feedback model is approximate, with a number of free parameters required to describe physics that is not fully modeled. In section 3.1, we explored the impact of some of these for halo H1. Remarkably, we find that the CR diffusion rate has a relatively small impact on the stellar mass produced (see Table 1), or the gas mass. The exception was the run with , which produced a substantially higher stellar mass and did not produce a disk (this combination does not generate winds).

We also carried out an additional run of halo H1 with the standard diffusion rate (and other parameters), but a reduced CR energy efficiency (10% instead of 30%) and find that relative to the fiducial run, it produces 50% more stars. The resulting disk rotates marginally faster such that the system stays close to the observed baryonic Tully-Fisher relation. However, the HI to gas mass ratio is even lower the fiducial run, resulting in an an even stronger disagreement with the observed relation. In principle, a higher CR efficiency might bring this into agreement with observations, but that would disagree with local observational and theoretical work on CR acceleration efficiencies.

Another parameter which can be varied is the star formation efficiency per free-fall time, which we took to be 1%. We carried out another simulation much like H1 but with a star formation efficiency four times lower (with cm s; this run was not included in figures and tables above). This had almost no effect on the resulting stellar mass or gas mass, consistent with the idea that the disk is in a state of self-regulation – indeed, we verified visually that the lower efficiency resulted in a higher density gas disk and hence essentially the same net star formation rate at (since the star formation rate increases as ).

There are, of course, other parameters to vary, such as the IMF, the fraction of energy in SN, as well as other uncertainties in the subgrid model. In addition, there is some indication from H2 and the high-mass halo explored in Salem, Bryan & Hummels (2014) that an additional form of feedback is required in high-density regions – possibly radiative heating or pressure.

Finally, we briefly discuss our results in comparison to (Uhlig et al., 2012), who recently carried out simulations of dwarf galaxies with cosmic rays. The implementation differs somewhat in that they adopted a model in which CRs stream along magnetic field lines and heat the gas (rather than isotropically diffuse, as in our model). Although they do not simulate cosmological volumes, their halo mass which most closely compares to ours produces approximately similar results: robust cosmic-ray driven outflows with a modest mass-loading factor. Although we are unable to compute a mass-loading factor, the large stellar mass to halo mass we find implies a relatively low mass-loading factor.

4.1 Feedback efficiency of the purely thermal run

Figure 10: The distribution of SN heating and radiative cooling times for all gas particles within 30 kpc of the center of our thermal-only simulation (H1-0) at . The black line indicates equality; since nearly all gas cells lie above this line, we see that pure thermal heating operates more rapidly than cooling. See text for the definition of the SN heating time.

Finally, we turn to a discussion of the thermal-only feedback run (H1-0). Interestingly, the feedback in this case was actually more effective than the CR+thermal feedback, in contrast to the result for the higher mass halo with the same code and CR parameters discussed in Salem, Bryan & Hummels (2014). One possibility is that the purely thermal feedback is more efficient than in more massive galaxies because the shallower potential well results in low gas densities, lowering the radiative cooling rate and permitting blast waves from SN feedback to accelerate the gas.

We explore the relative importance of cooling and heating in Figure 10, which shows a two-dimensional distribution of radiative cooling times and a measure of the possible supernova heating rate for all gas within 30 kpc of the center of this halo as . We compute the SN heating rate in a way similar to that in Simpson et al. (2013), where we estimate the timescale required for our thermal energy injection mechanism to heat a cell to (the precise temperature selected does not affect our conclusions):


where is the gas mass cell in the cell, , M, our typical stellar particle mass, Myr, the timescale over which a star particle returns its supernova energy to the ISM in our model, and , as defined in the methods section. We stress that not all cells are being heated at this rate, but it is the rate that a single star particle would typically heat its local gas and so is a good measure of the timescale over which feedback heating acts.

From this figure, we see that nearly all the gas in the simulation lies above the black line indicating equality. This demonstrates that the gas densities are sufficiently low that pure thermal feedback is efficient in heating the gas and driving outflows – we focus on density because the cooling rate is proportional to the density squared. A close examination of this figure indicates that a substantial amount of gas lies close to the line of equality (in fact, this gas is in the disk), and so we expect that the larger potential depths of higher mass halos will push gas beyond the line into a regime where purely thermal supernovae feedback cannot heat the gas. This is consistent with the finding that our dwarf galaxy simulations do drive outflows while higher mass (e.g. Milky Way mass) systems do not. Note that this statement depends on details of the simulations (e.g. see definition of ) and so may be a reflection of numerical models rather than physical quantities.

This demonstrates that purely thermal feedback is sufficient to drive feedback but it does not answer the question why adding cosmic rays decreases the bursty outflows, as seen in Figure 1. The answer appears to be that cosmic rays actually increase the pressure in the diffuse circumgalactic gas surrounding galaxies. The purely thermal run has a (thermal only) CGM pressure at a radius of 10 kpc of approximately erg cm, while even the CR run with the highest diffusion coefficient (H1), has a (mostly CR) pressure of nearly erg cm. This larger pressure resists the impulsive outflows driven by bursts of supernovae in the disk (although it does not appear to eliminate the slower gas outflows driven by the CR pressure gradient).

Finally, we note a consequence of the effective thermal only feedback, which is that it resulted in bursty star formation at a low-level; and was so effective that the resulting stellar system was “hot” – that is, it was supported by velocity dispersion rather than rotation (an issue that has been raised in the context of other simulations – see the discussion in Wheeler et al. (2015).

5 Summary

We simulate a series of dwarf galaxies with zoom-in cosmological simulations, and study the influence of CR feedback. First we compare the outcome by varying the CR physics in different runs of the same galaxy halo. The different CR physics include feedback without CR, CR feedback with zero diffusion, CR feedback with equal to , , and cm s. Then we fixed the CR diffusion at cm s and simulated five dwarf galaxies with different masses, ranging from 8 to 30 M. We summarize the results of this work as follows:

  1. Adding CRs and some realistic level of diffusion consistently produced thin, extended disk galaxies with nearly flat rotation curves, in contrast to a case with purely thermal feedback, which produced a hot stellar systems. The star formation rate in the CR systems was relatively constant, compared to the lower (and burstier) SFR in the purely thermal feedback model.

  2. Simulations of our five dwarf galaxies with different masses but the same CR model matched well with the observed baryonic Tully-Fisher relation. However, the five galaxies’ neutral hydrogen fractions were always lower than what is observed, by a factor of about five and the stellar content was larger than observed (but the total baryonic content matched observations).

  3. Simulations with CR feedback produce cuspy dark matter profiles, while our purely thermal feedback case resulted in a cored dark matter profile.

These results indicate that the impact of cosmic rays in dwarf galaxies is potentially quite important, a conclusion which is in general agreement with previous work (e.g., Jubelgas et al., 2008). In detail, the results are somewhat dependent on the physical model, but broadly speaking CRs with diffusion produce dwarf galaxies which are rotationally supported and have relatively flat rotation curves. This result appears to be quite robust (and occurs even though 70% of the feedback energy is in the form of thermal energy in these models) and is probably due to the smooth pressure provided by CRs (along with the outflows they drive). The baryonic content of the dwarf systems with CRs is in broad agreement with observations.

On the other hand, there are two areas of disagreement with observations in these models: (a) the detailed census of baryons in each system (stars vs. gas) is much too heavily weighted to the stellar side, and (b) the presence of dark matter cores may be a problem (although CR pressure may systematically support the gas at small radius and so bias the observed measures of dark matter on small scales). Clearly, more work is required to explore and refine these results. For example, it would be useful to explore more realistic models including magnetic fields, anisotropic CR diffusion, and the impact of field-CR streaming instabilities.


We thank the anonymous referee for a useful report which improved the clarity of the paper. We would like to thank Cameron Hummels, Mordecai-Mark Mac Low and Mary Putman for useful discussion related to this work. We acknowledge financial support from NSF grants AST-1312888, and NASA grants NNX12AH41G and NNX15AB20G, as well as computational resources from NSF XSEDE, and Columbia University’s Yeti cluster. Computations described in this work were performed using the publicly-available Enzo code (, which is the product of a collaborative effort of many independent scientists from numerous institutions around the world. Their commitment to open science has helped make this work possible.


  • Ackermann et al. (2012) Ackermann M. et al., 2012, \apj, 750, 3
  • Behroozi, Wechsler & Conroy (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, \apj, 770, 57
  • Boylan-Kolchin, Bullock & Kaplinghat (2011) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, \mnras, 415, L40
  • Boylan-Kolchin, Bullock & Kaplinghat (2012) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2012, \mnras, 422, 1203
  • Bradford, Geha & Blanton (2015) Bradford J. D., Geha M. C., Blanton M. R., 2015, \apj, 809, 146
  • Brooks & Zolotov (2014) Brooks A. M., Zolotov A., 2014, \apj, 786, 87
  • Bryan & Norman (1997) Bryan G. L., Norman M. L., 1997, ArXiv e-print, astro-ph/971087
  • Bryan et al. (2014) Bryan G. L. et al., 2014, \apjs, 211, 19
  • Cen & Ostriker (1992) Cen R., Ostriker J. P., 1992, \apjl, 399, L113
  • de Blok & Bosma (2002) de Blok W. J. G., Bosma A., 2002, \aap, 385, 816
  • de Blok et al. (2008) de Blok W. J. G., Walter F., Brinks E., Trachternach C., Oh S.-H., Kennicutt, Jr. R. C., 2008, \aj, 136, 2648
  • Dekel & Silk (1986) Dekel A., Silk J., 1986, \apj, 303, 39
  • Di Cintio et al. (2014) Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., Dutton A. A., Wadsley J., 2014, \mnras, 437, 415
  • Drury (1985) Drury L. O., 1985, in Cosmical Gas Dynamics, Kahn F. D., ed., pp. 131–136
  • Drury & Falle (1986) Drury L. O., Falle S. A. E. G., 1986, \mnras, 223, 353
  • Ellison et al. (2010) Ellison D. C., Patnaude D. J., Slane P., Raymond J., 2010, \apj, 712, 287
  • Evoli et al. (2011) Evoli C., Salucci P., Lapi A., Danese L., 2011, \apj, 743, 45
  • Geha et al. (2012) Geha M., Blanton M. R., Yan R., Tinker J. L., 2012, \apj, 757, 85
  • Girichidis et al. (2016) Girichidis P. et al., 2016, \apjl, 816, L19
  • Governato et al. (2010) Governato F. et al., 2010, \nat, 463, 203
  • Governato et al. (2007) Governato F., Willman B., Mayer L., Brooks A., Stinson G., Valenzuela O., Wadsley J., Quinn T., 2007, \mnras, 374, 1479
  • Governato et al. (2012) Governato F. et al., 2012, \mnras, 422, 1231
  • Guedes et al. (2011) Guedes J., Callegari S., Madau P., Mayer L., 2011, \apj, 742, 76
  • Guo et al. (2010) Guo Q., White S., Li C., Boylan-Kolchin M., 2010, \mnras, 404, 1111
  • Haardt & Madau (2012) Haardt F., Madau P., 2012, \apj, 746, 125
  • Hahn & Abel (2011) Hahn O., Abel T., 2011, \mnras, 415, 2101
  • Hanasz et al. (2013) Hanasz M., Lesch H., Naab T., Gawryszczak A., Kowalik K., Wóltański D., 2013, \apjl, 777, L38
  • Hopkins et al. (2014) Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, \mnras, 445, 581
  • Huang et al. (2012) Huang S., Haynes M. P., Giovanelli R., Brinchmann J., Stierwalt S., Neff S. G., 2012, \aj, 143, 133
  • Hummels & Bryan (2012) Hummels C. B., Bryan G. L., 2012, \apj, 749, 140
  • Jubelgas et al. (2008) Jubelgas M., Springel V., Enßlin T., Pfrommer C., 2008, \aap, 481, 33
  • Jun, Clarke & Norman (1994) Jun B.-I., Clarke D. A., Norman M. L., 1994, \apj, 429, 748
  • Kang & Jones (2006) Kang H., Jones T. W., 2006, Astroparticle Physics, 25, 246
  • Klypin et al. (2014) Klypin A., Yepes G., Gottlober S., Prada F., Hess S., 2014, ArXiv e-print, arXiv:1411.4001
  • Komatsu et al. (2009) Komatsu E. et al., 2009, \apjs, 180, 330
  • Kulsrud (2005) Kulsrud R. M., 2005, Plasma physics for astrophysics
  • Leauthaud et al. (2012) Leauthaud A. et al., 2012, \apj, 744, 159
  • McGaugh et al. (2010) McGaugh S. S., Schombert J. M., de Blok W. J. G., Zagursky M. J., 2010, \apjl, 708, L14
  • Moore (1994) Moore B., 1994, \nat, 370, 629
  • Moore et al. (1999) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, \apjl, 524, L19
  • Navarro, Eke & Frenk (1996) Navarro J. F., Eke V. R., Frenk C. S., 1996, \mnras, 283, L72
  • Navarro & White (1994) Navarro J. F., White S. D. M., 1994, \mnras, 267, 401
  • Oñorbe et al. (2015) Oñorbe J., Boylan-Kolchin M., Bullock J. S., Hopkins P. F., Kereš D., Faucher-Giguère C.-A., Quataert E., Murray N., 2015, \mnras, 454, 2092
  • Ogiya & Mori (2014) Ogiya G., Mori M., 2014, \apj, 793, 46
  • Oppenheimer & Davé (2008) Oppenheimer B. D., Davé R., 2008, \mnras, 387, 577
  • O’Shea et al. (2004) O’Shea B. W., Bryan G., Bordner J., Norman M. L., Abel T., Harkness R., Kritsuk A., 2004, ArXiv eprint, astro-ph/0403044
  • Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, \aap, 571, A16
  • Pontzen & Governato (2012) Pontzen A., Governato F., 2012, \mnras, 421, 3464
  • Ptuskin et al. (2006) Ptuskin V. S., Moskalenko I. V., Jones F. C., Strong A. W., Zirakashvili V. N., 2006, \apj, 642, 902
  • Rahmati et al. (2013) Rahmati A., Schaye J., Pawlik A. H., Raičevic̀ M., 2013, \mnras, 431, 2261
  • Salem & Bryan (2014) Salem M., Bryan G. L., 2014, \mnras, 437, 3312
  • Salem, Bryan & Corlies (2016) Salem M., Bryan G. L., Corlies L., 2016, \mnras, 456, 582
  • Salem, Bryan & Hummels (2014) Salem M., Bryan G. L., Hummels C., 2014, \apjl, 797, L18
  • Salucci & Burkert (2000) Salucci P., Burkert A., 2000, \apjl, 537, L9
  • Salucci et al. (2012) Salucci P., Wilkinson M. I., Walker M. G., Gilmore G. F., Grebel E. K., Koch A., Frigerio Martins C., Wyse R. F. G., 2012, \mnras, 420, 2034
  • Sawala et al. (2011) Sawala T., Guo Q., Scannapieco C., Jenkins A., White S., 2011, \mnras, 413, 659
  • Sawala et al. (2010) Sawala T., Scannapieco C., Maio U., White S., 2010, \mnras, 402, 1599
  • Schaye et al. (2015) Schaye J. et al., 2015, \mnras, 446, 521
  • Simpson et al. (2015) Simpson C. M., Bryan G. L., Hummels C., Ostriker J. P., 2015, \apj, 809, 69
  • Simpson et al. (2013) Simpson C. M., Bryan G. L., Johnston K. V., Smith B. D., Mac Low M.-M., Sharma S., Tumlinson J., 2013, \mnras, 432, 1989
  • Smith, Sigurdsson & Abel (2008) Smith B., Sigurdsson S., Abel T., 2008, \mnras, 385, 1443
  • Sparre et al. (2015) Sparre M. et al., 2015, \mnras, 447, 3548
  • Springel et al. (2008) Springel V. et al., 2008, \mnras, 391, 1685
  • Stinson et al. (2007) Stinson G. S., Dalcanton J. J., Quinn T., Kaufmann T., Wadsley J., 2007, \apj, 667, 170
  • Strong & Moskalenko (1998) Strong A. W., Moskalenko I. V., 1998, \apj, 509, 212
  • Tabatabaei et al. (2013) Tabatabaei F. S. et al., 2013, \aap, 552, A19
  • Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, \apjs, 192, 9
  • Uhlig et al. (2012) Uhlig M., Pfrommer C., Sharma M., Nath B. B., Enßlin T. A., Springel V., 2012, \mnras, 423, 2374
  • Walker & Peñarrubia (2011) Walker M. G., Peñarrubia J., 2011, \apj, 742, 20
  • Walter et al. (2008) Walter F., Brinks E., de Blok W. J. G., Bigiel F., Kennicutt, Jr. R. C., Thornley M. D., Leroy A., 2008, \aj, 136, 2563
  • Wefel (1987) Wefel J. P., 1987, Origin and transport of high energy particles in the galaxy. Tech. rep.
  • Wheeler et al. (2015) Wheeler C., Oñorbe J., Bullock J. S., Boylan-Kolchin M., Elbert O. D., Garrison-Kimmel S., Hopkins P. F., Kereš D., 2015, \mnras, 453, 1305
  • White & Frenk (1991) White S. D. M., Frenk C. S., 1991, \apj, 379, 52
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description