Impact of Cosmic-Ray Transport on Galactic Winds

Galactic winds are observed ubiquitously in star-forming galaxies and significantly affect their chemical and dynamical evolution (Veilleux et al. 2005). Galactic winds redistribute angular momentum, aiding in the formation of extended disks (Brook et al. 2011; Übler et al. 2014); help produce large-scale magnetic fields in dwarf galaxies (Moss & Sokoloff 2017); and pollute the intergalactic medium with metals (Steidel et al. 2010; Booth et al. 2012). Additionally, most galaxies are missing a large fraction of baryons compared to the cosmological average (Bell et al. 2003). Models matching observed luminosity functions to simulated halo mass functions find that 20% of the baryons are accounted for in L_* galaxies, and that this fraction decreases rapidly for both more and less massive galaxies (Guo et al. 2010). This suggests that the efficiency of converting baryons into stars is a strong function of halo mass.The discrepancies between halo and stellar properties, "the missing baryons problem," constitutes an outstanding challenge in galaxy formation. Galactic winds can possibly solve the missing baryons problem by ejecting baryons out of galaxies. For galaxies more massive than L_*, active galactic nuclei likely dominate the energetics of the outflows (e.g., Croton et al. 2006), while in less massive galaxies galactic winds are likely driven by stellar feedback (Larson 1974, Chevalier & Clegg 1985; Dekel & Silk 1986).

In the standard model of supernova-driven galactic winds (Chevalier & Clegg 1985), thermal energy is injected into the gas, launching it ballistically and entraining denser gas as it is flung out of the galaxy. Thermally driven winds may explain the superwinds observed in starburst galaxies, such as M82 (Bustard et al. 2016). However, results from high-resolution simulations demonstrated that such models might inject significant amounts of energy and launch metals out of galaxies but fail to expel a significant amount of mass into the intergalactic medium (Mac Low & Ferrara 1999; Melioli et al. 2013). Additionally, Steidel et al. (2010) found that the kinematic features of Lyman-break galaxies best match models in which the gas velocity increases with distance to at least 100 kpc. This result is also difficult to reconcile with the thermal feedback model. The insufficiencies of purely thermally driven winds hint at the importance of additional stellar feedback processes, such as cosmic rays (Boulares & Cox 1990; Breitschwerdt et al. 1993; Uhlig et al. 2012).

Cosmic rays can be accelerated by means of the diffusive acceleration mechanism operating in the shocks of supernova remnants (Blandford & Eichler 1987; Caprioli 2015) and in the winds from massive stars (Bykov 2014). Cosmic rays exert pressure that is in rough equipartition with magnetic and dynamical pressures in the interstellar medium (ISM; Zweibel & Heiles 1997; Beck 2001; Cox 2005), suggesting their dynamical importance. In particular, cosmic rays can provide pressure support against self-gravitating clouds, suppressing star formation (Jubelgas et al. 2008; Pfrommer et al. 2017). Additionally, Fermi γ-ray observations of the starburst galaxies M82 and NGC 253 imply cosmic-ray energy densities roughly two magnitudes higher than in the Milky Way (Paglione & Abrahams 2012; Yoast-Hull et al. 2013, 2014).

Cosmic rays escape the Galactic disk in ∼10 Myr (Strong et al. 2007). Compared to the thermal gas, cosmic rays can be relatively free of energy losses, which together with their fast escape from the disk, suggests that cosmic rays may efficiently transport supernova energy to regions occupied by tenuous gas above the disk, which they can accelerate into a wind (Hanasz et al. 2013).

Early work by Ipavich (1975) considered the emission of magnetohydrodynamic (MHD) waves by super-Alfvénic cosmic rays. These waves enable cosmic rays to be coupled to the thermal gas. A steady-state, spherically symmetric, hydrodynamic treatment by Ipavich (1975) suggested that cosmic rays could drive outflows at rates ≳1 M_⊙ yr⁻¹ from a typical galaxy. Breitschwerdt et al. (1991) extended the work by Ipavich (1975) by including the streaming of cosmic rays along large-scale magnetic fields and found that mass outflow rates of ∼1 M_⊙ yr⁻¹ are possible in Milky Way-like galaxies. Everett et al. (2008) further extended this work by combining cosmic-ray and thermal pressure under Milky Way conditions and found that both thermal and cosmic-ray pressures were essential for wind driving in the Milky Way.

Recently, 3D numerical studies of cosmic-ray winds have found that wind properties depended sensitively on the details of cosmic-ray transport. This was demonstrated in both Eulerian grid hydrodynamic (Uhlig et al. 2012; Booth et al. 2013; Salem & Bryan 2014) and MHD simulations (Hanasz et al. 2013; Ruszkowski et al. 2017) as well as unstructured moving mesh simulations (Pakmor et al. 2016a, 2016b; Simpson et al. 2016; Pfrommer et al. 2017; Jacob et al. 2018).

In predominantly cold, neutral gas, cosmic rays can propagate faster than in the ionized medium, and the effective transport can be substantially larger; i.e., cosmic rays can decouple from the gas. In this work, we study the consequences of this decoupling and show that it has a significant impact on the properties of cosmic-ray-driven galactic winds. In Section 2, we delineate the numerical methods and treatment of physics in our simulations. In Section 3, we present our results, and in Section 4, we conclude.

