First Results from the TNG50 Simulation: Galactic outflows driven by supernovae and black hole feedback

TL;DR

This paper introduces the TNG50 simulation, revealing detailed insights into galactic outflows driven by supernovae and black hole feedback, including their velocities, mass loading, and complex behaviors across galaxy types.

Contribution

First results from the high-resolution TNG50 simulation demonstrate detailed modeling of galactic outflows and their dependence on galaxy mass and black hole activity.

Findings

Outflow mass loading varies non-monotonically with galaxy mass.

Outflow velocities increase with stellar mass and redshift.

High velocity, multi-phase outflows can reach over 3000 km/s.

Abstract

We present the new TNG50 cosmological, magnetohydrodynamical simulation -- the third and final volume of the IllustrisTNG project. This simulation occupies a unique combination of large volume and high resolution, with a 50 Mpc box sampled by 2160^3 gas cells (baryon mass of 8x10^4 Msun). The median spatial resolution of star-forming ISM gas is ~100-140 parsecs. This resolution approaches or exceeds that of modern 'zoom' simulations of individual massive galaxies, while the volume contains ~20,000 resolved galaxies with M*>10^7 Msun. Herein we show first results from TNG50, focusing on galactic outflows driven by supernovae as well as supermassive black hole feedback. We find that the outflow mass loading is a non-monotonic function of galaxy stellar mass, turning over and rising rapidly above 10^10.5 Msun due to the action of the central black hole. Outflow velocity increases with stellar mass, and at fixed mass is faster at higher redshift. The TNG model can produce high velocity, multi-phase outflows which include cool, dense components. These outflows reach speeds in excess of 3000 km/s out to 20 kpc with an ejective, BH-driven origin. Critically, we show how the relative simplicity of model inputs (and scalings) at the injection scale produces complex behavior at galactic and halo scales. For example, despite isotropic wind launching, outflows exhibit natural collimation and an emergent bipolarity. Furthermore, galaxies above the star-forming main sequence drive faster outflows, although this correlation inverts at high mass with the onset of quenching, whereby low luminosity, slowly accreting, massive black holes drive the strongest outflows.

Tables

Table 1. Table 1: Key details of the TNG50 simulation in comparison to its larger volume siblings. These are: the simulated volume and box side-length (both comoving), the number of initial gas cells, dark matter particles, and Monte Carlo tracers. The mean baryon and dark matter particle mass resolutions, in solar masses. The minimum allowed adaptive gravitational softening length for gas cells (comoving Plummer equivalent), the redshift zero softening of the collisionless components in physical parsecs, and CPU time.

Run Name		TNG50	TNG100	TNG300
Volume	[ Mpc3 ]	${51.7}^{\mathrm{𝟑}}$	${110.7}^3$	${302.6}^3$
$L_{\mathrm{box}}$	[ Mpc/$h$ ]	$\mathrm{𝟑𝟓}$	75	205
$N_{\mathrm{GAS}}$	-	${\mathrm{𝟐𝟏𝟔𝟎}}^{\mathrm{𝟑}}$	${1820}^3$	${2500}^3$
$N_{\mathrm{DM}}$	-	${\mathrm{𝟐𝟏𝟔𝟎}}^{\mathrm{𝟑}}$	${1820}^3$	${2500}^3$
$N_{\mathrm{TR}}$	-	${\mathrm{𝟐𝟏𝟔𝟎}}^{\mathrm{𝟑}}$	$2\times {1820}^3$	${2500}^3$
$m_{\mathrm{baryon}}$	[ M⊙ ]	$8.5\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟒}}$	$1.4\times {10}^6$	$1.1\times {10}^7$
$m_{\mathrm{DM}}$	[ M⊙ ]	$4.5\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟓}}$	$7.5\times {10}^6$	$5.9\times {10}^7$
${ϵ}_{\mathrm{gas},\min }$	[ pc ]	$\mathrm{𝟕𝟒}$	185	370
${ϵ}_{\mathrm{DM},\mathrm{stars}}$	[ pc ]	$\mathrm{𝟐𝟖𝟖}$	740	1480
CPU Time	[ Mh ]	$\mathrm{𝟏𝟑𝟎}$	18	35

Equations

$$t_{\star }(\rho )=t_0^{\star }{\left(\frac{\rho }{{\rho }_{\mathrm{th}}}\right)}^{-\alpha }$$

$$\dot{M}=\left.\frac{\partial M}{\partial t}\right\rvert_{\rm\,rad}=\frac{1}{\Delta r}\sum_{\overset{i=0}{|r_{i}-r_{0}|<\Delta r/2}}^{N}\left(\frac{\mbox{\boldmath$v_{i}$}\cdot\mbox{\boldmath$r_{i}$}}{|r_{i}|}m_{i}\right)$$

$$\dot{M}=\dot{M}\left(z,M_{\star };r,v_{\mathrm{rad}},\rho ,T,Z\right).$$

$${\eta }_{\mathrm{M}}^{\mathrm{SN}}=\frac{{\dot{M}}_{\mathrm{out}}}{{\dot{M}}_{\star }}.$$

$$L_{\mathrm{bol}}=\left\{\begin{array}{rlr}& {ϵ}_r\dot{M}{c}^2& ;{\lambda }_{\mathrm{edd}}\ge 0.1\\ & 10{\lambda }_{\mathrm{edd}}^2{L}_{\mathrm{Edd}}& ;{\lambda }_{\mathrm{edd}}<0.1\end{array}\right\},$$

$${\dot{M}}_{\mathrm{out}}=\Omega r^2\mu (r)n_{\mathrm{H}}(r)v_{\mathrm{out}}(r)$$

$${\dot{M}}_{\mathrm{out}}(r)=\Omega r\mu (r)N_{\mathrm{H}}{v}_{\mathrm{out}}(r).$$

$${\dot{M}}_{\mathrm{out}}=M_{\mathrm{gas}}{v}_{\mathrm{out}}{r}_{\mathrm{out}}^{-1}.$$

Full text

††volume: submitted

1 Introduction

Energetic outflows driven out of galaxies signpost the resistance of astrophysical feedback processes against the inevitability of gravitational collapse. By itself, the accretion of gas into dark matter halos and its subsequent cooling onto a centrally forming galaxy (Silk, 1977; White & Rees, 1978) leads to an excess of baryons across the mass scale (Thoul & Weinberg, 1995). This is evidenced most clearly by the differential shape of the galaxy stellar mass function with respect to the underlying dark matter halo mass function from which it must emerge (Kauffmann et al., 1993; Cole et al., 1994).

To successfully regulate the stellar mass content of galaxies gas must be prevented from joining the star-forming phase, efficiently ejected out of this phase, or both. Preventive feedback alone is not thought to be sufficient (Benson et al., 2003), while galactic outflows are both theoretically expected (Tomisaka & Ikeuchi, 1988) and observationally detected (Heckman et al., 1990). Physically, they are expected to be driven by the energy released from stars and supernovae (Chevalier & Clegg, 1985), as well as from supermassive black holes (Begelman et al., 1991).

No model for galaxy formation currently succeeds without incorporating both mechanisms (Springel & Hernquist, 2003; Croton et al., 2006). Outflows and feedback are therefore a fundamental aspect of galaxy formation. However, the underlying physical mechanisms are complex (Murray et al., 2005) and operate at extremely small scales – a supermassive black hole accretion disk (10*-5* pc) is roughly ten orders of magnitude smaller than its host galaxy or the intracluster medium (105 pc) whose thermodynamics it can dominate (Tabor & Binney, 1993; McNamara et al., 2000; Fabian et al., 2003). As feedback-driven outflows are both difficult to simulate and difficult to observe, the topic invites study (reviewed in Veilleux et al., 2005; King & Pounds, 2015; Heckman et al., 2017).

1.1 Observations of galactic outflows

Pioneering studies have demonstrated the existence of outflows in the local universe (Heckman et al., 1990; Lehnert & Heckman, 1996) and at high redshifts $z>2$ (Pettini et al., 2001). Galactic outflows are understood to be ‘ubiquitous’ in star-forming galaxies (Shapley et al., 2003; Martin, 2005; Weiner et al., 2009; Zhu et al., 2015), implying that normal main-sequence galaxies continuously drive winds. Observational probes for the incidence and properties of outflows extend from $z\sim 0$, including analyses of SDSS (Chen et al., 2010) through $z\sim 0.5-1$ (Rubin et al., 2010, 2011) to the peak epoch of cosmic star formation at $z\sim 2-3$ in both absorption (Steidel et al., 2010) and emission, focusing on star formation (SF)-driven outflows (Newman et al., 2012; Genzel et al., 2014) and blackhole/active galactic nuclei (AGN) sources (Harrison et al., 2016; Leung et al., 2017; Circosta et al., 2018) even up to $z\gtrsim 5$ quasars (e.g. Bischetti et al., 2018).

Gas which makes up such a galactic-scale wind can be identified in a number of different phase tracers (Cicone et al., 2018): cold gas (e.g. $\lesssim 10^{4}$ K) including molecules such as CO (Fluetsch et al., 2018) as well as neutral HI and metals such as NaD (Concas et al., 2017; Bae & Woo, 2018). Cool gas (e.g. $\sim 10^{4}$ K) including metal ion tracers such as FeII, MgII and OIII (Martin & Bouché, 2009; Rubin et al., 2014), as well as hydrogen in the form of H$\alpha$ (Shapiro et al., 2009). Warm gas (e.g. $\sim 10^{5}$ K - $10^{6}$ K) can also be traced by ionized metals such as OVI (Kacprzak et al., 2015; Nielsen et al., 2017). Finally, hot gas (e.g. $\gtrsim 10^{6}$ K) as observed at x-ray wavelengths (Lehnert et al., 1999; Strickland & Heckman, 2009).

Except in the most local examples, observations of more than one of these outflow phases in the same system are rare. When available, such data imply that outflows are multi-phase, having several co-spatial, possibly kinematically coherent components with a large range in density and temperature (Heckman et al., 2017). Different phases can have significantly different inferred outflow velocities (Grimes et al., 2009). Typically, information is only available on one phase, so an inferred outflow rate of mass, energy, or momentum may only be a fraction of the total (Förster Schreiber et al., 2018a). This is particularly problematic when comparing data and theoretical models, for example of mass loading factors $\eta=\dot{M}_{\rm out}/\dot{M}_{\star}$ or outflow velocities, as we explore below.

Only in the past few years have we begun to spatially resolve outflows with integral field units (IFU) surveys (Rupke et al., 2017). Recent explorations with MUSE have demonstrated the power of blind, deep fields out to high-z (Finley et al., 2017; Feltre et al., 2018) as well as resolved detail at low-z (Venturi et al., 2018). Data from the KCWI instrument has revealed E+A galaxies hosting strong conical outflows, and shown how MgII-traced outflows can be mapped in spatially resolved emission (Baron et al., 2018, \textcolorblueBurchett et al. in prep). Studies using SINFONI have shown how AO-assisted ground-based IFU data can connect galaxy and outflow properties on resolved, kpc scales at $z\sim 2$ (Förster Schreiber et al., 2018b; Davies et al., 2018; Circosta et al., 2018).

From the observational point of view, even identifying a given data feature as an outflow can be challenging. For direct ‘down the barrel’ spectral observations of a host galaxy, two orthogonal techniques are at play: absorption and emission. A blue-shifted absorption line can be interpreted as an intervening mass of the given gas tracer, flowing towards the observer and so away from the galaxy, blocking some of its background continuum light. While the line center offset from systemic gives an indication of the velocity of the bulk of the material, the shape (equivalent width) simultaneously encodes the full velocity distribution of the outflowing material, the covering factor of the observed phase, and the optical depth (Chen et al., 2010). Optical emission-line profiles are often asymmetric in shape, plausibly consisting of superimposed narrow and broad components, the latter preferentially blueshifted. These can be understood as a galactic-scale outflow, where obscuration by the interstellar medium (ISM) of the host galaxy hides the receding, redshifted gas (Armus et al., 1989).

Studies in absorption as well as emission must make a large number of assumptions to convert observed quantities to physical values. There are enough caveats that final uncertainties are large (Chisholm et al., 2017). Given the breadth of techniques, definitions, galaxy selection functions, mass and redshift coverage, and physical assumptions, even fundamental scaling relations between outflow properties (e.g. velocity, mass/energy loading) and host galaxy properties (e.g. stellar mass, star formation rate, black hole luminosity) are often revised (Rupke et al., 2005; Perna et al., 2017).

As a result, and despite the abundance of evidence for the existence of outflowing gas around galaxies, questions of its origin (launch mechanism, energy source) as well as its eventual fate (ability to escape the galaxy, dark matter halo, or to return as a recycled fountain flow) remain largely open topics. Of particular interest is our ability to discriminate winds with stellar versus black hole feedback origins. In particular, the correlations of outflow properties with the star formation and black hole activity of the galaxy itself can shed light on this dichotomy (Förster Schreiber et al., 2018a), ultimately providing a crucial view of, and constraint on, astrophysical feedback processes.

1.2 Cosmological hydrodynamical simulations

One promising avenue to study this interplay of feedback, outflows, and the subsequent cycle of baryons in and around galaxies is through cosmological hydrodynamical simulations (see review in Somerville & Davé, 2015). Recent large-volume efforts have clarified the need for efficient outflows in order to make realistic galaxy populations whose characteristics and integral properties are roughly in agreement with observational constraints (Vogelsberger et al., 2014; Genel et al., 2014; Crain et al., 2015; Schaye et al., 2015; Dubois et al., 2016). These simulations have started to disentangle the dominant physics driving galaxy formation, understanding not only mean relations and ‘typical’ galaxies (Genel et al., 2018), but also outliers and unique examples (Zhu et al., 2018) – providing, as a result, interpretation for a diverse array of observations.

The IllustrisTNG project (Pillepich et al., 2018b; Naiman et al., 2018; Nelson et al., 2018a; Marinacci et al., 2018; Springel et al., 2018) is a recent addition to this class of simulation. The first two simulations of this effort, TNG100 and TNG300, realize volumes of $\sim 100^{3}$ and $\sim 300^{3}$ comoving Mpc3, respectively. This allows them to sample galaxy populations which include hundreds of thousands of objects with $M_{\star}\geq 10^{9}$ M*⊙* by $z=0$. However, the unavoidable trade off in such population-level statistical power is – by construction – limited resolution.

