Cosmic rays can drive strong outflows from gas-rich high-redshift disk galaxies
We present simulations of the magnetized interstellar medium (ISM) in models of massive star forming () disk galaxies with high gas surface densities () similar to observed star forming high-redshift disks. We assume that type II supernovae deposit 10 per cent of their energy into the ISM as cosmic rays and neglect the additional deposition of thermal energy or momentum. With a typical Galactic diffusion coefficient for CRs () we demonstrate that this process alone can trigger the local formation of a strong low density galactic wind maintaining vertically open field lines. Driven by the additional pressure gradient of the relativistic fluid the wind speed can exceed , much higher than the escape velocity of the galaxy. The global mass loading, i.e. the ratio of the gas mass leaving the galactic disk in a wind to the star formation rate becomes of order unity once the system has settled into an equilibrium. We conclude that relativistic particles accelerated in supernova remnants alone provide a natural and efficient mechanism to trigger winds similar to observed mass-loaded galactic winds in high-redshift galaxies. These winds also help explaining the low efficiencies for the conversion of gas into stars in galaxies as well as the early enrichment of the intergalactic medium with metals. This mechanism can be at least of similar importance than the traditionally considered momentum feedback from massive stars and thermal and kinetic feedback from supernova explosions.
Subject headings:galaxies: general — galaxies: ISM — ISM: magnetic fields — cosmic rays
The universal rate for the conversion of gas into stars in galaxies peaks at redshifts z (e.g. Hopkins & Beacom, 2006). At this epoch star forming galaxies drive powerful galactic winds which can transport a significant fraction of the gas away from the central galaxy making it temporarily unavailable for star formation (e.g. Pettini et al., 2000; Shapley et al., 2003; Weiner et al., 2009; Steidel et al., 2010; Genzel et al., 2011; Newman et al., 2012a, b). Spatially resolved high-redshift observations indicate that these winds are launched directly from the sites of the - typically strongly clustered - star formation (Genzel et al., 2011; Newman et al., 2012a, b). The estimated outflow rates can be several times higher than the star formation rates (SFR). The resulting high mass loading indicates that - even at the peak epoch of cosmic star formation - the amount of gas expelled from the galaxies is comparable to the amount of gas converted into stars inside the galaxies.
The direct observational evidence for inefficient conversion of gas into stars is supported by indirect constraints from halo abundance matching techniques. Here the galaxy formation efficiency can be defined as the fraction of the stellar mass of a galaxy to the total available baryonic mass of its host dark matter halo. In a concordance -CDM cosmology this efficiency peaks - almost independent of redshift - for galaxies in halos of about and never exceeds (e.g. Moster et al., 2010; Behroozi et al., 2010; Guo et al., 2010; Moster et al., 2013). Therefore at least (significantly more in halos of higher as well as lower mass than ) of the baryonic material is never converted into stars, eventually due to powerful galactic winds. It is plausible that the main physical processes responsible for driving the outflows also regulate the efficiency with which the available gas is converted into the stellar components of galaxies in the Universe.
In connection to the evolution of stellar populations a number of physical processes are - in principle - energetic enough to expel gas from star forming galactic disks. Besides AGN for high mass galaxies (Croton et al., 2006), type II supernovae have long been considered the most promising candidates, in particular for lower mass galaxies (Larson, 1974; Dekel & Silk, 1986). Although the amount of energy per event is significant the thermal energy is mainly deposited at the sites of star formation, i.e. dense molecular clouds. Here the cooling times are very short and the energy can be efficiently radiated away making it difficult but not impossible to drive large scale galactic winds (for recent discussions see e.g. Brook et al., 2011; Dalla Vecchia & Schaye, 2012). However, even before the supernova explosions, the momentum and energy input from massive stars in the form of stellar winds and stellar luminosity is significant and might support the wind driving (Murray et al., 2005; Hopkins et al., 2012; Agertz et al., 2013).
In this paper we focus on a separate mechanism: the formation of large scale magnetized galactic winds driven by cosmic rays. As supernovae drive strong shocks into the interstellar medium some fraction of the explosion energy is consumed to accelerate ionized particles to relativistic energies which are then injected into the ISM as cosmic rays (CR) (Krymskii, 1977; Bell, 1978; Blandford & Ostriker, 1978). This relativistic fluid is coupled to the galactic magnetic field and - in particular the hadronic component - is less prone to energy losses than the gaseous component of the ISM. Analytic estimates and numerical experiments without or only approximate inclusions of galactic magnetic fields clearly indicate that CRs can help driving large scale galactic winds (Breitschwerdt et al., 1991, 2002; Strong et al., 2007; Everett et al., 2008; Uhlig et al., 2012; Dorfi & Breitschwerdt, 2012). However, CRs are strongly coupled to magnetic fields whose evolution should be followed in a self-consistent way. Hanasz et al. (2004, 2009); Siejkowski et al. (2010); Kulpa-Dybeł et al. (2011) have shown that CRs promote buoyancy effects in the interstellar medium, leading to the break-out of magnetic fields from galactic disks (Parker, 1992) and, at the same time, to magnetic field amplification by CR-driven dynamo action. Plausibly, such processes are also relevant for star forming galaxies at high redshift which are observed to have significant magnetic fields at the level of tens of . (Bernet et al., 2008). Recent observations even demonstrate the existence of large magnetic fields up 50 kpc away from the galaxy indicating strong large-scale magnetized winds (Bernet et al., 2013). In this letter we present a three-dimensional full MHD simulation of a massive gas-rich disk galaxy (section 2) and follow the formation of large scale magnetized winds as a dynamic response to the injection of CRs (section 3). This is considered as a proof of principle for the importance of this physical process. The most important implications are presented in section 4.
2. Numerical setup
For the simulations we use the PIERNIK MHD code (Hanasz et al., 2010), a grid-MHD code based on the Relaxing TVD (RTVD) scheme by *jin-xin-95 and *2003ApJS..149..447P. The induction equation, including the Ohmic resistivity term, is integrated with a constraint transport (CT) algorithm (Evans & Hawley, 1988). The original scheme is extended to include dynamically independent, but interacting fluids: thermal gas and a diffusive CR gas, described within the fluid approximation. (Hanasz & Lesch, 2003). We incorporate selfgravity of interstellar gas and gravitational potential is obtained by solving the Poisson equation inside the computational domain with an iterative, multi-grid solver (Huang & Greengard, 1999) combined with a multipole solver (James, 1977) to properly treat the gravitational potential at ’isolated’ boundaries.
We assume a fixed gravitational field due to the stellar disk and the dark matter halo and compute the gravitational potential using the model of Allen & Santillan (1991) with within and . We neglect the contribution of a central bulge. Fresh gas is supplied to the disk at a fixed rate of following the initial gas density distribution. To simplify the setup the gas is added directly at the disk plane. This is clearly a simplified model. However, another simplified way of treating the gas supply by spherical accretion is much more difficult to control. Even in this case we expect that the highly collimated winds forming in our simulation would punch through spherically accreting gas not changing our conclusions. The more realistic alternative of gas accretion along filaments would provide fresh gas, carrying high angular momentum, at disk peripheries and is not expected to change the wind properties significantly.
We construct a three-dimensional gas distribution (e.g. Ferriere, 1998). At every time step we add gas at the given inflow rate within the disk volume. The disk collects gas until it becomes gravitationally unstable. We assume that star formation is controlled by a star formation efficiency parameter . In the actual models we assume . The star formation rate per unit volume is computed as
provided that gas density exceeds some threshold density , which we treat as a free parameter. Its value () is adopted to obtain highly localized star formation, and to regulate the overall SFR of the galaxy. We locally deplete the gaseous ISM at the same rate at every timestep.
We assume that one supernova occurs per of gas forming new stars, and that of the explosion energy is used to accelerate CRs. Individual CR particles propagate at relativistic speeds, however fast streaming of CRs along magnetic field lines leads to streaming instabilities (Kulsrud & Pearce, 1969), the generation of small-scale turbulence, and subsequently the scattering of CRs on self-excited turbulence. This implies that a CR pressure gradient term has to be taken into account in the gas equation of motion. The bulk motion of CRs is considered as a combination of diffusive and advective propagation and can be described by the diffusion-advection transport equation. In our model CRs diffuse preferentially along magnetic field lines, while diffusion perpendicular to the magnetic field is significantly much less efficient (see Berezinskii et al., 1990, Chap. 9).
A more elaborated picture of CR propagation should include additional effects, such as energy conversion from CRs to waves, energy sinks for MHD waves due to ion-neutral collisions and nonlinear Landau damping (for a detailed discussion see e.g. Breitschwerdt et al. (2002); Dorfi & Breitschwerdt (2012) and references therein). The above mentioned authors adopt constant CR diffusion coefficients ranging from up to in their 1-D numerical models. The assumption of constant diffusion coefficient implies that the large-scale diffusion velocity may become larger than the Alfvén speed in regions of steep gradients of the CR energy density. This can occur in the disk, and especially around CR production regions, where molecular and neutral gas components dominate. Since Alfvén waves are efficiently dissipated there by ion-neutral damping, the most important contribution to the random magnetic field in the disk is induced by supernova explosions. Therefore ion-neutral collisions should reduce the amplitude of Alfvén waves, leading to enlargement of the mean free path of the CR particles, and therefore to higher diffusion coefficients.
The numerical algorithm of the anisotropic CR propagation, within the framework of staggered mesh MHD code, has been described in Hanasz & Lesch (2003). The values of CR diffusion coefficients, parallel and perpendicular to the magnetic fields adopted for the actual models are , . Initially the toroidal magnetic field pervading the disk has a strength of and a uniform magnetic diffusivity , corresponding to a standard value of turbulent diffusivity of the ISM.
The simulations have been performed at a resolution of grid cells, distributed among equal-sized MPI blocks, in the Cartesian domain spanning a volume of . The disk is placed at the centre of the domain, and the disk plane is parallel to plane of the coordinate system. We impose outflow boundary conditions for the gas component at all domain boundaries. Fixed boundary conditions () on external domain boundaries are assumed for the CR component.
Initially the gaseous disk collects gas at the presumed global infall rate until it becomes locally gravitationally unstable. Supernovae start to explode and deposit CRs in the ISM after the gas density exceeds the critical value. After about the disk reaches an equilibrium state with a star formation rate at a level of . A typical snapshot of the system after of evolution is shown in Fig. 1. Most of the supernovae activity is confined to isolated regions in kpc-sized dense gas clouds (upper right panel). These regions can be also identified as spots of high CR energy density apparent as dark brown and black patches in the face-on map (lower right panel of Fig. 1). One can identify about discrete star formation regions with CR energy densities exceeding dropping to at larger distances away from the disk (lower panels of Fig. 1). The distribution of the CR energy density in the galactic halo is highly non-uniform. Sharp edges of CR-populated regions can be identified with similar structures in the maps of vertical mass flux and vertical magnetic field component shown in Fig. 2.
The vertical streams of rarefied gas visible in gas density distribution (upper panels of Figs. 1 and 2) are accelerated, by CRs, to high velocities (several ). The streams can extend several tens of above and below the disk plane and significant fraction of the outflowing gas has velocities above escape velocity and will be able to leave the galaxy altogether.
Maps of the mass flux (mid panels of Fig. 2) show the bimodal nature of the outflow perpendicular to the disk plane with peak values up to (the color scale of the mass flux panels is saturated at only to show the wind structure far from the disk plane).
Streams of gas emanating from a single star forming region have a large cross-section, visible at the horizontal slice of at . Individual SF regions generate outflows of on average and form streams of about in radius on both sides of the disk. In some cases the streams from neighboring SF regions merge. The horizontal area of a single stream is a few . This area multiplied by gives a number consistent with through the stream cross-section area, even though the average flux at may be half of the peak value for two outflows on both sides of the disk. The mass loaded wind is enriched with fresh CRs and is highly magnetized. Field strengths exceeding can be reached naturally at distances of away from the central galaxy. Magnetic flux tubes coherent over such large distances would be directly detectable with Faraday rotation measurement techniques.
To quantify the vertical structure of the CR-driven wind we plot in Fig. 3 the total vertical mass flux as a function of distance from the galactic disk from to . In general the flow patterns are symmetric and outflow dominates inflow by an order of magnitude. Total outflow rates decline from close to the disk plane to at . A more detailed inspection of wind density-velocity distribution (Fig. 4) shows that the density of the CR driven wind gas is typically below and the high-velocity gas () escapes at densities H atoms .
In Fig. 5 we show the evolution of the star formation rate and the integrated outflow rates at different altitudes at , , and above the disk plane. After about 400 Myr the star formation rate settles to a value of . The mean surface density of disk gas within the radius 10 kpc reaches an equilibrium value of and SFR surface density is . These values are in good agreement with nearby highly star forming galaxies and typical massive high-redshift disks (Kennicutt, 1998; Tacconi et al., 2013). At about away from the disk plane the mass outflow rate matches the star formation rate and the galaxy average mass-loading becomes of order unity. Closer to the disk plane the mass-loading is higher and further away from the disk plane it is still significant. It only decreases to values of .
4. Discussion and conclusions
We have demonstrated as a proof of principle that the injection of only of SN energy in the form of CRs and neglecting the thermal and kinetic energy input is sufficient to drive a large scale galactic wind in a gas-rich and highly star forming disk with properties similar to typical star forming high-redshift galaxies. The additional pressure gradient of the relativistic fluid - which in contrast to heated dense gas cannot easily dissipate its energy away - drives the formation of a strong bi-polar galactic wind with velocities exceeding . CRs can easily escape far from dense regions with almost negligible energy losses and deposit their energy and momentum in rarefied medium. This process is supported by the CR driven break-out of field lines whose vertically open structure is maintained by the wind. Cosmic rays can rapidly diffuse along these field lines far into the galactic halo.
To check consequences of our assumption of constant diffusion coefficients we examined the ratio of CR diffusion speed to the Alfvén speed over the simulation volume. As shown in Fig.1 the CR energy distribution is very smooth everywhere except the close to the SF regions and some apparent shock-like structures, especially when compared to the gas density distribution. Vertical stratification of the CR energy distribution can be observed only far from star formation regions and the corresponding scale height is of the order of a few kpc. The smooth distribution of CRs is a consequence of a highly irregular vertical magnetic field component, guiding the field aligned CR diffusion in the vertical direction. We have found that in a major part of the simulation volume. Consequently, in regions of galactic disk with high gas density. As we mentioned already, we consider the limitation of the diffusion velocity by Alfvén waves less restrictive in the disk volume, because neutrals are supposed to reduce the level of CR self-excited Alfvénic turbulence. Moreover, we note that steep gradients of CR energy density around star formation regions provide a specific feedback mechanism. High production rates of CRs implies fast expansion of overpressured regions, and subsequently enhanced CR advective expansion leads to a reduction of the CR energy density gradient. The fast expansion of CR overpressured bubbles increases the Alfvén speed locally, leading to a reduction of the ratio .
The CR driving is so significant that the mass outflow rate can become of the same order as the star formation rate in the galactic disk, even in our simplified setup where the disk plane is more or less treated as an inner boundary condition and thermal as well as kinetic feedback from stellar evolution and supernovae have been neglected entirely. Based on our simulations we can conclude that relativistic particles accelerated in supernova remnants in combination with a strong magnetic fields (typical for high-redshift galaxies (Bernet et al., 2008)) provide a natural and efficient mechanism to help explaining the ubiquitously observed mass-loaded galactic winds in high-redshift galaxies (e.g. (Shapley et al., 2003; Newman et al., 2012a)) as well as the highly magnetized medium surrounding these galaxies (Bernet et al., 2013). The efficiency - in terms of mass loading - of this wind driving process appears to be comparable to momentum and energy driving from stellar evolution and supernovae explosions (e.g. Hopkins et al., 2012; Dalla Vecchia & Schaye, 2012; Agertz et al., 2013) and requires further investigation.
- Agertz et al. (2013) Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25
- Allen & Santillan (1991) Allen, C., & Santillan, A. 1991, Revista Mexicana de Astronomia y Astrofisica, 22, 255
- Behroozi et al. (2010) Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379
- Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147
- Berezinskii et al. (1990) Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin, V. S. 1990, Astrophysics of cosmic rays (Amsterdam: North-Holland, 1990, edited by Ginzburg, V.L.)
- Bernet et al. (2013) Bernet, M. L., Miniati, F., & Lilly, S. J. 2013, ArXiv e-prints
- Bernet et al. (2008) Bernet, M. L., Miniati, F., Lilly, S. J., Kronberg, P. P., & Dessauges-Zavadsky, M. 2008, Nature, 454, 302
- Blandford & Ostriker (1978) Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29
- Breitschwerdt et al. (2002) Breitschwerdt, D., Dogiel, V. A., & Völk, H. J. 2002, A&A, 385, 216
- Breitschwerdt et al. (1991) Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, A&A, 245, 79
- Brook et al. (2011) Brook, C. B., et al. 2011, MNRAS, 415, 1051
- Croton et al. (2006) Croton, D. J., et al. 2006, MNRAS, 365, 11
- Dalla Vecchia & Schaye (2012) Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, 426, 140
- Dekel & Silk (1986) Dekel, A., & Silk, J. 1986, ApJ, 303, 39
- Dorfi & Breitschwerdt (2012) Dorfi, E. A., & Breitschwerdt, D. 2012, A&A, 540, A77
- Evans & Hawley (1988) Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659
- Everett et al. (2008) Everett, J. E., Zweibel, E. G., Benjamin, R. A., McCammon, D., Rocks, L., & Gallagher, III, J. S. 2008, ApJ, 674, 258
- Ferriere (1998) Ferriere, K. 1998, ApJ, 497, 759
- Genzel et al. (2011) Genzel, R., et al. 2011, ApJ, 733, 101
- Guo et al. (2010) Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010, MNRAS, 404, 1111
- Hanasz et al. (2004) Hanasz, M., Kowal, G., Otmianowska-Mazur, K., & Lesch, H. 2004, ApJ, 605, L33
- Hanasz et al. (2010) Hanasz, M., Kowalik, K., Wóltański, D., & Pawłaszek, R. 2010, in EAS Publications Series, Vol. 42, EAS Publications Series, ed. K. Goździewski, A. Niedzielski, & J. Schneider, 275–280
- Hanasz & Lesch (2003) Hanasz, M., & Lesch, H. 2003, A&A, 412, 331
- Hanasz et al. (2009) Hanasz, M., Wóltański, D., & Kowalik, K. 2009, ApJ, 706, L155
- Hopkins & Beacom (2006) Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142
- Hopkins et al. (2012) Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522
- Huang & Greengard (1999) Huang, J., & Greengard, L. 1999, SIAM Journal on Scientific Computing, 21, 1551
- James (1977) James, R. A. 1977, Journal of Computational Physics, 25, 71
- Jin & Xin (1995) Jin, S., & Xin, Z. 1995, Comm. Pure Appl. Math., 48, 235
- Kennicutt (1998) Kennicutt, R. C. 1998, ApJ, 498, 541
- Krymskii (1977) Krymskii, G. F. 1977, Akademiia Nauk SSSR Doklady, 234, 1306
- Kulpa-Dybeł et al. (2011) Kulpa-Dybeł, K., Otmianowska-Mazur, K., Kulesza-Żydzik, B., Hanasz, M., Kowal, G., Wóltański, D., & Kowalik, K. 2011, ApJ, 733, L18
- Kulsrud & Pearce (1969) Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445
- Larson (1974) Larson, R. B. 1974, MNRAS, 169, 229
- Moster et al. (2013) Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121
- Moster et al. (2010) Moster, B. P., Somerville, R. S., Maulbetsch, C., van den Bosch, F. C., Macciò, A. V., Naab, T., & Oser, L. 2010, ApJ, 710, 903
- Murray et al. (2005) Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569
- Newman et al. (2012a) Newman, S. F., et al. 2012a, ApJ, 752, 111
- Newman et al. (2012b) —. 2012b, ApJ, 761, 43
- Parker (1992) Parker, E. N. 1992, ApJ, 401, 137
- Pen et al. (2003) Pen, U.-L., Arras, P., & Wong, S. 2003, ApJS, 149, 447
- Pettini et al. (2000) Pettini, M., Steidel, C. C., Adelberger, K. L., Dickinson, M., & Giavalisco, M. 2000, ApJ, 528, 96
- Shapley et al. (2003) Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65
- Siejkowski et al. (2010) Siejkowski, H., Soida, M., Otmianowska-Mazur, K., Hanasz, M., & Bomans, D. J. 2010, A&A, 510, A97
- Steidel et al. (2010) Steidel, C. C., Erb, D. K., Shapley, A. E., Pettini, M., Reddy, N., Bogosavljević, M., Rudie, G. C., & Rakic, O. 2010, ApJ, 717, 289
- Strong et al. (2007) Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Annual Review of Nuclear and Particle Science, 57, 285
- Tacconi et al. (2013) Tacconi, L. J., et al. 2013, ApJ, 768, 74
- Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., Skory, S., Skillman, S. W., Abel, T., & Norman, M. L. 2011, ApJS, 192, 9
- Uhlig et al. (2012) Uhlig, M., Pfrommer, C., Sharma, M., Nath, B. B., Enßlin, T. A., & Springel, V. 2012, MNRAS, 423, 2374
- Weiner et al. (2009) Weiner, B. J., et al. 2009, ApJ, 692, 187