We model cosmic rays with a two-fluid model (e.g., Salem & Bryan 2014; Ruszkowski et al. 2017), in which cosmic rays take the form of an ultrarelativistic ideal fluid with an adiabatic index γ_cr = 4/3, and the thermal gas is characterized by an adiabatic index γ = 5/3. We include the advection of cosmic rays, dynamical coupling between cosmic rays and the gas, and model the transport of cosmic rays relative to the gas via anisotropic diffusion of cosmic rays along magnetic field lines rather than via the streaming instability (see Section 2.3). We model the effect of cosmic rays decoupling from the cold, neutral ISM via a temperature-dependent diffusion coefficient (see Section 2.3). Additionally, we include star formation and feedback, self-gravity of the gas, and radiative cooling. We solve the following equations:

where ρ is the gas density, u is the gas velocity, is a density sink term representing star formation, is a density source term representing stellar winds and supernovae (see Section 2.4), B is the magnetic field, p_tot is the sum of the gas (p_th), magnetic, and cosmic-ray (p_cr) pressures, g = −∇ϕ + g_NFW is the gravitational acceleration (including contributions from self-gravity of gas and stars: −∇ϕ and dark matter: g_NFW; see Section 2.1), is the momentum injection from stellar winds and supernovae, e = 0.5ρu² + e_g + e_cr + B²/8π is the total energy density (where e_g is the thermal energy density), is the temperature-dependent diffusion coefficient (see Section 2.3), e_cr is the cosmic-ray energy density, C is the radiative cooling term (see Section 2.2), H_SN is the supernova heating term (see Section 2.4), and ρ_b is the total (gas and stars) baryon density.

We use the adaptive-mesh refinement MHD code FLASH4.2 (Fryxell et al. 2000; Dubey et al. 2008) extended to include cosmic rays (Yang et al. 2012, 2013; Ruszkowski et al. 2017; Yang & Ruszkowski 2017) to solve the above equations. We utilize the directionally unsplit staggered mesh (USM) solver (Lee & Deane 2009; Lee 2013). The USM solver is a finite-volume, high-order Godunov scheme that utilizes constrained transport to ensure divergence-free magnetic fields.

Due to computational constraints, we employ sound-speed limiting. That is, we set a ceiling on the thermal and cosmic-ray energy so that the time step does not become unfeasibly small. Specifically, we impose an upper limit on the generalized sound speed,

In our fiducial runs, we limit c_s to km s⁻¹, but we have additionally run test cases for km s⁻¹ and found no significant differences in the star formation rates (c.f., Salem & Bryan 2014).

For the most expensive simulation, which employs the decoupling mechanism (Run DEC; cf., Sections 2.3 and 3), we use a raised density floor of 10⁻³ cm⁻³, which effectively limits the Alfvén speed. In test cases that use reduced diffusion coefficients (which are computationally easier), we found no difference in the results between runs that used the raised density floor of 10⁻³ cm⁻³ compared to the fiducial density floor of 5 × 10⁻⁷ cm⁻³.

2.1. Gravity

We include self-gravity of baryons (gas and stars) and solve the Poisson equation using the Barnes–Hut tree solver (Barnes & Hut 1986) implemented in FLASH4.2 by Richard Wünsch (Wünsch et al. 2018). This solver allows us to use mixed boundary conditions (see Section 2.5).

In addition to self-gravity, we include acceleration in the z-direction due to dark matter. This component of the gravitational field assumes that dark matter is distributed according to the Navarro–Frenk–White profile (Navarro et al. 1997) and has the following form:

where z is the height from the midplane, G is the universal gravitational constant, M₂₀₀ is the virial mass of the halo, , c is the halo concentration parameter, and r₂₀₀ is the virial radius defined as the radius of a sphere within which the average density exceeds the critical density at redshift zero by a factor of 200, i.e., . See Table 1 for the parameter values used in our simulations.

Table 1.  Model Parameters

Halo
	1012M⊙
c(2)	12
Disk
ρo(3)	5.24 g cm−3
	0.325 kpc
	100 M⊙ pc−2
	104 K
	1 μG
Star Formation
	10 cm−3
	300 K
	105 M⊙
	0.05
Stellar Feedback
	0.25
	0.1
	1051 erg/(Msfc2)
	100 M⊙

Note. From top to bottom, the rows contain the (1) halo mass, (2) concentration parameter, (3) initial midplane density, (4) initial scale height of the gas disk, (5) initial gas surface density, (6) initial temperature, (7) initial magnetic field strength, (8) gas density threshold for star formation, (9) floor temperature, (10) minimum stellar population particle mass, (11) star formation efficiency, (12) fraction of stellar mass returned to the ISM, (13) fraction of supernova energy bestowed upon cosmic rays, (14) supernova energy per rest-mass energy of newly formed stars, and (15) rest-mass energy of newly formed stars per supernova.

Download table as:  ASCIITypeset image

2.2. Radiative Cooling

We implemented the Townsend (2009) exact cooling method using the Rosen & Bregman (1995) piecewise power-law form of the cooling function, which extends down to a floor temperature of 300 K. The Rosen & Bregman (1995) cooling function is an approximation to the Dalgarno & McCray (1972) and Raymond et al. (1976) radiative cooling functions and is given by

