Introducing the NEWHORIZON simulation: Galaxy properties with resolved internal dynamics across cosmic time

Yohan Dubois; Ricarda Beckmann; Frédéric Bournaud; Hoseung Choi; Julien Devriendt; Ryan Jackson; Sugata Kaviraj; Taysun Kimm; Katarina Kraljic; Clotilde Laigle; Garreth Martin; Min-Jung Park; Sébastien Peirani; Christophe Pichon; Marta Volonteri; Sukyoung K. Yi

doi:10.1051/0004-6361/202039429

Home

All issues

Volume 651 (July 2021)

A&A, 651 (2021) A109

Full HTML

Open Access

Issue		A&A Volume 651, July 2021


Article Number		A109
Number of page(s)		29
Section		Cosmology (including clusters of galaxies)
DOI		https://doi.org/10.1051/0004-6361/202039429
Published online		28 July 2021

A&A 651, A109 (2021)

Introducing the NEWHORIZON simulation: Galaxy properties with resolved internal dynamics across cosmic time

Yohan Dubois¹, Ricarda Beckmann¹, Frédéric Bournaud²^,3, Hoseung Choi⁴, Julien Devriendt⁵, Ryan Jackson⁶, Sugata Kaviraj⁶, Taysun Kimm⁴, Katarina Kraljic⁷^,8, Clotilde Laigle¹, Garreth Martin⁹^,10, Min-Jung Park⁴, Sébastien Peirani¹¹^,1, Christophe Pichon¹^,12^,13, Marta Volonteri¹ and Sukyoung K. Yi⁴

¹ Institut d’Astrophysique de Paris, UMR 7095, CNRS, UPMC Univ. Paris VI, 98bis boulevard Arago, 75014 Paris, France
e-mail: dubois@iap.fr
² AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France
³ IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France
⁴ Department of Astronomy and Yonsei University Observatory, Yonsei University, Seoul 03722, Republic of Korea
e-mail: yi@yonsei.ac.kr
⁵ Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
⁶ Centre for Astrophysics Research, University of Hertfordshire, College Lane, Hatfield, Herts AL10 9AB, UK
⁷ Aix Marseille Université, CNRS, CNES, UMR 7326, Laboratoire d’Astrophysique de Marseille, Marseille, France
⁸ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
⁹ Steward Observatory, University of Arizona, 933 N. Cherry Ave, Tucson, AZ 85719, USA
¹⁰ Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea
¹¹ Université Côte d’ Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Nice, France
¹² IPHT, DRF-INP, UMR 3680, CEA, Orme des Merisiers, Bât 774, 91191 Gif-sur-Yvette, France
¹³ Korea Institute of Advanced Studies (KIAS) 85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic of Korea

Received: 14 September 2020
Accepted: 5 May 2021

Abstract

Hydrodynamical cosmological simulations are increasing their level of realism by considering more physical processes and having greater resolution or larger statistics. However, usually either the statistical power of such simulations or the resolution reached within galaxies are sacrificed. Here, we introduce the NEWHORIZON project in which we simulate at high resolution a zoom-in region of ∼(16 Mpc)³ that is larger than a standard zoom-in region around a single halo and is embedded in a larger box. A resolution of up to 34 pc, which is typical of individual zoom-in, up-to-date resimulated halos, is reached within galaxies; this allows the simulation to capture the multi-phase nature of the interstellar medium and the clumpy nature of the star formation process in galaxies. In this introductory paper, we present several key fundamental properties of galaxies and their black holes, including the galaxy mass function, cosmic star formation rate, galactic metallicities, the Kennicutt–Schmidt relation, the stellar-to-halo mass relation, galaxy sizes, stellar kinematics and morphology, gas content within galaxies and its kinematics, and the black hole mass and spin properties over time. The various scaling relations are broadly reproduced by NEWHORIZON with some differences with the standard observables. Owing to its exquisite spatial resolution, NEWHORIZON captures the inefficient process of star formation in galaxies, which evolve over time from being more turbulent, gas rich, and star bursting at high redshift. These high-redshift galaxies are also more compact, and they are more elliptical and clumpier until the level of internal gas turbulence decays enough to allow for the formation of discs. The NEWHORIZON simulation gives access to a broad range of galaxy formation and evolution physics at low-to-intermediate stellar masses, which is a regime that will become accessible in the near future through surveys such as the LSST.

Key words: galaxies: general / galaxies: evolution / galaxies: stellar content / galaxies: kinematics and dynamics / methods: numerical

© Y. Dubois et al. 2021

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

1. Introduction

The origin of the various physical properties of galaxies, such as their mass content, size, kinematics, or morphology, emerges from the complex multi-scale and highly non-linear nature of the problem. It involves a strong connection between the small-scale star formation embedded in large molecular complexes and the gas that is accreted from the intergalactic medium and ejected into large-scale galactic outflows. To draw a theoretical understanding of the process of galaxy formation and evolution, it is necessary to connect cosmological structure formation – which leads to gas accretion into galaxies, that is the fuel of star formation – to the relevant small-scale processes that lead to the formation of the stars. Therefore, cosmological simulations are now a key tool in this theoretical understanding by allowing us to track the anisotropic non-linear cosmic accretion (which spectacularly results in filamentary gas accretion; e.g., Kereš et al. 2005; Dekel & Birnboim 2006; Ocvirk et al. 2008) in a self-consistent fashion.

Important challenges exist in the field of galaxy formation that need to be addressed, such as the global inefficiency of the star formation process on galactic scales (e.g., Moster et al. 2013; Behroozi et al. 2013), the morphological diversity of galaxies across the whole mass range (e.g., Conselice 2006; Martin et al. 2020), and the important evolution of the nature of galaxies over time; galaxies are more gas rich (e.g., Daddi et al. 2010a) and turbulent (e.g., Kassin et al. 2007), clumpy and irregular (e.g., Genzel et al. 2011), and star forming (e.g., Elbaz et al. 2007) at early time than they are in the local Universe.

High-redshift galaxies substantially differ in nature from low-redshift galaxies because cosmic accretion is more efficiently funnelled to the centre of dark matter (DM) halos owing to higher large-scale densities (Dekel et al. 2009), bringing gas into galaxies with lower angular momentum, higher surface densities, and, hence, more efficient star formation. However, for this high-redshift Universe that is naturally more efficient at feeding intergalactic gas into structures, a significant amount of galactic-scale feedback has to regulate the gas budget. On the low-mass end, it is generally accepted that stellar feedback as a whole, and more likely feedback from supernovae (SNe), is able to efficiently drive large-scale galactic winds (e.g., Dekel & Silk 1986; Springel & Hernquist 2003; Dubois & Teyssier 2008; Dalla Vecchia & Schaye 2008), although the exact strength of that feedback, and hence, how much gas is driven in and out of galaxies is still largely debated and relies on several important physical assumptions (e.g., Hopkins et al. 2012; Agertz et al. 2013; Kimm et al. 2015; Rosdahl et al. 2017; Dashyan & Dubois 2020). On the high-mass end, because of deeper potential wells, stellar feedback remains largely inefficient and gas regulation relies on the activity of central supermassive black holes (e.g., Silk & Rees 1998; Di Matteo et al. 2005; Croton et al. 2006; Dubois et al. 2010, 2012; Kaviraj et al. 2017; Beckmann et al. 2017).

Low-mass and low surface-brightness regimes are becoming important frontiers for the study of galaxy evolution (e.g., Martin et al. 2019) as surveys such as the LSST will allow us to observe very faint structures such as tidal streams and, for the first time, thousands of dwarfs at cosmological distances (mostly at z < 0.5). Complementary high-resolution cosmological simulations and deep observational datasets will enable us to start addressing the considerable tension between theory and observations in the dwarf regime (e.g., Boylan-Kolchin et al. 2011; Pontzen & Governato 2012; Naab & Ostriker 2017; Silk 2017; Kaviraj et al. 2019; Jackson et al. 2021a) as well as in the high-mass regime, where faint tidal features encode information that can aid in understanding the role of galaxy mergers and interactions in the formation, evolution, and survival of discs (Jackson et al. 2020; Park et al. 2019) and spheroids (Toomre & Toomre 1972; Bournaud et al. 2007; Naab et al. 2009; Kaviraj 2014; Dubois et al. 2016; Martin et al. 2018a).

Owing to their modelling of the most relevant aspects of feedback, SNe and supermassive black holes, which occur at the two mass ends of galaxy evolution, respectively, and thanks to their large statistics, large-scale hydrodynamical cosmological simulations with box sizes of ∼50 − 300 Mpc have made a significant step towards a more complete understanding of the various mechanisms (accretion, ejection, and mergers) involved in the formation and evolution of galaxies; these large-scale simulations include Horizon-AGN (Dubois et al. 2014a), Illustris (Vogelsberger et al. 2014), EAGLE (Schaye et al. 2015), IllustrisTNG (Pillepich et al. 2018), SIMBA (Davé et al. 2019), Extreme-Horizon (Chabanier et al. 2020a), and Horizon Run 5 (Lee et al. 2021). However, as a result of their low spatial resolution in galaxies (typically of the order of 1 kpc), and therefore owing to their intrinsic inability to capture the multi-phase nature of the interstellar medium (ISM), their sub-grid models for star formation or the coupling of feedback to the gas has had to rely on cruder effective approaches than what a higher-resolution simulation might allow. A couple of simulations with an intermediate volume and a better mass and spatial resolution stand out; these are the TNG50 simulation (Pillepich et al. 2019) from the IllustrisTNG suite and the Romulus25 simulation (Tremmel et al. 2017), which offer sub-kiloparsec resolution of 100 and 250 pc, respectively.

An important aspect of the evolution of galaxies is that rather than occurring in a homogeneous medium of diffuse interstellar gas, star formation proceeds within clustered molecular complexes; these range from pc to 100 pc in size and have properties that vary from one galaxy to another (e.g., Hughes et al. 2013; Sun et al. 2018). This has several important consequences. A clumpier star formation affects the stellar distribution via a more efficient migration of stars; it can be locally efficient while globally inefficient, and it can also enhance the effect of stellar feedback by driving more concentrated input of energy. Therefore, the necessity of capturing this minimal small-scale clustering of gas in galaxies has constrained numerical simulations to either rely on isolated set-ups (i.e., an isolated disc of gas and stars or isolated spherical collapsing halos; see e.g., Dobbs et al. 2011; Bournaud et al. 2014; Semenov et al. 2018) or on zoomed-in cosmological simulations with a handful of objects (e.g., Ceverino et al. 2010; Hopkins et al. 2014, 2018; Dubois et al. 2015; Nuñez-Castiñeyra et al. 2021; Agertz et al. 2021); this is because of the strong requisite on spatial resolution, that is typically below the 100 pc scale. Since star formation occurs in molecular clouds that are gravitationally bound or marginally bound with respect to turbulence, a consistent theory of a gravo-turbulence-driven star formation efficiency can be built considering that this shapes the probability density function (PDF) of the gas density within the cloud (see e.g., Federrath & Klessen 2012, and references therein). Such a theory can only be used in simulations in which the largest-scale modes of the interstellar medium turbulence are captured (Hopkins et al. 2014; Kimm et al. 2017; Nuñez-Castiñeyra et al. 2021). Similarly, less ad hoc models for SN feedback can be used to accurately reproduce the distinct physical phases of the blown-out SN bubbles (the so-called Sedov and snowplough phases; e.g., Kimm & Cen 2014), depending on the exact location of these explosions in the multi-phase ISM.

Our approach in this new numerical hydrodynamical cosmological simulation called NEWHORIZON, which we introduce in this work¹, is to provide a complementary tool between these two standard techniques, that is between the few well-resolved objects vs. a large ensemble of poorly resolved galaxies. The NEWHORIZON tool is designed to capture the basic features of the multi-scale, clumpy, ISM with a spatial resolution of the order of 34 pc in a large enough high-resolution, zoomed-in volume of (16 Mpc)³. This is larger than a standard zoomed-in halo, has a standard cosmological mean density, and is embedded in the initial lower-resolution (142 Mpc)³ volume of the Horizon-AGN simulation (Dubois et al. 2012); at z = 0.25, the mass density in that zoom-in region is 1.2 times that of the cosmic background density. Although still limited in terms of statistics over the entire range of galaxy masses (in particular galaxies in clusters are not captured), this volume offers sufficient enough statistics – in an average density region – to meaningfully study the evolution of galaxy properties at a resolution sufficient to apply more realistic models of star formation and feedback.

This paper introduces the NEWHORIZON simulation with its underlying physical model and reviews the main fundamental properties of the simulated galaxies, including their mass budget, star formation rate (SFR), morphology, kinematics, and the mass and spin properties of the hosted black holes in galaxies.

The paper is organised as follows. Section 2 presents the numerical technique, resolution, and physical models implemented in NEWHORIZON. Section 3 presents the various results of the properties of the galaxies in the simulation and their evolution over time. Finally, we conclude in Sect. 4.

2. The NEWHORIZON simulation: Prescription

We describe the NEWHORIZON, simulation employed in this work², which is a sub-volume extracted from its parent HORIZON-AGN simulation (Dubois et al. 2014a)³, and the procedure we use to identify halos and galaxies. A number of physical sub-grid models have been substantially modified compared to the physics implemented in HORIZON-AGN (see e.g., Volonteri et al. 2016; Kaviraj et al. 2017), in particular regarding the models for star formation, feedback from SNe and from active galactic nuclei (AGN). A comparison with simulated galaxies in HORIZON-AGN within the same sub-volume will be the topic of a dedicated paper. Nonetheless, we describe the corresponding differences with HORIZON-AGN at the end of each of the subsections of the sub-grid model.

2.1. Initial conditions and resolution

The NEWHORIZON, simulation is a zoom-in simulation from the 142 Mpc size HORIZON-AGN simulation (Dubois et al. 2014a). The HORIZON-AGN simulation initial conditions had 1024³ DM particles, a 1024³ minimum grid resolution, and a ΛCDM cosmology. The total matter density is Ω_m = 0.272, dark energy density Ω_Λ = 0.728, amplitude of the matter power spectrum σ₈ = 0.81, baryon density Ω_b = 0.045, Hubble constant H₀ = 70.4 km s⁻¹ Mpc⁻¹, and n_s = 0.967 is compatible with the WMAP-7 data (Komatsu et al. 2011). Within this large-scale box, we define an initial spherical patch of 10 Mpc radius, which is large enough to sample multiple halos at a 4096³ effective resolution, that is with a DM mass resolution of M_DM, hr = 1.2 × 10⁶ M_⊙. The high-resolution initial patch is embedded in buffered regions with decreasing mass resolution of 10⁷ M_⊙, 8 × 10⁷ M_⊙, 6 × 10⁸ M_⊙ for spheres of 10.6 Mpc, 11.7 Mpc, and 13.9 Mpc radius, respectively, and a resolution of 5 × 10⁹ M_⊙ in the rest of the simulated volume. In order to follow the Lagrangian evolution of the initial patch, we fill this initial sub-volume with a passive colour variable with values of 1 inside and zero outside, and we only allow for refinement when this passive colour is above a value of 0.01. Within this coloured region, refinement is allowed in a quasi-Lagrangian manner down to a resolution of Δx = 34 pc at z = 0: refinement is triggered if the total mass in a cell becomes greater than eight times the initial mass resolution. The minimum cell size is kept roughly constant by adding an extra level of refinement every time the expansion factor is doubled (i.e., at a_exp = 0.1, 0.2, 0.4 and 0.8); the minimum cell size is thus between Δx = 27 and 54 pc. We also added a super-Lagrangian refinement criterion to enforce the refinement of the mesh if a cell has a size shorter than one Jeans’ length wherever the gas number density is larger than 5 H cm⁻³.

The NEWHORIZON simulation is run with the adaptive mesh refinement RAMSES code (Teyssier 2002). Gas is evolved with a second-order Godunov scheme and the approximate Harten-Lax-Van Leer-Contact (HLLC, Toro 1999) Riemann solver with linear interpolation of the cell-centred quantities at cell interfaces using a minmod total variation diminishing scheme. Time steps are sub-cycled on a level-by-level basis, that is each level of refinement has a time step that is twice as small as the coarser level of refinement, following a Courant–Friedrichs–Lewy condition with a Courant number of 0.8. The simulation was run down to z = 0.25 using a total amount of 65 single core central processing unit (CPU) million hours. The simulation contained typically 0.5–1 billion of leaf cells in total. With 30–100 millions of leaf cells per level of refinement in the zoom-in region from level 12 to level 22, the region had a total of 3.3 × 10⁸ star particles formed and completed 4.7 × 10⁶ fine time steps (the number of time steps of the maximum level of refinement), thus, corresponding to an average fine time step of size Δt ≃ 2.3 kyr), by z = 0.25.

Figure 1 shows a projection of the high-resolution region. Figure 2 illustrates the typical structure of the gas density achieved in one of the massive galaxies at z = 1 and the corresponding gas resolution. The diffuse ISM (0.1–1 cm⁻³) is resolved with a ∼100 pc resolution or such, while the densest clouds reach the maximum level of refinement corresponding to 34 pc and the immediate galactic corona is resolved with cells of size 500 pc. In terms of mass and spatial resolution, NEWHORIZON is comparable to TNG50 (Pillepich et al. 2019; 4.5 × 10⁵ M_⊙ DM mass resolution and a spatial resolution in galaxies of 100 pc) or zoomed-in cosmological simulations (such as for the most massive galaxies of the FIRE-2 runs; Hopkins et al. 2018).

Fig. 1.

Sequential zoom (clockwise from top left) over the projected density (silver blue colours) and projected temperature (red) of the NEWHORIZON simulation at redshift z = 2. The dashed white circles encompass the initial high-resolution volume. Each panel is a zoomed-in version of the previous panel (identified by the white square in the previous panel) with the panel sizes of 142, 18, 4.4, and 1.1 comoving Mpc width, respectively. Two top panels: encompass the zoom-in region, with its network of filaments. Two bottom panels: how narrow filaments break up and mix once they connect to one of the most massive galaxies of that zoom-in region.

Fig. 2.

Illustration of the structure of the gas density (left panels) and the corresponding spatial resolution (right panels) in a massive galaxy of M_s = 6 × 10¹⁰ M_⊙ at z = 1 seen edge-on (top panels) or face-on (bottom panels).

2.2. Radiative cooling and heating

We adopt the equilibrium chemistry model for primordial species (H and He) assuming collisional ionisation equilibrium in the presence of a homogeneous UV background. The primordial gas is allowed to cool down to ≈10⁴ K through collisional ionisation, excitation, recombination, Bremsstrahlung, and Compton cooling. Metal-enriched gas can cool further down to 0.1 K using rates tabulated by Sutherland & Dopita (1993) above ≈10⁴ K and those from Dalgarno & McCray (1972) below ≈10⁴ K. The heating of the gas from a uniform UV background takes place after redshift z_reion = 10, following Haardt & Madau (1996). Motivated by the radiation-hydrodynamic simulation results that the UV background is self-shielded in optically thick regions (n_H ≳ 0.01 H cm⁻³; Rosdahl & Blaizot 2012), we assume that UV photo-heating rates are reduced by a factor exp(−n_H/n_shield), where n_shield = 0.01 H cm⁻³.

Compared to HORIZON-AGN, the gas can now cool below 10⁴ K.

2.3. Star formation

Star formation occurs in regions with hydrogen gas number density above n₀ = 10 H cm⁻³ (the stellar mass resolution is n₀m_pΔx³ = 1.3 × 10⁴ M_⊙) following a Schmidt law: ${\dot{ρ}}_{⋆} = ϵ_{⋆} ρ_{g} / t_{ff}$ $\dot \rho_\star= \epsilon_\star {\rho_{\mathrm{g}} / t_{\mathrm{ff}}}$ , where ${\dot{ρ}}_{⋆}$ $\dot \rho_\star$ is the SFR mass density, ρ_g the gas mass density, $t_{ff} = \sqrt{3 π / (32 G ρ_{g})}$ $t_{\mathrm{ff}}=\sqrt{3\pi/(32G\rho_{\mathrm{g}})}$ the local free-fall time of the gas, G the gravitational constant, and ϵ_⋆ ia varying star formation efficiency (Kimm et al. 2017; Trebitsch et al. 2017, 2020).

The current theory of star formation provides a framework for working out the efficiency of the star formation where the gas density PDF is well approximated by a log-normal PDF (Krumholz & McKee 2005; Padoan & Nordlund 2011; Hennebelle & Chabrier 2011; Federrath & Klessen 2012). This PDF is related to the star-forming cloud properties through the cloud turbulent Mach number ℳ = u_rms/c_s, where u_rms is the root mean square velocity, c_s the sound speed. The virial parameter α_vir = 2E_kin/E_grav and the efficiency is fully determined by integrating how much mass passes above a given density threshold using the multi-free fall approach of Hennebelle & Chabrier (2011) as follows:

$\begin{matrix} ϵ_{⋆} = \frac{ϵ}{2 ϕ_{t}} exp (\frac{3}{8} σ_{s}^{2}) [1 + \erf (\frac{σ_{s}^{2} - s_{crit}}{\sqrt{2 σ_{s}^{2}}})], \end{matrix}$ $\begin{aligned} \epsilon _{\star } = \frac{\epsilon }{2\phi _{t}} \exp \left(\frac{3}{8}\sigma _{\rm s}^2\right)\left[1 + \mathrm{erf} \left(\frac{\sigma _{\rm s}^2 - s_{\rm crit}}{\sqrt{2\sigma _{\rm s}^2}}\right)\right], \end{aligned}$ (1)

where s = ln(ρ/ρ₀) is the logarithmic density contrast of the PDF with mean ρ₀ and variance $σ_{s}^{2} = \ln (1 + b^{2} M^{2})$ $\sigma_{\rm s}^2=\ln (1+b^2\mathcal{M}^2)$ . In this expression b = 0.4 conveys the fractional amount of solenoidal to compressional modes of the turbulence. The critical density contrast s_crit is determined by Padoan & Nordlund (2011) as follows:

$\begin{matrix} s_{crit} = ln (\frac{0.067}{θ^{2}} α_{vir} M^{2}) . \end{matrix}$ $\begin{aligned} s_{\rm crit} = \ln \left(\frac{0.067}{\theta ^2} \alpha _{\rm vir}\mathcal{M} ^2\right). \end{aligned}$ (2)

In the NEWHORIZON simulation, the turbulent Mach number is given by the local three-dimensional instantaneous velocity dispersion σ_g (obtained by computing $σ_{g}^{2} = sum {(\nabla \otimes u d x)}^{2}$ $\sigma_{\rm g}^2={\rm sum} (\nabla \otimes u {\rm d}x)^2$ ), and the virial parameter also takes the thermal pressure support $α_{vir, 0} = 5 (σ_{g}^{2} + c_{s}^{2}) / (π ρ_{g} G Δ x^{2})$ $\alpha_{\rm vir,0}=5(\sigma_{\rm g}^2+c_{\rm s}^2)/(\pi\rho_{\rm g}G\Delta x^2)$ into account. In this case, $ϕ_{t}^{- 1} = 0.57$ $\phi_{t}^{-1}=0.57$ and θ = 0.33 are empirical parameters of the model determined by the best-fit values between the theory and the numerical experiments (Federrath & Klessen 2012). The different values of $ϕ_{t}^{- 1}$ $\phi_{t}^{-1}$ and θ we use compared to those given in Federrath & Klessen (2012) arise from the difference between the definition of α_vir (measured over time, which are the values given in Federrath & Klessen 2012) and α_vir, 0 (the homogeneous cloud initial conditions). As our measurements of the virial parameter are meant to correspond to the initial cloud value α_vir, 0, that is to the virial parameter of a spherical gas cloud with the same mass, radius, and thermo-turbulent velocity dispersion (Bertoldi & McKee 1992; Krumholz & McKee 2005) of the gas cell, we use the best-fit values from Federrath & Klessen (2012) corresponding to this definition of the virial parameter (Fedderath, private communication). We ignore the role of the magnetic field in this model despite the effect it has on the critical density and variance of the density PDF due to its large pressure with respect to the thermal pressure in the cold neutral medium (e.g., Heiles & Troland 2005; Crutcher 2012). In Eq. (1) ϵ = 0.5 is a proto-stellar feedback parameter that controls the actual amount of gas above s_crit that is able to form stars (typical estimates of ϵ are around 0.3 − 0.5; see Matzner & McKee 2000; Alves et al. 2007; André et al. 2010).

Such a star formation law shows a significantly different behaviour on galactic scales with respect to simulations with constant (usually low) efficiencies since the efficiency can now vary by orders of magnitude. For instance, for gravitationally bound (α_vir < 1) and highly turbulent regions (ℳ > 1), the efficiency can go well above 1, while regions that are marginally bound have an efficiency that quickly drops to very low values. Star formation efficiency, in conjunction with stellar feedback, plays a key role in shaping galaxy properties (e.g., Agertz et al. 2011; Nuñez-Castiñeyra et al. 2021), and such potentially higher and more bursty star formation participates in driving stronger outflows and self-regulation of galaxy properties. We note that our gravo-turbulent model of SF is somewhat reminiscent of those adopted in Hopkins et al. (2014, 2018) or Semenov et al. (2018). In those models, α_vir is used as a criterion to trigger star formation (gas needs to be sufficiently bound), but star formation proceeds with a constant efficiency in contrast to our model.

Compared to HORIZON-AGN, star formation in NEW HORIZON occurs at above a hundred times larger gas density, and a varying gravo-turbulent-based star formation efficiency is used instead of assuming a constant 2% efficiency.

2.4. Feedback from massive stars

We include feedback from Type II SNe assuming that each explosion initially releases the kinetic energy of 10⁵¹ erg. Because the minimum mass of a star particle is 10⁴ M_⊙, each particle is assumed to represent a simple stellar population with a Chabrier initial mass function (IMF; Chabrier 2005) where the lower (upper) mass cut-off is taken as M_low = 0.1 (M_upp = 150) M_⊙, respectively. We further assume that the minimum mass that explodes is 6 M_⊙ in order to include electron-capture SNe (Chiosi et al. 1992, see also Crain et al. 2015). The corresponding specific frequency of SN explosion is $0.015 M_{⊙}^{- 1}$ $0.015 \,M_\odot^{-1}$ . We increase this number by a factor of 2 ( $0.03 M_{⊙}^{- 1}$ $0.03 \,M_\odot^{-1}$ ) because multiple clustered SN explosions can increase the total radial momentum, with respect to the total momentum predicted by the accumulation of individual SNe (Thornton et al. 1998), by decreasing the ambient density into which subsequent SNe explode (Kim et al. 2017; Gentry et al. 2019; Na et al., in prep.). Supernovae are assumed to explode instantaneously when a star particle becomes older than 5 Myr. The mass loss fraction of a stellar particle from the explosions is 31% and has a metal yield (mass ratio of the newly formed metals over the total ejecta) of 0.05.