With baryon mass resolutions of $\sim$ 106 M*⊙* (corresponding roughly to the canonical ‘1 kpc’ spatial resolution) of the current generation of cosmological volumes, we can realistically study the structure of galaxies with $M_{\star}\gtrsim 10^{9}$ M*⊙*. This hinders investigation of many interesting scientific questions, particularly those related to the structural properties of the stellar and gaseous components of galaxies, the internal structure of disks and their ISM, the co-evolution of black holes and galactic nuclei, feedback, and the resulting baryon cycle of inflows and outflows.

These small-scale phenomena are frequently investigated with ‘zoom’ simulations of individual galaxies (Hopkins et al., 2014), with project campaigns of order 10 to 100 galaxies in total (Wang et al., 2015), often oriented around a focus such as Milky Way mass halos (Grand et al., 2017), Local Group analogs (Sawala et al., 2016) or the challenging regime of galaxies in high-density cluster environments (Bahé et al., 2017; Tremmel et al., 2018). However, zoom simulations immediately sacrifice one of the greatest advantages of cosmological simulations – namely, the ability to make statistically robust statements about large, unbiased populations and thereby comment on galaxy formation and evolution in full generality.

Together with the companion paper (\textcolorbluePillepich et al. 2019) we here present the new TNG50 simulation, the third and final volume of the IllustrisTNG project. TNG50 has been designed to overcome the resolution/volume limitation by simulating a large, cosmological volume at a resolution which approaches or even exceeds that of modern ‘zoom’ simulations of individual galaxies. It enables an insightful view into the structure, chemodynamical evolution, and small-scale properties of galaxies and their halos.

1.3 Simulation work on outflows

At small scales there are many idealized studies aimed at understanding the details of supernova-driven wind launching mechanisms (Girichidis et al., 2016; Fielding et al., 2017; Kim & Ostriker, 2017) or detailed phase-interaction physics (Richings & Faucher-Giguère, 2018; Schneider & Robertson, 2018). General relativistic (GR)MHD simulations explore black hole outflow phenomena including both winds and jets (Blandford & Begelman, 1999; McKinney et al., 2014), recently reviewed in Yuan & Narayan (2014).

One of the long-term goals of such efforts is to capture and parameterize the emergent behavior at larger scales in order to develop effective models which can be used in full galaxy or cosmological simulations. The sophistication of such effective models is presently limited by several factors, most critically the finite numerical resolution available. As it will be impossible to ever resolve the energy injection scales of critical astrophysical processes, sub-resolution (or sub-grid) approximations are unavoidable.

As a result, cosmological simulations have rarely if ever been compared directly to outflow observations (i.e. the baryon cycle), but rather to indirect properties such as the amount of metals left in the ISM (Torrey et al., 2017), the ionized metal content of the circumgalactic medium (Nelson et al., 2018b), or the stellar formation efficiency of the galaxy itself (Pillepich et al., 2018b).

Notably, Oppenheimer et al. (2010) studied SN wind-driving models in early cosmological simulations, emphasizing the importance of recycling and the dependence of recycled accretion on the balance between mass loading and outflow velocity. In the context of the more realistic FIRE model, Muratov et al. (2015) presented measurements of the scaling of mass loading (and outflow velocity) which decrease (increase) strongly with $M_{\star}$, such that lower-mass progenitors at early times would have had much higher mass loadings than their present day descendants. Similarly, Christensen et al. (2016) evaluated outflow scalings (and metal content; Christensen et al., 2018) over a similar mass range (halos $<10^{12}$ M*⊙*) and zoom sample size ($\sim$10 - 20) using a delayed cooling blastwave SN feedback model, emphasizing the importance of ejective feedback (i.e. high mass loadings) towards low $M_{\star}$.

Non-cosmological, idealized simulations of dwarf galaxy disks enable the possibility of higher resolution and sophisticated physical modeling, and have started to explore the phase-structure and feedback dependencies of galactic outflows. Emerick et al. (2018) highlight how different treatments of the radiative feedback from young stars can modify the mass loading and temperature of outflows. Hu (2018) demonstrate how different SN energy injection schemes similarly lead to different outflow properties – in particular, that a terminal momentum injection approach does not produce the same pressure-driven wind as when residual thermal energy is also accounted for (see also Smith et al., 2018).

None of these works have emphasized any direct observational signatures of their outflows, i.e. forwarding modeling. Recently, Tescari et al. (2018) made a first comparison between EAGLE ‘warm gas’ and spatially resolved SAMI/IFU H$\alpha$ emission to relate gas kinematics to wind signatures. The notable exception is Ceverino et al. (2016), who generated synthetic H$\alpha$ emission line profiles of $z\sim 2$ zoom galaxies and correlated their information content with intrinsic outflow properties. Any robust comparison between outflow observations and hydrodynamical simulations will require such progress in the future, as we discuss below.

All of these investigations have focused exclusively on SN-driven outflows, mostly in models specifically neglecting BH feedback, which is however thought to play a dominant role in outflow production in massive galaxies. For example, the Muratov et al. (2015) measurement of $\eta_{\rm M}$ scalings in FIRE fails to capture the onset of the BH mechanism at $M_{\star}\gtrsim 10^{10.5}$ M*⊙* and therefore misses the inversion of this relation towards higher masses, as we demonstrate below. A notable exception is Brennan et al. (2018), which explores a kinetic BH feedback model in a suite of 24 zooms (halos $10^{12}-10^{13.5}$ M*⊙*), showing how black holes can produce faster outflows (up to 1000 km/s) which travel farther and are thus less likely to be subsequently recycled (see also Nelson et al., 2015).

The detailed properties (as opposed to the consequences) of SN and BH-driven outflows, as well as the interplay between these two origin mechanisms, have not yet been studied with cosmological simulations, nor across a full sampling of the galaxy mass range, from $M_{\star}=10^{7}$ M*⊙* dwarfs to $M_{\star}=10^{11.5}$ M*⊙* BCGs.

1.4 The present work

We first present the TNG50 simulation, a new class of cosmological volume simulation realized at a resolution of modern ‘zoom’ simulations of individual galaxies. Leveraging this new tool, we describe first results focusing on the properties of galactic outflows in the cosmological setting and with respect to the galaxies from which they arise. In this context we are not trying to understand the launching mechanism or the relevant microphysics, but rather the way in which outflows shape the galaxy population as a whole, modulate galaxy evolution, and generate associated observational signatures. As emphasized above, the feedback models in simulations such as TNG are effective in nature, meaning that the outflow properties of the model at injection are an input rather than an output of the model. We therefore contrast, throughout this work, our results on the emergent properties of galactic outflows with the assumptions of the underlying physical models.

TNG50 incorporates unchanged a robust scheme for cosmological magnetohydrodynamics coupled to a comprehensive and well-validated model for galaxy formation physics (Weinberger et al., 2017; Pillepich et al., 2018a). This motivates our study of the resulting outflows, as they reflect the net outcome of a model which produces a realistic galaxy population matching many key properties and scaling relations of observed galaxies (see \textcolorbluePillepich et al. (2019) in the context of TNG50 galactic structural properties). The TNG model was designed via calibration against a number of observations, mainly at $z=0$ (see Pillepich et al., 2018a, for a detailed discussion of our calibration/tuning), and parameters related to feedback processes at the smallest scales are arguably the most important ingredient of such a model. We focus here on intermediate and high redshift, namely $z=1$, $z=2$, and above, as particularly interesting and observationally rich epochs of galaxy assembly and galactic-scale outflows.

The paper is organized as follows: in Section 2.1 we introduce the TNG50 simulation and the TNG model for galaxy formation, while Section 2.3 details the methodology and analysis details for our study of outflows. Section 3 presents our results: the high redshift galaxy population of TNG50 (Section 3.1); measurement of mass outflow rates and mass loading factors (Section 3.2); outflow velocities and the ability of BHs to generate $\gtrsim$ 3000 km/s flows (Section 3.3); the multi-phase character of our outflows (Section 3.4); the angular dependence of winds and their hydrodynamic collimation (Section 3.5); and correlations between outflows and the properties of their host galaxy and central black hole, in comparison to observed relations (Section 3.6). Section 4 discusses future directions and comparison with observations, while Section 5 summarizes our conclusions. Finally, Appendix A presents halo-normalized outflow velocity measurements for comparison.

2 Methods

2.1 The TNG Simulations

The IllustrisTNG project111http://www.tng-project.org (Nelson et al., 2018a; Naiman et al., 2018; Pillepich et al., 2018b; Marinacci et al., 2018; Springel et al., 2018) is a series of three large cosmological volumes, simulated with gravo-magnetohydrodynamics (MHD) and incorporating a comprehensive model for galaxy formation physics (Weinberger et al., 2017; Pillepich et al., 2018a). All aspects of the model, including parameter values and the simulation code, are described in these two methods papers and remain in practice unchanged for our production simulations, and we give here only a brief overview.

The TNG project includes three distinct simulation volumes: TNG50, TNG100, and TNG300. The larger two have previously been presented: these are TNG100, which includes 2$\times$18203 resolution elements in a $\sim$ 100 Mpc (comoving) box, and TNG300, which includes 2$\times$25003 resolution elements in a $\sim$ 300 Mpc box. Here we present the third, final, and by far most computationally demanding simulation of the IllustrisTNG project: TNG50, our high-resolution effort. This run includes 2$\times$21603 resolution elements in a $\sim$ 50 Mpc (comoving) box. In TNG50 the baryon mass resolution is $8.5\times 10^{4}$ M*⊙*, the gravitational softening lengths are $290$ parsecs at $z$ = 0 for the stars and dark matter, while the gas softening is adaptive with a minimum of $74$ comoving parsecs. TNG50 has roughly fifteen times better mass resolution, and two and a half times better spatial resolution, than TNG100. This represents an unprecedented combination of resolution and volume for a cosmological hydrodynamical simulation, which approaches or exceeds that of modern ‘zoom’ simulations of individual galaxies. The main parameters of the three TNG volumes are compared in Table 1.

TNG uses the Arepo code (Springel, 2010) which solves for the coupled evolution under (self-)gravity and ideal, continuum MHD (Pakmor et al., 2011; Pakmor & Springel, 2013). Gravity employs the Tree-PM approach, and the fluid dynamics use a Godunov type finite-volume scheme where an unstructured, moving, Voronoi tessellation provides the spatial discretization. The simulations include a physical model for the most important processes relevant for the formation and evolution of galaxies. Specifically: (i) gas radiative processes, including primordial/metal-line cooling and heating from the background radiation field, (ii) star formation in the dense ISM, (iii) stellar population evolution and chemical enrichment following supernovae Ia, II, as well as AGB stars, with individual accounting for the nine elements H, He, C, N, O, Ne, Mg, Si, and Fe, plus NS-NS byproducts, (iv) supernova driven galactic-scale outflows or winds, (v) the formation, coalescence, and growth of supermassive black holes, (vi) and dual-mode BH feedback operating in a thermal ‘quasar’ state at high accretion rates and a kinetic ‘wind’ state at low accretion rates. Note that TNG does not yet include an explicit treatment of radiation nor cosmic rays, both of which could plausibly affect the phase structure of outflows.

Black holes are seeded in massive halos and then accrete nearby gas at the Eddington limited Bondi rate. Based on this accretion rate, their feedback mode is determined. When the Eddington ratio exceeds a threshold of $\chi=\rm{min}[0.002(M_{\rm BH}/10^{8}\rm{M}_{\odot})^{2},0.1]$ thermal energy is injected continuously into the surrounding gas. The rate is $\Delta E_{\rm high}=\epsilon_{\rm f,high}\epsilon_{\rm r}\dot{M}_{\rm BH}c^{2}$ where $\epsilon_{\rm f,high}$ is the high-state coupling efficiency, $\epsilon_{\rm r}$ is the radiative accretion efficiency, and $\epsilon_{\rm f,high}\epsilon_{\rm r}=0.02$. Below this threshold, kinetic energy is injected as a time-pulsed, oriented ‘wind’ with a random direction which reorients for each event. The rate is $\Delta E_{\rm low}=\epsilon_{\rm f,low}\dot{M}_{\rm BH}c^{2}$ with $\epsilon_{\rm f,low}\leq 0.2$ (its typical value at onset), the efficiency decreasing at low environmental density (see Weinberger et al., 2017).

Galactic-scale outflows generated by stellar feedback are modeled using a kinetic wind approach, whereby the available energy from SNII is used to stochastically eject star-forming gas cells from galaxies. The injection velocity of such wind particles is $v_{\rm w}\propto\sigma_{\rm DM}$ where $\sigma_{\rm DM}$ is the local dark matter velocity dispersion, subject to a minimum $v_{\rm w,min}=350$ km/s. The mass loading of the winds is then $\eta=\dot{M}_{\rm w}/\dot{M}_{\rm SFR}=2(1-\tau_{\rm w})e_{\rm w}/v_{\rm w}^{2}$ where $\tau_{\rm w}=0.1$ is the thermal energy fraction and $e_{\rm w}$ is a metallicity dependent modulation of the canonical $10^{51}$ erg available per SNII. Wind particles are hydrodynamically decoupled from surrounding gas until they exit the dense, star-forming environment. This occurs when the wind reaches a background density lower than 0.05 times the star-formation threshold, or after 0.025 of the current Hubble time elapses, and typically within a few kiloparsecs. The total energy available to drive winds from a gas cell therefore depends on its instantaneous star formation rate as $\dot{E}_{\rm w}=e_{\rm w}\dot{M}_{\rm SFR}$ which is roughly $\simeq 10^{41}-10^{42}$ erg/s $\dot{M}_{\rm SFR}/(M_{\odot}/\rm{yr})$ depending on the local gas metallicity (see Pillepich et al., 2018a).

In Figure 1 we offer a comparison between recent cosmological hydrodynamical simulations of galaxy formation, putting TNG50 in context given its combination of volume and resolution, which uniquely positions it in this space. TNG50 is by far the most computationally demanding run of the TNG suite, and one of the most ambitious cosmological simulations to date.

We adopt a cosmology consistent with the Planck Collaboration (2016) results, namely $\Omega_{\Lambda,0}=0.6911$, $\Omega_{m,0}=0.3089$, $\Omega_{b,0}=0.0486$, $\sigma_{8}=0.8159$, $n_{s}=0.9667$ and $h=0.6774$.

