Galactic winds driven by cosmic-ray streaming

Galactic winds driven by cosmic-ray streaming

Abstract

Galactic winds are observed in many spiral galaxies with sizes from dwarfs up to the Milky Way, and they sometimes carry a mass in excess of that of newly formed stars by up to a factor of ten. Multiple driving processes of such winds have been proposed, including thermal pressure due to supernova-heating, UV radiation pressure on dust grains, or cosmic ray (CR) pressure. We here study wind formation due to CR physics using a numerical model that accounts for CR acceleration by supernovae, CR thermalization by Coulomb and hadronic interactions, and advective CR transport. In addition, we introduce a novel implementation of CR streaming relative to the rest frame of the gas. Streaming CRs excite Alfvén waves on which they scatter, thereby limiting the CR’s effective bulk velocity. We find that CR streaming drives powerful and sustained winds in galaxies with virial masses . In dwarf galaxies () the winds reach a mass loading factor of , expel of the initial baryonic mass contained inside the halo’s virial radius and suppress the star formation rate by a factor of . In dwarfs, the winds are spherically symmetric while in larger galaxies the outflows transition to bi-conical morphologies that are aligned with the disc’s angular momentum axis. We show that damping of Alfvén waves excited by streaming CRs provides a means of heating the outflows to temperatures that scale with the square of the escape speed, . In larger haloes (), CR streaming is able to drive fountain flows that excite turbulence, providing another means of heating the halo gas. For halo masses , we predict an observable level of H and X-ray emission from the heated halo gas. We conclude that CR-driven winds should be crucial in suppressing and regulating the first epoch of galaxy formation, expelling a large fraction of baryons, and – by extension – aid in shaping the faint end of the galaxy luminosity function. They should then also be responsible for much of the metal enrichment of the intergalactic medium.

keywords:
cosmic rays – galaxies: formation – galaxies: evolution – galaxies: starburst – galaxies: dwarf – intergalactic medium
23

1 Introduction

Large-scale galactic outflows constitute an important process in the evolution of galaxies. They influence the gas content and therefore the star formation rate in galaxies. Moreover, outflows from galaxies also affect the surrounding medium by depositing mass, metals and energy (see Veilleux et al. 2005, for a recent review).

The earliest works on galactic outflows were motivated by the observations of the starburst galaxy M82 (Lynds & Sandage 1963; Burbidge et al. 1964), and the detection of broad optical emission lines in some elliptical galaxies (Osterbrock 1960) that implied winds (Burke 1968). It was suggested that galactic winds cause ellipticals to shed their interstellar medium (ISM; Mathews & Baker 1971; Johnson & Axford 1971). Galactic outflows have since been frequently observed in both local and high redshift galaxies. Multi-wavelength studies for local galaxies have demonstrated the ubiquity of the phenomenon of outflows (e.g., Heesen et al. 2009), while spectroscopic studies of Lyman Break Galaxies at have also revealed the presence of winds (e.g., Adelberger et al. 2003; Shapley et al. 2003). Large mass-loading factors (defined as mass loss rate divided by star formation rate, SFR) with values of a few are estimated for local galaxies and massive star-forming regions at (Martin 1999; Heckman et al. 2000; Sato et al. 2009; Steidel et al. 2010; Coil et al. 2011). The use of classic ionization diagnostic diagrams allows to cleanly separate winds that are driven by either starbursts or active galactic nuclei (AGN; Sharp & Bland-Hawthorn 2010): while AGN photoionization always dominates in the wind filaments, this is not the case in starburst galaxies where shock ionization dominates, suggesting that the latter is an impulsive event. Moreover, after the onset of star formation, there is a substantial delay before a starburst wind develops (Sharp & Bland-Hawthorn 2010). This may suggest that different physical processes, exhibiting different timescales, are responsible for driving winds in starbursts and AGNs, and reinforces the importance of studying different physical mechanisms that could drive winds.

Theoretically, it has proven to be challenging to launch those powerful winds with mass-loading factors of a few from the galactic disc. The first models of outflows invoke the heating of the ISM by repeated supernovae (Larson 1974; Dekel & Silk 1986), in which thermal pressure of the hot gas drives the outflow. Early simulation studies (Mac Low & Ferrara 1999) that placed supernova explosions in simple models for dwarf galaxies predicted comparatively little mass loss. Recently, simulations have begun to approach the necessary resolution to explicitly simulate the formation of a multiphase medium with a dense, cold molecular phase that harbours star formation in full three-dimensional models. The collective action of stellar winds and supernovae drives hot chimneys and superbubbles that start to launch galactic outflows (Ceverino & Klypin 2009), which may be important for low-redshift galaxies (Hopkins et al. 2011). However, these thermal pressure-driven winds suffer from their tendency to destroy the cool gas through evaporation which limits the mass-loading factor and calls for alternative mechanisms (Nath & Silk 2009; Murray et al. 2011). Observed correlations between wind speed and rotation speed of galaxies, as well as the SFR, may be easier to explain with radiation pressure driving on dust grains embedded in outflows (Martin 2005; Murray et al. 2005; Nath & Silk 2009). This appears to be one promising mechanism of explaining outflows in high-redshift massive starbursts (Hopkins et al. 2011).

CRs can also drive a large-scale outflow if the coupling between high energy particles and thermal gas is strong enough (Ipavich 1975; Breitschwerdt et al. 1991; Zirakashvili et al. 1996; Ptuskin et al. 1997). This is expected given the near equipartition of energy in the Milky Way disc between magnetic fields, CRs, and the thermal pressure (Zweibel & Heiles 1997; Beck 2001; Cox 2005). Fast streaming CRs along the magnetic field excite Alfvén waves through the ‘streaming instability’ (Kulsrud & Pearce 1969). Scattering off this wave field limits the CRs’ bulk speed. These waves are then damped, effectively transferring CR energy and momentum to the thermal gas; hence CRs exert a pressure on the thermal gas by means of scattering off Alfvén waves.

Ipavich (1975) studied this process assuming a spherical geometry and that the waves were completely damped away. Later, Breitschwerdt et al. (1991) considered a disc geometry, and calculated the effect of both small and large damping. They found solutions of the outflow equations with realistic ( M yr) mass loss rates from Milky Way-type galaxies. Such winds can explain the small gradient in -ray emission as a function of galacto-centric radius (Breitschwerdt et al. 2002). Similarly, an outflow from the Milky Way driven by CRs explains the observed synchrotron emission as well as the diffuse soft X-ray emission towards the Galactic bulge region much better than a static atmosphere model (Everett et al. 2008, 2010). This match is particularly striking in the hard X-ray band because of the CR Alfvén-wave heating that partially counteracts the adiabatic cooling of the thermal gas and makes a strong case in favour of CR physics being a necessary ingredient in understanding galactic winds. CR-driven outflows, in which the CR fluid provides an additional source of pressure on thermal gas, may eject substantial amounts of gas from spherically symmetric galaxies with a mass outflow rate per unit SFR of order for massive galaxies (Samui et al. 2010). Hence, CR pressure in starburst galaxies would provide a negative feedback to star formation and eventually limit the luminosity (Socrates et al. 2008).

These general ideas and analytical arguments have been supported by detailed numerical simulations employing an implementation of CR physics in the Gadget code (Jubelgas et al. 2008; Wadepuhl & Springel 2011). According to these simulations, CRs that are accelerated by supernova (SN) shocks can exert enough pressure to significantly influence star formation, particularly in low mass galaxies so that the observed faint-end of the satellite luminosity function in the Milky Way can be successfully reproduced. These simulations assumed the CR fluid to be perfectly coupled to the thermal gas and did not allow for CR streaming in the rest frame of the gas. In this paper, we improve upon these earlier studies and investigate the effect of CR streaming in hydrodynamical simulations of galaxy formation. We will show that CR streaming is potentially the dominant energy transport process powering outflows from dwarf galaxies. In Section 2, we introduce the concept of CR streaming and work out a spherically symmetric wind model driven by CR streaming in Section 3. In Section 4, we introduce our simulation setting. In Section 5, we study how CR streaming can drive outflows and assess their specific properties including mass loss rates of galaxy haloes as well as the CR streaming-induced heating of the halo gas. We conclude in Section 6. Supporting material on the numerical implementation, code and resolution tests are shown in the Appendices A and B.

2 Cosmic ray streaming

CR proton populations have several properties which could help galaxies to launch strong winds. First of all, supernovae, being one of the main energy sources of the ISM, are believed to convert a significant fraction () of their released kinetic energy into CRs via diffusive shock acceleration in SN remnants (Kang & Jones 2005; Helder et al. 2009). The CR population forms a very light fluid in the galaxy, with a significant pressure and the tendency to buoyantly escape more rapidly than any other of the ISM components from the galactic disc. Since the CR fluid is coupled via magnetic fields to the thermal ISM components, the CR pressure gradient helps to lift thermal gas out of galactic discs.

Second, the energy loss time-scales of trans- and ultra-relativistic protons due to Coulomb and hadronic interactions with the thermal ISM are long. The CR loss times are longer than typical residence times of these particles in the denser galactic discs, where the target densities are highest (and therefore loss times are shortest), and they are substantially longer than the cooling time of the thermal gas.