We employ the mechanical SN feedback scheme (Kimm & Cen 2014; Kimm et al. 2015), which ensures the transfer of a correct amount of radial momentum to the surroundings. Specifically, the model examines whether the blast wave is in the Sedov–Taylor energy-conserving or momentum-conserving phase (Chevalier 1974; Cioffi et al. 1988; Blondin et al. 1998) by calculating the mass swept up by SN. If the SN explosion is still in the energy-conserving phase, the assumed specific energy is injected into the gas since hydrodynamics naturally capture the expansion of the SN and imparts the correct amount of radial momentum. However, if the cooling length in the neighbouring regions is under-resolved owing to finite resolution, radiative cooling takes place rapidly, thereby suppressing the expansion of the SN bubble. This leads to an under-estimation of the radial momentum, hence weaker feedback. In order to avoid this artificial cooling, the mechanical feedback model directly imparts the radial momentum expected during the momentum-conserving phase if the mass of the neighbouring cell exceeds some critical value. This is done by first measuring the local ratio of the swept-up gas mass over the ejecta mass and examining whether the ratio is greater than the critical ratio corresponding to the energy-to-momentum phase transition. That is to say $70 E_{51}^{- 2 / 17} n_{1}^{- 4 / 17} Z'^{- 0.28}$ $70 \,E_{\mathrm{51}}^{-2/17}\,n_{\mathrm{1}}^{-4/17}\,Z{{\prime}}^{-0.28}$ , where E₅₁ is the total energy released in units of 10⁵¹ erg, n₁ is the hydrogen number density in units of cm⁻³, and Z′=max[Z/Z_⊙, 0.01] is the metallicity, normalised to the solar value (Z_⊙ = 0.02). The final momentum in the snowplough phase per SN explosion is taken from Thornton et al. (1998) as

$\begin{matrix} q_{SN} = 3 \times 10^{5} km s^{- 1} M_{⊙} E_{51}^{16 / 17} n_{1}^{- 2 / 17} Z'^{- 0.14} . \end{matrix}$ $\begin{aligned} q_{\rm SN}=3\times 10^5 \, \mathrm{km\,s^{-1}}\,M_\odot \, E_{51}^{16/17} n_{\rm 1}^{-2/17} Z{\prime }^{-0.14}. \end{aligned}$ (3)

We further assume that the UV radiation from the young OB stars over-pressurises the ambient medium near to young stars and increases the total momentum per SN to

$\begin{matrix} q_{SN + PH} = 5 \times 10^{5} km s^{- 1} M_{⊙} E_{51}^{16 / 17} n_{1}^{- 2 / 17} Z'^{- 0.14}, \end{matrix}$ $\begin{aligned} q_{\rm SN+PH}=5\times 10^5 \, \mathrm{km\,s^{-1}}\,M_\odot \, E_{51}^{16/17} n_{\rm 1}^{-2/17} Z{\prime }^{-0.14} , \end{aligned}$ (4)

following Geen et al. (2015).

It is worth noting that the specific energy used for SN II explosion in this study is larger than previously assumed. A Chabrier (2003) IMF with a low- to high-mass cut-off of M_low = 0.1 and M_upp = 100 M_⊙ and an intermediate-to-massive star transition mass at M_IM = 8 M_⊙ gives $e_{SN} = 1.1 \times 10^{49} erg M_{⊙}^{- 1}$ $e_{\rm SN}=1.1\times 10^{49}\,{\rm erg}\,M_\odot^{-1}$ . However, e_SN can be increased up to $3.6 \times 10^{49} erg M_{⊙}^{- 1}$ $3.6\times 10^{49}\,{\rm erg}\, M_\odot^{-1}$ if a non-negligible fraction (f_HN = 0.5) of hypernovae (with E_HN ≃ 10⁵² erg for stars more massive than 20 M_⊙; e.g., Iwamoto et al. 1998; Nomoto et al. 2006) is taken into account. This is necessary to reproduce the abundance of heavy elements, such as zinc (Kobayashi et al. 2006), or if a lower transition mass M_IM = 6 M_⊙ and a shallower (Salpeter) slope of −2.1 at the high-mass end (reflecting that early star formation should lead to a top-heavier IMF; e.g., Treu et al. 2010; Cappellari et al. 2012; Martín-Navarro et al. 2015) are assumed. Furthermore, various sources of stellar feedback that would contribute to the overall formation of large-scale outflows including type Ia SNe, stellar winds, shock-accelerated cosmic rays (e.g., Uhlig et al. 2012; Salem & Bryan 2014; Dashyan & Dubois 2020), multi-scattering of infrared photons with dust (e.g., Hopkins et al. 2011; Roškar et al. 2014; Rosdahl & Teyssier 2015), or Lyman-α resonant line scattering (Kimm et al. 2018; Smith et al. 2017) are neglected. In addition runaway OB stars (Ceverino & Klypin 2009; Kimm & Cen 2014; Andersson et al. 2020) or the unresolved porosity of the medium (Iffrig & Hennebelle 2015) are also ignored. In this regard, the NEWHORIZON simulation is unlikely to overestimate the effects of stellar feedback, as described in Sect. 3.

Unlike HORIZON-AGN, feedback from stars in NEWHORIZON only includes Type II SNe and ignores stellar winds and Type Ia SNe. In addition, NEWHORIZON adopts a mechanical scheme for SNe instead of a kinetic solution (Dubois & Teyssier 2008). The assumed IMF is also changed from the Salpeter IMF to a Chabrier type, and thus the mass loss, energy, and yield are all increased.

2.5. MBHs and AGN

We now briefly describe the models corresponding to massive black hole (MBH) formation and their AGN feedback.

2.5.1. Formation, growth, and dynamics of MBH

In NEWHORIZON, MBHs are assumed to form in cells that have gas and stellar densities above the threshold for star formation, a stellar velocity dispersion larger than 20 km s⁻¹, and that are located at a distance of at least 50 comoving kpc from any pre-existing MBH.

Once formed, the mass of MBHs grows at a rate Ṁ_MBH = (1−ϵ_r)Ṁ_Bondi, where ϵ_r is the spin-dependent radiative efficiency (see Eq. (7)) and Ṁ_Bondi is the Bondi–Hoyle–Lyttleton rate, that is

$\begin{matrix} \frac{d M_{Bondi}}{d t} = 4 π \bar{ρ} \frac{{(G M_{MBH})}^{2}}{{({\bar{u}}^{2} + {\bar{c_{s}}}^{2})}^{3 / 2}}, \end{matrix}$ $\begin{aligned} \frac{\mathrm{d}M_{\rm Bondi}}{\mathrm{d}t}=4\pi \bar{\rho }\frac{(GM_{\rm MBH})^2}{(\bar{u}^2+\bar{c_{\rm s}}^2)^{3/2}}, \end{aligned}$ (5)

where $\bar{u}$ $\bar u$ is the average MBH-to-gas relative velocity, $\bar{c_{s}}$ $\bar{c_{\mathrm{s}}}$ the average gas sound speed, and $\bar{ρ}$ $\bar \rho$ the average gas density. All average quantities are computed within 4Δx of the MBH, using mass weighting and a kernel weighting as specified in Dubois et al. (2012). We do not employ a boost factor in the formulation of the accretion rate, as is commonly done in cosmological simulations, because we have sufficient spatial resolution to model part of the multi-phase structure of the ISM of galaxies directly.

The Bondi–Hoyle–Lyttleton accretion rate is capped at the Eddington luminosity rate for the appropriate ϵ_r

$\begin{matrix} \frac{d M_{Edd}}{d t} = \frac{4 π G M_{MBH} m_{p}}{ϵ_{r} σ_{T} c}, \end{matrix}$ $\begin{aligned} \frac{\mathrm{d}M_{\rm Edd}}{\mathrm{d}t}=\frac{4\pi G M_{\rm MBH} m_{\rm p}}{\epsilon _{\rm r}\sigma _{\rm T}c}, \end{aligned}$ (6)

where σ_T is the Thompson cross-section, m_p the proton mass, and c the speed of light.

To avoid spurious motions of MBHs around high-density gas regions as a result of finite force resolution effects, we include an explicit drag force of the gas onto the MBH, following Ostriker (1999). This drag force term includes a boost factor with the functional form α = (n/n₀)² when n > n₀, and α = 1 otherwise. The use of a sub-grid drag force model is justified by our larger-than-Bondi-radius spatial resolution (Beckmann et al. 2018). We also enforce maximum refinement within a region of radius 4Δx around the MBH, which improves the accuracy of MBH motions (Lupi et al. 2015).

The MBHs are allowed to merge when they get closer than 4Δx (∼150 pc) and when the relative velocity of the pair is smaller than the escape velocity of the binary. A detailed analysis of MBH mergers in NEWHORIZON is presented in Volonteri et al. (2020).

2.5.2. Spin evolution of MBH

The evolution of the spin parameter a is followed on-the-fly in the simulation, taking the effects of gas accretion and MBH–MBH mergers into account. The model of MBH spin evolution is introduced in Dubois et al. (2014c), and technical details of the model are detailed in that paper. The only change is that we now use a different MBH spin evolution model at low accretion rates: χ = Ṁ_MBH/Ṁ_Edd < χ_trans, where χ_trans = 0.01. At high accretion rates (χ ≥ χ_trans), a thin accretion disc solution is assumed (Shakura & Sunyaev 1973), as in Dubois et al. (2014c). The angular momentum direction of the accreted gas is used to decide whether the accreted gas feeds an aligned or misaligned Lense–Thirring disc precessing with the spin of the MBH (King et al. 2005), thereby spinning the MBH up or down for co-rotating and counter-rotating systems, respectively (see top panel of Fig. 3). At low accretion rates (χ < 0.01), we assume that jets are powered by energy extraction from MBH rotation (Blandford & Znajek 1977) and that the MBH spin magnitude can only decrease. The change in the spin magnitude da/dM follows the results from McKinney et al. (2012), where we fitted a fourth-order polynomial to their sampled values; from their Table 7, s_H for AaN100 runs, where a is the value of the MBH spin. The functional form of the spin evolution as a function of MBH spin at low accretion rates is represented in the top panel of Fig. 4, where the dimensionless spin-up parameter s ≡ d(a/M)/dt is shown, where if s and a have opposite signs the black hole spins down.

Fig. 3.

Spin-up rate (top panel) and radiative efficiency ϵ_rad/f_att (bottom panel) as a function of the MBH spin for the thin disc solution (Shakura & Sunyaev 1973) applied to the quasar mode. At negative values of the MBH spin, the gas accreted from the thin accretion disc decreases the MBH spin, while for a positive MBH spin, the gas increases. For the thin disc solution, the radiative efficiency is an increasing function of the MBH spin with a sharp increase (by 4) between a MBH spin of 0.7 and 0.998.

In addition, MBH spins change in magnitude and direction during MBH–MBH coalescences, with the spin of the remnant depending on the spins of the two merging MBHs and the orbital angular momentum of the binary, following analytical expressions from Rezzolla et al. (2008).

The evolution of the spin parameter is a key component of the AGN feedback model because it controls the radiative efficiency of the accretion disc and the jet efficiency. Therefore, the Eddington mass accretion rate, used to cap the total accretion rate, and the AGN feedback efficiency in the jet and thermal modes vary with spin values. The spin-dependent radiative efficiency (see bottom panel of Fig. 3) is defined as

$\begin{matrix} ϵ_{r} = f_{att} (1 - e_{isco}) = f_{att} (1 - \sqrt{1 - 2 / (3 r_{isco})}), \end{matrix}$ $\begin{aligned} \epsilon _{\rm r}=f_{\rm att}\left(1-e_{\rm isco} \right)=f_{\rm att}\left(1-\sqrt{1-2/(3r_{\rm isco})}\right) ,\end{aligned}$ (7)

where e_isco is the energy per unit rest mass energy of the innermost stable circular orbit (ISCO), r_isco = R_isco/R_g is the radius of the ISCO in reduced units, and R_g is half the Schwarzschild radius of the MBH. The parameter R_isco depends on spin a. For the radio mode, the radiative efficiency used in the effective growth of the MBH is attenuated by a factor f_att = min(χ/χ_trans, 1) following Benson & Babul (2009). The MBH seeds are initialised with a zero spin value and a maximum value of the BH spin at a_max = 0.998 (due to the emitted photons by the accretion disc captured by the MBH; Thorne 1974) is imposed.

2.5.3. Radio and quasar modes of AGN feedback

Active galactic nuclei feedback is modelled in two different ways depending on the Eddington rate (Dubois et al. 2012): below χ < χ_trans the MBH powers jets (a.k.a. radio mode) continuously releasing mass, momentum, and total energy into the gas (Dubois et al. 2010), while above χ ≥ χ_trans the MBH releases only thermal energy back into the gas (a.k.a. quasar mode, Teyssier et al. 2011). The AGN releases a power that is a fraction of the rest-mass accreted luminosity onto the MBH, L_AGN,R,Q = η_R,QṀ_MBHc², where the subscripts R and Q stand for the radio jet mode and quasar heating mode, respectively.

For the jet mode of AGN feedback, the efficiency η_R is not a free parameter. This value scales with the MBH spin, following the results from magnetically chocked accretion discs (MCAD) of McKinney et al. (2012), where we fitted a fourth-order polynomial to the sampled values of jet plus wind efficiencies of this work (from their Table 5, η_j plus η_w, 0 for runs AaN100). This fit is shown in the bottom panel of Fig. 4. When active in our simulation, the bipolar AGN jet deposits mass, momentum, and total energy within a cylinder of size Δx in radius and semi-height, centred on the MBH, whose axis is (anti)aligned with the MBH spin axis (zero opening angle). Jets are launched with a speed of 10⁴ km s⁻¹, whose exact value has little impact on MBH growth or galaxy mass content (Dubois et al. 2012).

Fig. 4.

Spin-up rate (top panel), jet (red plus signs), wind (blue diamonds), and total (black) efficiencies (bottom panel) as a function of the MBH spin for the MCAD solution applied to the radio mode. The symbols represent the results from the simulations of McKinney et al. (2012) and the solid lines indicate the interpolated functions used in the NEWHORIZON simulation. As opposed to the thin disc solution (Fig. 3), gas accreted from the thick accretion disc always decreases the MBH spin. The MCAD feedback efficiency is an increasing function of the absolute value of the BH spin with a minimum efficiency for non-rotating MBHs.

The quasar mode of AGN feedback deposits internal energy into its surrounding within a sphere of radius Δx, within which the specific energy is uniformly deposited (uniform temperature increase). Because only a fraction of the AGN-driven wind is expected to thermalise and only some of the multiwavelength radiation emitted from the accretion disc couples to the gas on ISM scales (Bieri et al. 2017), we scale the feedback efficiency in quasar mode by a coupling factor of η_c = 0.15, which is calibrated on the local M_MBH − M_s in lower resolution (∼kpc) simulations (Dubois et al. 2012). The effective feedback efficiency in quasar mode is therefore η_Q = ϵ_rη_c.

Compared to HORIZON-AGN, NEWHORIZON now includes MBH spin evolution, which affects several compartments of MBH mass growth and feedback. The MBH accretion is changed owing to the spin-dependent radiative efficiency, thereby changing the maximum Eddington accretion rate. The AGN feedback is also changed by the spin-dependent radiative efficiency in the quasar mode. For the radio mode, the jet closely follows the spin-dependent mechanical efficiency of the MCAD model instead of a constant efficiency of 1, and the jet direction is now along the BH spin axis instead of along the accreted gas angular momentum.

2.6. Identification of halos and galaxies

Halos are identified with the AdaptaHOP halo finder (Aubert et al. 2004). The density field used in AdaptaHOP is smoothed over 20 particles. The minimum number of particles in a halo is 100 DM particles. We only consider halos with an average overdensity with respect to the critical density ρ_c, which is larger than δ_t = 80 and which overcomes the Poissonian noise filtering density threshold at $(1 + 5 / \sqrt{N}) δ_{t} ρ_{c}$ $(1+5/\sqrt{N})\delta_{t}\rho_{\mathrm{c}}$ (where N is the number of particles in the (sub)structure; seeAubert et al. 2004, for details). For a substructure, it is only kept if the maximum density is 2.5 times its mean density. The centre of the halo is recursively determined by seeking the centre of mass in a shrinking sphere, while decreasing its radius by 10% recurrently down to a minimum radius of 0.5 kpc (Power et al. 2003). The maximum DM density in that radius is defined as the centre of the halo. The shrinking sphere approach is used since strong feedback processes can significantly flatten the central DM density and smaller, but denser, substructures can be misidentified as being the centre of the main halo.

We run the same identification technique, using either AdaptaHOP or HOP, on stars to identify the galaxies in the simulation, except that we only consider galaxies with more than 50 star particles and a value of δ_t twice as large. The AdaptaHOP tool separates substructures that include in situ star-forming clumps as well as satellites already connected to a galaxy, while HOP keeps all substructures connected to the main structure (i.e., it does not detect substructures). Appendix B shows examples of how using HOP or AdaptaHOP affect the segmentation of galaxies. Both tools can be employed depending on context, as indicated in the corresponding text. For the centring of the galaxies at the low-mass end, particular attention has to be taken, since these galaxies tend to be extremely turbulent structures where bulges cannot be easily identified.

Since the NEWHORIZON simulation is a zoom simulation embedded in a larger cosmological volume filled with lower DM resolution particles, we also need to remove halos of the zoom regions polluted with low-resolution DM particles. To that end, we only consider halos as well as the embedded galaxies and MBHs encompassed in their virial radius, which are found devoid of low-resolution DM particles up to some threshold (see Appendix A for the halo mass function for different purity levels). With 100% purity, there are, respectively, 626, 245, 53, and 5 main galaxies (which are not substructures in the sense of AdaptaHOP) at z = 2 with stellar mass above 10⁷, 10⁸, 10⁹, and 10¹⁰ M_⊙; 403, 191, 70, and 12 at z = 1; and 276, 145, 58, and 16 at z = 0.25. For comparison, considering a contamination lower than 1 per cent in number of DM, the number of galaxies typically doubles at z = 0.25 (see Table 1 for detailed numbers). We note that the most massive unpolluted halo obtained at z = 0.25 has a DM virial mass of 8 × 10¹² M_⊙.

Table 1.

Number of galaxies for different stellar mass thresholds.

3. Cosmic evolution of baryons

In this section we present several standard properties of the simulated galaxies including their stellar and gas mass content, SFR, morphological and structural properties, and kinematics. We also present their hosted MBHs and compare these to observational relations down to the lowest redshift reached out by the simulation (z = 0.25).

3.1. Synthetic galaxy morphology

In order to qualitatively illustrate the variety of galaxy properties simulated in NEWHORIZON, we show in Fig. 5 a couple of galaxies at z = 4, z = 2 and z = 0.25 with their gas density and stellar emission. The 15 panels on the left show the images of a massive galaxy (stellar mass M_s = 3.0 × 10⁸ M_⊙ at z = 4, 8.2 × 10⁹ M_⊙ at z = 2 and 5.5 × 10¹⁰ M_⊙ at z = 0.25) and the 15 panels on the right represent a less massive galaxy (9.7 × 10⁷ M_⊙ at z = 4, 1.4 × 10⁹ M_⊙ at z = 2 and 4.0 × 10⁹ M_⊙ at z = 0.25). While the first, second, and fourth rows show their gas density maps, the third and fifth rows show the mock images; the second and third rows are shown with a face-on view (with respect to the stellar angular momentum of the galaxy) and the fourth and fifth rows an edge-on view. The mock images are in SDSS g − r − i bands and are generated via the SKIRT9 code (Camps & Baes 2020), which computes radiative transfer effects based on the properties and positions of stars and the dusty gas assuming a dust fraction f_dust = 0.4 following Saftly et al. (2015). The high resolution of NEWHORIZON (34 pc) reveals the detailed structure of the cosmologically simulated galaxies, and it is clearly evident that star formation (highlighted by the young blue region in the stellar maps) proceeds in clustered regions of dense gas. The massive galaxy settles its disc around z ≈ 2.5 and appears as a regular disc galaxy with well-defined spiral arms and a central bulge if witnessed at z = 0.25. We used the visual inspection as well as (V/σ)_gas > 3 (Kassin et al. 2012) for disc settling criteria; the calculation of the kinematics is detailed in Sect. 3.13. The less massive galaxy, on the other hand, exhibits an extremely irregular morphology at z = 2 with strong asymmetries in both gas and stars, and prominent off-centred (blue) star-forming clusters. This low-mass galaxy, which only grows moderately by z = 0.25, develops a galactic-scale disc at z ≈ 1.0 and maintains the marginally stable disc for the rest of the cosmic history.

Fig. 5.

Projection of the gas density and mock observations. Panels on the left are images of a massive galaxy at different epochs and the panels on the right are for a less massive galaxy. First row: gas density projections a of 1.4 Mpc at different epochs, while the second and fourth rows are zoomed-in gas density projections with face-on and edge-on views of the galaxy, respectively. Third and fifth rows are SKIRT mock observations in face-on and edge-on direction. Stellar mass and halo mass of each galaxy at each epoch are given in the second row in log scale. For a galaxy at a given epoch, the second to fifth panels are on the same scale and the white bar in third row indicates 5 kpc. Gas density maps share the same colour scheme as given in the colour bar.

3.2. Galaxy mass function

We compare the z = 0.25 mass function obtained from NEWHORIZON with the mass function obtained from an equivalent volume in the HSC-SSP survey (Aihara et al. 2019). In order to do this we take 100 random pointings from the HSC-SSP deep layer (encompassing the SXDS, COSMOS, ELIAS, and DEEP-2 fields), where each pointing has an equivalent volume to the NEWHORIZON box. The central redshift of each volume is varied by up to 0.02 around a central redshift of z = 0.25 for each random pointing. Since the photometric redshift errors are typically larger than the 20 Mpc box length, it is likely that we do not capture the full variance in the mass function since cosmic variance would be underestimated along the radial axis.

To infer the stellar masses of the HSC-SSP sample, we use the spectral energy distribution (SED) fitting code LEPHARE (Arnouts et al. 1999; Ilbert et al. 2006; Arnouts & Ilbert 2011) with the Bruzual & Charlot (2003) (BC03 here and after) templates to estimate galaxy stellar masses from the g, r, i, and z cModel magnitudes. We then use the luminosity function tool ALF (Ilbert et al. 2005) to construct galaxy stellar mass functions for each pointing using the method of Sandage et al. (1979). Galaxies are selected in the r band with an apparent magnitude cut of 26. We first constrain the knee of the mass function (M_s) by computing the mass function for each pointing in a larger redshift slice (0.1 < z < 0.4) before re-fitting the mass function with M_s fixed for the smaller volume. For the simulated sample we follow a similar procedure, first obtaining dust-attenuated g, r, i, and z magnitudes for galaxies identified with HOP using SUNSET (see Martin et al. 2021, Sect. 2.2.1). To approximate the selection effects present in real data, we select galaxies by their effective surface brightness, where the probability of selecting a galaxy is proportional to the surface brightness completeness of the HSC survey; this value is estimated by assuming that the true number of objects continues to rise exponentially as a function of effective surface brightness after the turnover in the number of galaxies observed by HSC. We again use LEPHARE and ALF to construct the galaxy stellar mass function using the same 26 mag cut in the r band. Because of the more limited volume of NEWHORIZON, the number of galaxies that are considerably more massive than the knee of the mass function is too small to effectively constrain this value, thereby leading to unrealistic fits. We therefore fix M_s at a value of 10^10.8 M_⊙, which is calculated from the full volume of the HORIZON-AGN simulation. While varying M_s also necessarily affects the slope at the low-mass end, this is not significant enough to qualitatively alter our comparison to the observed mass functions within a reasonable range of values (e.g., 10^10.6 M_⊙–10¹¹ M_⊙).

The galaxy mass function, which is a volume-integrated quantity poses a conceptual challenge to a zoom-in simulation. Indeed, galaxies within halos polluted with low-resolution DM particles continue to form stars, and it is questionable whether or not their contribution to the overall cosmic star formation should be taken into account. In addition, we have to determine the actual corresponding volume of the zoom-in region, which can expand or contract over time. For the volume entering the calculation of the galaxy mass function (and in other volume-integrated quantities measured in this work), we take the entire initial volume of the zoom-in region of the simulation, hence, (16 Mpc)³. We could alternatively use the sum of each individual leaf cell that passes a given threshold value of the passive scalar colour value (see Sect. 2.1). The corresponding initial volume can be reduced by 20–40% for a threshold value of resp. 0.1–0.9, depending on redshift. We decided to simplify the problem by taking the initial zoom-in volume, but we note that the presented volume-integrated quantities are only a lower limit and can be a few tens of per cent higher.

Figure 6 shows the galaxy stellar mass function from NEWHORIZON, HSC-SSP, and from the literature (Sedgwick et al. 2019a; Davidzon et al. 2017; Tomczak et al. 2016; Song et al. 2016; Grazian et al. 2015; D’Souza et al. 2015; Fontana et al. 2014; Bernardi et al. 2013; Ilbert et al. 2013; Baldry et al. 2012; Bielby et al. 2012; González et al. 2011; Pozzetti et al. 2007). We note that Sedgwick et al. (2019a) includes only star-forming galaxies. The light blue squares with error bars represent the NEWHORIZON stellar mass function with Poisson errors for all galaxies. The purple circles show the same, but include only galaxies whose halos are not contaminated by low-resolution particles from outside of the highest resolution zoom region – a simple correction is made to account for the smaller effective volume by dividing the mass function by the fraction of uncontaminated galaxies. Additionally, the mass function (with selection effects) for NEWHORIZON that is constructed using the Sandage et al. (1979) method (STY) is shown as a thick purple dotted line. The black line indicates the median galaxy stellar mass function from the 100 random pointings from the HSC-SSP deep layer. Various other mass functions from the literature are also indicated as thin coloured lines in each panel.

Fig. 6.

Galaxy stellar mass function at z = 0.25 in NEWHORIZON and the HSC-SSP survey. The light blue squares and dark blue circles with error bars indicate the NEWHORIZON stellar mass function for all galaxies (offset by 0.15 dex) and only uncontaminated galaxies (with a volume correction applied). The black line indicates the median galaxy stellar mass function from 100 random pointings from a volume of the HSC-SSP deep layer with the same volume as NEWHORIZON and the black error bar indicates the 90% scatter in the mass function at the low-mass end. A comparison is made with additional mass functions from literature (Sedgwick et al. 2019a; Davidzon et al. 2017; Tomczak et al. 2016; Song et al. 2016; Grazian et al. 2015; D’Souza et al. 2015; Fontana et al. 2014; Bernardi et al. 2013; Ilbert et al. 2013; Baldry et al. 2012; Bielby et al. 2012; González et al. 2011; Pozzetti et al. 2007), which are shown as thin coloured lines. The thick purple dotted line indicates the NEWHORIZON mass function constructed using the Sandage et al. (1979) method (STY) and including selection effects.

