EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 638 (June 2020) A&A, 638 (2020) A123 Full HTML
Open Access
Issue
A&A
Volume 638, June 2020
Article Number A123
Number of page(s) 19
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201936339
Published online 23 June 2020

© G. Dashyan and Y. Dubois 2020

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

It has by now been well established that feedback – the processes by which star formation and supermassive black hole growth regulate themselves by deposition of energy and momentum in the ambient medium – is a key mechanism to regulate the amount of baryons in galaxies (see Somerville & Davé 2015; Naab & Ostriker 2017, for recent reviews). Additionally, outflows are ubiquitously observed in galaxies (see review by Veilleux et al. 2005) and moderate density regions of the Universe – where stars do not form – are metal-enriched, which suggests that metal-enriched winds can deposit metals into the intergalactic medium (e.g. Booth et al. 2012). However, the processes regulating star formation and driving winds in galaxies, as well as their relative importance, are still a topic of active investigation.

Due to their shallow potential wells and lack of galactic mergers during their evolution, dwarf galaxies (with halo virial mass of Mvir ≤ 1011 M⊙) can act as both important probes for cosmological structure formation, and as astrophysical test-beds for feedback processes. As the least massive and least luminous galaxies, they are both difficult to resolve in cosmological simulations and under-observed. Several puzzles surrounding dwarf galaxies, including the too-big-to-fail problem (Boylan-Kolchin et al. 2011), the missing satellites (Klypin et al. 1999; Moore et al. 1999a), and the cusp-core problems (Moore et al. 1999b), have posed important challenges to our understanding of galaxy formation and dark matter (DM) for DM-only models that can potentially be alleviated by a realistic treatment of baryonic feedback processes (e.g. Pontzen & Governato 2012; Zolotov et al. 2012; Teyssier et al. 2013; Brooks & Zolotov 2014; Chan et al. 2015; Oñorbe et al. 2015).

Supernova (SN) feedback remains the commonly believed dominant feedback mechanism in low mass halos (e.g. Dekel & Silk 1986). The problem is that it is partly implemented in terms of subgrid physics in numerical simulations, which prevents the capture of its exact nature and extent. Ideally, we would like to directly simulate the physics without tuning the parameters to observables, and rather use observations as powerful discriminants for theories. Recent work argues for the need for interstellar medium (ISM) physics beyond SN feedback alone to suppress star formation and generate outflows (Rosdahl et al. 2017; Smith et al. 2019). In addition, resolution is generally too low to capture the detailed evolution of SN explosions (Hu 2019), which has produced a diversity of subgrid workarounds such as delayed cooling (Stinson et al. 2006; Teyssier et al. 2013), kinetic winds (Springel & Hernquist 2003; Dubois & Teyssier 2008), stochastic feedback (Dalla Vecchia & Schaye 2012), or momentum-driven approaches (Kimm et al. 2015), with diverse successes at producing efficient outflows (see Rosdahl et al. 2017, for a comparison of such models). Feedback mechanisms other than SNe are now suspected to play a key role in shaping galactic winds and regulating the baryon content in galaxies (Hopkins et al. 2014). The attention has therefore turned to other processes such as stellar radiation (Hopkins et al. 2011; Agertz et al. 2013; Aumer et al. 2013; Rosdahl et al. 2015; Emerick et al. 2018), the momentum transferred from resonantly scattered Lyman-α photons (Kimm et al. 2018), or cosmic rays (CRs) released in SNe (Wadepuhl & Springel 2011; Booth et al. 2013; Salem & Bryan 2014).

CRs are particles accelerated to relativistic velocities at shock fronts of supernova remnants (SNR; see Caprioli 2015 for a recent review), active galactic nuclei (e.g. Berezinsky et al. 2006), winds from massive stars (Bykov 2014), or structure formation shocks (e.g. Miniati et al. 2001). Observations of local SN remnants suggest that of the order of 10% of the explosion energy is converted to CRs (Morlino & Caprioli 2012; Helder et al. 2013). In the solar neighbourhood, CRs have an observed energy density of ∼10−12 ergs cm−3 (Wefel 1987), which is comparable to the magnetic, turbulent and vertical gravitational energies. That equipartition shows that CRs are important in the local dynamics of the disc. Recent observations suggest that CR energy density exceeds that seen in the Milky Way by orders of magnitude in the starbust galaxy M 82 (Veritas & Acciari 2009; Paglione & Abrahams 2012). The close correlation observed between the radio continuum and the far-infrared emission of spiral galaxies (e.g. de Jong et al. 1985) hints at a correlation between star formation and CR energy density and therefore at a powerful self-regulation mechanism.

The protons of a few GeV, where most of CR energy density resides, are essentially collisionless. Their coupling with the gas is mediated by the ambient magnetic field and can be described by a fluid picture that macroscopically describes the energy and momentum exchange between CRs and thermal gas, mediated by the gyroresonant streaming instability (Kulsrud 2005; Zweibel 2017). CRs have several properties which make them a promising candidate for driving galactic winds: they provide an additional pressure gradient that helps lifting thermal gas; due to its relativistic nature, the CR component has a softer equation of state (γCR = 4/3) than the thermal component (γ = 5/3) and therefore suffers less losses upon adiabatic expansion; the CR energy has a longer cooling time than that of the thermal component; finally, CRs diffuse along magnetic field lines and thereby move relative to the thermal component and leak out from their production sites, that is from dense and efficiently cooling star-forming regions to more diffuse regions.

The role of CRs on galaxy evolution and their ability to drive winds have been a subject of interest in both analytical (Ipavich 1975; Breitschwerdt et al. 1991; Recchia et al. 2016; Mao & Ostriker 2018) and numerical studies. In hydrodynamical simulations without magnetic field evolution, Uhlig et al. (2012) found that CR streaming – assuming a streaming velocity proportional to the sound speed – drives powerful and sustained winds in galaxies with virial masses Mvir < 1011 M⊙. Booth et al. (2013), performing simulations of Small Magellanic Cloud and Milky Way sized disc galaxies, showed that the star-formation rate (SFR) was suppressed for both masses when accounting for isotropic CR diffusion (i.e. diffusion of CRs that is independent of the magnetic field lines). With similar assumptions on the propagation of CRs, Salem & Bryan (2014) found that CRs thicken gaseous disc, suppress the SFR and produce mass loading factors (i.e. the ratio of the outflow rate to the star-formation rate) above unity in a Milky Way mass galaxy.

In order to include the anisotropic nature of CR diffusion into simulations – that is the fact that CRs diffuse along magnetic field lines preferentially – the evolution of magnetic fields would have to be modelled since CRs gyrates around magnetic field lines, CR diffusion and streaming are parallel to them, and the streaming velocity depends on the local Alfvén velocity, which depends on the amplitude of the magnetic field. Using magneto-hydrodynamical (MHD) simulations, Hanasz et al. (2013) demonstrated that even with anisotropic diffusion and CR feedback alone, efficient large-scale winds can be driven out of galaxies. Pakmor et al. (2016) showed that winds develop significantly later with anisotropic diffusion compared to isotropic diffusion. In simulations of a stratified disc, Girichidis et al. (2016) showed that including CRs thickens the galactic disc and leads to the formation of a warm and neutral galactic atmosphere providing a mass reservoir for galactic winds and outflows, as opposed to feedback without CRs that produces hot and ionised outflows. Ruszkowski et al. (2017) included CR streaming in numerical galactic MHD simulations. They found that the inclusion of CR streaming and anisotropic diffusion can have a significant effect on wind launching and mass loading factors. They also found that the presence of moderately super-alfvénic CR streaming greatly enhances the efficiency of galactic wind driving. However Holguin et al. (2019), taking into account realistic streaming instability suppression in the presence of turbulence that results in higher streaming speeds, found that the mass loading factor drops significantly as the strength of turbulence increases. Butsky & Quinn (2018), simulating a suite of isolated Milky Way-type disc galaxies, found that CR feedback suppresses star formation by supporting thermal gas against collapse, that models with CR transport drive strong outflows and reproduce the multiphase temperature and ionisation structure of the circumgalactic medium (CGM), but that the simulated CGM strongly depends on CR transport model. Girichidis et al. (2018) studied the impact of CR-driven galactic outflows from a multiphase ISM, and found that simulations including CR have denser, colder and slower outflows. Performing computations of resonant Lyman-α radiation transfer through snapshots of a suite of stratified disc simulations, Gronke et al. (2018) found that the absence of CR feedback leads to spectra incompatible with observations due to the smoother neutral gas distribution of CR supported outflows. Jacob et al. (2018) performed a parameter study varying the halo mass, the diffusion coefficient and the CR injection fraction and found that CRs are only able to drive winds in halos with masses Mvir < 1012 M⊙ and that the mass loading factor drops with virial mass. Chan et al. (2019) presented isolated simulations of dwarf and Milky Way-type galaxies including CR feedback. They found that only constant isotropic diffusion coefficients of 3 × 1028 − 29 cm2 s−1 reproduce observations of γ-ray emission. With advection-only or streaming-only models (even when allowing for super-alfvénic streaming speeds) too large γ-ray luminosities are obtained. However, the effects of CRs on SFR and gas density distributions that they obtained are relatively modest.

In a cosmological context, Jubelgas et al. (2008) studied the impact of CR acceleration on structure formation shocks and showed that CRs can significantly reduce the star-formation efficiency of small galaxies. Wadepuhl & Springel (2011) performed hydrodynamical simulations of the formation of Milky Way-sized galaxies and found that simulations that include CRs match the faint-end of the galaxy luminosity function. Salem et al. (2014) performed cosmological zoom-in simulations of a 1012 M⊙ halo and found that CRs generate thin, extended discs with a significantly more realistic rotation curve. Chen et al. (2016) performed cosmological zoom in simulations of dwarf galaxies; their simulations retrieve the observed baryonic Tully-Fisher relation and produce smoother SFRs than the bursty SFRs obtained with a thermal component only, but lead to too many stars and to the presence of DM cusps. Recently, Hopkins et al. (2020) presented the study of a large suite of high-resolution cosmological zoom-in simulations. They found that CRs have relatively weak effects on dwarf galaxies at redshifts above 1 − 2 for a diffusion coefficient, which matches observations of γ-ray emission. Overall, those numerical studies that have included CR transport and studied wind properties agree that CRs increase the mass loading and the density of galactic winds.

In the present work, we study the impact of CR feedback released by SNe on galaxies and wind formation using high resolution MHD simulations of dwarf galaxies of halo mass of 1010 M⊙ and 1011 M⊙, and varying the physics of CRs (isotropic or anisotropic diffusion, with different values of the diffusion coefficient, and testing for the effect of the streaming instability).

This paper is structured as follows. In Sect. 2, we describe the setup of our isolated galaxy disc simulations. In Sect. 3, we present the results of our analysis on the most massive dwarf galaxy, focusing on the suppression of star formation, the density distribution in the ISM, the generation of outflows, the distribution of CR energy and the evolution of the magnetic field. We also discuss how these results depend on the CR transport model. In Sect. 4 we discuss the dependency of CR-driven winds on halo mass. In Sect. 5, we examine how these results converge with numerical resolution. In Sect. 6 we test the dependency on the initial topology of the magnetic field. We discuss and summarise our results in Sect. 7.

2. Simulations and methods

Using the adaptive mesh refinement code RAMSES (Teyssier 2002), we ran idealised simulations of rotating isolated disc galaxies, consisting of gas and stars, embedded in a DM halo. The equations of ideal magneto-hydrodynamics (MHD; Fromang et al. 2006) are computed using the Harten-Lax-van Lear-Discontinuities Riemann solver from Miyoshi & Kusano (2005) and the MinMod slope limiter to construct variables at cell interfaces from the cell-centred values. The induction equation evolving the magnetic field is solved using constrained transport on the adaptive grid (Teyssier et al. 2006), which guarantees that the divergence of the magnetic field is zero at machine precision. Here we provide a short description of the simulation setup as well as the physical modules used, described in detail in Rosdahl et al. (2017). We used the method introduced by Dubois & Commerçon (2016) for solving the anisotropic diffusion of CR energy, that is, we used an implicit finite-volume method in the RAMSES code. We detail the CR magnetohydrodynamics in Sect. 2.7.