Third, nearly all the energy lost by the CR population, be it via particle-particle interactions, via volume work, or via collective plasma effects (as discussed below) is delivered to the thermal plasma. Dissipating CR energy into the thermal plasma increases the total pressure of the composite fluid due to the harder equation of state of the thermal plasma. It can be easily verified that the increase of thermal pressure is larger than the decrease in CR pressure by evaluating the amount of energy transferred per volume element, . The pressure can double in the limiting case of a fluid initially composed of an ultrarelativistic CR population with which dissipated all its energy. More importantly, any energy injected from CRs into the thermal plasma above the galactic disc is protected against the severe radiative losses any thermal gas element would have suffered within the dense and high-pressure environment within the galactic disc.

Thus, CRs could play an important role as a relatively loss-less vehicle for the transfer of SN energy into the wind regions. They can energize the wind with a higher efficiency than the thermal gas, which radiatively looses a large fraction of the energy injected into it within the galactic discs. The mobility of CRs to travel rapidly along already opened magnetic field lines and their ability to open ISM magnetic field lines via the Parker instability (Parker 1966) render them ideal for this function.

The dominant CR transport process along such opened field lines should be streaming as we will explain now. This process appears naturally provided there exists a gradient of the CR number density along a field line. Since the individual CR particles travel with nearly the speed of light on spiral trajectories guided by the magnetic field direction, a gradient of CRs at a location leads immediately to a bulk motion of particles from the dense to the less dense region. This bulk motion is mainly limited by resonant scattering of CR particles on plasma waves. The scattering isotropizes the CRs’ pitch angles, and thereby reduces the CR bulk speed.

The amplitude of plasma waves which are able to resonate with CRs of a certain energy therefore determines the rate at which the streaming CR population is re-isotropized, and thus the limiting bulk speed. The plasma waves are certainly seeded by the unavoidable turbulent cascade of kinetic energy from large to small spatial scales in the violent ISM. They are further amplified by plasma instabilities, which are powered by the free energy of the anisotropic CR population. Thus the CRs, while streaming through the plasma, excite the waves which limit this streaming. The growth of the waves itself is limited by plasma physical damping processes. The energy of the CRs transferred to the waves is finally thermalized and heats the plasma.

The delicate question is at which level does this interplay of CR streaming due to a CR gradient, plasma-wave excitation due to the CR streaming, and wave-damping processes saturate? What is the resulting effective bulk velocity achieved by the CRs?

A unique answer to these questions is difficult. In a low-beta plasma, where the Alfvén speed exceeds the sound speed by a large factor, most works on this topic seem to agree that the effective streaming speed should be of the order of the Alfvén velocity (e.g., Wentzel 1968; Kulsrud & Pearce 1969; Kulsrud & Cesarsky 1971; Kulsrud 2005). In a high-beta plasma, where the Alfvén speed is strongly subsonic, the Alfvén speed can not be the limiting velocity for CRs. This should become clear by the following Gedankenexperiment. If magnetic fields become weaker and weaker in an otherwise unaltered plasma, the Alfvén velocity would approach zero and any CR streaming, which would be limited to this speed, would vanish. However, weaker magnetic fields imply weaker coupling of the CRs to the thermal plasma and therefore intuitively one would expect a larger, or at least constant transport velocity that is independent of the magnetic field strength in this limit.

Consequently, a number of authors have worked on the problem of CR streaming in high-beta plasmas (Holman et al. 1979; Achterberg 1981; Felice & Kulsrud 2001; Yan & Lazarian 2008), and a summary of these works can be found in Enßlin et al. (2011). The consensus among most of these authors seems to be that the limiting velocity is of the order of the sound speed in this situation. It can be several times this speed, or less, depending on details of the pre-existing wave level, thermal and CR energy density, magnetic topology and so forth.

Here, we adopt the pragmatic approach of Enßlin et al. (2011) to parametrize our lack of knowledge on the precise dependence of the magnitude of the CR streaming speed on the plasma parameters and assume that it is proportional to the local sound speed ,

(1)

with a proportionality constant .

An important aspect of this work is the energy transfer of the CRs to the thermal plasma. This energy transfer is able to repower the thermal wind continuously, in particular at large heights above the galactic disc. The CR stream has a significantly higher velocity than the thermal wind, since the CRs’ bulk speed with respect to the thermal gas is . Thus the dominant energy flow of a galactic wind region might be carried by the mobile CR fluid, and not by the much slower thermal gas. The transfer of energy from CRs to the gas is mediated by plasma waves. To the CR population, the process of wave excitation resembles adiabatic work done during the expansion of the CR fluid with respect to the rest frame of magnetic irregularities of interacting plasma waves. Those are generated by streaming CRs and provide efficient scattering partners which are pushed away from the regions with high CR pressure. If the plasma wave motion would be reversed, they would compress back the expanded CRs and thereby return exactly the transferred energy. However, since the CRs are driving these waves, this usually does not happen under the circumstances we are interested in. Only the wave damping itself is a dissipative process, which thermalizes this energy and therefore re-powers the plasma.

For these reasons, the CR energy loss term in the wind equations due to CR streaming formally resembles an adiabatic loss term, which would be reversible if the streaming velocity field converged. Just because the relevant streaming velocity is always diverging away from the region with high CR pressure, the corresponding energy change of the CR population is always negative and the thermal plasma is always heated by CR streaming.

3 Analytical estimates

We begin with analytical estimates of the outflows excited by the heating from CR streaming. For simplicity, we assume a steady state situation and assume that CRs first diffuse out to a height above the disc, comparable to the scale height. Assuming spherical symmetry, we can write the constant gas mass flux per unit solid angle, , where is the gas mass density, is the gas velocity, and is the radial distance. The mass density of CRs, , is fixed by the constant CR flux per unit solid angle in the steady state: , where is the CR speed (for in equation (1)).

We will first calculate the terminal wind speed using the Bernoulli theorem, assuming that there are streamlines along which gas can travel from the disc to a large distance. In the steady state and for spherical symmetry, the energy equation for a compressible fluid is given by

(2)

where F and Q represent the external force and energy flux respectively, is the specific internal energy (thermodynamic energy per unit mass), and is the pressure. Here we have a two-component fluid composed of gas and CRs. Using an adiabatic index of for the gas, and writing , we have the following gas energy equation,

(3)

where is the gravitational potential and is the heating term due to damping of CR-excited Alfvén waves. We neglect radiative cooling for the tenuous gas in the wind in this analytical treatment. Since , we can safely neglect the bulk kinetic energy of CRs as well as their gravitational attraction. Thus, we have (using as an approximation for ultra-high-energy particles),

(4)

Here, the negative sign of indicates the loss of CR energy due to wave excitation. Adding these two equations and integrating, we get the following equation for total energy,

(5)

where is a constant. Equating the values at the base and the end of a streamline, we have

(6)

The sum of the terms inside the parenthesis is times the square of the effective combined sound speed of gas and CRs at infinity, and is negligible for an adiabatic flow. We are therefore left with the following expression,

(7)

where .

We will use the gravitational potential due to a baryonic component and the dark matter halo, given by the Navarro-Frenk-White (NFW) profile (1996). Since , we obtain

(8)

Here, is the total baryonic mass and is the halo mass of a sphere enclosing a mean density that is 200 times the critical density of the universe, is the characteristic radius of the NFW profile, in which the density profile is given by . Assuming the initial gas speed , the gas sound speed at the base, we have an expression for the wind speed at large distance which depends on the CR parameters and the gas sound speed at the base. The condition for an unbound flow or the existence of the wind is given by,

(9)

Multiplying equation (7) with the mass density, we arrive at the energy budget of the problem,

(10)

which essentially states that the kinetic energy of the wind is the difference between the total energy gain and the loss of energy against gravity. In other words, we have,

(11)

Now, using the value of CR mass flux at the base , we have

(12)
(13)

In the next step, we want to assess how the wind speed and mass-loading factor scale with the halo mass and how this compares to the simulations that we will introduce in Section 4. To this end, we have to specify the wind parameters which should reflect realistic ISM values and also match the late-time behaviour of our simulations with an approximate steady state. We use , corresponding to an ISM temperature of 10 K and assume erg cm. We adopt a mass flux (per unit solid angle) of M yr, which is typical in our simulations that show an almost spherical outflow.

The scaling of the base height exhibits two regimes separated by a characteristic circular velocity of . At , the vertical scale height increases roughly as according to two-dimensional fits of edge-on exponential discs in the -band (Dalcanton et al. 2004).4 At , the vertical scale height appears to be independent of .

This behaviour can be understood theoretically by the following line of arguments. In shallow gravitational potentials, as found in dwarfs, the vertical scale height should scale with the system size. Hence, we assume the base height to scale linearly with the disc size (which is in turn a function of halo mass, according to the model by Mo et al. 1998). When the gravitational potential of the halo is deep enough, the vertical scale height is set by the effective equation of state of the ISM. For an effective equation of state of , the solution to the hydrostatic equation yields a vertical scale height that is independent of surface mass density (see, e.g., Springel 2000). It turns out, that the sub-resolution model employed by our simulations has such a stiff equation of state just above the star formation threshold (see Fig. 2 of Springel & Hernquist 2003). Hence, the vertical scale height should be set by radiative cooling (which breaks the scale-invariance of gravity seen in smaller systems) in haloes heavy enough to support a thin disc, i.e., for corresponding to . This justifies a constant vertical scale height for these large haloes. For our main models, we adopt the following simple prescription for the base height, for . We note that these considerations are reflected by our simulations where we find identical trends of with halo mass.