Once selection effects are included, the NEWHORIZON mass function lies within the upper range of the observational mass functions shown. The discrepancy between the raw mass function likely emerges from incompleteness in the observed data at low surface brightness, meaning the observed mass function may be underestimated towards lower mass. The effect of selection effects and environment on the galaxy stellar mass function will be explored in more detail in Noakes-Kettel et al. (in prep.).

The 90% variance in the low-mass end of the HSC mass function is indicated by a black error bar. Over such a limited volume, the normalisation of the galaxy stellar mass function varies significantly. We note that our estimate of the variance may be an underestimate as redshift uncertainties are significantly larger than the selected volume. Additionally the location of the four HSC-SSP deep fields were chosen to enable certain science goals and not necessarily to sample representative volumes of the Universe as a whole.

3.3. Surface brightness-to-galaxy mass relation

While many comparisons between theory and observation treat galaxies as one-dimensional points, the two-dimensional distribution of baryons (e.g., as summarised by the effective surface brightness and effective radius) are important points of comparison. Since high-resolution cosmological simulations now offer predictions of the distribution of baryons, it is worth comparing these predictions to observational data of galaxy surface brightnesses.

Here, we obtain the dust attenuated surface brightness for each NEWHORIZON galaxy using the intensity-weighted second order central moment of the stellar particle distribution (e.g., Bernstein & Jarvis 2002) and we compare these values with Sedgwick et al. (2019a).

For each star particle, we first obtain the full SED from a grid of dust attenuated BC03 simple stellar population models corresponding to the age and metallicity of the star particle. We redshift each BC03 template to match the overall redshift distribution of the observed sample to account for surface-brightness dimming. Since low-mass galaxies in this sample are biased towards lower redshift, we also account for this by drawing a redshift for each galaxy from the conditional probability distribution P_z(M_s) (i.e., the redshift probability distribution at a given stellar mass as found in Sedgwick et al. 2019a) so that they are redshifted and their flux and size are calculated according to a redshift whose probability is related to the stellar mass of the object. The SEDs are then convolved with the response curve for the SDSS r-band filter. We then weight by the particle mass to obtain the luminosity contribution of each star particle, and obtain the apparent r-band magnitude by converting the flux to a magnitude and adding the distance modulus and zero point.

The second moment ellipse is obtained by firstly constructing the covariance matrix of the intensity-weighted second order central moments (sometimes called the moment of inertia) for all the star particles, that is

$\begin{matrix} cov [I (x, y)] = [\begin{matrix} I x^{2} & Ixy \\ I x y & I y^{2} \end{matrix}], \end{matrix}$ $\begin{aligned} \mathrm{cov} [I(x,y)]=\begin{bmatrix} Ix^{2}&Ixy\\ Ixy&Iy^{2}\\ \end{bmatrix}, \end{aligned}$ (8)

where I is the flux, and x and y are the projected positions from the barycentre in arcseconds of each star particle in the galaxy. The major ( $α = \sqrt{λ_{1} / Σ I}$ $\alpha=\sqrt{\lambda_{1}/\Sigma I}$ ) and minor ( $β = \sqrt{λ_{2} / Σ I}$ $\beta=\sqrt{\lambda_{2}/\Sigma I}$ ) axes of the ellipse are obtained from the covariance matrix, where λ₁ and λ₂ are its eigenvalues and ΣI is the total flux. The scaling factor, R, scales the ellipse so that it contains half the total flux of the object. Finally, the mean surface brightness within the effective radius can be calculated from ⟨μ⟩_e, r = m − 2.5log₁₀(2)+2.5log₁₀(A), where A = R²αβπ and m is the r-band apparent magnitude of the object. We repeat this process in multiple orientations (xy, xz and yz), taking the mean surface brightness for each object. The method presented in this work is equivalent to the method employed in SEXTRACTOR (Bertin & Arnouts 1996) to derive basic shape parameters, which Sedgwick et al. (2019a) use to derive their measurements.

In Fig. 7, we show the evolution of the surface brightness ⟨μ⟩_e, r vs. stellar mass plane in NEWHORIZON. In the bottom panel we compare the predicted surface brightnesses to a recent work that uses the IAC Stripe 82 Legacy Survey project (Sedgwick et al. 2019a,b). This study is one of few that probes the surface brightnesses of galaxies down into the dwarf regime, which is only possible at low redshift using past and current surveys, which are typically very shallow. To probe galaxy surface brightness down to faint galaxies Sedgwick et al. introduce a novel technique with core-collapse SNe (CCSNe). Using custom settings in SExtractor (Bertin & Arnouts 1996) they extract the host galaxies of these CCSNe, including those that are not detected in the IAC Stripe 82 Legacy survey. The resultant sample is free of incompleteness in surface brightness in the stellar mass range M_s > 10⁸ M_⊙; a host is identified for all 707 CCSNe candidates at z < 0.2. Given the high completeness of the sample at low surface brightness and the relative ease with which we can model the selection function and apply it to our simulated data, this dataset is an ideal choice to compare to the NEWHORIZON data. More details on how the matching between the two datasets has been completed are available in Jackson et al. (2021a).

Fig. 7.

Surface brightness vs. stellar mass in the NEWHORIZON simulation for 3 redshifts. The grey points indicate the entire galaxy population of NEWHORIZON in all 3 panels with the median lines for the redshift in question and z = 0.25 shown in blue and orange, respectively. In the bottom panel the open blue points indicate galaxies from Sedgwick et al. (2019a) and black points are NEWHORIZON galaxies. The predicted surface brightness vs. stellar mass plane in NEWHORIZON corresponds well to that where the Sedgwick et al. galaxies are complete. The red dotted lines in the bottom panel indicate the 70% and 10% completeness limits from the SDSS (see e.g., Table 1 in Blanton et al. 2005). The overwhelming majority of galaxies in the Universe lie below the surface brightness thresholds of surveys such as the SDSS; only those galaxies that depart strongly from the typical surface brightness vs. stellar mass relation are likely to be detectable in these datasets.

Figure 7 shows that the surface brightness vs. stellar mass plane in NEWHORIZON corresponds well to Sedgwick et al. (2019a), where the observational data is complete; we note that the simulation is not calibrated to reproduce galaxy surface brightness. The flattening seen in the observations is due to high levels of incompleteness at M_s < 10⁸ M_⊙. The prediction for the evolution of this plane to higher redshifts shows that NEWHORIZON galaxies have increasing brightness at higher redshift for a fixed stellar mass (i.e., galaxies are more concentrated, see Sect. 3.7) can be tested using data from future instruments such as the LSST.

3.4. Star formation rates

Figure 8 shows the cosmic SFR density as a function of redshift in NEWHORIZON compared to observations. The cosmic SFR density is obtained by summing all star particles formed at a given time over the last 100 Myr within the entire volume of the simulation, which are associated with a galaxy (sub)structure, whether pure or not⁴. Since stellar particles lose mass owing to stellar feedback, we compensate for this mass loss when reconstructing the SFR. The observational data in Fig. 8 (Behroozi et al. 2013; Madau & Dickinson 2014; Novak et al. 2017; Driver et al. 2018) are scaled to a Chabrier IMF whenever necessary (i.e., cosmic SFR is decreased by a factor 1.6 when going from a Salpeter to Chabrier IMF). Several datasets were selected to illustrate the typical variation from inter-publication variance and IMF assumptions (see e.g., Behroozi et al. 2013; Madau & Dickinson 2014, for a discussion). The obtained NEWHORIZON cosmic SFR density is slightly above the observational values collected by Behroozi et al. (2013). For the NEWHORIZON SFR values, the maximum offset is at the lowest redshift, although the numerical sampling is worse in this case and concentrated over a few rather massive M_s ≥ 10¹⁰ M_⊙ objects. To estimate the effect of cosmic variance, we relied on the measurement of the cosmic variance on the SFR density in Illustris (∼0.2 dex for a 35 Mpc box length; Genel et al. 2014) rescaled to the corresponding smaller volume here (the large error bar on the left-hand side of the top panel of Fig. 8). This shows that, with this small simulated volume, the cosmic variance is the largest source of uncertainty on the simulated cosmic SFR density compared to observational datasets.

Fig. 8.

Cosmic SFR density (top panel) and stellar density (bottom panel) as a function of redshift in NEWHORIZON (solid black or red lines) compared to observations (Behroozi et al. 2013, B13; Madau & Dickinson 2014; Novak et al. 2017; Driver et al. 2018) as indicated in the panels. All observational quantities are shown for a Chabrier (2005) IMF except for the dashed blue line, which indicates the Madau & Dickinson (2014) fit for a Salpeter IMF as originally assumed in their analysis. For the stellar density, the reconstructed result from the SFR density for NEWHORIZON (magenta line) and for the fit from B13 (grey lines) using two different return fractions R are also shown. The red coloured lines indicate the cosmic SFR densities of the different galaxy stellar mass bins as indicated in the panel in log M_⊙ units. The large error bar in black in the SFR density panel corresponds to the estimate of the cosmic variance (see text for details). The cosmic SFR in NEWHORIZON shows the qualitative expected trend over time; however, there is a systematic offset by a factor 1.5 − 2 with respect to the observational data assuming standard Chabrier-like IMF. The cosmic stellar density in NEWHORIZON is broad agreement with the data with a slight overestimate at low redshift, although there is an important uncertainty on the simulated cosmic SFR (and thus stellar density) due to the large cosmic variance associated with the simulated volume.

We also show (in Fig. 8) the cosmic stellar density; this value is obtained by summing over the individual mass of all the star particles in the simulation, which are again associated with a galaxy (sub)structure, whether pure or not. Comparing the NEWHORIZON cosmic stellar density to that directly measured in Driver et al. (2018) produces a factor 2.5 difference at z = 0.25, whose cosmic SFR density lies on the low side of aggregated observational values; this agrees with the mismatch that is also observed in the cosmic SFR density with the same observations. The reconstructed cosmic stellar density is obtained from the time-integrated cosmic SFR density with an instantaneous stellar mass return of R = 0.31, which is the value used in the simulation for SN feedback. This reconstructed stellar density is in excellent agreement with the direct measurement It has to be noted that the two direct measurements of Driver et al. (2018) are self-consistent for R = 0.5 (see their Fig. 16). We also show the reconstructed cosmic stellar density from the best fit to the cosmic SFR density from the Behroozi et al. (2013) data (Fig. 8), for two values of the return rate: R = 0.31 and R = 0.5. The NEWHORIZON predicted cosmic SFR density and stellar density⁵ compared to Behroozi et al. (2013) are within a factor of 2 (for R = 0.31); thus these values are in reasonable agreement with this set of data.

The contribution to the cosmic SFR and stellar density is further subdivided into separate galaxy stellar mass bins as indicated in the panels of Fig. 8. Low-mass galaxies M_s < 10⁹ M_⊙ dominate the SFR and stellar mass budget in the early Universe, while intermediate mass galaxies 10⁹ ≤ M_s/M_⊙ < 10¹¹ take over below the peak epoch of star formation (typically at z ≃ 1.5); Milky Way-like galaxies dominate the cosmic SFR by the lowest redshift z = 0.25, which is in qualitative agreement with previous theoretical predictions (e.g., Béthermin et al. 2013; Vogelsberger et al. 2014; Genel et al. 2014). Although only a few galaxies contribute to this range of mass, the highest mass bin M_s seems to marginally contribute to the cosmic SFR at low redshift⁶; however, it represents nearly half of the cosmic stellar density, thereby highlighting the role played by satellite infall (stars formed ex situ) in the assembly of massive galaxies (De Lucia & Blaizot 2007; Oser et al. 2010; Dubois et al. 2013, 2016).

The specific SFR (sSFR) of individual galaxies can be computed by measuring the stars younger than 100 Myr within their effective radius R_eff (see Sect. 3.7 for the calculation of R_eff) and dividing by the current stellar mass within R_eff at the given redshift. Figure 9 shows the resulting mean sSFR as a function of galaxy stellar mass for different redshifts⁷. Galaxies are usually separated into a main sequence of star-forming (active) galaxies and quenched galaxies that are passively evolving, showing a significant evolution over mass and redshift of their their respective fraction (see e.g., Kaviraj et al. 2007; Muzzin et al. 2013; Wetzel et al. 2013; Furlong et al. 2015; Fossati et al. 2017); in particular, the bulk of the quenched galaxies found in central galaxies are more massive than a few 10¹⁰ M_⊙ or in satellites hosted by massive groups or clusters. Only active galaxies are selected based on their level of sSFR, that is with a sSFR above 0.01 Gyr⁻¹. Including inactive galaxies slightly changes (by up to 20%) the mean sSFR at z = 0 and z = 1 for galaxies with stellar mass M_s ≤ 10⁹ M_⊙. The exact criterion for the separation between the two population is somewhat arbitrary, but the idea is that quenched galaxies should clearly stand out from the main population of galaxies (see e.g., Donnari et al. 2020, for how distinguishing between active and passive galaxies affects the quenched fraction). The sSFR at z = 4 and 2 show no trend with stellar mass, while the low redshift relations at z = 1 and 0.25 show a significant decrease (quenching) at M_s > 10¹⁰ M_⊙. These simulated values are compared to the best-fit relation from Behroozi et al. (2013) obtained from a collection of observational data (see references therein). The simulation agrees fairly well with the data at z = 4, 2 and 0.25, but the values at z = 1 are significantly below those of the data; however, this corresponds to the decrease in the cosmic SFR density right before a peak that almost doubles the overall SFR in galaxies more massive than M_s = 10⁹ M_⊙. Nonetheless, it has to be noted that there is a systematic offset between the NEWHORIZON sSFR and the data. Simulated sSFR are on average systematically lower than the data. This leads to a slight inconsistency with the cosmic SFR density, which is higher than the values found in the data. Similar tensions are noted within the data themselves (see Appendix C.4 of Behroozi et al. 2019).

Fig. 9.

sSFR as a function of galaxy stellar mass at different redshifts (as indicated in the panel) in NEWHORIZON with solid lines. The error bars stand for the error around the mean. The symbols correspond to the best fit from Behroozi et al. (2013) from their collection of observational data; their uncertainty is shown in dashed lines. The simulated sSFR show very little evolution with stellar mass at the two highest redshifts, while there is a significant quenching at the lowest redshifts for stellar masses above M_s > 10¹⁰ M_⊙. Simulated sSFR show a fair level of agreement with the observations.

3.5. Kennicutt–Schmidt relation

Figure 10 shows the Kennicutt–Schmidt relation of surface density of SFR (Σ_SFR) as a function of surface density of total (HI+H₂) gas (Σ_gas) for galaxies in the NEWHORIZON simulation at redshifts z = 4, 2, 1, 0.25 and with M_s > 10⁷ M_⊙, compared to observations. The quantities Σ_SFR and Σ_gas are computed within R_eff. The SFR is obtained by summing over all stellar particles with stellar age below 10 Myr and the HI+H₂ gas is selected as gas with density n > 0.1 cm⁻³ and temperature T < 2 × 10⁴ K. The observational data shown in Fig. 10 include z ∼ 1.5 BzK-selected normal galaxies (Daddi et al. 2010a) and z = 1 − 2.3 normal galaxies (Tacconi et al. 2010), high-z submillimetre selected galaxies (SMGs; Bouché et al. 2007; Bothwell et al. 2009), IR-luminous galaxies (ULIRGs), and spiral galaxies are taken from the sample of Kennicutt (1998) as compiled in Daddi et al. (2010b) with a consistent choice of the conversion factor used in Fig. 10 to derive molecular gas masses from CO luminosities (α_CO) and a Chabrier (2003) IMF.

Fig. 10.

Kennicutt–Schmidt relation compared to observations. Filled coloured circles correspond to NEWHORIZON galaxies at various redshifts. The mean and error on the mean are shown with solid and dotted coloured lines, respectively. The empty black symbols correspond to observational datasets: triangles are local spirals (Kennicutt 1998), stars are high-z BzK/Normal galaxies (Daddi et al. 2010a; Tacconi et al. 2010), circles high-z mergers (Bouché et al. 2007; Bothwell et al. 2009), and squares low-z mergers (Kennicutt 1998). The solid black and dotted black lines represent the sequence of discs and starbursts, respectively, from Daddi et al. (2010b). The NEWHORIZON tool is able to capture the galaxy main sequence, although the SFR surface densities are slightly lower than in the observational data. There is a larger fraction of high SFR galaxies at higher redshift in qualitative agreement with observations.

With decreasing redshift, the population of simulated galaxies as a whole moves roughly along the sequence of discs (solid line, Daddi et al. 2010b) towards lower values of Σ_gas and Σ_SFR. Qualitatively, simulated galaxies occupy comparable regions of the Kennicutt–Schmidt parameter space, reproducing the observed diversity of star-forming galaxies; however there are some notable differences. At z ∼ 4, there are galaxies at low Σ_gas < 10 M_⊙ pc⁻² and high Σ_SFR > 10⁻² M_⊙ yr⁻¹ kpc⁻², which seem to be offset from the bulk of the population. The reason for this offset is their low gas fraction that is ≲0.1 regardless of their stellar mass and SFR. These galaxies cover the entire stellar mass range, having similar SFR and size to galaxies with comparable Σ_SFR and higher Σ_gas. By redshift z ∼ 3, there are no galaxies left in this region of the parameter space. Below z = 3, the number of galaxies on the canonical sequence of starbursts (Daddi et al. 2010b) decreases with decreasing redshift. The slope of the average Σ_gas − Σ_SFR relation does not evolve strongly between z = 1 and 3; however, the slope is offset from the sequence of discs by ∼0.5 dex. At z = 0.25, the lowest available redshift, the slope is in qualitative agreement with observations of local spirals, albeit it still has an offset. We have checked that when considering star-forming gas only, that is gas with density n > 10 cm⁻³ and temperature T < 2 × 10⁴ K, the average Σ_gas − Σ_SFR relation at z < 3 follows the sequence of discs (Kraljic et al., in prep.).

3.6. Galaxy-to-halo mass relation

Figure 11 shows the stellar mass of galaxies as a function of their host DM halo mass at redshifts z = 4, 2, 1, 0.25 compared to the semi-empirical relation from Moster et al. (2013) and Behroozi et al. (2013). The mass of the DM halo is obtained by taking the total (DM, gas, stars, and MBHs) mass enclosed within the radius corresponding to Δ_c times the critical density at the corresponding redshift, where the Δ_c of the analytical form is taken from Bryan & Norman (1998). This procedure is similar to that of Behroozi et al. (2013), but differs from Moster et al. (2013), where Δ_c is fixed at 200. We note that this is a ∼0.1 − 0.2 dex increase in halo mass compared to the value of the virial mass obtained by the AdaptaHOP halo mass decomposition. AdaptaHOP galaxies are considered in this work; hence, satellites (and in situ stellar clumps) are not considered when the total stellar mass is measured, but these clumps only constitute a small fraction of the total galaxy mass as discussed in Sect. 3.10 (typically 10% and 1% of the total stellar mass at the most extreme redshifts, resp. z = 4 and z = 0.25). Satellite galaxies are connected to their host subhalo mass at the current redshift; thus these galaxies should not be compared directly with the semi-empirical constraints, whereby they reconstruct the relation with the halo mass before it becomes a subhalo. Satellite galaxies (not shown here) have systematically a larger value with respect to their host subhalo mass because the DM particles are first stripped by gravitational interactions with the main halo well before the more concentrated stellar mass becomes affected (Peñarrubia et al. 2008; Smith et al. 2016).

Fig. 11.

Stellar-to-halo mass relation (top) and baryon conversion efficiency M_s/M_h at different redshifts (as indicated in the panels). The solid red lines represent the average with their error of the mean (dashed) with individual points as plus symbols. The solid black line indicates the semi-empirical relation from Moster et al. (2013) and the green line from Behroozi et al. (2013) at the indicated redshift and dot-dashed lines the respective relations at z = 0. The cyan lines stand for the constant star formation conversion efficiencies. The dotted lines correspond to the minimum stellar mass detected by the galaxy finder. The simulated relation between central stellar mass and halo mass is in fair qualitative agreement with the semi-empirical relations, where the increase in the baryon conversion efficiency with mass up to the peak at Milky Way halo mass is captured; however this efficiency is significantly overestimated below that peak.

There is a fairly good qualitative agreement of the stellar-to-halo mass relation with the general trends from Moster et al. (2013) and Behroozi et al. (2013) at all redshifts for the population of main halos. The baryon conversion efficiency, that is the ratio of M_s/M_h, shows a maximum at near to the Milky Way mass at a few 10¹¹ M_⊙, although this value is slightly below the expected peak of the semi-empirical relations. This ratio steeply decreases with the decreasing halo mass and the ratio plateaus around the Milky Way scale, while it is expected to decrease above this mass. However, see Kravtsov et al. (2018) for the underestimated stellar light component at those group and cluster scales together with IMF variation effects. The simulated stellar masses are still significantly above the relation; however there is a better agreement at the higher masses M_h > 10¹¹ M_⊙.

3.7. Size-to-galaxy mass relation

Galaxy effective radii are obtained by taking the geometric mean of the half-mass radius of the projected stellar densities along each of the Cartesian axis. For this measurement, we consider AdaptaHOP galaxies since HOP galaxies can largely overestimate the effective radius of galaxies at the high-mass end (see Appendix B) when satellites orbit around centrals and are connected by the diffuse stellar light. At the same time, star-forming clumps are also removed but as they only represent a small fraction of the total stellar mass, they do not have a significant impact on the determination of the effective radius. We note that this procedure tends to reduce the scatter of the relation, but in this work we are mostly interested in investigating to what extent the observed mean relation is reproduced.

Figure 12 shows the effective radius of NEWHORIZON galaxies for four different redshifts with simulated data points in grey plus signs and its average and error around the mean with black lines, compared to the observational relation obtained by Mowla et al. (2019) in purple at the corresponding redshift, which we have linearly interpolated from the two contiguous redshifts provided in Mowla et al. (2019). To guide the eye, and because there is no observational data available beyond z = 2.69, we also overplotted the observational relation at the lowest (z = 0.37 in blue) and highest redshift (z = 2.69 in red). There is overall good agreement between the simulated size-mass relation and the observations at all redshifts. We note that galaxy sizes around 1 kpc to a few kiloparsec are described with several resolution elements (resolution is 34 pc), that is more than what is typically achieved in the low-mass range for large-scale simulations such as Horizon-AGN (Dubois et al. 2014a), EAGLE (Schaye et al. 2015), or IllustrisTNG (Pillepich et al. 2018). Pillepich et al. (2018) show that with the IllustrisTNG model, the mean size obtained at their low-mass end of around a few kiloparsec is systematically scaled to 3–4 times their minimum stellar spatial resolution. Thus, the simulated low-mass galaxies (though above a few 10⁷ M_⊙) are comfortably resolved in NEWHORIZON in comparison. The NEWHORIZON tool recovers the size growth of galaxies with mass and the size growth of galaxies with redshift (at a given mass): as they grow in mass, galaxies tend to be more extended, and the size-mass relation produces more extended galaxies with time as measured in observations (e.g., Trujillo et al. 2006; van der Wel et al. 2014; Mowla et al. 2019).

Fig. 12.

Effective radius as a function of stellar mass in NEWHORIZON at different redshifts for individual galaxies (grey crosses) and for the average (solid line) with standard deviation (dashed lines). The sizes are computed with the geometric mean of the x, y, z effective radius using the AdaptaHOP classification. Observations from Mowla et al. (2019) are shown as a purple line at the current redshift when available and with their most extreme redshift fits at z = 0.37 and z = 2.69 in blue and red, respectively, which are extrapolated (triple dot-dashed line) beyond their range of available data (solid line) to guide the eye. Galaxies in NEWHORIZON are more compact at high redshift as opposed to low redshift in good agreement with observations.

The large spread of simulated data points below M_s < 10⁸ M_⊙ corresponds to embedded in situ stellar clumps for the most compact ones. Outliers with large sizes at relatively low stellar masses are due to satellite galaxies embedded in diffuse stellar light of a much more massive companion, where the galaxy finder has failed to extract them from the background correctly. There is formation of compact massive galaxies at high-mass end (M_s > 10¹⁰ M_⊙) and more compact galaxies at higher redshift, with sizes below R_eff ≲ 1 kpc. This feature is reminiscent of the nugget formation of massive galaxies that have endured a gas-rich compaction event triggering high levels of SFR before quenching (see e.g., Dekel & Burkert 2014; Zolotov et al. 2015; Lapiner et al. 2021). The NEWHORIZON galaxies however seem to fail to reproduce the rapid rise in galaxy sizes at the high-mass end. This might partially be a consequence of the volume that is simulated corresponds to a region of average density, and because we lack the formation of dense environments that are the main producers of such massive objects, where galaxies are more likely to be more passive at low redshift and, hence, built though mergers (Martin et al. 2018b) that lead to a large size increase (e.g., Naab et al. 2009; Oser et al. 2010; Dubois et al. 2013, 2016).

3.8. Stellar and gas metallicities

The mass-weighted stellar metallicity Z_s is computed for all the stars within R_eff for each galaxies. The value is renormalised to the solar metallicity Z_⊙ = 0.01345 (Asplund et al. 2009). Figure 13 shows the NEWHORIZON stellar metallicities as a function of each galaxy stellar mass at different redshifts z = 4, 2, 1, 0.25 and is compared to the (luminosity-weighted) observations by Gallazzi et al. (2014) at z = 0.7 and z = 0.1 and to the (mass-weighted) observations by Panter et al. (2008) at z = 0.1⁸. We note that observational estimates of metallicities barely differ whether are they mass-weighted or luminosity-weighted, as also seen in simulations (e.g., as is reported in De Rossi et al. 2017, for the EAGLE simulation). The NEWHORIZON galaxies show an increase in stellar metallicity with mass and with decreasing redshift as expected from the continuous release of metals from stars and their increased level of retention in more massive galaxies (e.g., Tremonti et al. 2004). Despite the extremely crude modelling of metal release used in NEWHORIZON (all metals are released at once after 5 Myr), the mass-metallicity relation at z = 0.25 is consistent with observations at z = 0.1 − 0.7. However, the mass-metallicity relation in NEWHORIZON shows a slower evolution over redshift than what is suggested by observations. We also show in Fig. 13 the metallicity of the cold gas phase, namely gas with density n > 0.1 cm⁻³ and temperature T < 2 × 10⁴ K. Gas metallicity is systematically larger by ≃50% than the stellar metallicity at any given redshift as an effect of the stellar metallicity being composed of very old poorly enriched stars.

Fig. 13.

