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?

Horizon simulation

Basic properties

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 σ8 = 0.81, baryon density Ωb = 0.045, Hubble constant H0 = 70.4  km s − 1 Mpc − 1, and ns = 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 Lbox = 100 h − 1Mpcwith 10243 dark matter (DM) particles, which results in a DM mass resolution of MDM, res = 8 × 107  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 aexp = 0.1, 0.2, 0.4 and 0.8).

Gas cooling and heating

Gas is allowed to cool by H and He cooling with a contribution from metals using a [31] model down to 104  K. Heating from a uniform UV background takes place after redshift zreion = 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.

Star formation and stellar feedback

The star formation process is modelled with a Schmidt law: ρ̇* = ϵ*ρ ⁄ tff, where ρ̇* is the star formation rate density, ϵ* = 0.02 [17, 25] the constant star formation efficiency, and tffthe local free-fall time of the gas. Star formation is allowed in regions which exceed a gas Hydrogen number density threshold of n0 = 0.1  H cm − 3following a Poissonian random process [24, 27]with a stellar mass resolution of M* = ρ0Δx3≃2 × 106  M. The gas pressure is artificially enhanced above ρ > ρ0 assuming a polytropic equation of state T = T0(ρ ⁄ ρ0)κ − 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].

Feedback from black holes

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 ρ > ρ0 with an initial seed mass of 105  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παG2M2 BHρ ⁄ (c2s + u2)3 ⁄ 2, where MBHis the BH mass, ρis the average gas density, csis the average sound speed, uis the average gas velocity relative to the BH velocity, and α is a dimensionless boost factor with α = (ρ ⁄ ρ0)2 when ρ > ρ0 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πGMBHmp ⁄ (ϵrσTc), where σTis the Thompson cross-section, c is the speed of light, mpis 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ϵrBHc2. 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 104  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.

Mock observations of galaxies

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 M20, and age.

Identifying and segmenting galaxies

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.

Synthetic colours

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.

Projected stellar kinematics

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): vpixel = (Σivlos, iIi, filter)/(ΣiIi, filter), where vlos, i is the velocity along the line of sight of the i-th star in the pixel considered, Ii, 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: σ2 pixel = (Σiv2 los, iIi, filter)/(ΣiIi, filter) − v2 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.

Specific star formation rate

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, Ms.


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

Spin of galaxies

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.

Tracing large-scale structures via the skeleton

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

Distribution of Matter

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

[1 ] Roberto G Abraham, Preethi Nair, Patrick J McCarthy, Karl Glazebrook, et al. The Gemini Deep Deep Survey. VIII. When Did Early-Type Galaxies Form?. \apj, 669(1):184—201, 2007.

[2 ] D Aubert, C Pichon, S Colombi. The origin and implications of dark matter anisotropic cosmic infall on =L haloes. \mnras, 352(2):376—398, 2004.

[3 ] J. M. Bardeen, J. R. Bond, N. Kaiser, A. S. Szalay. The statistics of peaks of Gaussian random fields. apj, 304:15-61, 1986.

[4 ] J. R. Bond, L. Kofman, D. Pogosyan. How filaments of galaxies are woven into the cosmic web. nat, 380:603-+, 1996.

[5 ] C. M. Booth, J. Schaye. Cosmological simulations of the growth of supermassive black holes and feedback from active galactic nuclei: method and tests. mnras, 398:53-74, 2009.

[6 ] G. Bruzual, S. Charlot. Stellar population synthesis at the resolution of 2003. mnras, 344:1000-1028, 2003.

[7 ] M. Colless, G. Dalton, S. Maddox, et al. The 2dF Galaxy Redshift Survey: spectra and redshifts. mnras, 328:1039-1063, 2001.

[8 ] D. J. Croton, V. Springel, S. D. M. White, G. De Lucia, C. S. Frenk, L. Gao, A. Jenkins, G. Kauffmann, J. F. Navarro, N. Yoshida. The many lives of active galactic nuclei: cooling flows, black holes and the luminosities and colours of galaxies. mnras, 365:11-28, 2006.

[9 ] Y. Dubois, J. Devriendt, A. Slyz, R. Teyssier. Jet-regulated cooling catastrophe. mnras, 409:985-1001, 2010.

[10 ] Y. Dubois, J. Devriendt, A. Slyz, R. Teyssier. Self-regulated growth of supermassive black holes by a dual jet-heating active galactic nucleus feedback mechanism: methods, tests and implications for cosmological simulations. mnras, 420:2662-2683, 2012.