Figure 1: Wind speed as a function of halo mass (equation (12)) for four models where we vary the input parameters, namely the halo concentration , the constant of proportionality connecting the base height and disc scale , and the maximum base height . The blue dotted curve shows the maximum rotation speed as a function of halo mass. The triangles are the results of the simulation: km s for M, while the wind in the last case (grey triangle) was only sustained for . Filled circles show the circular velocity measured in the simulations at the virial radius.

Figure 2: Mass loading factor (equations (15) and (16)) as a function of halo mass for two values of the gas fractions, and 0.01 (solid and dashed). The two different models both have a halo concentration , as in Fig. 1. The analytic models compare well with our simulated values for the “peak mass loading” (, filled circles).

In Fig. 1, we show the wind speed at a large distance as a function of halo mass. This is compared with the rotation speed of the disc implied by the halo mass, adopting the model of Mo et al. (1998). We note that in this model, the disc mass (which we equate with ) is a factor of the halo mass. In the low-mass regime, the wind speed increases with halo mass (in accordance with the simulations) up to wind velocities of , depending on the adopted parameters. At halo masses around , the CR-driven wind is not any more powerful enough to overcome the increasing potential difference of the halo and ceases. It is informative to contrast this behaviour of our main models (solid lines) with two limiting cases where we either adopt a constant (dash-dotted) or do not impose a cutoff value to the scaling (dashed).

The behaviour of our main models is in reassuring agreement with our simulations where the wind speed increases with halo mass up to a critical halo mass of . At this mass scale, we see an onset of a wind that stalls after , coinciding with the time of the formation of the thin disc. Here the increasing stellar mass deepens the central potential to the point where the CR-driving is not any more powerful enough to support a wind; instead we end up with a fountain flow. The normalization of the wind speed and the transition mass varies with changes in the base height, , and the halo concentration, i.e., the effective potential difference from the starting point of the wind to infinity. Interestingly, in the model with the high concentration in Fig. 1, the wind speed starts to flatten towards lower halo masses. This is because for these small haloes the (constant) sound speed starts to dominate over the terms accounting for the loss of energy against gravity and the energy gain due to CR pressure in equation (12).

To estimate the mass loading of the wind, we first estimate the disc surface density , where the disc mass is . The gas surface density can be written as , where is the gas mass fraction. To determine the SFR per unit area, , we use the Schmidt-Kennicutt relation (Kennicutt 1998),

(14)

Hence the SFR of the disc is obtained as

(15)

We can estimate the value of the mass flux per unit solid angle, , for escaping winds, by requiring that in the wind equation (equation (12)) and by inverting it to get (for a spherical outflow):

(16)

This is reasonable since the escape speed for a NFW halo is (Sharma & Nath 2012). We can then divide this by the corresponding SFR of equation (15) to obtain the mass loading factor .

In Fig. 2, we compare the estimated mass loading factor with the corresponding values obtained in the simulation. We find that our simple theoretical estimate nicely describes the mass loading observed in simulations, which scales as . Hence we conclude that the mass loading in CR-driven winds scales with the halo mass in exactly the way that is needed to explain the low-mass end of luminosity function according to semi-analytical models of galaxy formation (Bower et al. 2011) and phenomenological simulation models (Puchwein & Springel 2012). CR driven winds are thus also very promising candidates to help explaining the “missing satellites problem” in the Milky Way.

4 Numerical simulations

4.1 Simulated physics

We now turn to numerical simulations of galaxy formation that include CR physics. Our simulations were carried out with an updated and extended version of the distributed-memory parallel TreeSPH code Gadget-2 (Springel 2005). Gravitational forces were computed using a tree algorithm. Hydrodynamic forces are computed with a variant of the smoothed particle hydrodynamics (SPH) algorithm that conserves energy and entropy where appropriate, i.e., outside of shocked regions (Springel & Hernquist 2002).

Radiative cooling was computed assuming an optically thin gas of primordial composition (mass-fraction of for hydrogen and for helium) in collisional ionization equilibrium, following Katz et al. (1996). We account for star formation according to a sub-resolution model by Springel & Hernquist (2003) which assumes that star forming regions establishes a self-regulated regime, described by an effective equation of state. The self-regulation arises due to the interplay of cold, dense clouds that continuously form stars and supernova feedback that evaporates these clouds, creating a hot ambient medium. The rate of star formation is adjusted such as to reproduce the observed Schmidt-Kennicut law. As far as the important parameters of this model are concerned, we choose Gyr, thereby setting the characteristic time-scale of star formation, and a density threshold of star formation of (in units of hydrogen atoms per unit volume). This choice is motivated by Dubois & Teyssier (2008) who studied wind formation using similar initial conditions and a similar description of star formation. Their work differs from ours in that they considered superbubbles as the mechanism that drives outflows and they observed the strongest winds for the star formation time-scale stated above. Since CRs as well as superbubbles are powered by supernova explosions, adopting their time-scale should also give the most favourable results in our case.

We model CR physics with a formalism that follows the most important CR injection and loss processes self-consistently while accounting for the CR pressure in the equations of motion (Pfrommer et al. 2006; Enßlin et al. 2007; Jubelgas et al. 2008). In our methodology, the non-thermal CR population of each gaseous fluid element is approximated by a simple power law spectrum in particle momentum, characterized by an amplitude, a low-momentum cut-off, and a fixed slope . Adiabatic CR transport processes such as compression and rarefaction, and a number of physical source and sink terms, which modify the CR pressure of each particle, are included. As the most important source of galactic CRs, we consider CR acceleration by supernova remnants with a spectral index for the freshly accelerated CR population of and an acceleration efficiency, i.e., the fraction of supernova feedback energy directed into the CR population, ranging in between and 0.3 (where the latter value is our standard choice). As sinks of CRs, we consider thermalization by Coulomb interactions and catastrophic losses by hadronic interactions.

In addition to the advective CR transport, we have implemented CR streaming relative to the rest frame of the gas. The details of the numerical implementation in Gadget-2 are described in Appendix A. We set the streaming speed equal to the local sound speed (i.e., in equation (1)).

4.2 Initial conditions and settings

Resolution NGB
(DM/gas)
64
64
64
Table 1: Simulation parameters adopted in this study. Here, and denote the virial mass and radius, and denote the particle masses used to sample dark matter and baryonic gas, and are the gravitational softening lengths for the two components, regulates the size of the streaming time step (see equation (42)) and NGB is the number of SPH smoothing neighbours used in this study.

In our simulations, we consider disc galaxy formation in isolated dark matter haloes, according to the standard picture (Fall & Efstathiou 1980). Initially, the dark matter potential well is filled with a hydrostatic, slowly rotating, baryonic gas in equilibrium with the dark matter. Right after the start of the simulation, radiative cooling sets in, removing thermal pressure support of the gas atmosphere which will then collapse, subject to gravity, to form a rotationally supported galactic disc in the centre of our halo. The gas in the disc becomes colder and denser until it starts to form stars at a rate given by our sub-resolution model. This picture of galaxy formation in isolation is clearly oversimplified, since galaxy formation is a highly hierarchical process, shaped by frequent mergers of different dark matter haloes. However, isolated haloes enable us to study the impact of CRs in a comparatively clean and well-controlled environment, so that our simulations should be able to reveal the significance of CR streaming for driving outflows and to make some qualitative predictions on the size of the effects.

Both dark matter and baryonic gas are set up in such a way that they follow the same NFW density profile. Furthermore, the two components are initially in equilibrium, i.e., the gas atmosphere will be stable in the absence of radiative cooling when evolved over a Hubble time. In addition, the halo carries an initial amount of angular momentum described by the spin parameter , and with a radial distribution in agreement with results from full cosmological simulations Bullock et al. (2001). We adopt a concentration parameter of the NFW profile of in all our runs. The initial mass fraction of baryons is kept fixed in all cases at its universal value, , to ease comparison to the results of Jubelgas et al. (2008). Our simulations represent the baryons with SPH-particles and the dark matter with twice as many particles, which produces numerically converged results (see Appendix B.2). Further, we use SPH smoothing neighbours in all our simulations. An overview of the most important numerical parameters adopted in our simulations is given in Table 1.

5 Results

5.1 Launching of galactic winds

CR streaming and advection CR advection-only
Figure 3: Maps of the baryonic gas density and the gas velocity field (indicated by arrows) for a halo of total mass (top panels) and (bottom panels) at . The maps show an edge-on slice through the midplane of the galaxy that forms in the centre. The left panels correspond to the simulations with CR streaming and advection, while the right panels show the simulations with purely advective CR transport. While slow (thermal) gas motions dominate the velocity field in the advection-only run, the simulation employing CR streaming and advection shows a powerful outflow.

Since outflows cannot occur in simulations without CR feedback (or radiation pressure), at least for the specific star formation model that we use, we will in the following compare simulations with only advective CR transport to simulations with both, advective and streaming transport of CRs. Since the advection-only case was shown not to drive outflows (Jubelgas et al. 2008), this numerical setup enables us to assess the question whether CR streaming is able to launch galactic winds. In Figure 3, we show maps of the baryonic gas density in haloes of total mass (top panels) and (bottom panels), seen in an edge-on slice through the midplane of the galaxy that forms in the centre. Those maps are supplemented by velocity vectors of the gas and are shown at after the start of the simulation when the initial burst of star formation has ended and the galaxies have settled approximately into steady state. For reference, the virial radii of the two haloes are and , respectively.