2.2 Model update for TNG50

To execute the TNG50-1 (highest resolution level of the TNG50 volume) simulation a minor change has been made to the fiducial TNG model. With the inclusion of MHD, we found that extremely dense gas cells often have $P_{\rm B}\gg P_{\rm gas}$. The resulting Courant-like constraint on the hydrodynamical timestep, restricted by the magnetic signal velocity instead of the sound speed, requires physical timesteps as small as 10 years at the resolution of TNG50-1. Such high-density gas should naturally convert rapidly into stars, but in the fiducial configuration the numerical timestep decreases faster than the star formation probability increases (with increasing gas density), leading to an unfavorable race condition. We instead want that gas cells with extremely short star formation timescales convert into collisionless star particles in a reasonable number of numerical timesteps. We have therefore modified the base Springel & Hernquist (2003) star formation model, wherein the dependence of the star formation timescale on density $t_{\star}(\rho)$ is given by (Eqn. 21)

[TABLE]

where $t_{0}^{\star}$ is a model parameter which modulates the global gas consumption time-scale, and $\rho_{\rm th}$ is a model parameter setting the threshold density for star formation. In the fiducial configuration, used for all simulations employing this model, the choice of $\alpha=1/2$ has been made, corresponding to the canonical assumption that $t_{\star}$ is proportional to the local dynamical time of the gas.

Instead, for TNG50-1, we postulate a somewhat steeper dependence of the star formation rate on gas density only for the densest gas. We take $\alpha=1$ for gas above the runaway threshold, which is defined as the density where the slope of the effective pressure versus density curve $n_{\rm eff}<4/3$. At this point the effective pressure can no longer support gas against dynamical instabilities, and this occurs at 24.4 cm*-3*, which is $\sim$ 230 $\rho_{\rm th}$. Since the star formation rate goes as $\dot{M}_{\star}\propto\rho^{\alpha}$, this implies the physical ansatz of more rapid star formation in the densest environments – the nuclear starburst regime of galactic disks, for instance.

In practice, this change was taken for numerical reasons, and we have verified with test simulations at full resolution ($1024^{3}$ particles in 25 Mpc/h volumes) that this change affects a negligible amount of gas, by mass or by number, and has no impact on galactic structure nor galaxy population statistics, at any redshift.

2.3 Measuring Outflow Rates and Velocities

To compute mass outflow rates two orthogonal techniques exist: (i) deriving instantaneous fluxes based on gas kinematics at one point in time (e.g. Ocvirk et al., 2008), and (ii) deriving mass fluxes by tracking the Lagrangian evolution of gas mass across two distinct points in time (e.g. Nelson et al., 2013). As only the former has the possibility for direct comparisons to observables, we here focus on these instantaneous rates. We compute the radial mass flux of all the gas cells in a thin shell centered on a galaxy as

[TABLE]

where the subscript $i$ enumerates all the gas cells with mass $m_{i}$ in a particular volume of space, taken as a spherical shell with some thickness $\Delta r$ and mean distance $r_{0}$ from a central galaxy. We calculate $\dot{M}$ for a range of $r_{0}$ and corresponding $\Delta r$ values, 5 kpc thick in the inner halo coarsening outwards. The cell position $r_{i}$ is relative to the subhalo center, taken as the position of the particle/cell with the minimum gravitational potential, while the velocity $v_{i}$ is relative to the subhalo center of mass motion, accounting for local Hubble expansion. $v_{\rm rad}>0$ denotes outflow, and $v_{\rm rad}<0$ is inflow.

Herein we frequently decompose the outflowing gas by its instantaneous properties, whereby $\dot{M}$ has parameter dependencies:

[TABLE]

To compute this distribution for each galaxy (with a given redshift and stellar mass) we calculate the radial mass flux of each gas cell and take the 5D histogram of this mass flux, or mass outflow rate, binning in the gas properties of interest above: radius, radial velocity, density, temperature, and metallicity. A total mass outflow rate is derived by marginalizing over the state of the gas (i.e. $\rho,T,Z$) and satisfying a threshold in outflow velocity of $v_{\rm rad}>v_{\rm thresh}$. We frequently require only $v_{\rm rad}>0$ km/s to define outflowing material, though higher values are also considered as conservative limits and/or to mimic observational sensitivity limits.

In measuring total mass outflow rates, mass loading factors, or outflow velocities, we consider actual gas cells as well as wind-phase particles, the latter of which are generated by the TNG wind model (see Pillepich et al., 2018a) and are hydrodynamically (though not gravitationally) decoupled from their surrounding gas for the duration of their existence, which is typically restricted to small distances away from galaxies. At $r\gtrsim 10$ kpc, where we generally focus our analysis, wind-phase particles are a negligible contribution and outflows are hydrodynamically resolved. Excluding wind particles has only a minor effect on all our measurements, and we include them as a relevant component of our simulations, one currently modeled at a particularly simple level (i.e. with strong sub-grid assumptions about phase structure and interaction).

To explore trends of outflow velocity, we must reduce the complex distribution of outflow velocities around a galaxy down to a single number $v_{\rm out}$. We do so by taking mass outflow rate weighted percentiles of $v_{\rm rad}$. In particular, we define the quantity $v_{\rm out,N}$ as the radial velocity above which $(1-N)$% of the outflow is moving. For example, $v_{\rm out,90}=500$ km/s implies that ten percent of outflowing mass flux is moving 500 km/s or faster. The ‘median’ value $v_{\rm out,50}$ therefore gives the speed achieved by at least half of the total outflow mass flux. This quantity can be measured at a particular radial distance from the galaxy, $v_{\rm out,50,r=10kpc}$ for example.

We note that these theoretical velocity percentiles may be similar, but not directly comparable to, observational measurements of outflow velocity (e.g. ‘$v_{\rm 05}$’ or ‘$v_{\rm 90}$’ in the literature), because we consider 3D $v_{\rm rad}$ instead of 1D projected $v_{\rm LOS}$, our velocities are mass flux (not optical depth or emission) weighted, and we do not yet attempt any modeling of synthetic spectral features.

We define and measure a mass loading factor with respect to the star formation rate of the galaxy in the usual way,

[TABLE]

This measure is typically used in the context of supernovae-driven winds, to diagnose the relative efficiency of stellar feedback produced galactic-scale outflows. It has the same dependencies as the mass outflow rate, namely $\eta_{\rm M}=\eta_{\rm M}\left(z,M_{\star};r,v_{\rm rad},\rho,T,Z\right)$.

We typically drop the superscript ‘SN’, which emphasizes that the normalization (denominator) assumes that stellar feedback is the relevant energetic source. As BH energy injection is typically coexistent and always present, we always interpret its possible role in setting $\eta_{\rm M}$ or contributing to any measured $\dot{M}_{\rm out}$.

2.4 Galaxy Identification and General Properties

Subhalos are identified with the Subfind algorithm (Springel et al., 2001), and we exclusively consider central (non-satellite) galaxies with $7<\log(M_{\star}/\rm{M}_{\odot})<11$ at redshifts $z\geq 1$.

Throughout this work we take 30 physical kpc aperture values for the stellar mass $M_{\star}$ of a galaxy. We compute star formation rates $\dot{M}_{\star}$ within the same aperture and with a temporal smoothing of 100 Myr, as appropriate for observationally derived SFR indicators based on the infrared or 24$\mu m$ measurements (Calzetti, 2008), using an accounting of recently formed stars, rather than the instantaneous star formation rate of dense gas.

We compute the bolometric luminosities of accreting black holes under the assumption of a $\dot{M}$-dependent radiative efficiency (following Churazov et al., 2005; Weinberger et al., 2018) as

[TABLE]

where $\lambda_{\rm edd}=\dot{M}/\dot{M}_{\rm edd}$ is the Eddington ratio, $L_{\rm edd}=\dot{M}_{\rm edd}c^{2}$ is the Eddington luminosity, and $\dot{M}_{\rm edd}=4\pi GM_{\text{BH}}m_{p}/(\epsilon_{r}\sigma_{\rm T}c)$ the Eddington accretion rate, with the usual fundamental constants. $L_{\rm bol}$ is an intrinsic luminosity – obscuration is not included.

3 Results

3.1 TNG50: Resolving outflows from resolved galaxies

The fundamental means by which feedback regulates galaxy formation and evolution is by generating non-gravitational motions in the gas. Energy injection gives rise to outflows, the subject of this paper, which are particularly evident in the velocity structure of the cosmic gas surrounding galaxies. To set the stage, Figure 2 therefore visualizes the time evolution of a strong BH-driven outflow in TNG50 originating from a massive $M_{\star}\simeq 10^{11.4}$ M*⊙* galaxy at $z\sim 2$. We show gas velocity (left column) and gas temperature (right column), with time progressing downwards, each row roughly 100 Myr apart. Energy injection from the BH produces a high-velocity, large-scale, and highly collimated (directional) outflow. The gas reaches velocities $\gtrsim$ 2500 km/s even as it crosses the halo virial radius (shown as the circle in the last row).

This halo-scale outflow results from numerous, high time cadence energy inputs within a spatial region of $\lesssim 150$ pc: in the center of the upper right panel, five distinct events are already visible as loop-like features in temperature which have begun to flow outwards. The most recent (smallest) has already reached a distance of $\sim 4$ kpc from the SMBH, which has a mass of $10^{8.7}$ M*⊙* and is accreting gas at a rate of $\dot{M}_{\rm BH}\simeq 10^{-2}$ M*⊙*  yr*-1*, such that it is in the ‘low accretion state’ of the TNG model. It has an Eddington ratio of $\lambda_{\rm edd}\simeq 0.002$ and a bolometric luminosity of only $L_{\rm bol}\simeq 10^{42.3}$ erg/s (cf. the system of May et al., 2018). This initial ‘burst’ of energy inputs is followed by continual, though less frequent, feedback from the BH over the next $\sim$500 Myr. The central galaxy is the eighth most massive at this time in TNG50, with a SFR of $\simeq$ 50  M*⊙*  yr*-1*, log(sSFR) $\simeq$ -9.7 yr*-1*, and is in the process of quenching.

We see that the outflow heats two centrally offset lobes to temperatures $\gtrsim 10^{7.5}-10^{8}$ K. A succession of cocoon like shells pile up as they encounter resistance propagating through the halo gas. Although not shown here, energy dissipation due to hydrodynamical shocks traces the outer envelopes (i.e. density edges) of individual bubbles. The outflow geometry is more jet-like in its collimation than bipolar, being quite linear along the horizontal direction. The overall morphology is reminiscent of early simulations by Suchkov et al. (1994) of superwind structure, as well as radio observations of Cecil et al. (2001) who studied the coherence of the magnetic field inside similarly trailing lobes. That work suggested vortex-like dynamics at the edge of the bubble could loft material outwards, while inversion of the observed rotation measure across the boundary implied a separation of outflow and infall regions – signatures in the magnetic field topology which could be explored in TNG50.

The occasional tilt with respect to the host rotation axis (as seen in the first and last rows) could, if the outflow origin was stellar in nature, reflect asymmetries in the distribution of rapidly star-forming gas and young stars. Here it is instead due to a large momentum injection by the BH in a random direction (e.g. Schmitt et al., 2002; Liska et al., 2018) which then propagates along a preferred direction. We discuss this natural collimation in Section 3.5.

The structure is not bi-symmetric as is the case in idealized jet/wind models, and asymmetries of the outflow at larger, halo scales are expected in this more realistic setting (Nelson et al., 2016). In particular, the left lobe is more pronounced, sustaining higher velocities, than the right lobe in the last two time snapshots. This is plausibly due to density fluctuations in the halo gas of the host galaxy. In particular, we see how an intruding satellite at $\Delta t=270$ Myr is shredded as it passes through pericenter, populating the inner halo to the right of the central disk with dense gas structures which likely help inhibit outflow propagation in that direction.

Figure 3 shows the same time sequence, this time in projected gas column density (left) and gas-phase metallicity (right). At the final $\Delta t=370$ Myr the remnants of this merging satellite are clearly visible as a number of dense gas fragments flowing outwards – a feature which could easily be mistaken as a feedback-driven outflow. The most striking feature in the density projections is the bipolar cavity structure, which inflates from a few kpc to the scale of the virial radius over the span of a few hundred Myr. Along the edges of each cavity accumulated shells of swept up and ejected mass are clearly visible as strong density contrasts relative to the under dense, post-shock regions in their wakes. Even without further energy injection from the central engine these highly over pressurized bubbles will expand outwards as they dissipate their energy into the halo gas over time.

The evolving gas metallicity distribution over this time span directly shows how galactic-scale outflows enrich the circumgalactic and intergalactic medium with heavy elements (Tegmark et al., 1993). Although a considerable fraction of the virial volume is already enriched to $0.1Z_{\odot}$ by $z\sim 2$, prior to the onset of strong BH activity, the distribution of metals is inhomogeneous and regions of the halo and its immediate environment remain at $\lesssim 0.01Z_{\odot}$. Sitting already on the mass-metallicity relation, the central galaxy is, however, highly enriched, and this $\gtrsim Z_{\odot}$ gas is directly ejected out to scales of 100s of kpc by the BH-driven wind. This high metallicity material is preferentially cospatial with hot gas in the cavities, particularly at early times. A rich turbulent mixing structure develops, producing metallicity sub-structure within the bubbles and possibly across their boundaries. We discuss the multiphase nature and metal content of the outflows in Section 3.4.

Because the outflow is not isotropic, it allows predominantly metal-poor gas to remain along the major axes of the disk (i.e. in the vertical direction of this projection). Any ongoing accretion will lead to metal dilution over time. Evidence for this fueling of the gas reservoir of the central galaxy, particularly near the end of this strong feedback episode, is also apparent. Cold, overdense gas fragments begin to populate the inner halo at $r\lesssim 0.25r_{\rm vir}$. They originate predominantly along the major axis of the galaxy, though also appear to form directly in the wake of the outflow. These gas overdensities are clearly infalling, with mixing tails extending radially outward, and often not obviously associated with any substructure in the halo at earlier times. We speculate that they may arise from an in-situ, condensation-like process in the thermally unstable hot halo gas (Sharma et al., 2012; McCourt et al., 2012; Voit et al., 2015). Given the ability of TNG50 to resolve detailed structure in the lower density halo gas, a detailed study of precipitation and its ability to feed the central BH and/or regulate the feedback cycle will be an interesting topic for future work.

