A model for cosmological simulations of galaxy formation physics

Abstract

1 INTRODUCTION

In the standard model of cosmology, the dominant mass contribution to the Universe is in the form of cold dark matter (CDM). Together with a large dark energy (DE) component (Perlmutter et al. 1999; Riess et al. 1999), this is the defining feature of the concordance ΛCDM cosmology. Early models of structure formation have predicted the condensation of baryons via radiative cooling at the centres of a population of hierarchically growing dark matter (DM) haloes, forming galaxies (Rees & Ostriker 1977; Silk 1977; White & Rees 1978; Blumenthal et al. 1984). Galaxy formation theory aims to explain the observed galaxy population of our Universe through a self-consistent model which requires as input only the initial conditions left behind after the big bang, as inferred from the cosmic microwave background (CMB). Such a model needs to correctly predict the evolution of all constituents of the Universe: DM, DE and baryons. Despite the fact that it is still unknown what DM and DE are, numerical N-body simulations of these ‘dark components’ have reached a high degree of sophistication over the past decade, making possible accurate calculations far into the non-linear regime (e.g. Springel et al. 2005b, 2008; Boylan-Kolchin et al. 2010; Klypin, Trujillo-Gomez & Primack 2011; Wu et al. 2013).

However, to make direct connections with the large number of available and upcoming observations of the visible universe, simulations must also account for baryonic physics. Especially in the era of ‘precision cosmology’ based on large surveys like SDSS, LSST, etc., it is crucial to have predictive, accurate and reliable galaxy formation models to test our understanding of structure formation by confronting these models with a plethora of data. Including baryons in cosmological models is also a necessity for predicting the back-reaction of baryons on the DM distribution, an effect that is particularly important for observational searches of DM and for concluding whether the ΛCDM model is viable on small scales (Vogelsberger et al. 2009; Ling et al. 2010).

There are currently two main approaches for computing cosmological models of galaxy formation and evolution: so-called semi-analytic models and self-consistent hydrodynamical simulations. Semi-analytical models (e.g. White & Frenk 1991; Kauffmann, White & Guiderdoni 1993; Baugh 2006; Croton et al. 2006; Guo et al. 2011; Benson 2012) account for baryonic physics by adding galaxy formation as a post-processing step on top of the output of N-body simulations or on to Monte Carlo realizations of DM merger history trees. Specifically, this approach uses simple analytical prescription to describe the baryonic physics that shape galaxy populations. Typically, the free parameters of the models are tuned to reproduce certain key observations, like the local stellar mass function, and are then used to predict other characteristics of the observed galaxy population. By contrasting those predictions to observations, the model can then be refined. While this approach allows a fast parameter space exploration of the relevant physics at play, its predictive power for certain observables is limited. For example, semi-analytic models do not permit a detailed prediction for the gas properties of the intergalactic medium (IGM), let alone the circumgalactic medium in and around galaxy haloes. Moreover, since the gas dynamics in the semi-analytic approach is approximated only in a very crude way, it is difficult to correctly identify the primary physics driving the outcome.

The most general way to overcome these limitations is to couple the gas physics to the ‘dark components’ and perform fully self-consistent hydrodynamical simulations in cosmological volumes. This approach is significantly more expensive from a computational point of view than semi-analytic models. In addition, modelling the required physics in a numerically meaningful way is also challenging. The difficulties arise from two different considerations: (i) the underlying numerical technique must solve the hydrodynamical equations efficiently and accurately, and offer a sufficiently large dynamic range to describe both cosmological and galactic scales, and (ii) the relevant baryonic processes must be implemented in a physically and numerically meaningful way. Especially the latter point should not be viewed as an afterthought – owing to inevitable technical limitations, brute-force implementations of the physics are often ill-posed numerically, non-convergent and may produce spurious results. These limitations can be overcome only by employing subgrid prescriptions that provide numerical closure, albeit at the cost of some modelling uncertainty. Such subresolution models aim to link the unresolved small-scale physics within galaxies to the scales that can actually be resolved in large-scale cosmological simulations. We stress that any numerical model always requires closure at some scale, since cosmological simulations cannot resolve all scales. Schemes without closure naturally lack convergence. This can be seen for example in the recent comparison paper by Scannapieco et al. (2012) and in more detail in Marinacci, Pakmor & Springel (2013), where different numerical implementations were compared. Fig. 21 (bottom left) in Marinacci et al. (2013) demonstrates that the convergence properties of a numerically well-posed scheme (as presented below) are significantly better than those of other methods. We note that this does not mean that smaller scales can and should not be resolved in future simulations. But we stress that this has to be done in a numerically meaningful way such that convergence can be achieved.

Many previous studies have argued that accurate schemes for solving the inviscid Euler equations are crucial for reliably modelling galaxy formation (Frenk et al. 1999; Agertz et al. 2007; Mitchell et al. 2009; Sijacki et al. 2012). Recently, we have corroborated this point explicitly through a detailed comparison campaign between simulations carried out with the moving-mesh code arepo (Springel 2010) and those with smoothed particle hydrodynamics (SPH) with the gadget (last described in Springel, Di Matteo & Hernquist 2005a) code. The differences we identified are sufficiently large to qualitatively change the way galaxies accrete their gas (Nelson et al. 2013), and to quantitatively modify the strength of required feedback processes. We stress that these results were obtained with ‘standard SPH’ as implemented in the gadget code. We hence consider it essential to eliminate numerical artefacts and errors in galaxy formation simulations as far as possible. Without doing so, the physical models used to describe effects associated with star formation (SF), black hole (BH) growth and feedback will be corrupted by numerical limitations. Furthermore, galaxy formation simulations should ideally be performed at the highest possible resolution in order to faithfully model also small galaxies and to limit the amount of relevant physics modelled as ‘subgrid’. This clearly requires a very good scaling behaviour of the employed numerical scheme to allow the use of state-of-the-art supercomputers.

Over the last several years, a large body of results based on hydrodynamical simulations of galaxy formation have been obtained (e.g. Ocvirk, Pichon & Teyssier 2008; Crain et al. 2009; Schaye et al. 2010; Davé, Oppenheimer & Finlator 2011a; McCarthy et al. 2012; Kannan et al. 2013; Puchwein & Springel 2013). A detailed investigation of the impact of various physical mechanisms on cosmic gas can be found in Schaye et al. (2010) as well as in other recent studies (e.g. Haas et al. 2013a,b). A consensus from these works is that feedback processes are crucial for any successful model of galaxy formation. For example, Davé et al. (2011a) focused on the impact of stellar winds on the faint end of the stellar mass function. Puchwein & Springel (2013) tried to match the stellar mass function and the stellar mass–halo mass relation through stellar and active galactic nucleus (AGN) feedback, but their simulations did not include effects like metal-line cooling and gas recycling. Other studies focused specifically on low-mass haloes and showed that specific forms of stellar feedback may be successful in producing the correct amount of stellar mass in these systems (e.g., McCarthy et al. 2012; Kannan et al. 2013).

We note that these cosmological simulations of the formation of representative galaxy populations are markedly different from studies where the computational power is focused on a single galaxy (e.g. Governato et al. 2010; Agertz, Teyssier & Moore 2011; Guedes et al. 2011; Genel et al. 2012; Scannapieco et al. 2012; Aumer et al. 2013). In the latter case, a very high spatial and mass resolution can be achieved, allowing one to resolve the physics to smaller scales. However, one does not obtain a statistical sample of objects to compare to observations. Moreover, the impact of numerical artefacts inherent to some methods such as SPH varies significantly from one halo to another (e.g. Sales et al. 2012; Vogelsberger et al. 2012). This means that physics models calibrated against individual objects will be corrupted by numerical errors in unpredictable ways, rendering their applicability to galaxy populations problematic. These considerations limit the ability of such ‘zoom-in’ simulations by themselves to constrain galaxy formation models. Nevertheless, when used in conjunction with large-scale simulations, the zoom-in procedure will be essential for constructing more reliable subresolution models for characterizing entire galaxy populations.

Modern cosmological hydrodynamical simulations typically include a subresolution model for SF (e.g. Ascasibar et al. 2002; Springel & Hernquist 2003a; Dubois & Teyssier 2008; Schaye & Dalla Vecchia 2008; Few et al. 2012), radiative cooling (e.g. Katz, Weinberg & Hernquist 1996; Wiersma, Schaye & Smith 2009a) and chemical enrichment (e.g. Steinmetz & Mueller 1994; Mosconi et al. 2001; Lia, Portinari & Carraro 2002; Springel & Hernquist 2003a; Kobayashi 2004; Scannapieco et al. 2005; Tornatore et al. 2007; Oppenheimer & Davé 2008; Wiersma et al. 2009b; Few et al. 2012). Furthermore, some small-scale simulations also include magnetic fields (e.g. Teyssier, Fromang & Dormy 2006; Dolag & Stasyszyn 2009; Collins et al. 2010; Pakmor, Bauer & Springel 2011; Pakmor & Springel 2013), radiative transfer (e.g. Abel & Wandelt 2002; Cantalupo & Porciani 2011; Petkova & Springel 2011), cosmic ray physics (Jubelgas et al. 2008) and thermal conduction (Dolag et al. 2004; Jubelgas, Springel & Dolag 2004). Most importantly, hydrodynamic simulations of galaxy formation usually include some form of stellar feedback (e.g. Dekel & Silk 1986; Navarro & White 1993; Mihos & Hernquist 1994; Gerritsen & Icke 1997; Thacker & Couchman 2000; Kay et al. 2002; Kawata & Gibson 2003; Sommer-Larsen, Götz & Portinari 2003; Springel & Hernquist 2003a; Brook et al. 2004; Oppenheimer & Davé 2006; Stinson et al. 2006, 2013; Dalla Vecchia & Schaye 2008, 2012; Dubois & Teyssier 2008; Okamoto et al. 2010; Piontek & Steinmetz 2011), and more recently also AGN feedback (e.g. Di Matteo, Springel & Hernquist 2005; Kawata & Gibson 2005; Springel et al. 2005a; Sijacki & Springel 2006; Thacker, Scannapieco & Couchman 2006; Sijacki et al. 2007; Di Matteo et al. 2008; Okamoto, Nemmen & Bower 2008; Booth & Schaye 2009; Kurosawa & Proga 2009; Debuhr, Quataert & Ma 2011; Teyssier et al. 2011; Dubois et al. 2012).

In a recent series of papers (Kereš et al. 2012; Sijacki et al. 2012; Torrey et al. 2012; Vogelsberger et al. 2012; Bird et al. 2013; Nelson et al. 2013), we have demonstrated that a moving-mesh approach as implemented in the arepo code (Springel 2010) provides a promising numerical method for performing galaxy formation calculations. Indeed, such a scheme combines most of the advantages of previous Eulerian and pseudo-Lagrangian techniques used in the field while avoiding many of their weaknesses. We therefore concluded in our previous work that the quasi-Lagrangian finite-volume scheme of arepo is well suited for reliably solving the hydrodynamical equations posed by the galaxy formation problem. However, the simulations presented in Vogelsberger et al. (2012) accounted for only a rather limited set of physical processes in order to simplify comparisons with earlier numerical approaches. Most importantly, they did not include any form of explicit stellar feedback and also did not account for feedback from AGN. But it has become clear over the last decade that both forms of feedback are crucial for shaping, for example, the stellar mass function. The stellar mass function may be viewed as a convolution of the underlying DM halo mass function with the efficiency with which stars form in these haloes (Springel & Hernquist 2003b). Stellar feedback regulates SF in low-mass haloes, since their potential wells are sufficiently shallow such that energy supplied through supernovae (SN) can remove gas efficiently from the galaxy, yielding galactic winds and outflows. For more massive systems, this channel of feedback becomes inefficient. However, as was realized already many years ago, AGN feedback from BHs can supply sufficient energy to quench SF in these more massive systems. A combination of both feedback channels is the leading theoretical conjecture to explain the observed stellar mass function in the local Universe.

Besides these feedback mechanisms, other physical processes are known to be important but were not included in our first moving-mesh simulations of galaxy formation presented in Vogelsberger et al. (2012). In particular, these simulations considered cooling processes only through a primordial mixture of hydrogen and helium. But evolving and dying stars release significant amounts of heavier elements, which leads to an important source of additional atomic line cooling once these metals are mixed with ambient gas. Stellar evolution processes also lead to gas recycling such that baryonic matter is not completely locked up in stars in a permanent fashion; rather it is returned at some rate to the gas phase, even by old stellar populations. Both metal-line cooling and gas recycling elevate the supply of cold gas to galaxies and hence tend to increase the star formation rates (SFRs) in these systems.

In this paper, we present a new implementation of galaxy formation physics in the moving code arepo, suitable for large-scale hydrodynamical simulations of structure formation. Our goal is to account for all primary physical processes that are known to be crucial for shaping the galaxy population. While we aim to be reasonably comprehensive in this endeavour, we note that it is clear that we cannot be fully complete in this first instalment of our model. For example, we neglect some ‘non-standard’ physics that has been proposed to influence galaxy formation, such as cosmic rays and magnetic fields. Besides presenting our implementation we also demonstrate how different feedback effects influence the properties of galaxies and how the simulated galaxy population compares to current observations. We also identify the physical parameters of a fiducial reference model that leads to reasonable agreement with many key observables.

This paper is organized as follows. In Section 2, we present our galaxy formation model. Analysis and post-processing techniques that we have implemented in arepo are described in Section 3. The first cosmological simulations using the newly added physics, and a comparison to different observational data sets, are discussed in Section 4. Here we also explore different feedback settings and introduce a fiducial reference model which provides a good fit to a set of important key observations. Finally, we give our conclusions and summary in Section 5.

2 GALAXY FORMATION MODEL

In the following subsections, we describe the physics implementation of our galaxy formation model. We first give a very brief description of our numerical scheme for solving the hydrodynamical equations. We then present our subresolution model for SF followed by a discussion of our stellar evolution model. Next we discuss our approach for handling gas cooling. Finally, the second half of this section covers our methods for including BHs and feedback processes (stellar and AGN).

2.1 Hydrodynamical method