Average stellar (coloured solid) and gas (coloured dashed) metallicity as a function of stellar mass for different redshifts in NEWHORIZON as indicated in the panel. The error bars stand for the error around the mean. Observational fits for the stellar metallicity from Gallazzi et al. (2014; solid black at z = 0.1, and dashed black at z = 0.7) and from Panter et al. (2008; solid grey) are also shown. The metallicity of both the gas and stellar component are decreasing with increasing redshift owing to less chemically enriched galaxies.

The gas oxygen abundance ratio and its relation to galaxy mass (or the so-called mass–metallicity relation; MZR) exhibits an evolution with redshift (e.g., Erb et al. 2006; Maiolino et al. 2008; Zahid et al. 2011, 2014; Yabe et al. 2015; Sanders et al. 2018), which is reflective of the redshift evolution of the SFR–M_s that drives the scatter of the MZR (Mannucci et al. 2010). To qualitatively appreciate the evolution over redshift of the MZR, we show the data from Zahid et al. (2014) measured within fixed aperture of 10 kpc at z = 1.6 together with the NEWHORIZON data in Fig. 14. We also report the various fits of the MZR obtained from the spatially resolved data of the SAMI galaxy survey within R_eff for different oxygen abundance calibrators from emission lines, which is known to be a major source of uncertainty (see details in Sánchez et al. 2019), and the fitting relation from Curti et al. (2020) for the MZR of SDSS galaxies at z = 0.027. For NEWHORIZON, even though the tool does not self-consistently following the amount of oxygen (mainly produced by massive stars), we rescale the gas metallicity values by a factor corresponding to the fractional abundance of oxygen in the solar atmosphere and further assume that the fraction of hydrogen is solar (Asplund et al. 2009): the mass fraction of hydrogen of the gas is 73.4%, and oxygen represents 43% of the total mass of elements heavier than H and He. This is, obviously, a crude estimate of the actual oxygen abundance in the simulation since the fractional amount of oxygen amongst metals varies with the age of the galaxy, and in turn with metallicity (or galaxy mass), and for instance its [O/Fe] ratio is known to increase faster with decreasing metallicity than some other significantly abundant elements such as carbon or nitrogen (see e.g., Prantzos et al. 2018, for a compilation of observational data). Figure 14 shows that the increase in the gas oxygen abundance with galaxy mass and time is well captured by the simulation despite the simple scaling for converting the metallicity in NEWHORIZON into oxygen, and the choice of different apertures amongst data and our simulated galaxies. However, the highest redshift bin NEWHORIZONz = 4 seems to have a serious offset (∼0.6 dex) with respect to observations at z = 3.5 (Mannucci et al. 2009), even though there is a large spread and large uncertainties in the observational data at this redshift. This range of redshift where observed metallicities are a factor 10 below local values calls for a more accurate treatment of stellar yield release (see e.g., Mannucci et al. 2009; Prantzos et al. 2018).

Fig. 14.

Average abundance of oxygen (coloured solid) in the cold gas as a function of stellar mass for different redshifts in NEWHORIZON as indicated in the panel. The error bars stand for the error around the mean. Several observational relations are also shown from Mannucci et al. (2009; M09; in dotted lines with four different colours corresponding to z = 0.07, 0.7, 2.2, and 3.5), Zahid et al. (2014; Z14; dot-dashed line), Sánchez et al. (2019; S19; dashed lines with different lines corresponding to various oxygen calibrators), and Curti et al. (2020; C20; black solid).

3.9. Stellar kinematics

Stellar kinematics are obtained by first computing the angular momentum vector of the stars around the centre of the galaxy. This vector is then used to decompose the kinematics into a cylindrical frame of reference. The stellar rotation V is the average of the tangential component of velocities, while the 1D stellar dispersion σ is the dispersion obtained from the dispersion around each mean component, that is $σ^{2} = (σ_{r}^{2} + σ_{t}^{2} + σ_{z}^{2}) / 3$ $\sigma^2=(\sigma_{r}^2+\sigma_{t}^2+\sigma_{z}^2)/3$ . The kinematics are computed from the AdaptaHOP extracted stars within two different radii R_eff or 2R_eff to exemplify the effect of aperture on the measured kinematics.

Figure 15 shows the rotation/dispersion ratio for stars within different radii and various redshifts. Independently of radius, the ratio increases with decreasing redshift and with stellar mass (except at high redshift z = 4). The stellar component is more rotationally supported over time. We note that these relations have a significant scatter, illustrated by the distribution of individual points for z = 0.25 in Fig. 15, of around 0.3. Galaxies also exhibit a stronger rotational support with respect to dispersion when more distant stars are taken into account. Galaxy interiors are more supported by dispersion, while the outskirts are more rotationally dominated when V/σ is measured either in R_eff or 2R_eff. This is observationally confirmed (e.g., Emsellem et al. 2011; Naab et al. 2014; van de Sande et al. 2017) and expected since central regions of galaxies are probing the bulge component mostly supported by dispersion, while the outskirts of the galaxy have a significant rotating disc components in cases where the galaxy is strongly discy. Nonetheless, elliptical galaxies could eventually show a reverse trend, where central regions have a significantly large amount of rotation with kinematically decoupled cores (Krajnović et al. 2013, 2015; Coccato et al. 2015). For example these regions are rejuvenated by a recent episode of star formation fed by counter-rotating filaments (Algorry et al. 2014) or mergers (Bois et al. 2011; Moody et al. 2014), while the large stellar halo of the elliptical is more likely to be dispersion-dominated.

Fig. 15.

Ratio of stellar rotation over dispersion as a function of galaxy stellar mass for different redshifts (colours) as indicated in the panels. The kinematics are measured within R_eff (top) or 2R_eff (bottom). The solid lines indicate the average and the dashed lines are the error around the mean, with individual points shown at z = 0.25 only (purple plus signs) to appreciate the scatter of the distribution. Galaxies show an increase support of stellar rotation over dispersion with time and galaxy mass (except for z = 4) with stars outside R_eff having more rotation than inside.

Alternatively it is possible to compute the fraction of dispersion-supported, elliptical, or irregular galaxies f_{ell + irr} = 1 − f_disc, by positing that a galaxy is elliptical or irregular when V/σ < 0.5 (and conversely a disc when V/σ ≥ 0.5), where the exact threshold value is arbitrary and where in this case this value is chosen to best fit the observational data. It has to be noted that the classification with kinematics does not allow us to distinguish between irregulars and ellipticals (both sum up to f_{ell + irr}) but this is done using the asymmetry index (see below). This fraction of elliptical/irregular galaxies can be compared qualitatively with observational morphological (visual) classifications from Conselice (2006). The exact value of the fraction of each morphological type with mass depends on the adopted threshold in V/σ, nonetheless the obtained trends with mass are robust against reasonable variations in the adopted thresholds as show in Appendix C. Figure 16 shows the fraction of elliptical/irregular galaxies as a function of galaxy stellar mass at the lowest redshift of NEWHORIZON, that is z = 0.25 compared to the observations at z = 0 for kinematics measured in different radii. There is a fair agreement of the simulated data at R_eff with observations with a similar stellar mass trend. The agreement is weaker for 2R_eff, but the level of agreement also depends on the arbitrary threshold value on V/σ adopted (here 0.5). There are fewer elliptical/irregular galaxies as galaxies increase in mass up to the maximum galaxy masses probed here: M_s ≃ 10¹¹ M_⊙. For the most massive galaxies, owing to the limited volume of the simulation and the lack of groups of galaxies, it is impossible to conclude whether our obtained fraction of ellipticals is consistent with observations.

Fig. 16.

Fraction of galaxies classified by their morphological type: f_{ell + irr} for elliptical+irregulars (red solid line; galaxies with V/σ < 0.5) in NEWHORIZON based on the stellar kinematics measured within R_eff (top) or 2R_eff (bottom) and as a function of stellar mass at z = 0.25; conversely, the fraction of discs f_disc = 1 − f_{ell + irr} is compared to data from Conselice (2006), which is at z ≃ 0, and is represented by symbols with error bars. The fraction of irregulars in NEWHORIZON inferred through the asymmetry index, within the same aperture in both panels, is also shown as the magenta curve. The error bars and dashed line stand for the error on the mean. The trends of the morphological fractions with mass in NEWHORIZON is qualitatively consistent with observational data.

We further classify NEWHORIZON galaxies as irregulars, using the asymmetry index A_r from Conselice et al. (2000) on the rest frame r-band extracted image of each individual galaxy (see Martin et al. 2021 for details). Since dwarf irregulars have significantly higher asymmetries than other classes of dwarf (see Conselice 2014), we define irregular galaxies using an arbitrary cut of A_r > 0.3; we note that we use the regular, un-smoothed, definition of asymmetry rather that the shape asymmetry as described in Martin et al. (2021). The exact value of A_r to be used when compared to observations might differ since A_r is sensitive to the point spread function and resolution in the observations and simulations, respectively. However, the qualitative trend of the fraction of irregulars with stellar mass is robust against realistic variations in the threshold value of A_r (see Appendix C). The fraction of irregular galaxies in NEWHORIZON, shown as the light blue curve in Fig. 16, is consistent with the observational result from Conselice (2006) with more irregular galaxies at the low mass. This is the result of fewer star-forming regions, and thus it provide galaxies with more patchy and more irregular star formation and mass distribution (Faucher-Giguère 2018). The NEWHORIZON data have lately been analysed in terms of morphology by Park et al. (2019), who use the circularity parameter to decompose disc and dispersion components of stars and pin down the origins of the discs and spheroids of spiral galaxies.

3.10. Stellar clusters

As a result of the high spatial resolution of the simulation, clumpy star formation located in large gas complexes is naturally captured (see Fig. 5). Depending on the amount of gas in galaxies and the level of turbulence, star formation can proceed in massive clouds or in a more diffuse fashion within smaller mass clouds. In this work, we measure the fraction of stellar mass locked into stellar clusters using a catalogue of stellar substructures, which puts a minimum detectable stellar cluster mass at 5 × 10⁵ M_⊙. Therefore, it should be noted that a non-negligible fraction of the remaining ‘diffuse’ stellar mass can be contained in stellar clusters with lower mass, passing below our detection threshold of the structure finder. Also, the AdaptaHOP substructure finder considers (sub)structures solely based on the density distribution of stars; therefore stellar clusters can certainly be bound or not with some non-negligible level of contamination (stellar particles from different phase-space) from the background host. Nonetheless, since stars form from high-gas densities, at large efficiencies for bound (α_vir < 1) gas clumps, they form relatively bound until secular evolution (e.g., Gnedin & Ostriker 1997; Baumgardt & Makino 2003) or sudden gas removal (Yu et al. 2020) eventually disrupt them. For this analysis, we also used a lower number of star particles to detect substructures at 10 instead of 50; hence the minimum cluster mass is 10⁵ M_⊙ in that case. The main qualitative evolution with mass and redshift is not affected (see Appendix D for further comments). It also has to be noted that, by construction, the AdaptaHOP galaxy finder cannot remove the most massive structure, which can either be a bulge, but could also be an off-centred stellar cluster in an irregular galaxy. Finally, to minimise the contribution from satellites that are already connected to the main progenitor, only substructures within 2R_eff and with a mass lower than 20% that of the main galaxy are taken into account as a stellar cluster. We note that there is a mix of accreted and in situ formed stellar clusters that contribute to the overall mass of the galaxy and we do not distinguish between these (Mandelker et al. 2014).

Figure 17 shows the mass fraction of stars contained in those relatively massive clumps as a function of galaxy stellar mass and redshift. High-redshift galaxies (z = 4) contain a large fraction (10%) of their stellar mass content within these stellar clusters, while this fraction decreases by an order of magnitude below z < 2. In addition, the low-mass galaxies have a lower fraction of stellar mass within stellar clusters at all redshifts, and at the two lowest redshifts (z = 1 and z = 0.25) there is a trend for the most massive galaxies to decrease their cluster mass fraction with respect to intermediate mass galaxies. At the lowest redshift, many of the individual values of the cluster mass fraction at the low-mass end (M_s ≲ 10⁹ M_⊙) do not appear in the figure since their value is exactly zero. It should be recalled that the exact value of the cluster mass fraction is affected by the capability to capture the formation and survival (see e.g., Pfeffer et al. 2018, and references therein) of the lowest mass stellar clusters due to limited mass and spatial resolution of the simulation, and to detect the stellar overdensities with AdaptaHOP. Despite such resolution effects, the qualitative trend of increasing mass fraction of stars inside cluster with redshift is expected (Elmegreen & Elmegreen 2005; Genzel et al. 2011) as a result of more gas-rich, compact, and turbulent galaxies that are increasingly more gravitationally unstable (Cacciato et al. 2012; Inoue et al. 2016), thereby forming stars into more numerous massive clusters and more efficiently; this is confirmed in Sects. 3.11 and 3.13.

Fig. 17.

Fraction of stellar mass in stellar clusters as a function of galaxy stellar mass at different redshifts as indicated in the panel. The solid lines indicate the mean values, and the dashed lines stand for the error on the mean. The data points are overplotted for the most extreme redshifts to appreciate the scatter in the distribution of sampled values. Galaxies at high redshift are more clumpy than at low redshift at fixed stellar mass.

3.11. Baryonic content

We decompose the baryonic content of galaxies by measuring the amount of stars and gas within twice the effective radius of galaxies; these values are also obtained within R_eff or 0.1R_vir and show similar behaviour. The gas content is further decomposed into a cold and dense gas component (identified by a subscript ‘cd’ in gas quantities) by considering only the star formation gas with density n ≥ 0.1, 10 cm⁻³ and temperature T < 2 × 10⁴ K. In the following, we distinguish between the neutral HI+H₂ gas component with the density cut-off at 0.1 cm⁻³, and the H₂ molecular dense component at 10 cm⁻³ (e.g., Lupi et al. 2018; Nickerson et al. 2019). An exact match of the ionisation and molecular states of the gas however would require a detailed treatment of radiative transfer and molecular chemistry.

In Fig. 18, we show the gas and stellar surface densities as a function of stellar mass and redshift. Surface densities are obtained by dividing the corresponding mass content by π(2R_eff)². All surface densities – gaseous Σ_g, cold Σ_cd, and stellar Σ_s – increase with mass and decrease with redshift; the cold component of the gas surface density decreases faster with time than the total gas component. The decrease of gas and stellar surface densities are consistent with less concentrated galaxies (large effective radius) over time and galaxies with lower sSFR as shown in the previous sections.

Fig. 18.

Surface densities within 2R_eff as a function of galaxy stellar mass at different redshifts for stars (top), HI+H₂ gas (n > 0.1 cm⁻³ and T < 2 × 10⁴ K, middle), and H₂ molecular gas (n > 10 cm⁻³ and T < 2 × 10⁴ K, bottom). The solid lines represent the mean values, and the dashed lines stand for the error on the mean. Galaxies are denser at high redshift at fixed stellar mass and denser in more massive galaxies at fixed redshift.

Figure 19 shows the fraction of baryons locked into the cold gas as a function of galaxy stellar mass for different redshifts and for different gas density cut-offs. The total cold gas fraction (M_cd/(M_cd + M_s) with n > 0.1 cm⁻³; top left panel of Fig. 19) strongly decreases with galaxy mass but does not evolve significantly over redshift. The fraction of cold and denser (n > 10 cm⁻³; top right panel) gas shows a weaker variation with mass, in particular at low redshift, and its value is particularly low compared to the total amount of gas (M_cd/M_g, bottom right), although the fraction of gas made of dense (n > 10 cm⁻³) material increases with stellar mass. Nonetheless, the fraction of cold dense gas decreases with decreasing redshift at a given stellar mass.

Fig. 19.

Ratio of gas over baryon content (top panels) within 2R_eff as a function of galaxy stellar mass at different redshifts for the cold dense star-forming gas with density higher than 0.1 or 10 cm⁻³ (left and right, respectively) and temperature lower than 2 × 10⁴ K, and the ratio of this cold selected gas over total gas (bottom panels). The solid lines indicate the mean values, and the dashed lines stand for the error on the mean. The points with error bars correspond to the observations from Santini et al. (2014; sampled from their interpolated curves in their Fig. 12) at various redshifts and from Catinella et al. (2018) from local galaxies for the HI+H₂ gas, while the dotted line shows the data from Tacconi et al. (2018) at various redshifts for the H₂ molecular gas only. The NEWHORIZON galaxies are less gas-rich at the high-mass end and also at lower redshift for a given stellar mass, in qualitative agreement with the data. But they show significantly less gas at the low-mass end compared to observations, and especially at high redshift with the data from Tacconi et al. (2018), and in contrast they are in good agreement with the data from Santini et al. (2014).

Observations of star-forming H₂ gas, via CO measurements, at z ∼ 1 − 4 have shown that at a given stellar mass, galaxies at high redshifts tend to be significantly more gas rich compared to their present-day counterparts (e.g., Tacconi et al. 2006, 2008, 2010, 2013, 2018, 2020; Daddi et al. 2010a; Scoville et al. 2017). Reported baryonic gas fractions range from ∼20 to 80%, a factor of ∼2–3 more than typically found in cosmological models of galaxy formation (e.g., Popping et al. 2014; Lagos et al. 2015; Davé et al. 2017, 2019). This is confirmed in this work with the comparison of the NEWHORIZON molecular gas to Tacconi et al. (2018) at high redshifts (top right panel of Fig. 19). Narayanan et al. (2012) suggested that the observationally inferred gas fractions of star-forming galaxies at high redshift are overestimated because of the adoption of locally calibrated conversion factors α_CO. When applying a smoothly varying α_CO with the physical properties of galaxies (essentially gas-phase metallicity and CO surface brightness), these authors found gas fractions of ∼10–40% in galaxies of stellar masses in the range 10¹⁰ − 10¹² M_⊙ and a reduced scatter in stellar mass vs. gas fraction by a factor of ∼2, bringing models and observations into a much better agreement. Similarly lower values of gas fractions have been found by Santini et al. (2014), who study main-sequence, star-forming galaxies with M_s > 10¹⁰ M_⊙ at z < 2 and infer gas content from dust mass measurement, or by Conselice et al. (2013), who derive the cold gas mass fraction by inverting the global Kennicutt–Schmidt relation for massive galaxies (M_s > 10¹¹ M_⊙) at 1.5 < z < 3. In the range of comparable stellar masses, gas fractions of high-z NEWHORIZON galaxies agree broadly with these revised measurements, hosting ∼10–40% of star-forming gas that decreases with increasing stellar mass and evolves weakly with redshift. In Fig. 19, we show a few data points obtained from the sampling of the interpolated curves from Santini et al. (2014) (see their Fig. 12). While their dust method used to derive the gas content of galaxies in principle traces both atomic and molecular components, the authors suggest that the bulk of the gas in studied galaxies is in the molecular phase. We also show in Fig. 19, the local scaling relation from SDSS galaxies from Catinella et al. (2018). The NEWHORIZON cold gas fractions agree fairly well with both local and high-z gas scaling relations for galaxy masses above a few 10¹⁰ M_⊙, but significantly differ with the local scaling relation at the low-mass end.

3.12. Tully–Fisher relation

To measure the baryonic Tully & Fisher (1977) relation (e.g., Bell & de Jong 2001; McGaugh & Schombert 2015; Ponomareva et al. 2018; Lelli et al. 2019) in simulated galaxies, we use the decomposition of the gas kinematics in a cylindrical frame of reference. In this frame, rotation velocity profile is measured using only (cold) gas resolution elements with densities larger than 0.1 cm⁻³ and temperatures below 2 × 10⁴ K. The flat rotational velocity involved in the baryonic Tully–Fisher relation V_g, f is obtained by measuring the average gas rotational velocity within [1.8R_eff,2.2R_eff]. The total baryonic mass corresponds to the total stellar and cold gas mass M_b = M_s + M_g within 2R_eff. A small fraction of galaxies in the sample do not have cold gas that reaches the 2R_eff radius. In that case, the measured kinematics and baryonic mass are replaced by their value at R_eff. We select only disc galaxies with V/σ ≥ 0.5 that are hosted within DM halos more massive than 10⁹ M_⊙; however, considering non-disc galaxies leads to a similar relation.

Figure 20 shows the mean baryonic Tully–Fisher relation at different redshifts z = 4, 2, 1, 0.25 in NEWHORIZON, with the two observational fits by McGaugh & Schombert (2015) with m ∼ 4 ( $V_{g, f} \propto M_{b}^{1 / m}$ $V_{\mathrm{g,f}}\propto M_{\mathrm{b}}^{1/m}$ ) from their different combined samples of disc galaxies, and with the observational fit from Ponomareva et al. (2018) with m ∼ 3. The obtained NEWHORIZON baryonic Tully–Fisher relation shows very little evolution with redshift; the mean relations are almost indistinguishable in Fig. 20. Although the bulk of the relation is captured, there are a few different behaviours of the NEWHORIZON sample with the observations, in particular with respect to the McGaugh & Schombert (2015) relations: at masses between 10⁸ ≲ M_b/M_⊙ ≲ 10¹⁰ M_⊙, NEWHORIZON galaxies have velocities below the observed values, and the slope of the relation V_g, f–M_b becomes steeper for masses larger than 10¹⁰ M_⊙. A similar discrepancy (bending of the relation) with observations was also found in NIHAO simulations (Dutton et al. 2017) or APOSTLE/EAGLE (Sales et al. 2017), although with different overall normalisation of the relation relative to ours. The baryonic Tully–Fisher relation obtained from our NEWHORIZON data at z = 0.25 produces a slope of m = 2.8 ± 0.1 that agrees best with the m ∼ 3 slope from Ponomareva et al. (2018). It has to be noted that the observational slope of the relation varies with the adopted stellar mass-to-light ratio: whether a (molecular) gas component remains undetected, depending on the mass range of the sample (see e.g., Ponomareva et al. 2018, for a discussion along these lines), and how the velocity is measured. For instance in Lelli et al. (2019), this value is m = 3.14 if the velocity is measured at 2R_eff as we adopted in this work, while it is m = 3.85 for the velocity in the flat part of the radial velocity profile. For comparison, the baryonic Tully–Fisher in SIMBA (Glowacki et al. 2020) differs significantly from the value in this work with a value of m ≃ 4 rather than m ≃ 3; the SIMBA value has a larger normalisation of the relation than in McGaugh & Schombert (2015) as opposed to this work, where the normalisation is lower in most of the mass range. Glowacki et al. (2020) also show that the value of the slope is sensitive to the velocity estimator. In the EAGLE simulation, Ferrero et al. (2017) show that the stellar Tully–Fisher relation is in excellent agreement with the observational data, even though the stellar (using M_s), as opposed to the baryonic (using M_b), Tully–Fisher relation, provides a less tight relation between the mass and the rotational velocity of galaxies (McGaugh 2012).

Fig. 20.

Baryonic Tully–Fisher relation for disc galaxies in NEWHORIZON at different redshifts as labelled in the panel with the raw data points overplotted as crosses for redshift z = 0.25 (purple plus signs), compared to the two best fits from McGaugh & Schombert (2015; dashed black and dot-dashed black) and the best fit from Ponomareva et al. (2018; solid black; in thick their fit over the range of their data points, while the thin part is the extrapolated part). The fit to the NEWHORIZON data at z = 0.25 (thick dotted purple line) shows a better level of agreement with the data from Ponomareva et al. (2018; with a slope of one-third) over those of McGaugh & Schombert (2015; slope of one-fourth).

3.13. Gas kinematics

The kinematics of the ISM is computed by measuring gas velocities within twice the effective radius of each galaxy and for the cold non-star-forming gas only (i.e., 10⁻² ≤ n/cm⁻³ < 10 and T ≤ 2 × 10⁴ K), which has similar properties to that observed in Weiner et al. (2006) or Simons et al. (2017), who estimated the gas kinematics from nebular emission lines (amongst which Hα and [O III]λ5007). The kinematics of the gas are decomposed along the cylindrical system of coordinates, which for each galaxy is given by the angular momentum of the selected gas elements. Therefore, the velocity dispersion is decomposed along the three components σ_g, r, σ_g, z, and σ_g, t and is simply the mass-weighted dispersion of each velocity component around the mean value of that component; the total velocity dispersion of the gas is obtained from the three components $σ_{g}^{2} = (σ_{g, r}^{2} + σ_{g, z}^{2} + σ_{g, t}^{2}) / 3$ $\sigma_{\rm g}^2=(\sigma_{{\rm g},r}^2+\sigma_{{\rm g},z}^2+\sigma_{{\rm g},t}^2)/3$ . Similarly the rotational velocity is obtained by taking mass-weighted mean tangential velocity component $V_{g, rot} = {\bar{V}}_{g, t}$ $V_{\mathrm{g, rot}}=\bar{V}_{\mathrm{g}, t}$ . A proxy for the total kinematics support is built out of the dispersion and rotation using $S_{g, 0.5} = {(0.5 V_{g, rot}^{2} + σ_{g}^{2})}^{1 / 2}$ $S_{\rm g,0.5}=(0.5V_{\rm g,rot}^2+\sigma_{\rm g}^2)^{1/2}$ (e.g., Weiner et al. 2006; Kassin et al. 2007). Although synthetic mocks of nebular emission lines with instrumental effects from which observed kinematics are reconstructed would certainly show differences with the direct gas kinematics measured from the simulation, these raw measurements on NEWHORIZON galaxies have been tested against projection effects on σ_g (i.e., vs. σ_g, z), density cut-offs (i.e., vs. selected gas with n > 0.1 cm⁻³ and T ≤ 2 × 10⁴ K) or aperture (i.e., vs. within R_eff), which changes the estimated averaged values shown in Fig. 21 by at most 10–20% and does not affect any of the trends obtained with mass and redshift.

Fig. 21.

Gas kinematics as a function of galaxy stellar mass at different redshifts for gas density 0.01 < n < 10 cm⁻³, temperature lower than T < 2 × 10⁴ K and within 2R_eff. From left to right and top to bottom: rotational velocity, velocity dispersion, velocity dispersion over S_g, 0.5, and rotational velocity over S_g, 0.5. The error bars stand for the error on the mean. The dashed plus signs in the top right panel indicate the values with uncertainties from the fit to the observational data in Simons et al. (2017). The dotted lines on the bottom panels stand for the cases of pure rotational support or pure dispersion support, respectively. The NEWHORIZON galaxies decrease their amount of gas velocity dispersion over time for a given stellar mass, while keeping a similar amount of rotation.