where T is the gas temperature in K and Λ(T) is in units of erg cm³ s⁻¹. The above cooling function is approximately correct for gas of solar abundance, which is completely ionized at 8000 K. Unlike explicit or implicit solvers, the Townsend integration scheme is exact and does not impose restrictions on the cooling time step. Our tests confirm that the gas temperature evolution computed using this method follows, down to machine precision, the evolution predicted analytically (see Appendix A).

2.3. Cosmic-Ray Decoupling from the Cold ISM

Cosmic rays can be efficiently confined to hot plasmas by scattering off self-excited hydromagnetic waves (Kulsrud & Pearce 1969; Kulsrud 2005). The relative drift speed of the cosmic rays with respect to the plasma v_D occurs at the local Alfvén speed v_A, unless the hydromagnetic waves are damped. Kulsrud & Cesarsky (1971) have shown that ion–neutral damping can significantly boost the relative drift velocity of cosmic rays. Moreover, as the ionization fraction decreases, the Alfvén speed v_A ∝ (density of ionized particles)^−1/2 increases. Both of these effects lead to the decoupling of cosmic rays from the low-temperature ISM.

Note that turbulent damping (Farmer & Goldreich 2004; Lazarian 2016; P. Holguin et al. 2018, in preparation), linear Landau damping (Wiener et al. 2018), and nonlinear Landau damping (Kulsrud 2005) will also increase the relative drift velocity of cosmic rays. All of these damping mechanisms (including the ion–neutral damping explored in this work) dissipate cosmic-ray streaming energy into heat. This will affect the equation of state, which is different for the waves (depending on the Alfvén Mach number), cosmic rays, and gas, an effect we ignore in the present work.

We note that even in the presence of cosmic rays in the cold ISM, the ionization fraction in most of this phase should be low because it is mostly the low-energy (tens of MeV) cosmic rays that are responsible for the ionization of hydrogen. The ionization cross-section of these cosmic rays is very high (Draine 2011), and consequently, they typically do not travel far from the sites of their injection. This allows higher-energy cosmic rays, which carry most of the energy of the cosmic-ray fluid, to propagate in weakly ionized cold ISM where the coupling is relatively weak.

The dynamics of cosmic-ray decoupling from the gas is governed by kinetic theory, yet we model cosmic rays as a fluid in order to perform simulations of tractable duration. The expectation from kinetic theory is that cosmic-ray pressure and energy density tend toward constant values in space when decoupling operates (cf., Everett & Zweibel 2011, who also model decoupling via a large diffusion coefficient). We model the decoupling mechanism in the "extrinsic turbulence" framework (Zweibel 2013) in which cosmic rays scatter off waves generated by turbulence driven by external sources rather than off the waves generated by the streaming instability. In this model, cosmic-ray transport proceeds via diffusion rather than streaming, and cosmic rays exert pressure on the gas but do not heat it (Zweibel 2017).⁵

Cosmic-ray streaming and streaming heating were previously investigated in detail in Ruszkowski et al. (2017). They found that the streaming speed (i.e., the boost factor f > 1) significantly affects wind launching, while whether or not streaming heating is included does not have any noticeable impact on the results (M. Ruszkowski et al. 2017, private communication). In this work, we focus on the influence of decoupling on the transport speed.

In order to emulate the decoupling effect in the fluid model, we adopt a simple treatment in which the parallel diffusion coefficient is amplified to large values in low-temperature regions. This has the same effect as increasing the effective transport speed in the low-temperature gas as discussed above. As such, it will smooth the gradients in cosmic-ray pressure and energy density, matching the expectations from kinetic theory.

The parallel diffusion coefficient can be expressed as κ_Bohm/, where is the ratio of the scattering frequency of cosmic rays on the waves generated by external turbulence to the gyrofrequency c/r_g, and is the Bohm diffusion coefficient, r_g is the gyroradius, and v_p is the particle speed (Schlickeiser 1989; Enßlin 2003). The scattering frequency depends on the properties of the MHD turbulence on scales comparable to r_g. Since we are working in the "extrinsic turbulence" model, the source of the magnetic field perturbations capable of deflecting cosmic rays is most likely compressive waves generated by the Goldreich–Sridhar cascade down to this small scale (Yan & Lazarian 2002). Frequent cosmic-ray scattering on these perturbations reduces field-aligned diffusion. The scattering frequency is proportional to δB²/B², where δB² is the power in the magnetic fluctuations corresponding to the scale equal to r_g. Thus, , and we assume that when significant wave damping is present in weakly ionized regions, the amplitude of the magnetic field perturbations decreases, and the parallel diffusion coefficient is boosted.

To illustrate that increasing the diffusion coefficient has the desired effect of flattening the cosmic-ray energy density distribution, let us consider a one-dimensional and steady-state form of the cosmic-ray energy density equation with spatially constant velocity v₀. For the sake of simplicity, we also ignore supernova heating in this example. In this case, the integration of Equation (5) over z yields

Let . In this example, low-temperature regions correspond to large values of z. The solution of Equation (10), which satisfies the boundary condition e_cr = e_cr,0 at z = 0, is