Sitting at the center of this outflow, nearly invisible at this scale, is the stellar body and gaseous reservoir of the central galaxy. Figure 4 shows one such massive galaxy at $z\sim 2$ in a face-on gas density projection (upper left), and at later times an example of the stellar structure of a $z\sim 1$ descendant of such a system, traced through a mock of its optical stellar light (upper right). In both cases the edge-on view is shown immediately below. The redshift two system hosts a large, flocculent, clumpy spiral structure, roughly $10-20$ kpc in radial extent. This disk is relatively thin in its vertical scale-height at the center, but warped at the outskirts, with gas features coupling the outer disk to ongoing accretion from its larger scale environment. The line-of-sight velocity dispersion of the star-forming (i.e. H$\alpha$-traced) component varies across the disk; between $1-2r_{\rm 1/2,\star}$ the mean value is $\sim 70$ km/s, reaching as high as $170$ km/s – see \textcolorbluePillepich et al. (2019). The maximum rotational velocity is $\sim 350$ km/s, giving $V/\sigma\sim 2-5$ depending on measurement methodology. The redshift one descendant hosts a large stellar bulge surrounded by a thin stellar disk composed of irregular spiral-like features. A stellar clump is visible migrating inwards to eventually coalesce with the bulge. Bombardment by minor mergers at $z<1$ will convert this galaxy into the spheroidal BCG of a massive group by $z=0$.

By showing the highly resolved structure of an individual galaxy we emphasize that TNG50 enables the study of the connection between small-scale (i.e. few hundred pc) feedback and large-scale (i.e. few hundred kpc) outflows. Furthermore, we clarify that winds launched by stellar feedback originate from the entirety of the star-forming gaseous disk, on a cell-by-cell basis at the scale of the simulation resolution. As the local star formation properties vary substantially within a galaxy, the resulting SN-driven outflows have non-constant, and non-trivial, distributions of outflow velocity, mass outflow rate, and even temperature (all at launch). That is, galaxies of a given mass or circular velocity do not produce outflows at a single constant velocity or mass loading, as is sometimes mistakenly inferred for large-volume feedback models like TNG. Simultaneously in time and possibly coincident in space, the central black hole will also be injecting either thermal or kinetic energy into its immediate surroundings – contributing to, or possibly dominating the production of an outflow.

Figure 4 also highlights the statistical power of the galaxy population realized in TNG50 already by $z=2$ (lower half), visualizing the most massive $\sim$ 750 galaxies present in the volume, in the same face-on gas density (left) and face-on stellar light mocks (right) as above. They are sorted from the most massive, in the upper left corner, to the least massive, in the lower right, ranging from $8.0<\log(M_{\star}/\rm{M}_{\odot})<11.7$. With few exceptions, galaxies at this mass and redshift are gas rich, with rotationally supported gas reservoirs which are typically much more extended than the stellar component. The diversity present even within this sample highlights the importance of the large cosmological volume. The $\sim 1.4\times 10^{5}$ cMpc3 of TNG50 provides a sampling of rare, high-mass halos, together with the statistics needed to make robust statements about a representative galaxy population.

3.2 Outflow rates and mass loading factors vs. M⋆ and z

In Figure 5 we begin to explore the quantitative properties of TNG50 outflows. Here we measure the outflow mass loading $\eta_{\rm M}$, the ratio of the mass outflow rate and the galaxy star formation rate, as a function of galaxy stellar mass at $z=2$ (top panel). Considering a fixed distance of 10 kpc from the galaxy, we evaluate the mass outflow rates subject to three different minimum $v_{\rm rad}$ values to be classified as an outflow (blue, orange, and green).

Considering all $v_{\rm rad}>0$ km/s material in the mass outflow rate (blue), $\eta_{\rm M}$ drops rapidly with increasing stellar mass, from $\eta_{\rm M}\simeq 50$ at $M_{\star}=10^{7.5}$ M*⊙* to $\eta_{\rm M}\simeq 3-5$ at $M_{\star}=10^{10.5}$ M*⊙*. Enforcing a more stringent $v_{\rm rad}>150$ km/s cut as measurably outflowing material (orange) weakens the dependence and flattens the relation with stellar mass considerably – $\eta_{\rm M}$ decreases by only a factor of two over the same mass range. For even higher velocity thresholds (green), the slope of $\eta_{\rm M}\propto M_{\star}$ changes sign, because low mass galaxies drive very little mass at such high velocities.

The black line shows the value of $\eta_{\rm M}^{\,i}$ at the injection scale, as prescribed by the TNG model. We derive a single mean value for each galaxy from its ensemble of star-forming gas cells, then present an approximate fit through these mean values as a function of $M_{\star}$, as well as its scatter (gray band). The model $\eta_{\rm M}^{\,i}$ is nearly flat below $M_{\star}\lesssim 10^{8.5}$ M*⊙* and then decreases monotonically with increasing mass, dropping to unity by $10^{10.5}$ M*⊙* and $<1$ thereafter. However, this behavior is fundamentally different than what is seen in the resolved flows which result at 10 kpc scales. Here the trend of $\eta_{\rm M}$ reverses, reaching a minimum at $M_{\star}\simeq 10^{10.5}$ M*⊙* and before rising rapidly again towards higher masses. This corresponds to the mass scale for the onset of quenching as galaxies begin to transition out of the blue, star-forming population in the TNG model (Nelson et al., 2018a; Weinberger et al., 2018). Star formation rates are suppressed by the action of efficient BH feedback which simultaneously produces a large mass of outflowing gas through direct ejection of the nuclear (i.e. central) ISM, driving $\eta_{\rm M}=\dot{M}_{\rm out}/\dot{M}_{\star}$ up. As we move beyond $M_{\star}\gtrsim 10^{11}$ M*⊙*, mass loading factors approach similar values ($>10$) as in the lowest mass dwarfs. As a result, $\eta_{\rm M}(M_{\star})$ has a broken ‘v’ shape: the low-mass behavior regulated by stellar feedback, while the high-mass behavior is set by the energetics and coupling efficiency of BH feedback.

Even where $\eta_{\rm M}$ is established largely by the SN energy budget alone, its value and scaling with $M_{\star}$ differs from the TNG model input at the injection scale. In the lower left panel we clarify this difference, by showing the dependence on threshold $v_{\rm rad}$ (different line colors) as well as at different distances away from the galaxy (different line styles). The wind mass loading decreases strongly as a function of distance from 10 kpc outwards at $M_{\star}<10^{10.5}$ M*⊙*, dropping for low-mass galaxies in the $v_{\rm rad}>50$ km/s case from $\eta_{\rm M}\simeq 20-50$ at 10 kpc to only $\eta_{\rm M}\simeq 1$ at 50 kpc. Depending on velocity cut and radius, the actual mass outflow rate can exceed that prescribed at injection by the TNG model, presumably due to local entrainment as winds begin to interact in the outskirts of the disk. Declining $\eta_{\rm M}(r)$ is generically expected for a flow which is at least partially fountain, i.e. not entirely escaping (e.g. Sarkar et al., 2014). In Appendix A we discuss mass loading trends scaled by the virial velocity of the parent halo.

The radial dependence of $\eta_{\rm M}$ in TNG weakens towards high stellar mass, pointing towards a different launch mechanism. At $M_{\star}\sim 10^{10}$ M*⊙* we see that even for $v_{\rm rad}>150$ km/s the wind mass loading exceeds that produced by the SF-driven wind model alone. We speculate that the mutual, positive coupling between the two feedback mechanisms is roughly maximal at this mass scale - where black holes are able to contribute energy to drive outflows, but not yet severely impacting galactic star formation and so total available SN energy. As the majority of BHs transition into the low accretion state at $M_{\star}\gtrsim 10^{10.5}$ M*⊙* the radial dependence is entirely eliminated, even for fast $v_{\rm rad}>350$ km/s outflows. As we discuss later, the vast majority of gas driven by kinetic BH feedback exceeds this velocity to distances greater than $r>50$ kpc, leaving $\eta_{\rm M}$ invariant to distance in this regime.

In the lower right panel of Figure 5 we explore the redshift dependence of $\eta_{\rm M}$ in different bins of stellar mass, from $\log(M_{\star}/\rm{M}_{\odot})=10^{8}-10^{11}$. For low mass galaxies with $M_{\star}\lesssim 10^{10}$ M⊙ the wind mass loading is roughly constant as a function of redshift. This is because, although the actual outflow rates around such galaxies decrease with redshift at fixed mass, by roughly one order of magnitude from $z=6$ to $z=0.2$, their star formation rates do likewise. On the other hand, galaxies at the high-mass end exhibit a rapidly increasing $\eta_{\rm M}$ towards late times. For example, galaxies with $M_{\star}\sim 10^{10.5}$ M⊙ have $\eta_{\rm M}\sim 1$ at $z=6$, increasing to $\eta_{\rm M}\sim 100$ by $z\simeq 0$. In this case mass outflow rates are roughly constant with redshift, while star formation rates begin to decline after some characteristic redshift as a result of the quenching process. At redshift two we see the inversion of $\eta_{\rm M}$ noted above as the crossover of the high-mass curves; at $z\gtrsim 4$ it is clear that $\eta_{\rm M}$ is monotonic with $M_{\star}$, unlike at low redshift.

3.3 Outflow velocities: SN versus BH-driven winds

Moving to a second fundamental property beyond mass outflow rates, Figure 6 measures outflow velocities as a function of stellar mass. As before, we select outflowing material at some distance away from the galaxy. We also condense the full spectrum of the outflowing velocity distribution around a galaxy into a single number - mass outflow rate weighted percentiles $v_{\rm out,50}$ to $v_{\rm out,95}$ (as described in Section 2.3).222For example, $v_{\rm out,95}$ gives the velocity which only 5% of the outflow mass flux exceeds, i.e. the velocity envelope of 95% of the mass flux. In some observational studies, this value might be referred instead to as $v_{\rm out,05}$. As a characterization of the extreme tails of the outflow velocity distribution, it will approach the similar though less well defined $v_{\rm max}$ measure. In the top panel we show the 75th and 95th percentiles of $v_{\rm out}$ at a distance of 10 kpc at $z=2$. Both increase monotonically with $M_{\star}$ from $100-200$ km/s at stellar masses of $10^{7.5}$ M*⊙* up to $\sim$ 400 km/s ($v_{\rm 75}$) or $\gtrsim$ 800 km/s ($v_{\rm 95}$) for $M_{\star}\gtrsim 10^{10.5}$ M*⊙*. They are typically below the escape velocity $v_{\rm esc}=[-2\phi(r=\rm{10kpc})]^{1/2}$ at this distance. Given that the global potential also traces the large scale environment, particularly for smaller galaxies (Oppenheimer & Davé, 2008), we also include a second line for $v_{\rm esc}$ (dashed), computed using the potential difference between 10 kpc and $r_{\rm vir}$. This roughly captures, therefore, the energy required to move between these two locations, against the force of gravity. In this case the conclusion would be that $v_{\rm out}$ can commonly exceed the escape velocity, although we caution that in neither case does this actually mean that an outflow can necessarily escape the halo, as the comparison neglects hydrodynamical interactions, drag, and other similar effects. In Appendix A we discuss velocity trends scaled by the virial velocity of the parent halo.

As before, we include the prescription from the TNG model for the injection (launch) velocity $v_{\rm out}^{\,i}$ of winds generated by stellar feedback as the gray line and shaded band, which represent a fit to the mean value derived for each galaxy. We see that $v_{\rm out}^{\,i}$ represents more an upper envelope to the resulting flow at 10 kpc, rather than its characteristic velocity. We note however that $v_{\rm out,99}$ does exceed the injection velocity even at this distance. The minimum wind velocity of 350 km/s enforced by the TNG wind model leads to the constant $v_{\rm out}^{\,i}$ below $M_{\star}\lesssim 10^{9}$ M*⊙. Interestingly, we see that this does not mean that material flowing away from a galaxy at any resolved scale actually inherits such a velocity (nor even its trend with stellar mass). The $v_{\rm out,75}$ velocity percentile is actually below this ‘minimum’ for all $M_{\star}\lesssim 10^{10}$ M⊙*. We interpret this as an indication of halo drag, either gravitational or hydrodynamical (Silich & Tenorio-Tagle, 2001; Dalla Vecchia & Schaye, 2008), an important demonstration of how emergent outflow properties can differ from model inputs due to resolved physical processes.

In the lower left panel of Figure 6 we explore the dependence of $v_{\rm out}$ on distance (colors) and velocity percentile (linestyles). Outflow velocity notably decreases with distance for all $r>10$ kpc. At the same time, our velocity percentile ‘summary statistic’ will be influenced by all gas with $v_{\rm rad}>0$ km/s, i.e. it could also be pulled down by a progressively larger fraction of quasi-equilibrium, slowly moving halo gas relative to an actual, faster outflow component. Outflow velocity scales strongly with the velocity percentile, particularly at the massive end. The distributions of $v_{\rm out}$ clearly have extended tails which are not captured by $v_{\rm 50}$ for instance. At $M_{\star}=10^{11}$ M*⊙*, the 95th tail reaches $\sim$ 1000 km/s while $v_{\rm out,50}$ is only $\sim 200$ km/s. By measuring the maximal velocities achieved by gas shortly after being kicked by the low-state BH feedback, we determine a corresponding ‘injection velocity’ of $\gtrsim$ 12,000 km/s, applicable to the high-mass end, which shows how rapidly such nuclear-originating flows must decelerate.