When comparing the simulations with CR streaming and advection (left) to the advection-only case (right), it becomes clear that the additional mobility of the streaming CRs enables the driving of a powerful outflow. This manifests itself in a shallower density profile in the CR-streaming case, suggestive of a substantial mass transport beyond the virial radius. The outward pointing gas velocities with maximal values of demonstrate the presence of a wind driven by streaming CRs that extends well beyond , allowing the gas to escape from the gravitational influence of the halo. At the largest scales after  Gyrs, we obtain converged values for the wind velocities of and for our two haloes ( and ), respectively. (Note however that owing to our simplified initial conditions, our simulations cannot be taken as a realistic picture for radii larger than a few where the wind will encounter the anisotropic mass distribution of the cosmic web.)

Comparing the CR-streaming simulations of the different halo masses, we observe a quasi-spherical outflow in the dwarf halo with . In contrast, there is a bi-conical, hour-glass shaped wind in the halo within that starts to become more spherical on larger scales. This different wind morphology can be traced back to the density distribution of the star-forming ISM which determines the initial conditions at the base of the outflow. While the shallow potential of the dwarf halo is not able to confine the ISM into a thin disc and instead has an ISM with a quasi-spherical morphology, the larger halo shows a well-confined galactic disc. Such a configuration has a vertical stratification of CR pressure which determines the initial direction of the CR streaming and focuses the forming wind into azimuthal symmetry.

We now detail the physics of the launching mechanism of the wind. Initially, the CR pressure distribution peaks at the star forming regions in the galactic disc which accelerate them in SN remnant shocks. This causes a steep CR pressure gradient that implies the onset of CR streaming out of the galaxy. A fluid element above the disc which is in approximate hydrostatic equilibrium with the surrounding fluid elements receives the (instreaming) CR flux and becomes overpressured. CR loss processes (wave heating as a result of CR streaming instabilities or particle-particle interactions) transform a part of the CR pressure to the thermal plasma. Since the thermal plasma has a harder equation of state, it gains more in terms of pressure than the CRs have lost by this transfer (see Section 2). This increases the total pressure of this fluid element furthermore. This excess pressure causes a weak shock that accelerates the gas in the stratified atmosphere. Momentum conservation ensures that this weak shock moves in the same (outwards) direction as the original CR-streaming flow. The rarefaction wave, moving in the opposite direction, dilutes the gas in the tail and causes a buoyant rearrangement of the flow. We see that the streaming CRs are crucial in this picture for launching the outflow.

Figure 4: Radial velocity of the gas in a cone of opening angle , aligned with the angular momentum axis of the disc for our three haloes of masses , , and at . In the two lower-mass haloes (as well as for in the halo), we see an accelerating wind over the entire radial range which is due to a continuous CR momentum and energy deposition during the ascent of the wind in the gravitational potential. For in the halo, the wind stalls and falls back onto the disc (shown by the negative velocity values) which is the characteristics of a fountain flow.

A central question for winds is the evolution of their outflow velocity along the streamlines and whether they are energy or momentum driven. In the first case, it is known to be difficult to obtain the observed large mass loading factors. Additionally, if the (uncertain) clumping factor in the wind is larger than unity, the wind would radiate away more energy per unit time than in the smooth case which could potentially stall it. In Fig. 4 we show the radial velocity of the gas in a cone of opening angle , aligned with the angular momentum axis of the disc. In the two lower-mass haloes, we see an accelerating wind over the entire radial range. This is due to the particular property of CR streaming-driven winds, which represents neither a classic energy-driven nor a momentum-driven wind in a sense that they are not preloaded with those conserved quantities. Instead, they are continuously re-loaded with energy and momentum during their ascent in the gravitational potential through the Alfvén-wave heating term in the energy equation and the term in the momentum equation.

While there was initially a wind launched in the halo that propagates on scales (at ), it apparently stalled for and started to fall back as a fountain flow onto the galactic disc (see Fig. 4). We observe the same behaviour for later times in the halo. The reason for this is the increasing stellar mass and the formation of a thin galactic disc which causes the wind to be launched deeper in the gravitational potential (implying a small ) and to loose more binding energy during its ascent in the potential.5 Another effect that can partly affect the sustainability of winds in larger haloes are the increased CR losses due to Alfvén-wave heating in the dense ISM which is characterized by steep CR pressure gradients. For our current description of the ISM with an effective equation of state, we may overestimate the CR losses as CRs are kept for a longer time in the dense phase of the ISM. In reality, they may be able to escape more quickly into the warm and more dilute phase of the ISM that is characterized by a smoother CR pressure gradient so that less CR energy is lost at the initial phases of the wind launching. Without improving our modelling of the ISM or the CR transport within the ISM the resulting wind velocities could thus be underestimated. We will return to this point in Section 5.2 and estimate its effect on the cumulative mass loss in a halo.

halo halo
Figure 5: CR-to-thermal pressure ratio, , in an edge-on slice through the galactic disc. Shown is the simulation with CR streaming and advection in haloes of (left) and (right) at time . The CR-to-thermal pressure ratio is less than in the vicinity of the centre, because of loss processes that effectively transfer CR energy into the thermal reservoir, but becomes dominant at larger heights due to the softer adiabatic index of CRs.

To highlight the launching mechanism of the wind, we study the CR-to-thermal pressure ratio, , in an edge-on slice through the galactic disc in Fig. 5. Generally, we find an increasing CR pressure fraction, , at larger radii. As the wind propagates, the composite gas of CRs and thermal plasma experiences adiabatic expansion so that the pressure drops less quickly than the thermal pressure due to the composite’s softer equation of state. Thus, the CR-to-thermal pressure ratio rises to at large radii. Apparently, this effect wins over CR energy losses due to CR Alfvén-wave heating during the ascent of the wind in the halo potential.

The map of the halo shows a homogeneous morphology for with values around . At larger radii, the morphology of the map becomes patchier. It is interesting that CR streaming is unable to smooth this inhomogeneous CR pressure distribution. This is because of the large wind velocities which reach values of , whereas the sound speed (equal to the streaming speed in our model) is only around or less. Thus, advection dominates the transport on these scales which is not expected to result in a smooth distribution of CR pressure and apparently does not significantly mix. The map of the halo shows the hour-glass morphology and supports the picture of a bi-conical outflow driven by CR streaming. Clearly visible are the torus-shaped regions in blue on scales below the virial radius which indicate a low relative CR pressure. This is due to a converging vortex flow towards the midplane and a necessary consequence of the bipolar outflow. Here, the CR pressure component is disfavoured in comparison to the thermal component upon adiabatic compression.

5.2 Mass loss due to galactic winds

Figure 6: Time-evolution of the baryonic mass (in form of gas and stars) enclosed within the virial radius, , normalized to the initial mass at the beginning of the simulations. Left: we compare for four different scenarios (as indicated in the legends) in a halo of total mass . Right: we show for three different halo masses, , , and , always adopting the CR streaming and advection scenario with .
halo halo
Figure 7: Mass loss rates and SFRs as a function of time for the CR streaming and advection case, computed within the virial radius in haloes of total mass (left) and (right). The mass loss rate was calculated for two different acceleration efficiencies, and 0.3.

After studying the launching mechanism of galactic winds through CR streaming, we will now quantify the associated mass loss and the mass loading of the wind. In Fig. 6, we show the fractional mass as a function of time for a reference simulation without CR feedback, for a simulation with advective CR transport-only, and for two runs that additionally include CR streaming and differ only in the adopted CR acceleration efficiencies of and 0.3. We define fractional mass as the baryonic mass in form of stars and gas contained within a certain radius at a given time, normalized to the total gas mass within this radius at the beginning of the simulation.

The reference run without CR feedback exhibits no mass loss and so does the advection-only run, in accordance with the findings in Section 5.1. The fact that there is no wind in the former is connected to the implementation of stellar feedback that we use, which does not allow for outflows without additional adjustments. In contrast, we detected a strong wind in our simulation with CR streaming which results in a considerable mass loss (see Fig. 6). Depending on the assumed CR acceleration efficiency, a fraction of % of the original gas mass contained within the virial radius of the halo at the start of the simulation is expelled until . The non-zero slope of the mass loss history indicates a moderate continuation of the mass loss if we increased the run time of our simulation.

The mass loss history shows a strong dependence on halo mass with a fractional mass loss of only % from the virial radius in the case of the and haloes (see right panel of Fig. 6). Thus, these outflows are weaker in comparison to that in the dwarf halo of , and do not result in a severe mass loss. Since the mass loss history has saturated in the two high-mass haloes, further mass loss is not expected for increasing simulation time. Apparently, the gravitational attraction is too strong for the CRs to launch a strong wind in these cases, but certainly drives powerful fountain flows.

However, the mass loss history in Fig. 6 demonstrates that the outflow is not able to develop until which is due to the ram pressure from inflowing gas that the wind needs to overcome first. Ram pressure can therefore significantly delay wind formation and is probably more important in this respect than the shallow gravitational potential of this halo. This result is in agreement with previous studies of outflows (Fujita et al. 2004; Dubois & Teyssier 2008) and may be partly due to our simplified initial conditions. If most of the gas mass is accreted onto a galaxy along filaments that cover a small solid angle, the ram pressure in directions other than the filaments would be substantially reduced, possibly allowing for an earlier onset of a wind after the first burst of star formation.