2.1. Initial conditions

The main parameters for the simulated galaxies and their host DM haloes are presented in Table 1. Similar simulation settings were used in Rosdahl et al. (2015, 2017). The initial conditions are generated with the MAKEDISC code (Springel et al. 2005). We focus most of our analysis on the higher mass galaxy that we name G9. The rotating isolated galaxy of G9 has a baryonic mass of Mbar = Mdisc + Mbulge = 3.8 × 109 M⊙, with an initial gas fraction of fgas = 0.5, and is hosted by a DM halo of mass Mhalo = 1011 M⊙. We also compare, at the same spatial resolution, a less detailed set of results for feedback models in a ten times less massive galaxy, that we name G8, with a baryonic mass Mbar = 3.8 × 108 M⊙, fgas = 0.5 and Mhalo = 1010 M⊙. Each simulation is run for 250 Myr that is 2.5 orbital times at the scale radii of G8 and G9.

Table 1.

Simulation initial conditions and parameters for the disc galaxy modelled in this paper.

The DM halo follows the NFW density profile of Navarro et al. (1997), with a concentration parameter c = 10 and spin parameter λ = 0.04 (Bullock et al. 2001). The DM is modelled by 106 collisionless particles, hence the galaxy halo has a DM mass resolution of 105 M⊙ for G9 and 104 M⊙ for G8. The initial disc consists of stars and gas, both set up with density profiles that are exponential in radius and Gaussian in height from the mid-plane (scale radius of 0.7 kpc for the G8 galaxy and 1.5 kpc for the G9 galaxy; scale height of one-tenth of the scale radius). The G8 and G9 galaxies contain a stellar bulge with masses and scale radii both one-tenth that of the disc. The initial stellar particle number is 1.1 × 106, a million of which are in the disc and the remainder in the bulge. The mass of the initial stellar particles is 1.7 × 102 M⊙ for G8 and 1.7 × 103 M⊙ for G9 (close to the masses of stellar particles formed during the simulation for G9, but ten times smaller for G8). While contributing to the dynamical evolution and gravitational potential of the rotating galaxy disc, the initial stellar particles do not explode as SNe. The temperature of the gas discs is initialised to a uniform T = 104 K and the ISM metallicity Zdisc is set to 0.1 for the G8 and G9 galaxies. The CGM initially consists of a homogeneous hot and diffuse gas, with nH = 10−6 cm−3, T = 106 K and zero metallicity. The cutoff radii and heights for the gas discs are chosen to minimise the density contrast between the disc edges and the CGM. The square box widths for the G8 and G9 galaxies are 150 kpc and 300 kpc, respectively, and we use outflow (i.e. zero gradient) boundary conditions on all sides.

This particular choice of idealised setup for the galaxy and its CGM is a great simplification of the cosmological infall composed of anisotropic cold narrow streams of matter and a hot and diffuse corona, with respective mass fraction depending on mass and redshift (see e.g. Dekel et al. 2009). It is feasible, alternatively, to study an idealised problem with an initially spherically-symmetric halo, where the collapse of the gaseous halo interferes with the propagation of winds (e.g. Jacob et al. 2018), but where such idealised setups are also a very simplified representation. The advantage of such a simplified setup is to allow for the study of the launching of galactic winds without the complex interplay with an inhomogeneous cosmological infall, whose aspect is beyond the scope of this paper.

2.2. Magnetic field

The initial magnetic field vector B is obtained by defining an initial value of the magnetic potential vector, so that the reconstruction of B = ∇ × A ensures ∇ ⋅ B = 0 (Dubois & Teyssier 2010). A toroidal B field is obtained by setting the value of A to

A = 3 2 B 0 r 0 ( ρ ρ 0 ) 2 / 3 e z , $$ \begin{aligned} \boldsymbol{A}=\frac{3}{2}B_{0}r_{0}\left(\frac{\rho }{\rho _0}\right)^{2/3}\boldsymbol{e}_{z}\, , \end{aligned} $$(1)

where ez is the z-Cartesian unit vector, and ρ is the gas density profile as defined in the disc within the vertical and radial cutoffs, ρ0 its normalisation (∼15 cm−3 for both G8 and G9) and r0 its scale radius (1.5 kpc for G9 and 0.7 kpc for G8). B0 is set to 1 μG. With this parametrization, a toroidal magnetic field with only non-zero initial components Bx and By is produced with a density-dependency ρ2/3 that mimics the effect of pure compression of magnetic field lines in the disc.

In Sect. 6, we test the influence of the initial topology of the magnetic field by changing the initial topology from toroidal to poloidal.

2.3. Adaptive refinement

Each refinement level uses half the cell width of the next coarser level, starting at the box width at the first level. Our simulations start at level 7 in a box of size Lbox = 300 kpc for G9 and Lbox = kpc for G8, corresponding to a coarse resolution of 1283 cells, and adaptively refine up to a maximum level 15 for G9 and level 14 for G8, corresponding to a minimum cell width of Δxmin = 9 pc. A cell is refined (derefined) up to the maximum level (respectively down to level 7) if its total mass exceeds 8 times 103 M⊙ (respectively is below 103 M⊙), or if its width exceeds a quarter of the local Jeans length. In Sect. 5, we perform a convergence study for the G9 galaxy, where the finest level is set to 14.

2.4. Gas cooling

Gas can radiatively cool down to a minimum temperature of 103 K (but further down due to adiabatic expansion) following Sutherland & Dopita (1993) for metal-dependant cooling rates above 104 K, and Rosen & Bregman (1995) below. We started with an initial metallicity of 0.1 Z⊙ that we kept constant (no release of metals by stellar evolution). We included heating of diffuse gas from a redshift zero Haardt & Madau (1996) ultraviolet (UV) background, and enforce an exponential damping of the UV radiation above the self-shielding density of nH = 10−2 cm−3.

2.5. Star formation

Our star formation model follows a Kennicutt-Schmidt law, ρ ˙ ∗ = ϵ ∗ ρ / t ff $ \dot{\rho}_*= \epsilon_*\rho/t_{\mathrm{ff}} $ (Kennicutt 1998; Krumholz & Tan 2007) where ρ ˙ ∗ $ \dot{\rho}_* $ is the SFR density, ϵ* = 0.02 the constant star-formation efficiency, ρ the gas density and tff is the free fall time of the gas. Star particles of mass m* = 2 × 103 M⊙ are formed stochastically using a Poissonian distribution (see Rasera & Teyssier 2006 for details), in cells that exceed the star-formation number density threshold n0 = 80 cm−3 and with a temperature below 3 × 103 K.

To prevent numerical fragmentation of gas below the Jeans scale (Truelove et al. 1997), an artificial “Jeans pressure” is maintained in each gas cell in addition to the thermal pressure. In terms of an effective temperature, it can be written as TJ = T0nH/n0, where we have set T0 = 103 K (and n0 is the aforementioned star-formation threshold), to ensure that the Jeans length is resolved by a minimum number of 7 cell widths at any density (see Eq. (3) in Rosdahl et al. 2015). The pressure floor is non-thermal, in the sense that the gas temperature that is evolved in the thermochemistry is the difference between the total temperature and the floor – therefore, we can have T ≪ TJ.

2.6. Supernova feedback

SN feedback is performed with single and instantaneous injections per stellar population particle tSN = 5 Myr after it has formed. Each stellar particle has an energy (ESN) and mass (mej) injection budget of

E SN = 10 51 η SN m ∗ m SN erg , $$ \begin{aligned}&E_{\rm SN} = 10^{51} \eta _{\rm SN} \frac{m_{*}}{m_{\mathrm{SN}}}\,\mathrm{erg}\, ,\end{aligned} $$(2)

m ej = η SN m ∗ , $$ \begin{aligned}&m_{\rm ej} = \eta _{\rm SN} m_{*}\,, \end{aligned} $$(3)

respectively. Here, ηSN is the fraction of stellar mass that is recycled into SN ejecta; mSN is the average stellar mass of a type II SN progenitor, and m* is the mass of a stellar particle. Assuming a Chabrier (2003) initial mass function, we set ηSN = 0.2, and mSN = 20 M⊙.

At each SN explosion, 0.1ESN is released in the form of CR energy in the cell hosting the SNR (see Sect. 2.7 for a description of CR physics). For the energy associated to the thermal component (0.9ESN), we used the mechanical feedback prescription detailed in Kimm & Cen (2014) and Kimm et al. (2015). In this method, the momentum is deposited into the neighbour cells of a SN hosting cell, with the magnitude adaptively depending on whether the adiabatic phase of the SN remnant is captured by this bunch of cells and the mass within them, or whether the momentum-conserving (snow-plough) phase is expected to have already developed on this scale. If the energy-conserving phase is resolved, the momentum is given by energy conservation, while if it is unresolved, the final momentum, which depends via the cooling efficiency on the density and metallicity, is given by Blondin et al. (1998), Thornton et al. (1998). We injected CR energy into the cell hosting the SN and neighbouring cells using the same relative weights as for the mass and momentum, to determine the CR energy assigned to each neighbour.

We note that CR injection could affect the SNR evolution itself and change the subgrid prescription for mechanical feedback: as a relativistic fluid (γCR = 4/3), CRs suffer less adiabatic loss than the thermal gas, so that at late times they dominate thermal pressure. Most importantly, CR energy has longer cooling time scales than the thermal energy, and is not radiated away during the snow-plough phase, but rather continues to support SNR expansion (Diesing & Caprioli 2018). In future work, we will investigate how the momentum deposition of SNR is affected by the presence of CRs, but that is beyond the scope of this paper.

2.7. Cosmic ray magnetohydrodynamics

CRs are modelled as a fluid represented by one separate energy equation with one single bin of energy density eCR which contributes to the total energy of the fluid in the equations of ideal MHD:

∂ρ dt +∇⋅(ρ u)=0, $$ \begin{aligned}&\frac{\partial \rho }{\mathrm{d}t} + \nabla \cdot (\rho \boldsymbol{u}) = 0, \end{aligned} $$(4)

∂ρu dt +∇⋅( ρuu+ P tot + BB 4π )=ρg+ p ˙ SN , $$ \begin{aligned}&\frac{\partial \rho \boldsymbol{u}}{\mathrm{d}t} +\nabla \cdot \left( \rho \boldsymbol{u}\boldsymbol{u}+P_{\rm tot} + \frac{\boldsymbol{B} \boldsymbol{B}}{4 \pi }\right)=\rho \boldsymbol{g} +\dot{p}_{\rm SN}, \end{aligned} $$(5)

∂ e d t + ∇ · ( ( e + P tot ) u − B ( B · u ) 4 π ) + ∇ · F st = ρ u · g + Q CR + Q th − Λ rad − Λ CR − ∇ · ( − κ b ( b · ∇ ) e CR ) , $$ \begin{aligned}&\frac{\partial e}{\mathrm{d}t} +\nabla \cdot \left( (e+P_{\rm tot}) \boldsymbol{u} - \frac{\boldsymbol{B} (\boldsymbol{B}\cdot \boldsymbol{u})}{4 \pi }\right) + \nabla \cdot \boldsymbol{F}_{\rm st} = \rho \boldsymbol{u} \cdot \boldsymbol{g} + Q_{\rm CR} &\\ \quad\quad& + Q_{\rm th} - \Lambda _{\rm rad} - \Lambda _{\rm CR} -\nabla \cdot (-\kappa \boldsymbol{b}(\boldsymbol{b}\cdot \nabla ) e_{\rm CR}),\end{aligned} $$(6)

