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

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 $\sim(16\,\rm Mpc)^3$ 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 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.


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 smallscale 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 formationwhich 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 selfconsistent 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 Article number, page 1 of 29 arXiv:2009.10578v4 [astro-ph.GA] 28 Jun 2021 A&A proofs: manuscript no. article 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. 2010Dubois et al. , 2012Kaviraj 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 highmass 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 largescale simulations include Horizon-AGN (Dubois et al. 2014a), Illustris , EAGLE , 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. 2014Hopkins et al. , 2018Dubois 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-turbulencedriven 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 1 , is to provide a complementary tool between these two standard techniques, that is between the few well-resolved objects versus 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) 3 . 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) 3 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 Section 4.

The NewHorizon simulation: Prescription
We describe the NewHorizon simulation employed in this work 2 , which is a sub-volume extracted from its parent Horizon-AGN simulation (Dubois et al. 2014a) 3 , 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.

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 3 DM particles, a 1024 3 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 σ 8 = 0.81, baryon density Ω b = 0.045, Hubble constant H 0 = 70.4 km s −1 Mpc −1 , 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 3 effective resolution, that is with a DM mass resolution of M DM,hr = 1.2 × 10 6 M . The high-resolution initial patch is embedded in buffered regions with decreasing mass resolution of 10 7 M , 8 × 10 7 M , 6 × 10 8 M for spheres of 10.6 Mpc, 11.7 Mpc, and 13.9 Mpc radius, respectively, and a resolution of 5 × 10 9 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 −3 .
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 8 star particles formed and completed 4.7 × 10 6 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. Fig. 1 shows a projection of the high-resolution region. Fig. 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 −3 ) 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 5 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).

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 4 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 4 K and those from Dalgarno & McCray (1972) below ≈ 10 4 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 −3 ) (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 −3 .
Compared to Horizon-AGN, the gas can now cool below 10 4 K.

Star formation
Star formation occurs in regions with hydrogen gas number density above n 0 = 10 H cm −3 (the stellar mass resolution is n 0 m p ∆x 3 = 1.3 × 10 4 M ) following a Schmidt law:ρ = ρ g /t ff , whereρ is the SFR mass density, ρ g the gas mass density, t ff = 3π/(32Gρ g ) the local free-fall time of the gas, G the gravitational constant, and ia varying star formation efficiency Trebitsch et al. 2017Trebitsch et al. , 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 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. The two top panels encompass the zoom-in region, with its network of filaments. The two bottom panels illustrate how narrow filaments break up and mix once they connect to one of the most massive galaxies of that zoom-in region. cloud turbulent Mach number M = 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: where s = ln(ρ/ρ 0 ) is the logarithmic density contrast of the PDF with mean ρ 0 and variance σ 2 s = ln(1+b 2 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: (2) 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 10 M at z = 1 seen edge-on (top panels) or face-on (bottom panels).
In the NewHorizon simulation, the turbulent Mach number is given by the local three-dimensional instantaneous velocity dispersion σ g (obtained by computing σ 2 g = sum(∇ ⊗ udx) 2 ), and the virial parameter also takes the thermal pressure support α vir,0 = 5(σ 2 g + c 2 s )/(πρ g G∆x 2 ) into account. In this case, φ −1 t = 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 φ −1 t 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 protostellar 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 (M > 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 NewHorizon 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 per cent efficiency.

Feedback from massive stars
We include feedback from Type II SNe assuming that each explosion initially releases the kinetic energy of 10 51 erg. Because the minimum mass of a star particle is 10 4 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 . We increase this number by a factor of 2 (0.03 M −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 momentumconserving 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 −2/17 51 n −4/17 1 Z −0.28 , where E 51 is the total energy released in units of 10 51 erg, n 1 is the hydrogen number density in units of cm −3 , and Z = max[Z/Z , 0.01] is the metallicity, normalised to Article number, page 5 of 29 A&A proofs: manuscript no. article the solar value (Z = 0.02). The final momentum in the snowplough phase per SN explosion is taken from Thornton et al. (1998) as q SN = 3 × 10 5 km s −1 M E 16/17 51 n −2/17 1 Z −0.14 . (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 q SN+PH = 5 × 10 5 km s −1 M E 16/17 51 n −2/17 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. However, e SN can be increased up to 3.6 × 10 49 erg M −1 if a non-negligible fraction ( f HN = 0.5) of hypernovae (with E HN 10 52 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 , 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 Section 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.

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 −1 , and that are located at a distance of at least 50 comoving kpc from any preexisting 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 whereū is the average MBH-to-gas relative velocity,c s the average gas sound speed, andρ 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 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 0 ) 2 when n > n 0 , 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).

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 corotating 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) (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 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.

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ṀMBH c 2 , 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 4 km s −1 , whose exact value has little impact on MBH growth or galaxy mass content (Dubois et al. 2012).
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. 2016), 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.

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/ √ N)δ t ρ c (where N is the number of particles in the (sub)structure; see Aubert 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 per cent 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 Adap-taHOP 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 7 , 10 8 , 10 9 , and 10 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 12 M .

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).

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 8 M at z = 4, 8.2 × 10 9 M at z = 2 and 5.5 × 10 10 M at z = 0.25) and the 15 panels on the right represent a less massive galaxy (9.7 × 10 7 M at z = 4, 1.4 × 10 9 M at z = 2 and 4.0 × 10 9 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 Section 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 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. The first row shows 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.
z ≈ 1.0 and maintains the marginally stable disc for the rest of the cosmic history.

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 red-shift 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, section 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 to 10 11 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) 3 . We could alternatively use the sum of each individual leaf cell that passes a given threshold value of the passive scalar colour value (see Section 2.1). The corresponding initial volume can be reduced by 20-40 per cent 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 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  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 regiona 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) Fontana+04 (K20, 0.7 < z < 1) Pozzetti+07 (VVDS, 0.9 < z < 1.2) Bielby+12 (WIRDS, 1 < z < 1.2) Ilbert+13 (UltraVISTA, 0.8 < z < 1.1) Tomczak+14 (ZFOURGE, 1 < z < 1.5) Davidzon+17 (COSMOS, 0.8 < z < 1.1) 10 7 10 8 10 9 10 10 10 11 Ms/Msun Pozzetti+07 (VVDS, 1.6 < z < 2.5) Ilbert+13 (UltraVISTA, 2 < z < 2.5) Tomczak+14 (ZFOURGE, 2 < z < 2.5) Davidzon+17 (COSMOS, 2 < z < 2.5) 10 7 10 8 10 9 10 10 10 11 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 per cent 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.

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  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.
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 .
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 red-shift 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 ) 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 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 (α = √ λ 1 /ΣI) and minor (β = √ λ 2 /ΣI) axes of the ellipse are obtained from the covariance matrix, where λ 1 and λ 2 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 10 (2)+2.5log 10 (A), where A = R 2 αβπ 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  use to derive their measurements.
In Fig. 7, we show the evolution of the surface brightness µ e,r versus 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 ). 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 8 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. Figure 7 shows that the surface brightness versus stellar mass plane in NewHorizon corresponds well to , 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 8 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 Section 3.7) can be tested using data from future instruments such as the LSST.  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. 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 4 . Since stellar particles lose mass owing to stellar feedback, we compensate for this mass loss when recon- structing 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 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.