The lower right panel of Figure 6 shows the redshift dependence of outflow speed at fixed stellar mass, where we focus on $v_{\rm out,90}$ and $r=10$ kpc for simplicity. Winds are faster at higher redshift for a given $M_{\star}$, both in the SF-driven and BH-driven regimes. In the former, a plausible cause is that, at fixed stellar mass, the DM velocity dispersion will be higher at earlier times. We note an important new aspect of the TNG wind model is the explicit scaling $v_{\rm out}^{\,i}\propto[H_{0}/H(z)]^{1/3}$. This dependence produces a co-evolution of wind velocity and virial halo mass, i.e. redshift-independent velocity at a fixed host halo mass (Pillepich et al., 2018a). After convolution with the evolving stellar mass to halo mass relation, we see that in practice the model injection velocities are only higher at a given $M_{\star}$ for very high redshift ($z\sim 6$), but settle down to a roughly constant set of the overlapping black curves for all $z<4$. In practice, this scaling is not relevant for $M_{\star}\lesssim 10^{9}$ M*⊙* due to the imposition of the minimum wind injection velocity of 350 km/s. However, we see that the resolved wind velocities evolve with redshift even at these low masses. In this regime $v_{\rm out,90}\propto(1+z)^{0.2-0.4}$ or $v_{\rm out}\propto H(z)^{0.2-0.3}$ depending on $M_{\star}$. We note that the measured scaling with $H(z)$ has the opposite sign as the model injection scaling, emphasizing that the mapping between model inputs and resolved outputs can be non-trivial. At the high mass end, the redshift scaling is much stronger, $v_{\rm out,90}\propto(1+z)^{0.8-1.3}$ for $M_{\star}\simeq 10^{11}$ M*⊙*, indicating that the velocities of BH-driven outflows evolve differently with redshift.

Having presented the population-level view of outflow velocity, we go on to explore the detailed velocity structure of outflows. Figure 7 shows the phase space diagram of gas mass in the plane of $(v_{\rm rad},\rm{radius})$ around a single halo with total mass $\simeq 10^{12.8}$ M*⊙*, from just outside the direct region of influence of the black hole (1 kpc) to about a third of the host virial radius (100 kpc). The bulk of mass at almost all radii is sitting just below $v_{\rm rad}\sim 0$ indicating consistent inflow through the inner halo. Outflowing winds are evident as a succession of discrete features at positive radial velocities. These outflows launched by black hole feedback form a sawtooth-like pattern in $(v_{\rm rad},r)$, consisting of a series of triangular components – each results from a single energy injection event at a small scale off the left-edge of this panel. There are three such flows evident at $r<20$ kpc, and $\sim$ 5 visible in the inner halo.

The maximum radial velocity of each outflow component increases with distance, i.e $v_{\rm out}\propto r^{\,\alpha}\,(\alpha>0)$. This shape is a natural consequence if material has a spectrum of launch velocities and the most rapidly outflowing gas reaches, at any given time, the largest distances. The maximal outflow velocity at each radius is established by the upper envelope of the successive peaks, corresponding to the highest velocity tails of each flow. Over time, each individual feature collapses outward, shifting to larger distance and lower velocity. As a result, the envelope of outflowing radial velocities actually decreases with distance, i.e $v_{\rm out,max}\propto r^{\,\alpha}\,(\alpha<0)$. This velocity structure, coming from the superposition of the individual outflow components, is continually reinforced by new injection events at $r\sim 0$ and generically maintains a quasi-steady structure over time. The high velocity tails of the outflowing material lie above the escape velocity curve $v_{\rm esc}(r)=[-2\phi(r)]^{1/2}$ at all radii (black line).

Figure 8 presents the distribution functions of outflow velocity, decomposing the mass outflow rate $\dot{M}_{\rm out}$ as a function of $v_{\rm rad}$. We stack galaxies in bins of stellar mass, from $8.0<\log(M_{\star}/\rm{M}_{\odot})<11.5$ (different line colors), and focus on a distance of $r=20$ kpc at $z=2$ for simplicity. At lower stellar masses the distributions have a similar shape, being peaked at low velocities $\lesssim$ a few hundred km/s, which increases slowly with $M_{\star}$. The amount of outflowing gas with higher velocity drops rapidly beyond this central core, such that there is an effective maximum outflow velocity of $\sim$ 300 km/s for $M_{\star}=10^{8}$ M*⊙, $\sim$ 400 kms for $M_{\star}=10^{9}$ M⊙, and $\sim$ 500 km/s for $M_{\star}=10^{10}$ M⊙*. This is the regime of the stellar feedback driven wind model, with the approximate scaling of $v_{\rm out}\propto M_{\rm halo}^{1/3}$. As a result, typical wind velocities scale up with mass, but have bounded maximal values.

However, at stellar masses $M_{\star}\gtrsim 10^{10.5}$ M*⊙* a high-velocity tail begins to emerge. By $10^{11}$ M*⊙* this feature is fully developed and results in a roughly linear extension out to high $v_{\rm rad}$ (in this log-linear plane). As a result, we can start to find outflowing mass at $>$ 1000 km/s velocities, where the distributions are then composed of two distinct components: the SF-driven core at low velocity, and the tail at high velocity. In the TNG model, only the kinetic BH feedback mechanism can generate such high velocities, and this mode becomes efficient at $M_{\star}\gtrsim 10^{10.5}$ M*⊙* as black holes transition below the threshold accretion rate. The mass outflow rate decreases more slowly with $v_{\rm rad}$ than in the SF-driven wind case, such that there is no strong cutoff at some maximal velocity. However, the decline is still precipitous: for the highest mass bin considered here at $M_{\star}=10^{11.5}$ M*⊙*, for $v_{\rm out}\sim 1000\pm 50$ km/s, the total flux is on average $\dot{M}_{\rm out}\simeq 10\,\rm{M}_{\odot}\,\rm{yr}^{-1}$. However, for gas moving with $v_{\rm out}\sim 2000\pm 50$ we have only $\dot{M}_{\rm out}\simeq 1\,\rm{M}_{\odot}\,\rm{yr}^{-1}$, and the flux of gas moving at $\sim$ 3000 km/s is a further factor of ten less.

As we discuss further in Section 3.6, it is important to note that these highest velocity BH-driven outflows in the TNG model do not arise from ultra-luminous or ‘radio-loud’ quasars (e.g. with space densities $n\sim 10^{-7}$ Mpc*-3* or less, and/or $L_{\rm bol}\sim 10^{45}$ erg/s or greater), which are much too rare to arise in the volume of TNG50. Instead, they are generated from the much more common low luminosity, low-accretion rate ($\lambda_{\rm edd}\sim 10^{-2}$ or less) population of BHs for which we posit a radiatively inefficient flow which converts gravitational binding energy into a non-relativistic wind (Blandford & Begelman, 1999; Yuan & Narayan, 2014). The implication is that – in the TNG model – low luminosity, slowly accreting black holes can drive some of the most powerful outflows. Other models produce strong BH-driven outflows by invoking different physical mechanisms and/or numerical implementations (e.g. McCarthy et al., 2011; Henden et al., 2018), and differentiating between these scenarios will provide useful theoretical constraints.

3.4 Multiphase outflows: temperature, density, metallicity

We have so far considered only the kinematics of the outflow as a whole, i.e. phase-agnostic. Observationally there exists a broad correlation between outflow velocity and gas phase temperature, such that neutral phases are slower than ionized components (Heckman et al., 2001; Rupke et al., 2002), which are themselves slower than a hot wind fluid phase (Veilleux et al., 2005).

In Figure 9 we show that different outflow phases in TNG have different kinematics. We present the galaxy-stacked 2D histogram of mass outflow rate, as a function of temperature (x-axis) and instantaneous velocity (y-axis). We focus on a mass scale ($M_{\star}=10^{11}$ M*⊙*) where the action of the central black hole plays an important role in setting the properties of the galactic-scale outflows. For winds which reach an appreciable distance away from the galaxy we find a relative dichotomy of hot/cooler gas phases and rapid/slower outflow velocities. In particular, the only outflowing gas which achieves velocities of $>$ 1000 km/s is hot, with $T\gtrsim 10^{6}$ K. Within this phase there is a minor correlation between $v_{\rm out}$ and temperature: the hottest gas travels the fastest, with $v_{\rm out}\sim 2000$ km/s outflows having temperatures of order $10^{6.5}$ K.

On the other hand, cooler gas is slower, and so unlikely to propagate as far. The outflow component which would be visible in ionized H$\alpha$ or metal transition lines, between $10^{4}$ K and $10^{5}$ K, has typical velocities of a few hundred km/s and a maximal $v_{\rm out}$ of $\sim$ 600 km/s. The phase separation at high velocity implies that the temperature of a high-velocity component can be clearly deduced. The bulk of the outflowing gas, however, coexists at a common, lower velocity interval, preventing one from easily identifying the phase structure of this dominant, low-velocity component. We note that for lower stellar masses $M_{\star}\sim 10^{9-10}$ M*⊙* (not shown), where stellar rather than BH-driven winds are most relevant, this correlation of hotter outflow phases moving faster likewise holds.

Figure 10 quantifies the multiphase nature of TNG outflows in terms of temperature, density, and metallicity. In the main panel (top) we decompose the mass outflow rate as a function of gas temperature, where each line shows the mean stack within a given stellar mass bin, from $8.0<\log(M_{\star}/\rm{M}_{\odot})<11.5$ (different line colors). Two important trends emerge.

First, there exists a warm/hot outflowing component with a well-behaved, roughly gaussian distribution of mass flux weighted log temperature. The peak of this component shifts smoothly upwards with increasing galaxy stellar mass, from $\sim 10^{5.5}$ K at $M_{\star}=10^{8}$ M*⊙* to $\sim 10^{6}$ K by $10^{9.5}$ M*⊙, and all the way up to $\sim 10^{7}$ K at $10^{11}$ M⊙*. The width (in log K) is roughly constant, implying that this constituent of the outflow is set by the temperature structure of the virialized hot halo gas and tracks the increasing host halo virial temperature. Note that bulk (i.e. non-circular) motions in the inner halo gas which exceed the chosen velocity threshold would contribute here, so the hot outflow component could partially reflect increasing amounts of such non-equilibrium kinematics of virialized gas. The peak outflow rate of this component increases by roughly a factor of ten with each decade in stellar mass – i.e., it is a roughly constant ‘specific gas outflow rate’. In the integral it dominates, at all mass scales, the total outflowing mass flux.

Next, a second cool outflow component emerges towards higher stellar masses, resulting in a bimodal outflow temperature distribution. This bimodality is clear at all $M_{\star}>10^{10}$ M*⊙* and becomes stronger towards higher mass, where the peak outflow rates at $T\sim 10^{4}$ K and $T\sim 10^{7}$ K come into near equality.333Some outflowing gas at temperatures even lower than $10^{4}$ K is also evident. While it would be tempting to interpret this in the context of material which could be neutral or even partially molecular, we remind that the effective cooling floor of TNG is $\simeq$ 104 K due to the radiative cooling assumptions, which neglect metal fine-structure contributions. As a result, we interpret this gas as representative only of a cool ionized component. The origin is twofold - the hint of this secondary peak already developing even for galaxies as small as $M_{\star}=10^{9}$ M*⊙, together with the invariance of the peak temperature itself in relation to the shape of our adopted cooling curve (Wiersma et al., 2009) suggests a cooling process (possibly in the outflow itself, e.g. Thompson et al., 2016; Schneider et al., 2018). In the future we can use the Monte Carlo tracer particle scheme in TNG (Genel et al., 2013) to directly track the thermal history of outflowing gas and (dis)prove this point. On the other hand, the strong increase of outflowing mass at this temperature above $M_{\star}>10^{10.5}$ M⊙* indicates a relation to the onset of efficient black hole feedback. We observe that the kinetic BH mechanism in this regime directly evacuates the central ISM reservoir, leading to centrally suppressed gas densities (and star formation rates, as we discuss later). At stellar masses $\gtrsim 10^{11}$ M*⊙* this dichotomy in outflow temperature is therefore a direct consequence of ejective feedback launching formerly dense ISM out into the halo (e.g. Scannapieco & Brüggen, 2015).

Observationally, Chisholm et al. (2017) recently measured the mass loading of OVI absorbing gas in a single lensed $z\sim 3$ galaxy, concluding that the fastest outflowing material (with $v_{\rm out}>v_{\rm esc}$) has a mass outflow rate larger than in the cooler, photoionized phase. Similarly, from a sample of local starbursts Grimes et al. (2009) find that OVI-traced gas at $\sim 10^{5.5}$K has a higher mean outflow velocity than seen in lower ionization states. If the expectation presented here is roughly correct, then the most important outflow phase for $M_{\star}\gtrsim 10^{10}$ M*⊙* moves quickly away from available ionized tracers and even past OVI bearing ‘coronal’ gas into a much hotter x-ray traced plasma, where observational constraints are presently challenging.

In the lower left panel of Figure 10 we similarly decompose the mass outflow rate now as a function of gas density, in the same stacked distributions as a function of stellar mass. A two component structure is also evident: the low density peak is associated to the high temperature gas, with the high density peak corresponding to the low temperature component which prominently emerges only at $M_{\star}>10^{10}$ M*⊙. As before, the location of the lower density component tracks the increasing halo gas density with higher host mass. On the other hand, the high density ISM component has a roughly constant maximal outflow rate for densities of $n\sim 1$ cm-3*, although the most massive galaxies above $M_{\star}>10^{11}$ M*⊙* can launch outflows to 20 kpc which contain densities of $n>10$ cm*-3* and higher. This also indicates a direct ejective origin for this dense material, having originally compressed in the center of the galactic potential. Given our star formation threshold density, some of this ejecta will therefore be star-forming and on the effective equation of state.

The lower right panel of Figure 10 measures the metallicity of outflowing gas, in three stellar mass bins (colored lines) at $z=1$, $z=2$, and $z=4$ (solid, dashed, dotted). At all redshifts, the outflow metallicity increases with the stellar mass of the galaxy as expected. At redshift two, the dominant metallicity of wind material is $0.1Z_{\odot}$ for $M_{\star}=10^{8}$ M*⊙, increasing to $0.3Z_{\odot}$ at $M_{\star}=10^{9.5}$ M⊙* and reaching $0.7Z_{\odot}$ by $M_{\star}=10^{10.5}$ M*⊙*. The distributions are broad and gas at $\pm 0.5$ dex of these values is always also present at roughly 10% of the mass outflow rate when compared to the dominant metallicity. At low masses the redshift evolution is roughly independent of stellar mass and fairly weak – the average $Z_{\rm out}$ increases by $\sim$ 0.2 dex from $z=4$ to $z=2$ and by a further $\sim$ 0.1 dex to $z=1$. At high masses the redshift evolution is even shallower due to the flattening of the gas-phase mass metallicity relation. In both cases the evolution of outflow metallicity is approximately similar to the evolving MZR of TNG (Torrey et al., 2017), although we defer to future work a detailed comparison of wind metallicity relative to the galaxy ISM metallicity itself (e.g. Muratov et al., 2017).

