Theoretical models of structure formation by gravitational instability and numerical simulations have predicted that small fluctuations from the early Universe lead to the formation of a large-scale cosmic web made of clustered halos, filaments, sheets and voids [36, 18, 3, 4]. The resulting properties of the Universe’s large-scale structure is the interplay of the planar local collapse, as emphasized in [35] [27] and the inherent structure of the gaussian initial density and velocity shear fields, leading to the cosmic web picture of dense peaks connected by filaments, framing the honeycomb-like structure of walls [3, 4].

The extension of the Center for Astrophysics redshift survey [16] gave spectacular observational evidence [37, 12] for this picture, triggering a renewed interest for such large-scale galaxy surveys [7 , 13].

Modern simulations have established a tight connection between the geometry and dynamics of the large-scale structure of matter, on the one hand, and the evolution of the physical properties of forming galaxies, on the other. Observational information on the morphology of galaxies and its dependence on environment is routinely becoming available for galaxies up to redshift two and beyond [1, 15, 16]. Matched samples at low and high redshifts allow for the study of the evolution of many physical properties of galaxies for most of the history of our Universe in unprecedented detail. A key question formulated decades ago is nevertheless not satisfactorily answered: what properties of galaxies are driven by the cosmic environment?

1. Basic properties
   • Gas cooling and heating
   • Star formation and stellar feedback
   • Feedback from black holes
2. Mock observations of galaxies
   • Identifying and segmenting galaxies
   • Synthetics colours
   • Projected stellar kinematics
   • Specific stars formation rate
   • Ages
   • Spin of galaxies
   • Tracing large-scale structures via the skeleton
3. Distribution of Matter
4. References

Horizon-AGN adopts a standard ΛCDM cosmology with total matter density Ω_m = 0.272, dark energy density Ω_Λ = 0.728, amplitude of the matter power spectrum σ₈ = 0.81, baryon density Ω_b = 0.045, Hubble constant H₀ = 70.4  km s^− 1 Mpc^− 1, and n_s = 0.967 compatible with the WMAP-7 cosmology [19]. The values of this set of cosmological parameters are compatible with those of the recent Planck results within a ten per cent relative variation [22]. The size of the box is L_box = 100 h^− 1  Mpcwith 1024³ dark matter (DM) particles, which results in a DM mass resolution of M_DM, res = 8 × 10⁷  M_⊙. The initial conditions have been produced with the mpgrafic software [23]. The simulation was run down to z = 1.2 and used 4 million CPU hours.

The Horizon simulation is run with the Adaptive Mesh Refinement code ramses [32]. The evolution of the gas is followed using a second-order unsplit Godunov scheme for the Euler equations. The HLLC Riemann solver [33] with MinMod Total Variation Diminishing scheme is used to reconstruct the interpolated variables from their cell-centered values. Collisionless particles (DM and star particles) are evolved using a particle-mesh solver with a Cloud-In-Cell interpolation. The initial mesh is refined up to Δx = 1  kpc(7 levels of refinement). This is done according to a quasi-Lagragian criterion: if the number of DM particles in a cell is more than 8, or if the total baryonic mass in a cell is 8 times the initial DM mass resolution, a new refinement level is triggered. In order to keep the minimum cell size approximately constant in physical units, we allow a new maximum level of refinement every time the expansion scale factor doubles (i.e. at a_exp = 0.1, 0.2, 0.4 and 0.8).

Gas is allowed to cool by H and He cooling with a contribution from metals using a [31] model down to 10⁴  K. Heating from a uniform UV background takes place after redshift z_reion = 10following [14]. Metallicity is modelled as a passive variable for the gas and its amount is modified by the injection of gas ejecta during supernovae explosions and stellar winds. We also account for the release of various chemical elements synthesised in stars and released by stellar winds and supernovae: O, Fe, C, N, Mg and Si. However, they do not contribute separately to the cooling curve (the ratio between each element is taken to be solar for simplicity) but can be used to probe the distribution of the various metal elements. The gas follows an equation of state for an ideal monoatomic gas with an adiabatic index of γ = 5 ⁄ 3.