Most previous hydrodynamical simulations of galaxy formation have employed the SPH technique (Gingold & Monaghan 1977; Lucy 1977; Monaghan 1992, 2005), where gas is discretized into a set of particles for which appropriate equations of motion can be derived. The SPH method is well suited for cosmological applications owing to its pseudo-Lagrangian character that automatically adjusts its resolution when structures collapse gravitationally. Furthermore, SPH manifestly conserves energy, momentum, mass, entropy and angular momentum, a feature unmatched by competing techniques.¹ Other methods for solving the equations of hydrodynamics in a cosmological context have been employed as well, most of which are descendants of classic Eulerian methods (Stone & Norman 1992) implemented as adaptive mesh refinement (AMR) codes (Berger & Colella 1989; Teyssier 2002; O'Shea et al. 2004). Several studies have however pointed out that these widely used schemes can lead to significant differences in the results for galaxy formation (e.g. Frenk et al. 1999; Agertz et al. 2007; Mitchell et al. 2009).

In this paper, we use a new numerical approach where a moving mesh is used to solve the hydrodynamical equations (the arepo code; Springel 2010), based on a quasi-Lagrangian finite-volume method. The primary motivation is similar to earlier implementations of moving meshes by Gnedin (1995) and Pen (1998), but our method avoids the mesh-twisting problems that troubled these first attempts. To this end, arepo employs an unstructured Voronoi tessellation of the computational domain. The mesh-generating points of this tessellation are allowed to move freely, offering significant flexibility for representing the geometry of the flow. This mesh is then used to solve the equations of ideal hydrodynamics with a finite-volume approach using a second-order unsplit Godunov scheme with an exact Riemann solver. This technique avoids several of the weaknesses of the SPH and AMR methods while retaining their most important advantages (see the detailed discussion in Springel 2010; Kereš et al. 2012; Sijacki et al. 2012; Torrey et al. 2012; Vogelsberger et al. 2012). We note that a number of SPH extensions were proposed, which solve some of the defects of standard SPH (e.g. Price 2008; Wadsley, Veeravalli & Couchman 2008; Cullen & Dehnen 2010; Abel 2011; Read & Hayfield 2012; Saitoh & Makino 2013).

2.2 Star formation

Owing to the limited resolution in large-scale structure simulations, we must treat the star-forming interstellar medium (ISM) in a subresolution fashion. Here we follow Springel & Hernquist (2003a) and model the star-forming dense ISM gas using an effective equation of state (eEOS), where stars form stochastically above a gas density of ρ_sfr with a star formation time-scale of t_sfr. Such an approach for simulations with an unresolved ISM structure is very common (e.g. Ascasibar et al. 2002; Springel & Hernquist 2003a; Springel, Di Matteo & Hernquist 2005c; Robertson et al. 2006; Dubois & Teyssier 2008; Hopkins et al. 2008, 2009a,b; Schaye & Dalla Vecchia 2008; Few et al. 2012). The idea behind the use of an eEOS is that the ISM is believed to be governed by small-scale effects like supersonic turbulence, thermal instability, thermal conduction, and molecular cloud formation and evaporation processes. The combination of these processes is believed to quickly establish a self-regulated equilibrium state of the ISM. In that regime, the average ISM temperature may be approximated as a function of density only.

We introduce a few modifications of the original Springel & Hernquist (2003a) implementation. First, to seamlessly connect to our stellar evolution model (see below), we use a Chabrier (2003) initial mass function (IMF) instead of the Salpeter IMF. Note that this also implies a larger Type II supernova (SNII) energy input in the ISM per formed solar mass of stars. Second, we take stellar mass-loss processes into account when deriving the SFR based on the cold gas fraction in the ISM. Third, we prevent cells from being star forming (both in terms of their SFR and the eEOS) if their temperature is above the temperature corresponding to the eEOS at their density. This prevents spurious SF in diffuse hot gaseous haloes. Fourth, for determining the temperature of star-forming gas, we interpolate between the full Springel & Hernquist (2003a) eEOS and an isothermal EOS at 10⁴ K with an interpolation parameter q_EOS = 0.3, which avoids an overpressurization of the ISM (Springel et al. 2005a; Hopkins & Quataert 2010). We did not incorporate metal-dependent density thresholds nor did we take into account metal cooling within the two-phase model. Although adopting this would be straightforward, these effects do not lead to any significant changes in the functionality of our model; hence, we ignore them here in order to keep the model as simple as possible. The adopted parameter choices for our SF model are ρ_sfr = 0.13 cm⁻³ and t_sfr = 2.2 Gyr, in the notation of Springel & Hernquist (2003a). Our ISM subresolution model alone, by construction, does not drive galactic winds. We will discuss this below in more detail and present a model for the explicit generation of galactic outflows.

We note that once the numerical resolution is sufficient to resolve all the relevant ISM phases and processes, an eEOS is in principle not necessary anymore. In this case, a direct modelling of the relevant physical ISM processes, like momentum/energy injection and radiation pressure associated with SF, is expected to lead to a turbulent and self-regulated quasi-equilibrium structure of the ISM. Impressive progress in this direction has recently been achieved (e.g. Hopkins, Quataert & Murray 2012a,b). However, we point out that the resolution requirements for such an approach are significantly beyond what is currently feasible in large-scale cosmological simulations; hence, this is applicable only in ‘zoom-in’ simulations of individual galaxies of low to moderate mass or in idealized simulations of isolated galaxies (Hopkins et al. 2012c, 2013a; Hopkins, Keres & Murray 2013b). In a forthcoming study (Marinacci et al., in preparation), we will present a much more detailed ISM prescription on top of the other physics implementations discussed below.

2.3 Stellar evolution and chemical enrichment

Elemental abundance patterns are an important diagnostic tool for galaxy formation. Cosmological simulations track the evolution of baryons as they fall into galaxies, form stars and are recycled back into the ISM. Including estimates for the recycling of material from ageing stellar populations is important because it plays a central role in determining the evolution of a galaxy's gas content and because the metal-rich ejecta from evolved stars determine its heavy element content.

Handling the mass and metal return from massive stars can be well approximated by assuming instantaneous massive star evolution, which is implemented by lowering the effective SFR while the metal content of nearby gas is increased to account for metal production (e.g. Steinmetz & Mueller 1994; Springel & Hernquist 2003a). This instantaneous recycling approach is justified for massive stars because they have relatively short lifetimes ( ≲ 10⁷ yr). However, an instantaneous recycling model cannot be accurately applied to intermediate-mass stars, as the time-scale for stellar evolution, and therefore the time-scale for mass return, is comparable to or larger than galactic dynamical times. Most notably, asymptotic giant branch (AGB) stars return a substantial fraction of their mass to the ISM over expected lifetimes that can exceed 10⁹ yr. Similarly, Type Ia supernovae (SNIa) explode with substantial delay times and, although they do not return a very large amount of mass to the ISM, they produce substantial quantities of iron. This gives rise to a situation where both the mass return rate and the elemental composition of the returned material will change as a function of time, depending on which mechanisms are operating.

A more general enrichment approach can be developed by tracking the mass and metal return of a stellar population to the ISM as a function of time (e.g. Mosconi et al. 2001; Lia et al. 2002; Kobayashi 2004; Scannapieco et al. 2005; Tornatore et al. 2007; Oppenheimer & Davé 2008; Wiersma et al. 2009b; Few et al. 2012). In our model, we integrate the time evolution of stellar particles to determine their mass-loss and chemical enrichment as a continuous function of time. We use information from stellar evolution calculations to determine the expected main-sequence lifetime, mass return fraction, and heavy element production for a wide range of initial stellar masses and metallicities. As a result, our model calculates not only the appropriate time-delayed mass-loss rate, but also the chemical composition of the mass-loss to the ISM. In the following, we briefly outline our scheme, which is similar to the implementation presented in Wiersma et al. (2009b).

2.3.1 Stellar evolution

We adopt the lifetime function from Portinari, Chiosi & Bressan (1998), which gives the expected stellar lifetimes as a function of the initial stellar mass and metallicity. These lifetimes are calculated as the sum of the hydrogen- and helium-burning time-scales for stars as a function of mass and metallicity.

2.3.2 Mass and metal return

We use the elemental mass yields for AGB stars from Karakas (2010), which are calculated by dynamically evolving the thermally pulsating AGB stars and then inferring the nucleosynthetic yields using a full reaction network over a wide range in metallicity (0.0001 < Z < 0.02) and initial stellar masses (1 < M < 6 M_⊙). We adopt the elemental mass yields for core collapse SN from Portinari et al. (1998), which are calculated using the Padova stellar evolutionary tracks in conjunction with the explosive nucleosynthesis of Woosley & Weaver (1995). We also extract f_rec(m, Z) from these references.

2.3.3 SN rates

The majority of a stellar population's mass-loss and metal production comes from AGB stars or core collapse SN. However, SNIa are also important to consider because they produce a substantial amount of iron. Besides, the SNIa rate is crucial for cosmological studies.

The progenitor systems of SNIa are still uncertain. The two most commonly discussed progenitor channels are the single degenerate scenario, where an SNIa results from a white dwarf accreting material from a companion star, and the double degenerate scenario, where an SNIa results from the merging of two white dwarf stars. In the standard single degenerate scenario, nearly all Type Ia progenitors will have the same mass and will be comprised mostly of carbon and oxygen. Since the conditions leading to a Type Ia event are comparatively well posed in this picture, the chemical yields expected in each SNIa have been calculated (e.g. Thielemann et al. 2003).

What remains nevertheless highly uncertain even in this scenario is the rate we expect for Type Ia events after the birth of an SSP. In contrast, determining the expected number and rate of AGB stars and core collapse SN is relatively straightforward once we have made a choice for the IMF and lifetime functions. This is because we expect that all stars of a given mass will end their lives in the same way (e.g. we assume that all 20  M_⊙ stars end as core collapse SN). The same simple logic cannot be applied to SNIa because of the uncertainty that exists about their progenitor systems and the long time delay between the birth of a star and its eventual explosion in a Type Ia event. Our lack of precise knowledge about the true form of the IMF, stellar binary fractions and Type Ia progenitor systems prevents us from implementing a first-principles-based model for SNIa rates.

With the SNIa rate specified, determination of the mass return and metal production is easily done using the yield calculations of Thielemann et al. (2003) and Travaglio et al. (2004). Mass and metals produced in Type Ia events are then returned to the ISM in the same fashion as for core collapse SN and AGB winds. We have checked that the exact choice of the DTD (e.g. power law versus exponential) does not affect most of our results in any significant way. The reason for this is that the SNIa mainly affect the iron production, but they are not a significant source of stellar mass-loss, which is dominated by stars in the AGB phase and SNII events.

2.3.4 Gas enrichment

For each active stellar particle, we calculate the total mass, total metal mass and element return during its current time step following the scheme described above. These masses are then returned to the nearby Voronoi cells of the stellar particle, where we enforce momentum and energy conservation. We distribute the mass and metals over a set of nearest neighbouring cells using a tophat kernel enclosing a mass of approximately N_enrich × m_target, where m_target is the cell target mass according to our (de-)refinement scheme (see below). For convergence studies, we normally hold this ‘enrichment mass’ fixed if we change resolution, which means that we scale N_enrich such that the change of the target mass with resolution is compensated. However, we find that using a constant N_enrich produces almost identical results.

2.4 Cooling and heating

Energy loss via radiative cooling is a key physical process for galaxy formation. Early cosmological hydrodynamical simulations (e.g. Katz, Hernquist & Weinberg 1992; Katz et al. 1996) included cooling due to a primordial mixture of hydrogen and helium, where two-body processes (collisional excitation, collisional ionization, recombination, dielectric recombination and free–free emission) and inverse Compton cooling off the CMB (e.g. Ikeuchi & Ostriker 1986) were considered. Furthermore, energy input due to photoionization was injected as heat into the gas. Typically these simulations considered a spatially uniform, time-varying UVB radiation field. The presence of metals in the ISM and IGM due to enrichment processes changes the cooling rates by increasing the number of possible transitions. A primordial treatment is then not accurate anymore for such an enriched gas.

For high-temperature gas, Sutherland & Dopita (1993) published cooling rate tables for 14 heavy elements over a range of metallicities, using solar relative abundances. Their rates assume collisional ionization equilibrium and do not include effects due to background radiation fields. However, it is important to include background radiation when calculating the cooling rates, because it affects both the thermal and ionization state of the plasma (e.g. Efstathiou 1992; Gnedin & Hollon 2012) leading to a suppression of the cooling rate and increase of heating rates. Smith, Sigurdsson & Abel (2008) introduced a table-based method to include metal-line cooling, accounting for photoionization due to the UVB. Their technique involves the use of the photoionization code cloudy to construct a lookup table of metal cooling rates. This approach was also followed in other studies; for example, Wiersma et al. (2009a) even considered the cooling contribution from individual elements.

We adopt a spatially uniform time-dependent UVB and also implement metal-line cooling in arepo based on cloudy cooling tables. We assume ionization equilibrium such that we can calculate cooling rates a priori for a given portion of dust-free and optically thin gas, for a given background radiation field. We still carry out a self-consistent calculation of the primordial cooling following Katz et al. (1996), on top of which we add the cooling contribution of the metals. In principle, the metal-line cooling can be implemented on an element-by-element basis once the main coolants are tracked by the enrichment scheme. Wiersma et al. (2009a) showed that it is sufficient to follow 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) to recover the main metal-line cooling contribution. Although most of these elements are explicitly followed by our chemical enrichment model, we decided not to implement cooling on an element-by-element basis in our default setup because (i) individual element yields have large uncertainties, and element-wise cooling would directly couple those into the hydrodynamical evolution via the cooling contribution; (ii) departures of the relative abundances from solar lead only to small corrections compared to the overall photoionization modification (see Wiersma et al. 2009a for example); and (iii) the computational overhead for evaluating the cooling rate scales directly with the number of elements used in the cooling. We base our implementation of metal-line cooling therefore on the rates for a solar composition gas, and these rates are scaled linearly with the total metallicity Z.

For simplicity, we neglect the metal contribution to the mean molecular weight μ, i.e. we calculate temperatures based on the molecular weight corresponding to the ionization state of the primordial gas composition. This is well justified since in a fully ionized metal enriched plasma of metallicity Z, the mean molecular weight is approximately given by μ ≅ 4/(8X + 3Y + 2Z), such that even at supersolar composition the metal contribution leads only to minor changes in μ (see also Smith et al. 2008). We also neglect the metal contribution to the free electron density, which is also justified since those will only be relevant for Z ≫ Z_⊙ (e.g. Wiersma et al. 2009a). For solving the primordial network, we use the advected H and He mass fractions as an input. Together with the spatially uniform but temporally varying UVB, this gives us the primordial abundances and a self-consistent temperature based on μ. With this temperature, we calculate Λ_p, and look up Λ_m in the pre-calculated cloudy tables to arrive at the total cooling rate through equation (12).

For the construction of the lookup tables, we tabulated metal-line net cooling rates on a grid in log (T/K), log (n_H/cm⁻³) and redshift: 1 < log (T/K) < 9 at 200 equally spaced grid points, −8 < log (n_H/cm⁻³) < 2 at 50 equally spaced grid points and 0 < z < 10 at 50 grid points. During the simulation, we linearly interpolate the metal-line cooling rates in [z, log (n_H/cm⁻³), log (T/K)] space.

2.4.1 Self-shielding

So far we have described our treatment of cooling and heating in the presence of a spatially uniform UVB radiation under the assumption that the gas is optically thin to the radiation. However, this approximation breaks down at gas densities ρ ≳ 10⁻³ cm⁻³, depending on the details of the radiation field. Above such densities, the gas absorbs the radiation to such a degree that the radiation field is attenuated compared to the optically thin case. This modified radiation field in turn affects the heating rate and ionization state, and hence the cooling rate, differently than the unattenuated radiation field. As a result, the equilibrium gas temperature ends up being different, and hence the dynamical evolution can be affected.

2.5 Stellar feedback

Although our subgrid ISM model implicitly invokes thermal SNII feedback, this is not sufficient to avoid the well-known overcooling problem (e.g. Springel & Hernquist 2003a). Large-scale cosmological simulations therefore typically resort to some subgrid model for galactic winds and outflows, capable of efficiently removing baryonic material from the star-forming phase. Two approaches are typically followed in a regime of limited resolution: injection of SN energy in kinetic form (e.g. Navarro & White 1993; Mihos & Hernquist 1994; Kay et al. 2002; Springel & Hernquist 2003a; Oppenheimer & Davé 2006; Dalla Vecchia & Schaye 2008; Dubois & Teyssier 2008; Okamoto et al. 2010; Hopkins, Quataert & Murray 2012c) or suppression of radiative cooling after injection of thermal energy (e.g. Gerritsen & Icke 1997; Thacker & Couchman 2000; Kawata & Gibson 2003; Sommer-Larsen et al. 2003; Brook et al. 2004; Stinson et al. 2006; Piontek & Steinmetz 2011; Dalla Vecchia & Schaye 2012; Stinson et al. 2013).

The kinetic wind models are typically characterized by a mass loading factor η_w, which gives the ratio of the wind mass flux to SFR, and an initial wind velocity v_w. In the original kinetic wind scheme of Springel & Hernquist (2003a), star-forming SPH particles are converted into wind particles stochastically. In order to allow a precise specification of the wind mass flux, the implementation ensured that wind particles can leave the dense ISM gas without entraining additional star-forming particles. In order to technically realize this subgrid prescription, hydrodynamical interactions were temporarily disregarded for newly launched wind particles until they ‘recouple’ just outside of the star-forming phase, based on a density threshold criterion. The mass loading and wind velocity in Springel & Hernquist (2003a) were chosen to be constant. More recent studies propose variable wind models, however, where the wind velocity is not held fixed. For example, Oppenheimer & Davé (2006) argue for momentum-driven winds, where the wind speed scales as the galaxy velocity dispersion and the mass loading is inversely proportional to that (see also Murray, Quataert & Thompson 2005). Such an approach appears to be supported by observations. Martin (2005) and Weiner et al. (2009) found that the final wind velocity scales approximately linearly with the circular velocity. Okamoto et al. (2010) use the local velocity dispersion of DM to prescribe the wind velocity, although they consider energy-driven winds, which are generated locally around the stellar particles similar to the approach of Dalla Vecchia & Schaye (2008), which uses constant local winds. Puchwein & Springel (2013) reported a good match to the observed stellar mass function with an energy-driven variable wind model, where the wind velocity is derived from the escape velocity of the host halo.

We have implemented two alternative wind schemes in arepo. Our basic approaches for modelling winds and outflows are similar to those presented in previous work (Springel & Hernquist 2003a; Oppenheimer & Davé 2006, 2008; Dalla Vecchia & Schaye 2008; Okamoto et al. 2010; Puchwein & Springel 2013), but our technical implementation is slightly different to account for our mesh-based hydro scheme. Also, we note that the overcooling problem tends to be more severe in arepo compared to previous SPH studies (see Kereš et al. 2012; Vogelsberger et al. 2012 for details). The origin of this lies in part in the treatment of subsonic turbulence and in spurious viscous heating in SPH, which leads to an artificial heating of gaseous haloes in SPH simulations (Bauer & Springel 2012; Vogelsberger et al. 2012). This effect is strongest for massive haloes, and we will show below that our simulations require stronger AGN feedback because of this (see below).

2.5.1 Non-local SNII feedback

In our first stellar feedback implementation, which we call ‘non-local’, we launch winds directly from the star-forming ISM gas, which is part of the two-phase medium. This ties the wind rate directly to the SFR of star-forming gas. For large-scale simulations, the time steps are typically not much shorter than the SNII delay time, i.e. it is a valid approximation to neglect locality for SNII-driven winds in that case.

We realize the generation of the wind at a technical level by turning gas cells (or a fraction of their mass) for a short period of time into massive wind particles which interact gravitationally but do not couple to the hydrodynamical calculation. For these non-local winds, we probabilistically select a star-forming ISM gas cell, which we then convert into a wind particle. This particle is then allowed to travel, decoupled from hydrodynamical forces, until it reaches a certain density threshold or a maximum travel time has elapsed. Once either of these criteria is fulfilled, we dissolve the particle and deposit its mass, momentum, thermal energy and all tracked metals into the gas cell in which it is currently located. The net effect of this procedure mimics the subgrid wind prescription of Springel & Hernquist (2003a) and facilitates a precise control of the wind parameters independent of numerical resolution. We note that Dalla Vecchia & Schaye (2008) demonstrated that decoupling wind particles can have an effect on the galaxy properties. However, a decoupled scheme gives better control over mass loadings and typically converges better, which is why we favour such an approach.