[11 ] Y. Dubois, R. Gavazzi, S. Peirani, J. Silk. AGN-driven quenching of star formation: morphological and dynamical implications for early-type galaxies. mnras, 433:3297-3313, 2013.

[12 ] M. J. Geller, J. P. Huchra. Mapping the universe. Science, 246:897-903, 1989.

[13 ] L. Greggio, A. Renzini. The binary model for type I supernovae - Theoretical rates. aap, 118:217-222, 1983.

[14 ] F. Haardt, P. Madau. Radiative Transfer in a Clumpy Universe. II. The Ultraviolet Extragalactic Background. apj, 461:20-+, 1996.

[15 ] R. W. Hockney, J. W. Eastwood. Computer Simulation Using Particles. , 1981.

[16 ] J. Huchra, M. Davis, D. Latham, J. Tonry. A survey of galaxy redshifts. IV - The data. apjs, 52:89-119, 1983.

[17 ] Jr. Kennicutt. The Global Schmidt Law in Star-forming Galaxies. apj, 498:541-+, 1998.

[18 ] A. A. Klypin, S. F. Shandarin. Three-dimensional numerical model of the formation of large-scale structure in the Universe. mnras, 204:891-907, 1983.

[19 ] E. Komatsu, K. M. Smith, J. Dunkley, et al. Seven-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation. apjs, 192:18-+, 2011.

[20 ] C. Leitherer, D. Schaerer, J. D. Goldader, et al. Starburst99: Synthesis Models for Galaxies with Active Star Formation. apjs, 123:3-40, 1999.

[21 ] H. Omma, J. Binney, G. Bryan, A. Slyz. Heating cooling flows with jets. mnras, 348:1105-1119, 2004.

[22 ] Planck Collaboration, P. A. R. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barandela, J. Aumont, C. Baccigalupi, A. J. Banday, et al. Planck 2013 results. XVI. Cosmological parameters. ArXiv e-prints, 2013.

[23 ] S. Prunet, C. Pichon, D. Aubert, D. Pogosyan, R. Teyssier, S. Gottloeber. Initial Conditions For Large Cosmological Simulations. apjs, 178:179-188, 2008.

[24 ] Y. Rasera, R. Teyssier. The history of the baryon budget. Cosmic logistics in a hierarchical universe. aap, 445:1-27, 2006.

[25 ] Edwin E Salpeter. The Luminosity Function and Stellar Evolution. apj, 121:161, 1955.

[26 ] N. I. Shakura, R. A. Sunyaev. Black holes in binary systems. Observational appearance. aap, 24:337-355, 1973.

[27 ] S. F. Shandarin, Y. B. Zeldovich. The large-scale structure of the universe: Turbulence, intermittency, structures in a self-gravitating medium. Reviews of Modern Physics, 61:185-220, 1989.

[28 ] T. Sousbie, S. Colombi, C. Pichon. The fully connected N-dimensional skeleton: probing the evolution of the cosmic web. mnras, 393:457-477, 2009.

[29 ] T. Sousbie, C. Pichon, H. Kawahara. The persistent cosmic web and its filamentary structure - II. Illustrations. mnras, 414:384-403, 2011.

[30 ] V. Springel, L. Hernquist. Cosmological smoothed particle hydrodynamics simulations: a hybrid multiphase model for star formation. mnras, 339:289-311, 2003.

[31 ] R. S. Sutherland, M. A. Dopita. Cooling functions for low-density astrophysical plasmas. apjs, 88:253-327, 1993.

[32 ] R. Teyssier. Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES. aap, 385:337-364, 2002.

[33 ] E. F. Toro, M. Spruce, W. Speares. Restoration of the contact surface in the HLL-Riemann solver. Shock Waves, 4:25-34, 1994.

[34 ] D. Tweed, J. Devriendt, J. Blaizot, S. Colombi, A. Slyz. Building merger trees from cosmological N-body simulations. Towards improving galaxy formation models using subhaloes. aap, 506:647-660, 2009.

[35 ] Y. B. Zel'dovich. Gravitational instability: An approximate theory for large density perturbations. aap, 5:84-89, 1970.

[36 ] I. B. Zeldovich, J. Einasto, S. F. Shandarin. Giant voids in the universe. nat, 300:407-413, 1982.

[37 ] V. de Lapparent, M. J. Geller, J. P. Huchra. A slice of the universe. apjl, 302:L1-L5, 1986.