∂B dt −∇×(u×B)=0, $$ \begin{aligned}&\frac{\partial \boldsymbol{B}}{\mathrm{d}t} - \nabla \times (\boldsymbol{u}\times \boldsymbol{B}) = \boldsymbol{0}, \end{aligned} $$(7)

∂ e CR d t + ∇ · ( e CR u ) + ∇ · F st = − P CR ∇ · u + Q CR − Λ CR − ∇ · ( − κ b ( b · ∇ ) e CR ) − u st · ∇ P CR , $$ \begin{aligned}&\frac{\partial e_{\rm CR}}{\mathrm{d}t} +\nabla \cdot ( e_{\rm CR} \boldsymbol{u}) + \nabla \cdot \boldsymbol{F}_{\rm st} = -P_{\rm CR}\nabla \cdot \boldsymbol{u} + Q_{\rm CR} - \Lambda _{\rm CR} \nonumber \\& -\nabla \cdot (-\kappa \boldsymbol{b}(\boldsymbol{b}\cdot \nabla ) e_{\rm CR}) -\boldsymbol{u}_{\mathrm{st}} \cdot \nabla P_{\rm CR}, \end{aligned} $$(8)

where ρ is the density, u is the gas vertical velocity, and B the magnetic field. The total energy density includes the kinetic, internal, CR and magnetic energy densities, that is e = 0.5ρv2 + eth + eCR + B2/8π, where eth is the internal energy; the total pressure includes the magnetic, thermal (Pth) and CR (PCR) pressures, that is Ptot = Pth + PCR + B2/8π, where PCR and Pth are related to their energy density with eth = Pth/(γ − 1) and eCR = PCR/(γCR − 1), with γ = 5/3 and γCR = 4/3 for a purely monoatomic ideal thermal component and a fully relativistic population of CRs; g is the gravitational field; ṗSN is the source therm for momentum injection by SNe; Qth is the thermal energy source term and includes UV background heating as well as CR collisional heating due to Coulomb collisions (Guo & Oh 2008); QCR is the energy source term for CRs from SNe explosions; Λrad is the radiative cooling of the thermal component; ΛCR is the total CR energy loss rate due to Coulomb and hadronic collisions (Guo & Oh 2008): ΛCR = 7.51 × 10−16(ne/cm−3)(eCR/erg cm−3) erg s−1 cm−3, with a corresponding reinjection to the thermal component of 2.63 × 10−16(ne/cm−3)(eCR/erg cm−3) erg s−1 cm−3 included in Qth. CR anisotropic diffusion term −∇⋅(−κb(b.∇)eCR) (where b is the magnetic field unit vector) is solved with the implicit solver of Dubois & Commerçon (2016). The value of the constant diffusion coefficient κ is varied within the range [3 × 1027, 1029] cm2 s−1 (see Table 2). Typical inferred values of the (isotropic) diffusion coefficient in the Milky Way from observational constraint obtained with GALPROP lie between 3 − 8 × 1028 cm2 s−1 for particles at a few GeV (Strong et al. 2007; Trotta et al. 2011). Hydrodynamical simulations with a fluid treatment of the CR physics produce consistent scaling of the gamma-ray luminosity from CR hadronic losses with SFR (Ackermann et al. 2012) or gas surface density (Lacki et al. 2011) using a broad range of diffusion coefficient 1028 − 3 × 1029 cm2 s−1 (Salem et al. 2016; Pfrommer et al. 2017a; Chan et al. 2019).

Table 2.

Models of CR injection and transport.

As CRs interact with Alfvén waves, they stream along their own gradient of pressure. This streaming of CRs is represented by the advection term Fst = fbust(eCR + PCR) and the associated cooling (for CRs) and heating (for the thermal component) term −ust.∇PCR, where ust = sign(b.∇eCR)uA is the CR streaming velocity ( u A = B / 4 π ρ $ \boldsymbol{u}_{\mathrm{A}}=\boldsymbol{B}/\sqrt{4\pi\rho} $ is the Alfvén speed, thus, the streaming is also fully anisotropic along the magnetic field lines). fb is the boost factor of the streaming velocity and accounts for the fact that Alfvén waves might experience damping by various mechanisms, allowing for weaker confinement by streaming, and, hence, a faster streaming velocity.

The ISM includes a cold, dense phase, where the effect of ion-neutral damping should be accounted for. The imperfect ionisation of the thermal particles causes CRs to couple less efficiently to the gas in this cold phase and allows them to escape these regions at velocities faster than the Alfvén velocity, while at the same time limiting the local diffusion coefficient (Wiener et al. 2013; Nava et al. 2016; Brahimi et al. 2019). As CRs enter the completely ionised warm and hot phases of the ISM, they establish a tight coupling while the weaker non-linear Landau damping dominates the wave damping (Wiener et al. 2013, 2017; Thomas & Pfrommer 2019). Turbulent damping of waves, effective in the various phases of the ISM, can control both the amount of coupling of CRs with the plasma and the level of anisotropy of the diffusive flux depending on the nature of the magnetised turbulence (Yan & Lazarian 2002; Farmer & Goldreich 2004; Lazarian 2016; Holguin et al. 2019; Nava et al. 2019). Finally, the shape of the spectrum of CRs, harder (flatter) close to acceleration sites (i.e. SN remnants), is expected to play an important role in how well CRs are tight to the plasma by the spectrum-dependent excited waves. In order to simplify the complex mechanisms of wave damping that change with the different ISM phases, in this work: (i) we have assumed diffusion with an effective galactic-wide constant coefficient (and which degree of anisotropy does not vary either), (ii) we have assumed streaming with a uniform value of velocity, and (iii) we have tested the effect of diffusion and streaming separately (while in practice their are intertwined and spatially varying). In this “effective” approach, the streaming velocity entering the advection term is boosted by a uniform factor of fb = 4 (not the heating term of streaming), which is a representative value for CR-excited waves that are damped mostly by turbulence in the diffuse ISM (see Ruszkowski et al. 2017, for details). We keep a more in-depth numerical investigation of wave-coupling effects for future work (Brahimi et al., in prep.).

While the heating term of streaming is straightforward to implement with a finite-differentiation, the advection term of streaming, if solved with a standard explicit Godunov-like approach puts stringent constraints on the time step Δt ∝ Δx3 due to the non-continuous nature of the streaming velocity at energy extrema (see Sharma et al. 2010, for discussion). A manipulation of the different terms in the streaming flux allows one to express this term equivalently as a pure diffusion term, which can then be added to the constant diffusion term κ, and to use the implicit solver for anisotropic diffusion and streaming (see Dubois et al. 2019, for details and tests of the method).

For our fiducial galaxy G9, we ran nine different models: a run without CRs (“noCR”); a run with CR injection but no transport (“Advection”); a run with isotropic diffusion of CRs (“Isodiff”); 4 runs with anisotropic diffusion where we varied the value of the diffusion coefficient κ parallel to the magnetic field; two runs with streaming where we used a boost factor of either fboost = 1 (“Streaming”) or1 fboost = 4 (“Streaming boost”). We present in Table 2 the different CR models presented in this paper. The Isodiff model has been run with a diffusion coefficient of κ = 3 × 1027 cm2 s−1: since in that model, CR diffusion is identical in all directions, the amount of CR diffusion is expected to be theoretically equivalent to the anisotropic diffusion run with κ = 1 × 1028 cm2 s−1.

The original implementation of the implicit diffusion solver introduced in Dubois & Commerçon (2016) has been here extended to include a minmod slope limiter to the transverse gradient of the CR energy density for the centred asymmetric scheme as proposed by Sharma & Hammett (2007) to preserve the monotonicity of the flux. Since the limiter produces a non-symmetric matrix, the linear system is solved with a biconjugate gradient stabilised method. This procedure guarantees that the values of the CR energy density obtained by the anisotropic diffusion/streaming implicit solver are always positives1.

3. Case study of the G9 galaxy

We begin by comparing the different models of CR injection and transport on the wind and ISM morphology, star formation, outflow generation, and CR energy distribution of the G9 galaxy.

3.1. Galaxy and wind morphologies

Figures 1 and 2 show the density, the velocity, the sum of CR and thermal pressures and ratio of CR to thermal pressures in 15 kpc-wide slices of the G9 galaxy viewed edge-on. Without CR injection (noCR), a low density wind (≤10−3 cm−3) is generated with velocities of a few hundred km s−1. With the injection of CR energy without any transport (Advection), the disc is puffed up but the wind is much weaker. CR streaming allows for a mild wind formation, weaker than in the simulations with CR diffusion. When adding isotropic diffusion of CRs, the wind is 10–100 times denser with similar velocities to the noCR case. The sum of CR and thermal pressures at 1–5 kpc from the mid-plane increases by 2 orders of magnitude and mainly consists of CR pressure. This strong effect on wind density and pressure is slightly mitigated but still substantial when modelling anisotropic diffusion. As already pointed out, the amount of CR diffusion in the Isodiff run is roughly equal to that of the anisotropic run with κ = 1 × 1028 cm2 s−1, but one sees that the wind is weaker in the κ = 1 × 1028 cm2 s−1 compared to the Isodiff run, which suggests that the azimuthal component of the magnetic field is stronger than its vertical component. We confirm this later in Sect. 3.5. The velocity of the winds increases with the anisotropic diffusion coefficient, and the density in the first 5 kpc away from the disc plane decreases at higher diffusion coefficients. The Streaming run do not differ much from the Advection case, only when streaming is boosted by a factor of 4, do we see a larger wind density and velocity.

thumbnail Fig. 1.

Slices of the G9 galaxy, from left to right: gas density, gas vertical velocity, sum of thermal and CR pressures, and ratio of CR to thermal pressure, seen edge-on at 250 Myr for the different simulations as indicated in the panels (one model per row). See Fig. 2 for the four remaining runs. In the noCR case, a low density wind is generated with velocities of a few hundred km s−1; in the Advection model, the disc is puffed up but the wind is much weaker. When adding CR diffusion, the wind is 10 times denser with velocities similar to noCR case. The sum of CR and thermal pressures, PCR + Pth, increases by 2 orders of magnitude and mainly consists of CR pressure. With CR streaming, the morphology is similar to the Advection case, but with a mild wind.

thumbnail Fig. 2.

Slices of the G9 galaxy, from left to right: gas density, gas velocity, sum of CR and thermal pressures, and ratio of CR to thermal pressure, seen edge-on at 250 Myr for the different simulations as indicated in the panels (one model per row). See Fig. 1 for the other runs. In the noCR case, a low density wind is generated with velocities of a few hundred km s−1; in the Advection model, the disc is puffed up but the wind is much weaker. When adding CR diffusion, the wind is 10 times denser with velocities similar to noCR case. The sum of CR and thermal pressures increases by 2 orders of magnitude and mainly consists of CR pressure. With CR streaming, the morphology is similar to the Advection case, but with a mild wind.

Figures 3 and 4 show the differences between the different models for the G9 galaxy viewed face-on. Unlike in the edge-on view, the density distributions in the ISM are more similar visually, but smoother with CR injection, as opposed to the presence of large pockets of low density gas in the noCR simulation. The injection of CR smooths the density distribution (as it will be seen later more quantitatively). This effect is mitigated at higher diffusion coefficient, because CR energy diffuses out more quickly. The sum of CR and thermal pressures in the disc increases by 1–2 orders of magnitude when adding CRs, and is either dominated by CR pressure by large in the diffuse parts of the ISM – and increasingly with the increase of the value of the diffusion coefficient –, or equally distributed between the CR and the thermal pressure in the cold filamentary ISM.

thumbnail Fig. 3.

Slices of the G9 galaxy, from left to right: gas density, slices of gas density zoomed-in compared to the left hand panels, sum of CR and thermal pressures, and ratio of CR to thermal pressure, seen edge-on at 250 Myr for the different simulations as indicated in the panels (one model per row). See Fig. 4 for the four remaining runs. The density distributions in the disc is smoother with CR injection. The sum of CR and thermal pressures in the disc increases by 1–2 orders of magnitude when adding CRs, and is either dominated by CR pressure, or equally distributed between CR and thermal pressure.