For choosing the wind direction, we implemented two commonly used approaches: isotropic and bipolar winds. For the isotropic case, each wind particle is kicked in a completely random direction. In the bipolar case, each wind particle is kicked, with a random sign, in the direction |${\boldsymbol v}\times \nabla \phi$| (Springel & Hernquist 2003a) in the rest frame of the corresponding FoF group. Note that both schemes do not add total momentum on average. In the following, we will only use the bipolar implementation.

2.5.2 Local SNII feedback

In the wind scheme described above, wind particles are launched from the ISM and are therefore not spawned directly from stellar particles. Although slightly inconsistent with our stellar evolution and enrichment model, where we use a delayed SNII return, this model works very well for cosmological applications because of the relevant time-scales, as argued above. Nevertheless, it is desirable to have the possibility of using a more local SNII feedback implementation for comparison. For example, Okamoto et al. (2010) use an energy-driven local stochastic scheme, where wind particles are launched in the neighbourhood of stellar particles. This approach is similar to the implementation of Dalla Vecchia & Schaye (2008) with the exception of the way in which ‘excess energy/mass loading’ is handled. Another variation of this approach was presented in Dalla Vecchia & Schaye (2012).

To explore the effect of non-locality, we implemented a wind creation scheme similar to that of Okamoto et al. (2010) for purely energy-driven winds, whose velocity scales with the local DM velocity dispersion. Specifically, we modified our enrichment routines to also distribute SNII energy among the gas cells that receive metals and mass. During the enrichment step, each cell keeps track of how much energy it receives from stellar particles such that once the enrichment is finished, some cells have received the SNII energy ΔE_SNII. Each of these cells is then assigned a probability |$p=\Delta E_{\rm SNII} / (0.5\,M v_{\rm w}^2)$|⁠, where M is the mass of the gas cell and v_w = κ_wσ is the wind velocity assigned to this particle. Based on this probability, it is then decided whether this cell should be converted to a wind particle. Once a cell is converted, it is treated in exactly the same way as in the non-local case described above. We note that we do not require gas cells to be above the density threshold for SF for this to happen.

2.5.3 Wind metal loading

Wind models with an energy or momentum scaling are the two most commonly applied wind parametrizations in the literature, and both rely on mass loading factors that decrease with galaxy mass. In fact, the mass loading factors associated with low-mass galaxies can often be well in excess of unity, naturally indicating that low-mass galaxies do not efficiently hold on to their gaseous content. Moreover, such high mass loading factors appear to be necessary for obtaining a reasonable shape for the galaxy stellar mass function (e.g. Davé et al. 2011a; Puchwein & Springel 2013). However, calculations of the metal content of low-mass galaxies based on their observed gas fractions and metallicities indicate that low-mass galaxies nevertheless retain a substantial fraction of their metals (Zahid et al. 2012). Thus, there is a certain tension between the need to efficiently expel a lot of material to reduce the SF efficiency in low-mass galaxies and the need to retain a sizeable fraction of the metals in the same low-mass galaxies.

The origin of this tension lies in our construction of the subgrid wind model. Traditionally, it is assumed that the wind particles have the same metallicity as the ambient ISM from where they are launched (e.g. Springel & Hernquist 2003a; Davé et al. 2011a). However, while the actual wind material will be partially composed of the material from the site where the wind is launched, a substantial fraction of the wind's mass will be entrained as the wind vents out of the galaxy and blows through low-density regions of the ISM. This suggests a picture where the metallicity of the wind material will almost certainly be different – and likely lower – than the metallicity of the ISM region from where the wind originated.

We found that the traditional γ_w = 1 leads to an underproduction of the gas-phase metallicity and an overenrichment of halo gas, especially in low-mass galaxies. Both of these problems directly result from an overly efficient ejection of metals from galaxies via winds. Not only does this lead to a poor match to the observationally constrained mass–metallicity relation, but the overenrichment of halo gas also leads to enhanced metal-line cooling. The other extreme, when using γ_w = 0, leads to overly enriched galaxy gas-phase metallicities compared to the mass–metallicity relations, suggesting that the correct behaviour is bracketed by these two extremes.

We note that wind metal loadings were recently discussed in Peeples & Shankar (2011), where it was argued that winds should actually be ‘superenriched’ such that their metallicity should be higher than the one of the surrounding ISM gas. We note however that this conclusion was reached based on the very low mass loading factors discussed in Peeples & Shankar (2011), which then automatically require highly enriched winds to match the observed mass–metallicity relation. Such low mass loadings are however incompatible with results from cosmological simulations; the latter require rather much higher mass loadings for low-mass systems to match the stellar mass function (see below).

2.6 BH growth and AGN feedback

As with SF, it is currently computationally impossible to resolve the details of accretion flows around central galactic BHs in large-scale simulations of galaxy formation. While the gravitational radius of supermassive BHs starts to be resolved in zoom-in simulations, the Schwarzschild radius of these BHs is still many orders of magnitude smaller than the spatial resolution limits of these calculations. Describing BH growth and related feedback therefore requires a subresolution model (e.g. Di Matteo et al. 2005, 2008; Kawata & Gibson 2005; Springel et al. 2005a; Thacker et al. 2006; Sijacki et al. 2007; Okamoto et al. 2008; Booth & Schaye 2009; Kurosawa & Proga 2009; Debuhr et al. 2011; Teyssier et al. 2011; Dubois et al. 2012). Our BH and AGN feedback implementation follows closely previous studies (Di Matteo et al. 2005, 2008; Springel et al. 2005a; Sijacki et al. 2007), with modifications to adopt the scheme to the moving-mesh and finite-volume technique. In what follows, BHs are represented by collisionless, massive sink particles that are created at early times, and subsequently grow in mass by gas accretion or BH mergers. We include three different forms of back-reaction of the accretion: thermal, mechanical and electromagnetic (EM) feedback. Thermal and mechanical feedback have already been considered in some previous large-scale cosmological simulations, whereas EM radiative feedback has typically been neglected (except for a few cases, see below). We here present a novel phenomenological treatment of this form of feedback.

2.6.1 BH seeding

Our seeding strategy for BHs in cosmological simulations follows previous work. We regularly run an FoF group finder and assign seed BHs to massive enough FoF groups (see Sijacki et al. 2007; Di Matteo et al. 2008 for details). We set the seeding mass threshold for FoF groups to 5 × 10¹⁰ h⁻¹ M_⊙. For an FoF group that is selected to acquire a seed BH, we turn the gas cell with the highest density in the group into a BH sink particle. The dynamical mass of the particle is set to the cell mass. This dynamical mass is typically significantly higher than the real expected seed mass of the BH, which we set in the following to 10⁵ h⁻¹ M_⊙, as in many previous studies. To address this problem, we follow Springel et al. (2005a) and track the BH mass as a separate subgrid variable. This internal BH mass smoothly increases according to the estimated BH accretion rate |$\dot{M}_{\rm BH}$| (see below), such that the value of this internal mass represents the BH mass for the case where accretion can be fully resolved and mass is not discretized.

We allow the BH sink particle to drain mass continuously from the primary cell of the BH particle. Specifically, during each time step, an active BH particle removes the mass |${\Delta M = (1-\epsilon _\mathrm{r}) \, \dot{M}_{\rm BH} \, \Delta t}$| from its primary cell, where Δt is the current time step and ϵ_r is the radiative efficiency (see below). BH accretion rates are typically very small compared to the actual cell masses of the primary cell such that ΔM is usually significantly smaller than the mass of the primary cell. In the rare event that this is not the case, we added a ‘bucket mechanism’ which collects ΔM in an accretion mass bucket to be then handled with a small time delay. During subsequent time steps, this bucket is then emptied such that we always accrete the correct amount of mass. Specifically, if ΔM is larger than 90 per cent of the primary cell mass, we remove this 90 per cent of the cell mass and put the remaining mass of ΔM into the bucket. During the next accretion event of the BH, we then try to remove the current ΔM from the cell plus the mass that is currently in the bucket. Here again, we allow at maximum a removal of 90 per cent of the cell mass. This procedure guarantees a continuous accretion scheme, where the internal and dynamical masses of the BH particles grow at the same rate, which is different from the original stochastic accretion implementation used in previous SPH studies of this model (Springel et al. 2005a). We note that the draining scheme does not change the real dynamical mass of BH sink particles until the internal mass is equal to or larger than the initially assigned dynamical mass. Also, only at that point do we actually drain mass from the primary cell, thereby maintaining mass conservation during the whole evolution.

2.6.2 BH growth

However, when no star-forming gas is present in the immediate vicinity of the BH (and the quasar accretion is hence in a ‘low state’), this prescription is expected to significantly overestimate the true accretion rate. In particular, this will happen when the BH has grown to a large mass and is embedded in comparatively low-density background gas. Then the residual BH accretion we estimate with equation (21), together with our quasar-mode feedback scheme (see below), can create a hot, low-density bubble around the BH, with a pressure that matches the background gas pressure P_ext at the centre of the halo. The formation of such a bubble with density well below the SF density threshold would clearly be an unrealistic artefact of our subgrid model. We address this as follows.

BHs do not only grow through gas accretion, but also when they merge with other BHs. BHs merge once they are within their ‘feedback radius’ (see below). Here we also do not consider the relative velocities of BHs for the same reason that we neglect the velocity term in the accretion rate.

2.6.3 Quasar- and radio-mode feedback

For AGN feedback, we use a two-state model as used in previous studies (see Springel et al. 2005a; Sijacki et al. 2007 for details). This model distinguishes between quasar-mode feedback (high BH accretion rates) and radio-mode feedback (low BH accretion rates). Quasar-mode AGN feedback is modelled by assuming that a fraction (ϵ_f) of the radiative energy released by the accreted gas couples thermally to nearby gas within a radius that contains a pre-calculated mass scale. Significant quasar activity requires high densities of relatively cold gas around the BH, supplied through large-scale inflows during galaxy mergers. We will in the following assume a radiative efficiency of ϵ_r ∼ 0.1–0.2 (Shakura & Sunyaev 1973; Yu & Tremaine 2002). For quasar-mode AGN feedback, we follow the same approach as for the enrichment scheme presented above and distribute energy over a number of BH neighbour cells such that a pre-defined mass is enclosed within the thermal ‘feedback radius’. We set the mass over which to spread the BH thermal energy to the same mass that is also used for the enrichment scheme, which again is typically set to a multiple of the target mass (m_target) of our (de-)refinement scheme. For convergence studies, we hold this mass fixed if we change resolution, i.e. the feedback energy is always distributed over a fixed physical mass, following the same strategy as for the enrichment scheme.

For low-activity states of the BH, we consider a form of mechanical radio-mode AGN feedback following Sijacki et al. (2007). Quasar- and radio-mode feedback are distinguished based on the BH accretion rate. For BH accretion rates below a fraction χ_radio of the Eddington rate, we assume that feedback operates in radio mode, where AGN jets inflate hot, buoyantly rising bubbles in the surrounding halo atmosphere. The duty cycle of bubble injection, energy content of the bubbles as well as their initial size are estimated based on the BH accretion rate and the current BH mass. The duty cycle is coupled to the mass growth of the BH such that energy is released once the BH has increased its mass by a factor of δ_BH. For the bubble properties, we assume the scaling relations of Sijacki et al. (2007), with R_bub,0 = 50 kpc, E_bub,0 = 10⁶⁰ erg and |$\rho _\mathrm{ICM,0}=10^{4}{\,\rm M_{\odot }}\,{\rm kpc}^{-3}$|⁠, which we will keep fixed in the following. For the scaling of the distance over which bubbles are placed, we use a normalization value of 100 kpc, which we will also keep fixed in the following. We assume the radio-mode feedback efficiency provided by the bubbles to be ϵ_m × ϵ_r of the accreted rest mass energy.

2.6.4 Radiative AGN feedback

In Fig. 1, we show net cooling rates for gases of different metallicities which are exposed to different AGN radiation sources characterized by their ionization parameter |$\zeta ^{\rm AGN}=L_{\rm bol}^{\rm AGN}/(r^2 n_{\rm H})$|⁠, where |$L_{\rm bol}^{\rm AGN}=(1-\epsilon _\mathrm{f}) \, \tilde{\epsilon }_\mathrm{r} \, \dot{M}_{\rm BH} c^2$| is the bolometric luminosity and r is the distance from the plasma cell to the AGN radiation source. We note that |$J^{\rm AGN}\equiv L_{\rm bol}^{\rm AGN}/r^2$| is the bolometric intensity that we tabulate in the cooling table to look up the metal-line cooling. As Fig. 1 demonstrates, a simple Z-scaling of the rates does not work in this regime. We hence need to tabulate metal-line cooling rates also as a function of Z. Specifically, we extend the grid of net cooling rates by [log (Z/Z_⊙), log (4πJ^AGN/( erg s^{− 1} cm^{− 2}))], with − 5 < log (4πJ^AGN/( erg s^{− 1} cm^{− 2})) < 5 sampled equally spaced at 11 bins, and six different metallicities log (Z/Z_⊙) = −4, −3, −2, −1, 0, 1. At each redshift, we also add grid points for a pure UVB which we use for heating and cooling if there is no nearby AGN. Although Fig. 1 clearly demonstrates that the AGN radiation field can strongly influence net cooling rates, we note that this occurs only for a short interval of cosmic time for each BH once it is in its maximum accretion phase with an accretion rate close to Eddington. This implies that the impact of EM AGN feedback is strongest during quasar activity, but is limited when AGN enter the low accretion radio-mode state.

Figure 1.

View largeDownload slide

Influence of the AGN radiation field (‘EM feedback’) on the net cooling rate for gases of different metallicities, in the range log [Z/Z_⊙] = −2.0, −1.0, 0.0, 1.0 (lowest to highest curve). The AGN ionization parameter is negligibly small in the left-hand panel (⁠|$\zeta ^{\rm AGN}_{\rm thresh}=10^{-1}\,{\rm erg}\,{\rm s}^{-1}\,{\rm cm}$|⁠), but substantial in the right-hand panel (⁠|$\zeta ^{\rm AGN}_{\rm thresh}=10^{3}\,{\rm erg}\,{\rm s}^{-1}\,{\rm cm}$|⁠). Clearly, a high ionization parameter leads to a strong cooling suppression and significant additional heating, as seen in the net cooling rates in the right-hand panel. We note that such large ionization parameters require high BH accretion rates close to Eddington. Therefore, halo gas is exposed to such a radiation field typically only for a rather short period of cosmic time during quasar activity. For lower BH accretion rates, the radiation field is essentially too small to produce large ionization parameters for a significant fraction of the halo gas.

Figure 1.

View largeDownload slide

Influence of the AGN radiation field (‘EM feedback’) on the net cooling rate for gases of different metallicities, in the range log [Z/Z_⊙] = −2.0, −1.0, 0.0, 1.0 (lowest to highest curve). The AGN ionization parameter is negligibly small in the left-hand panel (⁠|$\zeta ^{\rm AGN}_{\rm thresh}=10^{-1}\,{\rm erg}\,{\rm s}^{-1}\,{\rm cm}$|⁠), but substantial in the right-hand panel (⁠|$\zeta ^{\rm AGN}_{\rm thresh}=10^{3}\,{\rm erg}\,{\rm s}^{-1}\,{\rm cm}$|⁠). Clearly, a high ionization parameter leads to a strong cooling suppression and significant additional heating, as seen in the net cooling rates in the right-hand panel. We note that such large ionization parameters require high BH accretion rates close to Eddington. Therefore, halo gas is exposed to such a radiation field typically only for a rather short period of cosmic time during quasar activity. For lower BH accretion rates, the radiation field is essentially too small to produce large ionization parameters for a significant fraction of the halo gas.

The intensity assignment is performed by looping over all BHs and calculating their bolometric luminosities based on the current mass accretion rates. We first calculate a search radius for each BH. For this we derive the BH's bolometric obscured luminosity and solve for the radius within which the ionization parameter is larger than |$\zeta ^{\rm AGN}_{\rm thresh}=10^{-3}\,{\rm erg}\,{\rm s}^{-1}\,{\rm cm}$|⁠, assuming gas with mean hydrogen number density. We then calculate for all cells within this radius the incoming AGN intensity. To avoid excessively large search radii, we impose an upper limit of three times the virial radius to this search radius. We have checked that our results are not sensitive to this radius once it is chosen to be of the order of the virial radius or larger. Overlapping AGN contributions of several BHs are handled by summing up their bolometric intensities at the cell's position. Once we know all local intensities, we can calculate the net cooling rate for each cell. We do not include speed-of-light delay effects in the propagation from the AGN, which is a reasonable approximation since we limit the sphere of influence of the radiation field as described above, and do not allow the radiation field to propagate over arbitrarily large distances through the simulation volume.

3 ANALYSIS TECHNIQUES

3.1 Tracer particles

The pseudo-Lagrangian nature of SPH is convenient for approximately tracing mass elements in simulations. [However, we emphasize that this ‘advantage’ of SPH comes at the expense of an inaccurate integration of the mass continuity equation (see, e.g., Vogelsberger et al. 2012).] On the other hand, the Eulerian nature of finite-volume schemes does not directly allow this because the mass continuity equation is integrated correctly on the resolution scale. One approach to address this ‘disadvantage’ of grid codes is to use an additional passive tracer particle species which is advected along with the gas by inferring a tracer's velocity by interpolating the calculated hydrodynamic velocity field.