For z ≫ L, e_cr approaches a constant required to match the behavior expected from kinetic theory. Moreover, if v₀L/κ₀ ≪ 1, the constant value is e_cr ≈ e_cr(0), as expected. On the other hand, in the limit of z ≪ L,

which, as expected, matches the solution of an advection–diffusion equation with a constant diffusion coefficient.

In order to capture the effect of decoupling, we implement in the code the following simple dependence of the diffusion coefficient on the gas temperature,

and set cm² s⁻¹ for all temperatures. Here, κ_∥ and κ_⊥ denote the diffusion coefficient for cosmic-ray transport parallel and perpendicular to the magnetic field, respectively.

The adopted cold gas value of the diffusion coefficient is representative of the values inferred for the Galaxy, but the high-temperature value of the coefficient is lower. However, our volume-weighted diffusion coefficient is expected to lie in between these limits and is comparable to that adopted by Booth et al. (2013), Salem & Bryan (2014), and Ruszkowski et al. (2017). This average value is somewhat smaller than the Galactic value (c.f., Strong & Moskalenko 1998). However, the level of diffusion inferred from observations depends on the assumptions of the models used to quantify it. Specifically, diffusion coefficients derived from the GALPROP propagation model assume a spatially constant diffusion coefficient and/or often assume the absence of winds. Interestingly, Ptuskin et al. (1997; see also Zirakashvili et al. 1996) consider an analytic model that includes both of these effects. They study cosmic-ray-driven winds in which they treat streaming in the diffusion approximation and include decoupling due to ion–neutral damping. They find that the level of diffusion required for consistency with the Galactic data is only ∼10²⁷ cm² s⁻¹ outside the regions close to the disk midplane, and significantly higher close to the disk where decoupling operates. The vertical velocity gradient of the cosmic-ray-accelerated wind in their model is large and similar to the values predicted by simulations (e.g., Salem & Bryan 2014). Furthermore, Jóhannesson et al. (2016); see also Trotta et al. 2011) demonstrate that propagation parameters derived from low-mass isotope data differ significantly from those based on light elements, e.g., B and C, suggesting that these species probe different locations of the ISM where cosmic-ray transport may occur at different rates, though their GALPROP models ignore winds. Individual supernova remnants have also been used to put constraints on the diffusion coefficient. These studies suggested that the locally measured diffusion coefficient can be around ∼10²⁶ cm² s⁻¹ to ∼10²⁷ cm² s⁻¹ when isotropic diffusion is assumed and a somewhat larger cm² s⁻¹ (at 1 GeV) when anisotropic diffusion is assumed (e.g., Nava & Gabici 2013 and references therein).

Our approximate treatment assumes that the gas is fully ionized above 10⁴ K. The ionization level changes dramatically near this temperature threshold, and we note that the results are only weakly sensitive to the exact choice of this threshold. Consequently, the exact form of the cooling function, and specifically its dependence on the gas ionization near this critical temperature, is not critical to our conclusions.

Although we implement decoupling by boosting κ_∥ by a factor of 30 in regions with T < 10⁴ K, we experimented with larger boost factors and found that they did not significantly affect the results. Since the computational time step is inversely proportional to the diffusion coefficient, we decided to use smaller boost factors to accelerate the computations. Tests of the anisotropic diffusion module are presented in Appendix B.

2.4. Star Formation and Feedback

2.5. Simulation Setup

We simulate a slab of ISM, with box dimensions of (2 kpc)² × 40 kpc. This domain shape and size were motivated by the results of Hill et al. (2012), who found that an extended height was crucial in establishing a realistic halo temperature distribution in their simulations. We employ periodic boundary conditions on boundaries perpendicular to the disk and "diode" boundary conditions on boundaries parallel to the disk. Diode boundary conditions permit material to outflow from the simulation box but prevent infall (cf., Sur et al. 2016). Note that we ignore the magnetic field amplification due to rotational shear, and our simulations somewhat enhance gravitational instability since we ignore differential rotation. These caveats are important to keep in mind; however, ignoring differential rotation greatly simplifies the calculation without sacrificing its main objectives.

We use static mesh refinement throughout the duration of the simulations to maximize the resolution near the disk. Although using static mesh refinement possibly underestimates the shock heating of the halo gas and thus the temperature of the wind, it does not affect our main conclusions (i.e., whether winds are launched since it mainly depends on what happens close to the disk, and the relative differences of wind properties among the three transport cases presented below). We achieve a maximum resolution of 31.25 pc in the disk. Beyond kpc, our resolution begins to degrade down to a minimum resolution beyond kpc of 250 pc. One of the factors limiting the simulation time step is the magnetic-field-aligned diffusion. Our maximum resolution is comparable to that achieved in Girichidis et al. (2016), who also included magnetic-field-aligned cosmic-ray diffusion but considered a lower maximum value of the diffusion coefficient.

We initialize a constant temperature of 10⁴ K and a constant magnetic field of strength 1.0 μG along the (horizontal) x-direction throughout the computational volume. The initial density distribution follows a vertical density profile given by