3.5 Angular dependence: emergent bipolarity

In Figure 11 we show the dependency of the mass outflow rate (as the background color) on both galactocentric angle and galactocentric distance at $z=1$. Within each panel the two horizontal lines indicate the direction aligned with the minor axis of the galaxy.444We rotate each galaxy to an edge-on orientation, such that the longest axis is aligned with $\hat{x}$ (i.e. the $y=0$ line) when viewed in projection in the $x-y$ plane, and the angle $\theta$ for every gas cell is calculated from its $(x,y)$ coordinates in this projection, such that $\theta=\{0,\pm\pi\}$ corresponds to the $\pm x$ (major) axes, respectively, while $\theta=\pm\pi/2$ corresponds to alignment along the $\pm y$ (i.e. minor) axes, respectively. The top panel stacks relatively low-mass galaxies with $M_{\star}\simeq 10^{10}$ M*⊙* ($\sim$ 280 systems), a mass regime where outflows are generated by stellar feedback. The bottom panel stacks high-mass galaxies with $M_{\star}\simeq 10^{11}$ M*⊙* ($\sim$ 70 systems), where higher-velocity outflows result almost exclusively from BH feedback.

In both cases the conclusion is the same: the mass outflow rate of winds is not directionally isotropic. Rather, strong $\dot{M}$ outflows are preferentially found aligned with the minor axes of galaxies. This is particularly intriguing in the context of TNG, as it represents an emergent quality of galactic-scale outflows. In the TNG model, wind launching and energy injection from stellar as well as black hole feedback is entirely isotropic, in the time-average, multiple event sense. For stellar feedback driven winds, the directionality of the wind-phase mass is entirely random at launch (see Pillepich et al., 2018a), while BH energy injection in the low accretion state kinetic mode, which drives the outflows seen here at high-mass, is also entirely random with each injection event (see Weinberger et al., 2017). The preferential propagation of outflows along the minor axis of a galaxy therefore represents a natural, hydrodynamical collimation of the flow. White lines demarcate the ‘half mass outflow rate angle’, quantifying a measure of the opening angle; this is $\theta_{\rm 1/2}\simeq 70^{\circ}$ at $r\sim 10$ kpc and $\theta_{\rm 1/2}\simeq 40^{\circ}-50^{\circ}$ at $r\sim 50$ kpc.

Figure 11 quantifies the visual impression provided by Figures 2 and 3 for the case of a strong AGN driven outflow. To complement this high mass, high outflow velocity regime, we visualize the velocity structure of a smaller galaxy with $M_{\star}=10^{10.3}$ M*⊙* ($M_{\rm halo}=10^{11.8}$ M*⊙*) in Figure 12. Such a system observed at a higher redshift would be similar to a Lyman-break galaxy. Here we show a gas density projection of the galactic and halo gas, alloyed with its instantaneous velocity field using the line-integral convolution technique. Small-scale transport and swirling, vortical motions are revealed in regions without strong bulk outflow or inflow. Overlaid on top, streamlines indicate flow direction and velocity, with outflows reaching up to $\sim$ 350 km/s as they establish a circulatory, galactic fountain like baryon cycle across $r_{\rm vir}/2$ (dotted circle).555In addition to this one halo, a comprehensive set of examples, spanning a range of galaxy masses and redshifts, is available on the TNG50 website at www.tng-project.org/explore/gallery/.

We speculate that this collimation is largely a ‘path of least resistance’ effect, whereby outflows directed into the plane of a disk, comprised of the densest material, experience enhanced resistance. Consequently, we find that a roughly bipolar or biconical outflow is typical behavior. Sculpting of the flow is evident already by $\sim 10$ kpc in the inner halo, and persists out as far as strong outflows are found ($\gtrsim 50$ kpc at least). The inferred opening angle is large, decreasing for stronger outflow thresholds which become more concentrated along the minor axes. We note that the degree of this collimation appears to evolve with redshift (not shown), becoming stronger at late times and emerging particularly from $z=2$ to $z=1$. It is difficult to identify such a clear signal even at $z=2$. The degree of collimation is likely to depend on galaxy morphology, particularly the amount of rotational support of the system.

Observationally, preferential detection of outflow signatures along the minor axes of galaxies is a common conclusion at low-z (Heckman et al., 2000) based on SDSS (Chen et al., 2010; Concas et al., 2017; Bae & Woo, 2018) and also towards $z\sim 1$ (Martin et al., 2012; Kornei et al., 2012; Bordoloi et al., 2014b; Rubin et al., 2014). At the same time, data from even higher redshifts $z\geq 2$ support a lack of collimation at such early times (Law et al., 2012). Typical evidence along these lines comes from cold gas-phase tracers such as MgII or NaD, or hydrogen, and the association is typically with stellar feedback driven outflows, although asymmetric outflows from low-luminosity black holes may also be common (Cheung et al., 2016; Wylezalek et al., 2017; Roy et al., 2018). For hotter phases, Kacprzak et al. (2015) find OVI around $z\leq 0.7$ galaxies preferentially along either the major or minor axes (e.g. inflows and outflows, respectively), but not in between; similarly so for colder MgII absorbers at $z\geq 0.3$ (Nielsen et al., 2015).

In addition to spanning both feedback regimes, we also find no clear temperature dependence of this signal. We therefore expect (i) hotter outflowing gas phases to be similarly collimated, (ii) bipolar geometry outflows to be present across a broad stellar mass range, and (iii) a lack of similar signatures at high redshift ($z>2$; not shown). We speculate that this collimation arises only towards lower redshift due to the rise of ordered rotation and settling of gaseous galactic disks which helps to sculpt the flows.

Theoretically such a natural collimation was anticipated already by Tomisaka & Ikeuchi (1988) whose numerical simulations demonstrated galactic-scale superwind emergence preferentially elongated along the direction of the minor axis of a launching disk, caused by propagation in the direction of the maximum pressure gradient of the confining gaseous halo. Early wind models for global simulations sometimes including such hard-coded directionality by hand (Springel & Hernquist, 2003), including in the original Illustris model (Vogelsberger et al., 2013), but this was removed in TNG mainly as an opportunity to simplify an unnecessarily complex model feature (Pillepich et al., 2018a).

3.6 Relating outflow and host galaxy properties

With the exception of global stellar mass trends, we have not yet explored the dependence of outflows on additional properties of the galaxies from which they are launched. Here we specifically consider star formation activity (SFR) as well as black hole activity ($L_{\rm bol}$, $\lambda_{\rm edd}$). As both types of feedback activity are a strong function of galaxy mass, we naively expect a strong scaling of most outflow properties such as $v_{\rm out}$ with $M_{\star}$, and move beyond such zeroth order expectations by predominantly investigating trends at fixed mass.

In Figure 13 we consider the dependence of outflow velocity on the star formation rate of a galaxy, cast in terms of $\Delta$SFMS, the deviation from the median star-forming main sequence relation. To isolate out the general trend of increasing velocity with $M_{\star}$, we further normalize $v_{\rm out}$ by its median value in each stellar mass bin. Color then indicates higher (red), lower (blue), or typical (light gray) outflow velocity relative to galaxies at that same mass. As our operating definition of outflow velocity we take $v_{\rm out,50,<5kpc}$, the 50th velocity percentile of outflowing mass near the galaxy. The trends we discuss do not depend strongly on this choice.

First, we find a clear correlation between $v_{\rm out}$ and $\Delta$SFMS for star-forming galaxies, whereby outliers above the main sequence launch faster winds, whereas outliers below the MS drive slower outflows. This statement is true at fixed stellar mass, and holds from $z=1$ up to $z=6$ if not higher.

However, this correlation actually inverts at $M_{\star}\gtrsim 10^{10.5}$ M*⊙*. Galaxies which have grown above this mass scale drive faster outflows if they are below the main sequence, i.e. if they are in the process of quenching towards the red population. In the main panel, this is visible as an excess of red pixels for low $\Delta$SFMS at high $M_{\star}$. We quantify this anti-correlation in the inset panel, which plots relative $v_{\rm out}$ as a function of $\Delta$SFMS for high-mass galaxies (orange) versus low-mass star forming galaxies (blue). The two slopes have opposite sign, and the strength of this relation is such that $v_{\rm out}$ varies by $30-40$% for outliers offset by $\Delta$SFMS = $\pm 0.5$ dex.

Observationally, Chisholm et al. (2015) correlate outflow velocity and SFR in nearby star-forming galaxies, finding that mergers identified via optical morphology drive faster outflows, and that merging systems have the fastest outflows overall. If mergers are preferentially found at $\Delta$SFMS $>$ 0, this may be a similar signal as discussed here. Cicone et al. (2016) explicitly compare outflow velocity measurements with $\Delta$SFMS (their Fig. 21) for a ‘normal’ SDSS star-forming population, finding a clear positive trend above the main-sequence, with a weak or flat trend below. Towards higher masses, Sato et al. (2009) use a SF-agnostic sample from EGS to uncover a continuation of detectable NaI outflows into the red population, concluding that outflows outlive the star-formation phase and contribute to quenching galaxies en route to the quiescent population. Our analysis here implies that these ‘relic’ outflows may also be mixed with more recent outflows driven by black holes rather than supernovae.

We note that for star-forming galaxies a large number of properties correlate strongly with $\Delta$SFMS. That is, (a) outliers above the main sequence could be outliers because they drive faster outflows, or (b) they could simply have faster winds because they are already part of the $\Delta$SFMS $>$ 0 population for another reason. In this regime, relative to those systems below the median main-sequence at fixed stellar mass, galaxies have: faster outflows (as noted), lower gas-phase metallicities, higher gas fractions, lower central stellar surface densities, lower stellar metallicities, bluer (g-r) colors, larger stellar sizes (half mass radii), lower outflow mass loadings, less rotationally supported gas disks, higher halo masses, and lower mass BHs (all not shown, but see \textcolorbluePillepich et al. 2019). Oscillation and quasi-equilibria about the evolving SFMS leads to an interrelation between many of these galaxy properties. Although the root cause of the trend of $v_{\rm out}$ with $\Delta$SFMS is difficult to determine, and to some degree both (a) and (b) must apply, we note that the injection scaling of $v_{\rm out}\propto\sigma_{\rm DM}$ in the TNG model is sensitive to a spatially localized measurement of the dark matter velocity dispersion, which would naturally respond to merger activity (i.e. where $\sigma_{\rm gas}$ would also be enhanced) in contrast to more isolated, dynamically quiet environments. Likewise, the positive correlation of $M_{\rm halo}$ with $\Delta$SFMS implies larger $\sigma_{\rm DM}$ at fixed $M_{\star}$ (similar to the finding of Matthee et al., 2017, with the EAGLE model).

In order to reincorporate the dominant trend of increasing $v_{\rm out}$ with stellar mass, Figure 14 presents a view on the fraction of galaxies with ‘fast’ outflows, with respect to their position on the $\Delta$SFMS$-M_{\star}$ plane. We do this by defining $v_{\rm out,90,10kpc}>300$ km/s as a somewhat arbitrary though useful threshold. For example, if observational detectability of an outflow depended mainly on the existence of a component traveling above 200 km/s (i.e. above an instrumental limit or intrinsic level), our choice would map to an outflow detection fraction. The color of each pixel shows the fraction of galaxies in that bin which satisfy the velocity cut, from none (dark purple) to 100% (yellow). As a common feature, we find a diagonal boundary across the $\Delta$SFMS$-M_{\star}$ plane, with fast outflow galaxies residing towards the upper right. This can be understood as the superposition of the $v_{\rm out}-\Delta$SFMS secondary correlation on top of the dominant $v_{\rm out}-M_{\star}$ relation. A higher threshold velocity moves this boundary to the right, with roughly constant slope.

This trend can be compared in spirit to Fig. 17 of Rubin et al. (2014) at $z\sim 0.5$ based on MgII absorption modeling, or Fig. 3 of Förster Schreiber et al. (2018a) at $z=0.6-2.7$ based on broad emission line decomposition. Both hint towards a tilted trend in $\Delta$SFMS at fixed $M_{\star}$, in the sense that galaxies above the main sequence have a larger detected outflow fraction, while galaxies approaching $10^{11}$ M*⊙* and higher, even those below the extension of the main-sequence, do likewise. Similarly, Fig. 18 of Cicone et al. (2016) based on OIII emission from ionized outflows at $z\lesssim 0.7$ also highlights higher outflow velocities at positive $\Delta$SFMS. We emphasize such similarities are largely speculative and hindered by the current inability to make quantitative comparisons to these observables, as discussed below in Section 4.2.

Given our current inability to do a full forward-modeling, and the uncertainties inherent with inversion in the opposite direction starting from the observational quantities, our ability to make a robust comparison between the simulations and observations is presently limited. However, we would like to give a sense of the overlap between the outflow properties expected from the TNG model and available observational datasets. We therefore provide a series of suggestive plots, which are intended to be qualitative only in nature, and from which no direct statements can be made as to the level of (dis)agreement between the simulations and observations. Further modeling efforts from the theory side are needed before we can contrast against, or provide interpretation of, the wealth of high-quality observational data on outflows (see Rupke, 2018).

In Figure 15 we show a number of correlations between outflow properties and central galaxy or supermassive black hole properties. In these various spaces we add a diverse set of observations, with citations indicated in the caption. These span all redshifts, at least as high as $z\sim 4$, and include many local universe data points from $z\sim 0$. They also span a range of sample selections – in particular, a significant number specifically target some of the most energetic nearby ULIRGs or the most luminous observable starbursts. All detectable gas phases are represented, including molecular, neutral, and ionized tracers. Importantly, we combine datasets where outflows are believed to arise from stellar feedback as well as BH feedback – their separation being often difficult.

In addition, outflow properties including velocity and mass loading factors are measured in a diversity of ways and with a number of definitions across the observational works, and no homogenization has been undertaken. For example, outflow velocities may be derived from line center offsets from systemic, line wings below some percentile of the continuum, or based in some way on line width, with or without inclination corrections. Error bars are excluded in every case for clarity.