It is interesting to compare the mass loss rate in the wind with the SFR. These are expected to be roughly proportional to each other, as simple theoretical predictions and observations suggest. Thus, in Fig. 7, we show the time evolution of the SFR and the mass loss rate within the virial radius for two CR streaming simulations that differ only in the acceleration efficiency of CRs. The influence of ram pressure is apparent, which delays a significant escape of gas until .6 Once the CR-driven wind overcomes the ram pressure, the mass loss rate increases sharply.

halo halo
Figure 8: Top panels: SFR as a function of time in a halo of total mass (left) and (right). We compare five different simulations: a reference simulation without CRs (black line) and four runs with CR feedback, differing only in the acceleration efficiency and the transport scheme, i.e., advective transport-only (orange line: , green line: ) and CR streaming and advection (blue line: , red line: ), respectively. CR feedback leads to a strong suppression of the SFR in comparison to the reference simulation. Bottom panels: Time-evolution of the total mass lost in the wind out of the virial radius (dashed red line), and of the total mass in stars (solid red line) for the CR streaming run. These are compared to the total mass in stars in the advection-only run (green line) and in the reference simulation without CRs (black line). Both runs with CRs adopt . In the halo, the amount of mass removed from the halo in the CR streaming run is considerably larger than the stellar mass that forms in this simulation, in agreement with the large mass loading of the outflow. Moreover, the total mass in stars is significantly less than in the reference simulation.

In case of the high acceleration-efficiency run for the dwarf haloes of , the mass loss rate is about times the SFR and about times the SFR for the lower acceleration efficiency (we always consider the “peak mass loading”, i.e., we divide the respective maxima of the mass loss and star formation histories). For dwarfs, the mass loading of these winds is therefore considerable. In the end of our high acceleration-efficiency simulation, we find that the wind has expelled a gas mass fraction of while of the initial gas mass ended up in stars (which was considerably smaller in comparison to the stellar mass fraction of for our run without CR feedback, see Fig. 8). In contrast, the wind in the halo only has a small mass loading of about for a high acceleration efficiency. While only of the initial gas mass was expelled by the wind in this halo, the stellar mass fraction dropped from down to due to feedback by CR streaming as well as by advected CRs (Fig. 8). Apart from that, in the low acceleration-efficiency CR streaming run there is no mass loss out of the virial radius.

The energy transfer from the CRs to the thermal plasma following the excitation and damping of hydromagnetic waves within the dense ISM, that has a comparatively short cooling time, causes a large fraction of this energy to be lost to cooling radiation. If CRs are injected into the warm phase of the ISM from where they would drive a wind, then our simplified treatment of the multiphase structure of the ISM, which does not explicitly treat the warm and hot gas in superbubbles, underestimates the effect of CR streaming and hence the associated mass loss rates. To estimate the importance of this effect, we switch off the conversion of CR to thermal pressure in our CR streaming implementation and find a fraction of % of the original gas mass contained within the virial radius of the halo is already expelled by . For the reminder of this work, we always include CR-wave heating but bear in mind that the results may be too conservative as a result of too efficient CR energy losses to radiation in the dense, cool phase of the ISM.

5.3 Suppression of star formation rates

To study how CR feedback influences the star formation process, we plot the star formation history in Fig. 8. We vary the simulated physics and compare a model without CR feedback, one with advective CR transport only, and one that additionally includes CR streaming. For each of our CR models, we employ two CR acceleration efficiencies of and 0.3, respectively. Obviously, if feedback by CRs is taken into account, the SFR is suppressed in comparison to the reference simulation without CRs. In case of the CR streaming model with a high acceleration efficiency, the SFR is suppressed by about () at the peak of the star formation history for our halo with (). In our model that only accounts for advective CR transport, this suppression is even increased and amounts to more than () for our halo with (); in accordance with the findings of Jubelgas et al. (2008). As laid out there in detail, the suppression of the SFR by CR feedback is less effective in higher mass haloes which attain higher central gas densities due to their deeper potential wells. This is because adiabatically compressing a composite of CRs and thermal gas disfavours the CR pressure relative to the thermal pressure due to the softer equation of state of CRs. Hence, for the higher central densities present in larger haloes the relative importance of the CR pressure as well as their modulating effect on the SFR decreases.

There are two reasons why the inclusion of CR streaming suppresses the SFR in comparison to the reference model without CRs. First, CR pressure that quickly builds up as stars form and eventually explode will hinder the gas from collapsing to densities much larger than the threshold of star formation. This is due to the inefficient cooling processes of CRs as opposed to the thermal gas that cools comparatively quickly. Consequently, the SFR, that scales with the gas surface density as , will be reduced. As soon as CRs are removed by cooling processes, collapse can start again and the SFR increases. The overall effect is a suppression of star formation as well as an oscillatory behaviour of the SFR, that can be appreciated in Fig. 8. If CR transport via streaming is included, these oscillations basically vanish because freshly accelerated CRs continuously stream away from the star forming regions, allowing the gas to reach higher densities and consequently enhancing the SFR in comparison to the advection-only run.

Second, the wind launched by CR streaming prevents a large fraction of the infalling gas from accreting onto the galactic disc. Hence, the amount of gas accumulated in the galactic disc is reduced. This lowers the available amount of gas that may be turned into stars as well as the surface gravity. The resulting smaller central gas densities imply a lower SFR.

The runs adopting a low CR acceleration efficiency of for the CRs show less suppression of the SFR in comparison to their corresponding high-efficiency counterparts. This is due to the lower contribution of CR pressure to the total pressure in the star forming regions, weakening the influence of CRs on the process of star formation, as outlined above. Additionally, the mass loading of the winds is also reduced in this case. Nevertheless, the suppression is still significant and amounts to () at the maximum of the star formation history in the CR streaming model (advection-only model) for the halo.

The bottom panels of Fig. 8 show the time-evolution of the total amount of mass lost from the halo’s virial radius for our CR streaming simulation with as well as the total amount of mass formed in stars at a given time in this simulation. We contrast this to the cumulative stellar mass for a model without CRs and our model with advective CR transport-only. This shows again that the mass loading of the outflow in the CR streaming case is large, with the total mass loss being about times as large as the mass in stars at the end of the simulation. Furthermore, the final stellar mass at the end of the CR feedback runs is reduced by relative to the final stellar mass in the reference simulation (considering the halo). This is in good agreement with the suppression factor of the SFR at the peak of star formation history. As expected, the reduction is strongest in the model that exhibits only advective CR transport which results in the least amount of stars formed at the end of the simulation.

In contrast, the baryonic mass lost in the galactic wind of the larger halo of is only % of the total mass in stars at , reflecting the small mass loading of the wind (see Section 5.2). In this larger halo, the final stellar mass at the end of both CR feedback runs is only reduced by relative to the final stellar mass in the reference simulation. The large reduction of the SFR at the peak of in the advection-only run is therefore not seen in the final stellar mass. This is because of the different star formation history which shows more continuous rather than a bursty behaviour of the CR streaming models for this halo.

5.4 How CR streaming heats the halo gas

Figure 9: Heating mechanism of the halo gas in a halo of total mass for the CR streaming model. Top panels: Time evolution of the gas temperature and the gas velocity field in an edge-on slice through the galactic disc. Results are shown at two different epochs in the simulation, at (left panels) and (right panels). The maximum velocities encountered in these runs are 26 (left panels) and 70 (right panels). Collimated cavities of hot gas form below and above the disc and propagate further up in the course of time. Middle panels: Corresponding plots for the time-evolution of the ratio of the wave heating rate due to the damping of self-excited waves to the gas cooling rate, , demonstrating that wave heating is responsible for heating the cavities. Bottom panels: Gas density in a cutting plane parallel to the disc at a height of (to assess the density distribution in the hot cavities). At , an under-dense channel forms in the centre where the cooling rate is consequently lowered so that wave heating can dominate over the reduced cooling rate.
halo halo halo
Figure 10: Temperature distribution at the time of maximum CR Alfvén-wave heating in an edge-on slice through the galactic disc. We compare three different haloes of mass , , and (left to right) in our CR streaming model with an acceleration efficiency of . Note that the resulting halo temperatures roughly scale as . The temperature structure resembles that of the wind, which implies that with increasing halo size, the morphology of the hot patches becomes more conical. The broadening of the hot regions in our halo is associated with the inability of CR streaming to drive a sustained wind that escapes from such a halo. Hence, the kinetic energy of the fountain flow drives turbulence which dissipates energy and thereby heats larger regions of the halo gas.

In this section, we assess the role of the wave heating that is intimately connected to the physics of CR streaming. To this end, we only concentrate on our CR streaming model that successfully launches winds. In Fig. 9 we plot the gas temperature in an edge-on slice at two different epochs, and , about after the outflow started. At we can see collimated flows through chimneys of hot gas, extending up to above and below the disc. In these chimneys, the gas temperatures reach values greater than  K. In the centre of the maps, we can see that the disc is also heated up to  K owing to supernova feedback.7 These hot regions are surrounded by a torus-shaped region, filled with comparatively cool gas at  K.