where ρ₀ is the midplane density, z₀ is the scale height of the gas disk, and ρ_crit is the critical density of the universe. The normalization constant ρ₀ is obtained from the disk surface gas density . See Table 1 for the adopted parameter values. For a radially exponential gas distribution of 4.1 × 10¹⁰M_⊙ in baryons in a Milky Way-type galaxy with a scale factor of 3.6 kpc (Booth et al. 2013), the adopted gas surface density Σ₀ corresponds to the gas surface density averaged within a radius of ∼10 kpc from the Galactic center.

The initial setup is in rough hydrostatic equilibrium; that is, the density distribution is such that there would be hydrostatic equilibrium but for the gravity due to the NFW halo. However, gravity due to the NFW halo is small compared to self-gravity near the midplane. Consequently, gas rapidly accretes onto the midplane early in the simulation. We show results for times after 40 Myr, after memory of the initial condition has been forgotten.

We present results from three simulations, which include self-gravity, radiative cooling, magnetic fields, star formation and feedback, and cosmic-ray pressure forces, but differ in their treatment of cosmic-ray transport: in run ADV, cosmic rays are advected by gas motions and diffusive cosmic-ray transport is ignored; in run DIF, cosmic rays are advected by gas motions and additionally, the anisotropic diffusion of cosmic rays along magnetic field lines with cm² s⁻¹ and cm² s⁻¹ is included; run DEC includes cosmic-ray decoupling from low-temperature ISM treated via a temperature-dependent diffusion coefficient as described in Section 2.3.

In Figure 1, we show the projected gas densities in the central at ∼170 Myr for each run. In the ADV run (left panel), cosmic rays injected by supernovae are unable to efficiently drive the dense gas away from the midplane. In this case, cosmic rays remain trapped within the disk, provide additional pressure support, and thus puff up the disk. In contrast, in the DIF run (middle panel), cosmic rays are able to diffuse away from the midplane and drive the more tenuous gas located immediately above the midplane away from the disk, thus producing a wind. In the DEC run (right panel), the effects of stellar feedback are even stronger than in the DIF run—faster cosmic-ray transport out of the ISM reduces cosmic-ray pressure support near the midplane, which leads to enhanced star formation rate, and thus, stronger feedback. This stronger feedback, together with the fact that cosmic rays avoid dense gas in the disk due to their decoupling from the cold ISM phase and preferentially acting on more tenuous gas, leads to the formation of hot, low-density bubbles extending up to ±4 kpc away from the midplane (see the right panel in Figure 1). The timing of the snapshots shown in Figure 1 corresponds to the period following the peak in star formation rate and was chosen specifically to reveal the maximum spatial extent of the hot bubbles.

Zoom In Zoom Out Reset image size

Previous simulations also found that cosmic rays cannot drive winds without diffusive transport (Jubelgas et al. 2008, Uhlig et al. 2012; Simpson et al. 2016). However, the formation of small-scale structure is resolution dependent. We performed an additional ADV simulation at 15.625 pc resolution and still found that no wind was produced. Nevertheless, it is possible that at sufficiently high resolution, the thermal energy from supernovae would carve out channels in the gas density. The channels would allow the cosmic rays to escape rather than puffing up the disk. This hypothesis will be examined by future higher-resolution simulations, which can better resolve high-density clumps within the disk (but this is beyond the scope of the current work).

To better understand the impact of cosmic-ray decoupling on the properties of the ISM, we next consider temperature–density phase plots. Figure 2 shows these phase plots for kpc (bottom row) and kpc (top row). The left column shows results from the DIF run and the right from the DEC run. All panels correspond to 170 Myr. One of the most striking differences between these phase plots is the difference between the phase plots corresponding to kpc. These plots reveal that the low-density bubbles formed in the DEC case (see the right panel in Figure 1) are very hot (10⁶–10⁷ K). This feature is absent in the DIF case. We verified that the hot underdense bubbles are present throughout the kpc region rather than being confined to smaller distances from the midplane. These results also show that the gas in the DEC case is on average hotter than that in the DIF case (i.e., even gas with temperatures below 10⁴ K is less abundant). As mentioned above, this is a consequence of enhanced stellar feedback in the DEC run. At kpc (top panels), the maximum gas density in the DEC case is lower compared to that in the DIF run. This is consistent with the gas surface density distributions shown in Figure 1. This could again be understood as being due to faster cosmic-ray transport in the DEC case. This enhanced transport forces cosmic rays to interact with relatively less dense disk gas, which results in an outflow characterized by lower density.

Zoom In Zoom Out Reset image size

Figure 3 shows the evolution of the vertical wind velocity (top row), cosmic-ray number density (middle row), and magnetic field strength (bottom row) as a function of height above the midplane. All three variables are volume-weighted. The cosmic-ray number density is computed from the simulation output cosmic-ray energy density through n_cr =[(n − 4)/(n − 3)]e_cr/E_min, where n = 4.5 is the slope of the cosmic-ray distribution function in momentum, and E_min = 1 GeV is the minimum cosmic-ray energy.

Zoom In Zoom Out Reset image size

From left to right, the columns show the results for the ADV, DIF, and DEC cases. The profiles are shown from the beginning of star formation at 50 Myr to quiescence at 200 Myr.