Figure 21 shows the amount of dispersion and rotational support of the cold gas as a function of galaxy mass for different redshifts. The dispersion, rotation, and total support all increase with galaxy mass at any given redshift; however, the behaviour of these quantities over time differ significantly. The velocity dispersion of the cold gas decreases with decreasing redshift at fixed galaxy mass, such as for galaxies with stellar mass M_s = 10⁹ M_⊙, values as high as 60 km s⁻¹ (and up to 100 km s⁻¹ for the most massive galaxies at high redshift) at z = 4 and 20 km s⁻¹ at z = 0.25, which agrees well with the observations (Simons et al. 2017). In agreement with the observations, the rotation of cold gas in galaxies as a function of mass does not significantly evolve over time, although the small evolution with time is opposite to the small evolution in observations. In turn, the cold gas in galaxies is proportionally more supported by rotation (V_g, rot/S_g, 0.5) over time as an effect of reduced gas velocity dispersion that is qualitatively consistent with the observations (Simons et al. 2017). The ratio of V_g, rot/S_g, 0.5 increases over time (while σ_g/S_g, 0.5 decreases) and also has a sharp transition, that is an increase and a decrease, respectively, at stellar masses above 10¹⁰ M_⊙ (Simons et al. 2015). The obtained values of these ratios are in good quantitative agreement with Simons et al. (2017) for the high-mass range (M_s > 10¹⁰ M_⊙), but V_g, rot/S_g, 0.5 at M_s < 10¹⁰ M_⊙ is two standard deviations larger the observational relation at z = 1 and 2, and lower at z = 0.25. The turbulent support (σ_g/S_g, 0.5) agrees well with the observations except for the low mass range M_s < 10¹⁰ M_⊙ at z = 0.25, where the amount of turbulent support is significantly larger than observed. This decrease of the gas dispersion and hence of the ratio of dispersion to rotation over time and mass triggers the settling of galactic discs (Kassin et al. 2007, 2014; Ceverino et al. 2017; Simons et al. 2017; Pillepich et al. 2019) as galaxies become more quiescent (see Dubois et al., in prep.). Compared to TNG50 (Pillepich et al. 2019), NEWHORIZON galaxy gas kinematics produce smaller rotational velocities and higher turbulent velocities at all redshifts: the rotational velocity in NEWHORIZON is about half of that in TNG and the turbulent velocities about twice higher than in TNG. Interestingly, TNG50 gas fractions are also a factor two smaller than in NEWHORIZON, with galaxies that have a similar size-mass relation. Therefore, the gas Toomre parameter Q_g, which becomes simply proportional to σ_g/Σ_g, is of the same value, pointing towards a common saturation mechanism of the (gas) Toomre parameter despite the different models of feedback.

Turbulence can be sustained by several sources including stellar feedback (Joung et al. 2009; Ostriker & Shetty 2011; Hopkins et al. 2014), gravitational instabilities (Agertz et al. 2009; Bournaud et al. 2010), or cosmological infall (Elmegreen & Burkert 2010; Klessen & Hennebelle 2010) ranging from anisotropic gas infall or mergers. It is not clear which process dominates the driving of the turbulence. However Krumholz & Burkhart (2016) and Krumholz et al. (2018) show that observations favour a gravity-driven source for gas-rich (or high-z) galaxies and a stellar-driven source of turbulence for gas-poor (low-z) galaxies, although the role of cosmic infall is not tested in those models.

3.14. MBH-to-galaxy mass relation

In NEWHORIZON, MBH formation commonly occurs even in low-mass galaxies, as can be seen by the fact that about half of galaxies with a stellar mass 10⁷ − 10⁸ M_⊙ at redshift z = 0.25 host at least one BH. Observationally, a statistical analyses based on X-ray observations (Miller et al. 2015; She et al. 2017) suggest that between 10% and 100% of galaxies with mass ∼10⁹ M_⊙ at z = 0 host a MBH. While using dynamical masses or upper limits in nucleated local galaxies within 4 Mpc, Nguyen et al. (2018) find that between 50 and 100% of galaxies with mass ∼10⁹ − 10¹⁰ M_⊙ contain a MBH. By redshift z = 0.25, NEWHORIZON contains 583 MBHs that are located within two effective radii of a galaxy included in the catalogue. Each MBH can only be associated with a single galaxy at a given point in time, but a galaxy can contain multiple black holes and this particularly occurs for the most massive galaxies in NEWHORIZON. This is reflected in the occupation fraction (Volonteri et al. 2008), which increases with galaxy stellar mass from 0.28 at M_s = 10⁶ M_⊙ to unity at M_s = 10^8.5 M_⊙. At higher galaxy masses there is on average more than one MBH per galaxy, reaching an average of two MBHs per galaxy at M_s = 10¹⁰ M_⊙. The existence of additional ‘off-centre’ MBHs in massive galaxies is a standard prediction of models and simulations that study the cosmic evolution of MBHs in galaxies (e.g., Islam et al. 2003; Bellovary et al. 2010; Volonteri et al. 2016; Tremmel et al. 2018). The relatively high occupation fraction in low-mass galaxies allows for MBH mergers to also occur in low-mass galaxies, following galaxy mergers. This is studied in further detail in Volonteri et al. (2020).

Not all MBHs are located close to the centre of their host galaxy. Overall, 38% are located within 2.5 kpc of the centre, while the rest can be found on larger orbits. This agrees with other recent simulations of MBHs in dwarf galaxies (Bellovary et al. 2019; Pfister et al. 2019). In Bellovary et al. (2019), it is reported that 50% of MBHs in dwarf galaxies wander up to several kiloparsec from the galaxy centre. Wandering MBHs have also been found observationally, both in dwarf galaxies (Reines et al. 2020) and in more massive galaxies (Menezes et al. 2014, 2016; Shen et al. 2019). Maximum offset distances reported in observations are smaller than in simulations because they preferentially detect active MBHs, which biases their sample to centrally located MBHs that are fed by the gas supply of their host galaxy. In contrast, our sample includes all MBHs, whether active or not.

As previously shown in Dubois et al. (2015; see also Habouzit et al. 2017; Bower et al. 2017; Prieto et al. 2017; Anglés-Alcázar et al. 2017; McAlpine et al. 2018; Trebitsch et al. 2018, 2020), we confirm that MBH growth is regulated by SN feedback for galaxies with stellar masses below M_s < 5 × 10⁹ M_⊙, which make up the bulk of the stellar population in NEWHORIZON. As shown in Fig. 22, as a consequence of this most of the MBH population grows little over the course of the simulation. Across all stellar mass bins, MBHs in NEWHORIZON are under-massive in comparison to observations, lying up to two orders of magnitude below the bulk of the observed MBH masses. However, MBHs in more massive galaxies (> 10¹¹ M_⊙) lie closer to the locus of observational data than those in low-mass galaxies, suggesting that MBHs in NEWHORIZON are still in the process of catching up with the observed relation.

Fig. 22.

MBH mass as a function of galaxy mass for redshifts z = 4, 2, 1 and 0.25 as circles for all MBHs contained within 2 effective radii of their host galaxy. For the lowest redshift, z = 0.25, big circles (dark red) highlight the most massive MBH for a given galaxy, while the small circles (dark red) show all the secondary MBHs within the galaxy. The black line denotes the mean M_BH for all primary BH within a given M_S bin at z = 0.25 for NEWHORIZON. Also shown are observations of MBH vs. stellar mass for z ∼ 0 − 0.3 from Reines & Volonteri (2015; RV15; green triangles), Greene et al. (2020; Greene20; blue markers) and Baron & Ménard (2019; BM19; grey contours). The error bars in RV15 were omitted for clarity. On average, central MBHs in NEWHORIZON grow significantly only above a stellar mass threshold of a few 10⁹ M_⊙, with non-central MBHs generally failing to grow even in galaxies above this mass threshold.

By redshift z = 0.25, 79% of MBHs have retained a mass that remains within a factor 2 of their seed mass (chosen to be 10⁴ M_⊙ in NEWHORIZON). It is only once their host galaxy reaches the transition stellar mass of M_s > 5 × 10⁹ M_⊙ that the main MBH of a galaxy (big circles, here defined to be the most massive MBH for each galaxy) grows efficiently and soon reaches the observed M_BH − M_s relation (triangles and contours) because the deepening gravitational potential of their host galaxies allows for nuclear inflows to start efficiently feeding the central MBH. All secondary MBHs associated with a galaxy (small dark red circles for z = 0.25) continue to struggle to grow as they are too far from the centre of their galaxy to efficiently access this increased gas supply.

3.15. MBH spins

In NEWHORIZON, MBH spin is followed on the fly, taking into account the effects of gas accretion and MBH–MBH mergers. As shown in Fig. 23, even though MBHs are initially formed with zero spin, they spin up quickly as their mass doubles, which is in line with earlier findings that gas accretion leads to maximally spinning MBH (Dubois et al. 2014b,c; Bustamante & Springel 2019). From there, spin evolution becomes more complex because mass growth is accomplished through MBH–MBH mergers and gas accretion and has the broadest scatter in the mass range M_BH = 2 × 10⁴ − 10⁵ M_⊙. The MBHs are formed with zero spin and initial mass M_BH = 10⁴ M_⊙, and for masses up to 2 × 10⁴ M_⊙ the only way spin can change is by gas accretion. We note that at low accretion rates χ < 0.01 jets spin MBHs down, while at higher accretion rates gas accretion can either spin MBHs up or down, depending on the relative orientation of the spin and the gas feeding the MBH. As shown in Fig. 23, for MBHs with mass up to 2 × 10⁴ M_⊙, spin increases with mass, implying that gas accretion with χ ≥ 0.01 is responsible for increasing MBH spins. Most of accretion occurs in prograde fashion, meaning that gas and MBH spin are initially aligned or align during the accretion episode, which is in line with earlier findings that gas accretion leads to maximally spinning MBH (Dubois et al. 2014c). Above M_BH = 2 × 10⁴ M_⊙ spin can be affected by the combination of accretion and MBH–MBH mergers, which can modify abruptly both spin magnitude and orientation. The MBH–MBH mergers are responsible for the large scatter in the mass range M_BH = 2 × 10⁴ − 10⁵ M_⊙. As shown in the coloured evolution histories, individual MBHs can repeatedly spin up and down, but as they grow more massive all MBHs in our sample spin up to spin values larger than 0.7. All MBHs with mass > 10⁶ M_⊙ have acquired more than 80% of their mass through accretion, which is once again responsible for increasing MBH spins. This is in good agreement with observations, which find high spin values for the mass range M_BH = 10⁶ − 5 × 10⁷ M_⊙ (Reynolds 2013). For higher MBH masses, spins are expected to turn over again as MBH mergers become increasingly important in MBH growth, but no MBHs in our sample grows sufficiently massive to probe this regime. The impact of mergers and gas accretion on the mass and spin evolution of the MBH population will be explored further in a companion paper (Beckmann et al., in prep.).

Fig. 23.

MBH spin evolution with MBH mass for all MBH at redshift z = 0.25 (black markers), and the spin evolution history with mass for the six most massive MBH (coloured lines). Shown for comparison are observations from Reynolds (2013; grey triangles) and Soares & Nemmen (2020; grey crosses). MBHs experience a fast initial spin growth phase followed by a more stochastic evolution.

4. Conclusions

The NEWHORIZON tool provides a compromise to alternative strategies, namely large-volume cosmological simulations with poor spatial resolution within galaxies or selected high-resolution zoomed-in halos, by simulating an intermediate sub-volume (16 Mpc)³ of average cosmic density within a larger box that is resolved at a spatial resolution good enough (34 pc) to capture the turbulence and multi-phase structure of the ISM. Despite a still important effect of the cosmic variance on the various estimated properties, we believe this compromise between sampling a diversity of environments and high resolution is an important step to faithfully capture the relevant physical processes involving star formation and to improve our understanding of the various mechanisms at work that shape galaxies over cosmic times. We have shown that NEWHORIZON reproduces with a good level of accuracy several of the key properties that define galaxies, which we briefly review:

Galaxies simulated in NEWHORIZON naturally exhibit a multi-phase structure with dense cold star-forming clumps embedded in warm and hot diffuse gas, where turbulence shape star-forming properties of galaxies.
The galaxy mass function with the characteristic turnover at 10¹⁰ − 10¹¹ M_⊙ is well reproduced at low redshift when surface brightness limitations are accounted for.
The cosmic SFR density and stellar density compare reasonably well with observations, although they are higher than in observations (up to a factor of 2 depending on the observations), and with a significant uncertainty in the predicted value of the cosmic SFR (0.7 dex) owing to cosmic variance. The sSFR as a function of galaxy stellar mass was measured at various redshifts and a comparison with the observations shows that the simulated sSFRs are in a fair agreement with the observations, showing that galaxies were more active in the past.
The simulated Kennicutt–Schmidt relation points are in fair agreement with the observed main sequence although NEWHORIZON galaxies are a factor ≃0.3 dex below the observations. High-redshift galaxies show a larger amount of starburst galaxies above the main sequence in qualitative agreement with the data.
The stellar-to-halo mass relations show that galaxies are relatively too massive in NEWHORIZON at the low-mass end (M_h< a few 10¹¹ M_⊙) compared to observational data. The relation is in good agreement at the Milky Way mass scale.
The size-mass relation of galaxies falls well within the range of observations with a significant positive evolution of galaxy sizes over decreasing redshift for a given stellar mass. At the high-mass end, NEWHORIZON produces more compact galaxies than expected and no extended ellipticals, although there is a clear lack of statistics in this mass range.
The metallicities of stars and gas agree reasonablly well with the observations in mass and redshift, despite a crude model for chemical evolution. However, the stellar metallicities in NEWHORIZON show less variation with redshift than in the observations.
NEWHORIZON galaxies display an increase in stellar rotation over dispersion over mass and decreasing redshift. The fraction of ellipticals based on their stellar kinematics compares well with observations, although the exact fraction is sensitive to the radius within which the kinematics are measured along with the adopted threshold.
Simulated galaxies at high redshift have a larger fraction of their stellar mass enclosed in stellar clumps than at low redshift for a given stellar mass, as a result of more gas-rich, turbulent, and concentrated galaxies.
The NEWHORIZON galaxies have lower stellar and gas surface densities over decreasing redshift. The fraction of cold star-forming gas is decreasing over decreasing redshift at a given stellar mass, which makes galaxies less gravitationaly unstable.
The bulk of the Tully–Fisher relation is fairly well captured, where a slope of one-third is preferred over a slope of one-fourth.
The gas in NEWHORIZON galaxies settles into discs over decreasing redshift as a result of the decaying level of gas turbulence, while gas rotation remains similar (at constant mass).
The MBH-to-galaxy mass relation shows a sharp increase at a stellar mass of a few 10⁹ M_⊙. Observations are broadly reproduced, although MBH masses fall short.
Once MBHs manage to grow in mass their spin is also large (a > 0.7 − 0.8) and has a sharp rise during the initial growth phase with a more chaotic evolution during self-regulation of the MBH mass, leading to a fair agreement of the MBH spin-mass relation with the scarce data in this mass range.

This paper is the first step of an introduction to NEWHORIZON by reviewing the bulk properties of galaxies (stars and gas) and of their hosted MBHs. In particular, we have not investigated the properties of the circum-galactic medium of galaxies (including hot diffuse and cold flows) and how it interacts with the galactic outflows, the shape of the DM distribution within halos and galaxies, inner substructures within galaxies (bars, bulges and pseudo-bulges, thin and thick discs, and properties of star-forming clouds of gas). In addition, the NEWHORIZON set-up limits our study to a limited volume within a fairly average cosmic density, and similar studies should extend this work to various environments such as dense galaxy clusters (see e.g., Trebitsch et al. 2020, for a first attempt in a proto-cluster), voids, or within large-scale tens-of-Mpc wide cosmic filaments. Since the properties of galaxies in NEWHORIZON are expected to depend on the adopted sub-grid models (see for instance, Dubois et al. 2012, 2016; Kimm et al. 2015; Chabanier et al. 2020b; Nuñez-Castiñeyra et al. 2021) and resolutions (on at least some aspects, numerical resolution even in the diffuse large-scale medium can be more important than the parametrisation of sub-grid models, Chabanier et al. 2020a), it would be ideal to explore the effects of critical parameters systematically in order to assess the robustness of the results. However, it is beyond the available computing resources considering that the combination of the volume and resolution of NEWHORIZON is already pushed to the limit.

¹

See Park et al. (2019, 2021), Volonteri et al. (2020), Martin et al. (2021), Jackson et al. (2021a), and Jackson et al. (2021b) for early results on the origin of discs and spheroids, the thickness of discs, the mergers of black holes, the role of interactions in the evolution of dwarf galaxies, the DM deficient galaxies, and low-surface brightness dwarf galaxies, respectively.

²

http://new.horizon-simulation.org

³

http://horizon-simulation.org

⁴

To avoid arbitrary volume correction to select for purity of the zoom-in region, we use the whole galaxy sample and the initial volume of the simulation, that is (16 Mpc)³ to measure the volumetric quantities of this section.

⁵

The cosmic matter density of the simulated zoom-in region is a factor 1.2 in excess to the average cosmic matter density, which contributes to the total excess in cosmic SFR and stellar density.

⁶

The sharp increase at log(1 + z)≃0.25 in the cosmic SFR density of the M_s/M_⊙ ≥ 10¹¹ mass bin corresponds to a galaxy moving from the 10¹⁰ ≤ M_s/M_⊙ < 10¹¹ mass bin to the previous mass bin; there is a dominant contribution from this single galaxy to the cosmic SFR density in that mass bin.

⁷

Changing the timescale over which SFRs are measured to 10 Myr increases the uncertainty in the mean relation of the simulated points without changing the trends. Changing the measurement radius to 2R_eff or 3R_eff decreases the mean sSFRs by a few tens of per cent.

⁸

The observational values of Gallazzi et al. (2014) and Panter et al. (2008) are given for an assumed solar metallicity of Z_⊙ = 0.02 (Anders & Grevesse 1989), while more self-consistent calculations of the composition of the solar atmosphere give a significantly lower value of Z_⊙ = 0.01345 (Asplund et al. 2009). We have, therefore, scaled up their fitting relations accordingly.

⁹

Although a haze of stars from the companion is still seen and associated with the main object. This component could be removed by using a more robust structure finder that uses information from the velocity distribution (Cañas et al. 2019).

Acknowledgments