arepo employs a quasi-Lagrangian scheme, i.e. it is neither Eulerian nor strictly Lagrangian since mass can be advected through cell boundaries as in an ordinary finite-volume approach. (We note here that SPH is not actually Lagrangian either, at best ‘pseudo-Lagrangian’, as the remapping that implicitly occurs when local neighbourhood relations of particles change is ignored in this approach.) This quasi-Lagrangian nature of arepo implies that tracing the evolution of ‘closed-box’ fluid elements is not possible a priori, similar to the situation in other finite-volume schemes. We have therefore realized two different passive tracer techniques in arepo (Genel et al. 2013). First, we implemented the velocity field approach, which ties the tracers to the reconstructed and interpolated fluid velocity field. Second, we developed a new scheme, which determines the exchange of tracers between cells in a Monte Carlo fashion based on the actual mass exchanges between them, which is given by the Riemann solver. We briefly describe the main characteristics of these schemes in the following two subsections.

3.1.1 Velocity field tracer particles

Velocity field tracer particles move according to the piece-wise linear reconstruction of the velocity field. At each time step, we make a lookup of the closest mesh-generating point for every active tracer particle, which provides the Voronoi cell that this tracer resides in. We then use the fluid velocity field gradients of that cell to interpolate the velocity field to the tracer position. The gradient information is readily available for the cells, because it is already calculated for the MUSCL-Hancock step in the finite-volume solver. Once a new velocity is assigned to each active tracer particle, we drift them according to their individual time steps. Tracer particles inherit the time step from the cells they fall in, which makes their time integration adaptive and consistent with the dynamics of the underlying cell time step hierarchy.

Velocity field tracer particles try to follow the stream lines of the underlying flow field and are essentially noise free. Moving them is computationally very cheap; the only moderate computational effort lies in the required closest cell lookup, which however can be realized very efficiently through a tree search with an initial search radius guessed based on the nearest cell distance of the last search. However, we have found that this tracer particle approach does not follow the gas mass properly. Since the fluid is evolved based on the local solutions to the Riemann problems across each cell interface, simply using the reconstructed velocity field to advect the tracers does not lead to consistent evolution between the tracers and the fluid. Moreover, the fluid is averaged in the cell and reconstructed after it is evolved, while those steps are not performed on the tracers.

The inconsistency between the flow of the fluid and that of the tracers leads in cosmological simulations to large biases (Genel et al. 2013). For example, the tracer density profiles of haloes deviate significantly from the actual gas density profiles. This effect was already found in other studies, but was not interpreted appropriately. For example, Price & Federrath (2010) studied turbulence with the velocity field tracer implementation of the flash code and found that tracer particles tend to clump on subresolution scales in turbulence simulations.

3.1.2 Monte Carlo tracers

Instead of relating the tracer evolution directly to the velocity field, we can also link them to the mass exchange between cells. The basic idea of this approach, presented and studied in Genel et al. (2013), is to attach a population of N_i ≥ 0 Monte Carlo tracers to a computational Voronoi cell i. Based on the finite-volume update solution for this cell, we know the mass fluxes through each cell face during the time step. We can then probabilistically sample the transfer of tracer particles from one cell to the other. This results in a Monte Carlo sampling of the underlying gas mass fluxes over the computational domain, which, by construction, does not suffer from the bias effect. Such a scheme can be easily inlined in the finite-volume calculations, which loop over all Voronoi cell faces. This face list is constructed by the tessellation engine. During the interface loop, we keep track of the current number of tracers per cell and of the current total mass in each cell. We only consider outgoing mass fluxes ΔM_{i, j} < 0 from cell i to cell j. Furthermore, we keep track only of the reduced mass |$\widetilde{M}_i$| of each cell, which is updated for each outgoing flux, but not for in-going fluxes since the tracer particle exchange is done for all outgoing fluxes starting from a given cell. In-going fluxes into cell i are symmetrically treated by the outgoing fluxes of cell j. The probability for a tracer to leave cell i and go into cell j is then given by |$p_{i, j}^{\rm flux} = \Delta M_{i, j}/\widetilde{M}_i$|⁠. To decide whether a tracer should leave a cell, we draw a random number x_α ∈ U(0, 1) for each tracer α ≤ N_i of cell i. The tracer is put into cell j if |$x_\alpha < p_{i, j}^{\rm flux}$|⁠.

In such a Monte Carlo-based approach, tracers have no phase-space coordinates within the cell, corresponding to the assumption that they are always uniformly mixed within a cell. We therefore store them simply in a globally distributed linked list where each list entry has a tracer ID, a set of tracked fluid properties, and a pointer to the next and the previous tracer in the list. Each fluid cell then needs only a pointer to the first tracer associated with it. Tracer exchanges between cells can be implemented as operations on this global linked list.

Another advantage of Monte Carlo tracers is that any sort of mass transfer in or out of the gas cells can be modelled probabilistically in that scheme. This can be used, for example, to follow mesh (de-)refinement operations, but also for treating the formation of stars, BHs or wind particles. We exploit this in our implementation and have added, in addition to the tracers associated with gas cells, star, wind and BH tracers to keep track of all baryonic mass exchanges. For example, this allows us to quantify precisely (albeit with some Monte Carlo noise that can be reduced by averaging over many cells) what fraction of gas in a cell has previously been part of a star or a wind, or exactly which gas was swallowed by BHs.

3.2 On-the-fly volume rendering

Efficient visualization of complex simulation data sets becomes increasingly important since it can provide insights into the physical mechanisms at play. Both for SPH- and AMR-based simulations, a large number of visualization toolkits exist. For example, splash (Price 2007) or splotch (Dolag et al. 2011) provides easy-to-use tools to visualize various aspects of astrophysical SPH simulation data. AMR data can be efficiently visualized by many different software packages, including mayavi, visit, paraview and yt (Turk et al. 2011). Tools for visualizing an unstructured Voronoi mesh like the one used by arepo are much less common and hardly exist in comparison.

The actual ray traversing is done by exploiting mesh connectivity information which is available through the Voronoi mesh implementation in arepo. We also use this connectivity for a smoothing operation that we optionally apply during the render process. Especially in regions of low density, the number of Voronoi cells sampling the gas distribution can be quite low. Visually, this can then lead to discrete transitions from one cell to the other, which can show up in the rendered images, similar to the familiar block structure showing up in many images of AMR simulations. To reduce this effect, we have implemented a smoothing procedure, which averages cell values based on the values of the nearest cells which are connected to the current ray-casting cell. Specifically, we first make one loop over all neighbouring cells and calculate the maximum distance of their mesh-generating points to the current ray location in the primary cell. We then take this radius as a smoothing length for a tophat-kernel-based cell value estimate, where we weight the cell values that are going to be rendered according to the cell volume.

3.3 Stellar population synthesis models

To compare our galaxy formation model predictions with observations, we need to transform stellar information of our simulated galaxies into photometric properties. Stellar population synthesis models (e.g. Leitherer et al. 1999; Bruzual & Charlot 2003; Le Borgne et al. 2004) provide a convenient way to associate the stars in our simulations with observable spectra or broad-band luminosities. We include this in our simulations by assigning to each star particle a series of broad-band luminosities based on the Bruzual & Charlot (2003) catalogues after taking into account the star particle's age, mass and metallicity. Currently, we choose to include the U, V, B, K, g, r, i and z bands; however, this could be extended easily to include any other spectral information which is accurately tabulated in the stellar population synthesis models.

The assigned stellar luminosities are not currently used in any part of the dynamical evolution, nor do we include radiative transfer or dust attenuation. However, adding together the luminosity contributions from all star particles in a given subhalo allows us to construct estimates of the galaxy luminosities in several bands for direct comparison against observations. We construct broad-band luminosities for all galaxies and store the values in the group catalogues to simplify basic post-processing analysis.

4 COSMOLOGICAL SIMULATIONS

4.1 Simulations

We will present in this section the first cosmological arepo simulations including the newly added physics described above. In the following, we adopt the cosmological parameters Ω_m0 = 0.27, Ω_Λ0 = 0.73, Ω_b0 = 0.0456, σ₈ = 0.81 and H₀ = 100 h km s⁻¹ Mpc⁻¹ = 70.4 km s⁻¹ Mpc⁻¹ (h = 0.704). We create realizations of this cosmology in periodic boxes with a side length of 25 h⁻¹ Mpc. Initial conditions are generated at z = 127 based on a linear power spectrum made by camb, with gas particles/cells added to the initial conditions by splitting each original particle into a DM and gas particle/cell pair, displacing them with respect to each other such that two interleaved grids are formed, keeping the centre of mass of each pair fixed.

Our fiducial parameters for the different physical processes are summarized in Table 1. These parameters are physically plausible and were chosen such that they provide a reasonable fit to most key observables at z = 0, as we discuss below. Most importantly, the feedback parameters of our fiducial model are chosen such that they reasonably well reproduce the stellar mass function and the stellar mass–halo mass relation. We also explore a few modifications of this fiducial model, where we focus on differences in the stellar and AGN feedback processes. All our simulations are summarized in Table 2. We specify two Plummer-equivalent gravitational softening length values. For DM particles, we use a fixed comoving softening length (second value). For baryonic collisionless particles (stars and BHs), we further assume a maximum physical softening length (first value), which limits the growth of the physical softening length. Gas cells use an adaptive softening length tied to the cell radius, limited by a floor. This floor is set to the same value as for the other baryonic particles (i.e. stars and BHs). In the table, we also specify the simulation volume and the mass resolution for DM and baryons. Finally, the last column describes the physics characteristics of the specific model with respect to our fiducial model. We employ a (de-)refinement scheme which keeps the cell masses close to a specified target mass m_target, which results in a nearly constant number of cells. This scheme is identical to that used in Vogelsberger et al. (2012), and we also follow the same mesh regularization strategy.

Table 2 includes three simulations with the fiducial setup presented in Table 1: L25n128, L25n256 and L25n512. These three simulations serve as a resolution study of our fiducial physics setup to estimate convergence. The other simulations explore modifications of the fiducial setup. Each of these modified runs typically varies only one of the feedback parameters of our fiducial setup and explores it at the intermediate resolution equivalent to the L25n256 simulation, which is sufficient to obtain reasonably converged results and computationally efficient enough to explore different feedback settings quickly. Each of the modified parameters is increased or decreased by a factor of 2 compared to the fiducial value. We explore here only a rather limited number of parameters. For example, we do not change the quasar-mode feedback strength since we are mainly interested in the stellar properties of haloes, and quasar-mode feedback mainly regulates BH growth (see below). Furthermore, we do not alter the wind model, i.e. we only vary the parametrization of the non-local energy-driven winds presented in Table 1, but we do not explore, for example, local or momentum-driven winds. For the sake of brevity, we here also do not explore a large number of unrealistic setups like models without metal-line cooling or without stellar mass-loss. While such test simulations can be interesting for studying the specific implications of these physical effects, our goal is here to bracket the feedback parameters required to obtain an acceptable physical match to observations; hence, we explore only models with complete physics in this subsection. However, we include one simulation (‘no feedback’), which does not include stellar and AGN feedback except for the implicit ISM pressurization through SNII feedback as implemented in our eEOS for star-forming ISM gas. We consider this ‘no feedback’ simulation to demonstrate that such a model gives a very poor fit to observations, highlighting the difference with our fiducial model which includes the full feedback physics. We stress that this ‘no feedback’ simulation contains a weak form of stellar feedback in the form of the ISM pressurization. To understand the importance of the different physical processes in our model, we will present unrealistic simulations in Section 4.3 below. There we deactivate certain physical processes to understand better how they shape our z = 0 results. Clearly, such simulations are not expected to match observational constraints, which is why we do not discuss them here.

We will show below that the fiducial parameters presented in Table 1 produce a galaxy population in good agreement with various observational probes at z = 0 in the local Universe. We stress that some of these feedback parameters differ significantly from those previously used in SPH simulations with similar physics models. Most importantly, the AGN radio-mode feedback is more energetic than in our previous simulations (e.g. Sijacki et al. 2007; Puchwein, Sijacki & Springel 2008). This can be seen directly by the value for the radio-mode feedback energy fraction, which is about a factor of 2 larger than in previous SPH studies. Also, the accretion rate threshold for switching between radio- and quasar-mode feedback is five times larger than that in previous simulations. This means that more of the BH accretion energy goes into radio-mode feedback compared to previous work, where χ_radio = 0.01 was typically adopted. As a result, more feedback energy is deposited or displaced from the central BH in a bursty fashion rather than smoothly heating the central region around the BH. This automatically leads to a more efficient suppression of SF in massive haloes. We found that this combination of an increased radio-mode feedback factor and a higher accretion rate threshold for the quasar mode is necessary to reproduce the observed stellar masses in massive haloes. We also note that the amount of SN energy for stellar feedback is quite high in our fiducial setup. In fact, it is three times larger than what is nominally assumed for SNII events (10⁵¹ erg). We argue that this is still physically plausible for a couple of reasons. First, 10⁵¹ erg is just a canonical value which is in fact uncertain. Second, recent work has shown that accounting for early stellar feedback seems crucial for regulating SF in low-mass systems (Kannan et al. 2013; Stinson et al. 2013). In our model, we include this contribution in the wind energy budget for launching kinetic wind particles.

The revised settings of the feedback model are required to match the z = 0 observations reasonably well. Less effective feedback leads to a significant overproduction of stellar mass at both the faint and bright ends. The need for even stronger feedback than that in SPH was already pointed out in Vogelsberger et al. (2012), where it was demonstrated that arepo shows enhanced gas cooling compared to previous SPH simulations, which we primarily attribute to spurious viscous heating effects in SPH and an inaccurate treatment of subsonic turbulence (Bauer & Springel 2012). This effect is strongest for high-mass systems where a quasi-hydrostatic hot atmosphere can form. SPH simulations show here effectively a form of ‘numerical quenching’, reducing the amount of AGN feedback required to regulate the amount of stars forming in these systems. Since arepo does not suffer from this problem (Vogelsberger et al. 2012), we are not helped by ‘numerical quenching’ and instead need to offset the enhanced cooling in large haloes at late times physically by stronger forms of stellar and AGN feedback. We note that this form of quenching can also be seen in Springel & Hernquist (2003a,b), where the cosmic SFR density could be reproduced without the inclusion of AGN feedback. This is possible because numerical quenching prevents SF in massive haloes.

4.2 Properties of the simulated galaxy population

In the following subsections, we present some first comparisons of our simulation results to observational data, where we mainly focus on the local Universe (z = 0). This serves primarily as a validation of our galaxy formation model and an identification of acceptable parameter settings. An initial study of the redshift evolution of the predicted galaxy properties is given in Torrey et al. (2013), while an in-depth analysis will be provided in forthcoming studies based on larger and higher resolution simulations.

We identify haloes, subhaloes and galaxies based on the subfind algorithm (Springel et al. 2001; Dolag et al. 2009). In the following, we often refer to the stellar galaxy mass of haloes and subhaloes. For definiteness, we define this mass as the total gravitationally bound stellar mass that is contained within twice the stellar half-mass radius of each subfind (sub)halo. This definition of stellar mass does not significantly differ from the total stellar mass for low-mass systems, but it excludes some of the intracluster light stars for massive systems. We have checked this definition against surface brightness cuts in different bands and find it to give very similar results as such more detailed methods of excluding intracluster light.

4.2.1 General gas structure

A first qualitative view of the simulations is presented in Fig. 2, where we show gas density (left-hand panels), gas temperature (middle panels) and gas metallicity (right-hand panels) projections of the highest resolution L25n512 simulation with our fiducial physics setup at three redshifts z = 2, 1, 0 (top to bottom). Each panel is 25 h⁻¹ Mpc on a side and has a projection thickness of 1 h⁻¹ Mpc. At z = 2, some haloes show outflows generated mainly by winds associated with stellar feedback leading to IGM enrichment. The more dramatic heating effects at late times are largely caused by strong radio-mode AGN feedback in massive systems. Both stellar and AGN feedback lead to a significant enrichment of the IGM, as can be seen in the right-hand panels. Furthermore, the AGN feedback also alters the density structure of the gas at z = 0.

Figure 2.

View largeDownload slide

Gas density (left-hand panel), gas temperature (middle panel) and gas metallicity (right-hand panel) projections of the L25n512 simulation. Each panel is 25 h⁻¹ Mpc on a side and has a thickness of 1 h⁻¹ Mpc. We show the fields at three different redshifts z = 0, 1 and 2. At z = 2 some haloes show outflows generated mainly by strong winds through stellar SNII feedback. The more dramatic heating effects at late times (z = 0, 1) are largely caused by strong radio-mode AGN feedback. Both stellar and AGN feedback lead to a significant enrichment of the IGM as can be seen in the metallicity projections. Strong AGN feedback also alters the density structure of the gas at z = 0 (see the left-hand panel).

Figure 2.

View largeDownload slide