We begin the discussion of Figure 3 by considering the vertical gas velocity. In the ADV case, inflows (positive/negative "wind" velocity at negative/positive z) are present for most of the simulation time. Note that regions at large heights above the disk contain very little gas in this case. This has to be contrasted with the DIF and DEC runs that are dominated by outflows. Interestingly, the wind in the DEC case is faster than that in the DIF case (e.g., compare the green-yellow curves near 140 Myr when the wind is about twice as fast in the DEC run) and lasts longer. Both of these effects are the result of stronger stellar feedback and the fact that it is easier to accelerate lower-density gas in the DEC case.

Let us now consider the spatial distribution of cosmic rays and magnetic fields (shown in the middle and bottom rows, respectively). Cosmic rays are most tightly confined to the midplane in the ADV case. This is consistent with our analysis of the gas density projections (see the left panel in Figure 1). On the other hand, in both DIF and DEC runs, the cosmic rays are much more dispersed than in the ADV case; i.e., at large , the cosmic-ray number density is much greater for most of the simulated time. Moreover, the DEC run exhibits a much wider distribution of cosmic rays than that seen in the DIF case.

The trends seen in the evolution of the cosmic-ray distribution are generally reflected in the evolution of the magnetic field profiles. Specifically, the magnetic field distribution is much broader in the simulations that include cosmic-ray transport (DIF and DEC) compared to the ADV case. However, the magnetic field strength is much stronger in DEC than in DIF (except at very late times). In the DEC run, cosmic-ray feedback is substantially more explosive than in the DIF run, launching a strongly magnetized outflow.

While we initialize a unidirectional 1 μG magnetic field, that initial field is quickly erased. In the inflow case (initial stage of the ADV case), the field is simply accreted onto the midplane. As we use diode boundary conditions (inflow through the top and bottom boundaries is not permitted), accretion leads to the reduction of the magnetic field in the regions away from the midplane. Because no wind develops in this case, the gas density above and below the midplane is also very low, as can be seen in the left panel in Figure 1. Thus, we expect the results in the ADV case to be unaffected by our assumption of a spatially constant magnetic field, and an initial field decaying with the distance from the midplane should lead to very similar results. The bulk of the magnetic field amplification occurs only near the midplane as a result of turbulent motions associated with star formation and feedback. In the cases that include transport (DIF and DEC), the field is also amplified in the disk due to star formation and feedback, but following its amplification, it is expelled from the disk. During the outflow, the initial field is swept out of the simulation volume through the outer boundaries. We do not expect our assumption regarding the initial magnetic field to affect our conclusions in these cases either. Thus, the system quickly loses memory of the initial field configuration and strength. It is only in the cases that include cosmic-ray transport that we expect significant magnetization of the gas at large distances from the disk at late times. This magnetization process occurs earlier in the DEC case compared to the DIF case, because the wind speed is larger in the former case.

The above differences between the evolution of the wind velocity, gas density, temperature, magnetic field strength, and cosmic-ray number density will lead to different observational signatures. For example, we expect a stronger, spatially extended soft X-ray emission when the decoupling mechanism operates. The presence of such emission may mitigate the problem reported by Peters et al. (2015), who found that cosmic-ray-driven outflows eject too little hot gas to match the soft X-ray background. Note that if streaming were additionally included, the coupled regions would be collisionlessly heated by cosmic rays, possibly producing even higher temperatures and stronger soft X-ray emission.

Furthermore, elevated cosmic-ray number densities and magnetic field strengths in the halo, combined with advection times shorter than synchrotron cooling times, suggest a more extended radio emission in the DEC case than in the DIF case. However, we note that our simulations reflect the cosmic-ray proton rather than the cosmic-ray electron distribution; cosmic-ray electrons are subject to energy-dependent losses (synchrotron and inverse Compton cooling) which dominate the nonthermal radio emission. We will investigate radio spectra in future work (e.g., via a Lagrangian tracer particle approach to follow the synchrotron aging of electrons comoving with the wind).

Finally, as the decoupling reduces the amount of time cosmic rays spend in the cold ISM phase, we expect that this mechanism would have implications for the γ-ray emission due to hadronic processes. We defer the study of these effects, and the other observational signatures, to a future publication.

In Figure 4, we quantify the properties of the mass flow and compare them to the star formation rates in the ADV (left column), DIF (middle column), and DEC cases (right column). The top row shows the evolution of the star formation rates (solid red lines) and the mass outflow rates computed by integrating mass fluxes through three different pairs of surfaces parallel to the disk. These planes are positioned at ±1 kpc (dashed green lines), ±2.5 kpc (dotted–dashed blue lines), and ±5 kpc (dotted black lines). We find that there is essentially no outflow in the ADV case, and the star formation in this case is very weak. This is consistent with the findings of a number of authors (e.g., Salem & Bryan 2014; Girichidis et al. 2016; Simpson et al. 2016; Ruszkowski et al. 2017). In this case, cosmic rays are confined to the dense disk and the pressure forces they exert are too weak to expel the dense gas from the galaxy. The enhanced pressure support in the disk inhibits the collapse of cold gas clumps and thus significantly suppresses star formation.

Zoom In Zoom Out Reset image size