thumbnail Fig. 4.

Slices of the G9 galaxy, from left to right: gas density, slices of gas density zoomed-in compared to the left hand panels, sum of CR and thermal pressures, and ratio of CR to thermal pressure, seen edge-on at 250 Myr for the different simulations as indicated in the panels (one model per row). See Fig. 3 for the other runs. The density distributions in the disc is smoother with CR injection. The sum of CR and thermal pressures in the disc increases by 1–2 orders of magnitude when adding CRs, and is either dominated by CR pressure, or equally distributed between CR and thermal pressure.

3.2. Density distribution and star formation

Figure 5 shows the SFR and stellar mass of the G9 galaxy as a function of time. Adding CR injection from SNe reduces the SFR by a factor of 2–3 after 250 Myr. The simulation without any CR transport (Advection) is the most efficient at reducing the SFR because the CR pressure is only advected with the gas and therefore remains in the star-forming regions where it smooths out the highest densities (see Fig. 6) and regulates gas infall in dense clouds. In accordance with recent work from Chan et al. (2019) and Hopkins et al. (2020), we find that overall, CR injection has a modest effect on star formation. The simulations with the higher diffusion coefficient is where the SFR is the least affected because, as already pointed out, CR energy escapes from the star-forming regions and the pressure support against gravitational collapse in the ISM is weaker than for lower diffusion coefficients (see Sect. 3.4).

thumbnail Fig. 5.

Top panel: SFR as a function of time for the G9 galaxy with and without CR injection, with different CR transport models. Bottom panel: stellar mass as a function of time. Without CR injection (grey thick line), the SFR is 2–3 times greater than with CR injection. The simulation with advection of the CRs is where the SFR is the most affected because CR energy is trapped in the regions of star formation and does not escape.

Figure 6 shows the instantaneous mass-weighted probability density function (PDF) of the density in the ISM (typical wind densities are below n < 10−2 H cm−3 as it will be shown in Sect. 3.3) within a disc of 10 kpc radius and 4 kpc height at t = 100 and t = 250 Myr, for the nine different models2. The injection of CR energy systematically reduces the fraction of gas above the density threshold for star formation, which explains the decrease of the SFR observed in Fig. 5. The Advection model is the one where that fraction is the most suppressed, the reason being that the CR-pressure support in the star-forming regions is higher because the injected energy does not escape: we will see in Sect. 3.4 that the confinement of CRs in higher density regions also leads to more CR-energy loss since CR cooling is higher at higher densities. At early times (t = 100 Myr), the larger the diffusion coefficient is, the closer is the density PDF – and, hence, the SFR – to the run without CRs, which is reminiscent of a more uniform CR pressure in the disc, and Streaming runs are close to the result without any transport. Though, here, CR radiative losses and energy-injection by SNe complexifies the whole picture, the qualitative behaviour of the density PDF is consistent with simulations of turbulence in the cold neutral medium with CRs from Commerçon et al. (2019): when CR diffusion is large enough, CR pressure cannot be trapped in high gas density regions and the density PDF becomes similar to that without CRs. At late times (t = 250 Myr), the picture is more complex between the different diffusion coefficients but the qualitative behaviour as opposed to the Advection run is similar to that at early times, and the Streaming runs are now comparable to the runs with anisotropic diffusion. Conversely, CR injection increases the amount of gas at lower densities, and this effect is larger with diffusion, and CR streaming is the transport model for which this effect is the weakest at early times (t = 100 Myr) and becomes comparable to models with large diffusion coefficient at late times (t = 250 Myr).

thumbnail Fig. 6.

Mass-weighted probability density function of the gas density within a disc of 4 kpc height and 10 kpc radius, centered on the G9 galaxy, for the different runs, at t = 100 Myr (top panel) and t = 250 Myr (bottom panel) throughout the simulations as indicated. The vertical black dashed line shows the density threshold for star formation. Because it adds an additional pressure support against gravity that cools less efficiently than thermal pressure, CR injection reduces the high-density tail of the probability density function and increases the low-density tail.

3.3. Outflows and mass loading

In Fig. 7, we compare the time-evolution of outflows from the G9 galaxy in the different models. We measure the galactic-wide outflow properties across planes parallel to the galaxy disc, at a distance of 2 kpc in the left-hand panels and further out at 10 kpc in the right-hand panels. We show the mass outflow rate, the mass loading factor (the mass outflow rate divided by the SFR), the mass-weighted average velocity and the mass-weighted average density. To measure the mass outflow rate, we select the cells that cross the plane across which we compute the outflow rate, and sum the outflow rates ṁcell, obtained in individual cells as: m ˙ cell = ρ cell u z , cell Δ x , cell 2 $ \dot{m}_{\mathrm{cell}}=\rho_{\mathrm{cell}} u_{z,\mathrm{cell}} \Delta_{x,\mathrm{cell}}^2 $, where ρcell, uz, cell and Δx, cell are the cell gas density, gas vertical velocity and size. We select only the gas that is outflowing: that has a vertical velocity that is positive above the plane of the disc and negative below the plane of the disc. All quantities are very similar in the first 50 Myr, but the outflow rate quickly drops without CR injection or with CR injection but without transport, in the noCR and Advection models. In the Advection model, the lack of CR transport prevents the generation of outflows: the amount of outflowing gas is weaker than in the absence of CR energy, perhaps because the SN bubbles have more difficulty breaking out from the disc that is puffed up by CR pressure. The outflow rate in the Streaming simulation is very similar to the Advection case due to the low values of the Alfvén velocity compared to the SN shock velocity (see discussion in Sect. 3.3). For the Streaming Boost simulation, the outflow rate is similar to the Advection case at early times, but it rises to higher values at later times, although the outflow rate is sill weaker than with the other transport models. The mass outflow rate is higher when allowing for CR diffusion: CR energy is allowed to escape high density regions and accelerate more diffuse gas, that is accelerated at higher velocities than dense gas. The isotropic diffusion model and anisotropic diffusion with κ = 3 × 1028 − 1 × 1029 cm2 s−1 are the most efficient at driving winds: more that 10 times more gas is driven out of the 2 kpc and 10 kpc planes compared to the noCR case. The mass-averaged outflow velocity is less affected by the injection of CRs: it is increased at most by a factor of 2 for the highest diffusion coefficient, and for the lowest diffusion coefficient, it is even lower than in the noCR model. In the Isodiff run, the amount of CR diffusion is expected to be equivalent to the anisotropic diffusion run with κ = 1028 cm2 s−1 for a randomly oriented magnetic field. However, the wind mass loading is stronger in the Isodiff run, indicating that field lines in the disc are, indeed, preferentially within the plane of the galaxy suppressing the effective vertical diffusion of CRs (see Sect. 3.5). The wind after 250 Myr is ∼5 times denser with CR injection and diffusion with high diffusion coefficients than without CR injection.

thumbnail Fig. 7.

Mass outflow rates of the outflowing gas, mass loading factors, mass-weighted average velocities and mass-weighted average density of the outflowing gas across planes at 2 (left column) and 10 kpc (right column) above the disc plane of the G9 galaxy. The outflow rate quickly drops without CR injection and transport, in the noCR and Advection models. In the Advection-only model, the lack of CR transport prevents the generation of sustained outflows. The outflow rate in the Streaming Boost simulation is very similar to the Advection case at early times and rises to higher values at later times, but the outflow rate is sill weaker than with the other transport models. The isotropic diffusion model and anisotropic diffusion with κ = 3 × 1028 − 1 × 1029 cm2 s−1 are more than 10 times more efficient at driving winds compared to the noCR case. The mass average outflow velocity is less affected by the injection of CRs: it is increased at most by a factor of 2 for the highest diffusion coefficient, and is even lower for κ = 3 × 1027 cm2 s−1.

Figure 8 shows edge-on slices of the ratio of the gradient of the thermal (panels a) and CR (panels b) pressures over the vertical gravitational force in the different runs. A positive (red) ratio corresponds to a pressure gradient that pushes outwards. In the absence of CR injection or transport (noCR and Advection), the pressure above and below the mid-plane is dominated by the thermal pressure gradient, which is maximal in shocked regions. When including CR transport, the gradient of the CR pressure is much higher than that of the thermal pressure, and the gradient of the thermal pressure is lower than without CR transport because less shocks form (see Appendix A): the reason is that above and below the plane, the pressure is dominated by CR pressure (as shown in Fig. 1), of which the distribution is very smooth thanks to the transport of CR energy. One sees that with anisotropic diffusion, unlike with isotropic diffusion, the gradient of the thermal pressure is slightly higher in the wind, which stems from the fact that with anisotropic diffusion, CR transport is suppressed in the directions perpendicular to the magnetic field: this slightly mitigates the smoothness of the total pressure, and therefore shocks are more frequent. The ratio of the thermal pressure gradient to gravity exceeds unity only within shocks whereas the CR pressure gradient in the wind becomes uniformly greater (2 to 10 times) than the gravitational pull for κ ≥ 3 × 1028 cm2 s−1, and increases with the increasing value of the diffusion coefficient. Our results agree with that of Girichidis et al. (2018) (see also Simpson et al. 2016), where it is found that, when including CR diffusion, above ∼1 kpc away from the mid-plane, CRs provide the dominant gas acceleration mechanism. Gradients of CR pressure in the wind are higher for larger diffusion coefficient as more CR energy is allowed to escape the disc without hadronic and Coulomb losses as will be discussed in Sect. 3.4. The Streaming boost case shows a moderate amount of CR gradient of pressure in the wind, lower than that of the lowest anisotropic run studied, κ = 3 × 1027 cm2 s−1, though, even in that case, the CR gradient of pressure overwhelms that of thermal pressure.

thumbnail Fig. 8.

Thermal pressure gradient (a, two upper rows), and CR pressure gradients (b, two bottom rows) relative to the gravitational vertical pull, in 20 kpc slices of the G9 galaxy viewed edge-on, for the different runs, as indicated on the panels. A positive (red) ratio corresponds to a pressure gradient that pushes outwards. A linear threshold of 1 is used for the symmetric logarithmic colourbar. In the absence of CR injection or transport (noCR and Advection), the pressure above and below the mid-plane is dominated by thermal pressure gradient that is maximal in shocked regions. When including CR transport, the gradient of the CR pressure is much higher than that of the thermal pressure, and the gradient of the thermal pressure is lower and less shocks form. The ratio of the thermal pressure gradient to gravity exceeds unity only within shocks whereas the CR pressure gradient in the wind becomes uniformly greater (2 to 10 times) than the gravitational pull for κ ≥ 3 × 1028 cm2 s−1.

Figure 9 shows the temperature-density phase diagrams in the outflowing gas after 250 Myr for a representative subset of simulations (noCR, Isodiff, κ = 3 × 1028 cm2 s−1 and Streaming boost). We use a lower threshold of 10 km s−1 to detect the outflowing gas, selected only above 2 kpc from the disc plane. Without CR injection, the densest (10−4 − 10−3 cm−3) component of the outflowing gas has temperatures of 105 − 105.5 K. When including CR injection with transport, that component is colder (104 K) and denser (10−4 − 10−2 cm−3), in agreement with recent work from Girichidis et al. (2018) who simulated detailed slabs of the ISM. Interestingly, the temperature of the wind in the isotropic case is much more uniform than when modelling anisotropic diffusion: again, this is due to the fact that diffusion is suppressed in the directions perpendicular to the magnetic field, which reduces the smoothness of the total pressure, and shocks are more frequent, which complicates the temperature distribution of the outflow compared to the isotropic case. With CR injection but no CR transport (Advection), the wind component is almost non-existent (not shown here), as previously described, and the wind component is very weak for CR streaming with boost but has the same temperatures as with CR diffusion.

thumbnail Fig. 9.