Gas density (left-hand panel), gas temperature (middle panel) and gas metallicity (right-hand panel) projections of the L25n512 simulation. Each panel is 25 h⁻¹ Mpc on a side and has a thickness of 1 h⁻¹ Mpc. We show the fields at three different redshifts z = 0, 1 and 2. At z = 2 some haloes show outflows generated mainly by strong winds through stellar SNII feedback. The more dramatic heating effects at late times (z = 0, 1) are largely caused by strong radio-mode AGN feedback. Both stellar and AGN feedback lead to a significant enrichment of the IGM as can be seen in the metallicity projections. Strong AGN feedback also alters the density structure of the gas at z = 0 (see the left-hand panel).

To demonstrate the impact of stellar and AGN feedback more clearly, we compare in Fig. 3 the L25n256 simulation with the ‘no feedback’ simulation, which does not include stellar and AGN feedback. We note that the colour scales in Figs 2 and 3 are not the same and have been adapted for each resolution. Fig. 3 clearly demonstrates that stellar and AGN feedback strongly affect the thermodynamic properties of the gas as can be seen from the density and temperature maps, which differ significantly between the two runs. Furthermore, the metal distribution is very different. In fact, the simulation without any feedback does not lead to any significant IGM enrichment as metals are locked up efficiently in galaxies at the centres of haloes.

Figure 3.

View largeDownload slide

Gas density (left-hand panel), gas temperature (middle panel) and gas metallicity (right-hand panel) projections of the L25n256 simulation (top panels) and ‘no feedback’ simulation (bottom panels) at z = 0. Each panel is 25 h⁻¹ Mpc on a side and has a thickness of 1 h⁻¹ Mpc (same as in Fig. 2). Clearly, the addition of more efficient cooling and stellar/AGN feedback strongly affects the thermodynamic properties of the gas. Most dramatically, the metal distribution is altered through the additional feedback enriching the IGM. We note that the colour scale used here is different from the one in Fig. 2 as it has been adapted to the different simulation resolution in both figures.

Figure 3.

View largeDownload slide

Gas density (left-hand panel), gas temperature (middle panel) and gas metallicity (right-hand panel) projections of the L25n256 simulation (top panels) and ‘no feedback’ simulation (bottom panels) at z = 0. Each panel is 25 h⁻¹ Mpc on a side and has a thickness of 1 h⁻¹ Mpc (same as in Fig. 2). Clearly, the addition of more efficient cooling and stellar/AGN feedback strongly affects the thermodynamic properties of the gas. Most dramatically, the metal distribution is altered through the additional feedback enriching the IGM. We note that the colour scale used here is different from the one in Fig. 2 as it has been adapted to the different simulation resolution in both figures.

To give a better impression of the three-dimensional distribution of gas in the simulation volume, we show in Fig. 4 a volume-rendered gas temperature field of a large fraction of the L25n512 simulation volume.² Here, the newly implemented, inlined volume renderer has been used, as described above. We do not employ the smoothing procedure for this rendering. The render frustum has a near-field plane distance of 2 h⁻¹ Mpc and the far plane is at 20 h⁻¹ Mpc. The width of the near-field plane is 2 h⁻¹ Mpc with a height of 1.5 h⁻¹ Mpc. Hot regions are coloured in red and cold regions in blue. Galaxies are shown in bright white. Parts of the cosmic web and the hot halo atmosphere of one of the most massive haloes in the simulation volume can be seen. In Fig. 5, we show the same region of the simulation volume but render the gas metallicity instead. The hot gas atmosphere of the cluster is clearly significantly enriched. We note that the filaments, which are visible in the temperature rendering, are not clearly visible in Fig. 5.

Figure 4.

View largeDownload slide

Rendered gas temperature field of a large region of the simulation volume of L25n512 at z = 0. The render frustum has a near-field plane distance of 2 h⁻¹ Mpc and the far plane is at 20 h⁻¹ Mpc. The width of the near-field plane is 2 h⁻¹ Mpc with a height of 1.5 h⁻¹ Mpc. Hot regions are shown in red and cold regions in blue. Galaxies are shown in bright white. Part of the cosmic web and hot halo atmosphere of one of the most massive haloes in the simulation volume can be seen.

Figure 4.

View largeDownload slide

Rendered gas temperature field of a large region of the simulation volume of L25n512 at z = 0. The render frustum has a near-field plane distance of 2 h⁻¹ Mpc and the far plane is at 20 h⁻¹ Mpc. The width of the near-field plane is 2 h⁻¹ Mpc with a height of 1.5 h⁻¹ Mpc. Hot regions are shown in red and cold regions in blue. Galaxies are shown in bright white. Part of the cosmic web and hot halo atmosphere of one of the most massive haloes in the simulation volume can be seen.

Figure 5.

View largeDownload slide

Rendered gas metallicity field of a large region of the simulation volume of L25n512 at z = 0. The render frustum has a near-field plane distance of 2 h⁻¹ Mpc and the far plane is at 20 h⁻¹ Mpc. The width of the near-field plane is 2 h⁻¹ Mpc with a height of 1.5 h⁻¹ Mpc. The hot halo gas shown in Fig. 4 is significantly enriched.

Figure 5.

View largeDownload slide

Rendered gas metallicity field of a large region of the simulation volume of L25n512 at z = 0. The render frustum has a near-field plane distance of 2 h⁻¹ Mpc and the far plane is at 20 h⁻¹ Mpc. The width of the near-field plane is 2 h⁻¹ Mpc with a height of 1.5 h⁻¹ Mpc. The hot halo gas shown in Fig. 4 is significantly enriched.

4.2.2 Evolution of the cosmic SFR density

Observationally, it is now well established that the cosmic SFR density increased significantly from z ∼ 10 to z ∼ 2 and reached a peak value at around z ∼ 2 followed by a rapid decline (Lilly et al. 1996; Madau, Pozzetti & Dickinson 1998; Schiminovich et al. 2005; Hopkins & Beacom 2006; Bouwens et al. 2008). In Fig. 6, we show the simulated cosmic SFR density as a function of redshift. To compare our models against observational constraints, we have included a set of cosmic SFR density measurements from the literature, as compiled in Hopkins & Beacom (2006) and Behroozi et al. (2013). We show a convergence study in the left-hand panel, and an exploration of the different feedback models summarized in Table 2 in the right-hand panel. We will in the following subsections follow the same format and always explore numerical convergence and physical model variations side by side.

Figure 6.

View largeDownload slide