This picture is significantly altered when transport processes are included in the simulations. We find that in the DIF case, star formation is enhanced compared to the ADV case, and gas is displaced from the midplane. We observe significant gas outflow followed by some inflow (the signature of a fountain flow). When the decoupling physics is included, the star formation rate is enhanced further and an outflow is launched, but there appears to be no inflow. Notice that not only is the star formation peak highest in this case, but the duration of the star formation episode is the longest. In the DIF and DEC cases, there is a delay between the onset of star formation and the outflow.

The star formation rate increases from ADV to DIF to DEC due to decreasing cosmic-ray pressure support in the cold ISM phase caused by faster cosmic-ray escape from the disk. However, including energy-dependent losses of cosmic rays could decrease the cosmic-ray pressure inside dense regions in the ADV case. In such a case, cosmic rays would additionally enhance the pressure in the ambient medium, boosting the star formation rate in this case relative to the DIF or DEC cases since the Bonnor–Ebert mass goes as , where P₀ is the ambient pressure (Ebert 1955; Bonnor 1956). Thus, it is possible that a treatment including energy-dependent losses of cosmic rays would find a boosted star formation rate in the ADV case relative to the DIF or DEC cases. We will investigate the effect of energy-dependent cosmic-ray transport in future work.

In the bottom row, we present the evolution of the integrated mass-loading factor: the ratio of the integrated mass in the wind m_w to that in the stars m_*. In the early stages of disk evolution, the gas cools and quickly settles very close to the disk midplane. Therefore, in measuring the wind mass, we delay the integration of the disk mass until 40 Myr. This allows us to exclude accretion through the set of planes positioned closest to the disk (±1 kpc), which would appear as a negative wind mass. Thus, this approach allows us to better quantify the true amount of gas expelled to large distances from the midplane over time.

In agreement with the findings presented in the first row, we observe that in the ADV case, the integrated mass-loading factor is very low when the mass flux is measured at ±5 kpc from the disk. For smaller heights (±1 kpc), the integrated mass loading is higher, and this simply reflects the fact that the disk puffs up due to the increased pressure caused by the inability of the cosmic rays to leave the disk.

As expected, in the cases that include cosmic-ray transport, the integrated mass-loading factors are much larger compared to the no-transport case. Surprisingly, the level of the integrated mass-loading factors is roughly comparable in both the DIF and DEC cases despite the differences in the transport physics that lead to a number of important differences in the properties of the outflows and their observational signatures. This can be understood by the faster wind speed partially compensating for the much lower gas density in the wind for DEC compared to DIF. The DIF run exhibits integrated mass loadings of ∼0.3, which is comparable to that found in other papers (e.g., Booth et al. 2013), while the DEC run is somewhat smaller, ∼0.2.

We perform simulations of cosmic-ray feedback and its impact on the launching of galactic winds. In our simulations, we find that cosmic-ray transport is essential for driving galactic winds by cosmic rays—in the absence of transport effects, cosmic rays alone fail to drive winds, as found in previous studies (Jubelgas et al. 2008; Uhlig et al. 2012; Simpson et al. 2016).

However, a novel element of our simulations is that they incorporate the effect of the ISM temperature on cosmic-ray propagation. At low temperatures, when the gas ionization fraction is low, ion–neutral friction can damp waves generated by the cosmic-ray streaming instability, and cosmic rays can propagate unimpeded through the ISM rather than scatter off these waves; i.e., cosmic rays are said to be decoupled from the ISM. Furthermore, low gas ionization leads to larger ion Alfvén speed and consequently faster cosmic-ray transport. Both of these effects result in faster cosmic-ray transport in the cold ISM. We model this transport phenomenon by introducing enhanced diffusion in low-temperature regions.

We note here that part of our motivation to treat decoupling via diffusion rather than streaming is to enable comparisons to previous work that also considered diffusion. Our work generalizes previously obtained results by investigating the impact of environment-dependent diffusion. It represents the first attempt to approximate the effects of the decoupling of cosmic rays from the low-temperature plasma on wind launching and the properties of the outflow.

Our simulations focus on a patch of the galactic disk to achieve a resolution higher than would be otherwise possible. These simulations address a well-posed question of how the star formation rates and wind properties are affected by the decoupling of cosmic rays from the low-ionization phase of the ISM. Our specific conclusions can be summarized as follows.

• 1.  
  We observe the formation of low-density and high-temperature bubbles in the simulation that includes cosmic-ray decoupling from the ISM. This raises the possibility that this case may result in enhanced soft X-ray emission from edge-on galaxies undergoing intense stellar feedback. We suggest that the formation of these structures is due to a combination of the increase in (i) the star formation rate and (ii) the effective cosmic-ray transport speed. Our simulations show that faster transport leads to the expulsion of more tenuous gas from the galaxy because cosmic rays avoid cold and dense gas clouds in the disk and preferentially act on lower-density ISM.
• 2.  
  Our simulations corroborate earlier findings that cosmic-ray feedback reduces star formation rates. We emphasize that this does not occur as a result of wind launching. In fact, the impact of cosmic rays on star formation is strongest when transport processes are ignored and no wind is present. Our results are consistent with other studies in that they demonstrate a monotonic trend for the star formation to increase with the average cosmic-ray transport speed. While the star formation rate could be moderated by a number of model parameters, cosmic rays play a very important role in shaping the properties of the outflows and controlling star formation rates.