The only conceivable process which could provide a source of heating several kiloparsecs above the disc, where the gas density is low, must be the damping of CR-excited waves. In Sec. 2 we have seen that the hydromagnetic waves are created by streaming CRs in the background plasma. Those waves are damped on short time-scales, thereby transferring energy to the plasma. This heating process is in operation as long as CRs are streaming through the rest frame of the gas, i.e., as long as there is a gradient in the CR pressure. Thus, unlike Coulomb or hadronic losses, which scale with gas density, wave heating is also important in low-density regions. In order to demonstrate that this process is responsible for the hot gas chimneys, we show an edge-on slice through the galaxy with the ratio of the radiative cooling rate of the gas to the wave-heating rate of equation (49) due to the damping of self-excited waves, (middle panels of Fig. 9).

At , wave heating is nowhere able to overcome the radiative cooling of the gas. This is in particular true for the galactic disc, where the density is so large that the cooling dominates the heating by more than two orders of magnitudes. At , however, we can see the structure of the hot gas chimneys again, traced by the evidence that the wave heating dominates the cooling there by about a factor of ten, showing that heating via wave damping is indeed the process that creates the hot cavities in our simulations.

To better illustrate the reasons for the dominance of wave heating over gas cooling above and below the disc, we take a slice parallel to the disc that cuts the outflowing gas stream at a height of above the disc, which is near the base of the upper chimney. We plot the gas density in this slice in the bottom panels of Fig. 9. At , when the cavities have not yet formed, we see that the outflow covers a circular area in the plane and features a denser core in the centre, with decreasing density towards the outskirts. This indicates that the outflow is very collimated. Later, at , an under-dense hole has appeared in the very centre of the core. Thus, the gas cooling rate will drop there, so that the CR-wave heating can heat up the gas.

What is the reason for this under-dense channel in the centre of the outflowing material? Interestingly, the hot chimneys first occur about after the outflow started. The likely reason for this is that the CR-driven wind becomes more and more collimated as time progresses, owing to the disc that forms in the centre of the simulation box. When the wind first occurs, the disc formation is not yet finished, so that it still has an approximately spherical shape. At these early times, CRs can stream at a large angle with respect to the -axis without encountering too much disc material, resulting in an outflow with a wide opening angle. Later, when the disc has formed, the inertia of the dense and cool star forming gas prohibits a long pathway of the CRs through the disc and forces an outflow with a smaller opening angle. When the outflow is collimated enough, it will then quickly dig a diluted channel in the old ejecta above the disc with a correspondingly lower cooling rate.

This argument is also in agreement with the fact that we do not observe chimneys of hot gas for our lower mass halo. The weaker gravity of those dwarf haloes does not support a thin disc forming at the centre, but a rather spherical density distribution which is unable to collimate the wind. This is highlighted in Fig. 10 which shows temperature maps for haloes with , , and . These show a clear dependence of the wave heating on the halo mass. The resulting halo temperatures roughly scale as with the largest halo in our simulations reaching temperatures in excess of  K so that the outflow is expected to emit thermal bremsstrahlung emission. Besides, the larger SFR obtained for the higher mass halo and correspondingly the higher CR energy densities could contribute to the preferred creation of the chimneys in the higher mass halo. We also see in our largest halo that the heated regions are less collimated in comparison to the intermediate-mass haloes: the wind starts deeper in the gravitational potential of this halo (see Section 3) and looses a good fraction of its kinetic energy in climbing up the greater potential difference so that it is unable to counteract the ram pressure of the infalling gas (that reaches also bigger infall velocities in assuming the larger gravitational binding energy of this halo). However, the interaction of the outflow with the infalling gas results in violent turbulence that also dissipates by cascading down in length scale, thereby additionally heating the halo gas. Potentially, the additional momentum deposition from radiation pressure may help in circumventing the stalling of the CR-driven wind seen for this halo.

Another process that certainly enhances the collimation of the wind into a narrow channel is directly connected to the physics of CR streaming. Recall that the magnitude of wave heating due to CR-excited hydromagnetic waves is proportional to the streaming speed which is equal to the sound speed in our model. Thus, CRs will stream fastest and deposit most of their energy in the high-temperature gas inside the channel, thereby further increasing the temperature there and more importantly the outflow velocity because of enhanced thermal pressure. This is illustrated by the increasing magnitudes of the velocity vectors in the outflow (see top panels of Fig. 9). We find maximum velocities of inside the hot channel, providing the gas with sufficient momentum to overcome the gravitational attraction of the halo.

Moreover, one can see vortex like structures in the velocity field at the sides of the chimneys, where gas is stripped off from the flow and returns to the disc. The vortical motions are presumably created by shear along the sides of the outflow and they converge towards the base of the cavities, where gas is then compressed to high densities, facilitating efficient radiative cooling (similar effects also occur in the case of hot, rising supernova bubbles, see Robinson et al. 2004). This converging gas motion therefore explains the cool torus structure surrounding the disc. As time progresses, the vortical motions will rip the warm chimneys more and more apart so that a warm shell of gas is created around torus and the galaxy.

5.5 Maps of H emission

In the previous section we have seen that the wave heating owing to CR streaming results in characteristic, high-temperature structures. In these structures, the hydrogen gas is highly ionized and therefore allows for optical line emission of recombining electrons. Important in this context is the H line (the dominant line of the Balmer series corresponding to the transition in the hydrogen atom), because hydrogen is the most abundant element in the Universe and thus ubiquitous in astronomical spectra. H-line emission from hot gas has been detected in many galaxies that exhibit a wind, e.g., in the well-known starburst galaxy M82 (Bland & Tully 1988). Thus, we calculate a map the expected H emission in our CR streaming and advection simulation to see how well they compare to these observations.

The H emissivity, , is given by (Spitzer Jr. 1978)

(17)

where eV is the energy of a H photon, is the effective recombination coefficient and denotes the product of the (free) electron number density and the proton number density, respectively. The product gives the number of H photons emitted per second and per cubic centimetre. The temperature dependence of the effective recombination coefficient can be fitted,

(18)

with K (Pequignot et al. 1991, assuming case B). Using equation (18) we can then integrate the emissivity, equation (17), along the line-of-sight.

Figure 11: H emission (edge-on slice) in our CR streaming simulation at for a total halo mass of . We see diffuse H emission in the halo (), which is explained by ionized gas, entrained in the wind.

The result is shown in Fig. 11, where we show the H emissivity in an edge-on view of our simulation at , when the wind is already well-developed. The disc dominates the surface brightness map due to the supernova feedback which heats the gas to sufficiently high temperatures so that the gas becomes ionized. As a result, the product is large and consequently also the recombination rate. Interestingly, there is a halo of H emission with a characteristic surface brightness of that extends to scales of . Such a diffuse emission might be detectable, given reported detections at comparably low levels (Young et al. 1996). This emission is caused by the outflow that transports ionized gas into the halo and beyond. Especially inside the cavities and in their surroundings, where the streaming CRs heat the gas via wave excitation, the gas has a high degree of ionization, ranging from to up to in some places. Unfortunately, the overall gas density in the chimneys is very low, so that their clear structure that we saw in the temperature maps is not visible in H. There is in addition even lower surface brightness diffuse emission with an elliptical morphology as a result of the galactic wind, extending up to .

5.6 Limitations of the model

We need to caution that some of our conclusions depend on certain assumptions used in this work and may not be fully generic to the physics of winds driven by CR streaming. (1) The sub-resolution model for star formation that we adopt here is arguably not well-suited for studying outflows. In particular, the scheme is designed to produce stars in a quiescent mode, i.e., star formation in general only takes place at an average rate. In reality, however, winds are typically observed for starburst galaxies (cf. Veilleux et al. 2005, and references therein), which exhibit high rates of star formation and therefore inject a larger amount of thermal and CR pressure. Apart from that, the star formation scheme does not allow for an explicit treatment of supernova remnants and of the hot, over-pressurized gas bubbles they contain. These bubbles are expected to rise buoyantly in the disc potential of their host galaxy which should open channels along which the wind can escape much easier. (2) The use of such a sub-resolution model with an effective equation of state of the ISM may overestimate the energy transfer from the CRs to the thermal plasma following the excitation and damping of hydromagnetic waves within the dense ISM, that has a comparatively short cooling time, causing a large fraction of this energy to be lost to cooling radiation (as shown in Section 5.2). This overestimate of CR streaming losses may be particularly important in larger haloes with where larger ISM overdensities are reached due to the larger gravity. When explicitly modelling these two phases, it has been shown that CRs do not effectively couple to gas within cool clouds since they do not exert forces inside of cool clouds (Everett & Zweibel 2011). Hence, if CRs were injected into the warm phase of the ISM (for an improved ISM modelling in the simulations), CR streaming could more efficiently drive winds, possibly producing larger mass loading factors for haloes with . (3) To isolate the physics of CR streaming, we have here used an extremely simplified model for a forming and evolving galaxy. Our initial conditions assume spherical symmetry which causes the gas accretion and the associated ram pressure to be initially spherically symmetric. However, a galaxy that forms in a cosmological environment accretes the gas in part along filaments that have a much smaller angular covering factor, hence implying a ram pressure that varies with angle and could potentially modify the resulting wind morphology. This effect is particularly important for collimated winds in larger galaxy haloes with . Those should feel a substantially reduced amount of ram pressure in the wind direction since filaments are generally not aligned with the disc’s angular momentum axis. (4) In our approach, we use a simplified treatment of the CR streaming speed and equate it to the local sound speed. However, a realistic value of should depend on the pre-existing wave level, thermal and CR energy density, magnetic topology, and the damping mechanism of the Alfvén waves, to name a few. (5) We do not follow MHD; in particular, we assume the existence of locally open field lines along which CRs can escape into the halo. While only 10% of the galactic SN remnants will create flux tubes that reach heights above the disc, SN remnant lobes are expected to overlap, suggesting that at every location there is sufficient supply of open field lines (Breitschwerdt et al. 1993). This lets this assumption appear plausible (see also Appendix A). (6) We neglect diffusive shock acceleration of CRs in winds from high-mass Wolf-Rayet stars. The recent detection of TeV-gamma rays from the young stellar cluster Westerlund 2 is compatible with the hypothesis of pion-decay emission of relativistic nuclei interacting with a close-by molecular cloud complex (HESS Collaboration et al. 2011). If confirmed, that would imply that CR feedback could already be an important mechanism in dispersing molecular clouds when the most massive stars turn on, long before the first star became a supernova, making it possibly a competing mechanism to radiation pressure in launching outbursts.