Y.D., S.K.Y., J.D., T.K., C.P., S.P., and M.V. conceived and planned the NEWHORIZON simulation project. Y.D., J.D. and T.K. prepared the preparatory runs and Y.D., S.K.Y., H.C., K.K. and F.B. ran the simulation. Y.D., J.D., T.K. and M.V. developed the models for galaxy formation. C.P. produced the high-resolution initial conditions and managed the website. R.B. and M.V. contributed to the analysis of MBH properties. K.K. contributed to the Kennicutt–Schmidt and gas fraction analysis. H.C. generated mock images of galaxies. M.-J.P. and G.M. contributed to kinematic analyses on galaxies. G.M. and S.K. contributed to the analysis on galaxy mass functions. R.J., G.M. and S.K. contributed to the analysis on galaxy surface brightness. C.L. contributed to the cosmic SFR analysis. Y.D. produced halo and galaxy catalogues and contributed to all the analysis of the paper. The manuscript has been written by Y.D. with contributions from all co-authors. We thank the Camille Noûs collective, and thus declare our support for the activities of this embodiment of the academic community as a whole. The collective gathers scientists from all horizons wishing to symbolically recall that the scientific edifice is a collegial construction by nature, i.e. that any publication is deeply impregnated with the great body of human knowledge generated by our past and todays colleagues. We thank Olivier Ilbert for assistance using LEPHARE. We thank the anonymous referee for his/her constructing comments that improved the clarity of this paper. This work was granted access to the HPC resources of CINES under the allocations c2016047637, A0020407637 and A0070402192 by Genci, KSC-2017-G2-0003 by KISTI, and as a “Grand Challenge”? project granted by GENCI on the AMD Rome extension of the Joliot Curie supercomputer at TGCC. This research is part of the Spin(e) ANR-13-BS05-0005 (http://cosmicorigin.org), Segal ANR-19-CE31-0017 (http://secular-evolution.org) and Horizon-UK projects. This work has made use of the Horizon cluster on which the simulation was post-processed, hosted by the Institut d’Astrophysique de Paris. We warmly thank S. Rouberol for running it smoothly. The large data transfer was supported by KREONET which is managed and operated by KISTI. S.K.Y. acknowledges support from the Korean National Research Foundation (NRF-2020R1A2C3003769). TK was supported in part by the National Research Foundation of Korea (NRF-2017R1A5A1070354 and NRF-2020R1C1C1007079) and in part by the Yonsei University Future-leading Research Initiative (RMS2-2019-22-0216). This work is partly based on tools and data products produced by GAZPAR operated by CeSAM-LAM and IAP. This paper is based in part on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center (ADC) at National Astronomical Observatory of Japan. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org Catalogues extracted from the simulation will be made available at this URL https://new.horizon-simulation.org/data.html.

References

Agertz, O., Lake, G., Teyssier, R., et al. 2009, MNRAS, 392, 294 [NASA ADS] [CrossRef] [Google Scholar]
Agertz, O., Teyssier, R., & Moore, B. 2011, MNRAS, 410, 1391 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25 [NASA ADS] [CrossRef] [Google Scholar]
Agertz, O., Renaud, F., Feltzing, S., et al. 2021, MNRAS, 503, 5826 [CrossRef] [Google Scholar]
Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, PASJ, 71, 114 [CrossRef] [Google Scholar]
Algorry, D. G., Navarro, J. F., Abadi, M. G., et al. 2014, MNRAS, 437, 3596 [NASA ADS] [CrossRef] [Google Scholar]
Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
Andersson, E. P., Agertz, O., & Renaud, F. 2020, MNRAS, 494, 3328 [CrossRef] [Google Scholar]
André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Anglés-Alcázar, D., Faucher-Giguère, C.-A., Quataert, E., et al. 2017, MNRAS, 472, L109 [NASA ADS] [CrossRef] [Google Scholar]
Arnouts, S., & Ilbert, O. 2011, Astrophysics Source Code Library [record ascl:1108.009] [Google Scholar]
Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540 [NASA ADS] [CrossRef] [Google Scholar]
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376 [NASA ADS] [CrossRef] [Google Scholar]
Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621 [NASA ADS] [Google Scholar]
Baron, D., & Ménard, B. 2019, MNRAS, 487, 3404 [CrossRef] [Google Scholar]
Baumgardt, H., & Makino, J. 2003, MNRAS, 340, 227 [NASA ADS] [CrossRef] [Google Scholar]
Beckmann, R. S., Devriendt, J., Slyz, A., et al. 2017, MNRAS, 472, 949 [NASA ADS] [CrossRef] [Google Scholar]
Beckmann, R. S., Slyz, A., & Devriendt, J. 2018, MNRAS, 478, 995 [CrossRef] [Google Scholar]
Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [Google Scholar]
Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143 [Google Scholar]
Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212 [NASA ADS] [CrossRef] [Google Scholar]
Bellovary, J. M., Governato, F., Quinn, T. R., et al. 2010, ApJ, 721, L148 [NASA ADS] [CrossRef] [Google Scholar]
Bellovary, J. M., Cleary, C. E., Munshi, F., et al. 2019, MNRAS, 482, 2913 [NASA ADS] [Google Scholar]
Benson, A. J., & Babul, A. 2009, MNRAS, 397, 1302 [CrossRef] [Google Scholar]
Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, MNRAS, 436, 697 [Google Scholar]
Bernstein, G. M., & Jarvis, M. 2002, AJ, 123, 583 [Google Scholar]
Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140 [NASA ADS] [CrossRef] [Google Scholar]
Béthermin, M., Wang, L., Doré, O., et al. 2013, A&A, 557, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bielby, R., Hudelot, P., McCracken, H. J., et al. 2012, A&A, 545, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bieri, R., Dubois, Y., Rosdahl, J., et al. 2017, MNRAS, 464, 1854 [NASA ADS] [CrossRef] [Google Scholar]
Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [Google Scholar]
Blanton, M. R., Lupton, R. H., Schlegel, D. J., et al. 2005, ApJ, 631, 208 [Google Scholar]
Blondin, J. M., Wright, E. B., Borkowski, K. J., & Reynolds, S. P. 1998, ApJ, 500, 342 [NASA ADS] [CrossRef] [Google Scholar]
Bois, M., Emsellem, E., Bournaud, F., et al. 2011, MNRAS, 416, 1654 [NASA ADS] [CrossRef] [Google Scholar]
Bothwell, M. S., Kennicutt, R. C., & Lee, J. C. 2009, MNRAS, 400, 154 [NASA ADS] [CrossRef] [Google Scholar]
Bouché, N., Cresci, G., Davies, R., et al. 2007, ApJ, 671, 303 [NASA ADS] [CrossRef] [Google Scholar]
Bournaud, F., Jog, C. J., & Combes, F. 2007, A&A, 476, 1179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bournaud, F., Elmegreen, B. G., Teyssier, R., Block, D. L., & Puerari, I. 2010, MNRAS, 409, 1088 [NASA ADS] [CrossRef] [Google Scholar]
Bournaud, F., Perret, V., Renaud, F., et al. 2014, ApJ, 780, 57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bower, R. G., Schaye, J., Frenk, C. S., et al. 2017, MNRAS, 465, 32 [NASA ADS] [CrossRef] [Google Scholar]
Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40 [CrossRef] [Google Scholar]
Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]
Bustamante, S., & Springel, V. 2019, MNRAS, 490, 4133 [CrossRef] [Google Scholar]
Cañas, R., Elahi, P. J., Welker, C., et al. 2019, MNRAS, 482, 2039 [CrossRef] [Google Scholar]
Cacciato, M., Dekel, A., & Genel, S. 2012, MNRAS, 421, 818 [NASA ADS] [Google Scholar]
Camps, P., & Baes, M. 2020, Astron. Comput., 31, 100381 [CrossRef] [Google Scholar]
Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2012, Nature, 484, 485 [NASA ADS] [CrossRef] [Google Scholar]
Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875 [NASA ADS] [CrossRef] [Google Scholar]
Ceverino, D., & Klypin, A. 2009, ApJ, 695, 292 [Google Scholar]
Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 404, 2151 [NASA ADS] [Google Scholar]
Ceverino, D., Primack, J., Dekel, A., & Kassin, S. A. 2017, MNRAS, 467, 2664 [NASA ADS] [CrossRef] [Google Scholar]
Chabanier, S., Bournaud, F., Dubois, Y., et al. 2020a, A&A, 643, L8 [CrossRef] [EDP Sciences] [Google Scholar]
Chabanier, S., Bournaud, F., Dubois, Y., et al. 2020b, MNRAS, 495, 1825 [CrossRef] [Google Scholar]
Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
Chabrier, G. 2005, in The Initial Mass Function: From Salpeter 1955 to 2005, eds. E. Corbelli, F. Palla, & H. Zinnecker, Astrophys. Space Sci. Lib., 327, 41 [Google Scholar]
Chevalier, R. A. 1974, ApJ, 188, 501 [NASA ADS] [CrossRef] [Google Scholar]
Chiosi, C., Bertelli, G., & Bressan, A. 1992, ARA&A, 30, 235 [NASA ADS] [CrossRef] [Google Scholar]
Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252 [NASA ADS] [CrossRef] [Google Scholar]
Coccato, L., Fabricius, M., Morelli, L., et al. 2015, A&A, 581, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Conselice, C. J. 2006, MNRAS, 373, 1389 [NASA ADS] [CrossRef] [Google Scholar]
Conselice, C. J. 2014, ARA&A, 52, 291 [NASA ADS] [CrossRef] [Google Scholar]
Conselice, C. J., Bershady, M. A., & Jangren, A. 2000, ApJ, 529, 886 [NASA ADS] [CrossRef] [Google Scholar]
Conselice, C. J., Mortlock, A., Bluck, A. F. L., Grützbauch, R., & Duncan, K. 2013, MNRAS, 430, 1051 [Google Scholar]
Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937 [NASA ADS] [CrossRef] [Google Scholar]
Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11 [NASA ADS] [CrossRef] [Google Scholar]
Crutcher, R. M. 2012, ARA&A, 50, 29 [NASA ADS] [CrossRef] [Google Scholar]
Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944 [Google Scholar]
Daddi, E., Bournaud, F., Walter, F., et al. 2010a, ApJ, 713, 686 [Google Scholar]
Daddi, E., Elbaz, D., Walter, F., et al. 2010b, ApJ, 714, L118 [Google Scholar]
Dalgarno, A., & McCray, R. A. 1972, ARA&A, 10, 375 [NASA ADS] [CrossRef] [Google Scholar]
Dalla Vecchia, C., & Schaye, J. 2008, MNRAS, 387, 1431 [Google Scholar]
Dashyan, G., & Dubois, Y. 2020, A&A, 638, A123 [EDP Sciences] [Google Scholar]
Davé, R., Rafieferantsoa, M. H., Thompson, R. J., & Hopkins, P. F. 2017, MNRAS, 467, 115 [NASA ADS] [Google Scholar]
Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827 [NASA ADS] [CrossRef] [Google Scholar]
Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2 [NASA ADS] [CrossRef] [Google Scholar]
De Rossi, M. E., Bower, R. G., Font, A. S., Schaye, J., & Theuns, T. 2017, MNRAS, 472, 3354 [Google Scholar]
Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2 [NASA ADS] [CrossRef] [Google Scholar]
Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870 [Google Scholar]
Dekel, A., & Silk, J. 1986, ApJ, 303, 39 [Google Scholar]
Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451 [Google Scholar]
Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Dobbs, C. L., Burkert, A., & Pringle, J. E. 2011, MNRAS, 413, 2935 [Google Scholar]
Donnari, M., Pillepich, A., Nelson, D., et al. 2020, MNRAS, submitted [arXiv:2008.00004] [Google Scholar]
Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891 [NASA ADS] [CrossRef] [Google Scholar]
D’Souza, R., Vegetti, S., & Kauffmann, G. 2015, MNRAS, 454, 4027 [Google Scholar]
Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2010, MNRAS, 409, 985 [Google Scholar]
Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662 [NASA ADS] [CrossRef] [Google Scholar]
Dubois, Y., Gavazzi, R., Peirani, S., & Silk, J. 2013, MNRAS, 433, 3297 [Google Scholar]
Dubois, Y., Pichon, C., Welker, C., et al. 2014a, MNRAS, 444, 1453 [NASA ADS] [CrossRef] [Google Scholar]
Dubois, Y., Volonteri, M., & Silk, J. 2014b, MNRAS, 440, 1590 [Google Scholar]
Dubois, Y., Volonteri, M., Silk, J., Devriendt, J., & Slyz, A. 2014c, MNRAS, 440, 2333 [Google Scholar]
Dubois, Y., Volonteri, M., Silk, J., et al. 2015, MNRAS, 452, 1502 [NASA ADS] [CrossRef] [Google Scholar]
Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948 [Google Scholar]
Dutton, A. A., Obreja, A., Wang, L., et al. 2017, MNRAS, 467, 4937 [Google Scholar]
Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Elmegreen, B. G., & Burkert, A. 2010, ApJ, 712, 294 [Google Scholar]
Elmegreen, B. G., & Elmegreen, D. M. 2005, ApJ, 627, 632 [Google Scholar]
Emsellem, E., Cappellari, M., Krajnović, D., et al. 2011, MNRAS, 414, 888 [Google Scholar]
Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, ApJ, 644, 813 [NASA ADS] [CrossRef] [Google Scholar]
Faucher-Giguère, C.-A. 2018, MNRAS, 473, 3717 [Google Scholar]
Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156 [NASA ADS] [CrossRef] [Google Scholar]
Ferrero, I., Navarro, J. F., Abadi, M. G., et al. 2017, MNRAS, 464, 4736 [Google Scholar]
Fontana, A., Dunlop, J. S., Paris, D., et al. 2014, A&A, 570, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fossati, M., Wilman, D. J., Mendel, J. T., et al. 2017, ApJ, 835, 153 [Google Scholar]
Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486 [Google Scholar]
Gallazzi, A., Bell, E. F., Zibetti, S., Brinchmann, J., & Kelson, D. D. 2014, ApJ, 788, 72 [Google Scholar]
Geen, S., Rosdahl, J., Blaizot, J., Devriendt, J., & Slyz, A. 2015, MNRAS, 448, 3248 [Google Scholar]
Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175 [Google Scholar]
Gentry, E. S., Krumholz, M. R., Madau, P., & Lupi, A. 2019, MNRAS, 483, 3647 [Google Scholar]
Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101 [Google Scholar]
Glowacki, M., Elson, E., & Davé, R. 2020, MNRAS, 498, 3687 [Google Scholar]
Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223 [NASA ADS] [CrossRef] [Google Scholar]
González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJ, 735, L34 [Google Scholar]
Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [CrossRef] [Google Scholar]
Haardt, F., & Madau, P. 1996, ApJ, 461, 20 [NASA ADS] [CrossRef] [Google Scholar]
Habouzit, M., Volonteri, M., & Dubois, Y. 2017, MNRAS, 468, 3935 [NASA ADS] [CrossRef] [Google Scholar]
Heiles, C., & Troland, T. H. 2005, ApJ, 624, 773 [NASA ADS] [CrossRef] [Google Scholar]
Hennebelle, P., & Chabrier, G. 2011, ApJ, 743, L29 [Google Scholar]
Hopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, 417, 950 [NASA ADS] [CrossRef] [Google Scholar]
Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522 [Google Scholar]
Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581 [Google Scholar]
Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800 [NASA ADS] [CrossRef] [Google Scholar]
Hughes, A., Meidt, S. E., Colombo, D., et al. 2013, ApJ, 779, 46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Iffrig, O., & Hennebelle, P. 2015, A&A, 576, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Inoue, S., Dekel, A., Mandelker, N., et al. 2016, MNRAS, 456, 2052 [NASA ADS] [CrossRef] [Google Scholar]
Islam, R. R., Taylor, J. E., & Silk, J. 2003, MNRAS, 340, 647 [NASA ADS] [CrossRef] [Google Scholar]
Iwamoto, K., Mazzali, P. A., Nomoto, K., et al. 1998, Nature, 395, 672 [Google Scholar]
Jackson, R. A., Martin, G., Kaviraj, S., et al. 2020, MNRAS, 494, 5568 [CrossRef] [Google Scholar]
Jackson, R. A., Kaviraj, S., Martin, G., et al. 2021a, MNRAS, 502, 1785 [CrossRef] [Google Scholar]
Jackson, R. A., Martin, G., Kaviraj, S., et al. 2021b, MNRAS, 502, 4262 [CrossRef] [Google Scholar]
Joung, M. R., Mac Low, M.-M., & Bryan, G. L. 2009, ApJ, 704, 137 [NASA ADS] [CrossRef] [Google Scholar]
Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L35 [NASA ADS] [CrossRef] [Google Scholar]
Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, ApJ, 758, 106 [NASA ADS] [CrossRef] [Google Scholar]
Kassin, S. A., Brooks, A., Governato, F., Weiner, B. J., & Gardner, J. P. 2014, ApJ, 790, 89 [NASA ADS] [CrossRef] [Google Scholar]
Kaviraj, S. 2014, MNRAS, 440, 2944 [Google Scholar]
Kaviraj, S., Schawinski, K., Devriendt, J. E. G., et al. 2007, ApJS, 173, 619 [NASA ADS] [CrossRef] [Google Scholar]
Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739 [NASA ADS] [Google Scholar]
Kaviraj, S., Martin, G., & Silk, J. 2019, MNRAS, 489, L12 [Google Scholar]
Kennicutt, R. C., Jr 1998, ApJ, 498, 541 [Google Scholar]
Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2 [Google Scholar]
Kim, C.-G., Ostriker, E. C., & Raileanu, R. 2017, ApJ, 834, 25 [NASA ADS] [CrossRef] [Google Scholar]
Kimm, T., & Cen, R. 2014, ApJ, 788, 121 [NASA ADS] [CrossRef] [Google Scholar]
Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015, MNRAS, 451, 2900 [NASA ADS] [CrossRef] [Google Scholar]
Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826 [NASA ADS] [Google Scholar]
Kimm, T., Haehnelt, M., Blaizot, J., et al. 2018, MNRAS, 475, 4617 [NASA ADS] [CrossRef] [Google Scholar]
King, A. R., Lubow, S. H., Ogilvie, G. I., & Pringle, J. E. 2005, MNRAS, 363, 49 [NASA ADS] [CrossRef] [Google Scholar]
Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145 [NASA ADS] [CrossRef] [Google Scholar]
Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [NASA ADS] [CrossRef] [Google Scholar]
Krajnović, D., Alatalo, K., Blitz, L., et al. 2013, MNRAS, 432, 1768 [Google Scholar]
Krajnović, D., Weilbacher, P. M., Urrutia, T., et al. 2015, MNRAS, 452, 2 [NASA ADS] [CrossRef] [Google Scholar]
Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, Astron. Lett., 44, 8 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M. R., & Burkhart, B. 2016, MNRAS, 458, 1671 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250 [NASA ADS] [CrossRef] [Google Scholar]
Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716 [NASA ADS] [CrossRef] [Google Scholar]
Lagos, C. D. P., Crain, R. A., Schaye, J., et al. 2015, MNRAS, 452, 3815 [CrossRef] [Google Scholar]
Lapiner, S., Dekel, A., & Dubois, Y. 2021, MNRAS, 505, 172 [CrossRef] [Google Scholar]
Lee, J., Shin, J., Snaith, O. N., et al. 2021, ApJ, 908, 11 [CrossRef] [Google Scholar]
Lelli, F., McGaugh, S. S., Schombert, J. M., Desmond, H., & Katz, H. 2019, MNRAS, 484, 3267 [NASA ADS] [CrossRef] [Google Scholar]
Lupi, A., Haardt, F., & Dotti, M. 2015, MNRAS, 446, 1765 [NASA ADS] [CrossRef] [Google Scholar]
Lupi, A., Bovino, S., Capelo, P. R., Volonteri, M., & Silk, J. 2018, MNRAS, 474, 2884 [CrossRef] [Google Scholar]
Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]
Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Mandelker, N., Dekel, A., Ceverino, D., et al. 2014, MNRAS, 443, 3675 [NASA ADS] [CrossRef] [Google Scholar]
Mannucci, F., Cresci, G., Maiolino, R., et al. 2009, MNRAS, 398, 1915 [NASA ADS] [CrossRef] [Google Scholar]
Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115 [Google Scholar]
Martin, G., Kaviraj, S., Devriendt, J. E. G., Dubois, Y., & Pichon, C. 2018a, MNRAS, 480, 2266 [NASA ADS] [CrossRef] [Google Scholar]
Martin, G., Kaviraj, S., Devriendt, J. E. G., et al. 2018b, MNRAS, 474, 3140 [CrossRef] [Google Scholar]
Martin, G., Kaviraj, S., Laigle, C., et al. 2019, MNRAS, 485, 796 [NASA ADS] [CrossRef] [Google Scholar]
Martin, G., Kaviraj, S., Hocking, A., Read, S. C., & Geach, J. E. 2020, MNRAS, 491, 1408 [Google Scholar]
Martin, G., Jackson, R. A., Kaviraj, S., et al. 2021, MNRAS, 500, 4937 [CrossRef] [Google Scholar]
Martín-Navarro, I., Vazdekis, A., La Barbera, F., et al. 2015, ApJ, 806, L31 [NASA ADS] [CrossRef] [Google Scholar]
Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364 [NASA ADS] [CrossRef] [Google Scholar]
McAlpine, S., Bower, R. G., Rosario, D. J., et al. 2018, MNRAS, 481, 3118 [NASA ADS] [CrossRef] [Google Scholar]
McGaugh, S. S. 2012, AJ, 143, 40 [Google Scholar]
McGaugh, S. S., & Schombert, J. M. 2015, ApJ, 802, 18 [CrossRef] [Google Scholar]
McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083 [Google Scholar]
Menezes, R. B., Steiner, J. E., & Ricci, T. V. 2014, ApJ, 796, L13 [CrossRef] [Google Scholar]
Menezes, R. B., Steiner, J. E., & da Silva, P. 2016, ApJ, 817, 150 [CrossRef] [Google Scholar]
Miller, B. P., Gallo, E., Greene, J. E., et al. 2015, ApJ, 799, 98 [CrossRef] [Google Scholar]
Moody, C. E., Romanowsky, A. J., Cox, T. J., Novak, G. S., & Primack, J. R. 2014, MNRAS, 444, 1475 [NASA ADS] [CrossRef] [Google Scholar]
Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121 [Google Scholar]
Mowla, L., van der Wel, A., van Dokkum, P., & Miller, T. B. 2019, ApJ, 872, L13 [NASA ADS] [CrossRef] [Google Scholar]
Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18 [Google Scholar]
Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59 [Google Scholar]
Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJ, 699, L178 [NASA ADS] [CrossRef] [Google Scholar]
Naab, T., Oser, L., Emsellem, E., et al. 2014, MNRAS, 444, 3357 [NASA ADS] [CrossRef] [Google Scholar]
Narayanan, D., Bothwell, M., & Davé, R. 2012, MNRAS, 426, 1178 [NASA ADS] [CrossRef] [Google Scholar]
Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2018, ApJ, 858, 118 [NASA ADS] [CrossRef] [Google Scholar]
Nickerson, S., Teyssier, R., & Rosdahl, J. 2019, MNRAS, 484, 1238 [Google Scholar]
Nomoto, K., Tominaga, N., Umeda, H., Kobayashi, C., & Maeda, K. 2006, Nucl. Phys. A, 777, 424 [NASA ADS] [CrossRef] [Google Scholar]
Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Nuñez-Castiñeyra, A., Nezri, E., Devriendt, J., & Teyssier, R. 2021, MNRAS, 501, 62 [Google Scholar]
Ocvirk, P., Pichon, C., & Teyssier, R. 2008, MNRAS, 390, 1326 [Google Scholar]
Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312 [NASA ADS] [CrossRef] [Google Scholar]
Ostriker, E. C. 1999, ApJ, 513, 252 [NASA ADS] [CrossRef] [Google Scholar]
Ostriker, E. C., & Shetty, R. 2011, ApJ, 731, 41 [NASA ADS] [CrossRef] [Google Scholar]
Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40 [Google Scholar]
Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2008, MNRAS, 391, 1117 [Google Scholar]
Park, M.-J., Yi, S. K., Dubois, Y., et al. 2019, ApJ, 883, 25 [Google Scholar]
Park, M. J., Yi, S. K., Peirani, S., et al. 2021, ApJS, 254, 2 [Google Scholar]
Peñarrubia, J., Navarro, J. F., & McConnachie, A. W. 2008, ApJ, 673, 226 [Google Scholar]
Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., & Bastian, N. 2018, MNRAS, 475, 4309 [NASA ADS] [CrossRef] [Google Scholar]
Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101 [CrossRef] [Google Scholar]
Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077 [Google Scholar]
Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196 [NASA ADS] [CrossRef] [Google Scholar]
Ponomareva, A. A., Verheijen, M. A. W., Papastergis, E., Bosma, A., & Peletier, R. F. 2018, MNRAS, 474, 4366 [Google Scholar]
Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464 [Google Scholar]
Popping, G., Somerville, R. S., & Trager, S. C. 2014, MNRAS, 442, 2398 [Google Scholar]
Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14 [NASA ADS] [CrossRef] [Google Scholar]
Pozzetti, L., Bolzonella, M., Lamareille, F., et al. 2007, A&A, 474, 443 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Prantzos, N., Abia, C., Limongi, M., Chieffi, A., & Cristallo, S. 2018, MNRAS, 476, 3432 [Google Scholar]
Prieto, J., Escala, A., Volonteri, M., & Dubois, Y. 2017, ApJ, 836, 216 [CrossRef] [Google Scholar]
Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82 [NASA ADS] [CrossRef] [Google Scholar]
Reines, A. E., Condon, J. J., Darling, J., & Greene, J. E. 2020, ApJ, 888, 36 [CrossRef] [Google Scholar]
Reynolds, C. S. 2013, Class. Quant. Grav., 30, 244004 [NASA ADS] [CrossRef] [Google Scholar]
Rezzolla, L., Barausse, E., Dorband, E. N., et al. 2008, Phys. Rev. D, 78, 044002 [CrossRef] [Google Scholar]
Rosdahl, J., & Blaizot, J. 2012, MNRAS, 423, 344 [NASA ADS] [CrossRef] [Google Scholar]
Rosdahl, J., & Teyssier, R. 2015, MNRAS, 449, 4380 [NASA ADS] [CrossRef] [Google Scholar]
Rosdahl, J., Schaye, J., Dubois, Y., Kimm, T., & Teyssier, R. 2017, MNRAS, 466, 11 [NASA ADS] [CrossRef] [Google Scholar]
Roškar, R., Teyssier, R., Agertz, O., Wetzstein, M., & Moore, B. 2014, MNRAS, 444, 2837 [CrossRef] [Google Scholar]
Saftly, W., Baes, M., De Geyter, G., et al. 2015, A&A, 576, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312 [NASA ADS] [CrossRef] [Google Scholar]
Sales, L. V., Navarro, J. F., Oman, K., et al. 2017, MNRAS, 464, 2419 [CrossRef] [Google Scholar]
Sánchez, S. F., Barrera-Ballesteros, J. K., López-Cobá, C., et al. 2019, MNRAS, 484, 3042 [CrossRef] [Google Scholar]
Sandage, A., Tammann, G. A., & Yahil, A. 1979, ApJ, 232, 352 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2018, ApJ, 858, 99 [NASA ADS] [CrossRef] [Google Scholar]
Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521 [Google Scholar]
Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150 [NASA ADS] [CrossRef] [Google Scholar]
Sedgwick, T. M., Baldry, I. K., James, P. A., & Kelvin, L. S. 2019a, IAU Symposium 355: The Realm of the Low Surface Brightness Universe [Google Scholar]
Sedgwick, T. M., Baldry, I. K., James, P. A., & Kelvin, L. S. 2019b, MNRAS, 484, 5278 [Google Scholar]
Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2018, ApJ, 861, 4 [NASA ADS] [CrossRef] [Google Scholar]
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
She, R., Ho, L. C., & Feng, H. 2017, ApJ, 842, 131 [CrossRef] [Google Scholar]
Shen, Y., Hwang, H.-C., Zakamska, N., & Liu, X. 2019, ApJ, 885, L4 [CrossRef] [Google Scholar]
Silk, J. 2017, ApJ, 839, L13 [NASA ADS] [CrossRef] [Google Scholar]
Silk, J., & Rees, M. J. 1998, A&A, 331, L1 [NASA ADS] [Google Scholar]
Simons, R. C., Kassin, S. A., Weiner, B. J., et al. 2015, MNRAS, 452, 986 [NASA ADS] [CrossRef] [Google Scholar]
Simons, R. C., Kassin, S. A., Weiner, B. J., et al. 2017, ApJ, 843, 46 [NASA ADS] [CrossRef] [Google Scholar]
Smith, R., Choi, H., Lee, J., et al. 2016, ApJ, 833, 109 [NASA ADS] [CrossRef] [Google Scholar]
Smith, A., Bromm, V., & Loeb, A. 2017, MNRAS, 464, 2963 [NASA ADS] [CrossRef] [Google Scholar]
Soares, G., & Nemmen, R. 2020, MNRAS, 495, 981 [CrossRef] [Google Scholar]
Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5 [NASA ADS] [CrossRef] [Google Scholar]
Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [Google Scholar]
Sun, J., Leroy, A. K., Schruba, A., et al. 2018, ApJ, 860, 172 [NASA ADS] [CrossRef] [Google Scholar]
Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 [NASA ADS] [CrossRef] [Google Scholar]
Tacconi, L. J., Neri, R., Chapman, S. C., et al. 2006, ApJ, 640, 228 [Google Scholar]
Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246 [Google Scholar]
Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781 [Google Scholar]
Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74 [Google Scholar]
Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179 [NASA ADS] [CrossRef] [Google Scholar]
Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157 [CrossRef] [Google Scholar]
Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Teyssier, R., Moore, B., Martizzi, D., Dubois, Y., & Mayer, L. 2011, MNRAS, 414, 195 [CrossRef] [Google Scholar]
Thorne, K. S. 1974, ApJ, 191, 507 [Google Scholar]
Thornton, K., Gaudlitz, M., Janka, H. T., & Steinmetz, M. 1998, ApJ, 500, 95 [NASA ADS] [CrossRef] [Google Scholar]
Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2016, ApJ, 817, 118 [NASA ADS] [CrossRef] [Google Scholar]
Toomre, A., & Toomre, J. 1972, ApJ, 178, 623 [Google Scholar]
Toro, E. 1999, Riemann Solvers and Numerical Methods for Fluid Dynamics (Berlin: Springer-Verlag) [CrossRef] [Google Scholar]
Trebitsch, M., Blaizot, J., Rosdahl, J., Devriendt, J., & Slyz, A. 2017, MNRAS, 470, 224 [NASA ADS] [CrossRef] [Google Scholar]
Trebitsch, M., Volonteri, M., Dubois, Y., & Madau, P. 2018, MNRAS, 478, 5607 [CrossRef] [Google Scholar]
Trebitsch, M., Dubois, Y., Volonteri, M., et al. 2020, A&A, submitted [arXiv:2002.04045] [Google Scholar]
Tremmel, M., Karcher, M., Governato, F., et al. 2017, MNRAS, 470, 1121 [CrossRef] [Google Scholar]
Tremmel, M., Governato, F., Volonteri, M., Pontzen, A., & Quinn, T. R. 2018, ApJ, 857, L22 [Google Scholar]
Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898 [NASA ADS] [CrossRef] [Google Scholar]
Treu, T., Auger, M. W., Koopmans, L. V. E., et al. 2010, ApJ, 709, 1195 [Google Scholar]
Trujillo, I., Förster Schreiber, N. M., Rudnick, G., et al. 2006, ApJ, 650, 18 [NASA ADS] [CrossRef] [Google Scholar]
Tully, R. B., & Fisher, J. R. 1977, A&A, 500, 105 [Google Scholar]
Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, MNRAS, 423, 2374 [NASA ADS] [CrossRef] [Google Scholar]
van de Sande, J., Bland-Hawthorn, J., Brough, S., et al. 2017, MNRAS, 472, 1272 [NASA ADS] [CrossRef] [Google Scholar]
van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28 [NASA ADS] [CrossRef] [Google Scholar]
Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177 [NASA ADS] [CrossRef] [Google Scholar]
Volonteri, M., Lodato, G., & Natarajan, P. 2008, MNRAS, 383, 1079 [NASA ADS] [CrossRef] [Google Scholar]
Volonteri, M., Dubois, Y., Pichon, C., & Devriendt, J. 2016, MNRAS, 460, 2979 [Google Scholar]
Volonteri, M., Pfister, H., Beckmann, R. S., et al. 2020, MNRAS, 498, 2219 [Google Scholar]
Weiner, B. J., Willmer, C. N. A., Faber, S. M., et al. 2006, ApJ, 653, 1027 [NASA ADS] [CrossRef] [Google Scholar]
Wetzel, A. R., Tinker, J. L., Conroy, C., & van den Bosch, F. C. 2013, MNRAS, 432, 336 [NASA ADS] [CrossRef] [Google Scholar]
Yabe, K., Ohta, K., Akiyama, M., et al. 2015, PASJ, 67, 102 [CrossRef] [Google Scholar]
Yu, S., Bullock, J. S., Wetzel, A., et al. 2020, MNRAS, 494, 1539 [CrossRef] [Google Scholar]
Zahid, H. J., Kewley, L. J., & Bresolin, F. 2011, ApJ, 730, 137 [NASA ADS] [CrossRef] [Google Scholar]
Zahid, H. J., Kashino, D., Silverman, J. D., et al. 2014, ApJ, 792, 75 [NASA ADS] [CrossRef] [Google Scholar]
Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Purity of halos

Figure A.1 shows the halo DM mass function at redshift z = 1 for different levels of purity of the host: for 90%, 99%, or 100% of purity, where the levels of purity are computed in terms of the number fraction of high-resolution DM particles over the total number of DM particles in the halo. The amount of massive halos can increase by up to 30% if halos other than perfectly pure halos are considered. However, in this study we only consider halos with a 100% purity.

Fig. A.1.

Halo mass function (bottom panels) and fraction of non-pure halos (top panels) within the zoom simulation for halos with different levels of pollution (as indicated in the panel) for redshift z = 1. Only halos with a 100% purity are considered in this study.

Appendix B: Comparison of HOP vs. AdaptaHOP sizes

To illustrate how well the AdaptaHOP can separate galaxy substructures, in Fig. B.1 we show two stellar density maps of a galaxy at z = 1 (the most massive galaxy), extracted from the identified structures with HOP (including the stellar clumps) or with AdaptaHOP (which removes sub-structures). With the HOP finder, multiple stellar clumps can be identified within the galaxy, while AdaptaHOP has removed all substructures. In Fig. B.2, we show two stellar density maps of a merging galaxy at z = 0.4 as extracted from either HOP or AdaptaHOP. Indeed, AdaptaHOP efficiently identifies stellar clumps and removes them from the main structures and also removes companion satellite galaxies⁹ when their stellar distribution connects to the main object. Since HOP is not able to remove any substructure, evaluating the size of a galaxy becomes a severe issue for galaxies in a significant interaction; see the two bottom panels of Fig. B.2, which leads to an incorrect evaluation of the size of the main galaxy.

Fig. B.1.

Stellar density distribution for the most massive galaxy M_s = 1.5 × 10¹¹ M_⊙ at z = 1 as identified by HOP (left panel) or AdaptaHOP (right panel). Image size is 30 kpc across.

Fig. B.2.