• 3.  
  Simulations with decoupling exhibit significantly elevated cosmic-ray number densities in the halo at all times compared to the other cases. Combined with the fact that the wind speed is generally faster in this case, and that the advection times are shorter while the synchrotron cooling times may be comparable to those observed in the case without decoupling, we speculate that the wider spatial distribution of cosmic rays may result in broader radio halos when the decoupling physics is taken into consideration.
• 4.  
  Cosmic-ray decoupling reduces pressure support near the galactic midplane due to faster cosmic-ray escape from the cold ISM regions. This faster transport consequently leads to increased star formation rates and further injection of cosmic rays. This may have implications for hadronic losses and associated γ-ray emission.
• 5.  
  Compared to the simulations without temperature-dependent cosmic-ray transport, in the simulation including decoupling, the wind speeds are larger and the wind duration is longer.
• 6.  
  The magnetic fields amplified near the midplane, which subsequently reach large distances away from the midplane, are much stronger in magnitude in the simulation including decoupling than in the case with diffusion. Additionally, since the winds are faster in the decoupling case, the dispersal of the magnetic fields occurs at earlier times in this case.
• 7.  
  While the simulations with decoupling exhibit faster wind speeds compared to the case with diffusion, their winds are characterized by lower gas density. Consequently, wind mass-loading factors—quantified in terms of the ratio of the integrated wind mass to the cumulative mass in newly formed stars—appear to be roughly insensitive to the physics of cosmic-ray transport.
• 8.  
  Simulations with cosmic-ray transport included reveal the presence of both fountain flows and net mass loss in the wind.

We thank the anonymous referee for a very thoughtful report. R.F. gratefully acknowledges Hui Li, Zhijie Qu, Paco Holguin, and Nathan Goldbaum for helpful discussions. M.R. acknowledges support from NASA grant NASA ATP 12-ATP12-0017 and NSF grant AST 1715140. H.Y.K.Y. acknowledges support from NSF grant AST 1713722 and the Einstein Postdoctoral Fellowship by NASA (grant number PF4-150129). E.Z. happily acknowledges support from NSF Grant AST 1616037, the WARF Foundation, and the Vilas Trust. M.R. thanks the University of Maryland College Park, where part of this work was completed, for hospitality. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. Resources supporting this work were provided by the NASA High-End Computing Program through the NASA Advanced Supercomputing Division at Ames Research Center. Data analysis presented in this paper was performed with the publicly available yt visualization software (Turk et al. 2011). We are grateful to the yt community for their support.

To test our implementation of the Townsend (2009) exact integration radiative cooling scheme, we initialized gas of constant density and temperature in a cubical box with periodic boundary conditions. This setup ensured that, despite the fact that the simulation included hydrodynamics, no gas flow developed, and the only quantity with time dependence was the temperature, which decreased due to radiative cooling.

To compare our implementation to an analytic solution, we simplified the problem and considered only one branch of the cooling function,

which corresponds to the uppermost temperature regime of the Rosen & Bregman (1995) cooling curve.⁷ In this case, the thermal energy equation , where e_g is the thermal energy density and n_H is the hydrogen number density, has the following simple solution

where T_f is the temperature after an elapsed time t, T_i is the initial temperature, and k_B is the Boltzmann constant.

In the numerical test, we set K and n_H = 1 cm⁻³, and evolved the simulation for a few cooling times . The agreement between the analytical result and the simulated solution obtained using the Townsend method implemented in the code was excellent (see Figure 5). The simulated temperature evolution agreed with the analytic result to machine precision. Importantly, the temperature did not overshoot the lowest temperature below which the cooling function was set to zero. Other methods may suffer from the overshooting problem in regions where cooling is fast and gas temperatures are low.

Zoom In Zoom Out Reset image size

In order to validate the cosmic-ray decoupling module that we implemented, we performed the following test. We set up a simulation with a temperature-dependent diffusion coefficient κ such that κ(T) = T and T(x) = 1 − x². We used constant gas density and vanishing gas velocity throughout the computational domain (note that all quantities are in code units). Our spatial resolution was 64 zones in the x-direction, and the temporal resolution was 10⁻⁵ time units. The initial cosmic-ray energy density was initialized according to

where e_cr,0 is an arbitrary constant background cosmic-ray energy density, and the last term in Equation (17) is the fourth-order Legendre polynomial. These initial conditions are not isobaric. Since we were interested in testing the implementation of the diffusion module alone, we switched off hydrodynamics in this test, which allowed us to ignore pressure forces.

Using the above initial conditions, we solved the diffusion equation

in a domain that was periodic in the x-direction. The domain consisted of 64 zones and extended from to . These endpoints were chosen to ensure that the initial distribution of cosmic-ray energy density had a vanishing slope at the boundaries of the periodic domain. This allowed us to eliminate jumps in the energy density of cosmic rays at the boundaries. We then compared our solution to the analytic solution of Equation (18) given by

Our simulated solution was in perfect agreement with the analytical solution (see Figure 6), validating the implementation of the decoupling module.

Zoom In Zoom Out Reset image size