Star formation rates
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 figure 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 5 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 9 M dominate the SFR and stellar mass budget in the early Universe, while intermediate mass galaxies 10 9 ≤ M s /M < 10 11 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 6 ; 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. 2013Dubois et al. , 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 Section 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 7 . 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 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 −1 . 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 9 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 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 9 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 C4 of Behroozi et al. 2019). Figure 10 shows the Kennicutt-Schmidt relation of surface density of SFR (Σ SFR ) as a function of surface density of total (HI+H 2 ) gas (Σ gas ) for galaxies in the NewHorizon simulation at redshifts z = 4, 2, 1, 0.25 and with M s > 10 7 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 2 gas is selected as gas with density n > 0.1 cm −3 and temperature T < 2×10 4 K. The observational data shown in Fig. 10 include z ∼ 1.5 BzKselected 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), IRluminous 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.

Kennicutt-Schmidt relation
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 −2 and high Σ SFR > 10 −2 M yr −1 kpc −2 , 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 − 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 −3 and temperature T < 2 × 10 4 K, the average Σ gas -Σ SFR relation at z < 3 follows the sequence of discs (Kraljic et al., in prep.).
Article number, page 13 of 29 A&A proofs: manuscript no. article 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.
3.6. Galaxy-to-halo mass relation Fig. 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. Adapta-HOP 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 Section 3.10 (typically 10 per cent and 1 per cent 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).   Behroozi et al. (2013) at the indicated redshift and dotdashed 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 stellarto-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 11 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 11 M .

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 , 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 7 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).
The large spread of simulated data points below M s < 10 8 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 highmass end (M s > 10 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. 2020). 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. 2013Dubois et al. , 2016.

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 (massweighted) observations by Panter et al. (2008) at z = 0.1. 8 We note that observational estimates of metallicities barely dif- fer 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 −3 and temperature T < 2 × 10 4 K. Gas metallicity is systematically larger by 50 per cent than the stellar metallicity at any given redshift as an effect of the stellar metallicity being composed of very old poorly enriched stars. 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. 2011Zahid et al. , 2014Yabe 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 Z = 0.01345 (Asplund et al. 2009). We have, therefore, scaled up their fitting relations accordingly. 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 per cent, and oxygen represents 43 per cent 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 NewHorizon z = 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).

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 = (σ 2 r +σ 2 t +σ 2 z )/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  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. decoupled cores (Krajnović et al. 2013(Krajnović et al. , 2015Coccato 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.
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 ro-bust 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 11 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.
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.

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 5 M . Therefore, it should be noted that a nonnegligible 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 Adapta-HOP 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 substruc- tures at 10 instead of 50; hence the minimum cluster mass is 10 5 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 Adap-taHOP 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 9 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 gasrich, compact, and turbulent galaxies that are increasingly more gravitationally unstable (Cacciato et al. 2012;Inoue et al. 2016), . Surface densities within 2R eff as a function of galaxy stellar mass at different redshifts for stars (top), HI+H 2 gas (n > 0.1 cm −3 and T < 2 × 10 4 K, middle), and H 2 molecular gas (n > 10 cm −3 and T < 2 × 10 4 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.
thereby forming stars into more numerous massive clusters and more efficiently; this is confirmed in sections 3.11 and 3.13.

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 −3 and temperature T < 2 × 10 4 K. In the following, we distinguish between the neutral HI+H 2 gas component with the density cut-off at 0.1 cm −3 , and the H 2 molecular dense component at 10 cm −3 (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 ) 2 . 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. 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 −3 ; 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 −3 ; 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 −3 ) material increases with stellar mass. Nonetheless, the fraction of cold dense gas decreases with decreasing redshift at a given stellar mass.
Observations of star-forming H 2 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. 2006Tacconi et al. , 2008Daddi et al. 2010a;Tacconi et al. 2010Tacconi et al. , 2013Scoville et al. 2017;Tacconi et al. 2018Tacconi et al. , 2020. Reported baryonic gas fractions range from ∼ 20-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. 2017Davé et al. , 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 − 10 12 M and a reduced scatter in stellar mass versus 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 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 11 M ) at 1.5 < z < 3. In the range of comparable stellar masses, gas fractions of highz 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 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 −3 (left and right, respectively) and temperature lower than 2 × 10 4 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 figure 12) at various redshifts and from Catinella et al. (2018) from local galaxies for the HI+H 2 gas, while the dotted line shows the data from Tacconi et al. (2018) at various redshifts for the H 2 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). 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 10 M , but significantly differ with the local scaling relation at the low-mass end.

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 −3 and temperatures below 2 × 10 4 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 9 M ; however, considering non-disc galaxies leads to a similar relation.  Fig. 21. Gas kinematics as a function of galaxy stellar mass at different redshifts for gas density 0.01 < n < 10 cm −3 , temperature lower than T < 2 × 10 4 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 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 ∝ M 1/m b ) 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 8 M b /M 10 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 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 (Mc-Gaugh 2012).

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 −2 ≤ n/cm −3 < 10 and T ≤ 2 × 10 4 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 σ 2 g = (σ 2 g,r + σ 2 g,z + σ 2 g,t )/3. Similarly the rotational velocity is obtained by taking mass-weighted mean tangential velocity component V g,rot =V g,t . A proxy for the total kinematics support is built out of the dispersion and rotation using S g,0.5 = (0.5V 2 g,rot + σ 2 g ) 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. versus σ g,z ), density cut-offs (i.e. versus selected gas with n > 0.1 cm −3 and T ≤ 2 × 10 4 K) or aperture (i.e. versus 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. 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 9 M , values as high as 60 km s −1 (and up to 100 km s −1 for the most massive galaxies at high redshift) at z = 4 and 20 km s −1 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 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 10 M ), but V g,rot /S g,0.5 at M s < 10 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 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(Kassin et al. , 2014Ceverino 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  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 9 M , with non-central MBHs generally failing to grow even in galaxies above this mass threshold. 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.

MBH-to-galaxy mass relation
In NewHorizon, MBH formation commonly occurs even in lowmass galaxies, as can be seen by the fact that about half of galaxies with a stellar mass 10 7 − 10 8 M at redshift z = 0.2 5 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 percent and 100 percent of galaxies with mass ∼ 10 9 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 percent of galaxies with mass ∼ 10 9 − 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   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.
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 6 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 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(Menezes et al. , 2016Shen 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. 2018Trebitsch et al. , 2020, we confirm that MBH growth is regulated by SN feedback for galaxies with stellar masses below M s < 5 × 10 9 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 11 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.
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 4 M in NewHorizon). It is only once their host galaxy reaches the transition stellar mass of M s > 5 × 10 9 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.

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 4 − 10 5 M . The MBHs are formed with zero spin and initial mass M BH = 10 4 M , and for masses up to 2 × 10 4 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 4 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 4 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 4 − 10 5 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 6 M have acquired more than 80% of the 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 6 − 5 × 10 7 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.).

Conclusions
The NewHorizon tool provides a compromise to alternative strategies, namely large-volume cosmological simulations with poor spatial resolution within galaxies or selected highresolution zoomed-in halos, by simulating an intermediate subvolume (16 Mpc) 3 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 multiphase structure with dense cold star-forming clumps embedded in warm and hot diffuse gas, where turbulence shape starforming properties of galaxies. -The galaxy mass function with the characteristic turnover at 10 10 − 10 11 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 11 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 starforming 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 onefourth. -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 9 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 starforming 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 largescale 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. 2012Dubois et al. , 2016Kimm 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.