The star formation process is modelled with a Schmidt law: ρ̇_* = ϵ_*ρ ⁄ t_ff, where ρ̇_* is the star formation rate density, ϵ_* = 0.02 [17, 25] the constant star formation efficiency, and t_ffthe local free-fall time of the gas. Star formation is allowed in regions which exceed a gas Hydrogen number density threshold of n₀ = 0.1  H cm^− 3following a Poissonian random process [24, 27]with a stellar mass resolution of M_* = ρ₀Δx³≃2 × 10⁶  M_⊙. The gas pressure is artificially enhanced above ρ > ρ₀ assuming a polytropic equation of state T = T₀(ρ ⁄ ρ₀)^κ − 1 with polytropic index κ = 4 ⁄ 3 to avoid excessive gas fragmentation and mimic the effect of stellar heating on the mean temperature of the interstellar medium [30]. Feedback from stars is explicitly taken into account assuming a [25]initial mass function with a low-mass (high-mass) cut-off of 0.1  M_⊙(100  M_⊙), as described in detail in Kimm et al. (in prep.). Specifically, the mechanical energy from supernovae type II and stellar winds is taken from starburst99 [20, 31], and the frequency of supernovae type Ia explosions is computed following [13].

The same “canonical” Active Galactic Nuclei (AGN) feedback modelling employed in [10] is used here. Black holes (BHs) are created where the gas mass density is larger than ρ > ρ₀ with an initial seed mass of 10⁵  M_⊙. In order to avoid the formation of multiple BHs in the same galaxy, BHs are not allowed to form at distances less than 50 kpc from each other. The accretion rate onto BHs follows the Bondi-Hoyle-Lyttleton rate Ṁ_BH = 4παG²M²_BHρ ⁄ (c²_s + u²)^3 ⁄ 2, where M_BHis the BH mass, ρis the average gas density, c_sis the average sound speed, uis the average gas velocity relative to the BH velocity, and α is a dimensionless boost factor with α = (ρ ⁄ ρ₀)² when ρ > ρ₀ and α = 1 otherwise [5]in order to account for our inability to capture the colder and higher density regions of the inter-stellar medium. The effective accretion rate onto BHs is capped at the Eddington accretion rate: Ṁ_Edd = 4πGM_BHm_p ⁄ (ϵ_rσ_Tc), where σ_Tis the Thompson cross-section, c is the speed of light, m_pis the proton mass, and ϵ_r is the radiative efficiency, assumed to be equal to ϵ_r = 0.1for the [26] accretion onto a Schwarzschild BH.

The AGN feedback is a combination of two different modes, the so-called radio mode operating when χ = Ṁ_BH ⁄ Ṁ_Edd < 0.01and the quasar mode active otherwise. The quasar mode consists of an isotropic injection of thermal energy into the gas within a sphere of radius Δx, and at an energy deposition rate: Ė_AGN = ϵ_fϵ_rṀ_BHc². In this equation, ϵ_f = 0.15is a free parameter chosen to reproduce the scaling relations between BH mass and galaxy properties (mass, velocity dispersion) and BH density in our local Universe (see [10]). At low accretion rates, the radio mode deposits AGN feedback energy into a bipolar outflow with a jet velocity of 10⁴  km s^− 1. The outflow is modelled as a cylinder with a cross-sectional radius Δx and height 2 Δx following [21] (more details are given in [9]). The efficiency of the radio mode is larger than the quasar mode with ϵ_f = 1.

One can discern the large-scale pattern of the cosmic web, with filaments and walls surrounding voids and connecting halos. Massive halos are filled with hot gas, and feedback from supernovae and AGN pours warm and metal-rich gas in the diffuse inter-galactic medium. As demonstrated in [11], the modelling of AGN feedback is critical to create early-type galaxies and provide the sought morphological diversity in hydrodynamical cosmological simulations [8] for semi-analytical models.

We describe how we produce various observables that can be compared qualitatively with data from modern observational surveys. In this paper we focus on observables which are known to correlate with the Hubble type of galaxies, namely mass, V ⁄ σ, colour, morphological parameters like Gini and M₂₀, and age.

Galaxies are identified with the AdaptaHOP finder ([2], updated to its recent version by [34] for building merger trees) which directly operates on the distribution of star particles. A total of 20 neighbours are used to compute the local density of each particle, a local threshold of ρ_t = 178times the average total matter density is applied to select relevant densities, and the force softening (minimum size below which substructures are considered irrelevant) is    2 kpc. Only galactic structures identified with more than 50 particles are considered. It allows for a clear separation of galaxies (defined as sets of star particles segmented by AdaptaHOP), including those in the process of merging. Catalogues of around    150 000 galaxies are produced for each redshift analysed in this paper from z = 1.2 to z = 3.