Total mass of the outflowing gas showed in bins of temperature-density, after 250 Myr for the noCR, Isodiff, κ = 3 × 1028 cm2 s−1 and Streaming boost runs as indicated in the corresponding panels. We use a lower threshold of 10 km s−1 to detect the outflowing gas, selected only above 2 kpc from the disc plane. Without CR injection, the densest (10−4 − 10−3 cm−3) component of the outflowing gas has temperatures of 105 − 105.5 K. When including CR injection with transport, that component is colder (104 K) and denser (10−4 − 10−2 cm−3). In the Advection case, the wind component is almost nonexistent. The wind component is very weak with CR streaming but with similar temperatures to simulations with CR diffusion. The upper left part of the phase diagram corresponds to gas of the intergalactic medium at the outer edges of the disc.

3.4. Fate of the cosmic ray energy

Figure 10 shows the total CR energy in the simulated box for the different runs. The amount of CR energy in the box at a given time increases with the diffusion coefficient, the reason being that it escapes more efficiently the high density star-forming regions where it is injected, and where the cooling rate is higher. The total energy in the Advection run is lower for two reasons: first because the SFR is lower and therefore the total CR injection is weaker, second because the injected CR energy stays confined in the higher density regions from where it can escape only by advection with the gas, and therefore it looses energy to CR cooling (cf. Fig. 11). In the Isodiff run, CR energy rises faster because it does not have to follow the magnetic field lines, thereby avoiding more rapid cooling in the higher density disc. Conversely, in the anisotropic diffusion, the energy injected in the disc stays trapped in it for longer times. At later times, however, the amount of CR energy in the anisotropic diffusion simulation with κ = 3 × 1028 − 1 × 1029 cm2 s−1 is higher than in the Isodiff run because the SFR in that simulation being higher, CR energy injection is higher too, which compensates for the higher cooling rates.

thumbnail Fig. 10.

Total CR energy in the box the different runs of the G9 galaxy. The higher the diffusion coefficient, the higher the CR energy: the reason is that the more CRs can diffuse, the more they can escape high gas density and CR-energy-density regions where the CR cooling rate is higher. At the beginning of the simulation, with isotropic diffusion, CR energy quickly escapes in directions perpendicular to the disc, whereas it remains trapped in the disc for anisotropic diffusion.

Figure 11 displays the sum of hadronic and Coulomb CR energy losses relative to the total energy in the box (upper panel), as well as adiabatic losses3 (bottom panel), that is CR energy exchanges with the gas via compression and expansion (the volumetric rate of CR adiabatic losses expresses as −PCR∇⋅u). In the bottom panel, we also show the losses due to Alfén wave heating in the Streaming boost simulations −ust ⋅ ∇PCR (recall that the loss term is not boosted). We find that Alfvén wave heating is subdominant and that adiabatic losses are mainly subdominant except in the Isodiff and anisotropic diffusion runs with high diffusion coefficient simulations. In accordance with previous work (Pfrommer et al. 2017a; Girichidis et al. 2018; Chan et al. 2019; Hopkins et al. 2020), we find that the hadronic and Coulomb cooling rate is higher for lower diffusion and the two streaming solutions. The reason is that CR energy spends more time in dense regions where the hadronic and Coulomb CR cooling rate is higher (it scales linearly with gas density), as shown in Fig. 12, which displays, at different times, the instantaneous PDF of the density in the disc, weighted by CR energy: it shows the densities at which bulk of CR energy lies. One sees that in the Advection model, CR energy is concentrated in regions with densities that are higher than when including diffusion, and CR diffusion allows CR energy to populate lower density regions, where the hadronic and Coulomb cooling time is longer.

thumbnail Fig. 11.

Top panel: hadronic and Coulomb cooling rate of CR energy inside the box for the different runs. Bottom panel: CR adiabatic losses (−PCR∇⋅u <  0). We find that the cooling rate is higher for lower diffusion. We find that adiabatic losses are subdominant except in the isotropic and anisotropic diffusion simulations with the highest diffusion coefficients. The black dotted line corresponds to the streaming losses (−ust ⋅ ∇PCR) in the “Streaming Boost” simulation.

thumbnail Fig. 12.

CR energy-weighted density PDF a disc centred on the G9 galaxy (4 kpc height and 10 kpc radius), at 250 Myr. The upper axis shows the corresponding hadronic and coulomb cooling time. The vertical black dashed line is the density threshold for star formation.

3.5. Magnetic field evolution

We show in Fig. 13 the evolution of the vertical component – perpendicular to the disc – of the magnetic field for the different runs, in the disc (panels a,b) and in the wind (panels c,d), and the norm of the magnetic field vector in the wind (panel e). This is particularly relevant to CR transport since CR streaming and anisotropic diffusion are along magnetic field lines, and additionally the efficiency of streaming scales with the strength of the magnetic field. If the vertical component is strong, CR energy can escape from the disc to more diffuse and higher altitude regions. Conversely, if the magnetic field lines are parallel to the disc, CR energy is confined at lower altitudes close to the disc. In the disc (panel a), the amplification of the vertical component of the magnetic field is very similar in the different runs, although the noCR runs stands out with a slightly higher amplification. As shown in panel b, the vertical component is on average 2–3 times lower than the value for a randomly oriented field – where one would have 1 / 3 × | B | $ 1/\sqrt{3} \times |\boldsymbol{B}| $ –, which shows that in the disc, the magnetic field lines are dominated by their azimuthal component. Therefore, when the transport is anisotropic, CRs are more confined in the disc because they follow the magnetic field lines, as opposed to isotropic diffusion. This explains the fact that isotropic diffusion with κ = 3 × 1027 cm2 s−1 is more efficient at generating winds than the anisotropic run with κ = 1 × 1028 cm2 s−1 since field lines are more likely to be within the plane of the disc than along the direction of propagation of the galactic wind.

thumbnail Fig. 13.

Evolution of the (a) mass averaged magnetic field B in the disc (of 4 kpc height and 10 kpc radius); (b) mass averaged ratio of Bz over B in the disc (c) mass averaged Bz in the wind (more than 2 kpc away from the disc) (d) mass averaged ratio of Bz over B in the disc wind; mass average norm of B in the wind. The mean magnetic field strength in the disc saturates at ⟨B⟩≃5 μG with a similar vertical structure in all of the simulated cases. In the wind, a clear difference appears when adding diffusion (isotropic or anisotropic), and the vertical component is more than ten times greater with CR injection and diffusion than without CR injection.

In the wind (panel c), a clear difference appears when adding diffusion (isotropic or anisotropic): the vertical component is more than ten times greater with CR injection and diffusion than without CR injection. As shown in panel d, this difference is partly due to the fact that the vertical component is greater (relative to the total amplitude of the magnetic field) in the runs with diffusion, but the main reason is that the total amplitude of the magnetic field is greater in the wind, as shown in panel e. One sees that in the anisotropic diffusion runs, the vertical component is roughly 3 to 4 times smaller than the norm of the magnetic field, which means that the field in the wind is also dominated by its azimuthal component: this also explains the fact that isotropic diffusion is more efficient at generating winds than the anisotropic run with κ = 1028 cm2 s−1 because more CR energy is transported vertically. Nonetheless, the magnetic field lines are more vertical in the wind than they are in the disc indicating that the launching of the wind is able to open up the strongly azimuthal field lines from the disc. The higher values of the magnetic field in the wind in the runs with diffusion come from the largest densities in the wind, as magnetic field scales with density (not shown here, but the scaling of the magnetic field shows a linear dependency with density with a large scatter) in agreement with Pakmor et al. (2016).

The magnetic field in the disc quickly increases – as a result of both differential rotation (Dubois & Teyssier 2010), turbulent dynamo (Rieder & Teyssier 2016) or in combination with CRs (Hanasz et al. 2004) – in the first 100 Myr until it saturates at a value of ⟨B⟩mw, disc between 3–6 μG depending on the simulation. Those final values of the saturated magnetic field are in agreement with observations (see e.g. Beck 2000), and can reach values as high as 10–20 μG in the dense star-forming gas, consistent with the magnetic field observed in cold clouds (Crutcher 2012). The saturated values of the magnetic obtained in the simulated disc are broadly consistent with those obtained in Pakmor et al. (2016) or Pfrommer et al. (2017b) for a similar halo mass, despite the more resolved ISM structure, due to higher resolution, obtained here. This sheds light onto the delay for the Streaming run to operate: the characteristic wave speed of the streaming instability is the Alfvén velocity (here boosted by a factor of 4), and until the magnetic field reaches rough equipartition with the kinetic and thermal energy in the disc, the Alfvén velocity cannot compete with others characteristic velocities of the flow. It is possible to estimate crudely the expected effect of streaming in comparison to diffusion. Considering star-forming clouds of size of Lc ≃ 50 pc, the characteristic diffusion velocity is uD = κ/Lc, hence, for κ = 3 × 1027 cm2 s−1, we obtain uD = 195 km s−1. The boosted streaming velocity in comparison is of the order 4 × uA = 50 km s−1 for a magnetic field of 5 μG and a gas density of 10 cm−3 (in practice the boosted streaming velocity in the ISM varies from a few to 100 km s−1), well below the velocity of the smallest diffusion coefficient used here. Indeed, the streaming model in which the streaming velocity equals the Alfvén velocity produces galactic-wide properties more similar to that of the Advection case than they already are in the boosted case, that is it further weakens the mass loading of the galactic wind.

4. Galaxy/Halo mass dependency

We perform simulations of a subset of the models presented in Table 2 for the smaller G8 galaxy. We use the same spatial resolution of 9 pc, and the same resolution for the formed stellar particle mass, although the mass of the initial stellar particles (that do not explode as SNe) is ten times smaller than for the simulations of the G9 galaxy. The DM particle mass is also ten times smaller.

In Fig. 14, we show how the mass-loading (top panel), outflow rate (middle panel), and SFR (bottom panel) scale with the circular velocities of the G8/G9 galaxies and how these scalings compare with recent observations from Heckman et al. (2015) and Chisholm et al. (2017). We compute the circular keplerian velocity at the galactic scale radii of the G8 and G9 galaxies, and we average the outflow rate, mass loading factor and SFR between 75 and 250 Myr. Overall the SFR of the galaxies for any of the models are very similar and fall within the observed SFR. As a result of CR feedback, the SFR for the low-mass G8 galaxy are reduced by a factor ∼1.5 − 2 and the impact of CRs on the SFR is reduced with larger diffusion coefficients in a similar fashion to G9. Therefore, as in the G9 galaxy, the injection of CRs in the G8 galaxy has a mild impact on star formation.

thumbnail Fig. 14.

Markers show the time averaged (between 75 – 250 Myr) outflow rates (top) mass loading factors (middle) and SFR (bottom) as a function of the circular keplerian velocity of G8 (41 km s−1) and G9 (90 km s−1) at the galactic scale radii. For G9, some of the markers are shifted for clarity. The data points are taken from Heckman et al. (2015) and Chisholm et al. (2017). The labels with the suffix “-conv” are convergence study runs and are discussed in Sect. 5. The poloidal (blue diamond) run is discussed in Sect. 6.