6 Discussion and Conclusions

6.1 Properties of winds driven by CR streaming

Using numerical simulations and analytical arguments, we show in this paper that CR feedback is able to power large-scale galactic outflows. These are not only able to escape the galactic disc, but also its host dark matter halo. In particular, we demonstrate that CR streaming relative to the rest frame of the gas is necessary to drive winds through CRs. Solely considering advective CR transport, i.e., employing the flux-freezing approximation would not result in galactic outflows. In this advective-only case, the CR pressure does volume work on the gaseous discs (in the low-mass haloes), causing them to become more dilute which increases the thermal cooling time. As a result, the SFR declines strongly. After some time, the CR pressure is dissipated such that the gas can collapse again, triggering another “cycle” of star formation and causing an oscillatory mode of star formation with substantial suppression of the stellar mass formed by a factor of 8 for a halo of mass (Jubelgas et al. 2008). In this case, most of the CR energy is eventually transferred to the thermal gas, which radiates it away. In contrast, allowing for CR streaming in the rest frame of the gas down their pressure gradient enables CRs to quickly escape their acceleration sites in the dense, star forming regions, and to move into the CR-dilute, low density gas above the disc (or in between the cold phase of the ISM). There, the CR cooling rates (due to Coulomb or hadronic interactions) are significantly reduced since all these processes scale with the gas density. Thus, a big fraction of their initial energy can be used to power the wind. Moreover, launching a wind from a low-density gas is considerably easier than from the dense phase, because the smaller column density of the material that needs to be removed. Star formation is still suppressed in comparison to the case without CR feedback, but less in comparison to the advection-only case. Remarkably, the star formation rate is much smoother and does not show the oscillatory behaviour discussed above.

We find that CR streaming powers winds most effectively in small dwarf haloes where the wind speed exceeds the escape velocity for host halo masses . Below this mass threshold, we find that the wind velocity increases as a function of host halo mass. We can trace this behaviour back to the potential difference the wind has to climb until it can escape, and which is effectively set by the scale height of the disc. On dwarf scales, the disc scale height is increasing with circular velocity so that the system is approximately scale invariant. In larger galaxies with , the disc height is set by the effective equation of state of the ISM and roughly constant with halo mass. This breaks scale invariance, causes an increasing potential difference for larger haloes, and renders CR streaming not to be powerful enough to be solely responsible for galactic winds. While CR-driven winds are approximately spherical in dwarf haloes (), the denser discs in larger haloes are able to shield the solid angles subtended by the disc and hence focus the wind into bi-conical outflows around the angular momentum axis of the disc.

In a halo of mass , the relative baryonic mass loss amounts to () of the initial gas mass contained inside the halo’s virial radius, assuming a CR acceleration efficiency of (). Similarly, the mass loading factor of the wind increases towards low-mass haloes, reaching values as great as for a halo of . The reason for these increased mass losses towards smaller halo masses are the shallower gravitational potential wells of these systems. However, despite the inability of CR streaming to drive winds in large galaxies that escape the halo, we nevertheless find powerful fountain flows. While we expect a smaller impact of CR streaming towards larger halo masses , we caution that our simplified simulation setting may artificially weaken the CR streaming efficiency of driving winds in large galaxies. First, our employed sub-resolution model with an effective equation of state of the ISM may overestimate the CR energy losses in the ISM due to Alfvén wave excitation since larger ISM overdensities are reached in larger galaxy haloes due to their greater potential depths. Second, a cosmologically forming galaxy accretes the gas in part along filaments that have a much smaller angular covering factor. Since filaments are generally not aligned with the disc’s angular momentum axis, this reduces the accretion ram pressure that collimated winds (aligned with the angular momentum axis) have to exert -work on. Future work that realistically models the ISM in a galaxy forming in a cosmological setting is required to quantify these effects.

Another unique property of CR-driven winds is the effective heating of the outflow by dissipating Alfvén waves that are excited by streaming CRs. Because the wind speeds increase as , we also see an increased rate of heating with increasing halo mass that results in temperatures of the halo gas of . We predict outflows in haloes with masses to reach temperatures in excess of  K so that the outflow is expected to emit thermal bremsstrahlung emission.

We would like to emphasize that neither the classic energy-driven nor the classic momentum-driven wind fully captures the physics of CR streaming-driven winds. The classical picture assumes a pre-loading of the wind with energy or momentum. Instead, the physical mechanism of CR-driven winds implies wave excitation by the CR streaming instability that transfers both, energy and momentum to the gas. This manifests itself through the continuous Alfvén-wave heating term in the energy equation and the term in the momentum equation. CR-driven winds are a two-stream phenomenon composed of the gas flow (carrying momentum) and the fast CR stream (exerting -work on the gas) which continuously re-powers the wind with energy and momentum during its ascent in the gravitational potential.

6.2 Comparing different physical mechanisms for galactic winds

Figure 12: Contours of for cold ( K) wind models driven by a combination of radiation and ram pressure (equation (12) of Sharma & Nath 2012). Different symbols indicate observations of winds with (partly) considerable mass loadings (Martin 1998; Genzel et al. 2001; Schwartz & Martin 2004; Rupke et al. 2005; Martin 2005; Weiner et al. 2009). In these models, the region to the lower right of these contour lines allows for powerful winds driven by ram and radiation pressure, ejecting the gas from the halo. In the region to the upper left, these two processes are not able to drive galactic winds, leaving the observational data in this region unexplained. The range of SFRs encountered in our CR streaming simulations is shown with red lines and indicates the importance of this process in driving winds, especially for galaxies showing smaller values of and SFRs.

We now scrutinize different physical mechanisms responsible for driving galactic winds and put CR streaming – the mechanism studied here – into the context of other plausible processes, namely winds driven by radiation or ram pressure. Generally, it is believed that galaxies have winds when their star formation per unit area exceeds a threshold value of . In our simulation of CR-driven winds, the SFR densities are lower than this threshold. Therefore, this process can address winds with low SFRs.

In Fig. 12, we show contours of in the parameter space spanned by and SFR for cold ( K) wind models driven by a combination of ram pressure and radiation pressure (Sharma & Nath 2012). According to these models, the region of parameter space to the lower right of these contour lines allows for powerful winds driven by ram and radiation pressure that are able to eject the gas from the halo, i.e., it is easier to drive winds for a smaller gravitational potential well (smaller ) and a greater total power input to the system (larger SFR). In contrast, in the region to the upper left, these two processes are not able to drive galactic winds and demonstrate the failure of ram and radiation pressure driven wind models to explain the wind observations in starburst dwarfs (Martin 1998; Schwartz & Martin 2004). Our simulation results for winds driven by CR streaming are shown with red lines (indicating the range of SFRs) and highlight the importance of CR streaming in driving winds, especially for galaxies showing small values of and SFRs.

Good examples are the winds observed in the dwarf starburst galaxies NGC 4214 ( and ; Walter et al. 2001; Martin 2005) and NGC 1569 ( and ; Martin et al. 2002; Stil & Israel 2002). The latter shows direct evidence for a metal-enriched wind from a dwarf starburst galaxy with a mass loading factor of 9 that carries nearly all the metals ejected by the starburst. Both galaxies populate the no-wind region considering only ram and radiation pressure. These dwarf starbursts are close to our halo in the parameter space shown in Fig. 12. The CR-driven bi-conical outflow in this galaxy resembles the observed winds in NGC 4214 and NGC 1569, consistent with a similar physical origin.

While wind velocities and masses can be reliably determined for neutral phases (as it is the case for the two examples discussed), and for H-emitting warm ionized gas (as done in Martin 1998), the velocity of the hot gas ( K) is more difficult to determine. However, it is believed to be higher than that of the neutral gas, with velocities of order its own thermal sound speed of several hundred to a thousand kilometres per second. How does a neutral cloud (or some of its material) get entrained into the wind, in particular for a CR-streaming driven wind? Such a wind is necessarily magnetized with contributions from open field lines (extending from the disc into the halo), flux frozen disc fields, as well as self-amplified field by means of CR streaming. This magnetized super-Alfvénic and supersonic wind (with respect to the gas in the cloud that exhibits comparably low temperatures and sound speeds) impacts the neutral cloud and forms a contact discontinuity. Independent of the degree of magnetization, the cloud sweeps up enough wind magnetic field at the interface to the wind to build up a dynamically important sheath (Gregori et al. 2000; Lyutikov 2006; Asai et al. 2007; Dursi & Pfrommer 2008; Pfrommer & Dursi 2010).