We compute the absolute AB magnitudes and rest-frame colours of galaxies using single stellar population models from [6] assuming a Salpeter Initial Mass Function (IMF). Each star particle contributes to a flux per frequency that depends on its mass, age, and metallicity. The sum of the contribution from all stars is passed through the u, g, r, and i filters from the SDSS. Fluxes are expressed as rest-frame quantities (i.e. that do not take into account the red-shifting of spectra). We also neglect the contribution to the reddening of spectra from internal (interstellar medium) or external (intergalactic medium) dust extinction. Once the flux in each waveband is obtained for a star particle, we build two-dimensional projected maps from single galaxies (satellites are excised with the galaxy finder), and we can sum up the total contribution of their stars to the total luminosity.

For each galaxy, we build a field of view centred on the galaxy, which is made of 256 × 256 pixels over 100 kpc size (corresponding to a pixel size of 0.4 kpc, or 0.05 arcsec at z = 1.83). We compute the luminosity-weighted velocity along the line of sight (arbitrary defined as the x-axis of the simulation): v_pixel = (Σ_iv_los, iI_i, filter)/(Σ_iI_i, filter), where v_los, i is the velocity along the line of sight of the i-th star in the pixel considered, I_i, filteris the intensity in the corresponding filter bandwidth (u, g, r, i) of the i-th star in the pixel considered. Then, the velocity dispersion along the line of sight is: σ²_pixel = (Σ_iv²_los, iI_i, filter)/(Σ_iI_i, filter) − v²_pixel , The velocity maps are then smoothed with a gaussian kernel of 15 pixels. The position of the fastest (respectively slowest) pixel, which defines V for that galaxy, is then identified automatically and a 0.75 arcsec “slit” is put across so as to interpolate through the kinematic major axis of the galaxy. The smoothed velocity dispersion map is also interpolated along the same axis and the maximum of that curve defines σ. V ⁄ σ is then straightforwardly the corresponding ratio.

The calculation of the star formation rate (SFR) is done on stars as identified by the galaxy finder that belong to a given galaxy. To compute the SFR, we compute the amount of stars formed over the last 100 Myr. The choice of 100  Myrcorresponds to a minimum measurable SFR of M_* ⁄ 100  Myr = 0.02  M_⊙ yr^− 1. The specific star formation rate (sSFR) is then calculated by diving the SFR by the galaxy stellar mass, M_s.

The mean ages of galaxies are obtained through the summation of the mass-weighted age of star particles belonging to the galaxy.

To compute the spin of galaxies, we compute the total angular momentum of their stars with respect to the particle of maximum density (centre of the galaxy) from the smoothed stellar density constructed with the AdaptaHOP algorithm.

We have also tested the effect of grid-locking on the cartesian axes of the box (a common issue of cartesian-based Poisson solvers for which a numerical anisotropy in the force calculation arises, [15]) for galaxies and filaments.

In order to quantify the orientation of galaxies relative to the cosmic web, we use a geometric three-dimensional ridge extractor well suited to identify filaments, called the “skeleton”. A gas density cube of 512³ pixels is drawn from the simulation and gaussian-smoothed with a length of 3 h^− 1  Mpccomoving chosen so as to trace large-scale filamentary features. Two implementations of the skeleton, based on “watershed” [28] and “persistence” [29] were implemented, without significant difference for the purpose of this investigation. The first method identifies ridges as the boundaries of walls which are themselves the boundaries of voids. The second one identifies ridges as the “special” lines connecting topologically robust (filament-like) saddle points to peaks.

The clustering of the galaxies follows quite closely the skeleton of the gas, i.e. the cosmic filaments. Note that, on large-scales, the skeleton built from the gas is equivalent to that built from the DM as the gas and DM trace each other closely. The rest of the paper is devoted to studying the orientation of the spin of these galaxies relative to the direction of the nearest skeleton segment. In practice, an octree is built from the position of the mid-segment of the skeleton to speed up the association of the galaxy position to its nearest skeleton segment. It was checked that our results were not sensitive to how many such segments were considered to define the local direction of the skeleton. The orientation of the segment of the skeleton is used to define the relative angle between the filament and the spin of the galaxy.

The segments are also tagged with their curvilinear distance to the closest node (where different filaments merge), which allows us to study the evolution of this (mis)alignment along the cosmic web.

Baryonic processes such as the cooling of gas, star formation and feedback from AGN, can have an impact on the distribution of matter. This is of particular relevance to cosmological studies with large-scale galaxy surveys, which look into extracting cosmological information from observables such as weak gravitational lensing. From the simulation, it is possible to measure the distribution of matter in a wide range of scales and make predictions for those observables. Complementary runs to Horizon-AGN (our full baryonic physics run), such as Horizon-DM (dark matter-only) and Horizon-noAGN (without AGN feedback), give us clues into the role of baryonic physics in shaping the distribution of matter. Our results are made available in this webpage