Cosmic SFR density as a function of redshift (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data are taken from Hopkins & Beacom (2006) and Behroozi, Wechsler & Conroy (2013). The left-hand panel demonstrates that the overall shape of the SFR evolution is well converged, with a minor normalization change towards higher SFRs with increasing resolution, but the peak location does not change with resolution. Doubling the wind velocity (‘faster winds’) leads to a strong suppression of the overall SFR. Increasing the overall wind energy (‘stronger winds’) leads to a clear overproduction of stars at late times. Variations in the AGN radio-mode feedback do not change the SFR density nearly as much as changes in the stellar feedback. Stronger radio-mode feedback leads to less SF at late times. The same is true for an increased radio-mode threshold, where more AGN feedback energy is channelled into the radio mode. Since this mode is more bursty, it is more efficient at reducing SF compared to the continuous quasar-mode AGN feedback, which mainly regulates the growth of BHs.

Figure 6.

View largeDownload slide

Cosmic SFR density as a function of redshift (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data are taken from Hopkins & Beacom (2006) and Behroozi, Wechsler & Conroy (2013). The left-hand panel demonstrates that the overall shape of the SFR evolution is well converged, with a minor normalization change towards higher SFRs with increasing resolution, but the peak location does not change with resolution. Doubling the wind velocity (‘faster winds’) leads to a strong suppression of the overall SFR. Increasing the overall wind energy (‘stronger winds’) leads to a clear overproduction of stars at late times. Variations in the AGN radio-mode feedback do not change the SFR density nearly as much as changes in the stellar feedback. Stronger radio-mode feedback leads to less SF at late times. The same is true for an increased radio-mode threshold, where more AGN feedback energy is channelled into the radio mode. Since this mode is more bursty, it is more efficient at reducing SF compared to the continuous quasar-mode AGN feedback, which mainly regulates the growth of BHs.

The left-hand panel of Fig. 6 demonstrates that the shape of the SFR evolution is well preserved for all three resolutions, with a minor normalization change towards higher SFRs with increasing resolution. The magnitude of this normalization offset between our two highest resolution runs is ≲1.5, which we consider to be reasonably well converged. The overall shape of the cosmic SFR density is largely determined by stellar feedback at high redshifts and radio-mode AGN feedback towards lower redshifts. We find that the rapid decline of the SFR towards lower redshifts can only be achieved through strong radio-mode AGN feedback (see also Schaye et al. 2010; van de Voort et al. 2011; Bower, Benson & Crain 2012, for example). In fact, in the absence of AGN feedback, but with the addition of stellar feedback, the cosmic SFR density would continue to rise beyond z = 0, since gas recycling and metal-line cooling can support a large amount of SF at late times (see also Schaye et al. 2010). We note in particular that the location of the peak of the SFR density does not change significantly with numerical resolution. The intermediate-resolution run L25n256 is already sufficient to give a reasonably reliable estimate of the SFR in our simulation volume. We can hence study the effects of different physics parametrizations by running simulations at this intermediate resolution, as listed in Table 2.

The comparison among these different feedback models is shown in the right-hand panel of Fig. 6. Our fiducial model (L25n256) provides clearly the best match to the observational data compared to any of the other model variations presented in Table 2. The figure also highlights that the impact of changes in the wind model can be rather dramatic. For example, doubling the wind velocity (‘faster winds’) leads to a strong suppression of the overall SFR, especially at late times. Such fast winds lead to a significant reduction of SF also for massive haloes as we show below (see also Schaye et al. 2010 for a similar finding). Increasing the overall wind energy (‘strong winds’) leads to a clear overproduction of stars at late times, because a very substantial fraction of the gas that is blown out of galaxies at early times falls back in at late times – fuelling further SF along a ‘wind accretion’ channel (see also Oppenheimer et al. 2010). The peak of the SFR density is also sensitive to the details of the wind and AGN parametrization. Interestingly, the changes of the SFR density with respect to variations in the AGN radio mode are less dramatic. As expected, ‘stronger radio’ feedback leads to more efficient suppression of SF at late times. Likewise, ‘weaker radio’ mode feedback increases SF at late times. Similar trends can be seen for changes in the radio threshold. A ‘higher radio threshold’ puts more energy into the radio mode and is therefore more efficient in suppressing SF in massive systems at late times. A ‘lower radio threshold’ on the other hand puts more AGN feedback energy into the quasar mode, which is less efficient in suppressing SF. Consequently, such a simulation produces more stars at late times. In fact, we find that quasar-mode feedback is mainly responsible for establishing the BH scaling relations, but does not significantly affect SF in massive systems. The ‘no feedback’ simulation, as expected, strongly overproduces the amount of stars at all times. However, the peak of the SFR occurs still at about the right time even for this simulation.

4.2.3 Stellar mass–halo mass relation

Figure 7.

View largeDownload slide

Stellar mass of central galaxies as a function of the total halo mass (M_{200, crit}) at z = 0 (left-hand panel: resolution study; right-hand panel: various feedback models). Different black lines show abundance matching results (Guo et al. 2010; Moster et al. 2010, 2013; Behroozi et al. 2013), which are extrapolated beyond the constrained regime. The solid coloured lines mark the median relations of the simulations, whereas the two-dimensional histograms indicate the distribution of the L25n512 (left-hand panel) and L25n256 (right-hand panel) results. The left-hand panel demonstrates that the stellar content of low-mass galaxies (i.e. M_⋆ < 10⁹ M_⊙) is not yet fully converged. But at somewhat higher galaxy masses (i.e. M_⋆ > 10⁹ M_⊙), we find convergence in the stellar mass–halo mass relationship. The overall shape of the relation and the turnover mass agree reasonably well with the results derived based on abundance matching techniques. The right-hand panel demonstrates that both stellar feedback at the faint end and AGN feedback at the massive end shape the stellar mass content. The wind speed has the most dramatic impact on the stellar mass content of haloes. ‘Faster winds’ reduce the amount of stellar mass substantially over a large range of halo masses. ‘Weaker winds’ clearly do not suppress SF efficiently enough at the faint end. The massive end is most sensitive to changes in the radio-mode accretion rate threshold, where more radio-mode AGN feedback leads to more efficient quenching of massive systems. The same effect can be achieved through a larger radio-mode feedback strength.

Figure 7.

View largeDownload slide

Stellar mass of central galaxies as a function of the total halo mass (M_{200, crit}) at z = 0 (left-hand panel: resolution study; right-hand panel: various feedback models). Different black lines show abundance matching results (Guo et al. 2010; Moster et al. 2010, 2013; Behroozi et al. 2013), which are extrapolated beyond the constrained regime. The solid coloured lines mark the median relations of the simulations, whereas the two-dimensional histograms indicate the distribution of the L25n512 (left-hand panel) and L25n256 (right-hand panel) results. The left-hand panel demonstrates that the stellar content of low-mass galaxies (i.e. M_⋆ < 10⁹ M_⊙) is not yet fully converged. But at somewhat higher galaxy masses (i.e. M_⋆ > 10⁹ M_⊙), we find convergence in the stellar mass–halo mass relationship. The overall shape of the relation and the turnover mass agree reasonably well with the results derived based on abundance matching techniques. The right-hand panel demonstrates that both stellar feedback at the faint end and AGN feedback at the massive end shape the stellar mass content. The wind speed has the most dramatic impact on the stellar mass content of haloes. ‘Faster winds’ reduce the amount of stellar mass substantially over a large range of halo masses. ‘Weaker winds’ clearly do not suppress SF efficiently enough at the faint end. The massive end is most sensitive to changes in the radio-mode accretion rate threshold, where more radio-mode AGN feedback leads to more efficient quenching of massive systems. The same effect can be achieved through a larger radio-mode feedback strength.

The left-hand panel of Fig. 7 demonstrates that low-mass galaxies (i.e. M_⋆ < 10⁹ M_⊙) experience a factor of ≲2 increase in their stellar masses as we increase the resolution. Although we do not expect similarly sized changes to the stellar masses of these galaxies with further increased resolution, the stellar masses of these low-mass galaxies are clearly not yet fully converged. Galaxies with somewhat higher masses (i.e. M_⋆ > 10⁹ M_⊙) are better resolved and show reasonable convergence in their stellar mass–halo mass relationship. As a result, the overall shape of the relation and the turnover mass are quite well converged and agree with the results derived based on abundance matching techniques. We note that the break in the stellar mass-to-halo mass ratio around |$M_\mathrm{200,crit}\sim 10^{12}{\,\rm M_{\odot }}$| can be reproduced only through a combination of stellar and AGN feedback. In fact, the main challenge in reproducing the abundance matching results is to correctly match the change of slope towards more massive systems. Reproducing the faint-end or massive-end slope alone is typically easier to achieve. Our simulation volume is too small to properly sample massive systems, but based on selected test calculations of clusters of galaxies using a zoom-in procedure, we have verified that our feedback model also reduces the stellar masses of these more massive systems to the observed levels.

The right-hand panel of Fig. 7 demonstrates that both stellar feedback at the faint end and AGN feedback at the massive end are needed to produce a stellar mass–halo mass relation that matches the derived relation from abundance matching reasonably well. In agreement with our findings for the cosmic SFR density, the wind speed has the most dramatic impact on the stellar mass content of haloes, with ‘faster winds’ being capable of substantially suppressing the stellar mass content of haloes over a wide mass range, also affecting more massive systems. ‘Weaker winds’ in contrast clearly do not reduce SF efficiently enough at the faint end. This can be seen from the ‘weaker wind’ curve which significantly overproduces the stellar mass for low-mass systems. On the other hand, ‘stronger winds’ lead to an excessive suppression of SF towards the low-mass end, and also an undershoot around the turnover point of the observed relation. The massive end is very sensitive to changes in the radio-mode accretion rate threshold. Lowering this value, and thereby putting more AGN feedback energy into the quasar mode, leads to an overproduction of stars in massive systems as can be seen from the ‘lower radio threshold’ curve. The same is true if we decrease the radio feedback factor by a factor of 2 (‘weaker radio’ curve). These findings agree with the conclusions drawn from the cosmic SFR density plots. As for the cosmic SFR density, the fiducial L25n256 model provides the best fit to the observational data. The ‘no feedback’ simulation clearly overproduces stars at all halo masses. Furthermore, the number of stars scales nearly linearly with the halo mass in this case (above |$M_\mathrm{200,crit}\sim 10^{10.5}{\,\rm M_{\odot }}$|⁠), and there is clearly no turnover towards higher masses. This implies that radio-mode AGN feedback is crucial for the quenching of SF at the massive end (as has been previously found by Bower et al. 2006; Cattaneo et al. 2006; Croton et al. 2006; Puchwein et al. 2008; Somerville et al. 2008; McCarthy et al. 2010; Teyssier et al. 2011; Dubois et al. 2013, for example).

We stress again that arepo does not experience the ‘numerical quenching’ of cooling in large haloes present in SPH simulations. This artificial quenching that occurs mainly in massive haloes helps to reduce the amount of SF in these systems (see Vogelsberger et al. 2012), which implies that even without explicit AGN feedback the stellar masses of massive haloes are pushed towards the observed relationship. Even though this is helpful in this sense, it is an undesirable effect as it is purely of numerical origin. In fact, it will then lead one to adopt physically incorrect settings for the radio-mode AGN feedback parameters.

4.2.4 Stellar mass function

The shape of the stellar mass function is determined by the underlying halo mass function convolved with the efficiency at which stars can form in these haloes. As shown in the previous subsection, SF in our simulations is most efficient for haloes with masses around |$M_\mathrm{200,crit}\sim 10^{12}{\,\rm M_{\odot }}$|⁠, in agreement with abundance matching results. SF in lower mass haloes is suppressed through stellar feedback, whereas SF in more massive systems is quenched by AGN feedback, as found in previous semi-analytic models (mainly Bower et al. 2006; Croton et al. 2006) and more recently in self-consistent hydrodynamical simulations (e.g. Puchwein & Springel 2013). Especially, the exponential suppression at higher masses requires strong AGN feedback.

In Fig. 8, we show the simulated z = 0 stellar mass function and compare it to observational data of Baldry et al. (2008) and Pérez-González et al. (2008), with correction factors to account for our Chabrier IMF. Here we include all subhaloes and treat the stellar component of each as a galaxy contributing to the stellar mass function. As discussed above, the stellar mass is measured within twice the stellar half-mass radius. The left-hand panel of Fig. 8 shows that the number density of galaxies at any given stellar mass increases marginally with resolution. This is directly tied to the increase in stellar mass for low-mass haloes, as discussed in the previous subsection. Consistently with our previous discussion, we find that the stellar mass function shows reasonable convergence for higher mass systems (i.e. M_* > 10⁹ M_⊙). But increasing the resolution allows an examination of progressively smaller galaxy masses. With our highest resolution simulation, we are able to compare to the observed stellar mass functions over quite a large range of masses. Although both functions differ in their detailed shape, we find that our simulated stellar mass functions are remarkably similar to the observations. This includes the slope of the low-mass end of the stellar mass function, which is shaped by our stellar wind model, as well as the sharp drop towards the massive end, which is due to the radio-mode AGN feedback. We emphasize that the energy scaling of our stellar wind model is crucial to produce a very flat stellar mass function towards lower mass systems. Note that the mass loading for low-mass systems can achieve very high values if an energy scaling is adopted, and it grows more rapidly towards small galaxy sizes than for a momentum scaling. A similar conclusion was reached by Puchwein & Springel (2013), who also argued for energy-driven winds to obtain a sufficiently flat stellar mass function for low-mass systems (but see Oppenheimer et al. 2010 for arguments in favour of momentum-driven winds).

Figure 8.

View largeDownload slide

Stellar mass function at z = 0 including all subhaloes (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data points are taken from Baldry, Glazebrook & Driver (2008) and Pérez-González et al. (2008). The left-hand panel shows that the stellar mass function is reasonably converged for higher mass systems (i.e. M_* > 10⁹ M_⊙). In the right-hand panel, the ‘faster wind’ simulation leads to a strong reduction of stellar mass in agreement with similar findings for other diagnostics like the SFR density and the stellar mass–halo mass relation. The ‘weaker winds’ simulation overproduces the number of galaxies at the faint end significantly, whereas the ‘strong winds’ suppress SF too much in these systems. A weaker radio-mode feedback leads to a less sharp drop-off of the stellar mass function towards the massive end. The same is true for a lower radio-mode accretion threshold since more energy goes then into the non-bursty quasar-mode AGN feedback, which is less effective at suppressing SF.

Figure 8.

View largeDownload slide

Stellar mass function at z = 0 including all subhaloes (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data points are taken from Baldry, Glazebrook & Driver (2008) and Pérez-González et al. (2008). The left-hand panel shows that the stellar mass function is reasonably converged for higher mass systems (i.e. M_* > 10⁹ M_⊙). In the right-hand panel, the ‘faster wind’ simulation leads to a strong reduction of stellar mass in agreement with similar findings for other diagnostics like the SFR density and the stellar mass–halo mass relation. The ‘weaker winds’ simulation overproduces the number of galaxies at the faint end significantly, whereas the ‘strong winds’ suppress SF too much in these systems. A weaker radio-mode feedback leads to a less sharp drop-off of the stellar mass function towards the massive end. The same is true for a lower radio-mode accretion threshold since more energy goes then into the non-bursty quasar-mode AGN feedback, which is less effective at suppressing SF.

In the right-hand panel of Fig. 8, we explore variations around our fiducial model. The ‘faster wind’ simulation leads to a strong reduction of stellar mass in agreement with similar findings for other diagnostics presented above. Such winds are capable of suppressing SF even in more massive haloes. The other models behave largely very similarly. The ‘weaker winds’ simulation overproduces the number of galaxies at the faint end significantly, whereas ‘stronger winds’ suppress SF too much in these systems. The ‘stronger winds’ also lead to an undershoot of the stellar mass function around the cut-off scale. These findings are consistent with the ones presented above for the stellar mass–halo mass relationship. A ‘weaker radio’ mode feedback leads to a less sharp drop-off of the stellar mass function towards the massive end. The same is true for a lower radio-mode accretion threshold since in that case more energy goes then into the non-bursty quasar-mode AGN feedback. The ‘no feedback’ simulation overproduces the number of galaxies at all stellar masses. Here the stellar mass function closely follows the shape of the halo mass function since most of the baryons in a halo can efficiently cool and form stars if not moderated by feedback. Note that the simulation without feedback does not show the strong exponential drop towards high galaxy masses. As for the stellar mass–halo mass relation, we hence argue that a quenching of SF at the massive halo end requires AGN feedback.

4.2.5 Stellar mass density

Integrating the cosmic SFR density leads to the overall stellar mass density in the Universe. In Fig. 9, we show the stellar mass density of the total simulation volume as a function of redshift, comparing it to estimates based on observational data (Dickinson et al. 2003; Fontana et al. 2006; Pozzetti et al. 2007; Pérez-González et al. 2008; Marchesini et al. 2009; Caputi et al. 2011; González et al. 2011; Mortlock et al. 2011). We note that we here include all stellar mass formed in the simulation volume, i.e. we do not restrict ourselves to stellar mass associated with (sub)haloes or that has formed within a certain radial range in a halo. The reason for this is to make this diagnostic consistent with the cosmic SFR density (see Fig. 6), which also includes all reservoirs of stellar mass formed without restriction.

Figure 9.

View largeDownload slide

Stellar mass density as a function of redshift for the various models (left-hand panel: resolution study; right-hand panel: different feedback models). We compare the simulation results to observational data taken from Dickinson et al. (2003), Fontana et al. (2006), Pozzetti et al. (2007), Pérez-González et al. (2008), Marchesini et al. (2009), Mortlock et al. (2011), González et al. (2011) and Caputi et al. (2011). The convergence is reasonable although the L25n512 simulation is clearly not yet fully converged (left-hand panel). But both L25n256 and L25n512 have a stellar mass density which is consistent with the observations. The ‘no feedback’ simulation strongly overproduces the amount of stellar mass, whereas the ‘faster winds’ strongly underproduce it. All the other feedback variations are quite close to the fiducial model, which provides the best fit to the data. Changes in the radio-mode AGN feedback only affect the late-time evolution of the stellar mass density.

Figure 9.

View largeDownload slide

Stellar mass density as a function of redshift for the various models (left-hand panel: resolution study; right-hand panel: different feedback models). We compare the simulation results to observational data taken from Dickinson et al. (2003), Fontana et al. (2006), Pozzetti et al. (2007), Pérez-González et al. (2008), Marchesini et al. (2009), Mortlock et al. (2011), González et al. (2011) and Caputi et al. (2011). The convergence is reasonable although the L25n512 simulation is clearly not yet fully converged (left-hand panel). But both L25n256 and L25n512 have a stellar mass density which is consistent with the observations. The ‘no feedback’ simulation strongly overproduces the amount of stellar mass, whereas the ‘faster winds’ strongly underproduce it. All the other feedback variations are quite close to the fiducial model, which provides the best fit to the data. Changes in the radio-mode AGN feedback only affect the late-time evolution of the stellar mass density.

Our fiducial model reproduces the observations reasonably well at all redshifts. As already seen above, the ‘faster winds’ suppress SF too strongly and therefore also significantly underpredict the total stellar mass density at nearly all redshifts but especially towards lower redshifts. Changes in the radio mode lead, as expected, to differences in the late-time behaviour since the radio mode becomes dominant only at late times once more massive haloes with lower BH accretion rates form. For the radio mode, we also find the anticipated behaviour: a ‘stronger radio’ mode reduces the total stellar mass density at late times, whereas a ‘weaker radio’ mode leads to an increase. A similar trend can be seen for variations of the radio-mode accretion threshold. A higher threshold suppresses late-time SF compared to our fiducial model, which is in agreement with the findings above. The ‘no feedback’ simulation strongly overproduces the number of stars.

The left-hand panel of Fig. 9 demonstrates that the convergence of the overall amount of stellar mass is reasonable although the L25n512 simulation is clearly not yet fully converged. But both L25n256 and L25n512 have a stellar mass density at z = 0 which is consistent with the observations. The non-converged L25n128 simulation clearly fails in this respect and has far too few stars at z = 0 since it does not have adequate mass resolution. We note that such a low resolution does not even resolve the halo mass function properly for our simulation volume. It is therefore not surprising that the stellar mass is substantially underpredicted in that case.

4.2.6 Tully–Fisher relation

Figure 10.

View largeDownload slide

Tully–Fisher relation at z = 0 (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data are taken from Bell & de Jong (2001) and Reyes et al. (2012). The solid lines mark the simulation median relations, whereas the two-dimensional histograms indicate the distribution of the L25n512 result (left-hand panel) and L25n256 result (right-hand panel). The runs converge towards a relation which is very similar to the marked observed relation with the slope and normalization of the simulated Tully–Fisher relation falling almost directly on top of the two marked observational Tully–Fisher relations. Most variations of the feedback model parameter choices lead to only minor changes in the resulting Tully–Fisher relation. The exception is the ‘faster wind’ model which gives rise to the largest change in the Tully–Fisher relation by increasing the rotational velocity at a fixed stellar mass. This result is a consequence of the efficiency with which the fast wind model suppresses SF.

Figure 10.

View largeDownload slide

Tully–Fisher relation at z = 0 (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data are taken from Bell & de Jong (2001) and Reyes et al. (2012). The solid lines mark the simulation median relations, whereas the two-dimensional histograms indicate the distribution of the L25n512 result (left-hand panel) and L25n256 result (right-hand panel). The runs converge towards a relation which is very similar to the marked observed relation with the slope and normalization of the simulated Tully–Fisher relation falling almost directly on top of the two marked observational Tully–Fisher relations. Most variations of the feedback model parameter choices lead to only minor changes in the resulting Tully–Fisher relation. The exception is the ‘faster wind’ model which gives rise to the largest change in the Tully–Fisher relation by increasing the rotational velocity at a fixed stellar mass. This result is a consequence of the efficiency with which the fast wind model suppresses SF.

Examining the different resolutions shown in the left-hand panel of Fig. 10, we find that there is good convergence for high-mass systems, consistent with our previous findings. Furthermore, our model converges towards a relation which is very similar to the observed relation, with the slope and normalization of the simulated Tully–Fisher relation falling almost directly on top of the two observational Tully–Fisher determinations. We note that we did not select the galaxies in Fig. 10 according to their type. This means that especially at the more massive end our result will be influenced by a certain fraction of elliptical galaxies, implying that we should not expect perfect agreement with the observations in this regime.

In the right-hand panel of Fig. 10, we find that most variations in the feedback parameters do not induce significant changes in the resulting Tully–Fisher relation. The exception is the ‘faster wind’ model which gives rise to the largest change in the Tully–Fisher relation by increasing the rotational velocity at a fixed stellar mass. This result is a consequence of the efficiency with which the fast wind model suppresses SF, as seen previously in Fig. 7. Aside from the ‘faster wind’ model, we find that our feedback models are able to produce reasonable normalizations and slopes for the Tully–Fisher relation independent of the detailed feedback parameter choices. We therefore conclude that the Tully–Fisher relation is not very sensitive to the details of our physics parametrization. In fact, among all the observables considered in this paper, the Tully–Fisher relation is the most stable against variations of the underlying physics. However, the ‘no feedback’ simulation clearly fails to reproduce this key observable.

4.2.7 Mass–metallicity relation

Our stellar evolution and enrichment implementation allows us to study chemical abundances of galaxies in detail and confront those with observations. In Fig. 11, we show the z = 0 mass–metallicity relation and compare it to SDSS DR7 data (Abazajian et al. 2009). Two distinct mass–metallicity relations are shown. Although both use the SDSS data, one is the KK04 (Kobulnicky & Kewley 2004) nebular emission line diagnostic relation, while the other is the PP04 (Pettini & Pagel 2004) diagnostic. To facilitate this comparison, we have taken the metallicity values as tabulated in Zahid et al. (2012) using the KK04 diagnostic and converted them to the PP04 diagnostic using the empirical relations of Kewley & Ellison (2008). The main reason for using this procedure is to include the relevant information about the known uncertainty in the normalization of nebular emission line metallicity indicators, and to avoid comparing our models to only one diagnostic method. We calculate the metallicity value for each simulated galaxy by finding the average oxygen-to-hydrogen abundance, weighted by the SFRs of gas and excluding non-star-forming gas from this calculation in order to provide a more straightforward comparison to observations, where the nebular emission lines naturally probe star-forming regions.

Figure 11.

View largeDownload slide

Mass–metallicity relation as a function of stellar mass (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data are taken from SDSS DR7 (Abazajian et al. 2009). The solid lines mark the simulation median relations, whereas the two-dimensional coloured histograms indicate the distribution of the L25n512 (left-hand panel) and L25n256 (right-hand panel) results. The metallicity values are well converged, but the relationship which they converge to does not perfectly match the observations. Specifically, the observed mass–metallicity slope is shallower than the simulated slope, and there is no apparent turnover in the simulated mass–metallicity relation at high masses. However, the simulated mass–metallicity relation falls between the PP04 and KK04 relations indicating that galaxies are retaining a reasonable fraction of their metals. The right-hand panel shows that the simulated mass–metallicity relation is quite sensitive to the adopted physics parameters since feedback variations strongly affect the amount of metals retained in galaxies and expelled to haloes.

Figure 11.

View largeDownload slide

Mass–metallicity relation as a function of stellar mass (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data are taken from SDSS DR7 (Abazajian et al. 2009). The solid lines mark the simulation median relations, whereas the two-dimensional coloured histograms indicate the distribution of the L25n512 (left-hand panel) and L25n256 (right-hand panel) results. The metallicity values are well converged, but the relationship which they converge to does not perfectly match the observations. Specifically, the observed mass–metallicity slope is shallower than the simulated slope, and there is no apparent turnover in the simulated mass–metallicity relation at high masses. However, the simulated mass–metallicity relation falls between the PP04 and KK04 relations indicating that galaxies are retaining a reasonable fraction of their metals. The right-hand panel shows that the simulated mass–metallicity relation is quite sensitive to the adopted physics parameters since feedback variations strongly affect the amount of metals retained in galaxies and expelled to haloes.

In the left-hand panel of Fig. 11, we find that the metallicity values are well converged for our runs with varying resolution. In fact, the convergence is slightly better than that in most of the other relations which we examine in this paper. However, the relationship which they converge to does not perfectly reproduce the observations. There are two main features that make the simulation results distinct from the observations. The first is that the observed mass–metallicity slope is shallower than the simulated slope – regardless of whether we adopt the PP04 or KK04 diagnostic. The second issue is that there is no apparent turnover in the simulated mass–metallicity relation at high masses, unlike observed. Although there is a flattening visible for the two highest resolutions (L25n256 and L25n512), this occurs at too large stellar masses and at too high metallicities to account for this observational signature. However, despite these shortcomings, the simulated mass–metallicity relation falls between the PP04 and KK04 relations for most systems, which indicates that the simulated galaxies retain a reasonable fraction of their metals, in acceptable agreement with observations, on average (similar level agreement was found in Davé, Finlator & Oppenheimer 2011b).

We note that this agreement has become possible only through the independent treatment of mass and metal loading of stellar winds. As described above, we have adopted a wind metal loading factor (γ_w), which determines the metallicity of the wind material relative to the local ISM metallicity from where the wind is launched. Without such a scheme our galaxies would retain too few metals due to the large wind mass loadings required to match the observed stellar masses of lower mass systems. This creates an interesting tension between the need for high wind efficiencies to reduce the buildup of stellar mass in low-mass galaxies and the need for low-mass galaxies to be relatively efficient at retaining their metal content (Zahid et al. 2012). The solution to this tension likely either lies in (i) fundamentally reducing the accretion efficiency of low-mass galaxies via some feedback mechanism that is able to disrupt and heat the IGM around these systems or (ii) in a more detailed wind model that self-consistently handles the entrainment of low-metallicity material as winds propagate out of galaxies (e.g. Hopkins et al. 2012b).

An interpretation of the differences between the simulated mass–metallicity relation and the observed relation is guided by the right-hand panel of Fig. 11. We find that the simulated mass–metallicity relation is quite sensitive to the adopted feedback physics parameters. The reason for this behaviour is that the level of central galactic metallicity is easily reduced if we encourage mixing between the central galaxy gas and the halo gas, which is mediated by AGN or stellar feedback. Our winds, which by construction transport highly enriched disc gas into the halo, are very efficient at mixing metals. Thus, when we use a ‘weaker wind’ prescription, we find that the average metal content of the central dense galaxy gas goes up. On the other hand, if we employ ‘stronger winds’, we find that the average metal content goes down. The need to regulate the growth of low-mass galaxies successfully while at the same time avoiding overly depleting galaxies of their metal content prompted us to implement a specific wind metal loading, as described above. We note that variations in the AGN radio mode alone are not able to remove the metallicity tension at the high-mass end. Most of the simulations overshoot the amount of metals in massive galaxies by 0.3-0.5 dex.

4.2.8 Luminosity function

Figure 12.

View largeDownload slide

SDSS-band (g, r, i, z) luminosity functions for our simulated galaxy population (resolution study). We compare to double Schechter fits of the observational data, taken from Blanton et al. (2005). We find reasonable agreement of L25n256 and L25n512 with these fits. The faint-end slope agrees very well with the slope inferred from observational data, although the highest resolution simulations have a slightly too high normalization. Also, the exponential suppression is reasonably well reproduced but starts at slightly too bright galaxies. This discrepancy is largest for the g-band luminosity function.

Figure 12.

View largeDownload slide

SDSS-band (g, r, i, z) luminosity functions for our simulated galaxy population (resolution study). We compare to double Schechter fits of the observational data, taken from Blanton et al. (2005). We find reasonable agreement of L25n256 and L25n512 with these fits. The faint-end slope agrees very well with the slope inferred from observational data, although the highest resolution simulations have a slightly too high normalization. Also, the exponential suppression is reasonably well reproduced but starts at slightly too bright galaxies. This discrepancy is largest for the g-band luminosity function.

We compare the luminosity functions of the different feedback models to the observational fits in Fig. 13. Our fiducial L25n256 model gives a good fit to the observed luminosity function, whereas a model with ‘faster winds’ clearly gives the worst agreement. Weakening the radio-mode AGN feedback leads to a significant overshoot of the luminosity function at the bright end. ‘Weaker winds’ cause an overproduction of faint galaxies. Finally, the ‘no feedback’ clearly fails dramatically in reproducing the observed luminosity function, as it did for the stellar mass function. Interestingly, the ‘higher radio threshold’ simulation produces luminosity functions which are in better agreement with the observations than our fiducial model. This is true for all bands, where the exponential drop at the bright end is systematically in better agreement with the observational data for the ‘higher radio threshold’ simulation. Also the model with ‘stronger radio-mode’ feedback leads to better agreement of the luminosity functions in all bands for the bright end. However, although these two models agree better with the observed luminosity functions, their stellar mass functions do not agree as well with the observations as our model with the fiducial feedback settings.

Figure 13.

View largeDownload slide

SDSS-band luminosity (g, r, i, z) functions for our simulated galaxy population (different feedback models). We compare to observational double Schechter fits taken from Blanton et al. (2005). Our fiducial model reproduces the observational fit reasonably well. Models with faster winds clearly give the poorest agreement. Weakening the radio-mode AGN feedback leads to a significant overshoot of the luminosity function at the bright end. ‘Weaker winds’ lead to an overproduction of faint galaxies. The ‘no feedback’ simulation clearly fails dramatically in reproducing the observed luminosity function. Interestingly, a model with a higher radio-mode threshold leads to better agreement for the bright end in all bands.

Figure 13.

View largeDownload slide

SDSS-band luminosity (g, r, i, z) functions for our simulated galaxy population (different feedback models). We compare to observational double Schechter fits taken from Blanton et al. (2005). Our fiducial model reproduces the observational fit reasonably well. Models with faster winds clearly give the poorest agreement. Weakening the radio-mode AGN feedback leads to a significant overshoot of the luminosity function at the bright end. ‘Weaker winds’ lead to an overproduction of faint galaxies. The ‘no feedback’ simulation clearly fails dramatically in reproducing the observed luminosity function. Interestingly, a model with a higher radio-mode threshold leads to better agreement for the bright end in all bands.

4.2.9 BH mass–stellar mass relation

The BH mass growth is mainly regulated by quasar-mode feedback which together with the BH accretion yields a tightly self-regulated feedback loop. As a result, a BH mass–stellar mass relation is expected to be produced through these two processes (see also Di Matteo et al. 2005; Springel et al. 2005a). Here we compare the BH mass–stellar mass relation to observational data of Häring & Rix (2004). As the left-hand panel of Fig. 14 demonstrates, the implemented quasar-mode model leads to the correct growth of BHs in our simulation. Most importantly, we reproduce the correct slope of the relation. This demonstrates that the modifications of the original Springel et al. (2005a) model described above do not affect the BH mass–stellar mass relation in any significant way. Reproducing this relation is also crucial to inject the correct amount of radio-mode AGN feedback generated by these BHs and to create the appropriate radiation field used for the radiative (EM) AGN feedback. We note that we do not plot the stellar bulge mass in Fig. 14 but rather the stellar mass within twice the stellar half-mass radius, as discussed above. Correcting for the bulge mass does however not affect the conclusions of this subsection. We will in forthcoming work study this relation in more detail using extracted bulge masses (Sijacki et al., in preparation).

Figure 14.

View largeDownload slide

BH mass–stellar mass relation (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data are taken from Häring & Rix (2004). The BH mass–stellar mass relation is well converged and agrees with the observations both in slope and in normalization. Also the amount of scatter is in reasonable agreement with the observational data. Normalization changes in this relation are mainly due to changes in the stellar mass caused by the different feedback choices. For example, ‘fast winds’ strongly reduce the number of stars forming in haloes, as demonstrated above. We note that the right-hand panel does not include the ‘no feedback’ simulation since this run did not include any BHs.

Figure 14.

View largeDownload slide

BH mass–stellar mass relation (left-hand panel: resolution study; right-hand panel: different feedback models). Observational data are taken from Häring & Rix (2004). The BH mass–stellar mass relation is well converged and agrees with the observations both in slope and in normalization. Also the amount of scatter is in reasonable agreement with the observational data. Normalization changes in this relation are mainly due to changes in the stellar mass caused by the different feedback choices. For example, ‘fast winds’ strongly reduce the number of stars forming in haloes, as demonstrated above. We note that the right-hand panel does not include the ‘no feedback’ simulation since this run did not include any BHs.

The left-hand panel of Fig. 14 also shows that our BH mass–stellar mass relation is well converged and agrees with the observations both in slope and in normalization. Also, the amount of scatter is in reasonable agreement with the observational data. The right-hand panel of Fig. 14 demonstrates that the slope of this relation does not change significantly if we alter the stellar or radio-mode AGN feedback. We note that this plot does not include the ‘no feedback’ simulation since this run did not include BHs due to the absence of any BH growth regulating feedback mechanism. Normalization offsets in this relation are mainly due to the changes in the stellar masses caused by the different feedback choices. For example, ‘faster winds’ strongly reduce the amount of stars forming in haloes. This implies that the low-mass stellar systems end up with too massive BHs at their centres, as can be seen in the right-hand panel. Overall we find that our fiducial model reproduces the BH mass–stellar mass relation reasonably well.

4.3 Disentangling the impact of different physical processes

So far we have explored feedback variations around our fiducial model, varying the strength of stellar and AGN feedback. In this subsection, we briefly discuss the impact of the different physical effects on the cosmic SFR density, the stellar mass function, the stellar mass–halo mass relation and the stellar mass density. Most of the simulations presented here are unrealistic in the sense that they deliberately ignore important processes required to shape the galaxy population. But they are useful for demonstrating the impact of these physical processes and how they influence the galaxy population. In particular, we discuss the importance of our newly implemented radiative AGN feedback scheme and demonstrate the impact of our wind metal loading scheme.

It is instructive to consider a simulation set with increasing model complexity, starting from a physics setup close to the simulations presented in Vogelsberger et al. (2012) and ending at the full physics implementation discussed above. The different simulations of this series are summarized in Table 3. The ‘plain’ simulation differs only in terms of the IMF (Chabrier instead of Salpeter), a softer eEOS (q = 0.3 instead of q = 1.0) and self-shielding corrections from the simulations presented in Vogelsberger et al. (2012). We then add more physical processes on top of this plain setup, where the parameters of the individual physics modules are chosen according to our fiducial setup discussed above (see Table 1). We start with the simulation MeGa, which includes on top of the ‘plain’ setup metal-line cooling and gas recycling. This simulation is identical to the ‘no feedback’ simulation presented in the previous section. The last simulation listed (MeGaWiMlQuEmRa) includes all physical effects discussed above and is identical to the fiducial L25n256 simulation. Table 3 also includes a simulation which is similar to L25n256, i.e. our full fiducial model, but without stellar feedback (NoWindFeed). This simulation is useful for studying the detailed impact of stellar feedback, but still considering all the other physical processes included in our model. We note that MeGaWiMl explores the situation with full stellar feedback, but without AGN feedback.

The upper-left panel of Fig. 15 shows the cosmic SFR density for the different simulations, whereas the upper-right panel presents the stellar mass function. We do not include observational data in these plots for graphical clarity and because a detailed comparison to the observational data was already presented above. Going from the ‘plain’ model to the model which includes in addition metal-line cooling and gas recycling (MeGa = ‘no feedback’) strongly increases the SFRs since net cooling rates are increased and recycled gas can form additional stars. However, such a model still correctly produces an SFR peak although at a slightly different redshift compared to the ‘plain’ and full MeGaWiMlQuEmRa setup. Including stellar winds (MeGaWi) regulates this behaviour at early times by strongly suppressing the SF in low-mass systems. This implies that a significant amount of gas is not turned into stars at high redshifts in this simulation. This gas is then available for SF at late times, which results in a rising SFR density towards lower redshift. Metal-line cooling and gas recycling are so efficient then that the SFR density almost does not decline. In fact, a substantial decline can be achieved to a sufficient degree only by including radio-mode AGN feedback as demonstrated by the MeGaWiMlQuEmRa curve.

Figure 15.

View largeDownload slide

Cosmic SFR density (upper-left panel), stellar mass function (upper-right panel), stellar mass–halo mass relation (lower-left panel) and stellar mass density (lower-right panel) of simulations with increasing modelling complexity (see Table 3). The ‘plain’ setup does not include any explicit stellar and AGN feedback. The MeGaWiMlQuEmRa simulation includes all physical processes discussed above with parameters set as for our fiducial model (see Table 1). Although most of these simulations are unrealistic since they miss important physical ingredients, they are helpful in disentangling the impact of various physical processes on the galaxy population. Simulations including stellar winds, but no radio-mode AGN feedback, lead to a significant overproduction of stars at late times. This can only be regulated by strong radio-mode AGN feedback. Quasar-mode feedback alone has essentially no effect on the SF (see MeGaWiMl versus MeGaWiMlQu). However, radiative (EM) AGN feedback leads to a minor reduction of the SFR. Changing the metal loading of stellar winds (MeGaWi versus MeGaWiMl) has about the same effect on the SFR. As expected, stellar feedback regulates SF in low-mass systems, whereas the stellar content of higher mass haloes is mainly set through radio-mode AGN feedback.

Figure 15.

View largeDownload slide

Cosmic SFR density (upper-left panel), stellar mass function (upper-right panel), stellar mass–halo mass relation (lower-left panel) and stellar mass density (lower-right panel) of simulations with increasing modelling complexity (see Table 3). The ‘plain’ setup does not include any explicit stellar and AGN feedback. The MeGaWiMlQuEmRa simulation includes all physical processes discussed above with parameters set as for our fiducial model (see Table 1). Although most of these simulations are unrealistic since they miss important physical ingredients, they are helpful in disentangling the impact of various physical processes on the galaxy population. Simulations including stellar winds, but no radio-mode AGN feedback, lead to a significant overproduction of stars at late times. This can only be regulated by strong radio-mode AGN feedback. Quasar-mode feedback alone has essentially no effect on the SF (see MeGaWiMl versus MeGaWiMlQu). However, radiative (EM) AGN feedback leads to a minor reduction of the SFR. Changing the metal loading of stellar winds (MeGaWi versus MeGaWiMl) has about the same effect on the SFR. As expected, stellar feedback regulates SF in low-mass systems, whereas the stellar content of higher mass haloes is mainly set through radio-mode AGN feedback.

The MeGaWiMl simulation employs our new metal loading scheme for outflows. Reducing the metal loading of winds from γ_w = 1 (MeGaWi) to our fiducial γ_w = 0.4 (MeGaWiMl) decreases the SFR slightly towards lower redshifts. This is because less enriched winds lead to less enriched gaseous haloes. As a result, these haloes then experience less metal-line cooling yielding a lower gas supply to central galaxies such that the SFR decreases. Although this effect is clearly visible, we note however that it does not represent a substantial change in the SFRs, which is clearly largely shaped by feedback processes. As we discussed above, such a metal loading scheme is required to retain enough metals in galaxies and to obtain a reasonable match to the galaxy mass–metallicity relation. In fact, γ_w = 1 would lead to a far too small normalization of the mass–metallicity relation.

Comparing the MeGaWiMl simulation with a run that additionally includes quasar-mode AGN feedback (MeGaWiMlQu) demonstrates that continuous thermal AGN feedback has no significant effect on the global SFRs. The MeGaWiMlQu simulation differs from the MeGaWiMlQuEm only by radiative EM AGN feedback, which is included in the latter. Most importantly, this causes quasar-mode feedback to have a slightly more noticeable effect on the SFRs, i.e. our new radiative feedback is stronger than the continuous thermal AGN feedback of our quasar-mode feedback implementation alone. Radiative AGN feedback leads to a suppression of atomic cooling and to an increase of the heating rates, so that SF is further suppressed by this feedback channel. This can be seen in the cosmic SFR density, where we find a lower SFR once quasars become active. Although the effect is visible in the SFRs, it is not particularly strong compared to radio-mode AGN feedback, which is implemented in the MeGaWiMlQuEmRa simulation. The reason is that EM AGN feedback is effective when the accretion rate is near Eddington in the quasar regime, which typically happens only for short periods of cosmic time for any given BH. This questions previous claims suggesting that such radiative feedback is a highly important and very efficient feedback channel compared to thermal and mechanical AGN feedback (e.g. Gnedin & Hollon 2012). We stress however that we focus here only on the effect of EM feedback on the amount of stellar mass formed, and not do analyse the impact of the AGN radiation field on the surrounding IGM, for example.

Finally, the MeGaWiMlQuEmRa simulation has an equivalent setup as the L25n256 simulation presented above. It includes all its processes and uses the same parametrization as our fiducial model (see Table 1). Clearly, radio-mode feedback leads to a strong suppression of SF at late times in massive haloes. This induces a steep decline of the cosmic SFR density and an exponential drop-off of the stellar mass function towards massive haloes. The NoWindFeed run, which does not include winds, but is otherwise identical to L25n256, does not show any SF suppression at early times such that the overall normalization of the SFR density is significantly too high. We conclude that late-time SF is largely shaped by metal-line cooling, gas recycling and radio-mode AGN feedback. The high-redshift SFR density, on the other hand, is essentially only determined through stellar feedback. Metal-line cooling and gas recycling have essentially no effect in that regime. This is expected since not many stars formed at higher redshifts such that no significant amount of metals is available and also mass return of stars is not important at these early times.

The upper-right panel of Fig. 15 demonstrates the effect of our physics implementation on the stellar mass function. Clearly, models without any feedback (‘plain’ and MeGa) lead to a substantial overproduction of stars. This is also true for the simulation without any explicit stellar feedback (NoWindFeed). The incorporation of winds is the main driver for the low-mass slope of the stellar mass function. The separate metal loading of these winds has essentially no impact on the stellar mass function (MeGaWi versus MeGaWiMl). The bright end is largely shaped by radio-mode AGN feedback, which starts to dominate around |$\mathrm{M_\star } \sim 10^{11} {\,\rm M_{\odot }}$|⁠. Also radiative (EM) AGN feedback has only a rather minor impact on the massive end of the stellar mass function.

The lower-left panel of Fig. 15 shows the stellar mass–halo mass relation. For clarity, we only show the median relations here. As expected, simulations without stellar feedback (‘plain’, MeGa and NoWindFeed) lead to a large overproduction of stellar mass for lower- and intermediate-mass systems. The addition of stellar winds leads to about an order of magnitude suppression of stellar mass. The inclusion of quasar-mode and radiative (EM) AGN feedback does not significantly affect the massive end. Only the simulation with radio-mode AGN feedback (MeGaWiMlQuEmRa) leads to a reduction of stellar mass for higher mass systems, resulting in the correct turnover behaviour of the stellar mass–halo mass relation.

The lower-right panel of Fig. 15 shows the stellar mass density as a function of redshift, i.e. the integrated version of the upper-left panel. Here we find a behaviour similar to what we found for the other three diagnostics discussed so far. The lack of stellar feedback in the ‘plain’, MeGa and NoWindFeed simulations strongly overproduces the stellar mass density already at high redshifts. Radio-mode AGN feedback on the other hand strongly reduces the number of stars formed at later times (MeGaWiMlQuEmRa). All other curves follow in a narrow range between these extremes.

We conclude that the stellar content of haloes is largely dominated by stellar feedback and radio-mode AGN feedback. Other model details, like the independent metal loading of winds or radiative (EM) AGN feedback, lead only to very minor changes in the stellar mass, although they can have important consequences for other observables, such as the galaxy mass–metallicity relation in the case of wind metal loading.

5 SUMMARY AND CONCLUSIONS

Matching the enormous amount of low-redshift observational data is a great challenge for any galaxy formation model that tries to explain the galaxy population from first principles. Especially in the era of ‘precision cosmology’ based on large surveys like SDSS, LSST, etc., it is crucial to have predictive, accurate and reliable galaxy formation models to test our understanding of structure formation by confronting these models with a plethora of data. The most general self-consistent way to study structure formation is through hydrodynamical simulations which follow the coupled evolution of DM, DE and baryonic physics. This task requires two main ingredients: (i) a highly reliable and computationally efficient scheme to solve the basic hydrodynamical and gravity equations, and (ii) a well-motivated and numerically meaningful model for the required physics of galaxy formation.

We have recently demonstrated that a moving-mesh approach for solving the underlying gas dynamics is very attractive for cosmic structure formation and offers advantages compared to other numerical schemes for galaxy formation simulations (Springel 2010; Bauer & Springel 2012; Kereš et al. 2012; Sijacki et al. 2012; Torrey et al. 2012; Vogelsberger et al. 2012; Nelson et al. 2013). This approach hence addresses requirement (i). The natural next step is to equip this new numerical scheme with the relevant physical processes to fulfil requirement (ii), with the goal of carrying out studies of galaxy formation in large volumes with higher accuracy and fidelity to the physics than possible thus far.

Along these lines, we have presented in this paper the first galaxy formation implementation for the moving-mesh code arepo specifically aimed towards large-scale cosmological hydrodynamical simulations. Our objective has been to construct a model which results in a realistic galaxy population, similar to what the best semi-analytic models achieve, yet realized in a self-consistent hydrodynamical simulation of galaxy formation. This is an important next step in galaxy formation theory and promises a much higher predictive power from the theoretical models.

Our SF prescription (see Section 2.2) models star-forming dense ISM gas with an eEOS, where stars form stochastically above a certain density threshold. To reach these densities, gas cools through primordial and metal-line cooling in the presence of a photoionizing UV background with self-shielding corrections (see Section 2.4). Primordial cooling is treated with a self-consistent chemical network, whereas metal-line cooling rates are tabulated based on cloudy calculations. Both processes include self-shielding corrections from the otherwise spatially uniform UV background radiation. We follow stellar evolution and chemical enrichment processes based on AGB, SNIa and SNII yields and account for stellar mass-loss processes (see Section 2.3). For the chemical enrichment, we follow nine elements (H, He, C, N, O, Ne, Mg, Si, Fe). Stellar feedback effects are based on a kinetic outflow prescription, where stellar winds are either non-locally driven from the ISM gas directly or locally around evolving stellar particles (see Section 2.5). Our implementation supports energy- and momentum-driven winds or a combination of both, although in practice we prefer energy-driven winds since they provide a better fit to the faint end of the stellar mass function. We introduce a novel wind metal mass loading scheme, which allows us to decouple mass and metal loading for gas put into the outflowing wind (see Section 2.5.3), which is required to match the mass–metallicity relation and stellar mass function for low-mass galaxies simultaneously. We also follow BH growth and include quasar- and radio-mode feedback processes (see Section 2.6). Furthermore, we introduce a novel model for radiative (EM) AGN feedback (see Section 2.6.4), which is implemented assuming an average AGN SED and a bolometric, luminosity-dependent scaling to take into account obscuration effects. Finally, we introduced inlined analysis techniques for gas tracking in the form of classical Lagrangian velocity field tracer particles and a new Monte Carlo-based scheme, and we added on-the-fly visualization routines to allow for direct volume rendering through a ray-casting technique while the simulation is running.

The scope of the paper is to present our galaxy formation implementation, to identify a fiducial model with promising settings of the feedback parameters and to show some basic results for the newly added physics that are compared with some key observations at z = 0. To this end, we have explored three different sets of simulations: one resolution study using our fiducial feedback parameters with a maximum resolution of 2 × 512³ resolution elements in a (25 h⁻¹ Mpc)³ simulation volume, a set of simulations where we slightly modified the strength of stellar and radio-mode AGN feedback to explore changes in the feedback prescriptions, and simulations where we systematically explored the impact of various physical processes by specifically deactivating certain parts of our model.

Our fiducial feedback parameters were chosen to be physically plausible and set such that they reasonably well reproduce key observations at z = 0. We note that the stellar and radio-mode AGN feedback settings of this fiducial model differ significantly from those previously used in SPH simulations. Most important, the AGN radio-mode feedback is more energetic than adopted earlier. These energetic feedback processes are required to match the z = 0 stellar mass function and stellar mass–halo mass relation. Such a need for stronger feedback was already pointed out in Vogelsberger et al. (2012) and Kereš et al. (2012) where it was demonstrated that arepo shows significantly more efficient gas cooling compared to previous SPH simulations, which is mainly due to spurious viscous heating effects in SPH and the lack of a proper cascade of subsonic turbulent energy to small scales. This in turn results in an overall increase of SF by a factor of 2 towards lower redshifts (Vogelsberger et al. 2012) and also a stellar mass function which significantly overpredicts the number of massive systems (Kereš et al. 2012). Stronger feedback needs to compensate for this increase in cooling and SF. The artificial heating effect of SPH is largest for high-mass systems where a quasi-hydrostatic hot atmosphere can form. SPH simulations suffer here a form of ‘numerical quenching’, reducing the amount of AGN feedback required in this technique to regulate the amount of stars forming in these systems. Since arepo does not suffer from this unphysical effect, we need to invoke a corresponding increased radio-mode feedback strength.

The feedback variations explored in this paper are summarized in Table 2. We only varied the strength of stellar and radio-mode feedback since those influence SF most significantly. These modifications also affect the times at which stars form in the simulation (see Fig. 16), and therefore many of the key observables studied above are significantly altered by these feedback changes. The most extreme feedback models, in terms of stellar formation times, are our simulations with ‘stronger winds’, which produces too much late-time SF, and the simulation with ‘faster winds’, which leads to almost no late-time SF. All other feedback variations fall more or less between these two extreme setups. ‘Weaker winds’ lead, as expected, to more SF at higher redshifts. Weakening the radio-mode AGN feedback results in more late-time SF, whereas a stronger radio-mode feedback has the opposite effect. Interestingly, a simulation without any feedback leads to a lower late-time SFR than the model with very strong winds. The reason is that strong winds very efficiently suppress SF at early times such that significant quantities of gas remain for late-time SF.

Figure 16.

View largeDownload slide

Distribution of SF times for all stars at z = 0 (left-hand panel: resolution study; right-hand panel: different feedback models). Stellar formation times are well converged. The feedback variations strongly affect when stars form in the simulation. The most extreme feedback models are the simulation with ‘strong winds’, which leads to a great deal of late-time SF, and the simulation with ‘faster winds’, which yields almost no late-time SF. All other feedback variations fall more or less between these two extreme setups.

Figure 16.

View largeDownload slide

Distribution of SF times for all stars at z = 0 (left-hand panel: resolution study; right-hand panel: different feedback models). Stellar formation times are well converged. The feedback variations strongly affect when stars form in the simulation. The most extreme feedback models are the simulation with ‘strong winds’, which leads to a great deal of late-time SF, and the simulation with ‘faster winds’, which yields almost no late-time SF. All other feedback variations fall more or less between these two extreme setups.

Before discussing how our simulations compare to observations and how changes in the stellar and AGN feedback strength affect this comparison, we will first briefly elaborate on the importance of the two most novel aspects of our galaxy formation model, namely the incorporation of radiative AGN feedback and the special treatment of metals for stellar outflows.

5.1 Impact of radiative AGN feedback

Our galaxy formation model includes a new phenomenological model for radiative (EM) AGN feedback. This feedback acts on halo gas by changing the net cooling rates due to the radiation coming from AGN, which is particularly strong during quasar activity. Our model assumes a fixed SED with a simple and rather conservative prescription for obscuration effects. We explored for the first time the impact of this radiative AGN feedback in fully self-consistent large-scale cosmological simulations with a focus on its consequences for the host galaxy in terms of the stellar mass content and overall SFRs. We find that this form of feedback is more efficient in regulating SF than purely thermal, continuous quasar-mode AGN feedback, which is mainly relevant to regulate BH growth. However, this radiative form of feedback is only efficient if the BH accretion rates are close to Eddington where the radiation field is strongest. Since this happens only for a rather short period of cosmic time, the impact of mechanical radio-mode feedback is typically significantly stronger than radiative feedback. However, we did not explore possible effects on the ionization state of the IGM, which could be more substantially influenced by radiative AGN feedback.

5.2 Metal loading of stellar winds

Our stellar feedback model allows the metal loading of winds to be adjusted separately from the actual mass loading. We found that such a treatment is necessary to reproduce the stellar mass and metal content of galaxies simultaneously. Specifically, we have introduced a wind metal loading factor which is set independently from the wind mass loading factor. This wind metal loading factor defines the relationship between the metallicity of newly created wind particles and the metallicity of the ambient ISM. We found that the traditional wind metal loading leads to an underproduction of the gas-phase metallicity and an overenrichment of halo gas. Both of these directly result from the efficient ejection of metals from galaxies via winds. Not only does this lead to poor matches to the observationally constrained mass–metallicity relation, but the overenrichment of halo gas leads to enhanced metal-line cooling as well. We can remove this tension with our metal loading scheme for stellar winds. In fact, such a scheme allows us to simultaneously match the mass–metallicity relation and stellar mass function reasonably well. While this wind treatment significantly affects the mass–metallicity relation, it only slightly modifies other observables like the stellar mass function. We note that metal loading factors can be constrained empirically (Zahid et al. 2013). Such studies point to rather low values for γ_w.

5.3 Comparison to observations

Our fiducial model reproduces many key observables over a wide range of halo masses, covering low-mass and high-mass systems through the combination of strong stellar and AGN feedback. Most importantly, our model correctly describes the transition region around |$M_\mathrm{200,crit}\sim 10^{12}{\,\rm M_{\odot }}$|⁠, where SF is maximally efficient. We briefly summarize our main findings for each of the key characteristics presented in this paper.

• Cosmic SFR density. Our fiducial model reproduces the cosmic SFR density as a function of redshift reasonably well, although our highest resolution simulation is not yet fully converged. The wind velocity of the stellar feedback has a significant impact on the cosmic SFR density. Fast winds lead to a strong suppression of the SFR, and such fast winds can reduce the SF even in more massive haloes. The impact of variations of the AGN radio mode is less dramatic and mainly affects the late-time decline of the cosmic SFR density, but no radio mode altogether gives a far too high late-time SFR.

• Stellar mass–halo mass relation. Our model predicts a stellar mass–halo mass relation that is in good agreement with results based on abundance matching. At higher galaxy masses (i.e. M > 10⁹ M_⊙), we find clear convergence in the stellar mass–halo mass relation. Most importantly, our combination of stellar feedback and radio-mode AGN feedback leads to a shape in the stellar mass–halo mass relation which is in good agreement with abundance matching results, including a correct turnover at the massive end due to efficient quenching through AGN radio-mode feedback in massive haloes. We need strong stellar and AGN feedback to reproduce the stellar content of low- and high-mass haloes.

• Stellar mass function. Our model reproduces the stellar mass function well, although we find that we have slightly too many low-mass galaxies in our highest resolution simulation. We find a similar problem for lower mass systems for the stellar mass–halo mass relation. Here we also slightly overproduce the amount of stars in these systems. Both findings indicate the need for yet stronger stellar feedback, which we partially explored in this work by increasing the SNII energy factor by a factor of 2. But we note that our fiducial simulation setup includes already rather strong SNII feedback, and increasing this further appears energetically problematic. Radio-mode AGN feedback successfully produces the sharp exponential drop-off of the stellar mass function towards larger stellar masses. We also find that energy-driven winds are required to suppress the low-mass end of the stellar mass function sufficiently since they allow for large mass loading factors for low-mass systems compared to momentum-driven winds. This is essentially required to obtain a shallow slope of the stellar mass function at the low-mass end.

• Stellar mass density. Our fiducial model reproduces the observed stellar mass density reasonably well. The ‘no feedback’ simulation strongly overproduces the amount of stellar mass, whereas the ‘faster winds’ strongly underproduce it. All the other feedback variations are quite close to the fiducial model, which provides the best fit to the data. Changes in the radio-mode AGN feedback only affect the late-time evolution of the stellar mass density.

• Tully–Fisher relation. We find good convergence for the Tully–Fisher relation, and our highest resolution fiducial model agrees very well with observations both in terms of normalization and slope. This is an important result of our model since it implies that our galaxies, in general, have the correct dynamical structure. Variations of the feedback model parameter choices lead to only negligible changes in the resulting Tully–Fisher relation, implying that the Tully–Fisher relation is a rather stable outcome of our model. The exception is the ‘fast wind’ model which gives rise to the largest change in the Tully–Fisher relation by increasing the rotational velocity at a fixed stellar mass. Also a model without any feedback does not reproduce the Tully–Fisher relation.

• Mass–metallicity relation. Among all the relations that we examined in this paper, the mass–metallicity relation shows the largest tension between observations and our simulations. Although we find a positive correlation between the stellar mass and [O/H], the overall amplitude of this relation is not fully correctly reproduced. But given the observational uncertainties, we do not consider this a very severe problem at this point. The main issue is that the simulations overpredict the metallicity of higher mass systems. Here our simulations do not predict the turnover in the mass–metallicity relation correctly, and therefore lead to too many metals in massive galaxies. For lower- to intermediate-mass systems, we can only achieve reasonable agreement with the observations through our novel wind metal loading scheme, which decouples the actual mass loading from the metal loading of winds. Without such a scheme, low-mass galaxies particularly would lose far too many metals owing to the high mass loadings of our energy-driven wind scalings. However, such high mass loadings are required to reduce SF sufficiently in these low-mass systems.

• Luminosity function. We used stellar population synthesis models to derive SDSS-band luminosity functions (g, r, i, z bands). Our highest resolution simulations produce luminosity functions at z = 0 which are in reasonable agreement with observations. Especially, the faint-end slope and the exponential drop towards brighter systems are well reproduced. Our highest resolution simulation has a slightly too high faint-end normalization compared to the observations. Also the exponential drop happens at slightly too high halo masses. Interestingly, this effect is strongest for the g-band luminosity function, whereas the agreement with the other bands is typically better. Feedback models with faster winds strongly disagree with the observations. Weakening the radio-mode AGN feedback leads to a significant overshoot of the luminosity function at the bright end while weak winds cause an overproduction of faint galaxies.

• BH mass–stellar mass relation. BH growth is mainly regulated by quasar-mode feedback in our model, which results together with BH accretion in a tightly self-regulated feedback loop. Our simulated BH mass–stellar mass relation is well converged and agrees with the observations both in slope and in normalization. This demonstrates that our modifications to the BH accretion procedure do not affect BH growth in any significant way. Furthermore, the slope of the BH mass–stellar mass relation does not change strongly if we alter the stellar or radio-mode AGN feedback parameters. Normalization changes in the relation are mainly due to changes in the stellar mass caused by the different feedback choices.

Overall, we find that our fiducial model reproduces a significant number of key observations of the low-redshift Universe. As demonstrated in Torrey et al. (2013), our model also reproduces many high-redshift observables. Furthermore, Marinacci et al. (2013) show that our galaxy formation implementation also produces close analogues of Milky Way-like disc galaxies in high-resolution simulations of Milky Way-sized DM haloes. We therefore conclude that we have implemented a successful galaxy formation model in the new moving-mesh code arepo, which is suitable for large-scale galaxy formation simulations. Clearly, the simulations presented in this paper are not sufficient to fully explore all regimes of galaxy formation. Most importantly, our simulation volume is too small to sample very high mass systems and ideally we would like to have slightly higher mass and spatial resolution. The convergence studies presented in this paper show that our model requires a mass resolution of at least 10⁶ h⁻¹ M_⊙ for the baryonic component to achieve reasonably converged results. Therefore, predictive large-scale structure simulations with our model require a simulation volume of the order of (100 h⁻¹ Mpc)³ with a mass resolution of the order of 10⁶ h⁻¹ M_⊙, which implies that the simulation has to follow the evolution of many billions of resolution elements. In forthcoming work, we will present such simulations based on the methodology presented here (Vogelsberger et al., in preparation).

The simulations in this paper were run on the Odyssey cluster supported by the FAS Science Division Research Computing Group at Harvard University, the Stampede supercomputer at the Texas Advanced Computing Center, the Magny Cluster at HITS and the CURIE supercomputer at CEA/France as part of PRACE project RA0844. We thank Simeon Bird, Laura Blecha, Charlie Conroy, Daniel Eisenstein, Claude-Andre Faucher-Giguere, Chris Hayward, Dusan Keres, Lisa Kewley, Bence Kocsis, Federico Marinacci, Diego Munoz, Dylan Nelson, Rüdiger Pakmor, Christoph Pfrommer, Ewald Puchwein, Vicente Rodriguez-Gomez, Laura Sales, Rob Simcoe, Gregory Snyder, Joshua Suresh, Rob Wiersma, Dandan Xu and Jabran Zahid for useful discussions. MV acknowledges support from NASA through Hubble Fellowship grant HST-HF-51317.01. VS acknowledges support through SFB 881 ‘The Milky Way System’ of the DFG, and through the European Research Council under ERC-StG grant EXAGAL-308037.

REFERENCES