These observational studies by Heckman et al. (2015) and Chisholm et al. (2017) use UV absorption lines to measure outflow velocities and determine mass loss rates in nearby star-forming galaxies. The inferred mass loading factors are rather uncertain. These quantities can be directly determined in numerical simulations whereas observationally, they rely on challenging measurements and estimates, and caution is therefore needed in the interpretation of this comparison. Namely, the radius of the outflow measurement is rather uncertain. We measure the outflow at 1kpc for the G8 galaxy close to the value of 2r* where r* is the half-light radius of the galaxy, following Heckman et al. (2015). The mass loading factor and outflow rates fall significantly below the observations, especially for the G8 galaxy, despite the significant improvement due to the presence of CR feedback and transport. Therefore, it suggests that even though the energy injection from the SFR is similar, the outflow generation is still too weak. We note that our results differ from Jacob et al. (2018) who find ten times higher mass loading factors for this mass range. This difference might stem from the differences in the initial set-up: Jacob et al. (2018) use an initial rotating gaseous halo in a static Hernquist potential, which might produce higher outflow rates – compared to our simulations where the gas is already settled in a disc – since the gas initially distributed in a sphere can be swept along with the outflows. However, Jacob et al. (2018) also find some disagreement with observations: they find that CR-driven winds drop rapidly with halo mass, much faster than in observations. In addition, one should be aware that there are some important uncertainties in the fraction of the CR energy released into SN explosions (or, equivalently, the CR acceleration efficiency at shocks), and, though, our fiducial value of CR energy fraction of 10 per cent is standard, larger values can be considered and would increase the strength of galactic winds (see Jacob et al. 2018). Another interesting aspect is the different hierarchy of the diffusion coefficient found for G8 and G9. For instance, in the G8 galaxy, the highest diffusion coefficient is the least efficient at launching winds, whereas it is amongst the most efficient ones in G9. Conversely, the lowest diffusion coefficient is the most efficient at driving winds in G8 and the least efficient in G9. This might stem from the fact that if CR energy escapes too quickly from the galaxy, it cannot help generating winds, and at fixed diffusion coefficient, CR energy escapes faster from a smaller size galaxy.

5. Convergence study

To study the convergence, we ran a subset of the simulations presented in Table 2 of the G9 galaxy with lower spatial resolution. Since the spatial resolution is twice lower, we reduced the density threshold for star formation from 80 cm−1 for the higher resolution fiducial runs to 10 cm−1 for the lower resolution runs. The masses of the formed stellar particles and the DM particle mass are the same as in our fiducial runs. We show the results in Fig. 14. The SFRs are slightly higher in the higher resolution runs, and therefore the mass loading factor are slightly higher in the lower resolution run. The values of the outflow rate agree very well for the two resolutions.

6. Influence of the initial magnetic field topology

We test the influence of the initial magnetic field topology by running one of the models (κ = 1 × 1028 cm2 s−1) with a poloidal topology rather than toroidal. A poloidal B field is obtained by setting the value of A to

A = B 0 ( ρ ρ 0 ) 2 / 3 ( − y x 0 ) , $$ \begin{aligned} \boldsymbol{A}=B_{0}\left(\frac{\rho }{\rho _0}\right)^{2/3} \begin{pmatrix} -{ y} \\ x \\ 0 \end{pmatrix} \,, \end{aligned} $$(9)

where x and y, are the Cartesian coordinates, and ρ is the gas density profile as defined in the disc within the vertical and radial cutoffs, ρ0 its normalisation (∼15 cm−3 for both G8 and G9) and r0 its scale radius (1.5 kpc for G9 and 0.7 kpc for G8). B0 is set to 1 μG. We compare the SFR and the outflows in the κ = 1 × 1028 cm2 s−1 simulations run with the two different magnetic field initialisations: toroidal and poloidal. We show the result for the SFR, outflow rate and mass loading in Fig. 14 (blue diamonds). The SFR in the case of poloidal initial topology is similar to the SFR in the case of toroidal initial topology. The outflow rate is slightly higher: the reason is that the vertical component is strong in the initial poloidal topology but nonexistent in the toroidal topology and therefore CR energy escapes more easily in the vertical directions and generates more winds. Figure 13 shows the evolution of the vertical component – perpendicular to the disc – of the magnetic field for the different runs, in the disc (panels a,b) and in the wind (panels c,d), and the norm of the magnetic field vector in the wind (panel e). In the disc, the vertical component evolution is very similar in all cases. In the wind, the vertical component of the magnetic field is greater with a poloidal initialisation (panel c), because the norm is greater (panel e), but relative to the norm of the magnetic field in the wind, it converges to the same value (panel d).

7. Summary and discussion

We performed high resolution simulations of isolated disc galaxies embedded in halos of 1010 M⊙ and 1011 M⊙, with and without CR injection, and with different CR transport models. We summarise below our main results:

  • The amount of star formation is reduced by a factor of ∼2, which is due to the suppression of the high density tail of the density PDF by CR pressure support. Overall, the effect of CR injection has a soft effect on star formation.

  • The injection of CRs significantly increases the wind generation only when diffusive transport is included. The isotropic diffusion model (κ = 3 × 1027 cm2 s−1) and anisotropic diffusion with κ = 3 × 1028 − 1 × 1029 cm2 s−1 are more than 10 times more efficient at driving winds than without CR injection. Outflows with CR diffusion are ten times colder and up to a hundred times denser than without CRs. The gradient of CR pressure provides the dominant gas acceleration mechanism above and below the plane of the disc.

  • The injection and diffusion of CRs increases the total pressure gradient, which exceeds the gravitational pull by a factor of 2 to 10 in the wind. Conversely, the ratio of the thermal pressure gradient to gravity does not exceed unity in the wind except within shocks.

  • The energy losses for CRs are higher for lower diffusion. The reason is that CR energy spends more time in dense star-forming regions where the CR cooling rate is higher. Adiabatic losses and gain are mostly subdominant, and are higher for higher diffusion coefficients.

  • The streaming of CRs only plays a minor role: the star formation, outflow rates and morphology are very similar to the case without transport at all. This suggests that CR streaming is not an efficient mechanism to generate winds. Our results on CR streaming differ from Ruszkowski et al. (2017) but agree with Chan et al. (2019) and Hopkins et al. (2020).

  • Although the injection and transport of CRs increases the mass outflow rate and brings our simulated mass loading factors closer to observations and in much better agreement than without CRs, the outflow rates in our simulations are still on the weak part of observed outflow rates.

Although this idealised set up of initial conditions allows for a good understanding of the physical impact of the modelled physics, such an isolated disc set-up – meaning that we do not have a dense pre-existing circum-galactic medium nor realistic anisotropic gas accretion – is a limitation to how well our simulations can capture the rich interactions of galactic winds with their environment. The next step would be to run cosmological zoom-in simulations of CR feedback in dwarf galaxies, that allow to simulate dwarf galaxies with enough resolution to capture the scales relevant to stellar feedback and enough volume to take into account environmental effects and gas inflow.

There are as additional forms of stellar feedback processes at play in galaxies that we have neglected here, such as the stellar radiation from the UV producing photo-ionised bubbles or from the IR interacting with the dusty gas (Hopkins et al. 2014; Rosdahl et al. 2015), that need to be combined with the SN release of CRs (see Chan et al. 2019). We have also simplified the process of star formation by using a constant star-formation efficiency, while recent models of simulated galaxies now scale the efficiency with turbulent properties of the star-forming gas (Semenov et al. 2016; Kimm et al. 2017; Trebitsch et al. 2018; Lupi et al. 2018) calibrated on the multi-free fall models of Federrath & Klessen (2012) (based on Krumholz & McKee 2005; Padoan & Nordlund 2011; Hennebelle & Chabrier 2011). Such coupling of varying star-formation efficiency with the release of CRs in SNe remains to be explored.

Another shortcoming of our model is that the adopted CR physics model is simplified since it considers one single CR momentum-energy bin, and, hence, uniform energy losses and a single diffusion coefficient across the entire spectrum of CR momenta. However, in reality the diffusion coefficient and the amount of energy losses through hadronic and Coulomb losses depend on the considered CR energy. Future improvements will include several CR energy bins to account for the CR spectrum, especially since it represents a crucial observable that one can use in order to constrain wind models (Recchia et al. 2017).

Finally, the injection of CRs at shocks, with CR acceleration efficiencies depending on the structure of the background magnetic field (Caprioli & Spitkovsky 2014), might change how CRs are released into the ISM (Pais et al. 2018; Dubois et al. 2019), and, therefore, how they propel the large-scale galactic winds. We defer this to future investigations.

1

All the simulations have also been run with the original unlimited version, which does not guarantee the positivity of the solution. Thus, in that case, the CR energy can become negative, and has to be readjusted. This translates into a spurious CR energy injection upon SN explosions, that depends on the diffusion coefficient: namely at each SN explosion, an additional 3%, 7%, 12%, 15% of CR energy (relative to the expected CR energy to be injected) is injected for diffusion coefficients of respectively 3 × 1027 cm2 s−1, 1 × 1028 cm2 s−1, 3 × 1028 cm2 s−1, 1 × 1029 cm2 s−1. The results in terms of, for instance, mass outflow rates or star-formation rates do not differ by more than ∼10% with the current positivity-preserving version of the solver.

2

Although it is a small correction, all PDF are renormalised to that of the noCR run, so that their integrated value leads to the same amount of total mass.

3

We found that the adiabatic gain was nonexistent or negligible depending on the run so that we decided not to show it.

Acknowledgments

We thank Corentin Cadiou for his tremendous help with YT and other fruitful scientific discussions. We thank Joakim Rosdahl for his help in setting up the initial conditions and valuable advice and comments on the paper. We thank A. Marcowith for enlightening discussions. This work was granted access to the HPC resources under the allocation A0060406955. This work has made use of the Horizon Cluster hosted by the Institut d’Astrophysique de Paris. We thank S. Rouberol for running smoothly this cluster for us.