Stellar density maps of a galaxy and companion satellite at z = 0.4 as obtained by HOP (top left panel) and AdaptaHOP (top right panel). Stellar mass and effective radius are M_s = 2.14 × 10¹¹ M_⊙ and R_eff = 39 kpc, and M_s = 1.25 × 10¹¹ M_⊙ and R_eff = 3.3 kpc for HOP and AdaptaHOP respectively. The resulting size-mass relations (plotted as in Fig. 12) are affected by the galaxy finder as shown in the two bottom panels (left: HOP; right: AdaptaHOP).

Appendix C: Fraction of morphological types: Changing the ad hoc thresholds for classification

The exact value of the fraction of galaxies per morphological type depends on the adopted thresholds for their kinematic V/σ or their asymmetry A_r classifier. The adopted values of (V/σ)_c = 0.5 and A_r, c = 0.3 were chosen to best fit the fractions from Conselice (2006). We show in Fig. C.1 how changing those values between (V/σ)_c = 0.3 − 0.7 and A_r, c = 0.2 − 0.4 affect the fraction of each morphological type. These fractions change with adopting different values for those thresholds, as expected. There are a higher (lower) fraction of ellipticals and irregulars when increasing (decreasing) (V/σ)_c (and conversely for the fraction of kinematically classified discs) and decreasing (increasing) A_r, c, but the qualitative trends with mass are preserved.

Fig. C.1.

Same as Fig. 16: Fraction of each morphological type as a function of stellar mass with different thresholds for (V/σ)_c and A_r, c, with 0.7 and 0.2 (top panel), and 0.3 and 0.4 (bottom panel), respectively.

Appendix D: Fraction of low-mass stellar clusters

The minimum number of stellar particles used to detect substructures with AdaptaHOP can affect the amount of stellar mass contained in stellar clusters. We evaluate the effect by lowering our fiducial number of 50 stellar particles down to 10 on the fraction of stellar mass contained in clusters in Fig. D.1. The mean relation are relatively similar at all considered redshifts, except for the most massive galaxies M_s ≥ 10¹¹ M_⊙ at the lowest redshift z = 0.25 considered (to be compared with Fig. 17).

Fig. D.1.

Same as Fig. 17: fraction of stellar mass in stellar clusters as a function of galaxy stellar mass at different redshifts as indicated in the panel, except that a lower number of stellar particles are used to detect stellar clumps (10 instead of 50) in AdaptaHOP. The solid lines indicate the mean values, and dashed lines stand for the error on the mean. The data points are overplotted for the most extreme redshifts to appreciate the scatter in the distribution of values.

All Tables

Table 1.

Number of galaxies for different stellar mass thresholds.

In the text

All Figures

Fig. 1.

Sequential zoom (clockwise from top left) over the projected density (silver blue colours) and projected temperature (red) of the NEWHORIZON simulation at redshift z = 2. The dashed white circles encompass the initial high-resolution volume. Each panel is a zoomed-in version of the previous panel (identified by the white square in the previous panel) with the panel sizes of 142, 18, 4.4, and 1.1 comoving Mpc width, respectively. Two top panels: encompass the zoom-in region, with its network of filaments. Two bottom panels: how narrow filaments break up and mix once they connect to one of the most massive galaxies of that zoom-in region.

In the text

	Fig. 2. Illustration of the structure of the gas density (left panels) and the corresponding spatial resolution (right panels) in a massive galaxy of M_s = 6 × 10¹⁰ M_⊙ at z = 1 seen edge-on (top panels) or face-on (bottom panels).
In the text

Fig. 3.

Spin-up rate (top panel) and radiative efficiency ϵ_rad/f_att (bottom panel) as a function of the MBH spin for the thin disc solution (Shakura & Sunyaev 1973) applied to the quasar mode. At negative values of the MBH spin, the gas accreted from the thin accretion disc decreases the MBH spin, while for a positive MBH spin, the gas increases. For the thin disc solution, the radiative efficiency is an increasing function of the MBH spin with a sharp increase (by 4) between a MBH spin of 0.7 and 0.998.

In the text

Fig. 4.

Spin-up rate (top panel), jet (red plus signs), wind (blue diamonds), and total (black) efficiencies (bottom panel) as a function of the MBH spin for the MCAD solution applied to the radio mode. The symbols represent the results from the simulations of McKinney et al. (2012) and the solid lines indicate the interpolated functions used in the NEWHORIZON simulation. As opposed to the thin disc solution (Fig. 3), gas accreted from the thick accretion disc always decreases the MBH spin. The MCAD feedback efficiency is an increasing function of the absolute value of the BH spin with a minimum efficiency for non-rotating MBHs.

In the text

Fig. 5.

Projection of the gas density and mock observations. Panels on the left are images of a massive galaxy at different epochs and the panels on the right are for a less massive galaxy. First row: gas density projections a of 1.4 Mpc at different epochs, while the second and fourth rows are zoomed-in gas density projections with face-on and edge-on views of the galaxy, respectively. Third and fifth rows are SKIRT mock observations in face-on and edge-on direction. Stellar mass and halo mass of each galaxy at each epoch are given in the second row in log scale. For a galaxy at a given epoch, the second to fifth panels are on the same scale and the white bar in third row indicates 5 kpc. Gas density maps share the same colour scheme as given in the colour bar.

In the text

Fig. 6.

Galaxy stellar mass function at z = 0.25 in NEWHORIZON and the HSC-SSP survey. The light blue squares and dark blue circles with error bars indicate the NEWHORIZON stellar mass function for all galaxies (offset by 0.15 dex) and only uncontaminated galaxies (with a volume correction applied). The black line indicates the median galaxy stellar mass function from 100 random pointings from a volume of the HSC-SSP deep layer with the same volume as NEWHORIZON and the black error bar indicates the 90% scatter in the mass function at the low-mass end. A comparison is made with additional mass functions from literature (Sedgwick et al. 2019a; Davidzon et al. 2017; Tomczak et al. 2016; Song et al. 2016; Grazian et al. 2015; D’Souza et al. 2015; Fontana et al. 2014; Bernardi et al. 2013; Ilbert et al. 2013; Baldry et al. 2012; Bielby et al. 2012; González et al. 2011; Pozzetti et al. 2007), which are shown as thin coloured lines. The thick purple dotted line indicates the NEWHORIZON mass function constructed using the Sandage et al. (1979) method (STY) and including selection effects.

In the text

Fig. 7.

Surface brightness vs. stellar mass in the NEWHORIZON simulation for 3 redshifts. The grey points indicate the entire galaxy population of NEWHORIZON in all 3 panels with the median lines for the redshift in question and z = 0.25 shown in blue and orange, respectively. In the bottom panel the open blue points indicate galaxies from Sedgwick et al. (2019a) and black points are NEWHORIZON galaxies. The predicted surface brightness vs. stellar mass plane in NEWHORIZON corresponds well to that where the Sedgwick et al. galaxies are complete. The red dotted lines in the bottom panel indicate the 70% and 10% completeness limits from the SDSS (see e.g., Table 1 in Blanton et al. 2005). The overwhelming majority of galaxies in the Universe lie below the surface brightness thresholds of surveys such as the SDSS; only those galaxies that depart strongly from the typical surface brightness vs. stellar mass relation are likely to be detectable in these datasets.

In the text

Fig. 8.

Cosmic SFR density (top panel) and stellar density (bottom panel) as a function of redshift in NEWHORIZON (solid black or red lines) compared to observations (Behroozi et al. 2013, B13; Madau & Dickinson 2014; Novak et al. 2017; Driver et al. 2018) as indicated in the panels. All observational quantities are shown for a Chabrier (2005) IMF except for the dashed blue line, which indicates the Madau & Dickinson (2014) fit for a Salpeter IMF as originally assumed in their analysis. For the stellar density, the reconstructed result from the SFR density for NEWHORIZON (magenta line) and for the fit from B13 (grey lines) using two different return fractions R are also shown. The red coloured lines indicate the cosmic SFR densities of the different galaxy stellar mass bins as indicated in the panel in log M_⊙ units. The large error bar in black in the SFR density panel corresponds to the estimate of the cosmic variance (see text for details). The cosmic SFR in NEWHORIZON shows the qualitative expected trend over time; however, there is a systematic offset by a factor 1.5 − 2 with respect to the observational data assuming standard Chabrier-like IMF. The cosmic stellar density in NEWHORIZON is broad agreement with the data with a slight overestimate at low redshift, although there is an important uncertainty on the simulated cosmic SFR (and thus stellar density) due to the large cosmic variance associated with the simulated volume.

In the text

Fig. 9.

sSFR as a function of galaxy stellar mass at different redshifts (as indicated in the panel) in NEWHORIZON with solid lines. The error bars stand for the error around the mean. The symbols correspond to the best fit from Behroozi et al. (2013) from their collection of observational data; their uncertainty is shown in dashed lines. The simulated sSFR show very little evolution with stellar mass at the two highest redshifts, while there is a significant quenching at the lowest redshifts for stellar masses above M_s > 10¹⁰ M_⊙. Simulated sSFR show a fair level of agreement with the observations.

In the text

Fig. 10.

Kennicutt–Schmidt relation compared to observations. Filled coloured circles correspond to NEWHORIZON galaxies at various redshifts. The mean and error on the mean are shown with solid and dotted coloured lines, respectively. The empty black symbols correspond to observational datasets: triangles are local spirals (Kennicutt 1998), stars are high-z BzK/Normal galaxies (Daddi et al. 2010a; Tacconi et al. 2010), circles high-z mergers (Bouché et al. 2007; Bothwell et al. 2009), and squares low-z mergers (Kennicutt 1998). The solid black and dotted black lines represent the sequence of discs and starbursts, respectively, from Daddi et al. (2010b). The NEWHORIZON tool is able to capture the galaxy main sequence, although the SFR surface densities are slightly lower than in the observational data. There is a larger fraction of high SFR galaxies at higher redshift in qualitative agreement with observations.

In the text

Fig. 11.

Stellar-to-halo mass relation (top) and baryon conversion efficiency M_s/M_h at different redshifts (as indicated in the panels). The solid red lines represent the average with their error of the mean (dashed) with individual points as plus symbols. The solid black line indicates the semi-empirical relation from Moster et al. (2013) and the green line from Behroozi et al. (2013) at the indicated redshift and dot-dashed lines the respective relations at z = 0. The cyan lines stand for the constant star formation conversion efficiencies. The dotted lines correspond to the minimum stellar mass detected by the galaxy finder. The simulated relation between central stellar mass and halo mass is in fair qualitative agreement with the semi-empirical relations, where the increase in the baryon conversion efficiency with mass up to the peak at Milky Way halo mass is captured; however this efficiency is significantly overestimated below that peak.

In the text

Fig. 12.

Effective radius as a function of stellar mass in NEWHORIZON at different redshifts for individual galaxies (grey crosses) and for the average (solid line) with standard deviation (dashed lines). The sizes are computed with the geometric mean of the x, y, z effective radius using the AdaptaHOP classification. Observations from Mowla et al. (2019) are shown as a purple line at the current redshift when available and with their most extreme redshift fits at z = 0.37 and z = 2.69 in blue and red, respectively, which are extrapolated (triple dot-dashed line) beyond their range of available data (solid line) to guide the eye. Galaxies in NEWHORIZON are more compact at high redshift as opposed to low redshift in good agreement with observations.

In the text

Fig. 13.

Average stellar (coloured solid) and gas (coloured dashed) metallicity as a function of stellar mass for different redshifts in NEWHORIZON as indicated in the panel. The error bars stand for the error around the mean. Observational fits for the stellar metallicity from Gallazzi et al. (2014; solid black at z = 0.1, and dashed black at z = 0.7) and from Panter et al. (2008; solid grey) are also shown. The metallicity of both the gas and stellar component are decreasing with increasing redshift owing to less chemically enriched galaxies.

In the text

Fig. 14.

Average abundance of oxygen (coloured solid) in the cold gas as a function of stellar mass for different redshifts in NEWHORIZON as indicated in the panel. The error bars stand for the error around the mean. Several observational relations are also shown from Mannucci et al. (2009; M09; in dotted lines with four different colours corresponding to z = 0.07, 0.7, 2.2, and 3.5), Zahid et al. (2014; Z14; dot-dashed line), Sánchez et al. (2019; S19; dashed lines with different lines corresponding to various oxygen calibrators), and Curti et al. (2020; C20; black solid).

In the text

Fig. 15.

Ratio of stellar rotation over dispersion as a function of galaxy stellar mass for different redshifts (colours) as indicated in the panels. The kinematics are measured within R_eff (top) or 2R_eff (bottom). The solid lines indicate the average and the dashed lines are the error around the mean, with individual points shown at z = 0.25 only (purple plus signs) to appreciate the scatter of the distribution. Galaxies show an increase support of stellar rotation over dispersion with time and galaxy mass (except for z = 4) with stars outside R_eff having more rotation than inside.

In the text

Fig. 16.

Fraction of galaxies classified by their morphological type: f_{ell + irr} for elliptical+irregulars (red solid line; galaxies with V/σ < 0.5) in NEWHORIZON based on the stellar kinematics measured within R_eff (top) or 2R_eff (bottom) and as a function of stellar mass at z = 0.25; conversely, the fraction of discs f_disc = 1 − f_{ell + irr} is compared to data from Conselice (2006), which is at z ≃ 0, and is represented by symbols with error bars. The fraction of irregulars in NEWHORIZON inferred through the asymmetry index, within the same aperture in both panels, is also shown as the magenta curve. The error bars and dashed line stand for the error on the mean. The trends of the morphological fractions with mass in NEWHORIZON is qualitatively consistent with observational data.

In the text

Fig. 17.

Fraction of stellar mass in stellar clusters as a function of galaxy stellar mass at different redshifts as indicated in the panel. The solid lines indicate the mean values, and the dashed lines stand for the error on the mean. The data points are overplotted for the most extreme redshifts to appreciate the scatter in the distribution of sampled values. Galaxies at high redshift are more clumpy than at low redshift at fixed stellar mass.

In the text

Fig. 18.

Surface densities within 2R_eff as a function of galaxy stellar mass at different redshifts for stars (top), HI+H₂ gas (n > 0.1 cm⁻³ and T < 2 × 10⁴ K, middle), and H₂ molecular gas (n > 10 cm⁻³ and T < 2 × 10⁴ K, bottom). The solid lines represent the mean values, and the dashed lines stand for the error on the mean. Galaxies are denser at high redshift at fixed stellar mass and denser in more massive galaxies at fixed redshift.

In the text

Fig. 19.

Ratio of gas over baryon content (top panels) within 2R_eff as a function of galaxy stellar mass at different redshifts for the cold dense star-forming gas with density higher than 0.1 or 10 cm⁻³ (left and right, respectively) and temperature lower than 2 × 10⁴ K, and the ratio of this cold selected gas over total gas (bottom panels). The solid lines indicate the mean values, and the dashed lines stand for the error on the mean. The points with error bars correspond to the observations from Santini et al. (2014; sampled from their interpolated curves in their Fig. 12) at various redshifts and from Catinella et al. (2018) from local galaxies for the HI+H₂ gas, while the dotted line shows the data from Tacconi et al. (2018) at various redshifts for the H₂ molecular gas only. The NEWHORIZON galaxies are less gas-rich at the high-mass end and also at lower redshift for a given stellar mass, in qualitative agreement with the data. But they show significantly less gas at the low-mass end compared to observations, and especially at high redshift with the data from Tacconi et al. (2018), and in contrast they are in good agreement with the data from Santini et al. (2014).

In the text

Fig. 20.

Baryonic Tully–Fisher relation for disc galaxies in NEWHORIZON at different redshifts as labelled in the panel with the raw data points overplotted as crosses for redshift z = 0.25 (purple plus signs), compared to the two best fits from McGaugh & Schombert (2015; dashed black and dot-dashed black) and the best fit from Ponomareva et al. (2018; solid black; in thick their fit over the range of their data points, while the thin part is the extrapolated part). The fit to the NEWHORIZON data at z = 0.25 (thick dotted purple line) shows a better level of agreement with the data from Ponomareva et al. (2018; with a slope of one-third) over those of McGaugh & Schombert (2015; slope of one-fourth).

In the text

Fig. 21.

Gas kinematics as a function of galaxy stellar mass at different redshifts for gas density 0.01 < n < 10 cm⁻³, temperature lower than T < 2 × 10⁴ K and within 2R_eff. From left to right and top to bottom: rotational velocity, velocity dispersion, velocity dispersion over S_g, 0.5, and rotational velocity over S_g, 0.5. The error bars stand for the error on the mean. The dashed plus signs in the top right panel indicate the values with uncertainties from the fit to the observational data in Simons et al. (2017). The dotted lines on the bottom panels stand for the cases of pure rotational support or pure dispersion support, respectively. The NEWHORIZON galaxies decrease their amount of gas velocity dispersion over time for a given stellar mass, while keeping a similar amount of rotation.

In the text

Fig. 22.

MBH mass as a function of galaxy mass for redshifts z = 4, 2, 1 and 0.25 as circles for all MBHs contained within 2 effective radii of their host galaxy. For the lowest redshift, z = 0.25, big circles (dark red) highlight the most massive MBH for a given galaxy, while the small circles (dark red) show all the secondary MBHs within the galaxy. The black line denotes the mean M_BH for all primary BH within a given M_S bin at z = 0.25 for NEWHORIZON. Also shown are observations of MBH vs. stellar mass for z ∼ 0 − 0.3 from Reines & Volonteri (2015; RV15; green triangles), Greene et al. (2020; Greene20; blue markers) and Baron & Ménard (2019; BM19; grey contours). The error bars in RV15 were omitted for clarity. On average, central MBHs in NEWHORIZON grow significantly only above a stellar mass threshold of a few 10⁹ M_⊙, with non-central MBHs generally failing to grow even in galaxies above this mass threshold.

In the text

	Fig. 23. MBH spin evolution with MBH mass for all MBH at redshift z = 0.25 (black markers), and the spin evolution history with mass for the six most massive MBH (coloured lines). Shown for comparison are observations from Reynolds (2013; grey triangles) and Soares & Nemmen (2020; grey crosses). MBHs experience a fast initial spin growth phase followed by a more stochastic evolution.
In the text

	Fig. A.1. Halo mass function (bottom panels) and fraction of non-pure halos (top panels) within the zoom simulation for halos with different levels of pollution (as indicated in the panel) for redshift z = 1. Only halos with a 100% purity are considered in this study.
In the text

	Fig. B.1. Stellar density distribution for the most massive galaxy M_s = 1.5 × 10¹¹ M_⊙ at z = 1 as identified by HOP (left panel) or AdaptaHOP (right panel). Image size is 30 kpc across.
In the text

Fig. B.2.

Stellar density maps of a galaxy and companion satellite at z = 0.4 as obtained by HOP (top left panel) and AdaptaHOP (top right panel). Stellar mass and effective radius are M_s = 2.14 × 10¹¹ M_⊙ and R_eff = 39 kpc, and M_s = 1.25 × 10¹¹ M_⊙ and R_eff = 3.3 kpc for HOP and AdaptaHOP respectively. The resulting size-mass relations (plotted as in Fig. 12) are affected by the galaxy finder as shown in the two bottom panels (left: HOP; right: AdaptaHOP).

In the text

	Fig. C.1. Same as Fig. 16: Fraction of each morphological type as a function of stellar mass with different thresholds for (V/σ)_c and A_r, c, with 0.7 and 0.2 (top panel), and 0.3 and 0.4 (bottom panel), respectively.
In the text

Fig. D.1.

Same as Fig. 17: fraction of stellar mass in stellar clusters as a function of galaxy stellar mass at different redshifts as indicated in the panel, except that a lower number of stellar particles are used to detect stellar clumps (10 instead of 50) in AdaptaHOP. The solid lines indicate the mean values, and dashed lines stand for the error on the mean. The data points are overplotted for the most extreme redshifts to appreciate the scatter in the distribution of values.

In the text

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

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

Initial download of the metrics may take a while.

[1] Agertz, O., Lake, G., Teyssier, R., et al. 2009, MNRAS, 392, 294 [NASA ADS] [CrossRef] [Google Scholar]

[2] Agertz, O., Teyssier, R., & Moore, B. 2011, MNRAS, 410, 1391 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[3] Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25 [NASA ADS] [CrossRef] [Google Scholar]

[4] Agertz, O., Renaud, F., Feltzing, S., et al. 2021, MNRAS, 503, 5826 [CrossRef] [Google Scholar]

[5] Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, PASJ, 71, 114 [CrossRef] [Google Scholar]

[6] Algorry, D. G., Navarro, J. F., Abadi, M. G., et al. 2014, MNRAS, 437, 3596 [NASA ADS] [CrossRef] [Google Scholar]

[7] Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[8] Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]

[9] Andersson, E. P., Agertz, O., & Renaud, F. 2020, MNRAS, 494, 3328 [CrossRef] [Google Scholar]

[10] André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[11] Anglés-Alcázar, D., Faucher-Giguère, C.-A., Quataert, E., et al. 2017, MNRAS, 472, L109 [NASA ADS] [CrossRef] [Google Scholar]

[12] Arnouts, S., & Ilbert, O. 2011, Astrophysics Source Code Library [record ascl:1108.009] [Google Scholar]

[13] Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540 [NASA ADS] [CrossRef] [Google Scholar]

[14] Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]

[15] Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376 [NASA ADS] [CrossRef] [Google Scholar]

[16] Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621 [NASA ADS] [Google Scholar]

[17] Baron, D., & Ménard, B. 2019, MNRAS, 487, 3404 [CrossRef] [Google Scholar]

[18] Baumgardt, H., & Makino, J. 2003, MNRAS, 340, 227 [NASA ADS] [CrossRef] [Google Scholar]

[19] Beckmann, R. S., Devriendt, J., Slyz, A., et al. 2017, MNRAS, 472, 949 [NASA ADS] [CrossRef] [Google Scholar]

[20] Beckmann, R. S., Slyz, A., & Devriendt, J. 2018, MNRAS, 478, 995 [CrossRef] [Google Scholar]

[21] Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [Google Scholar]

[22] Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143 [Google Scholar]

[23] Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212 [NASA ADS] [CrossRef] [Google Scholar]

[24] Bellovary, J. M., Governato, F., Quinn, T. R., et al. 2010, ApJ, 721, L148 [NASA ADS] [CrossRef] [Google Scholar]

[25] Bellovary, J. M., Cleary, C. E., Munshi, F., et al. 2019, MNRAS, 482, 2913 [NASA ADS] [Google Scholar]

[26] Benson, A. J., & Babul, A. 2009, MNRAS, 397, 1302 [CrossRef] [Google Scholar]

[27] Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, MNRAS, 436, 697 [Google Scholar]

[28] Bernstein, G. M., & Jarvis, M. 2002, AJ, 123, 583 [Google Scholar]

[29] Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[30] Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140 [NASA ADS] [CrossRef] [Google Scholar]

[31] Béthermin, M., Wang, L., Doré, O., et al. 2013, A&A, 557, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[32] Bielby, R., Hudelot, P., McCracken, H. J., et al. 2012, A&A, 545, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[33] Bieri, R., Dubois, Y., Rosdahl, J., et al. 2017, MNRAS, 464, 1854 [NASA ADS] [CrossRef] [Google Scholar]

[34] Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [Google Scholar]

[35] Blanton, M. R., Lupton, R. H., Schlegel, D. J., et al. 2005, ApJ, 631, 208 [Google Scholar]

[36] Blondin, J. M., Wright, E. B., Borkowski, K. J., & Reynolds, S. P. 1998, ApJ, 500, 342 [NASA ADS] [CrossRef] [Google Scholar]

[37] Bois, M., Emsellem, E., Bournaud, F., et al. 2011, MNRAS, 416, 1654 [NASA ADS] [CrossRef] [Google Scholar]

[38] Bothwell, M. S., Kennicutt, R. C., & Lee, J. C. 2009, MNRAS, 400, 154 [NASA ADS] [CrossRef] [Google Scholar]

[39] Bouché, N., Cresci, G., Davies, R., et al. 2007, ApJ, 671, 303 [NASA ADS] [CrossRef] [Google Scholar]

[40] Bournaud, F., Jog, C. J., & Combes, F. 2007, A&A, 476, 1179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[41] Bournaud, F., Elmegreen, B. G., Teyssier, R., Block, D. L., & Puerari, I. 2010, MNRAS, 409, 1088 [NASA ADS] [CrossRef] [Google Scholar]

[42] Bournaud, F., Perret, V., Renaud, F., et al. 2014, ApJ, 780, 57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[43] Bower, R. G., Schaye, J., Frenk, C. S., et al. 2017, MNRAS, 465, 32 [NASA ADS] [CrossRef] [Google Scholar]

[44] Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40 [CrossRef] [Google Scholar]

[45] Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]

[46] Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 [NASA ADS] [CrossRef] [Google Scholar]

[47] Bustamante, S., & Springel, V. 2019, MNRAS, 490, 4133 [CrossRef] [Google Scholar]

[48] Cañas, R., Elahi, P. J., Welker, C., et al. 2019, MNRAS, 482, 2039 [CrossRef] [Google Scholar]

[49] Cacciato, M., Dekel, A., & Genel, S. 2012, MNRAS, 421, 818 [NASA ADS] [Google Scholar]

[50] Camps, P., & Baes, M. 2020, Astron. Comput., 31, 100381 [CrossRef] [Google Scholar]

[51] Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2012, Nature, 484, 485 [NASA ADS] [CrossRef] [Google Scholar]

[52] Catinella, B., Saintonge, A., Janowiecki, S., et al. 2018, MNRAS, 476, 875 [NASA ADS] [CrossRef] [Google Scholar]

[53] Ceverino, D., & Klypin, A. 2009, ApJ, 695, 292 [Google Scholar]

[54] Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 404, 2151 [NASA ADS] [Google Scholar]

[55] Ceverino, D., Primack, J., Dekel, A., & Kassin, S. A. 2017, MNRAS, 467, 2664 [NASA ADS] [CrossRef] [Google Scholar]

[56] Chabanier, S., Bournaud, F., Dubois, Y., et al. 2020a, A&A, 643, L8 [CrossRef] [EDP Sciences] [Google Scholar]

[57] Chabanier, S., Bournaud, F., Dubois, Y., et al. 2020b, MNRAS, 495, 1825 [CrossRef] [Google Scholar]

[58] Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]

[59] Chabrier, G. 2005, in The Initial Mass Function: From Salpeter 1955 to 2005, eds. E. Corbelli, F. Palla, & H. Zinnecker, Astrophys. Space Sci. Lib., 327, 41 [Google Scholar]

[60] Chevalier, R. A. 1974, ApJ, 188, 501 [NASA ADS] [CrossRef] [Google Scholar]

[61] Chiosi, C., Bertelli, G., & Bressan, A. 1992, ARA&A, 30, 235 [NASA ADS] [CrossRef] [Google Scholar]

[62] Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252 [NASA ADS] [CrossRef] [Google Scholar]