The field strength of this magnetic draping layer is set by a competition between ‘ploughing up’ and slipping around of field lines, maintaining a magnetic energy density in steady state that is of order the wind ram pressure seen by the cloud. This is an inherently three-dimensional problem where the third dimension is essential in allowing for field lines to slip around the obstacle. This strongly magnetized layer modifies the dynamics of the cloud, potentially suppressing hydrodynamic instabilities and mixing (Dursi 2007), and changing the geometry of stripped material. Most importantly, magnetic draping accelerates the cloud into the direction of the wind through magnetic tension of draped field lines that are anchored in the hot phase of the wind, streaming already ahead of the cloud (Dursi & Pfrommer 2008). Note that these authors study the opposite problem of a cloud moving through a magnetized wind (which can be transformed by means of a Galilean transformation to our problem at hand) and find that the deceleration through magnetic tension is always more important, by a factor of , than for the case of highly turbulent (Reynolds number ) hydrodynamic drag. The terminal velocity that the cloud can reach should be a sizable fraction of the wind velocity, and may depend on magnetic geometry and the magnetization of the wind.

This emerging picture suggests that different physical processes dominate the launching and powering of the wind, with a relative contribution that depends on the SFR and (Hopkins et al. 2011; Sharma & Nath 2012). For increasing values of SFR and , we seem to have first CR streaming, then ram pressure, and finally radiation pressure as the dominant agent in powering winds. We conclude that (1) CR-driven winds are most effective in dwarfs, ejecting a large fraction () of metal-enriched baryons into the IGM, and (2) CR Alfvén-wave heating is especially efficient in heating a fountain flow in high-mass galaxies. While we expect that our results for this physical scenario are robust irrespective of redshift or environment, we caution that much more work is needed to firmly confirm this picture, especially given the limitations in modelling these processes at the present time (see Section 5.6).

6.3 Cosmological and observational implications of CR-driven winds

In hierarchical structure formation, dwarf galaxies represent the building blocks of larger galaxies. CR-driven winds appear to be very efficient in removing baryons from small haloes, hence this process should be the dominant feedback mechanism during the formation of these dwarfs, i.e., shortly before and after reionization. We expect the successive mergers of dwarfs to be less gas rich than without this feedback mechanism and to form less stars at high redshift. As galaxy systems grow larger, the ejected gas may be accreted later, but, due to the modified thermal history and angular momentum, may end up at different places within the larger galaxies. Hence, it is not inconceivable that CR-driven winds in small-scale haloes have an indirect effect up to scales of at the knee of the luminosity function; but this needs to be carefully studied in future work.

Additionally, we expect our CR-driven winds to influence the early history of metal enrichment during the epoch of dwarf formation because of their considerable mass loading factors. This important role of low-mass systems in the metal enrichment of the IGM is in line with expectations from other authors (Nath & Trentham 1997; Madau et al. 2001; Samui et al. 2010).

What about specific observational predictions of CR-driven winds? The bi-conical structure of the hot gas cavities that formed in the halo has some similarities with H-structures, e.g., in the well-known starburst galaxy M82 (Bland & Tully 1988) or with the dwarf starburst galaxy NGC 1569 (Martin et al. 2002). However, the bi-conical structure observed in the temperature field has no direct resemblance in our simulated H emission maps, due to the low density of the gas in these cavities. This is the same problem that commonly makes observations of winds difficult, which face the difficulty of detecting the weak emission from a wind that is in general hot and dilute (see Veilleux et al. 2005, and references therein). Our simulations are not capable of describing the full spectrum of fluid phenomena that are likely to be associated with these hot cavities. Fluid instabilities of Kelvin-Helmholtz and Rayleigh-Taylor type (Robinson et al. 2004) as well as thermal instabilities should be important because of the density contrast around the chimneys and the shear along their sides. In addition, there will be entrainment of denser clouds from the interstellar medium when the hot wind pushes on them (Mac Low & McCray 1988). All these phenomena will influence the morphology of the cavities and could enhance the density locally, so that their signature might indeed be visible in H. However, we lack the necessary resolution to see these effects in our simulations and the SPH method may also significantly contribute to this limitation, as recent studies found (Sijacki et al. 2011). Moreover, magnetic fields that are confined within the hot chimney gas and rise up along with it are of particular importance. They will support the chimneys against the hydrodynamical instabilities and shear which try to tear them apart (Robinson et al. 2004; Dursi 2007; Ruszkowski et al. 2007; Dursi & Pfrommer 2008). Nevertheless, our simulations predict a low but observable level of extended diffuse H emission in haloes of mass . In fact, extended H emission forming an H-halo, that extends several kiloparsecs away from the disc, has been reported for several galaxies. The strength and spatial extent of the emission seem to be positively correlated with the SFR in the galaxy (e.g., Rand 1996). This highlights the importance of studying CR-driven winds including a better description of star formation and stellar feedback.

Some galaxies also exhibit radio haloes with an overall correlation between H and radio emission. Such a correlation would naturally arise when CRs are convected into the halo, as is the case in our simulations. To discriminate a pure advective scenario from our model, the spectral age of CR electrons could potentially be decisive (see Heesen et al. 2009, for such a measurement). Including CR streaming should decrease the spectral age since streaming electrons escape faster into the halo. However, the adiabatic expansion history in both scenarios needs to be identical to allow for a fair comparison. Future work that includes MHD will address those questions and allow for more definite predictions of the expected radio synchrotron emission in this scenario.

Acknowledgments

We thank Crystal Martin and Peng Oh for useful discussions. We also thank the referee for a thorough reading of the manuscript and for his constructive comments. This research was supported in the framework of the DFG Forschergruppe 1254 “Magnetization of Interstellar and Intergalactic Media: The Prospects of Low-Frequency Radio Observations”. C.P. gratefully acknowledges financial support of the Klaus Tschira Foundation and would furthermore like to thank KITP for their hospitality during the conference and workshop on “First Light and Faintest Dwarfs” (which was supported by the National Science Foundation under Grant No. NSF PHY05-51164). B.B.N. thanks the Alexander-von-Humboldt foundation for supporting his visit to the Max-Planck-Institute for Astrophysics (MPA). B.B.N. also thanks the MPA for the hospitality during his visit. V.S. acknowledges support by the DFG Research Centre 881 “The Milky Way System”. The rendered plots in this work were produced using the visualization tool SPLASH by Price (2007). The simulations were performed at the Rechenzentrum Garching (RZG) of the Max Planck Society as well as using resources of the bwGRiD project8.

Appendix A Implementation of CR streaming

Here, we discuss our implementation of CR streaming into the Lagrangian SPH code Gadget-2. First, we translate the transport equations for CR energy and number density into Lagrangian form and then cast them into the language of SPH. Our starting point are the evolution equations for CR energy density and number density which can be derived from the CR transport equation (e.g., Skilling 1971, 1975),

(19)
(20)

where is the gas speed, is the streaming speed of CRs and is the CR pressure. Here and in the reminder of this derivation, we neglect all CR source terms as well as CR diffusion. We model CR sources and sinks explicitly in our numerical CR formalism (Pfrommer et al. 2006; Enßlin et al. 2007; Jubelgas et al. 2008). For the limiting case of strong CR-wave scattering (the Bohm limit), the effective diffusion velocity is usually much slower than the CR streaming velocity, justifying our approximation. As discussed in Section 2, the streaming speed is assumed to be proportional to the local sound speed, , and anti-parallel to a unit vector along the CR pressure gradient in our model, i.e.,

(21)

where is a proportionality constant for a plasma where the thermal pressure dominates over the magnetic pressure. The CR streaming velocity has been taken to be opposite to the CR pressure gradient rather than the gradient of the CR number density of a CR energy interval (which would be the formal criterion for evaluating the direction of the CR streaming velocity). For power-law momentum distributions, , and hence the gradients of and point into the same direction. If this was not the case, i.e., if different energy regimes dominated and , then we would be interested in the energy regime that dominates the CR pressure (which is responsible for driving the wind) and hence in the streaming direction given by the CR pressure gradient.

Since we do not follow MHD in our simulations, we are unable to project the CR pressure gradient onto the magnetic field lines to determine the formally correct streaming direction. However, as we will now argue, our approach should correctly capture the main aspects of CR streaming and the associated galactic mass loss even without magnetic fields. In hydrostatic equilibrium, the gradient of the sum of thermal and CR pressure is determined by the gradient of the gravitational potential,

(22)

Above the disc, when wave heating induced by streaming CRs provides substantial heat input, the CR pressure distribution determines the total pressure distribution. For CRs to escape the galactic disc via streaming requires locally open field lines (that are stretched to large heights above the disc). While only 10% of the galactic SN remnants will create flux tubes that reach those altitudes, SN remnant lobes are expected to overlap, suggesting that at every location there is sufficient supply of open field lines (Breitschwerdt et al. 1993). The field topology at the disc-halo interface and beyond will be determined by the buoyancy of the CR component that drive a large-scale upward convection and possibly a galactic mass loss (Parker 1966). Hence, the topology of the open field line structure is determined by the CR pressure stratification which itself is shaped by the gravitational potential gradient. This should establish a field topology consistent with our assumptions, at least when averaged over sufficiently large scales so that small-scale fluctuations are averaged out.

To derive the Lagrangian form of the CR evolution equations, we define a Lagrangian time derivative,

(23)

and introduce the specific CR energy, , and CR particle number, , through

(24)
(25)

After using the continuity equation for the gas, , we arrive at the Lagrangian form of the CR evolution equations,