References

  1. Ackermann, M., Ajello, M., Allafort, A., et al. 2012, ApJ, 755, 164 [NASA ADS] [CrossRef] [Google Scholar]
  2. Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25 [NASA ADS] [CrossRef] [Google Scholar]
  3. Aumer, M., White, S. D. M., Naab, T., & Scannapieco, C. 2013, MNRAS, 434, 3142 [NASA ADS] [CrossRef] [Google Scholar]
  4. Beck, R. 2000. Philos. Trans. R. Soc. London Ser. A, 358, 777 [Google Scholar]
  5. Berezinsky, V., Gazizov, A., & Grigorieva, S. 2006, Phys. Rev. D, 74, 043005 [NASA ADS] [CrossRef] [Google Scholar]
  6. Blondin, J. M., Wright, E. B., Borkowski, K. J., & Reynolds, S. P. 1998, ApJ, 500, 342 [NASA ADS] [CrossRef] [Google Scholar]
  7. Booth, C. M., Schaye, J., Delgado, J. D., & Dalla Vecchia, C. 2012, MNRAS, 420, 1053 [Google Scholar]
  8. Booth, C. M., Agertz, O., Kravtsov, A. V., & Gnedin, N. Y. 2013, ApJ, 777, L16 [NASA ADS] [CrossRef] [Google Scholar]
  9. Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40 [CrossRef] [Google Scholar]
  10. Brahimi, L., Marcowith, A., & Ptuskin, V. 2019, Int. Cosmic Ray Conf., 36, 40 [Google Scholar]
  11. Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, A&A, 245, 79 [NASA ADS] [Google Scholar]
  12. Brooks, A. M., & Zolotov, A. 2014, ApJ, 786, 87 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240 [NASA ADS] [CrossRef] [Google Scholar]
  14. Butsky, I. S., & Quinn, T. R. 2018, ApJ, 868, 108 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bykov, A. M. 2014, A&ARv., 22, 77 [NASA ADS] [CrossRef] [Google Scholar]
  16. Caprioli, D. 2015, Int. Cosmic Ray Conf., 34, 8 [Google Scholar]
  17. Caprioli, D., & Spitkovsky, A. 2014, ApJ, 783, 91 [NASA ADS] [CrossRef] [Google Scholar]
  18. Chabrier, G. 2003, PASP, 115, 763 [NASA ADS] [CrossRef] [Google Scholar]
  19. Chan, T. K., Kereš, D., Oñorbe, J., et al. 2015, MNRAS, 454, 2981 [NASA ADS] [CrossRef] [Google Scholar]
  20. Chan, T. K., Kereš, D., Hopkins, P. F., et al. 2019, MNRAS, 488, 3716 [Google Scholar]
  21. Chen, J., Bryan, G. L., & Salem, M. 2016, MNRAS, 460, 3335 [Google Scholar]
  22. Chisholm, J., Tremonti, C. A., Leitherer, C., & Chen, Y. 2017, MNRAS, 469, 4831 [NASA ADS] [CrossRef] [Google Scholar]
  23. Commerçon, B., Marcowith, A., & Dubois, Y. 2019, A&A, 622, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Crutcher, R. M. 2012, ARA&A, 50, 29 [NASA ADS] [CrossRef] [Google Scholar]
  25. Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, 426, 140 [NASA ADS] [CrossRef] [Google Scholar]
  26. de Jong, T., Klein, U., Wielebinski, R., & Wunderlich, E. 1985, A&A, 147, L6 [NASA ADS] [Google Scholar]
  27. Dekel, A., & Silk, J. 1986, ApJ, 303, 39 [NASA ADS] [CrossRef] [Google Scholar]
  28. Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  29. Diesing, R., & Caprioli, D. 2018, Phys. Rev. Lett., 121, 091101 [NASA ADS] [CrossRef] [Google Scholar]
  30. Dubois, Y., & Commerçon, B. 2016, A&A, 585, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Dubois, Y., & Teyssier, R. 2010, A&A, 523, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Dubois, Y., Commerçon, B., Marcowith, A. R., & Brahimi, L. 2019, A&A, 631, A121 [CrossRef] [EDP Sciences] [Google Scholar]
  34. Emerick, A., Bryan, G. L., & Mac Low, M.-M. 2018, ApJ, 865, L22 [Google Scholar]
  35. Farmer, A. J., & Goldreich, P. 2004, ApJ, 604, 671 [NASA ADS] [CrossRef] [Google Scholar]
  36. Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156 [NASA ADS] [CrossRef] [Google Scholar]
  37. Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Girichidis, P., Naab, T., Walch, S., et al. 2016, ApJ, 816, L19 [NASA ADS] [CrossRef] [Google Scholar]
  39. Girichidis, P., Naab, T., Hanasz, M., & Walch, S. 2018, MNRAS, 479, 3042 [NASA ADS] [CrossRef] [Google Scholar]
  40. Gronke, M., Girichidis, P., Naab, T., & Walch, S. 2018, ApJ, 862, L7 [CrossRef] [Google Scholar]
  41. Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251 [NASA ADS] [CrossRef] [Google Scholar]
  42. Haardt, F., & Madau, P. 1996, ApJ, 461, 20 [NASA ADS] [CrossRef] [Google Scholar]
  43. Hanasz, M., Kowal, G., Otmianowska-Mazur, K., & Lesch, H. 2004, ApJ, 605, L33 [NASA ADS] [CrossRef] [Google Scholar]
  44. Hanasz, M., Lesch, H., Naab, T., et al. 2013, ApJ, 777, L38 [NASA ADS] [CrossRef] [Google Scholar]
  45. Heckman, T. M., Alexandroff, R. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 809, 147 [NASA ADS] [CrossRef] [Google Scholar]
  46. Helder, E. A., Vink, J., Bamba, A., et al. 2013, MNRAS, 435, 910 [NASA ADS] [CrossRef] [Google Scholar]
  47. Hennebelle, P., & Chabrier, G. 2011, ApJ, 743, L29 [NASA ADS] [CrossRef] [Google Scholar]
  48. Holguin, F., Ruszkowski, M., Lazarian, A., Farber, R., & Yang, H. Y. K. 2019, MNRAS, 490, 1271 [NASA ADS] [CrossRef] [Google Scholar]
  49. Hopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, 417, 950 [NASA ADS] [CrossRef] [Google Scholar]
  50. Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hopkins, P. F., Chan, T. K., Garrison-Kimmel, S., et al. 2020, MNRAS, 492, 3465 [Google Scholar]
  52. Hu, C.-Y. 2019, MNRAS, 483, 3363 [Google Scholar]
  53. Ipavich, F. M. 1975, ApJ, 196, 107 [Google Scholar]
  54. Jacob, S., Pakmor, R., Simpson, C. M., Springel, V., & Pfrommer, C. 2018, MNRAS, 475, 570 [NASA ADS] [CrossRef] [Google Scholar]
  55. Jubelgas, M., Springel, V., Enßlin, T., & Pfrommer, C. 2008, A&A, 481, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Kennicutt, R. C., Jr 1998, ApJ, 498, 541 [NASA ADS] [CrossRef] [Google Scholar]
  57. Kimm, T., & Cen, R. 2014, ApJ, 788, 121 [NASA ADS] [CrossRef] [Google Scholar]
  58. Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015, MNRAS, 451, 2900 [NASA ADS] [CrossRef] [Google Scholar]
  59. Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826 [NASA ADS] [Google Scholar]
  60. Kimm, T., Haehnelt, M., Blaizot, J., et al. 2018, MNRAS, 475, 4617 [NASA ADS] [CrossRef] [Google Scholar]
  61. Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82 [Google Scholar]
  62. Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250 [NASA ADS] [CrossRef] [Google Scholar]
  63. Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304 [NASA ADS] [CrossRef] [Google Scholar]
  64. Kulsrud, R. M. 2005, Plasma Physics for Astrophysics [CrossRef] [Google Scholar]
  65. Lacki, B. C., Thompson, T. A., Quataert, E., Loeb, A., & Waxman, E. 2011, ApJ, 734, 107 [NASA ADS] [CrossRef] [Google Scholar]
  66. Lazarian, A. 2016, ApJ, 833, 131 [NASA ADS] [CrossRef] [Google Scholar]
  67. Lupi, A., Bovino, S., Capelo, P. R., Volonteri, M., & Silk, J. 2018, MNRAS, 474, 2884 [Google Scholar]
  68. Mao, S. A., & Ostriker, E. C. 2018, ApJ, 854, 89 [NASA ADS] [CrossRef] [Google Scholar]
  69. Miniati, F., Ryu, D., Kang, H., & Jones, T. W. 2001, ApJ, 559, 59 [NASA ADS] [CrossRef] [Google Scholar]
  70. Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
  71. Moore, B., Ghigna, S., Governato, F., et al. 1999a, ApJl, 524, L19 [NASA ADS] [CrossRef] [Google Scholar]
  72. Moore, B., Quinn, T., Governato, F., Stadel, J., & Lake, G. 1999b, MNRAS, 310, 1147 [NASA ADS] [CrossRef] [Google Scholar]
  73. Morlino, G., & Caprioli, D. 2012, A&A, 538, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59 [Google Scholar]
  75. Nava, L., Gabici, S., Marcowith, A., Morlino, G., & Ptuskin, V. S. 2016, MNRAS, 461, 3552 [NASA ADS] [CrossRef] [Google Scholar]
  76. Nava, L., Recchia, S., Gabici, S., et al. 2019, MNRAS, 484, 2684 [NASA ADS] [CrossRef] [Google Scholar]
  77. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
  78. Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., et al. 2015, MNRAS, 454, 2092 [NASA ADS] [CrossRef] [Google Scholar]
  79. Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40 [NASA ADS] [CrossRef] [Google Scholar]
  80. Paglione, T. A. D., & Abrahams, R. D. 2012, ApJ, 755, 106 [NASA ADS] [CrossRef] [Google Scholar]
  81. Pais, M., Pfrommer, C., Ehlert, K., & Pakmor, R. 2018, MNRAS, 478, 5278 [NASA ADS] [CrossRef] [Google Scholar]
  82. Pakmor, R., Pfrommer, C., Simpson, C. M., & Springel, V. 2016, ApJ, 824, L30 [NASA ADS] [CrossRef] [Google Scholar]
  83. Pfrommer, C., Pakmor, R., Simpson, C. M., & Springel, V. 2017a, ApJ, 847, L13 [NASA ADS] [CrossRef] [Google Scholar]
  84. Pfrommer, C., Pakmor, R., Schaal, K., Simpson, C. M., & Springel, V. 2017b, MNRAS, 465, 4500 [NASA ADS] [CrossRef] [Google Scholar]
  85. Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464 [NASA ADS] [CrossRef] [Google Scholar]
  86. Rasera, Y., & Teyssier, R. 2006, A&A, 445, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Recchia, S., Blasi, P., & Morlino, G. 2016, MNRAS, 462, 4227 [NASA ADS] [CrossRef] [Google Scholar]
  88. Recchia, S., Blasi, P., & Morlino, G. 2017, MNRAS, 470, 865 [NASA ADS] [CrossRef] [Google Scholar]
  89. Rieder, M., & Teyssier, R. 2016, MNRAS, 457, 1722 [NASA ADS] [CrossRef] [Google Scholar]
  90. Rosdahl, J., Schaye, J., Teyssier, R., & Agertz, O. 2015, MNRAS, 451, 34 [NASA ADS] [CrossRef] [Google Scholar]
  91. Rosdahl, J., Schaye, J., Dubois, Y., Kimm, T., & Teyssier, R. 2017, MNRAS, 466, 11 [NASA ADS] [CrossRef] [Google Scholar]
  92. Rosen, A., & Bregman, J. N. 1995, ApJ, 440, 634 [NASA ADS] [CrossRef] [Google Scholar]
  93. Ruszkowski, M., Yang, H.-Y. K., & Zweibel, E. 2017, ApJ, 834, 208 [NASA ADS] [CrossRef] [Google Scholar]
  94. Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312 [NASA ADS] [CrossRef] [Google Scholar]
  95. Salem, M., Bryan, G. L., & Hummels, C. 2014, ApJl, 797, L18 [NASA ADS] [CrossRef] [Google Scholar]
  96. Salem, M., Bryan, G. L., & Corlies, L. 2016, MNRAS, 456, 582 [Google Scholar]
  97. Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2016, ApJ, 826, 200 [NASA ADS] [CrossRef] [Google Scholar]
  98. Sharma, P., & Hammett, G. W. 2007, J. Comput. Phys., 227, 123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Sharma, P., Colella, P., & Martin, D. F. 2010, SIAM J. Sci. Comput., 32, 3564 [CrossRef] [Google Scholar]
  100. Simpson, C. M., Pakmor, R., Marinacci, F., et al. 2016, ApJ, 827, L29 [NASA ADS] [CrossRef] [Google Scholar]
  101. Smith, M. C., Sijacki, D., & Shen, S. 2019, MNRAS, 485, 3317 [NASA ADS] [CrossRef] [Google Scholar]
  102. Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [Google Scholar]
  104. Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776 [NASA ADS] [CrossRef] [Google Scholar]
  105. Stinson, G., Seth, A., Katz, N., et al. 2006, MNRAS, 373, 1074 [NASA ADS] [CrossRef] [Google Scholar]
  106. Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Annu. Rev. Nucl. Part. Sci., 57, 285 [NASA ADS] [CrossRef] [Google Scholar]
  107. Sutherland, R. S., & Dopita, M. A. 1993, ApJs, 88, 253 [NASA ADS] [CrossRef] [Google Scholar]
  108. Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44 [NASA ADS] [CrossRef] [Google Scholar]
  110. Teyssier, R., Pontzen, A., Dubois, Y., & Read, J. I. 2013, MNRAS, 429, 3068 [NASA ADS] [CrossRef] [Google Scholar]
  111. Thomas, T., & Pfrommer, C. 2019, MNRAS, 485, 2977 [NASA ADS] [CrossRef] [Google Scholar]
  112. Thornton, K., Gaudlitz, M., Janka, H.-T., & Steinmetz, M. 1998, ApJ, 500, 95 [NASA ADS] [CrossRef] [Google Scholar]
  113. Trebitsch, M., Volonteri, M., Dubois, Y., & Madau, P. 2018, MNRAS, 478, 5607 [CrossRef] [Google Scholar]
  114. Trotta, R., Jóhannesson, G., Moskalenko, I. V., et al. 2011, ApJ, 729, 106 [NASA ADS] [CrossRef] [Google Scholar]
  115. Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179 [NASA ADS] [CrossRef] [Google Scholar]
  116. Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, MNRAS, 423, 2374 [NASA ADS] [CrossRef] [Google Scholar]
  117. Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769 [Google Scholar]
  118. Veritas, C., Acciari, V. A., et al. 2009, Nature, 462, 770 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  119. Wadepuhl, M., & Springel, V. 2011, MNRAS, 410, 1975 [Google Scholar]
  120. Wefel, J. P. 1987, Origin and Transport of High Energy Particles in the Galaxy, Tech. rep. [Google Scholar]
  121. Wiener, J., Oh, S. P., & Guo, F. 2013, MNRAS, 434, 2209 [NASA ADS] [CrossRef] [Google Scholar]
  122. Wiener, J., Pfrommer, C., & Oh, S. P. 2017, MNRAS, 467, 906 [NASA ADS] [Google Scholar]
  123. Yan, H., & Lazarian, A. 2002, Phys. Rev. Lett., 89, 281102 [NASA ADS] [CrossRef] [Google Scholar]
  124. Zolotov, A., Brooks, A. M., Willman, B., et al. 2012, ApJ, 761, 71 [Google Scholar]
  125. Zweibel, E. G. 2017, Phys. Plasmas, 24, 055402 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Shocks galactic winds