All simulation results are restricted to a consistent minimum stellar mass of $M_{\star}\geq 10^{9}$ M*⊙*, but otherwise include galaxies of all stellar masses as representative of a mass complete volume-limited sample, and are shown at $z=1$ (colored lines). To highlight the definitional range, we include five velocity percentile ranges ($v_{25}$ to $v_{99}$) for $v_{\rm out}$, or three velocity cut criteria (0 - 350 km/s) for $\eta_{\rm M}$, as appropriate. Measurements are made at $r=10$ kpc; varying galactocentric distance would produce an additional fan of lines in each panel. Observational data points are shown as black symbols.

• •

The upper left panel shows the relationship between $v_{\rm out}$ and star formation rate of the central galaxy. The stellar mass of galaxies in this plane increases monotonically from left to right, while outliers which scatter towards high $v_{\rm out}$ at a given SFR are massive, low-sSFR, and black hole outflow dominated. Although outflow velocity increases with SFR, the slope is fairly shallow: from $0.1\,\rm{M}_{\odot}\,\rm{yr}^{-1}$ to $100\,\rm{M}_{\odot}\,\rm{yr}^{-1}$ the outflow velocity increases by a factor of $\sim 3-5$. Observational datasets which span only a small dynamic range in SFR are unlikely to robustly identify a trend and would conclude a flat $v_{\rm out}$ (as appreciated in e.g. Chen et al., 2010). We highlight the results of Cicone et al. (2016), where we have taken the LoSVD $v_{\rm 0.1}$ (OIII) values based on a consistent analysis of a large, low-z sample which spans nearly our entire SFR range and shows an encouraging similar scaling as $v_{\rm out,95}$ for instance. We note that due to the finite volume of TNG50, there is only a partial sampling of high $\Delta$SFMS outliers and very few galaxies with star formation rates above 100 $\rm{M}_{\odot}\,\rm{yr}^{-1}$.

• •

With the upper right panel we show the dependence on the mass loading factor $\eta_{\rm M}$ with the SFR of the galaxy. As above, the stellar mass of galaxies in this plane increases monotonically from left to right, while outliers which scatter towards high $\eta_{\rm M}$ at a given SFR are quenching (or quenched) and black hole outflow dominated. Note that many observations present $\eta_{\rm M}$ of a specific gas phase which is only a sub-component of the total. Many of these (e.g. Chisholm et al., 2015) imply to zeroth order a trend of decreasing $\eta_{\rm M}$ with SFR, but as we have shown in Section 3.4 higher SFR (and so $M_{\star}$) galaxies progressively shift the majority of their outflowing mass into a difficult to observe hot phase. Without a phase-specific modeling of $\eta_{\rm M}$ we cannot robustly compare to trends in SFR, other than to comment that a large fraction of observed points, which are predominantly at $1.0<\log(\rm{SFR}/\rm{M}_{\odot}\,yr^{-1})<2.5$ have lower inferred mass outflow rates than the phase-total rates of TNG.

• •

In the center left panel we show the dependence of outflow velocity on the bolometric luminosity of the black hole. As with extreme starbursts, the volume of TNG50 is also too small to host the brightest luminosity AGN. There are essentially no black holes with $L_{\rm bol}>10^{45}$ erg/s, while the majority of observations of BH-driven outflows focus on this regime of bright quasars which host $1000$ km/s or greater winds, which are offset from our simulated sample. When observations of ‘normal’ galaxies with black holes exist, inferred values for $v_{\rm out}$ cover the range of velocity percentiles from TNG (Fluetsch et al., 2018, where $v_{\rm out}$ = FWMHbroad/2 + $|v_{\rm broad}-v_{\rm narrow}|$). On this plane, stellar mass increases with $L_{\rm bol}$ from left to right, but also with $v_{\rm out}$ from bottom to top. In particular, the tails visible at low black hole luminosities are the most active BH-driven outflows in the simulation.

• •

The center right panel shows $\eta_{\rm M}$ as a function of $L_{\rm bol}$ of the central black hole. Depending on velocity threshold, the trend can be slightly decreasing, flat, or increasing with $L_{\rm bol}$, but is weak in any case. In comparing to observed trends, the previous caveats about phase specific mass outflow rate tracers apply. In TNG, the scatter in mass loading increases strongly below $L_{\rm bol}<10^{42}$ erg/s, where low and high mass galaxies are simultaneously present, with the latter driving high $\eta_{\rm M}$ outflows.

• •

The lower left panel highlights the relation between $\eta_{\rm M}$ and $\Sigma_{\rm SFR}$, an integrated star formation rate surface density, computed within one times the stellar half mass radius (roughly similar to $R_{\rm e}$). As before, depending on velocity threshold, the trend of mass loading with galaxy-integral $\Sigma_{\rm SFR}$ can be negative, flat, or even possibly positive. At high-z and based on kpc-scale correlations, Davies et al. (2018) find an increasing trend, suggestive of a small-scale relationship which warrants an in-depth comparison.

• •

Finally, the lower right panel shows the relationship between outflow velocity and mass loading factor. The usual scalings of $\eta_{\rm M}\propto v_{\rm out}^{-1}$ (’momentum’ driven) and $\eta_{\rm M}\propto v_{\rm out}^{-2}$ (’energy’ driven) are indicated by the dashed gray lines. The latter is prescribed by the TNG model at injection (Pillepich et al., 2018a), and we roughly recover this expectation between $\eta_{\rm M}$ of a few and a few tens. In this range, stellar mass increases leftwards, from $M_{\star}\sim 10^{8}$ M*⊙* at $\eta_{\rm M}$ of a few tens to $M_{\star}\sim 10^{10.5}$ M*⊙* at $\eta_{\rm M}$ $\sim$ a few. High mass, quenching or quenched systems populate the top of this plane, distributed across all values of $\eta_{\rm M}$ as high outliers to the median lines. For example, at $v_{\rm out}\sim 1000$ km/s and $\eta_{\rm M}\sim 50$. Many of the observational data points in this regime, for instance of the Fiore et al. (2017) compilation, are thought to be black hole driven outflows. The upturn of outflow velocity at $\eta_{\rm M}\gtrsim 10$ occurs when low-mass galaxies disappear and $M_{\star}>10^{10.5}$ M*⊙* systems set the median. The monotonic scaling behavior of $v_{\rm out}(\eta_{\rm M})$ breaks down when low-SFR galaxies begin to dominate at a given mass loading.

Outflow properties are clearly related to galaxy stellar mass, star formation rate, and black hole activity, as well as the distribution of star formation activity. To emphasize the co-existence of stellar driven outflows across the mass range where BH activity also becomes relevant, Figure 16 presents stacked radial profiles of star formation rate surface density $\Sigma_{\rm SFR}(r)$ at $z=1$ for galaxies of different stellar masses, from $10^{8}$ to $10^{11.5}$ M*⊙*.

Overall, the star formation scales up with increasing $M_{\star}$ at all radii. Except for the two highest mass bins, $\Sigma_{\rm SFR}$ profiles are always monotonically decreasing away from the galactic center. They have a characteristic break at $\sim 1-2\,r_{\rm 1/2,\star}$, outside of which profiles steepen noticeably. We draw attention to the behavior within the galaxy body itself at high-mass, where $\Sigma_{\rm SFR}$ flattens at $r<10$ kpc for $M_{\star}\gtrsim 10^{11}$ M*⊙*. For the most massive systems, the profile even inverts towards the center. This central inversion is also present in the total gas density profile as well as the gas-phase metallicity profiles (not shown), indicative of dilution of the remaining central gas content by metal-poor inflows.

The drop in star formation (and so gas mass) in the central ISM is another example of the ejective character of the low-accretion mode BH feedback in TNG, which starts to operate at $M_{\star}\gtrsim 10^{10.5}$ M*⊙*. We previously saw evidence for this in Figure 10 as the culprit for the cold phase of the bimodal temperature distribution of outflowing gas. Even though such galaxies are either in the progress of quenching for the first time, or maintaining their quenched state, Figure 16 shows that there can be non-negligible, residual star formation, particularly at characteristic larger radii.

4 Discussion

4.1 The dichotomy of SN versus BH-driven outflows and the transition to quiescence

As any star-forming gas is eligible to drive a wind-phase in TNG, galaxies do not by construction require a threshold in $\Sigma_{\rm SFR}$ such as the canonical 0.1 M*⊙* yr*-1* kpc*-2* suggested in Heckman (2002) as a requirement to launch a galactic wind (but see Rubin et al., 2014). In our case, outflowing mass flux will correlate with the spatial locations of star formation, where the most supernova energy is available. The suppression of central SF at high $M_{\star}$ implies that stellar feedback driven winds will then become less nuclear and more extended. At the same time, the co-existence of SF-driven and BH-driven outflows also becomes inevitable at $M_{\star}\gtrsim 10^{10.5}$ M*⊙* (for example, as seen in NGC 3079; Cecil et al., 2001).

In general, energy inputs from stellar and black hole feedback are essentially always both present, except in the case where the halo is not yet massive enough to host a BH seed. That is, one cannot absolutely separate out the energetic contribution from these two sources to the production of an outflow, neither in simulations nor in observations. If part or all of a measured $\dot{M}_{\rm out}$ is generated due to energy injection from black holes rather then supernovae, then the denominator in our expression for $\eta_{\rm M}^{\rm SN}=\dot{M}_{\rm out}/\dot{M}_{\star}$ (Eqn. 4) neglects some if not most of the relevant energetics, as implied by the ‘SN’ superscript.

It would be tempting to write $\eta_{\rm M}^{\rm tot}=\dot{M}_{\rm out}/(\dot{M}_{\star}+\dot{M}_{\rm BH})$ or $\eta_{\rm E}^{\rm tot}=\dot{E}_{\rm out}/(\dot{E}_{\rm SN}+\dot{E}_{\rm BH})$ and so measure the outflow properties relative to the total available energetics. However, SN versus BH feedback operate at different (spatial) scales, and have different effective (galactic or halo-scale) coupling efficiencies. For instance, although $\dot{E}_{\rm BH}$ might dominate $\dot{E}_{\rm SN}$ for some main-sequence galaxies with rapidly accreting black holes, the resulting thermal AGN feedback may have little dynamical impact and so not actually be the cause of any outflows. Conversely, even a quenching or almost entirely quiescent galaxy may have $\dot{E}_{\rm SN}\gg\dot{E}_{\rm BH}$; nonetheless, the relatively small amount of energy available from the low-accretion state BH feedback could in practice be solely responsible for a high mass flux, high velocity outflow, possibly due to the spatial (i.e. nuclear) or temporal (i.e. rapid) coherence of that energy injection, or due to the physical/numerical coupling mechanism.

These differences leave signatures in properties of contemporaneous outflows. In TNG, the process of reddening and migration away from the blue, star-forming population is tied closely to a threshold in black hole and so stellar mass (Nelson et al., 2018a). As a result, the galactic-scale winds which arise from galaxies in different locations of the SFR$-M_{\star}$ plane (and with different BH properties) directly constrain underlying assumptions of the feedback models. Emblematically, Figure 13 reveals the correlation of $v_{\rm out}$ with $\Delta$SFMS as well as its inversion into an anti-correlation above a particular mass scale. Feedback also leaves signatures in the residual star formation of the galaxy itself – Figure 16 demonstrates that the evolutionary process of quenching proceeds ‘inside-out’ in TNG, with star-formation truncated in the center of the galaxy first. IFU survey data at low redshift will unambiguously measure the spatially resolved progression of SFR through the transition to quiescence (Belfiore et al., 2017), providing stringent and complementary constraints to the galactic outflows themselves.

Observationally it has been noted that, at least in the local universe, ‘normal’ AGN-host galaxies do not show a clear difference in outflow prevalence or properties when compared to non-AGN control samples (Sarzi et al., 2016; Roberts-Borsani & Saintonge, 2018). This may come down to the timescale issue, given that the BH need not be ‘on’ and injecting energy at the moment of observation, although a recent (i.e. Myr-scale) outflow could still be present (Förster Schreiber et al., 2014; Gabor & Bournaud, 2014). As shown in Figures 2 and 3, Myr-timescale variability at the energy injection scale can produce a coherent, large-scale BH-driven outflow which lasts for at least several hundred million years. In Figure 15 we noted only weak dependencies of both $\eta_{\rm M}$ and $v_{\rm out}$ on instantaneous black hole luminosity, likely due in part to this issue.

At the same time, it is evident that the highest velocity BH-driven outflows in the TNG50 simulation do not arise from high luminosity quasars (e.g. with space densities $n\lesssim 10^{-7}$ Mpc*-3*), which are not present in the small volume. Instead, they are generated from low-accretion rate BHs for which the model posits a non-relativistic wind mechanism produced by a radiatively inefficient flow (Blandford & Begelman, 1999). In TNG black holes which are slowly accreting and have low bolometric luminosities can drive some of the most powerful outflows.

Observational evidence has emerged in this direction, probing low accretion rate black hole populations possibly underrepresented in previous surveys (Cheung et al., 2016; Penny et al., 2018). Observational constraints on the outflow properties as well as the mass distribution of low luminosity AGN, particularly as a function of location in the SFR$-M_{\star}$ plane and so in relation to the quenching state of the galaxy, will provide an important assessment of the TNG black hole model in the future.

4.2 Towards quantitative comparisons with observations

Outflows are usually observed through the Doppler shift of specific gas tracers, either through blueshifted features in absorption (e.g. Steidel et al., 2002), or through broad line components in emission (e.g. Heckman et al., 1990). In both cases several complexities arise which make the inference of intrinsic physical values difficult.

First, the outflow signature itself must be separated from the host. Studying kinematics using absorption spectra when the galaxy stellar continuum is the background source is complicated by existence of strong stellar absorption features in the same lines; Chen et al. (2010) find for instance that $\sim$ 80% of Na D in stacked SDSS spectra arises from stellar atmospheres – careful subtraction is required and the signal of interest represents a potentially small residual. Handling emission which partially fills in absorption signatures from the ISM or resonantly scattered emission similarly requires careful modeling to back out the intrinsic absorption (e.g. for FeII and MgII; Zhu et al., 2015). Outflows derived from emission lines also face difficulties, particularly in the decomposition of narrow and broad-line components, inclination corrections, the treatment of the complex underlying ISM kinematics, and the diverse excitation mechanisms for nebular lines (Newman et al., 2012). Significant assumptions are still required to convert observed line properties into e.g. mass outflow rates or mass loading factors; in particular, on the geometry, metal content, and ionization state of the outflowing gas. The uncertainties therein are significant (Murray et al., 2007; Chisholm et al., 2016), and ignoring or improperly accounting for galaxy-to-galaxy variation in assumed metallicity, ionization corrections, and outflow radii can introduce a factor of ten uncertainty in $\dot{M}_{\rm out}$, possibly hiding (or producing misleading) trends (Chisholm et al., 2017). For a mass conserving outflow with density varying only as a function of distance from the galaxy, the expression for mass outflow rate is