[63] Coccato, L., Fabricius, M., Morelli, L., et al. 2015, A&A, 581, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[64] Conselice, C. J. 2006, MNRAS, 373, 1389 [NASA ADS] [CrossRef] [Google Scholar]

[65] Conselice, C. J. 2014, ARA&A, 52, 291 [NASA ADS] [CrossRef] [Google Scholar]

[66] Conselice, C. J., Bershady, M. A., & Jangren, A. 2000, ApJ, 529, 886 [NASA ADS] [CrossRef] [Google Scholar]

[67] Conselice, C. J., Mortlock, A., Bluck, A. F. L., Grützbauch, R., & Duncan, K. 2013, MNRAS, 430, 1051 [Google Scholar]

[68] Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937 [NASA ADS] [CrossRef] [Google Scholar]

[69] Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11 [NASA ADS] [CrossRef] [Google Scholar]

[70] Crutcher, R. M. 2012, ARA&A, 50, 29 [NASA ADS] [CrossRef] [Google Scholar]

[71] Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944 [Google Scholar]

[72] Daddi, E., Bournaud, F., Walter, F., et al. 2010a, ApJ, 713, 686 [Google Scholar]

[73] Daddi, E., Elbaz, D., Walter, F., et al. 2010b, ApJ, 714, L118 [Google Scholar]

[74] Dalgarno, A., & McCray, R. A. 1972, ARA&A, 10, 375 [NASA ADS] [CrossRef] [Google Scholar]

[75] Dalla Vecchia, C., & Schaye, J. 2008, MNRAS, 387, 1431 [Google Scholar]

[76] Dashyan, G., & Dubois, Y. 2020, A&A, 638, A123 [EDP Sciences] [Google Scholar]

[77] Davé, R., Rafieferantsoa, M. H., Thompson, R. J., & Hopkins, P. F. 2017, MNRAS, 467, 115 [NASA ADS] [Google Scholar]

[78] Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827 [NASA ADS] [CrossRef] [Google Scholar]

[79] Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[80] De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2 [NASA ADS] [CrossRef] [Google Scholar]

[81] De Rossi, M. E., Bower, R. G., Font, A. S., Schaye, J., & Theuns, T. 2017, MNRAS, 472, 3354 [Google Scholar]

[82] Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2 [NASA ADS] [CrossRef] [Google Scholar]

[83] Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870 [Google Scholar]

[84] Dekel, A., & Silk, J. 1986, ApJ, 303, 39 [Google Scholar]

[85] Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451 [Google Scholar]

[86] Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]

[87] Dobbs, C. L., Burkert, A., & Pringle, J. E. 2011, MNRAS, 413, 2935 [Google Scholar]

[88] Donnari, M., Pillepich, A., Nelson, D., et al. 2020, MNRAS, submitted [arXiv:2008.00004] [Google Scholar]

[89] Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891 [NASA ADS] [CrossRef] [Google Scholar]

[90] D’Souza, R., Vegetti, S., & Kauffmann, G. 2015, MNRAS, 454, 4027 [Google Scholar]

[91] Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[92] Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2010, MNRAS, 409, 985 [Google Scholar]

[93] Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662 [NASA ADS] [CrossRef] [Google Scholar]

[94] Dubois, Y., Gavazzi, R., Peirani, S., & Silk, J. 2013, MNRAS, 433, 3297 [Google Scholar]

[95] Dubois, Y., Pichon, C., Welker, C., et al. 2014a, MNRAS, 444, 1453 [NASA ADS] [CrossRef] [Google Scholar]

[96] Dubois, Y., Volonteri, M., & Silk, J. 2014b, MNRAS, 440, 1590 [Google Scholar]

[97] Dubois, Y., Volonteri, M., Silk, J., Devriendt, J., & Slyz, A. 2014c, MNRAS, 440, 2333 [Google Scholar]

[98] Dubois, Y., Volonteri, M., Silk, J., et al. 2015, MNRAS, 452, 1502 [NASA ADS] [CrossRef] [Google Scholar]

[99] Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948 [Google Scholar]

[100] Dutton, A. A., Obreja, A., Wang, L., et al. 2017, MNRAS, 467, 4937 [Google Scholar]

[101] Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[102] Elmegreen, B. G., & Burkert, A. 2010, ApJ, 712, 294 [Google Scholar]

[103] Elmegreen, B. G., & Elmegreen, D. M. 2005, ApJ, 627, 632 [Google Scholar]

[104] Emsellem, E., Cappellari, M., Krajnović, D., et al. 2011, MNRAS, 414, 888 [Google Scholar]

[105] Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, ApJ, 644, 813 [NASA ADS] [CrossRef] [Google Scholar]

[106] Faucher-Giguère, C.-A. 2018, MNRAS, 473, 3717 [Google Scholar]

[107] Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156 [NASA ADS] [CrossRef] [Google Scholar]

[108] Ferrero, I., Navarro, J. F., Abadi, M. G., et al. 2017, MNRAS, 464, 4736 [Google Scholar]

[109] Fontana, A., Dunlop, J. S., Paris, D., et al. 2014, A&A, 570, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[110] Fossati, M., Wilman, D. J., Mendel, J. T., et al. 2017, ApJ, 835, 153 [Google Scholar]

[111] Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486 [Google Scholar]

[112] Gallazzi, A., Bell, E. F., Zibetti, S., Brinchmann, J., & Kelson, D. D. 2014, ApJ, 788, 72 [Google Scholar]

[113] Geen, S., Rosdahl, J., Blaizot, J., Devriendt, J., & Slyz, A. 2015, MNRAS, 448, 3248 [Google Scholar]

[114] Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175 [Google Scholar]

[115] Gentry, E. S., Krumholz, M. R., Madau, P., & Lupi, A. 2019, MNRAS, 483, 3647 [Google Scholar]

[116] Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101 [Google Scholar]

[117] Glowacki, M., Elson, E., & Davé, R. 2020, MNRAS, 498, 3687 [Google Scholar]

[118] Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223 [NASA ADS] [CrossRef] [Google Scholar]

[119] González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJ, 735, L34 [Google Scholar]

[120] Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[121] Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [CrossRef] [Google Scholar]

[122] Haardt, F., & Madau, P. 1996, ApJ, 461, 20 [NASA ADS] [CrossRef] [Google Scholar]

[123] Habouzit, M., Volonteri, M., & Dubois, Y. 2017, MNRAS, 468, 3935 [NASA ADS] [CrossRef] [Google Scholar]

[124] Heiles, C., & Troland, T. H. 2005, ApJ, 624, 773 [NASA ADS] [CrossRef] [Google Scholar]

[125] Hennebelle, P., & Chabrier, G. 2011, ApJ, 743, L29 [Google Scholar]

[126] Hopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, 417, 950 [NASA ADS] [CrossRef] [Google Scholar]

[127] Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522 [Google Scholar]

[128] Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581 [Google Scholar]

[129] Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800 [NASA ADS] [CrossRef] [Google Scholar]

[130] Hughes, A., Meidt, S. E., Colombo, D., et al. 2013, ApJ, 779, 46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[131] Iffrig, O., & Hennebelle, P. 2015, A&A, 576, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[132] Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[133] Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[134] Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[135] Inoue, S., Dekel, A., Mandelker, N., et al. 2016, MNRAS, 456, 2052 [NASA ADS] [CrossRef] [Google Scholar]

[136] Islam, R. R., Taylor, J. E., & Silk, J. 2003, MNRAS, 340, 647 [NASA ADS] [CrossRef] [Google Scholar]

[137] Iwamoto, K., Mazzali, P. A., Nomoto, K., et al. 1998, Nature, 395, 672 [Google Scholar]

[138] Jackson, R. A., Martin, G., Kaviraj, S., et al. 2020, MNRAS, 494, 5568 [CrossRef] [Google Scholar]

[139] Jackson, R. A., Kaviraj, S., Martin, G., et al. 2021a, MNRAS, 502, 1785 [CrossRef] [Google Scholar]

[140] Jackson, R. A., Martin, G., Kaviraj, S., et al. 2021b, MNRAS, 502, 4262 [CrossRef] [Google Scholar]

[141] Joung, M. R., Mac Low, M.-M., & Bryan, G. L. 2009, ApJ, 704, 137 [NASA ADS] [CrossRef] [Google Scholar]

[142] Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L35 [NASA ADS] [CrossRef] [Google Scholar]

[143] Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, ApJ, 758, 106 [NASA ADS] [CrossRef] [Google Scholar]

[144] Kassin, S. A., Brooks, A., Governato, F., Weiner, B. J., & Gardner, J. P. 2014, ApJ, 790, 89 [NASA ADS] [CrossRef] [Google Scholar]

[145] Kaviraj, S. 2014, MNRAS, 440, 2944 [Google Scholar]

[146] Kaviraj, S., Schawinski, K., Devriendt, J. E. G., et al. 2007, ApJS, 173, 619 [NASA ADS] [CrossRef] [Google Scholar]

[147] Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739 [NASA ADS] [Google Scholar]

[148] Kaviraj, S., Martin, G., & Silk, J. 2019, MNRAS, 489, L12 [Google Scholar]

[149] Kennicutt, R. C., Jr 1998, ApJ, 498, 541 [Google Scholar]

[150] Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2 [Google Scholar]

[151] Kim, C.-G., Ostriker, E. C., & Raileanu, R. 2017, ApJ, 834, 25 [NASA ADS] [CrossRef] [Google Scholar]

[152] Kimm, T., & Cen, R. 2014, ApJ, 788, 121 [NASA ADS] [CrossRef] [Google Scholar]

[153] Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015, MNRAS, 451, 2900 [NASA ADS] [CrossRef] [Google Scholar]

[154] Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826 [NASA ADS] [Google Scholar]

[155] Kimm, T., Haehnelt, M., Blaizot, J., et al. 2018, MNRAS, 475, 4617 [NASA ADS] [CrossRef] [Google Scholar]

[156] King, A. R., Lubow, S. H., Ogilvie, G. I., & Pringle, J. E. 2005, MNRAS, 363, 49 [NASA ADS] [CrossRef] [Google Scholar]

[157] Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[158] Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145 [NASA ADS] [CrossRef] [Google Scholar]

[159] Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [NASA ADS] [CrossRef] [Google Scholar]

[160] Krajnović, D., Alatalo, K., Blitz, L., et al. 2013, MNRAS, 432, 1768 [Google Scholar]

[161] Krajnović, D., Weilbacher, P. M., Urrutia, T., et al. 2015, MNRAS, 452, 2 [NASA ADS] [CrossRef] [Google Scholar]

[162] Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, Astron. Lett., 44, 8 [NASA ADS] [CrossRef] [Google Scholar]

[163] Krumholz, M. R., & Burkhart, B. 2016, MNRAS, 458, 1671 [NASA ADS] [CrossRef] [Google Scholar]

[164] Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250 [NASA ADS] [CrossRef] [Google Scholar]

[165] Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716 [NASA ADS] [CrossRef] [Google Scholar]

[166] Lagos, C. D. P., Crain, R. A., Schaye, J., et al. 2015, MNRAS, 452, 3815 [CrossRef] [Google Scholar]

[167] Lapiner, S., Dekel, A., & Dubois, Y. 2021, MNRAS, 505, 172 [CrossRef] [Google Scholar]

[168] Lee, J., Shin, J., Snaith, O. N., et al. 2021, ApJ, 908, 11 [CrossRef] [Google Scholar]

[169] Lelli, F., McGaugh, S. S., Schombert, J. M., Desmond, H., & Katz, H. 2019, MNRAS, 484, 3267 [NASA ADS] [CrossRef] [Google Scholar]

[170] Lupi, A., Haardt, F., & Dotti, M. 2015, MNRAS, 446, 1765 [NASA ADS] [CrossRef] [Google Scholar]

[171] Lupi, A., Bovino, S., Capelo, P. R., Volonteri, M., & Silk, J. 2018, MNRAS, 474, 2884 [CrossRef] [Google Scholar]

[172] Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [NASA ADS] [CrossRef] [Google Scholar]

[173] Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[174] Mandelker, N., Dekel, A., Ceverino, D., et al. 2014, MNRAS, 443, 3675 [NASA ADS] [CrossRef] [Google Scholar]

[175] Mannucci, F., Cresci, G., Maiolino, R., et al. 2009, MNRAS, 398, 1915 [NASA ADS] [CrossRef] [Google Scholar]

[176] Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115 [Google Scholar]

[177] Martin, G., Kaviraj, S., Devriendt, J. E. G., Dubois, Y., & Pichon, C. 2018a, MNRAS, 480, 2266 [NASA ADS] [CrossRef] [Google Scholar]

[178] Martin, G., Kaviraj, S., Devriendt, J. E. G., et al. 2018b, MNRAS, 474, 3140 [CrossRef] [Google Scholar]

[179] Martin, G., Kaviraj, S., Laigle, C., et al. 2019, MNRAS, 485, 796 [NASA ADS] [CrossRef] [Google Scholar]

[180] Martin, G., Kaviraj, S., Hocking, A., Read, S. C., & Geach, J. E. 2020, MNRAS, 491, 1408 [Google Scholar]

[181] Martin, G., Jackson, R. A., Kaviraj, S., et al. 2021, MNRAS, 500, 4937 [CrossRef] [Google Scholar]

[182] Martín-Navarro, I., Vazdekis, A., La Barbera, F., et al. 2015, ApJ, 806, L31 [NASA ADS] [CrossRef] [Google Scholar]

[183] Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364 [NASA ADS] [CrossRef] [Google Scholar]

[184] McAlpine, S., Bower, R. G., Rosario, D. J., et al. 2018, MNRAS, 481, 3118 [NASA ADS] [CrossRef] [Google Scholar]

[185] McGaugh, S. S. 2012, AJ, 143, 40 [Google Scholar]

[186] McGaugh, S. S., & Schombert, J. M. 2015, ApJ, 802, 18 [CrossRef] [Google Scholar]

[187] McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083 [Google Scholar]

[188] Menezes, R. B., Steiner, J. E., & Ricci, T. V. 2014, ApJ, 796, L13 [CrossRef] [Google Scholar]

[189] Menezes, R. B., Steiner, J. E., & da Silva, P. 2016, ApJ, 817, 150 [CrossRef] [Google Scholar]

[190] Miller, B. P., Gallo, E., Greene, J. E., et al. 2015, ApJ, 799, 98 [CrossRef] [Google Scholar]

[191] Moody, C. E., Romanowsky, A. J., Cox, T. J., Novak, G. S., & Primack, J. R. 2014, MNRAS, 444, 1475 [NASA ADS] [CrossRef] [Google Scholar]

[192] Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121 [Google Scholar]

[193] Mowla, L., van der Wel, A., van Dokkum, P., & Miller, T. B. 2019, ApJ, 872, L13 [NASA ADS] [CrossRef] [Google Scholar]

[194] Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18 [Google Scholar]

[195] Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59 [Google Scholar]

[196] Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJ, 699, L178 [NASA ADS] [CrossRef] [Google Scholar]

[197] Naab, T., Oser, L., Emsellem, E., et al. 2014, MNRAS, 444, 3357 [NASA ADS] [CrossRef] [Google Scholar]

[198] Narayanan, D., Bothwell, M., & Davé, R. 2012, MNRAS, 426, 1178 [NASA ADS] [CrossRef] [Google Scholar]

[199] Nguyen, D. D., Seth, A. C., Neumayer, N., et al. 2018, ApJ, 858, 118 [NASA ADS] [CrossRef] [Google Scholar]

[200] Nickerson, S., Teyssier, R., & Rosdahl, J. 2019, MNRAS, 484, 1238 [Google Scholar]

[201] Nomoto, K., Tominaga, N., Umeda, H., Kobayashi, C., & Maeda, K. 2006, Nucl. Phys. A, 777, 424 [NASA ADS] [CrossRef] [Google Scholar]

[202] Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[203] Nuñez-Castiñeyra, A., Nezri, E., Devriendt, J., & Teyssier, R. 2021, MNRAS, 501, 62 [Google Scholar]

[204] Ocvirk, P., Pichon, C., & Teyssier, R. 2008, MNRAS, 390, 1326 [Google Scholar]

[205] Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312 [NASA ADS] [CrossRef] [Google Scholar]

[206] Ostriker, E. C. 1999, ApJ, 513, 252 [NASA ADS] [CrossRef] [Google Scholar]

[207] Ostriker, E. C., & Shetty, R. 2011, ApJ, 731, 41 [NASA ADS] [CrossRef] [Google Scholar]

[208] Padoan, P., & Nordlund, Å. 2011, ApJ, 730, 40 [Google Scholar]

[209] Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2008, MNRAS, 391, 1117 [Google Scholar]

[210] Park, M.-J., Yi, S. K., Dubois, Y., et al. 2019, ApJ, 883, 25 [Google Scholar]

[211] Park, M. J., Yi, S. K., Peirani, S., et al. 2021, ApJS, 254, 2 [Google Scholar]

[212] Peñarrubia, J., Navarro, J. F., & McConnachie, A. W. 2008, ApJ, 673, 226 [Google Scholar]

[213] Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., & Bastian, N. 2018, MNRAS, 475, 4309 [NASA ADS] [CrossRef] [Google Scholar]

[214] Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101 [CrossRef] [Google Scholar]

[215] Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077 [Google Scholar]

[216] Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196 [NASA ADS] [CrossRef] [Google Scholar]

[217] Ponomareva, A. A., Verheijen, M. A. W., Papastergis, E., Bosma, A., & Peletier, R. F. 2018, MNRAS, 474, 4366 [Google Scholar]

[218] Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464 [Google Scholar]

[219] Popping, G., Somerville, R. S., & Trager, S. C. 2014, MNRAS, 442, 2398 [Google Scholar]

[220] Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14 [NASA ADS] [CrossRef] [Google Scholar]

[221] Pozzetti, L., Bolzonella, M., Lamareille, F., et al. 2007, A&A, 474, 443 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[222] Prantzos, N., Abia, C., Limongi, M., Chieffi, A., & Cristallo, S. 2018, MNRAS, 476, 3432 [Google Scholar]

[223] Prieto, J., Escala, A., Volonteri, M., & Dubois, Y. 2017, ApJ, 836, 216 [CrossRef] [Google Scholar]

[224] Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82 [NASA ADS] [CrossRef] [Google Scholar]

[225] Reines, A. E., Condon, J. J., Darling, J., & Greene, J. E. 2020, ApJ, 888, 36 [CrossRef] [Google Scholar]

[226] Reynolds, C. S. 2013, Class. Quant. Grav., 30, 244004 [NASA ADS] [CrossRef] [Google Scholar]

[227] Rezzolla, L., Barausse, E., Dorband, E. N., et al. 2008, Phys. Rev. D, 78, 044002 [CrossRef] [Google Scholar]

[228] Rosdahl, J., & Blaizot, J. 2012, MNRAS, 423, 344 [NASA ADS] [CrossRef] [Google Scholar]

[229] Rosdahl, J., & Teyssier, R. 2015, MNRAS, 449, 4380 [NASA ADS] [CrossRef] [Google Scholar]

[230] Rosdahl, J., Schaye, J., Dubois, Y., Kimm, T., & Teyssier, R. 2017, MNRAS, 466, 11 [NASA ADS] [CrossRef] [Google Scholar]

[231] Roškar, R., Teyssier, R., Agertz, O., Wetzstein, M., & Moore, B. 2014, MNRAS, 444, 2837 [CrossRef] [Google Scholar]

[232] Saftly, W., Baes, M., De Geyter, G., et al. 2015, A&A, 576, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[233] Salem, M., & Bryan, G. L. 2014, MNRAS, 437, 3312 [NASA ADS] [CrossRef] [Google Scholar]

[234] Sales, L. V., Navarro, J. F., Oman, K., et al. 2017, MNRAS, 464, 2419 [CrossRef] [Google Scholar]

[235] Sánchez, S. F., Barrera-Ballesteros, J. K., López-Cobá, C., et al. 2019, MNRAS, 484, 3042 [CrossRef] [Google Scholar]

[236] Sandage, A., Tammann, G. A., & Yahil, A. 1979, ApJ, 232, 352 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[237] Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2018, ApJ, 858, 99 [NASA ADS] [CrossRef] [Google Scholar]

[238] Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[239] Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521 [Google Scholar]

[240] Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150 [NASA ADS] [CrossRef] [Google Scholar]

[241] Sedgwick, T. M., Baldry, I. K., James, P. A., & Kelvin, L. S. 2019a, IAU Symposium 355: The Realm of the Low Surface Brightness Universe [Google Scholar]

[242] Sedgwick, T. M., Baldry, I. K., James, P. A., & Kelvin, L. S. 2019b, MNRAS, 484, 5278 [Google Scholar]

[243] Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2018, ApJ, 861, 4 [NASA ADS] [CrossRef] [Google Scholar]

[244] Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]

[245] She, R., Ho, L. C., & Feng, H. 2017, ApJ, 842, 131 [CrossRef] [Google Scholar]

[246] Shen, Y., Hwang, H.-C., Zakamska, N., & Liu, X. 2019, ApJ, 885, L4 [CrossRef] [Google Scholar]

[247] Silk, J. 2017, ApJ, 839, L13 [NASA ADS] [CrossRef] [Google Scholar]

[248] Silk, J., & Rees, M. J. 1998, A&A, 331, L1 [NASA ADS] [Google Scholar]

[249] Simons, R. C., Kassin, S. A., Weiner, B. J., et al. 2015, MNRAS, 452, 986 [NASA ADS] [CrossRef] [Google Scholar]

[250] Simons, R. C., Kassin, S. A., Weiner, B. J., et al. 2017, ApJ, 843, 46 [NASA ADS] [CrossRef] [Google Scholar]

[251] Smith, R., Choi, H., Lee, J., et al. 2016, ApJ, 833, 109 [NASA ADS] [CrossRef] [Google Scholar]

[252] Smith, A., Bromm, V., & Loeb, A. 2017, MNRAS, 464, 2963 [NASA ADS] [CrossRef] [Google Scholar]

[253] Soares, G., & Nemmen, R. 2020, MNRAS, 495, 981 [CrossRef] [Google Scholar]

[254] Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5 [NASA ADS] [CrossRef] [Google Scholar]

[255] Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [Google Scholar]

[256] Sun, J., Leroy, A. K., Schruba, A., et al. 2018, ApJ, 860, 172 [NASA ADS] [CrossRef] [Google Scholar]

[257] Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 [NASA ADS] [CrossRef] [Google Scholar]

[258] Tacconi, L. J., Neri, R., Chapman, S. C., et al. 2006, ApJ, 640, 228 [Google Scholar]

[259] Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246 [Google Scholar]

[260] Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781 [Google Scholar]

[261] Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74 [Google Scholar]

[262] Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179 [NASA ADS] [CrossRef] [Google Scholar]

[263] Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157 [CrossRef] [Google Scholar]

[264] Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[265] Teyssier, R., Moore, B., Martizzi, D., Dubois, Y., & Mayer, L. 2011, MNRAS, 414, 195 [CrossRef] [Google Scholar]

[266] Thorne, K. S. 1974, ApJ, 191, 507 [Google Scholar]

[267] Thornton, K., Gaudlitz, M., Janka, H. T., & Steinmetz, M. 1998, ApJ, 500, 95 [NASA ADS] [CrossRef] [Google Scholar]

[268] Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2016, ApJ, 817, 118 [NASA ADS] [CrossRef] [Google Scholar]

[269] Toomre, A., & Toomre, J. 1972, ApJ, 178, 623 [Google Scholar]

[270] Toro, E. 1999, Riemann Solvers and Numerical Methods for Fluid Dynamics (Berlin: Springer-Verlag) [CrossRef] [Google Scholar]

[271] Trebitsch, M., Blaizot, J., Rosdahl, J., Devriendt, J., & Slyz, A. 2017, MNRAS, 470, 224 [NASA ADS] [CrossRef] [Google Scholar]

[272] Trebitsch, M., Volonteri, M., Dubois, Y., & Madau, P. 2018, MNRAS, 478, 5607 [CrossRef] [Google Scholar]

[273] Trebitsch, M., Dubois, Y., Volonteri, M., et al. 2020, A&A, submitted [arXiv:2002.04045] [Google Scholar]

[274] Tremmel, M., Karcher, M., Governato, F., et al. 2017, MNRAS, 470, 1121 [CrossRef] [Google Scholar]

[275] Tremmel, M., Governato, F., Volonteri, M., Pontzen, A., & Quinn, T. R. 2018, ApJ, 857, L22 [Google Scholar]

[276] Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898 [NASA ADS] [CrossRef] [Google Scholar]

[277] Treu, T., Auger, M. W., Koopmans, L. V. E., et al. 2010, ApJ, 709, 1195 [Google Scholar]

[278] Trujillo, I., Förster Schreiber, N. M., Rudnick, G., et al. 2006, ApJ, 650, 18 [NASA ADS] [CrossRef] [Google Scholar]

[279] Tully, R. B., & Fisher, J. R. 1977, A&A, 500, 105 [Google Scholar]

[280] Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, MNRAS, 423, 2374 [NASA ADS] [CrossRef] [Google Scholar]

[281] van de Sande, J., Bland-Hawthorn, J., Brough, S., et al. 2017, MNRAS, 472, 1272 [NASA ADS] [CrossRef] [Google Scholar]

[282] van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28 [NASA ADS] [CrossRef] [Google Scholar]

[283] Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177 [NASA ADS] [CrossRef] [Google Scholar]

[284] Volonteri, M., Lodato, G., & Natarajan, P. 2008, MNRAS, 383, 1079 [NASA ADS] [CrossRef] [Google Scholar]

[285] Volonteri, M., Dubois, Y., Pichon, C., & Devriendt, J. 2016, MNRAS, 460, 2979 [Google Scholar]

[286] Volonteri, M., Pfister, H., Beckmann, R. S., et al. 2020, MNRAS, 498, 2219 [Google Scholar]

[287] Weiner, B. J., Willmer, C. N. A., Faber, S. M., et al. 2006, ApJ, 653, 1027 [NASA ADS] [CrossRef] [Google Scholar]

[288] Wetzel, A. R., Tinker, J. L., Conroy, C., & van den Bosch, F. C. 2013, MNRAS, 432, 336 [NASA ADS] [CrossRef] [Google Scholar]

[289] Yabe, K., Ohta, K., Akiyama, M., et al. 2015, PASJ, 67, 102 [CrossRef] [Google Scholar]

[290] Yu, S., Bullock, J. S., Wetzel, A., et al. 2020, MNRAS, 494, 1539 [CrossRef] [Google Scholar]

[291] Zahid, H. J., Kewley, L. J., & Bresolin, F. 2011, ApJ, 730, 137 [NASA ADS] [CrossRef] [Google Scholar]

[292] Zahid, H. J., Kashino, D., Silverman, J. D., et al. 2014, ApJ, 792, 75 [NASA ADS] [CrossRef] [Google Scholar]

[293] Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327 [NASA ADS] [CrossRef] [Google Scholar]