We measure the shock Mach number ℳ in the simulations following the procedure described in Dubois et al. (2019), where the Mach number is given by the modified Rankine-Hugoniot jump relation for a composite CR-thermal mixture (Pfrommer et al. 2017b)

M 2 = 1 γ e R P C C − [ ( γ 1 + 1 ) + ( γ 1 − 1 ) R P ] ( γ 2 − 1 ) , $$ \begin{aligned} \mathcal{M} ^2=\frac{1}{\gamma _{\rm e}} \frac{\mathcal{R} _{\rm P}\mathcal{C} }{\mathcal{C} -\left[(\gamma _1+1)+(\gamma _1-1)\mathcal{R} _{\rm P}\right](\gamma _2-1)} \, , \end{aligned} $$(A.1)

where ℛP = P2/P1 is the ratio of the total (thermal plus CR) downstream to upstream pressures (respectively indexed 2 and 1), 𝒞 = [(γ2 + 1)ℛP + (γ2 − 1)](γ1 − 1), γi = Pi/ei + 1 for i = {1, 2} and γe = (γPth, 2 + γCRPCR, 2)/P2 for the downstream region. Figure A.1 compares the Mach number for the simulations without CRs and with CR diffusion (κ = 3 × 1028 cm2 s−1). Without CRs, strong shock develop in the immediate circum-galactic medium (ℳ >  10) of the galaxy and in the large-scale galactic wind (ℳ ≃ 4). With CRs and their diffusive transport, most of those large-scale shocks are erased due to the smooth redistribution of the overall dominating CR pressure.

thumbnail Fig. A.1.

Edge-on view of the maximum value of the Mach number along the line of sight for the simulation without CRs (left panel) and the simulation with CR and anisotropic diffusion with κ = 3 × 1028 cm2 s−1 (right panel) for the G9 simulation at time t = 250 Myr.

All Tables

Table 1.

Simulation initial conditions and parameters for the disc galaxy modelled in this paper.

In the text
Table 2.

Models of CR injection and transport.

In the text

All Figures

thumbnail Fig. 1.

Slices of the G9 galaxy, from left to right: gas density, gas vertical velocity, sum of thermal and CR pressures, and ratio of CR to thermal pressure, seen edge-on at 250 Myr for the different simulations as indicated in the panels (one model per row). See Fig. 2 for the four remaining runs. In the noCR case, a low density wind is generated with velocities of a few hundred km s−1; in the Advection model, the disc is puffed up but the wind is much weaker. When adding CR diffusion, the wind is 10 times denser with velocities similar to noCR case. The sum of CR and thermal pressures, PCR + Pth, increases by 2 orders of magnitude and mainly consists of CR pressure. With CR streaming, the morphology is similar to the Advection case, but with a mild wind.
In the text
thumbnail Fig. 2.

Slices of the G9 galaxy, from left to right: gas density, gas velocity, sum of CR and thermal pressures, and ratio of CR to thermal pressure, seen edge-on at 250 Myr for the different simulations as indicated in the panels (one model per row). See Fig. 1 for the other runs. In the noCR case, a low density wind is generated with velocities of a few hundred km s−1; in the Advection model, the disc is puffed up but the wind is much weaker. When adding CR diffusion, the wind is 10 times denser with velocities similar to noCR case. The sum of CR and thermal pressures increases by 2 orders of magnitude and mainly consists of CR pressure. With CR streaming, the morphology is similar to the Advection case, but with a mild wind.
In the text
thumbnail Fig. 3.

Slices of the G9 galaxy, from left to right: gas density, slices of gas density zoomed-in compared to the left hand panels, sum of CR and thermal pressures, and ratio of CR to thermal pressure, seen edge-on at 250 Myr for the different simulations as indicated in the panels (one model per row). See Fig. 4 for the four remaining runs. The density distributions in the disc is smoother with CR injection. The sum of CR and thermal pressures in the disc increases by 1–2 orders of magnitude when adding CRs, and is either dominated by CR pressure, or equally distributed between CR and thermal pressure.
In the text
thumbnail Fig. 4.

Slices of the G9 galaxy, from left to right: gas density, slices of gas density zoomed-in compared to the left hand panels, sum of CR and thermal pressures, and ratio of CR to thermal pressure, seen edge-on at 250 Myr for the different simulations as indicated in the panels (one model per row). See Fig. 3 for the other runs. The density distributions in the disc is smoother with CR injection. The sum of CR and thermal pressures in the disc increases by 1–2 orders of magnitude when adding CRs, and is either dominated by CR pressure, or equally distributed between CR and thermal pressure.
In the text
thumbnail Fig. 5.

Top panel: SFR as a function of time for the G9 galaxy with and without CR injection, with different CR transport models. Bottom panel: stellar mass as a function of time. Without CR injection (grey thick line), the SFR is 2–3 times greater than with CR injection. The simulation with advection of the CRs is where the SFR is the most affected because CR energy is trapped in the regions of star formation and does not escape.
In the text
thumbnail Fig. 6.

Mass-weighted probability density function of the gas density within a disc of 4 kpc height and 10 kpc radius, centered on the G9 galaxy, for the different runs, at t = 100 Myr (top panel) and t = 250 Myr (bottom panel) throughout the simulations as indicated. The vertical black dashed line shows the density threshold for star formation. Because it adds an additional pressure support against gravity that cools less efficiently than thermal pressure, CR injection reduces the high-density tail of the probability density function and increases the low-density tail.
In the text
thumbnail Fig. 7.

Mass outflow rates of the outflowing gas, mass loading factors, mass-weighted average velocities and mass-weighted average density of the outflowing gas across planes at 2 (left column) and 10 kpc (right column) above the disc plane of the G9 galaxy. The outflow rate quickly drops without CR injection and transport, in the noCR and Advection models. In the Advection-only model, the lack of CR transport prevents the generation of sustained outflows. The outflow rate in the Streaming Boost simulation is very similar to the Advection case at early times and rises to higher values at later times, but the outflow rate is sill weaker than with the other transport models. The isotropic diffusion model and anisotropic diffusion with κ = 3 × 1028 − 1 × 1029 cm2 s−1 are more than 10 times more efficient at driving winds compared to the noCR case. The mass average outflow velocity is less affected by the injection of CRs: it is increased at most by a factor of 2 for the highest diffusion coefficient, and is even lower for κ = 3 × 1027 cm2 s−1.
In the text
thumbnail Fig. 8.

Thermal pressure gradient (a, two upper rows), and CR pressure gradients (b, two bottom rows) relative to the gravitational vertical pull, in 20 kpc slices of the G9 galaxy viewed edge-on, for the different runs, as indicated on the panels. A positive (red) ratio corresponds to a pressure gradient that pushes outwards. A linear threshold of 1 is used for the symmetric logarithmic colourbar. In the absence of CR injection or transport (noCR and Advection), the pressure above and below the mid-plane is dominated by thermal pressure gradient that is maximal in shocked regions. When including CR transport, the gradient of the CR pressure is much higher than that of the thermal pressure, and the gradient of the thermal pressure is lower and less shocks form. The ratio of the thermal pressure gradient to gravity exceeds unity only within shocks whereas the CR pressure gradient in the wind becomes uniformly greater (2 to 10 times) than the gravitational pull for κ ≥ 3 × 1028 cm2 s−1.
In the text
thumbnail Fig. 9.

Total mass of the outflowing gas showed in bins of temperature-density, after 250 Myr for the noCR, Isodiff, κ = 3 × 1028 cm2 s−1 and Streaming boost runs as indicated in the corresponding panels. We use a lower threshold of 10 km s−1 to detect the outflowing gas, selected only above 2 kpc from the disc plane. Without CR injection, the densest (10−4 − 10−3 cm−3) component of the outflowing gas has temperatures of 105 − 105.5 K. When including CR injection with transport, that component is colder (104 K) and denser (10−4 − 10−2 cm−3). In the Advection case, the wind component is almost nonexistent. The wind component is very weak with CR streaming but with similar temperatures to simulations with CR diffusion. The upper left part of the phase diagram corresponds to gas of the intergalactic medium at the outer edges of the disc.
In the text
thumbnail Fig. 10.

Total CR energy in the box the different runs of the G9 galaxy. The higher the diffusion coefficient, the higher the CR energy: the reason is that the more CRs can diffuse, the more they can escape high gas density and CR-energy-density regions where the CR cooling rate is higher. At the beginning of the simulation, with isotropic diffusion, CR energy quickly escapes in directions perpendicular to the disc, whereas it remains trapped in the disc for anisotropic diffusion.
In the text
thumbnail Fig. 11.

Top panel: hadronic and Coulomb cooling rate of CR energy inside the box for the different runs. Bottom panel: CR adiabatic losses (−PCR∇⋅u <  0). We find that the cooling rate is higher for lower diffusion. We find that adiabatic losses are subdominant except in the isotropic and anisotropic diffusion simulations with the highest diffusion coefficients. The black dotted line corresponds to the streaming losses (−ust ⋅ ∇PCR) in the “Streaming Boost” simulation.
In the text
thumbnail Fig. 12.

CR energy-weighted density PDF a disc centred on the G9 galaxy (4 kpc height and 10 kpc radius), at 250 Myr. The upper axis shows the corresponding hadronic and coulomb cooling time. The vertical black dashed line is the density threshold for star formation.
In the text
thumbnail Fig. 13.

Evolution of the (a) mass averaged magnetic field B in the disc (of 4 kpc height and 10 kpc radius); (b) mass averaged ratio of Bz over B in the disc (c) mass averaged Bz in the wind (more than 2 kpc away from the disc) (d) mass averaged ratio of Bz over B in the disc wind; mass average norm of B in the wind. The mean magnetic field strength in the disc saturates at ⟨B⟩≃5 μG with a similar vertical structure in all of the simulated cases. In the wind, a clear difference appears when adding diffusion (isotropic or anisotropic), and the vertical component is more than ten times greater with CR injection and diffusion than without CR injection.
In the text
thumbnail Fig. 14.

Markers show the time averaged (between 75 – 250 Myr) outflow rates (top) mass loading factors (middle) and SFR (bottom) as a function of the circular keplerian velocity of G8 (41 km s−1) and G9 (90 km s−1) at the galactic scale radii. For G9, some of the markers are shifted for clarity. The data points are taken from Heckman et al. (2015) and Chisholm et al. (2017). The labels with the suffix “-conv” are convergence study runs and are discussed in Sect. 5. The poloidal (blue diamond) run is discussed in Sect. 6.
In the text
thumbnail Fig. A.1.

Edge-on view of the maximum value of the Mach number along the line of sight for the simulation without CRs (left panel) and the simulation with CR and anisotropic diffusion with κ = 3 × 1028 cm2 s−1 (right panel) for the G9 simulation at time t = 250 Myr.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.