[TABLE]

where $n_{\rm H}(r)$ is the hydrogen number density, $v_{\rm out}(r)$ is the dependence of the outflow velocity on distance, $\Omega=4\pi$ for a unity covering factor, and $\mu(r)/m_{\rm p}$ the mean molecular weight. In terms of a column density $N_{\rm H}$ observed along the line of sight, the mass flow rate across a shell at radius $r$ is

[TABLE]

Alternatively, given a total, uniform mass $M_{\rm gas}$ at a constant outflow speed $v_{\rm out}$ within an outer radius of $r_{\rm out}$, dimensional analysis yields an even simpler expression for total outflow rate of

[TABLE]

In all cases, observations must assume values for $\Omega$ and $r_{\rm out}$, which are largely or entirely unconstrained (Harrison et al., 2018). The functional form of $v_{\rm out}(r)$ is also typically unknown and often assumed to be constant, taking some derivable characteristic velocity. Finally, whatever phase of the gas is observed represents only a fraction of the total outflow (Cicone et al., 2018), and for metal ion lines this can be a truly minuscule component. Inversion to the total hydrogen mass is therefore crucial, and yet also error prone due to the large correction factors. These complexities, which together accumulate into large uncertainties, motivates a comparison in the other direction. That is, forward modeling of the simulation outcome into the space of direct observables.

From the simulation side, there are several difficulties in modeling the actual gas phases. First, many important observational tracers are very cold, dense phases at $T\lesssim 10^{4}$ K, including neutral or molecular gas. This is an unresolved regime in many current cosmological models, due primarily to resolution limitations. Second, a standard ionization post-processing treatment, as commonly applied in studies in the low-density CGM for instance, is insufficient. This is because local radiation sources are important so close to the galaxy, from stars as well as possibly a central AGN, which are present in addition to the uniform, meta-galactic background. A treatment here would require an expensive, essentially radiative transfer calculation to derive the total incident radiation on each gas parcel, as well as assumptions on free parameters such as escape fractions from young stars. In addition, the ionization calculation is complicated by non-equilibrium effects due to rapid temporal evolution, as arise in shocks. This is a difficult problem even for specialized photoionization codes (e.g. MAPPINGS; Sutherland et al., 2018). It is also unclear how to efficiently capture the needed short time-scale information from a cosmological simulation.

Overcoming these challenges would imply a level of modeling detail which would also enable direct predictions of synthetic nebular emission lines from the ISM itself (Gutkin et al., 2016; Byler et al., 2017), as well as extended and extraplanar ionized gas reservoirs (e.g. Jones et al., 2017), in the context of cosmological galaxy simulations (Hirschmann et al., 2017). This is the inevitable direction of future modeling efforts and will be an important step to bridge the gap between hydrodynamical simulations and many of the most accessible and important observational probes of galaxies. For outflows, we will then be able to directly calculate phase-specific predictions – for example, $\eta_{\rm M,MgII}$ or $v_{\rm out,H\alpha}$. This will facilitate direct quantitative comparisons with observables, and allow us to assess the ‘tip of the iceberg’ effect and the importance of gas phases missing in current observations of galactic outflows.

5 Conclusions

In this work and together with the companion paper \textcolorbluePillepich et al. (2019) we have presented the new TNG50 simulation, the third and final volume of the IllustrisTNG project. TNG50 has been designed to overcome the resolution limitations inherent in cosmological simulations by sampling a large, statistically representative volume at unprecedented numerical resolution. Leveraging this new resource, we describe a first exploration of the properties of galactic outflows – driven by both supernovae and black hole feedback at $z>1$ – in the cosmological setting and with respect to the galaxies from which they arise. By resolving the internal structure of individual galaxies, TNG50 enables us to study the connection between small-scale (i.e. few hundred pc) feedback and large-scale (i.e. few hundred kpc) outflows. The diversity of galaxies and outflow properties realized in the fully representative galaxy sample of TNG50 highlights the importance of the large cosmological volume. We summarize our principal results:

• •

The new TNG50 simulation occupies a unique combination of large volume and high resolution, with a 50 Mpc box sampled by $2160^{3}$ gas cells. This provides a baryon mass resolution of $8\times 10^{4}$ M*⊙* and an average spatial resolution of star-forming ISM gas of $\sim 100-200$ parsecs. This resolution approaches or exceeds that of modern ‘zoom’ simulations of individual galaxies, while the volume contains $\sim$ 20,000 resolved galaxies with $M_{\star}\ga 10^{7}$ M*⊙*. It offers an unparalleled view into the small-scale structure of galaxies and the co-evolution of dark matter, gas, stars, and supermassive black holes from dwarfs to massive ellipticals. (§2.1 and §3.1)

• •

By measuring mass outflow rates around galaxies we quantify the mass loading factor $\eta_{\rm M}=\dot{M}_{\rm out}/\dot{M}_{\star}$ as a function of stellar mass, redshift, distance, and velocity. Whereas the TNG model input value for $\eta_{\rm M}$ at the injection scale of stellar feedback monotonically decreases with $M_{\star}$, this is not the case for the emergent outflows at 10 kpc. Instead, the mass loading inverts and rises rapidly for $M_{\star}\gtrsim 10^{10.5}$ M*⊙*. As a result, $\eta_{\rm M}(M_{\star})$ has a broken ‘v’ shape: the low-mass behavior is regulated by stellar feedback, while the high-mass behavior is set by the energetics and coupling efficiency of BH feedback, as it begins to drive strong outflows away from quenching, low-SFR systems. Strong outflows reach larger distances towards low redshift, and $\eta_{\rm M}$ decreases monotonically with increasing galactocentric distance for $r>10$ kpc. (§3.2)

• •

We extract the velocity of outflows as a function of $M_{\star}$, redshift, and distance, characterizing heterogeneous $v_{\rm out}$ distributions around individual galaxies through mass flux weighted velocity percentiles such as $v_{\rm out,95}$. In this case, median outflow velocities increase from $\sim$ 200 km/s at $M_{\star}=10^{7.5}$ M*⊙* to $\sim$ 1000 km/s at $M_{\star}=10^{11}$ M*⊙*. At fixed stellar mass, outflows are faster at higher redshift, by roughly a factor of two from $z=1$ to $z=6$. Despite a ‘minimum’ launch velocity of 350 km/s imposed by the TNG model, low-mass galaxies host much slower outflows due to resolved halo drag. At the Milky Way mass scale, despite high wind velocities of $\gtrsim$ 800 km/s prescribed by the model, the bulk of outflowing material escapes the galaxy at $\sim$ 150 km/s. The provided scalings $v_{\rm out}\propto M_{\star}^{\alpha}$ depend on redshift and velocity percentile, and $v_{\rm out}$ always declines with galactocentric distance. At the high-mass end, BH feedback produces high-velocity outflows at speeds exceeding $\sim 3000$ km/s to 20 kpc and beyond into the inner halo. (§3.3)

• •

Outflows from TNG galaxies are multiphase, and different phases have different kinematics. In the temperature distribution of outflowing gas, the dominant component is a ‘hot’ phase whose peak temperature increases with the virial temperature of the host halo gas. At $M_{\star}\gtrsim 10^{10.5}$ M*⊙*, however, a second cool outflow component emerges, resulting in a bimodal outflow temperature distribution. This colder gas is also denser, and is a direct consequence of ejective BH-driven feedback launching formerly nuclear ISM material out into the halo. The highest outflow velocities are achieved by the hottest outflow phases. (§3.4)

• •

Despite the directional isotropy of all energy inputs from both SN and BH feedback mechanisms, outflows are found preferentially aligned with the minor axes of galaxies. These bipolar outflows are therefore an emergent feature of our simulations, resulting from hydrodynamical collimation. Strong angular anisotropy is evident already by $\sim$ 10 kpc, and persists out as far as strong outflows are present ($\gtrsim$ 50 kpc at least). It exists at both low and high mass scales, where SN and BH feedback, respectively, dominate the production of galactic-scale outflows. Collimation increases with time, being ubiquitous by $z=1$ but difficult to convincingly demonstrate even at $z=2$. We speculate that this redshift trend is associated with the epoch of disk settling, where the rise of orderly rotating gas structures helps to sculpt outflow propagation. (§3.5)

• •

We demonstrate a correlation of $v_{\rm out}$ with $\Delta$SFMS, such that star-forming galaxies above the main sequence drive faster winds than otherwise. This relationship inverts for $M_{\star}\gtrsim 10^{10.5}$ M*⊙*, where quenching galaxies far below the main sequence launch the fastest outflows as a result of BH feedback. (§3.6)

• •

We measure the trends of outflow properties ($v_{\rm out}$ and $\eta_{\rm M}$) with galactic star formation rate, black hole bolometric luminosity, and the spatial distribution and concentration of star formation activity. Where available we compare against a large, heterogeneous mix of observations spanning all possible galaxy types and selections, redshifts, gas phases, and outflow tracers. While data occupies similar parameter ranges as our intrinsic predictions from TNG, a robust comparison with observations of outflow properties requires additional modeling. In TNG, low luminosity, slowly accreting black holes can drive some of the most powerful outflows. (§3.6)

In the future, similar measurements of outflow properties can be made in other hydrodynamical cosmological simulations, facilitating meaningful comparisons of the feedback physics implemented in different models. This will be particularly powerful because the gas dynamics of resolved outflows and halo-scale flows in general are followed faithfully without subgrid prescriptions. It remains to be seen if other models can produce similarly realistic galaxy populations with radically different outflow properties, or vice versa.

In this work we have shown that inputs and parameters of the feedback model at the injection scale differ from the emergent properties of galactic-scale outflows. Despite the relative simplicity of the physical assumptions invoked on the smallest scales, the resulting outflow properties are diverse and complex. That is, model parameterizations do not translate directly into observable outflow signatures. The rich phenomenology of outflows encodes the effective functioning of feedback physics and provides a unique way to study the impact of the baryon cycle on galaxy formation.

Acknowledgements

DN would like to thank Kate Rubin, Gwen Rudie, and Alice Shapley for insightful discussions and suggestions, as well as the anonymous referee for a constructive report. SG, through the Flatiron Institute, is supported by the Simons Foundation. The primary TNG simulations were realized with compute time granted by the Gauss Centre for Supercomputing (GCS): TNG50 under GCS Large-Scale Project GCS-DWAR (2016; PIs Nelson/Pillepich), and TNG100 and TNG300 under GCS-ILLU (2014; PI Springel) on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS). GCS is the alliance of the three national supercomputing centres HLRS (Universität Stuttgart), JSC (Forschungszentrum Jülich), and LRZ (Bayerische Akademie der Wissenschaften), funded by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Württemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF). Additional simulations were carried out on the Draco and Cobra supercomputers at the Max Planck Computing and Data Facility (MPCDF).

Appendix A Halo-normalized velocities

In this Appendix we revisit two primary results, namely the trends of mass loading and outflow velocity as a function of stellar mass. For $\eta_{\rm M}$ $(M_{\star})$ we change the $v_{\rm rad}$ threshold from a constant value in km/s to a threshold which scales with the halo mass. We likewise present the trends of $v_{\rm out}(M_{\star})$ with a halo-based normalization. In both cases, we take a halo virial velocity defined as $v_{\rm 200}$, using the spherical overdensity measurements for each halo. For reference, at $z=2$, $v_{\rm 200}$ increases from $\sim 25$ km/s at $M_{\star}=10^{8}\,\rm{M}_{\odot}$ to $\sim 100$ km/s at $M_{\star}=10^{11}\,\rm{M}_{\odot}$, reaching $\sim 150-200$ km/s only for the most massive $M_{\star}\gtrsim 10^{11.5}\,\rm{M}_{\odot}$ galaxies.

The top panel of Figure 17 shows $\eta_{\rm M}$ as a function of $M_{\star}$ for four different thresholds: $v_{\rm rad}>\{0,1,5,10\}v_{\rm 200}$. As with a constant threshold, increasingly restrictive choices lead to lower measured mass outflow rates and hence mass loading values. For reasonable choices of velocity threshold, the results of Figure 5 at the high-mass end are qualitatively and quantitatively similar. $\eta_{\rm M}$ still shows a minimum at $M_{\star}\sim 10^{10.5}\,\rm{M}_{\odot}$, rapidly increasing above this mass due to strong BH-driven winds.

The main difference arises at the low-mass end, where the use of a halo-normalized velocity threshold allows slow outflows driven by relatively shallow potential halos to still contribute. As a result, no threshold value causes a turnover and subsequent decline of $\eta_{\rm M}$ towards the smallest stellar masses, as was the case for the $v_{\rm rad}>150$ km/s and $v_{\rm rad}>250$ km/s choices of Figure 5. Based on this zeroth order comparison with halo circular velocity, it is clear that low mass galaxies continue to drive outflows which are ‘fast’, relatively speaking.

In the bottom panel of Figure 17 we show the trend of outflow velocity as a function of stellar mass, again at $z=2$ and for a fixed distance of 20 kpc from the galaxy. In contrast to Figure 6, velocities are in each case normalized to $v_{\rm 200}$ of the parent halo. Three different percentiles are contrasted: $v_{\rm out,50}$, $v_{\rm out,75}$, and $v_{\rm out,95}$. As expected, measurements further into the tails of the distribution reveal gas flowing at larger velocities relative to the halo velocity. The trend with $M_{\star}$ in the first two cases is roughly constant, and the typical outflow velocities are between $\sim 2-5v_{\rm 200}$ and $\sim 4-7v_{\rm 200}$, respectively. The high velocity tail starts to reveal different behavior, increasing relative to $v_{\rm 200}$ at both the low mass and high mass ends. The former can be understood in terms of the minimum injection velocity of 350 km/s for the TNG wind model, and the latter in terms of the rise of fast BH-driven outflows. In both cases the fastest moving components of the outflow can exceed $10v_{\rm 200}$ at these distances.